Merge pull request #9 from philipplndwhr/try-replacing-lambdas-with-defs

Try replacing lambdas with defs

src/doptics/functions.py

@@ -3,21 +3,23 @@
 import scipy as sp
 from typing import Callable, List
 from numpy.typing import ArrayLike
+from icecream import ic
 
 YL = 0
 YR = 2
 SIGMA = 1
 XL = 1
 XR = 9
 
-uniform = lambda x: 1
-triangle = lambda x: x
-normal_10 = lambda x: np.exp(-0.5*((x-0)/SIGMA)**2) + 0.1
+
+def uniform(x): return 1
+def triangle(x): return x
+def normal_source(x, xl=XL, xr=XR, sigma=SIGMA): return np.exp(-0.5*((x-((xl+xr)*0.5))/sigma)**2) + 0.01
+def normal_target(x, yl=YL, yr=YR, sigma=SIGMA): return np.exp(-0.5*((x-((yl+yr)*0.5))/sigma)**2) + 0.01
+normal_10 = lambda x: np.exp(-0.5 * ((x - 0) / SIGMA) ** 2) + 0.1
 G1 = lambda y1: 1
 G2 = lambda y2: 1 - np.abs(y2 - 12) if 11 < y2 < 13 else 0
 E = lambda x: 1 / (np.exp(10 * (x - 0.5)) + np.exp(-10 * (x - 0.5)))
-normal_target = lambda x: np.exp(-0.5*((x-((YR+YL)*0.5))/SIGMA)**2) + 0.01
-normal_source = lambda x: np.exp(-0.5*((x-((XR+XL)*0.5))/SIGMA)**2) + 0.01
 
 
 def cdf_sampling_source(source_density: Callable, linear_samples: ArrayLike, total_samples: int = 1000):
@@ -48,14 +50,13 @@ def cdf_sampling_source(source_density: Callable, linear_samples: ArrayLike, tot
         Phi_arr[i] = Phi_arr[i - 1] + elt  # abuses np.zeros()
     Phi_arr = np.roll(Phi_arr, 1)
     Phi_arr[0] = 0
-    # return Phi_arr
     print(Phi_arr)
 
     distr_samples = np.zeros(len(linear_samples))
-    inverted_Phi = reversed(Phi_arr)
+    reversed_Phi = reversed(Phi_arr)
     quantiles = np.linspace(0, 1, len(linear_samples))
     for i in range(len(distr_samples)):
-        # ladder_index = next(x[0] for x in enumerate(inverted_Phi) if x[1] < delta*i)
+        # ladder_index = next(x[0] for x in enumerate(reversed_Phi) if x[1] < delta*i)
         # ladder_index = len([x for x in Phi_arr if x <= + i*delta]) - 1
         ladder_index = len([x for x in Phi_arr if x <= + quantiles[i]]) - 1
         distr_samples[i] = fine_ladder[ladder_index]
@@ -95,25 +96,31 @@ def construct_target_density_intervals_from_angular(angle_density: Callable,
         np.array([np.cos(large_angle), np.cos(small_angle) * np.sin(abs_large_angle)/np.sin(abs_small_angle)])
     )
 
-    print(y1_span)
-    print(y2_span)
+    ic(y1_span)
+    ic(y2_span)
 
     if l1 == l2:
         l2 = 1.4 * l2
         y2_span[0] = 1.4 * y2_span[0]
         y2_span[1] = 1.4 * y2_span[1]
 
-    print(y2_span)
+    ic(y2_span)
+
+    def y1_density(y1):
+        return (angle_density(np.arccos(y1 / np.linalg.norm(np.array(
+            [y1 - center[0], l1 - center[1]], dtype=object)))) /
+                np.linalg.norm(np.array([y1 - center[0], l1 - center[1]], dtype=object))
+                )
+
+    def y2_density(y2):
+        return (angle_density(np.arccos(y2 / np.linalg.norm(np.array(
+            [y2 - center[0], l2 - center[1]], dtype=object)))) /
+                np.linalg.norm(np.array([y2 - center[0], l2 - center[1]], dtype=object))
+                )
 
-    y1_density = lambda y1: (angle_density(np.arccos(y1 / np.linalg.norm(np.array(
-        [y1-center[0], l1-center[1]], dtype=object)))) /
-                             np.linalg.norm(np.array([y1-center[0], l1-center[1]], dtype=object)))
-    y2_density = lambda y2: (angle_density(np.arccos(y2 / np.linalg.norm(np.array(
-        [y2-center[0], l2-center[1]], dtype=object)))) /
-                             np.linalg.norm(np.array([y2-center[0], l2-center[1]], dtype=object)))
+    # for i in np.linspace(small_angle, large_angle, 100):
+    #     print(f'y2({i}) = {y2_density(i)}')
 
-    for i in np.linspace(small_angle, large_angle, 100):
-        print(f'y2({i}) = {y2_density(i)}')
     y1_span = y1_span + center[0]
     y2_span = y2_span + center[0]
     for j in np.linspace(y2_span[0], y2_span[1], 100):
@@ -123,15 +130,13 @@ def construct_target_density_intervals_from_angular(angle_density: Callable,
     return y1_density, y2_density, y1_span, y2_span, l1 + center[1], l2 + center[1]
 
 
-def f(x):
-    return 1
+def f(x): return 1
 
 
-def g(x):
-    return x**2
+def g(x, mu=0): return (x-mu)**2
 
 
-def rescaling_target_distribution():
+def rescaling_target_distribution() -> None:
     xl = -10
     xr = 4
     # f is the density on x

src/doptics/solver.py

@@ -47,7 +47,7 @@ def solve_single_mirror_parallel_source(starting_distribution: Callable, target_
         # plt.savefig(f'solution-mirror-{FILENAME[result_type]}.png')
 
 
-def f(x:float, u:float) -> float:
+def f(x: float, u: float) -> float:
     return u**2
 
 

src/doptics/symbolic.py

@@ -19,4 +19,3 @@ def u_prime_mv():
     dudx = solve(Eq(diff(sol, x), 0), diff(u))[0].simplify()
 
     return lambdify([x, y, u, V, t1, t2, L1], dudx), lambdify([x, y, u, V, t1, t2, L1], sol)
-

src/doptics/two_mirrors.py

@@ -37,7 +37,7 @@ def solve_two_mirrors_parallel_source_two_targets(starting_density: Callable, ta
     """
     x_discrete = np.linspace(x_span[0], x_span[1], number_rays)
     ax = sp.integrate.quad(starting_density, x_span[0], x_span[1])[0]
-    starting_density_rescaled = lambda x: 1 / ax * starting_density(x)
+    def starting_density_rescaled(x): return 1 / ax * starting_density(x)
 
     res = []
 
@@ -156,7 +156,7 @@ def solve_two_mirrors_parallel_source_point_target(starting_density: Callable, t
         small_angle=-0.7 * np.pi,
         large_angle=-0.3 * np.pi)
     # new stuff ends
-    print(f'back in main program')
+    print('back in main program')
     print(G1(y1_span[1]))
     x_discrete = np.linspace(x_span[0], x_span[1], number_rays)
 
@@ -283,4 +283,3 @@ def solve_two_mirrors_parallel_source_point_target(starting_density: Callable, t
         plt.show()
 
     return None
-

src/ray/ray.py

@@ -4,4 +4,4 @@
 class Ray:
     def __init__(self, x: NDArray[float], y: NDArray[float]):
         self.x = x
-        self.y = y
+        self.y = y

src/ray/raytracing/raytracing.py

@@ -4,6 +4,7 @@
 import matplotlib.pyplot as plt
 import scipy as sp
 import matplotlib as mpl
+from functools import lru_cache
 mpl.rcParams['figure.dpi'] = 100
 COLORS = {'szegedblue': '#3f3fa6'}
 
@@ -27,44 +28,43 @@ def intersection(l1, l2):
         return np.nan
 
 
-
-
 """
 Example of ray tracing
 Variable names are same as in the report
 """
 
 L1 = 8
 L2 = 14
-y1_span = (10.5, 13.5)
-y2_span = (11, 13)
-x_span = (0, 3)
+y1_SPAN = (10.5, 13.5)
+y2_SPAN = (11, 13)
+x_SPAN = (0, 3)
 u0 = 4
 w0 = 6
 # E = lambda x: np.exp(-((x - 1.5) / 1)**2 / 2) / (2 * np.pi)**.5
-G1 = lambda y1: 1
-G2 = lambda y2: 1
+def G1(y1): return 1
+def G2(y2): return 1
+
+
 E = G2
-pm1 = -1
-pm2 = 1
+PM1 = -1
+PM2 = 1
+
+
+@lru_cache(maxsize=100)
+def integrate(f_, left: float, right: float) -> float:
+    return sp.integrate.quad(f_, left, right)[0]
 
 
 if __name__ == '__main__':
     # calculate theoretical target from the mapping between emitting density E and target density G1
-    G1_int = sp.integrate.quad(G1, y1_span[0], y1_span[1])[0]
-    G1_fixed = lambda y1: G1(y1) / G1_int
-    E_int = sp.integrate.quad(E, x_span[0], x_span[1])[0]
-    E_fixed = lambda x: E(x) / E_int
-    # a2 = G1_int / sp.integrate.quad(G2, y2_span[0], y2_span[1])[0]
-    # print(f'calculated a2 as {a2}')
-    # G2_fixed = lambda y2: a2 * G2(y2)
-    # a1 = G1_int / sp.integrate.quad(E, x_span[0], x_span[1])[0]
-    # E_fixed = lambda x: a1 * E(x)
-
-    m1_prime = lambda x, m: pm1 * E_fixed(x) / G1_fixed(m)
-    y1_0 = y1_span[1]
-    m1_solved = sp.integrate.solve_ivp(m1_prime, x_span, [y1_0], dense_output=1)
-    m1 = lambda x: m1_solved.sol(x)[0]
+    def G1_fixed(y1, y1_span=y1_SPAN): return G1(y1) / integrate(G1, y1_span[0], y1_span[1])
+    def E_fixed(x, x_span=x_SPAN): return E(x) / integrate(x, x_span[0], x_span[1])
+    def m1_prime(x, m, pm1=PM1): return pm1 * E_fixed(x) / G1_fixed(m)
+
+    y1_0 = y1_SPAN[1]
+    m1_solved = sp.integrate.solve_ivp(m1_prime, x_SPAN, [y1_0], dense_output=1)
+    def m1(x): return m1_solved.sol(x)[0]
+    # m1 = lambda x: m1_solved.sol(x)[0]
 
     # absue that solution was calculated already by importing mirror points
     reflectors = np.load("reflectors.npz")
@@ -79,7 +79,7 @@ def intersection(l1, l2):
     # spline_B = sp.interpolate.CubicSpline(np.flipud(B[:][0]), np.flipud(B[:][1]))
     # spline_B = sp.interpolate.CubicSpline(np.flipud(B[:][0]), B[:][1])
 
-    # x_arr_u = np.linspace(x_span[0], x_span[1], 1000)
+    # x_arr_u = np.linspace(x_SPAN[0], x_SPAN[1], 1000)
     # x_arr_w = np.linspace(B[0][0], B[-1][0], 1000)
     # plt.plot(x_arr_u, spline_A(x_arr_u), "r")
     # plt.plot(A[:, 0], A[:, 1], "k")
@@ -161,7 +161,6 @@ def intersection(l1, l2):
         # plt.plot([xis, xis, xis + 10], [np.zeros(len(xis)), hits_U, hits_U - 10 * U_reflection])
         # plt.show()
 
-
         # Plotting
         # for i in range(n_rays):
         #     plt.plot([Ps[i][0], As[i][0], Bs[i][0], Bs[i][0] + 4 * u2s[i][0]],
@@ -181,7 +180,7 @@ def intersection(l1, l2):
             ray_low = (Bs[i][0], Bs[i][1])
             ray_upp = (Bs[i][0] + 4 * u2s[i][0], Bs[i][1] + 4 * u2s[i][1])
             ray = line(ray_low, ray_upp)
-            target_line = line((y1_span[0], L1), (y1_span[1], L1))
+            target_line = line((y1_SPAN[0], L1), (y1_SPAN[1], L1))
             intersec = intersection(target_line, ray)
             # print(intersec)
             try:

test2.py

@@ -8,10 +8,6 @@
 import doptics.two_mirrors as tms
 from doptics.functions import cdf_sampling_source
 
-# Enable ctrl+c
-# import signal
-# signal.signal(signal.SIGINT, signal.SIG_DFL)
-
 if __name__ == '__main__':
     l1 = 5
     l2 = 10
@@ -35,20 +31,18 @@
                                                                number_rays=15
                                                                )
 
-    G1 = lambda y1: 1
-    G2 = lambda y2: 1 - np.abs(y2 - 12) if 11 < y2 < 13 else 0
-    E = lambda x: 1 / (np.exp(10 * (x - 0.5)) + np.exp(-10 * (x - 0.5)))
+    def G1(y1): return 1
+    def G2(y2): return 1 - np.abs(y2 - 12) if 11 < y2 < 13 else 0
+    def E(x): return 1 / (np.exp(10 * (x - 0.5)) + np.exp(-10 * (x - 0.5)))
 
     # new stuff
     # E = func.uniform
 
     angle_result = tms.solve_two_mirrors_parallel_source_point_target(starting_density=func.uniform,
-                                                                      target_distribution_1=lambda y1: 1,
-                                                                      target_distribution_2=lambda y2: 1 -
-                                                                                                       np.abs(y2 - 12)
-                                                                      if 11 < y2 < 13 else 0,
+                                                                      target_distribution_1=G1,
+                                                                      target_distribution_2=G2,
                                                                       x_span=x_span,
                                                                       y1_span=y1_span, y2_span=y2_span, u0=u0, w0=w0,
                                                                       l1=l1, l2=l2,
                                                                       # color='#a69f3f',
-                                                                      number_rays=16)
+                                                                      number_rays=16)