diff --git a/docs/source/tutorials/BasicCustomSimulator.ipynb b/docs/source/tutorials/BasicCustomSimulator.ipynb index 7c165220..ff4ae8eb 100644 --- a/docs/source/tutorials/BasicCustomSimulator.ipynb +++ b/docs/source/tutorials/BasicCustomSimulator.ipynb @@ -30,9 +30,9 @@ "### Part 1: The `__init__` function\n", "\n", "First, we begin by creating a new **class** for our simulator. \n", - "**For those new to object-oriented programming**: a class consists of a recipe for building an object. An object is a reusable container which contains functions and their parameters. \n", + "**For those new to object-oriented programming**: a class consists of a recipe for building an object. An object is a reusable container which contains functions and their parameters. Objects are also referred to as **instances** of a class, and the parameters which are assigned to a class are called **instance variables**. \n", "\n", - "We want our simulator to inherit from the **Module** class in Caustics, which is a basic framework for constructing simulator objects. In the `__init__` function, we need a few basic ingredients to create the simulator:\n", + "We want our simulator to inherit from the **Module** class in Caustics, which is a basic framework for constructing simulator objects. To create inheritance, we put the parent class as an argument (in parentheses) to the child class. In the `__init__` function, we need a few basic ingredients to create the simulator:\n", "1. A lens mass distribution\n", "2. A model for the lens light\n", "3. A model for the source light\n", @@ -46,12 +46,12 @@ "\n", "Within our `__init__` function, we need to provide instructions to construct the basic structure of the simulator object, which is done by calling the `__init__` function of the `super` class, which in this case is `Module` from Caustics.\n", "\n", - "Within `__init__` we also need to construct the components of our simulator. For components which are constructed once (lens mass model, lens light model, source light model, and telescope psf), we simply need to make them a part of the current object (`self`). We do the same for parameters whose value we wish to only set once, such as the coordinate grid, which we generate with the `meshgrid` function of caustics. For parameters which we wish to sample with our MCMC (which are not already parameters of any of the existing components), we need to register them as a `Param` object, which will allow our simulator to find them in the flattened vector of parameters which we will pass to the simulator. In this example, we register the source redshift `z_s` as a `Param`. " + "Within `__init__` we also need to construct the components of our simulator. For components which are constructed once (lens mass model, lens light model, source light model, and telescope psf), we simply need to make them instance variables of the current object being constructed (`self`). We do the same for parameters whose value we wish to only set once, such as the coordinate grid, which we generate with the `meshgrid` function of caustics. For parameters which we wish to sample with our MCMC (which are not already parameters of any of the existing components), we need to register them as a `Param` object, which will allow our simulator to find them in the flattened vector of parameters which we will pass to the simulator. In this example, we register the source redshift `z_s` as a `Param`. For more information on the underlying functionality of `Module`, `Param`, and related parameter handling capabilities in Caustics, see the underlying **caskade** package and associated documentation: https://caskade.readthedocs.io/en/latest/notebooks/BeginnersGuide.html" ] }, { "cell_type": "code", - "execution_count": 20, + "execution_count": 48, "id": "7e022279-469a-45e6-8920-c2f40ac88466", "metadata": {}, "outputs": [], @@ -128,14 +128,16 @@ "4. Downsample the image to the correct pixel scale\n", "5. Convolve with the PSF of the telescope\n", "\n", + "To ensure that all the Param parameters in the simulator are handled correctly, we need to add the `@forward` decorator from `Caustics` (which wraps the `@forward` decorator in `caskade`) to our `forward` function.\n", + "\n", "### Part 3: Instantiating our simulator\n", "\n", - "Now that we have completed our custom simulator, we need to **instantiate** the components of the simulator and the simulator itself. The instantiation process creates an object in memory from a class." + "Now that we have completed our custom simulator, we need to **instantiate** the components of the simulator and the simulator itself. The instantiation process creates an object in memory from a class. " ] }, { "cell_type": "code", - "execution_count": 21, + "execution_count": 49, "id": "4a50ae4d-3964-4dfa-803d-3dec4e506c00", "metadata": {}, "outputs": [], @@ -175,12 +177,12 @@ "id": "87ec40c1-3860-4f3a-a053-ba1a63d5ec5f", "metadata": {}, "source": [ - "Now that we have instantiated our simulator, we can visualize its structure using graphviz. The grayed out squares are parameters which are fixed, while the white squares are parameters which are registered as `Param` objects (known as **active parameters** in Caustics). The arrows indicate which object contains which component. " + "Now that we have instantiated our simulator, we can visualize its structure using graphviz. The grayed out squares are parameters which are fixed (known as **static parameters** in Caustics), while the white squares are parameters whose value will be set once the `forward` function is run (these are known as **dynamic parameters** in Caustics). The arrows indicate which object contains which component. " ] }, { "cell_type": "code", - "execution_count": 22, + "execution_count": 50, "id": "1180567a-c4f2-436d-821e-2d673bee0675", "metadata": {}, "outputs": [ @@ -198,357 +200,357 @@ "\n", "%3\n", "\n", - "\n", + "\n", "\n", - "134556280499904\n", + "140601116849776\n", "\n", - "Singlelens(sim_2)\n", + "Singlelens(sim_1)\n", "\n", - "\n", + "\n", "\n", - "134556276696032\n", + "140601092836656\n", "\n", "EPL(epl_2)\n", "\n", - "\n", + "\n", "\n", - "134556280499904->134556276696032\n", + "140601116849776->140601092836656\n", "\n", "\n", "\n", - "\n", + "\n", "\n", - "134556276702128\n", + "140601105280704\n", "\n", "Sersic(sourcelight_2)\n", "\n", - "\n", + "\n", "\n", - "134556280499904->134556276702128\n", + "140601116849776->140601105280704\n", "\n", "\n", "\n", - "\n", + "\n", "\n", - "134556276521472\n", + "140601105284928\n", "\n", "Sersic(lenslight1_2)\n", "\n", - "\n", + "\n", "\n", - "134556280499904->134556276521472\n", + "140601116849776->140601105284928\n", "\n", "\n", "\n", - "\n", + "\n", "\n", - "134556278498896\n", + "140601116845168\n", "\n", "z_s\n", "\n", - "\n", + "\n", "\n", - "134556280499904->134556278498896\n", + "140601116849776->140601116845168\n", "\n", "\n", "\n", - "\n", + "\n", "\n", - "134556280492272\n", + "140601105272784\n", "\n", "FlatLambdaCDM(cosmo_2)\n", "\n", - "\n", + "\n", "\n", - "134556276696032->134556280492272\n", + "140601092836656->140601105272784\n", "\n", "\n", "\n", - "\n", + "\n", "\n", - "134556280490928\n", + "140601116845216\n", "\n", "z_l\n", "\n", - "\n", + "\n", "\n", - "134556276696032->134556280490928\n", + "140601092836656->140601116845216\n", "\n", "\n", "\n", - "\n", + "\n", "\n", - "134556276696368\n", + "140601105272688\n", "\n", "x0\n", "\n", - "\n", + "\n", "\n", - "134556276696032->134556276696368\n", + "140601092836656->140601105272688\n", "\n", "\n", "\n", - "\n", + "\n", "\n", - "134556280501008\n", + "140601116851600\n", "\n", "y0\n", "\n", - "\n", + "\n", "\n", - "134556276696032->134556280501008\n", + "140601092836656->140601116851600\n", "\n", "\n", "\n", - "\n", + "\n", "\n", - "134556280494000\n", + "140601116851216\n", "\n", "q\n", "\n", - "\n", + "\n", "\n", - "134556276696032->134556280494000\n", + "140601092836656->140601116851216\n", "\n", "\n", "\n", - "\n", + "\n", "\n", - "134556280500816\n", + "140601116851264\n", "\n", "phi\n", "\n", - "\n", + "\n", "\n", - "134556276696032->134556280500816\n", + "140601092836656->140601116851264\n", "\n", "\n", "\n", - "\n", + "\n", "\n", - "134556280499856\n", + "140601116851360\n", "\n", "b\n", "\n", - "\n", + "\n", "\n", - "134556276696032->134556280499856\n", + "140601092836656->140601116851360\n", "\n", "\n", "\n", - "\n", + "\n", "\n", - "134556280491552\n", + "140601116850784\n", "\n", "t\n", "\n", - "\n", + "\n", "\n", - "134556276696032->134556280491552\n", + "140601092836656->140601116850784\n", "\n", "\n", "\n", - "\n", + "\n", "\n", - "134556280491072\n", + "140601105268992\n", "\n", "h0\n", "\n", - "\n", + "\n", "\n", - "134556280492272->134556280491072\n", + "140601105272784->140601105268992\n", "\n", "\n", "\n", - "\n", + "\n", "\n", - "134556276697376\n", + "140601093131952\n", "\n", "critical_density_0\n", "\n", - "\n", + "\n", "\n", - "134556280492272->134556276697376\n", + "140601105272784->140601093131952\n", "\n", "\n", "\n", - "\n", + "\n", "\n", - "134556276604256\n", + "140601105278016\n", "\n", "Om0\n", "\n", - "\n", + "\n", "\n", - "134556280492272->134556276604256\n", + "140601105272784->140601105278016\n", "\n", "\n", "\n", - "\n", + "\n", "\n", - "134556278900192\n", + "140601105280320\n", "\n", "x0\n", "\n", - "\n", + "\n", "\n", - "134556276702128->134556278900192\n", + "140601105280704->140601105280320\n", "\n", "\n", "\n", - "\n", + "\n", "\n", - "134556280492560\n", + "140601103245584\n", "\n", "y0\n", "\n", - "\n", + "\n", "\n", - "134556276702128->134556280492560\n", + "140601105280704->140601103245584\n", "\n", "\n", "\n", - "\n", + "\n", "\n", - "134556280495392\n", + "140601103240640\n", "\n", "q\n", "\n", - "\n", + "\n", "\n", - "134556276702128->134556280495392\n", + "140601105280704->140601103240640\n", "\n", "\n", "\n", - "\n", + "\n", "\n", - "134556280490352\n", + "140601103253456\n", "\n", "phi\n", "\n", - "\n", + "\n", "\n", - "134556276702128->134556280490352\n", + "140601105280704->140601103253456\n", "\n", "\n", "\n", - "\n", + "\n", "\n", - "134556280490448\n", + "140601103246928\n", "\n", "n\n", "\n", - "\n", + "\n", "\n", - "134556276702128->134556280490448\n", + "140601105280704->140601103246928\n", "\n", "\n", "\n", - "\n", + "\n", "\n", - "134556280500240\n", + "140601103239152\n", "\n", "Re\n", "\n", - "\n", + "\n", "\n", - "134556276702128->134556280500240\n", + "140601105280704->140601103239152\n", "\n", "\n", "\n", - "\n", + "\n", "\n", - "134556280491312\n", + "140601103239824\n", "\n", "Ie\n", "\n", - "\n", + "\n", "\n", - "134556276702128->134556280491312\n", + "140601105280704->140601103239824\n", "\n", "\n", "\n", - "\n", + "\n", "\n", - "134556280501056\n", + "140601116845264\n", "\n", "x0\n", "\n", - "\n", + "\n", "\n", - "134556276521472->134556280501056\n", + "140601105284928->140601116845264\n", "\n", "\n", "\n", - "\n", + "\n", "\n", - "134556280490304\n", + "140601116845744\n", "\n", "y0\n", "\n", - "\n", + "\n", "\n", - "134556276521472->134556280490304\n", + "140601105284928->140601116845744\n", "\n", "\n", "\n", - "\n", + "\n", "\n", - "134556280491504\n", + "140601116848048\n", "\n", "q\n", "\n", - "\n", + "\n", "\n", - "134556276521472->134556280491504\n", + "140601105284928->140601116848048\n", "\n", "\n", "\n", - "\n", + "\n", "\n", - "134556280501680\n", + "140601116849680\n", "\n", "phi\n", "\n", - "\n", + "\n", "\n", - "134556276521472->134556280501680\n", + "140601105284928->140601116849680\n", "\n", "\n", "\n", - "\n", + "\n", "\n", - "134556280490112\n", + "140601116844640\n", "\n", "n\n", "\n", - "\n", + "\n", "\n", - "134556276521472->134556280490112\n", + "140601105284928->140601116844640\n", "\n", "\n", "\n", - "\n", + "\n", "\n", - "134556280493136\n", + "140601116845840\n", "\n", "Re\n", "\n", - "\n", + "\n", "\n", - "134556276521472->134556280493136\n", + "140601105284928->140601116845840\n", "\n", "\n", "\n", - "\n", + "\n", "\n", - "134556280492896\n", + "140601116846176\n", "\n", "Ie\n", "\n", - "\n", + "\n", "\n", - "134556276521472->134556280492896\n", + "140601105284928->140601116846176\n", "\n", "\n", "\n", @@ -556,10 +558,10 @@ "\n" ], "text/plain": [ - "" + "" ] }, - "execution_count": 22, + "execution_count": 50, "metadata": {}, "output_type": "execute_result" } @@ -575,12 +577,12 @@ "source": [ "### Part 4: Passing parameters to our simulator\n", "\n", - "Now that we have " + "Now that we have designed our simulator class and instantiated our simulator object, we can use the `forward` method to run the simulator. Thanks to `caskade`, we can pass all of the dynamic parameters at once as a flattened Pytorch tensor. However, we need to know what order to put our parameters in the tensor. We can find the order by literally printing our simulator:" ] }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 51, "id": "379c21bd-a2cc-4a13-a882-a27bf8646ff2", "metadata": {}, "outputs": [ @@ -588,9 +590,9 @@ "name": "stdout", "output_type": "stream", "text": [ - "sim_0|module\n", - " epl_0|module\n", - " cosmo_0|module\n", + "sim_1|module\n", + " epl_2|module\n", + " cosmo_2|module\n", " h0|static\n", " critical_density_0|static\n", " Om0|static\n", @@ -601,7 +603,7 @@ " phi|dynamic\n", " b|dynamic\n", " t|dynamic\n", - " sourcelight_0|module\n", + " sourcelight_2|module\n", " x0|dynamic\n", " y0|dynamic\n", " q|dynamic\n", @@ -609,7 +611,7 @@ " n|dynamic\n", " Re|dynamic\n", " Ie|dynamic\n", - " lenslight1_0|module\n", + " lenslight1_2|module\n", " x0|dynamic\n", " y0|dynamic\n", " q|dynamic\n", @@ -627,26 +629,99 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 56, "id": "6e901bae-fb96-4eaa-8eea-a6263183e755", "metadata": {}, "outputs": [], "source": [ - "params_lens_mass = torch.tensor([\n", - " 0.25, # epl x0\n", - " 0.3, # epl y0\n", - " 1/1.14, # epl q\n", - " np.pi/2 + 1.6755160819145565, # epl phi\n", - " 0.97, # epl b\n", - " 1.04, # epl t\n", - " x0=0.25,y0=0.3,q=1-0.29,phi=-28/180*torch.pi,n=4,Re=0.84/6.646,Ie=36\n", + "params_for_simulator = torch.tensor([\n", + " #Lens redshift and mass\n", + " 1.5, #z_l\n", + " 0.25, # epl x0\n", + " 0.3, # epl y0\n", + " 1/1.14, # epl q\n", + " np.pi/2 + 1.6755160819145565, # epl phi\n", + " 0.97, # epl b\n", + " 1.04, # epl t\n", + " #Source light\n", + " 0.25, #x0\n", + " 0.3, #y0\n", + " 1-0.29, #q\n", + " -30/180*torch.pi, #phi\n", + " 4, #n\n", + " 0.1, #Re\n", + " 36, #Ie\n", + " #Lens light\n", + " 0.25, #x0\n", + " 0.3, #y0\n", + " 1-0.29, #q\n", + " -30/180*torch.pi, #phi\n", + " 4, #n\n", + " 0.1, #Re\n", + " 100, #Ie\n", + " #Source redshift\n", + " 3.5 #z_s\n", " ])" ] }, { "cell_type": "code", - "execution_count": null, - "id": "ac5c918c-2193-4b89-879d-885033a126c9", + "execution_count": 57, + "id": "adbfbd05-163d-440c-90d6-e079ad2f3710", + "metadata": {}, + "outputs": [], + "source": [ + "lensed_image = simulator.forward(params_for_simulator)" + ] + }, + { + "cell_type": "markdown", + "id": "87304927-ee07-4943-bcfc-94b9aadab936", + "metadata": {}, + "source": [ + "We can then view the lensed image output by our simulator (here we have created an \"Einstein cross\"):" + ] + }, + { + "cell_type": "code", + "execution_count": 58, + "id": "32756499-ea15-411c-ac14-4b5b776567d1", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAakAAAGhCAYAAADbf0s2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA1fklEQVR4nO3de3Bc5X3H/8/Zi1YXS/INSxa2wfyq/riYJMSmnhoaOwWcaYGEYZoQLgEm/QNqIChuY3CdNoYpErit62mcmDHTIbTEhekEEtqhrZUmMWHcFsfgBEx/kExccABFXGRdLGlX2n1+f1D2+R5Za8v2Cj1av18zmjk+e3Z1zsr2V+ez3+d5IuecEwAAAUpM9QkAAFAKRQoAECyKFAAgWBQpAECwKFIAgGBRpAAAwaJIAQCCRZECAASLIgUACBZFCgAQrCktUt/85je1ePFiVVdXa+nSpfrxj388lacDAAjMlBWpxx9/XG1tbdqwYYNeeOEF/c7v/I5+7/d+T6+//vpUnRIAIDDRVE0wu3z5cn384x/Xtm3bivvOOeccXXXVVero6DjqcwuFgt58803V19criqLJPlUAQJk559Tf36+WlhYlEqXvl1If4jkV5XI57d27V3fffXds/+rVq7V79+4jjs9ms8pms8U/v/HGGzr33HMn/TwBAJPr4MGDWrBgQcnHp6RIvfPOO8rn82pqaortb2pqUldX1xHHd3R06J577jli/8X6faWUnrTzBABMjlGN6Fk9rfr6+qMeNyVF6gNjozrn3Ljx3fr167V27drin/v6+rRw4UKllFYqokgBwLTzfx80HesjmykpUnPnzlUymTzirqm7u/uIuytJymQyymQyH9bpAQACMSXdfVVVVVq6dKk6Oztj+zs7O7VixYqpOCUAQICmLO5bu3atvvCFL2jZsmX67d/+bW3fvl2vv/66br311qk6JQBAYKasSF1zzTV69913de+99+qtt97SkiVL9PTTT+uMM86YqlMCAARmysZJnYy+vj41NjZqlT5D4wQATEOjbkQ/0vfU29urhoaGkscxdx8AIFgUKQBAsChSAIBgUaQAAMGiSAEAgkWRAgAEiyIFAAgWRQoAECyKFAAgWBQpAECwKFIAgGBRpAAAwaJIAQCCRZECAASLIgUACBZFCgAQLIoUACBYFCkAQLAoUgCAYFGkAADBokgBAIJFkQIABIsiBQAIFkUKABAsihQAIFgUKQBAsChSAIBgUaQAAMGiSAEAgkWRAgAEiyIFAAgWRQoAECyKFAAgWBQpAECwKFIAgGBRpAAAwaJIAQCCRZECAASLIgUACBZFCgAQLIoUACBYFCkAQLAoUgCAYFGkAADBokgBAIJFkQIABIsiBQAIFkUKABAsihQAIFgUKQBAsChSAIBglb1IdXR06MILL1R9fb3mzZunq666Sq+88krsGOecNm7cqJaWFtXU1GjVqlXav39/uU8FADDNlb1I7dq1S7fddpv+67/+S52dnRodHdXq1at1+PDh4jGbNm3S5s2btXXrVu3Zs0fNzc267LLL1N/fX+7TAQBMY5Fzzk3mN3j77bc1b9487dq1S5/4xCfknFNLS4va2tp01113SZKy2ayampr0wAMP6JZbbjnma/b19amxsVGr9BmlovRknj4AYBKMuhH9SN9Tb2+vGhoaSh436Z9J9fb2SpJmz54tSTpw4IC6urq0evXq4jGZTEYrV67U7t27x32NbDarvr6+2BcAoPJNapFyzmnt2rW6+OKLtWTJEklSV1eXJKmpqSl2bFNTU/GxsTo6OtTY2Fj8Wrhw4WSeNgAgEJNapG6//Xb97Gc/0z/+4z8e8VgURbE/O+eO2PeB9evXq7e3t/h18ODBSTlfAEBYUpP1wnfccYeeeuopPfPMM1qwYEFxf3Nzs6T376jmz59f3N/d3X3E3dUHMpmMMpnMZJ0qACBQZb+Tcs7p9ttv1xNPPKEf/OAHWrx4cezxxYsXq7m5WZ2dncV9uVxOu3bt0ooVK8p9OgCAaazsd1K33XabduzYoe9973uqr68vfs7U2NiompoaRVGktrY2tbe3q7W1Va2trWpvb1dtba2uu+66cp8OAGAaK3uR2rZtmyRp1apVsf0PP/ywbr75ZknSunXrNDQ0pDVr1qinp0fLly/Xzp07VV9fX+7TAQBMY5M+TmoyME4KAKa3YMZJAQBwoihSAIBgUaQAAMGiSAEAgkWRAgAEiyIFAAgWRQoAECyKFAAgWBQpAECwKFIAgGBRpAAAwaJIAQCCRZECAASLIgUACBZFCgAQLIoUACBYFCkAQLAoUgCAYFGkAADBokgBAIJFkQIABIsiBQAIFkUKABAsihQAIFgUKQBAsChSAIBgUaQAAMGiSAEAgkWRAgAEiyIFAAgWRQoAECyKFAAgWBQpAECwKFIAgGBRpAAAwaJIAQCCRZECAASLIgUACBZFCgAQLIoUACBYFCkAQLAoUgCAYFGkAADBokgBAIJFkQIABIsiBQAIVmqqTwCAEUVTfQYnxrmpPgNUKO6kAADBokgBAIJF3Ad8WI43yosC+R3SFY59zESujUgQJ2DS/xV0dHQoiiK1tbUV9znntHHjRrW0tKimpkarVq3S/v37J/tUAADTzKQWqT179mj79u36yEc+Etu/adMmbd68WVu3btWePXvU3Nysyy67TP39/ZN5OgCAaWbSitTAwICuv/56PfTQQ5o1a1Zxv3NOW7Zs0YYNG3T11VdryZIleuSRRzQ4OKgdO3ZM1ukAkyeKJviVGPcrSiaP/ZVOFb8SVemyfdnXnch5lLqG2NeJvE9ACZNWpG677TZdfvnluvTSS2P7Dxw4oK6uLq1evbq4L5PJaOXKldq9e/e4r5XNZtXX1xf7AgBUvklpnHjsscf0/PPPa8+ePUc81tXVJUlqamqK7W9qatJrr7027ut1dHTonnvuKf+JAgCCVvYidfDgQd15553auXOnqqurSx4XjbnFd84dse8D69ev19q1a4t/7uvr08KFC8tzwsBETSSWGhN3RYlo3MeipDkumfT70+afZCpl9qf9/iqzbZ7rUn47dq5juuqi0bz/w8ioPyyX8/tzI37/qD8mGjH786brz3QAukKJgOZoXYKl3ls6Ak95ZS9Se/fuVXd3t5YuXVrcl8/n9cwzz2jr1q165ZVXJL1/RzV//vziMd3d3UfcXX0gk8kok8mU+1QBAIEr+2dSl1xyiV588UXt27ev+LVs2TJdf/312rdvn8466yw1Nzers7Oz+JxcLqddu3ZpxYoV5T4dAMA0VvY7qfr6ei1ZsiS2r66uTnPmzCnub2trU3t7u1pbW9Xa2qr29nbV1tbquuuuK/fpACenVAxlozsb6Zn4TdL7HXEfbJv4TiYZiGx8l6kqbroaf0y+tsps+9cZrfHb+Yw/D5e0cV/81JNZH7ulBn30lxrwcV+if9hvD/ptN5z155312zYSVD4/7nbJGPD9B8fff5TYEqeGKZlxYt26dRoaGtKaNWvU09Oj5cuXa+fOnaqvr5+K0wEABCpybvr9etLX16fGxkat0meUitLHfgJwogK5kypM4Z1UVOJOSiXupFzJO6mj/FczkamXpt9/VTiKUTeiH+l76u3tVUNDQ8njApkcDACAIzHBLCCd1B1T7G6pasydvb1jqvbbzmznZ9htf8c0MsP/88w2+O+RbfTnkTPbo3X+TqPgX+aIO6nUoH+tql5/vtXv+SdVv1vjL6HH3zEle4f89Qz4bU3ksyrT7j72zmlCbeu0qZ+SuJMCAASLIgUACBZxH05dE4n4YrGe2V9l8jQb6dXGZ1mJxXr1/rFRG+vVm1iv0X+P7Eyz7edo1vBpvhkhMddHa6fN8qsIzKkZLG4XXPw6fz0wo7j93ju+o3awy59T7Vs+Bqzt9u9B9Tt+f1WPPz7RZ2LAQRN5DvumC0UlYkBJUeQju1jjhZ3Bgzb1UxJ3UgCAYFGkAADBIu7DqaVExGdjPZUY22QjvqjWd785E/EV6uJx32i9j/tGGvxr5eqTZtt06zX47exMH13l5voIrHbe4eL2OfN+Xdy+cKZfReA3Mn5/MorHZD/P+jky/7PxrOL2/lo/l+ZAVW1xO5/xv8uOmnFcNdX+GjLVPuJLHfLb0YB5XxPmfTUdgNKYyW1jD5SYuJbo75TBnRQAIFgUKQBAsIj7UPmiEms6lRqca2M9syaa7dxzM0wcZiO9GfHBvLZzb6TOf7+RWjsI12ybtLBgX6pEI6Lt3Ms7f21584S04tFYbcJHa7OqfFdeQ53vxHtnpn8PsiPmRMz3cwm7lpU/PpPyx6RS/pwSCb/tEvELiv3JTrFkBgBHCaK/UxF3UgCAYFGkAADBIu5DZZpAxFeqc081JuKr87Feod539MW69up9HDZaF/+9b7TaRHl2lnJzWMKvyK6UmQ4vKkTmGB+tDeX8YNyfDvrzfq3Xj/g9re43ittVSTM4VtLQqD/fdw/76+vr99vKmWjOPDdvFsgeMV2JUcFEf+a9z5jttNlOjOmytN8jGhrWeE4q+sO0xZ0UACBYFCkAQLCI+1A5Sg3UnUDEFxucazv3Gv12LOKbYRce9K9fSI+Jscw5Jc141VTWBFym8yyyu82vkIWU7QA0g2vr/HkPmWv4Zc0c/zrxdRhjIpMEpkb896gaGf8YG01asRiw1saA5v12/qDUmG67hPlzyeivYCK+WII5gbn+6PSbtriTAgAEiyIFAAgWcR+mr1JLbegoc/GVivjq64rbBRPxjTT4iGq0zr+O7dpzSRNvxRvplM76yCmZ8zFTImf3++1o1ER/BROB2Us13y+f9r9nFswce/kq0503Ju6zA2ljr2uyRtOsF7u+2H7zK24sprTnZ87JRqSxfFDx+C8W/RVKxKIyuWNsfr8JRHlEf9MKd1IAgGBRpAAAwSLuw/RylIgv1tlVIuKzA3VLRnyNtovPv06+yg4Qtl17PjJKZuODSVODPv9LDvnBqIkhH1dFOTNINWtiLNPNVjKWsvPhJRPj7lcy/ruoi0WhJvpL+/2FKrOdsTGneT9Ml6HtaiwVCdroL6qL/9cTmfgvGjXxZ968BwX/XroS700kewyDfCsBd1IAgGBRpAAAwSLuw/QVxX/HikqtqFs9/lx8+Qbf3WcjvpxZQddGfLaDzXbkpQbN9kB8tGtywK9AGw2a7SG/7XL+OW7Ej/h1NurKm7bBwnFGV8l4e1+s8zHtrzWR9nP6Jc0KvM5sp2v9+5Sv8/tH6+zg5vFjwNjg5Kr4zy5fa85jNGO2/XVHZgmP2PtRqhtQY1otx0OnX/C4kwIABIsiBQAIFnEfwjeBOfkkKTJdbPFBuz7uKzTYLj6/30Z8dqCujfgSposvNeCjpKpeH9El+sxaG5Kiw/7PbtjPQ1fImon8RkzcZyO+CXSk2Xhr7PtRZGMySc7GpFlzrSXiUmVM/DZkuvCG/PuXyJrtGSYGrDXdgKa7z85FOPaxhIn+ohHz/cxSHTIxoErEopFdRZhOv2mLOykAQLAoUgCAYBH3YXqxUdXY7r4S8/IVzKDdkZkm4pvp//rb5SXs8hqpYdvF56OkdJ+J+HoH/fftPxw7Jzdk4r5YF5+Jro537rkSJvpcu6qtbYCLPb9QorPQRIe22y5hjk/l7dIjZiC1+e+mUBWP+2wXoB0wHI34jsMoVyr6G7/r77g7/RAk7qQAAMGiSAEAgkXchzBNZJXd9Ji/vqYLzZkVa/Omiy87y8dHuRlmaQu/O7aCrp2XLz3gY6Vkv+/UsxGfO+yjP0lyOTM418RSJaO5k+k8iyb2O2ep721jwNjKtxMY5Gp/WrGziC3zYeYZTMQHGMfm/jMDffM1ZpBvrspsm0HTpjuyVPQ3oU6/sX/nGNwbBO6kAADBokgBAIJFkQIABIvPpBC+Em3nUVU6fpiZWcJOHpud5T/LyDaa5cxN23nCtE0nBs3MEkNmmXc7WWy//+wp3mZuPtBSOJ9DHTd7TuYzHDdm9ooi+1nhkP+8KWE+N0zG1quKn3fsMymzbdesSuRMO/qwmQVj2LznqfE/q+LTpemLOykAQLAoUgCAYBH3IRxHWxr+g0PsUuim5VyKt53nZvrHhmf7mCk708RK5m9/2kwUkRg1y8EPmhkWzHpQdrLY2EwS+Xh0d1IR30lEeSUnmz0RJaK/2OwOtiXctJdHQ6aFPOPjukTNmOXjzaS++bSJF228O+pfNznsX8sNm1kt7MS9dt2s2M+F2SemE+6kAADBokgBAIJF3Icw2ZjHdo7ZSWRrqmNPGZnp477huT4OGjrNPz/XYL6FSX1Sw6bTzzSwJYdN3GeWfLfrQdnJYl3+KFHSSUR8E4rvytnpN5G1rMyMDJG97hG7VpadJcJMSDsy5vWd6fyzk1GYdaci80Bq2MSIg+bvhO34zPrjnfmhss7U9MKdFAAgWBQpAECwiPsQPhtjmTjHdvNJUna27+g73OSfM9TkY6nRGT6WSg2YiKnXxn0mxjIRlbOdY/nYQkzHuoIjHW+sZ4+38aftiExM8HdOs/aTKzWJamECr2U7FyM7Oa15P+zkryN2Lar4e2aebpM/5avsOl9+fzJrBgkfTptt81+a6e6z79OEB/ba95bJZqcMd1IAgGBRpAAAwSLuw9Q6zgG8kRnAO9IQ7+4bmuvjnaH5Jp5Z4OfWq834+Gkw4ZeVL3Sbfwo2xbKx1OgEu/hKmUjEZ7sa7cBlG13ZQap2u1RUOHZAsT33EtvOnuoElqUv2elnl5u3cV++dESajy0lP/7g6+SIv+60ifsS/WZgb8o8YaJRKIIzKT+5N954QzfccIPmzJmj2tpafexjH9PevXuLjzvntHHjRrW0tKimpkarVq3S/v37J+NUAADTWNmLVE9Pjy666CKl02n967/+q15++WX99V//tWbOnFk8ZtOmTdq8ebO2bt2qPXv2qLm5WZdddpn6+/vLfToAgGms7HHfAw88oIULF+rhhx8u7jvzzDOL2845bdmyRRs2bNDVV18tSXrkkUfU1NSkHTt26JZbbin3KWEaikVgNqrJ+DhnpCG+VMfQXDP/2+k+4jv/9LeK26mEj6J+NtJS3C5UzTDf3LyoiaViy1SYCKzk/HxHUTLiM0tb2GgzqjbzFJpBzM7Mhxd7n0bN0ul2Xj1JGvaDkmWXtzfHRXaZi+OM/mz3YCxONPujMa9jv0fBXNJonT3IDuz1u9MD/j1LVZvoz3aCDo7/fjOPX/jKfif11FNPadmyZfrsZz+refPm6YILLtBDDz1UfPzAgQPq6urS6tWri/symYxWrlyp3bt3j/ua2WxWfX19sS8AQOUre5H65S9/qW3btqm1tVX//u//rltvvVVf+tKX9Pd///eSpK6uLklSU1NT7HlNTU3Fx8bq6OhQY2Nj8WvhwoXlPm0AQIDKHvcVCgUtW7ZM7e3tkqQLLrhA+/fv17Zt23TjjTcWj4vGdHU5547Y94H169dr7dq1xT/39fVRqCpRqbnnTAebqzZx34xk7LDcLB8hnT7vUHH7E3N+7p9jRor+qn9mcftQxsd9sXjLDuIsnNzcbhOK+EyUF80w3Ycz/fmNzPKDmEfr7Px0ZkByzp9rVW98teBkj19VOBrw2/EI00Rztuuv1FBYG98lo/H3HyUqLJg5+kZr/PZIvTkn++M2/1dUmUHZVT1mrsA03X2VoOw/ufnz5+vcc8+N7TvnnHP0+uuvS5Kam5sl6Yi7pu7u7iPurj6QyWTU0NAQ+wIAVL6yF6mLLrpIr7zySmzfq6++qjPOOEOStHjxYjU3N6uzs7P4eC6X065du7RixYpynw4AYBore9z35S9/WStWrFB7e7s+97nP6bnnntP27du1fft2Se/HfG1tbWpvb1dra6taW1vV3t6u2tpaXXfddeU+HUxXNg6z866lfbSTz8Tj4ZEZPhpaVN9T3F5SfbC43V/wUVl9xne59STtCN4S51Ri/jYb443t9Dvugbo1/vxsxDd0ut8emO//2WZnme9t4rC0Gc1R9+t4LFprzillu++ypuvPxmP2mESJKM+IzQdonhuZ/S4V//3Y/ixHTFCSPc1877SPIAsp37mXHvDPrX7HdPrZGNW838zCN72UvUhdeOGFevLJJ7V+/Xrde++9Wrx4sbZs2aLrr7++eMy6des0NDSkNWvWqKenR8uXL9fOnTtVX19f7tMBAExjkzIt0hVXXKErrrii5ONRFGnjxo3auHHjZHx7AECFYO4+fPgmMF9fTCwyM3O5Jce8TtLHQRmzEmtVNP6AzVGzHIUdHBo7/DiXaDjqCrolBihHZo45V2PmJpzto7+BFn9M31n+ZUbmmc49E1km3zFz2CXicV9qyD+W7DMDhk92aYsP2BjQvn/2Zzcm7hup849lZ/sfRsN8n1vWV/s48q3qxuL2cJ9/n0bq/fuUMQO/o3KuWowPFT85AECwKFIAgGAR92HaGjv/W5TzkdHbWd8N9/9l5xe33xv1+989XFvcTg6beeHsyrx27r4TOskSvwfayNPGfXZuwhnjd/HZiK/5dN/FmIz8Gb4ZzSpu594z8/5Jylebc0qN30VZKLlir43ySgxuLjko2+8vZMYMxG4w73+Tn3fxotMPFLdbMoeK2z9O/0Zx+xdvLyhu29jQme9x1BgWQeNOCgAQLIoUACBYxH2YXsxKualsPJJKD/jfuQ68N7u4vbvq/yluH8r5iK+/x283HPavkzTz3p3sfH1WbG5KO1g21vXnN+3g3Nh8giWSq7yzg4rHX9bi/T+PPx+hKxHxldo/Ifb9M9ccixwljZghkgtP8xHm5bP2FbfPTPn9dg7GVxubi9ujNWZ5jiS/g1cCfooAgGBRpAAAwSLuw7QSjfhBuqnBeI5V1esjoP5f+y6+F2S6v0b8Mcl3zfxvh32kZZe5KDUYNT6H3Un+rmfit8jEmckhu9yGPyb3tj/vLue7+OzKtem3/T/t6nfHxKJ9Zj69YbNqr115OF9ixdpSHX0TYOO3fFX8PRutNsus1PUWtz9a9U5xe0HK/0wXZ94ubier/XkXknalYjr6KgF3UgCAYFGkAADBokgBAILFZ1IIn/3MJvaZVPxzk0yPmaHhbf/Z04Dzn2VEebP20CH/O1rKfP5jPxeKKecS5LY1u2A/I/KzSVT1+AlVZ9TYyXD9teXMRLK21TxzyL9ndV3msyZJVe/5GR2iweHitjPvbfz8SrSml9gfJWzbufkvxn6mN/bjolhHvlk3yhySdf7zs0N5P3ygkPPvh5lXWFF+YkvXl3QyrfcoG+6kAADBokgBAIJF3IfgObsEuY37Do/Ejsv0+fbj6nfN5KIuPpnpB9IDfjuZG39SWSs2Y0RsvaWJxUJ25oaoxBLrGvLxW/KQP+9ac3i630d8sZkbzGmnB8371OPjPUlK9PrpNdyQf8yZFnRn3wPTdl4q4jteYycHTpgf5VuDfq2oPcMtxe36hD/Xn/Se6Z/Qa4YSDJr3OGte1L7HJ9FGjw8fd1IAgGBRpAAAwSLuQ5hciVkfRnyEEw3F4750v490MqZzz5mZB+xErakhM8tEvAHOsx19Znn12H53Al1kdmJX21UX+Y4+Gwkms77rL9Frloa3k6jaLsiceW+G/WtKkjN/dlmzbc/jZCK+qMR7Zg8Z834nB/3P6FfvzSxuP133kXGfv/eNhcXtml/775fpNX8/zGwabhJm0MCHgzspAECwKFIAgGAR9yF8NhozHWiRicAkKT3gH6vq9zFTIWUnNvWxUsI83UZrdiLU2JpEdjs6gclLTWxmY8fIRFEu508qdq3Dw+YJJX63LBHRHRF12e9Xpi6+2PLsZjsq8f6N7e5Lmcvrf7emuP1calFxO5fz/13lD/rBvDO7/WtVHbJxsIky7TWXqUMRHw7upAAAwaJIAQCCRdyH4MUGwdoILBfv7ksM+j+nB0w0ZCK+yCRfpTv6THRlutOiVGrc/bH4bIIDe2PRn32OvVY7YDh+qcfkjtZxWKKjrWQMZo8vFTUasUHPad+J6Ox7NiYttdFrygxiPpz368onBv33rnvTv0Dt2/4Hmew1cxHaLkbm4Zu2uJMCAASLIgUACBZxHz58pZZkj3WXmU4wkxLZbrSx3X0Js8xF6rCZz83ObzeB38sKadPdl/b/RKK0WZo8ab53qYG9Ujxqm8jAUXPdE44Oj+f1NcFY73jZGNC8T1GV3y5kbMdlPO9LjPpzSg/4x1LD/jnpfn98bZc/1+p3zEDngcHith2oPKFB1kSCQeJOCgAQLIoUACBYxH0In41qbCfdaLw9z65qmxzyXWWpwz4yckkfJdnIyUUl9lf7fyJu2MR9ORNp2cGyhXhkFovsbCRWKlo7zsjthAamHm+sV6KjL7KdjzYWNRGfq8kUtwsZf8zYuM++TSmzsoid48+uNlzzrlmK5JBZbsQsdRL7++HGH7SM8HEnBQAIFkUKABAs4j6Ez8YzNk4bE/fJDO5NmGU8UtWmw8x07qnGTqCncY+xEVUi4yNEu+prbG68MR1ikexxx46Zyjav3MlGWqUivhLz8tmIL6rxc+8Van3cN1pr4r50PO6LzOmmzOq6Nvqrede/z1U9ZkmTwybuyx7nfH109AWPOykAQLAoUgCAYBH3YVqJzcE2Ujrus8s0JDI+ikpV2WU4TEefiZ+c6TzLV5sOtjof9yVHfZQUFUqsIizJmcdswOXysYP8MXYV4eON/iYr4rNdfLGIz8SfJuJzdX47P8PHffkaO5g3/j3sPIoJs2xHVb+/pqpDdnViM2jXxn12PsdYDEtH33TFnRQAIFgUKQBAsIj7MLWOdx6/aPyBvZKkETOfnhlsmxj00V/Szh9n5+iLrcZrvoWZ9y8qmA42E8XZ3/Sio6zYa6OoeNefPegkor8JLKNx9KdH9g9+s8RA3VjEV19X3M43+v0j9f740erS700ya+buG/LbGbPSbuqQzwHp6Dt1cCcFAAgWRQoAECziPkwvNgZ0Y+Ij2+1nIiAbVyUG/XYyZQbtJu2SEqa7z3YDRuPP9ZcyMVkyGf+9L0qZ7NDMK6esOc6et42r7LIk5mXK2sVnB+fGVtQ1sZ4dxFxd7U/DdvHZiK/BHz9aO/68idGYpDY17CO4qj7/fqR7/HsWW4Zj2MzRlzMxLx19FYc7KQBAsChSAIBgEfchfBPp9JPiUY8d1Dlsoj8TvyXTJooynX4jSTPHnB+Lqny1HfBrugSrfASWysT/SSVrffSVGPTH2cHGkTlX26kWmbkJY/P+HWVpEP/N7LyEUcnHopQ5X7NtIz5X7bftXHx2oO7oDB8PjtTZDkoT8dn5+XLx8073+2tKmy6+RF+JQbvDtqPP/v2go6/ScCcFAAgWRQoAECyKFAAgWGX/TGp0dFQbN27Ut7/9bXV1dWn+/Pm6+eab9dWvflWJ/8vCnXO65557tH37dvX09Gj58uX6xje+ofPOO6/cp4PppNTsEyWPH7NUu21rtmtNZX2Lsv0MJmFa023reGw9KfO5VT7Wmu63R/1HTRqpi//elxo23y/rP8NJZv0MDYkRM1ntaGHcbfvZU5Q375OZ6NZ+VhXlS096a99b+9mazPtRqPLb+Zq02bb7/bWOVtsZO+zMIf57J3NmVomBeA96utdMCNxnPns6XKLtPPZ5HW3nlazsd1IPPPCAHnzwQW3dulX/8z//o02bNukv//Iv9fWvf714zKZNm7R582Zt3bpVe/bsUXNzsy677DL19/eX+3QAANNY2YvUf/7nf+ozn/mMLr/8cp155pn6gz/4A61evVo/+clPJL1/F7VlyxZt2LBBV199tZYsWaJHHnlEg4OD2rFjR7lPBwAwjZU97rv44ov14IMP6tVXX9Vv/uZv6qc//ameffZZbdmyRZJ04MABdXV1afXq1cXnZDIZrVy5Urt379Ytt9xS7lPCdHScE89KUpQwj8Xa0c2MBHY9JLOdNG3ZLjYThVlzyuy3M1HkfRp2ZEzp7AS1fnciNqmseXqpVmmz2669ZCO05HDB7C+Y48dEYDYJNNdnt+31FarGX3fLvjd2sSwbR9qZJFKHTZt5n4/3pPj6UFH/YX9OQycxswRt5xWh7EXqrrvuUm9vr84++2wlk0nl83ndd999uvbaayVJXV1dkqSmpqbY85qamvTaa6+N+5rZbFZZM36kr6+v3KcNAAhQ2eO+xx9/XI8++qh27Nih559/Xo888oj+6q/+So888kjsuLFLGjjnSi5z0NHRocbGxuLXwoULy33aAIAAlf1O6itf+Yruvvtuff7zn5cknX/++XrttdfU0dGhm266Sc3NzZJU7Pz7QHd39xF3Vx9Yv3691q5dW/xzX18fhQpHdvdNYN0pZzr9YmsmmV+QkmY7bddYSphJaO0S86brb7Qm/otW3sxYYZdML1TZY9y4+wsZm/GZFx01k9sO+gfSfb5TLz3gD08NxWOvpH0LzAwNrlRHZWzd+/Gfa2NH+/3SAz6bTNkOvn4f70lSZCePHbTLwY8f8bmxa4kVHyDiqzRlv5MaHBwstpp/IJlMqvB/LbSLFy9Wc3OzOjs7i4/ncjnt2rVLK1asGPc1M5mMGhoaYl8AgMpX9jupK6+8Uvfdd58WLVqk8847Ty+88II2b96sL37xi5Le/421ra1N7e3tam1tVWtrq9rb21VbW6vrrruu3KcDAJjGyl6kvv71r+vP/uzPtGbNGnV3d6ulpUW33HKL/vzP/7x4zLp16zQ0NKQ1a9YUB/Pu3LlT9fX15T4dVILjHeSr0nGQXbY9NqjVHGNzgFSJrMElbNefHfAbP86sOK+RBv9dRmaaTreZvoOtZZYfK3hmw7vF7dlVPg7Lmdzwl/1zituvvzeruN37Tq1//R67GJWUGjDRZtZsm2QtMeLPNWH2p8wy77Htwz56TfebJd/7bcQ3/iBdqXQXX2zQLpPHnpIi56bfT7ivr0+NjY1apc8oFaWP/QRUjqMVqWj8ihIlk2bbfA5VZT4AyvjqEpnF/AoN/j/70QY/tURuln/u8Ez/+rnG+PnlTDI9kSLVXKYiNTThIuX3H3eRGi5fkdLxFil3lBk1MC2MuhH9SN9Tb2/vUT/CYe4+AECwWE8K08vRoj/727W9q7IDgG0KaDvHSvw2njD702aQamwpK2fuyKL4XYsd/DpaZx5I+XOa3eDvKi6c68cKXta4v7j9/6b9XdWIabfb19hS3P6PGecWt/dULSpu96TiMbozXYrq969lBxhHZqV7e/eUPmy2+80y7wP+7inZ558cHS4xD5/tstRRuvi4ezrlcScFAAgWRQoAECziPkxfR1mCovQcf8fu+ou9qrNLZPjtVGH8/YnReHtflLfxn10CxEdu783wTQ5vNTYWt7vrfEy3MHXIvMr4UVfBzBOYt9c8Gv9dNNYI4dM4pfv961bFtu2ce6YpYsAsgTLgXygaMku7lxiYaxsi3j9hBupifNxJAQCCRZECAASLuA+VYwLLe9iuv5IrutrXMbGeTKwXmbgqlfURWGLYR3eSlBz2Y6tSwz7iS2R9DDiY821//zV8VnH79X4/7unMhveK22nThvfmYR8PvvauPz73a38eNd3xjsPMu2a7tzDudpUZ65Q0sV5iwEd50aAZ2zRsIj6zYkEs4svbLssxkR7LbaAE7qQAAMGiSAEAgkXch8o0kejPHm6X+VCJwaS2A83EfdGIifvGDFKtGvRTLKUG/HZVv48Bq3v8P8OhX/vuwHdnNxe3f13nl7FxCTtlkb+2qj6/XXvIH2NjvPf/PH63XvKwuY7DZkCu7dazsd6IifJypjtyxHTuxbosSwzMHYuIDwZ3UgCAYFGkAADBIu5D5ZvIUh8TGfzrxu/0k4n7bOwlSdGwj82SZh676l4f91W97WPAGQ1+HsCROv/PM19tuhLN6SVG/XkkY7OS+8gtORQ/p8SgGYRrojzZyC7WoWf221nJbYfeRObbs4j0MEHcSQEAgkWRAgAEi7gPp5bjjP7iu48zBlR8MGtkOuOiw4eL24lDPuJLZvx2VZUf/OvM6r/xbkWzfMioOSfbYTcSj/vcyPideK5E96JKxHcT7tYb51yBieJOCgAQLIoUACBYxH04dZWKn443BrTT0EVjXtM+ZuM00z2nyKyXYSXGP4+oxPkVSl1P4SgxW6n4rsQxE0KshzLiTgoAECyKFAAgWMR9wFgTiauOMxKUxsaCx14tuNTukw7Tjje+iz2XKA8fLu6kAADBokgBAIJFkQIABIvPpIATMdHPZk7gs6spw+dNCBB3UgCAYFGkAADBIu4DJhMRGnBSuJMCAASLIgUACBZFCgAQLIoUACBYFCkAQLAoUgCAYFGkAADBokgBAIJFkQIABIsiBQAIFkUKABAsihQAIFgUKQBAsChSAIBgUaQAAMGiSAEAgkWRAgAEiyIFAAgWRQoAECyKFAAgWBQpAECwjrtIPfPMM7ryyivV0tKiKIr03e9+N/a4c04bN25US0uLampqtGrVKu3fvz92TDab1R133KG5c+eqrq5On/70p/WrX/3qpC4EAFB5jrtIHT58WB/96Ee1devWcR/ftGmTNm/erK1bt2rPnj1qbm7WZZddpv7+/uIxbW1tevLJJ/XYY4/p2Wef1cDAgK644grl8/kTvxIAQMWJnHPuhJ8cRXryySd11VVXSXr/LqqlpUVtbW266667JL1/19TU1KQHHnhAt9xyi3p7e3XaaafpH/7hH3TNNddIkt58800tXLhQTz/9tD71qU8d8/v29fWpsbFRq/QZpaL0iZ4+AGCKjLoR/UjfU29vrxoaGkoeV9bPpA4cOKCuri6tXr26uC+TyWjlypXavXu3JGnv3r0aGRmJHdPS0qIlS5YUjxkrm82qr68v9gUAqHxlLVJdXV2SpKamptj+pqam4mNdXV2qqqrSrFmzSh4zVkdHhxobG4tfCxcuLOdpAwACNSndfVEUxf7snDti31hHO2b9+vXq7e0tfh08eLBs5woACFdZi1Rzc7MkHXFH1N3dXby7am5uVi6XU09PT8ljxspkMmpoaIh9AQAqX1mL1OLFi9Xc3KzOzs7ivlwup127dmnFihWSpKVLlyqdTseOeeutt/TSSy8VjwEAQJJSx/uEgYEB/eIXvyj++cCBA9q3b59mz56tRYsWqa2tTe3t7WptbVVra6va29tVW1ur6667TpLU2NioP/zDP9Qf//Efa86cOZo9e7b+5E/+ROeff74uvfTS8l0ZAGDaO+4i9ZOf/ESf/OQni39eu3atJOmmm27St771La1bt05DQ0Nas2aNenp6tHz5cu3cuVP19fXF5/zN3/yNUqmUPve5z2loaEiXXHKJvvWtbymZTJbhkgAAleKkxklNFcZJAcD0NiXjpAAAKCeKFAAgWBQpAECwKFIAgGBRpAAAwaJIAQCCRZECAASLIgUACBZFCgAQLIoUACBYFCkAQLAoUgCAYFGkAADBokgBAIJFkQIABIsiBQAIFkUKABAsihQAIFgUKQBAsChSAIBgUaQAAMGiSAEAgkWRAgAEiyIFAAgWRQoAECyKFAAgWBQpAECwKFIAgGBRpAAAwaJIAQCCRZECAASLIgUACBZFCgAQLIoUACBYFCkAQLAoUgCAYFGkAADBokgBAIJFkQIABIsiBQAIFkUKABAsihQAIFgUKQBAsChSAIBgUaQAAMGiSAEAgkWRAgAEiyIFAAgWRQoAECyKFAAgWBQpAECwKFIAgGClpvoEToRzTpI0qhHJTfHJAACO26hGJPn/z0uZlkWqv79fkvSsnp7iMwEAnIz+/n41NjaWfDxyxypjASoUCnrzzTflnNOiRYt08OBBNTQ0TPVpfSj6+vq0cOHCU+qaJa77VLruU/GapVPvup1z6u/vV0tLixKJ0p88Tcs7qUQioQULFqivr0+S1NDQcEr8UK1T8ZolrvtUcipes3RqXffR7qA+QOMEACBYFCkAQLCmdZHKZDL62te+pkwmM9Wn8qE5Fa9Z4rpPpes+Fa9ZOnWv+1imZeMEAODUMK3vpAAAlY0iBQAIFkUKABAsihQAIFjTtkh985vf1OLFi1VdXa2lS5fqxz/+8VSfUll1dHTowgsvVH19vebNm6errrpKr7zySuwY55w2btyolpYW1dTUaNWqVdq/f/8UnXH5dXR0KIoitbW1FfdV6jW/8cYbuuGGGzRnzhzV1tbqYx/7mPbu3Vt8vNKue3R0VF/96le1ePFi1dTU6KyzztK9996rQqFQPKYSrvmZZ57RlVdeqZaWFkVRpO9+97uxxydyjdlsVnfccYfmzp2ruro6ffrTn9avfvWrD/Eqppibhh577DGXTqfdQw895F5++WV35513urq6Ovfaa69N9amVzac+9Sn38MMPu5deesnt27fPXX755W7RokVuYGCgeMz999/v6uvr3Xe+8x334osvumuuucbNnz/f9fX1TeGZl8dzzz3nzjzzTPeRj3zE3XnnncX9lXjN7733njvjjDPczTff7P77v//bHThwwH3/+993v/jFL4rHVNp1/8Vf/IWbM2eO+5d/+Rd34MAB90//9E9uxowZbsuWLcVjKuGan376abdhwwb3ne98x0lyTz75ZOzxiVzjrbfe6k4//XTX2dnpnn/+effJT37SffSjH3Wjo6Mf8tVMjWlZpH7rt37L3XrrrbF9Z599trv77run6IwmX3d3t5Pkdu3a5ZxzrlAouObmZnf//fcXjxkeHnaNjY3uwQcfnKrTLIv+/n7X2trqOjs73cqVK4tFqlKv+a677nIXX3xxyccr8bovv/xy98UvfjG27+qrr3Y33HCDc64yr3lskZrINR46dMil02n32GOPFY954403XCKRcP/2b//2oZ37VJp2cV8ul9PevXu1evXq2P7Vq1dr9+7dU3RWk6+3t1eSNHv2bEnSgQMH1NXVFXsfMpmMVq5cOe3fh9tuu02XX365Lr300tj+Sr3mp556SsuWLdNnP/tZzZs3TxdccIEeeuih4uOVeN0XX3yx/uM//kOvvvqqJOmnP/2pnn32Wf3+7/++pMq85rEmco179+7VyMhI7JiWlhYtWbKkYt6HY5l2E8y+8847yufzampqiu1vampSV1fXFJ3V5HLOae3atbr44ou1ZMkSSSpe63jvw2uvvfahn2O5PPbYY3r++ee1Z8+eIx6r1Gv+5S9/qW3btmnt2rX60z/9Uz333HP60pe+pEwmoxtvvLEir/uuu+5Sb2+vzj77bCWTSeXzed1333269tprJVXuz9qayDV2dXWpqqpKs2bNOuKYSv3/bqxpV6Q+EEVR7M/OuSP2VYrbb79dP/vZz/Tss88e8VglvQ8HDx7UnXfeqZ07d6q6urrkcZV0zdL7S88sW7ZM7e3tkqQLLrhA+/fv17Zt23TjjTcWj6uk63788cf16KOPaseOHTrvvPO0b98+tbW1qaWlRTfddFPxuEq65lJO5Bor8X0oZdrFfXPnzlUymTzit4ju7u4jfiOpBHfccYeeeuop/fCHP9SCBQuK+5ubmyWpot6HvXv3qru7W0uXLlUqlVIqldKuXbv0t3/7t0qlUsXrqqRrlqT58+fr3HPPje0755xz9Prrr0uqzJ/1V77yFd199936/Oc/r/PPP19f+MIX9OUvf1kdHR2SKvOax5rINTY3NyuXy6mnp6fkMZVu2hWpqqoqLV26VJ2dnbH9nZ2dWrFixRSdVfk553T77bfriSee0A9+8AMtXrw49vjixYvV3Nwcex9yuZx27do1bd+HSy65RC+++KL27dtX/Fq2bJmuv/567du3T2eddVbFXbMkXXTRRUcML3j11Vd1xhlnSKrMn/Xg4OARC90lk8liC3olXvNYE7nGpUuXKp1Ox45566239NJLL1XM+3BMU9aycRI+aEH/u7/7O/fyyy+7trY2V1dX5/73f/93qk+tbP7oj/7INTY2uh/96EfurbfeKn4NDg4Wj7n//vtdY2Oje+KJJ9yLL77orr322mnXonsstrvPucq85ueee86lUil33333uZ///Ofu29/+tqutrXWPPvpo8ZhKu+6bbrrJnX766cUW9CeeeMLNnTvXrVu3rnhMJVxzf3+/e+GFF9wLL7zgJLnNmze7F154oThcZiLXeOutt7oFCxa473//++755593v/u7v0sL+nTwjW98w51xxhmuqqrKffzjHy+2ZlcKSeN+Pfzww8VjCoWC+9rXvuaam5tdJpNxn/jEJ9yLL744dSc9CcYWqUq95n/+5392S5YscZlMxp199tlu+/btsccr7br7+vrcnXfe6RYtWuSqq6vdWWed5TZs2OCy2WzxmEq45h/+8Ifj/ju+6aabnHMTu8ahoSF3++23u9mzZ7uamhp3xRVXuNdff30KrmZqsFQHACBY0+4zKQDAqYMiBQAIFkUKABAsihQAIFgUKQBAsChSAIBgUaQAAMGiSAEAgkWRAgAEiyIFAAgWRQoAECyKFAAgWP8/58w3WoOyZ94AAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.imshow(lensed_image)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "5da7997e-c1c5-4e99-bc0e-3f66454c9fff", + "metadata": {}, + "source": [ + "### Part 5: Customizing your simulator\n", + "\n", + "So far, we have focused on re-creating the LensSource simulator provided by default in Caustics, but the real power of the Caustics package is reflected by its extensibility. \n", + "\n", + "Suppose we want to have a single background light source and a single lens mass distribution, but instead of a single lens light source, we want two lens light sources (this could be a modeling choice for merging lensed galaxies).\n", + "\n", + "We can implement this by creating a new simulator class, which we will call Doublelens. This class is identical to Singlelens, except for two things: we add an extra `lens_light` to our `__init__`, and in the `forward` we add the second `lens_light` to the image." + ] + }, + { + "cell_type": "code", + "execution_count": 66, + "id": "d2b5c2fc-063b-42a7-a640-c6e90bf458ca", "metadata": {}, "outputs": [], "source": [ @@ -655,7 +730,7 @@ " self,\n", " lens,\n", " lens_light1,\n", - " lens_light2,\n", + " lens_light2, #NEW!\n", " source,\n", " psf,\n", " pixelscale,\n", @@ -668,11 +743,11 @@ "\n", " self.lens = lens\n", " self.src = source\n", - " self.psf = torch.as_tensor(psf, dtype=torch.float32)\n", + " self.psf = psf\n", " self.lens_light1 = lens_light1\n", - " self.lens_light2 = lens_light2\n", - " self.upsample_factor = Param(upsample_factor)\n", - " self.z_s = Param(z_s)\n", + " self.lens_light2 = lens_light2 #NEW!\n", + " self.upsample_factor = upsample_factor\n", + " self.z_s = Param(\"z_s\", z_s)\n", " \n", " # Create the high-resolution grid\n", " thx, thy = caustics.utils.meshgrid(\n", @@ -691,83 +766,697 @@ " bx, by = self.lens.raytrace(self.thx, self.thy, self.z_s)\n", " \n", " # Evaluate the lensed source brightness at high resolution\n", - " mu_fine = self.src.brightness(bx, by)\n", + " image = self.src.brightness(bx, by)\n", " \n", - " # Add the lens light contributions\n", - " mu_fine += self.lens_light1.brightness(self.thx, self.thy)\n", - " mu_fine += self.lens_light2.brightness(self.thx, self.thy)\n", + " # Add the lens light\n", + " image += self.lens_light1.brightness(self.thx, self.thy)\n", + " image += self.lens_light2.brightness(self.thx, self.thy) #NEW!\n", " \n", " # Downsample to the desired resolution\n", - " mu = avg_pool2d(mu_fine[None, None], self.upsample_factor)[0, 0]\n", + " image_ds = avg_pool2d(image[None, None], self.upsample_factor)[0, 0]\n", " \n", " \n", " # Convolve with the PSF using conv2d\n", " psf_normalized = (self.psf.T / self.psf.sum())[None, None]\n", - " mu = conv2d(\n", - " mu[None, None], psf_normalized, padding=\"same\"\n", + " image_ds = conv2d(\n", + " image_ds[None, None], psf_normalized, padding=\"same\"\n", " ).squeeze(0).squeeze(0)\n", " \n", - " return mu" + " return image_ds" ] }, { "cell_type": "code", - "execution_count": null, - "id": "a92fb895-8765-4c1c-93ff-a341ae44fa65", + "execution_count": 67, + "id": "00e39f80-128f-44ff-a831-71348870a656", "metadata": {}, "outputs": [], "source": [ - "class Lens_only(caustics.Simulator):\n", - " def __init__(\n", - " self,\n", - " lens_light,\n", - " psf,\n", - " pixelscale,\n", - " pixels_x,\n", - " upsample_factor,\n", - " name: str = \"sim\",\n", - " ):\n", - " super().__init__(name)\n", - "\n", - " self.psf = torch.as_tensor(psf, dtype=torch.float32)\n", - " self.lens_light = lens_light\n", - " self.upsample_factor = upsample_factor\n", - "\n", - " \n", - " # Create the high-resolution grid\n", - " thx, thy = caustics.utils.meshgrid(\n", - " pixelscale / upsample_factor,\n", - " upsample_factor * pixels_x,\n", - " dtype=torch.float32, device=\"cuda\"\n", - " )\n", - " #thx=thx.requires_grad_()\n", - " #thy=thy.requires_grad_()\n", - " \n", - " # CHANGE THIS TO self.thx = thx\n", - " #self.thx = thx\n", - " #self.thy = thy\n", - " self.add_param(\"thx\", thx)\n", - " self.add_param(\"thy\", thy)\n", - " \n", - " @forward\n", - " def forward(self, params):\n", - " # Unpack the parameters\n", - " thx, thy = self.unpack(params)\n", - " \n", - " # Add the lens light contributions\n", - " mu_fine = self.lens_light.brightness(thx, thy, params)\n", - " \n", - " # Downsample to the desired resolution\n", - " mu = avg_pool2d(mu_fine[None, None], self.upsample_factor).squeeze(0).squeeze(0)\n", - " \n", - " # Convolve with the PSF using conv2d\n", - " psf_normalized = ((torch.flip(self.psf, (0, 1))) / self.psf.sum())[None, None]\n", - " mu = conv2d(\n", - " mu[None, None], psf_normalized, padding=\"same\"\n", - " ).squeeze(0).squeeze(0)\n", - " \n", - " return mu" + "# Cosmology model\n", + "cosmology = caustics.FlatLambdaCDM(name=\"cosmo\")\n", + "# Source light model\n", + "source_light = caustics.Sersic(name=\"sourcelight\")\n", + "# Lens mass model\n", + "epl = caustics.EPL(name=\"epl\", cosmology=cosmology)\n", + "# Lens Light models\n", + "lens_light1 = caustics.Sersic(name=\"lenslight1\")\n", + "lens_light2 = caustics.Sersic(name=\"lenslight2\")\n", + "# PSF and image resolution\n", + "pixscale=0.11/2\n", + "fwhm = 0.269 # full width at half maximum of PSF\n", + "psf_sigma = fwhm / (2 * np.sqrt(2 * np.log(2)))\n", + "psf_image = torch.as_tensor(caustics.utils.gaussian(\n", + " nx=10,\n", + " ny=10,\n", + " pixelscale=pixscale,\n", + " sigma=psf_sigma,\n", + " upsample=1,\n", + "), dtype = torch.float32)\n", + "# Instantiate simulator\n", + "simulator = Doublelens(\n", + " lens=epl,\n", + " lens_light1=lens_light1,\n", + " lens_light2=lens_light2,\n", + " source=source_light,\n", + " psf=psf_image,\n", + " pixels_x=60*2,\n", + " pixelscale=pixscale,\n", + " upsample_factor=5)" + ] + }, + { + "cell_type": "code", + "execution_count": 68, + "id": "f4348821-5389-4574-b8f3-8f6d4fb25bfe", + "metadata": {}, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "%3\n", + "\n", + "\n", + "\n", + "140601093072176\n", + "\n", + "Doublelens(sim_3)\n", + "\n", + "\n", + "\n", + "140601103619248\n", + "\n", + "EPL(epl_3)\n", + "\n", + "\n", + "\n", + "140601093072176->140601103619248\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "140601093071504\n", + "\n", + "Sersic(sourcelight_3)\n", + "\n", + "\n", + "\n", + "140601093072176->140601093071504\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "140601093066272\n", + "\n", + "Sersic(lenslight1_3)\n", + "\n", + "\n", + "\n", + "140601093072176->140601093066272\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "140601093074864\n", + "\n", + "Sersic(lenslight2_1)\n", + "\n", + "\n", + "\n", + "140601093072176->140601093074864\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "140601093074960\n", + "\n", + "z_s\n", + "\n", + "\n", + "\n", + "140601093072176->140601093074960\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "140601093072416\n", + "\n", + "FlatLambdaCDM(cosmo_3)\n", + "\n", + "\n", + "\n", + "140601103619248->140601093072416\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "140601093065120\n", + "\n", + "z_l\n", + "\n", + "\n", + "\n", + "140601103619248->140601093065120\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "140601093077456\n", + "\n", + "x0\n", + "\n", + "\n", + "\n", + "140601103619248->140601093077456\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "140601093070880\n", + "\n", + "y0\n", + "\n", + "\n", + "\n", + "140601103619248->140601093070880\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "140601093068528\n", + "\n", + "q\n", + "\n", + "\n", + "\n", + "140601103619248->140601093068528\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "140601093070112\n", + "\n", + "phi\n", + "\n", + "\n", + "\n", + "140601103619248->140601093070112\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "140601093074624\n", + "\n", + "b\n", + "\n", + "\n", + "\n", + "140601103619248->140601093074624\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "140601093074720\n", + "\n", + "t\n", + "\n", + "\n", + "\n", + "140601103619248->140601093074720\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "140601093069824\n", + "\n", + "h0\n", + "\n", + "\n", + "\n", + "140601093072416->140601093069824\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "140601103333024\n", + "\n", + "critical_density_0\n", + "\n", + "\n", + "\n", + "140601093072416->140601103333024\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "140601094748384\n", + "\n", + "Om0\n", + "\n", + "\n", + "\n", + "140601093072416->140601094748384\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "140601093069872\n", + "\n", + "x0\n", + "\n", + "\n", + "\n", + "140601093071504->140601093069872\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "140601093078608\n", + "\n", + "y0\n", + "\n", + "\n", + "\n", + "140601093071504->140601093078608\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "140601093078896\n", + "\n", + "q\n", + "\n", + "\n", + "\n", + "140601093071504->140601093078896\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "140601093078464\n", + "\n", + "phi\n", + "\n", + "\n", + "\n", + "140601093071504->140601093078464\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "140601093077936\n", + "\n", + "n\n", + "\n", + "\n", + "\n", + "140601093071504->140601093077936\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "140601093068480\n", + "\n", + "Re\n", + "\n", + "\n", + "\n", + "140601093071504->140601093068480\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "140601093073040\n", + "\n", + "Ie\n", + "\n", + "\n", + "\n", + "140601093071504->140601093073040\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "140601093071216\n", + "\n", + "x0\n", + "\n", + "\n", + "\n", + "140601093066272->140601093071216\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "140601093070400\n", + "\n", + "y0\n", + "\n", + "\n", + "\n", + "140601093066272->140601093070400\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "140601093074336\n", + "\n", + "q\n", + "\n", + "\n", + "\n", + "140601093066272->140601093074336\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "140601093074384\n", + "\n", + "phi\n", + "\n", + "\n", + "\n", + "140601093066272->140601093074384\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "140601093073088\n", + "\n", + "n\n", + "\n", + "\n", + "\n", + "140601093066272->140601093073088\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "140601093074240\n", + "\n", + "Re\n", + "\n", + "\n", + "\n", + "140601093066272->140601093074240\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "140601093070784\n", + "\n", + "Ie\n", + "\n", + "\n", + "\n", + "140601093066272->140601093070784\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "140601093066704\n", + "\n", + "x0\n", + "\n", + "\n", + "\n", + "140601093074864->140601093066704\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "140601093067808\n", + "\n", + "y0\n", + "\n", + "\n", + "\n", + "140601093074864->140601093067808\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "140601093074096\n", + "\n", + "q\n", + "\n", + "\n", + "\n", + "140601093074864->140601093074096\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "140601093068960\n", + "\n", + "phi\n", + "\n", + "\n", + "\n", + "140601093074864->140601093068960\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "140601093074144\n", + "\n", + "n\n", + "\n", + "\n", + "\n", + "140601093074864->140601093074144\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "140601093074192\n", + "\n", + "Re\n", + "\n", + "\n", + "\n", + "140601093074864->140601093074192\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "140601093077504\n", + "\n", + "Ie\n", + "\n", + "\n", + "\n", + "140601093074864->140601093077504\n", + "\n", + "\n", + "\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 68, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "simulator.graphviz()" + ] + }, + { + "cell_type": "markdown", + "id": "2b3fc745-943d-4d48-9dd7-15f0f3a7f813", + "metadata": {}, + "source": [ + "When passing parameters to the `forward`, we need to use " + ] + }, + { + "cell_type": "code", + "execution_count": 69, + "id": "f66634dd-c33a-4358-8230-17885bee90dd", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "sim_3|module\n", + " epl_3|module\n", + " cosmo_3|module\n", + " h0|static\n", + " critical_density_0|static\n", + " Om0|static\n", + " z_l|dynamic\n", + " x0|dynamic\n", + " y0|dynamic\n", + " q|dynamic\n", + " phi|dynamic\n", + " b|dynamic\n", + " t|dynamic\n", + " sourcelight_3|module\n", + " x0|dynamic\n", + " y0|dynamic\n", + " q|dynamic\n", + " phi|dynamic\n", + " n|dynamic\n", + " Re|dynamic\n", + " Ie|dynamic\n", + " lenslight1_3|module\n", + " x0|dynamic\n", + " y0|dynamic\n", + " q|dynamic\n", + " phi|dynamic\n", + " n|dynamic\n", + " Re|dynamic\n", + " Ie|dynamic\n", + " lenslight2_1|module\n", + " x0|dynamic\n", + " y0|dynamic\n", + " q|dynamic\n", + " phi|dynamic\n", + " n|dynamic\n", + " Re|dynamic\n", + " Ie|dynamic\n", + " z_s|dynamic\n" + ] + } + ], + "source": [ + "print(simulator)" + ] + }, + { + "cell_type": "code", + "execution_count": 70, + "id": "5b34f5b4-16b1-4300-a99b-dcb1196c31e5", + "metadata": {}, + "outputs": [], + "source": [ + "params_for_simulator = torch.tensor([\n", + " #Lens redshift and mass\n", + " 1.5, #z_l\n", + " 0.25, # epl x0\n", + " 0.3, # epl y0\n", + " 1/1.14, # epl q\n", + " np.pi/2 + 1.6755160819145565, # epl phi\n", + " 0.97, # epl b\n", + " 1.04, # epl t\n", + " #Source light\n", + " 0.25, #x0\n", + " 0.3, #y0\n", + " 1-0.29, #q\n", + " -30/180*torch.pi, #phi\n", + " 4, #n\n", + " 0.1, #Re\n", + " 36, #Ie\n", + " #Lens light\n", + " 0.25, #x0\n", + " 0.1, #y0\n", + " 1-0.29, #q\n", + " -30/180*torch.pi, #phi\n", + " 4, #n\n", + " 0.1, #Re\n", + " 100, #Ie\n", + " #Lens light\n", + " 0.25, #x0\n", + " 0.6, #y0\n", + " 1-0.29, #q\n", + " -30/180*torch.pi, #phi\n", + " 4, #n\n", + " 0.1, #Re\n", + " 100, #Ie\n", + " #Source redshift\n", + " 3.5 #z_s\n", + " ])" ] + }, + { + "cell_type": "code", + "execution_count": 71, + "id": "631655f6-5516-4551-9cd9-d254fb77cb91", + "metadata": {}, + "outputs": [], + "source": [ + "lensed_image = simulator.forward(params_for_simulator)" + ] + }, + { + "cell_type": "code", + "execution_count": 72, + "id": "07ee42df-728c-458c-8842-f2811187613b", + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAakAAAGhCAYAAADbf0s2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA2R0lEQVR4nO3de3Bc5X3/8c/Zi1YXS/INSxa2wfyq/riYJMSmnhoaOwWcaYCEYZoQLgEm/QNqIChuY3CdNoYpErit62mcmDHTIbTEhekEEtqhrZUmMWHcFsfgBEx/kExcYy6KwMi6WNKuVvv8/qDs8z1rLch4hR7J79eMZo7Onl2ds758dT77fZ4ncs45AQAQoMRknwAAAOVQpAAAwaJIAQCCRZECAASLIgUACBZFCgAQLIoUACBYFCkAQLAoUgCAYFGkAADBmtQi9a1vfUuLFy9WdXW1li5dqp/85CeTeToAgMBMWpF69NFH1dbWpg0bNui5557T7/zO7+j3fu/39Morr0zWKQEAAhNN1gSzy5cv18c//nFt27atuO+ss87SFVdcoY6Ojvd8bqFQ0Ouvv676+npFUTTRpwoAqDDnnPr7+9XS0qJEovz9UupDPKeiXC6nvXv36s4774ztX716tXbv3n3M8dlsVtlstvj9a6+9prPPPnvCzxMAMLEOHTqkBQsWlH18UorUW2+9pdHRUTU1NcX2NzU1qaur65jjOzo6dNdddx2z/0J9WimlJ+w8AQATI68RPa0nVV9f/57HTUqReldpVOecGzO+W79+vdauXVv8vq+vTwsXLlRKaaUiihQATDn/+0HT+31kMylFau7cuUomk8fcNXV3dx9zdyVJmUxGmUzmwzo9AEAgJqW7r6qqSkuXLlVnZ2dsf2dnp1asWDEZpwQACNCkxX1r167VF7/4RS1btky//du/re3bt+uVV17RzTffPFmnBAAIzKQVqauuukqHDx/W3XffrTfeeENLlizRk08+qdNOO22yTgkAEJhJGyd1Ivr6+tTY2KhV+iyNEwAwBeXdiH6s76u3t1cNDQ1lj2PuPgBAsChSAIBgUaQAAMGiSAEAgkWRAgAEiyIFAAgWRQoAECyKFAAgWBQpAECwKFIAgGBRpAAAwaJIAQCCRZECAASLIgUACBZFCgAQLIoUACBYFCkAQLAoUgCAYFGkAADBokgBAIJFkQIABIsiBQAIFkUKABAsihQAIFgUKQBAsChSAIBgUaQAAMGiSAEAgkWRAgAEiyIFAAgWRQoAECyKFAAgWBQpAECwKFIAgGBRpAAAwaJIAQCCRZECAASLIgUACBZFCgAQLIoUACBYFCkAQLAoUgCAYFGkAADBokgBAIJFkQIABIsiBQAIFkUKABAsihQAIFgUKQBAsChSAIBgUaQAAMGqeJHq6OjQ+eefr/r6es2bN09XXHGFXnrppdgxzjlt3LhRLS0tqqmp0apVq7R///5KnwoAYIqreJHatWuXbrnlFv3nf/6nOjs7lc/ntXr1ah09erR4zKZNm7R582Zt3bpVe/bsUXNzsy655BL19/dX+nQAAFNY5JxzE/kD3nzzTc2bN0+7du3SJz7xCTnn1NLSora2Nt1xxx2SpGw2q6amJt1333266aab3vc1+/r61NjYqFX6rFJReiJPHwAwAfJuRD/W99Xb26uGhoayx034Z1K9vb2SpNmzZ0uSDhw4oK6uLq1evbp4TCaT0cqVK7V79+4xXyObzaqvry/2BQCY/ia0SDnntHbtWl144YVasmSJJKmrq0uS1NTUFDu2qamp+Fipjo4ONTY2Fr8WLlw4kacNAAjEhBapW2+9VT//+c/1D//wD8c8FkVR7Hvn3DH73rV+/Xr19vYWvw4dOjQh5wsACEtqol74tttu0xNPPKGnnnpKCxYsKO5vbm6W9M4d1fz584v7u7u7j7m7elcmk1Emk5moUwUABKrid1LOOd1666167LHH9MMf/lCLFy+OPb548WI1Nzers7OzuC+Xy2nXrl1asWJFpU8HADCFVfxO6pZbbtGOHTv0/e9/X/X19cXPmRobG1VTU6MoitTW1qb29na1traqtbVV7e3tqq2t1TXXXFPp0wEATGEVL1Lbtm2TJK1atSq2/8EHH9SNN94oSVq3bp2Ghoa0Zs0a9fT0aPny5dq5c6fq6+srfToAgClswsdJTQTGSQHA1BbMOCkAAD4oihQAIFgUKQBAsChSAIBgUaQAAMGiSAEAgkWRAgAEiyIFAAgWRQoAECyKFAAgWBQpAECwKFIAgGBRpAAAwaJIAQCCRZECAASLIgUACBZFCgAQLIoUACBYFCkAQLAoUgCAYFGkAADBokgBAIJFkQIABIsiBQAIFkUKABAsihQAIFgUKQBAsChSAIBgUaQAAMGiSAEAgkWRAgAEiyIFAAgWRQoAECyKFAAgWBQpAECwKFIAgGBRpAAAwaJIAQCCRZECAASLIgUACBZFCgAQLIoUACBYFCkAQLBSk30CACZYFMW/d25yzgP4ALiTAgAEiyIFAAgWcR8QktJobqr8DCJETBDupAAAwaJIAQCCRdwHTKQTidaiQH6HdIX3P2Y810kkiA9gwv8VdHR0KIoitbW1Ffc557Rx40a1tLSopqZGq1at0v79+yf6VAAAU8yEFqk9e/Zo+/bt+shHPhLbv2nTJm3evFlbt27Vnj171NzcrEsuuUT9/f0TeToAgClmworUwMCArr32Wj3wwAOaNWtWcb9zTlu2bNGGDRt05ZVXasmSJXrooYc0ODioHTt2TNTpABMnisp/xY5LjPkVJZPv/5VOFb8SVemKfdnXHdd5mK9y13PC7xNgTFiRuuWWW3TppZfq4osvju0/cOCAurq6tHr16uK+TCajlStXavfu3WO+VjabVV9fX+wLADD9TUjjxCOPPKJnn31We/bsOeaxrq4uSVJTU1Nsf1NTkw4ePDjm63V0dOiuu+6q/IkCAIJW8SJ16NAh3X777dq5c6eqq6vLHheV3OI7547Z967169dr7dq1xe/7+vq0cOHCypwwMF7jiaVK4q4oEY35WJQ0xyWTfn/a/JNMpcz+tN9fZbbNc13Kb8fOtaSrLsqP+m9G8v6wXM7vz434/fm8Od7v16jp+jMdgK5QJqB5ry7Bcu8tHYEnvYoXqb1796q7u1tLly4t7hsdHdVTTz2lrVu36qWXXpL0zh3V/Pnzi8d0d3cfc3f1rkwmo0wmU+lTBQAEruKfSV100UV6/vnntW/fvuLXsmXLdO2112rfvn0644wz1NzcrM7OzuJzcrmcdu3apRUrVlT6dAAAU1jF76Tq6+u1ZMmS2L66ujrNmTOnuL+trU3t7e1qbW1Va2ur2tvbVVtbq2uuuabSpwMcv+OM9WKRnonfJL3TBffutonvZJKBKFPl95soz9X4Y0Zrq8y2f518jd8ezfjzcEkb98VPPZn1sVtq0Ed/qQEf9yX6h/32oN92w1l/3lm/HYsER0fH3paJJgslJ1UuCiQGPOlNyowT69at09DQkNasWaOenh4tX75cO3fuVH19/WScDgAgUJFzU+9Xkr6+PjU2NmqVPqtUlH7/JwDHI8A7qcIk3klFZe6kVOZOypW9kzKnNN47qXKm3n9bKJF3I/qxvq/e3l41NDSUPS6QycEAADgWE8wCUvm7p3HcMcXulqpK7uztHVO133Zme3SG3fZ3TCMz/D/PbIP/GdlGfx45s52v83cXBXNzVnonlRr0r1XV68+3+m3/pOrDNf4SevwdU7J3yF/PgN/WOD6rsndYkeJ3WMfdtv4eLfaYXriTAgAEiyIFAAgWcR9OXuOJ+GKxntlfZfI0G+nVxmdZceb70Tp/XN7GevUm1mv0PyM702z7OZo1fIqPyhJzfbR2yiy/isCcmsHidsHFr/PXAzOK22+/5TtqB7v8OdW+4WPA2m7/HlS/5fdXHfHXk+j1Py8aNJHn8LDZNk0Xpe/9eGavIPo7KXEnBQAIFkUKABAs4j6cXKIyE76W69yzk7xW2049H+PZSK9QF4/78vX+OSMN/rVy9Umzbbr1Gvx2dqaPrnJzfcRXO+9ocfuseb8ubp8/068i8BsZvz8ZxWOyX2T9HJn/0XhGcXt/rZ9Lsz9TW9werfbvU96M6aqp9teQyfhrSx3xcV80YMaNJcy26QCUJJnJbZ1p/IsSRH8nO+6kAADBokgBAIJF3Ifp73gjPtO5Z2M927nnZpg4zEZ6M+KDeW3n3kid/3kjtXYQrtk2aWHBvlSZRkTbuTfqxv6dM6F4NFab8NHanIzvymuc4Qfnvtno34NszpxIwU69ZNey8sdnUv6YVMqfUyLht10ifkGx7+wAYLPeFdHfyYk7KQBAsChSAIBgEfdhehpHxBfr3LODc2tMxFfnY71CvZ/PLta1V+/jsHxd/Pe+fLWJ8uws5eawyCzFlDLT4UUmWkvkfbQ2lPODcX826M/7YK8f8XtK3W8Ut6uT5gdIOpr3zzl81F9fX7/fVs5Ec+a5o2aB7JEZ5r0smOjPvPcZs50224nSwbz2ezsA2CD6OzlxJwUACBZFCgAQLOI+TB8nEPFFtT7Ki3XuNfrtWMQ3wy486F+/kI7HWDb6SvqmOqWyJnIy8VNkd5tfIUfN68YG15rzHjLbv6qZ418nvg5jTGQGzqZG/M+oGhn7mIRJDu2UgKMmLbWdi1HBvN/Ov3+l//HEflu274ddHLFgIr7YSh+x0b/mdY5zIUUEiTspAECwKFIAgGAR92HqKrfUhsY5UNdGfPV1xe2CifhGGszyGnX+dWzXnkuaeCu+4KzSWR85JXM+ukrk/P7EiNk/4vfbqMtGa84MkC1UmW0bCdr9pf/KzfsWe92E+Xk2NTPXZ5r4VGbscOz4vIkmEyYi1Wi8wy5loryEjfUKZWJRmTzSlTm+HDr9phTupAAAwaJIAQCCRdyHqeU9Ir5YZ1eZiM8O1C0b8TXagbr+dWyEZtm4LjUc7yhLDvnvU4M+okoM+e0oZ1rm7IDVvMkObSwV62I00V3KLoth34v4edu40B7n0v75hbSJCzM25jTvh4nyRqts/Gl+mO0AzJiOy5I5Dm20GeVN5GlX7C3498MV7Eq+Nga0x9DpNx1wJwUACBZFCgAQLOI+TF1R/HesqOyKumPPxTfa4Lv7bMSXMyvo2hjLDrS1EV/6qI+YUgOm60xScsCvQBsNmu1hsxJt1ozyNctUFMy2ysRbZdlIMFH6PpnvzfuUSPsILmlW4HVmO1Xr36fROjtnoXnPqsfuOIwNTs7EY9tEjTmPEf8zEibyjOz7MTp2FBrv7hvHIF86/YLHnRQAIFgUKQBAsIj7EL4yHX1R6equJsaKDdqt87FeocF28fkYMNdo5uLLlIn4zGDcdL+PntJ9Pq5L9pq1NiRFg37ZCTfkH3M5Hws6uxKtjbHGMzC1nMR7dUHaKNCch5lcMBr20aSqTAw4bKK4Yf/+JbN+O2869/K1Y3dHlq7Mazv/ErUmqjXRX2Q6H2U7H20HoHn/IjNSeVyDfBEk7qQAAMGiSAEAgkXch6nFdmmVdveVmZevYJbeGJlpIr6Z/q+/XV7CLq9hB+emzMDcWMR3ZND/3P6jsXMqG/GVi/XKrixbud8nXWwpDDtY1i6FYWIzs22jSTvY2HbhpQr+vZdZnsMu8Vuoisd9tgvQDhiORnx0GOXKRX9jd/2V7fQrpzRWptsvCNxJAQCCRZECAASLuA9hGkdHX5Qu+eub8XGQMx19o6aLLzvLx0e5GWbuOTOVnF1BNxnr6PNxXbLfd+3ZiM/Ge5LkTJecs11o45lLrlzEZ9+D95rL8ETY6M/GaWUiMHsWCRuX2rkFE+bPJxFfLtjGfXb5kVE7yDdXZbbNoOkRs10m+ot3+jGn31TCnRQAIFgUKQBAsChSAIBg8ZkUwlem7TyqKlmTqNZ/9mQnj83O8p9lZBv9823becIsZ54YMutDmbbz2GSxA6a13H7ulItPMHtCn0PZz57M5Ln287rY/nKzTJTMtmBb0MtOXFumNT32Ouazu9g5mQltE2aNq2SVWa+qZG0u+5lUbFJf8/lRImfa0e3MF1nznpf5rIpm8qmLOykAQLAoUgCAYBH3IRzH23ZuWs6leNt5bqZ/bHi2j5myM02rs3mp1KCJ/kZMxDdoZlgw60HFZpLImv12JgTpuCM+e312TSx7rZFZ30n2/Si3NlLpOY2MHY/FWs2tcUR/sdex63oNmxbyIR/X2fWjJCmqNn8uZqLgvFmbKsqb6HDYv5Yz2zI/TzYKtbHreGafQDC4kwIABIsiBQAIFnEfwlRmtoXYsvA11bHHRmb6uG94ro+Ahk7xUVKuwbyWSX1SfgIJJUzqlRw2cZ9d8r3cZLHjncHARnzJsTsWo7o6/7Izxl72Pm+WcLdrNCVydkl724YX70yMjpoZMsx74GysZ5efH8+EtLbDLu23E7m82S55n5zp/LOTUaRsp5+J+7L+uqsGzcTCtuMz64935g913LNPsLR8ELiTAgAEiyIFAAgWcR8m13gmSLXRn4lzbDefJGVn+w64o03+OUNNPqrJz7AxmImYek1UljfrLZmIyuVM3Gc74ca7NHmZiE9pE/FV+wjTNc4obuea/PbgPDs42URXJiZLm2Wtag7HBz1Xd5vYzJ676fqzg4TLdv1ZNg6z8ad9rukyjEbj0VpkU0QT8RXMqZuUTkkT5aVq/fuRHLQzBZsuQzvR7ZgXgFBxJwUACBZFCgAQLOI+BC/W/WYGtY40xLv7hub6eGdovgl1FvgOttqMj7QGE757rtBt/imYp0Z5k0PZrrVyg1pLlZmLz3bMxQbtmmXvR2b7jr6BFh9p9Z/mnzt8iumwS/oTr+rx78VoVfyfeWLEv4fVgz7CjAbNAGUb040nki3T6acyy83H3tfSlzKnm68ee/B1csRfX/qoGTBsu/vs+5rg9/GpakL+5F577TVdd911mjNnjmpra/Wxj31Me/fuLT7unNPGjRvV0tKimpoarVq1Svv375+IUwEATGEVL1I9PT264IILlE6n9S//8i968cUX9Vd/9VeaOXNm8ZhNmzZp8+bN2rp1q/bs2aPm5mZdcskl6u/vr/TpAACmsIrHfffdd58WLlyoBx98sLjv9NNPL24757RlyxZt2LBBV155pSTpoYceUlNTk3bs2KGbbrqp0qeEKSgqE43JzFs30hDvWhuaa5Z7ONVHV+ee+kZxO5XwkdPPR1qK24Uq3z0XWwu9TNdabLmLcQ7gjS31brdNLOVqbJxpBySbbsVT/XnUtfhf7BKRP6e+an89qcH4+1Rz2EdlVSYKTCbjS7q/H1duKfnYezb24N+o5Km2M7FgpibM19mD7MBev9vGfakas8S87QQdtO+9/d2cefxCV/E7qSeeeELLli3T5z73Oc2bN0/nnXeeHnjggeLjBw4cUFdXl1avXl3cl8lktHLlSu3evXvM18xms+rr64t9AQCmv4oXqV/96lfatm2bWltb9W//9m+6+eab9eUvf1l/93d/J0nq6uqSJDU1NcWe19TUVHysVEdHhxobG4tfCxcurPRpAwACVPG4r1AoaNmyZWpvb5cknXfeedq/f7+2bdum66+/vnhcVNIx5Jw7Zt+71q9fr7Vr1xa/7+vro1BNR2Xm67ODMl21iftmxOOp3CyfIZ0670hx+xNzflHcHjVZ3qv9M4vbRzI+HnP2NMxg1/iKtrYF8Ph/17ODZW33YiFlts0KtaNmVRJX7SOqGdV+mZBY3Jf2nYGu5PRKv/c/cJwdi+/Hvo7t9Bt9j46+pL9W29GXazDvsz1v839FlRmUXVVjlgOxy57Yv0OKr56MsFX8Tmr+/Pk6++yzY/vOOussvfLKK5Kk5uZmSTrmrqm7u/uYu6t3ZTIZNTQ0xL4AANNfxYvUBRdcoJdeeim27+WXX9Zpp50mSVq8eLGam5vV2dlZfDyXy2nXrl1asWJFpU8HADCFVTzu+8pXvqIVK1aovb1dn//85/XMM89o+/bt2r59u6R3Yr62tja1t7ertbVVra2tam9vV21tra655ppKnw5CdJzz9cXmXUubQaqZ+OuMzPDR0KL6nuL2kupDxe3+gh8sW5/xUVlP0sZ35kVjEd8HiMPKdSmWY35eYsRvp8yKGok+/8/2zeqxU4VEj4+90gPxx1JDpstuxC4zUibOPIFlKmxEGntbk/E/O/tnOWIuKTfXnF/an3chZa/PP7f6LdPpZwf2Jsbxdw5BqniROv/88/X4449r/fr1uvvuu7V48WJt2bJF1157bfGYdevWaWhoSGvWrFFPT4+WL1+unTt3qr6+vtKnAwCYwiZkWqTLLrtMl112WdnHoyjSxo0btXHjxon48QCAaYK5+xA+E9U4E5kVSiIjJX0clDErsVZFYw/YzJtVWe3g0NjhrjKx13uJRWJZf97pAb9d85ZZViTtzzvXb+YvNINdM0f87trueExZ1WuW5Mia5Ufeo/turHONsV2QJtWMYiv82nn44tHnSJ1/LDvHP6dhvh+sXG86Gd+obixuD/eZ+Q7r/X9pGbvK8QfowEQY+JMDAASLIgUACBZxH6YUOy9cVLIibpTzkdGbWT849/9l5xe33877/YeP+gGvyWEzL1x+HBGf7RYb78q85ZgBrzZ+S/X6eKvWxGOJEf/PdqTWRJbmXNNHfWRWfTg+eDXZP+yfM2zivkKZuQnLGc9129cxXZqj1SUDsRvMysjzfCvjBaceKG63mAzzJ+nfKG7/8s0F/nVm+PfDZcyAabr7pizupAAAwaJIAQCCRdyHqcWs6JrKxuOm9ID/nevA27OL27ur/k9x+0jOR3z9PX67YdC/TjJnBruajjc3nkHI4xTr6LORmF29dshEf/0+urJjmFNDYy+vkRz2r5McjMd9toNwojoWx+JSpkMxU9LdZ4ZILjjFD8S+dNa+4vbpKb9/xKzt8XJjc3F7tNosz5Hkd/DpgD9FAECwKFIAgGAR92FKiUZ8VJUaLB2k6iOg/l/7Lr7n5Lu/Rkb8McnDdv43M2eeifticdh45t77IF1k45knzw5oNoOYC2m74qzZHDVdbumSSNDEbiq3WvCJsNdjf7Tt7quKv5f5arPMSl1vcfujVW8Vtxek/J/p4syb/kdU+78ThSTz9U033EkBAIJFkQIABIsiBQAIFp9JIXzmM474Z1LxiWOrjvi/zlVv+g9DBpz/LCMaNWsPHfG/o8XWWMqP/ZlUVObzm7LLsY+X/ezEtmnX+jXjczP99uA8f53ZxrHPw7bj11bFP5upMW31STPDhYb9DBexNbze7/zHy3ymd8x7Zr5PmNl+7aeOWedb6Y+M+uEDhZx/z8y8wopGy6yPhSmFOykAQLAoUgCAYBH3IXix2Rls3Hc0PpNCda9vPx45bCYXdWPPymCXVU/myswAYSO+ckvBj469XtU7Tx+7DTq2P+HPz2Wqitv5Rh/xHW32/1T7F5m1l+aZSWET/ryr3javmYpff3LY/4zqAf8zEgN+2o0TCscS79/WHpUsXZUwf5RvDPq1ovYMtxS36xN+4tmf9p7un2D+3NODdm0u86L2z8i9/7pZ7xxHRBgC7qQAAMGiSAEAgkXch/DZqGbELH0+FI/70v3+uMwR20lmZ2vwxyeHzSwT5RI7G1cl7UwNeXOIef0PEBFFZiLUgpkg1S6FPjTX/4zhU/11z27xszMkTdz3pl1e/aiP96T4UvSZdJnZJyokSo4dtSby8fcpOeh/9qtvzyxuP1n3kTGfv/e1hcXtml/79y/Ta/5+DPtt9x6RLMLGnRQAIFgUKQBAsIj7ECbbgVUwazrlTcxmB6JKSh31j1WZ9ZcKZun1UTOwNWmebpeit+tGxTrjEnZ77O6+ct18//ugeX6ZDjjzunby2IJJ7BK1/jpn1/qOt4TpyXs74wcwH9PcWO4Uy0WVhXF2w7378uUiUnPNUcng2pRf0V79h2uK28+kFhW37eTA+Vfqitszu/1rpW3cN+QHJ5eN+8bb6YdJw50UACBYFCkAQLCI+xC82GBeG/fl4t19STO4N33Ud8nZiC8yqY+d5y3285Jjz6UX2ejPdOQ5s3bT8UZjx/5w03E44rdTQ+aQPp/9vVo9c8yXKbztj0n3xx+z63Alcv5NcGZOv7JdiuOZA89GoSn/X4wr0+n3znmYpxzxxx0d9evKJwb969a97v+Mat80A7yP+NzQmbkI7fU45vGbUriTAgAEiyIFAAgWcR+mFBtJlcZ9iWGfGaUGzHxuGdvOZn4vM7vjHX1jL70eVZmlyXMmn0p+gIjPdhOaiNAOQE33+xir5k3b9efPKdtrliExKdaMt/311HXFz6+qx597NGjiMROl2o7F8QxQtgOSIxPxRWYuwkLGdlzGWwzt4N70gH8sNeyfY2PLWnNN1W/564nNP5j1cd+4Ykrm6gsSd1IAgGBRpAAAwSLuw4ev3FIYZmClK5j4yOZYNoYaicd9kYn7kkM+ZkodNctWmM690shprP2uynSnmbgvSpvoz55TaXffeCIkOxjYdKSl3/bR1QzbcDjkf3a+xr5/frNqwL9m9Zsm9pKUOuzXKHGDvm0wFvfZ6ygXlSXKDNo175Or8UuB2HkJ7UDl0nO3nYy2GzPT4w+qOWw7+sw1DJWJL+2gXQbwTincSQEAgkWRAgAEi7gP4Sszj5/yJaNxTbdfwizjkYrFTOb3spqxO/3sMaPVplPNrGgb6yy0HYcl0Vi5OePKDVC2A1Dtb5BVI/51UmYwr+0+tLFcwnQJJgbMxHiSNDh2PGaXQSnb0Wfn30vazkfz3tT4ufcKtT7uy9f440dL4j67Um9qcOxBzDVv+/egqse/T9FRcz3ZsefrYwDv1MWdFAAgWBQpAECwiPswpcRiqJHycZ/t9EsM+rgvVWUG6tpOv7QdzOu3R6t9RJWoNd19o9V+2w7MLekci/LmdUfLLD9ijy8TV8nEgMk+88/WdkeW6TJ0JYOej3vQbtmIz7wf1T7WczN83Dc6w+8frbFdliU/wpxSwiSQVf3+Ouwg5ESvGbRr475YDGvev3IdfQzgDR53UgCAYFGkAADBIu5DmI5zYK8kacTMp5f1UVRi2G8nh8z8cWkb/Zm58ewgXxP3RWZ5XHsa9je90uHBtlsvikwXX2ln4rv7bSQ4aubYs3MLjvnMkrjOdkGWRlqF91+So9zqurEuvlof67k6E/HV+yh0ZIb/LyZfXX7V4mTWzN035LczR0yXplmGg46+kwd3UgCAYFGkAADBIu7D5Co3j1/Z421UVXK87fYzEVA05P+aJ8zg16TZLiRtR5/p+jOr+uZlWtIiMzegeW6iZPXZhF3N10R/GjG/H+bHXh3XKjcouGx32jhXCLbdevb9j9Lmv4a0ifjKdfHZiK/BH5+vHXvexKjkclIm7qvqtfPymcHNdhmOYTNHn102hY6+aYc7KQBAsChSAIBgEfchfOPp9JPiUY8d1Gk77EyMZeM+O4A3NpjXrOprl8UomNcZNQOEU9Xxf1KpGtNlOGi7D03nXtYMQrbz59kOwNgyGnbwcJnoynbklcao5ZbYsCvqmi4+V21W1zVz8dmBurEuvlrbQTn2UiI23pOkqj5/feleM3C596h/uu3oG7YdffbvR5n3g4hvyuJOCgAQLIoUACBYFCkAQLAq/plUPp/Xxo0b9Z3vfEddXV2aP3++brzxRn3ta19TIvFOTXTO6a677tL27dvV09Oj5cuX65vf/KbOOeecSp8OppITaEeXJGfbmu1nOPbzH/O5i20Pj7Wjp+3sE2a/+XxqxH9Mo7yZlSJfF/+9L2k+q0nm/Gc4yaxZMj7nryORH3sS2mjUvDd5u99ctF3XKl+mZb2Esy3ysc/Z/Hahxuy329X+WvPVY8/YYT82TA3780sPxM/Pfg4Vnzy2TNt5rG1/HG3nmLIqfid133336f7779fWrVv13//939q0aZP+4i/+Qt/4xjeKx2zatEmbN2/W1q1btWfPHjU3N+uSSy5Rf39/pU8HADCFVbxI/cd//Ic++9nP6tJLL9Xpp5+u3//939fq1av105/+VNI7d1FbtmzRhg0bdOWVV2rJkiV66KGHNDg4qB07dlT6dAAAU1jF474LL7xQ999/v15++WX95m/+pn72s5/p6aef1pYtWyRJBw4cUFdXl1avXl18TiaT0cqVK7V7927ddNNNlT4lTCdl2tElKUqYx2wEZKM/ExlFZlLZZMJMNmtatAuxNadMy7qdicJPtqCRRElMaRPMUTtZrVmLyY19vBXZ5ebNpdlW7qSJ0xKxCDH+ova1nIlVbUwXizyr7P6xj4+9vokmU8Nm+6g/8XRfNvYcG/FF/abtfIiZJU52FS9Sd9xxh3p7e3XmmWcqmUxqdHRU99xzj66++mpJUldXlySpqakp9rympiYdPHhwzNfMZrPKmmlu+vr6Kn3aAIAAVTzue/TRR/Xwww9rx44devbZZ/XQQw/pL//yL/XQQw/FjisdYOicO3bQ4f/q6OhQY2Nj8WvhwoWVPm0AQIAqfif11a9+VXfeeae+8IUvSJLOPfdcHTx4UB0dHbrhhhvU3NwsScXOv3d1d3cfc3f1rvXr12vt2rXF7/v6+ihU093xdvopPttAbDYKO4uDfUJkZ6/wPyNpIr60ie9iy82n/HPtrBQjtfFzLWTMtvnXZpam0mjGmf1jb7ukuTYzOW1qyG+n+3ycmB6QOSYeeyVMahaVa4Ybx1tuY71kzsR65uelB8xksbaDr9/He5IU2cljB+1y8GNHfMwscfKo+J3U4OBgsdX8XclkUoX/baddvHixmpub1dnZWXw8l8tp165dWrFixZivmclk1NDQEPsCAEx/Fb+Tuvzyy3XPPfdo0aJFOuecc/Tcc89p8+bN+tKXviTpnd9Y29ra1N7ertbWVrW2tqq9vV21tbW65pprKn06AIAprOJF6hvf+Ib+9E//VGvWrFF3d7daWlp000036c/+7M+Kx6xbt05DQ0Nas2ZNcTDvzp07VV9fX+nTwXRQLvor7eoy8V2sC9A0gkUyE8+WWZLd5gCxfyD2Z0e+O69gJmm1nX5SPOLLNfifkp/pT6pqlu9gWzDLNwUtqn+7uH1Klc/vhkxW+D8Ds/32Yb/d+1ZtcTt9JL7GVarfRJumyS5pkrXEiD9XGw/aY2xnYeqoGajbb5Z877cRn4nxjsbjvnJdfLFBuzbis3/2RHzTWuTKTqMcrr6+PjU2NmqVPquU+c8CJ4H3+nzKfsZkW8HtftN2bmf6VsZ/eBTV+cX8Cg3+P/t8o9+fm+n/3g3P9EUgOzN+fiPm967xFKn5FSpSQx9CkUoP+UJxokVKFKmTTt6N6Mf6vnp7e9/zIxzm7gMABIv1pDB9lFt3KjF29CfzG7u9/ykb/Zlf3qNYd5lZeykRv2uxa1Pl68Z+sdkNfvDq+XP9WMGLGl4sbv/f9OHitg059zW2FLd31vm5L3+a8d2vb6fiv6U6mz6Yu6qEjUXNHZaN9dJHzXa/WQNqwN89JfvMgGm7BtSguVsyg6qlcXbxcfd0UuJOCgAQLIoUACBYxH2YWsY7yPcEor/YzzCvk7BLYZhIKootkWFG70pKjNr4zwwArvKR2+EZPgd8rWGm3183o7g9mDpiXsWf36gJKgvOzLdn5zXMx38XtY0QKZ/GKd3vX7cqtm3n3DNNEQMmLj1qIj4b65UZmGsbIt65ECI+jI07KQBAsChSAIBgEfdh6iqNf8rFf2WX9yiz3IN9XbM6bmzlWzMfYCrrtxPDfnySJKWGqs22j/gSWR8DDmZ93Pefw2cUt1/pn1XcPr3Bj5lKmza81482FrcPHvbH537tz6P6zXjHYbVvFFTmiL+mTK9dOdfMuWfHOg347VisN+z3O7NiQSziG7VRa8nKwSy3gTK4kwIABIsiBQAIFnEfpo/xdP7ZWMlMlxRb5kNlOs1sB5rpTrPRXyJn5gaUVDXoo6/UgJ9Wqarfx4DVPf6f4VC37w48PKu5uP3rOr+MjUvYKYvMsiJmYO7MHn+MjfHe+X7sbr3kUXMdtltvyMR3NtYbMVGeue5YlDeerr1SRHwwuJMCAASLIgUACBZxH6an443+YruPswPQrvxbEvdFQ34wa9LMY1fd6+O+qjd9DDjjdT8P4MgM/89zNGOiSXN6ibw/p+SwXTrDx5HJofg5JQbNIFwT5clGdrEOPbPfzkp+orGeRcSHMriTAgAEiyIFAAgWcR+mv3JR0onEgAUz95xdYHE0/lw7mDUynXHRUb88R+KIj/iSGb9dZeb3cykzIDe2OrHpSsyb+G3EnN9IPO5zNp4cGUd8Z68pNjD6OGM9Ij18ANxJAQCCRZECAASLuA8nr4mIAVUyJ5391kZopntOkVkvIzH2z47ea1mSd1/TbsfmH3yPmK3sNZV5DrEePmTcSQEAgkWRAgAEi7gPKFWhGPCYw2z0F5WJCEvSwuJzx/UTPoDxDrYd87nEeph43EkBAIJFkQIABIsiBQAIFp9JAeM1ns9gxtEq/s5rncBnQROFz5gQIO6kAADBokgBAIJF3AdUEpEZUFHcSQEAgkWRAgAEiyIFAAgWRQoAECyKFAAgWBQpAECwKFIAgGBRpAAAwaJIAQCCRZECAASLIgUACBZFCgAQLIoUACBYFCkAQLAoUgCAYFGkAADBokgBAIJFkQIABIsiBQAIFkUKABAsihQAIFjHXaSeeuopXX755WppaVEURfre974Xe9w5p40bN6qlpUU1NTVatWqV9u/fHzsmm83qtttu09y5c1VXV6fPfOYzevXVV0/oQgAA089xF6mjR4/qox/9qLZu3Trm45s2bdLmzZu1detW7dmzR83NzbrkkkvU399fPKatrU2PP/64HnnkET399NMaGBjQZZddptHR0Q9+JQCAaSdyzrkP/OQo0uOPP64rrrhC0jt3US0tLWpra9Mdd9wh6Z27pqamJt1333266aab1Nvbq1NOOUV///d/r6uuukqS9Prrr2vhwoV68skn9alPfep9f25fX58aGxu1Sp9VKkp/0NMHAEySvBvRj/V99fb2qqGhoexxFf1M6sCBA+rq6tLq1auL+zKZjFauXKndu3dLkvbu3auRkZHYMS0tLVqyZEnxmFLZbFZ9fX2xLwDA9FfRItXV1SVJampqiu1vamoqPtbV1aWqqirNmjWr7DGlOjo61NjYWPxauHBhJU8bABCoCenui6Io9r1z7ph9pd7rmPXr16u3t7f4dejQoYqdKwAgXBUtUs3NzZJ0zB1Rd3d38e6qublZuVxOPT09ZY8plclk1NDQEPsCAEx/FS1SixcvVnNzszo7O4v7crmcdu3apRUrVkiSli5dqnQ6HTvmjTfe0AsvvFA8BgAASUod7xMGBgb0y1/+svj9gQMHtG/fPs2ePVuLFi1SW1ub2tvb1draqtbWVrW3t6u2tlbXXHONJKmxsVF/8Ad/oD/6oz/SnDlzNHv2bP3xH/+xzj33XF188cWVuzIAwJR33EXqpz/9qT75yU8Wv1+7dq0k6YYbbtC3v/1trVu3TkNDQ1qzZo16enq0fPly7dy5U/X19cXn/PVf/7VSqZQ+//nPa2hoSBdddJG+/e1vK5lMVuCSAADTxQmNk5osjJMCgKltUsZJAQBQSRQpAECwKFIAgGBRpAAAwaJIAQCCRZECAASLIgUACBZFCgAQLIoUACBYFCkAQLAoUgCAYFGkAADBokgBAIJFkQIABIsiBQAIFkUKABAsihQAIFgUKQBAsChSAIBgUaQAAMGiSAEAgkWRAgAEiyIFAAgWRQoAECyKFAAgWBQpAECwKFIAgGBRpAAAwaJIAQCCRZECAASLIgUACBZFCgAQLIoUACBYFCkAQLAoUgCAYFGkAADBokgBAIJFkQIABIsiBQAIFkUKABAsihQAIFgUKQBAsChSAIBgUaQAAMGiSAEAgkWRAgAEiyIFAAgWRQoAECyKFAAgWBQpAECwKFIAgGBRpAAAwUpN9gl8EM45SVJeI5Kb5JMBABy3vEYk+f/Py5mSRaq/v1+S9LSenOQzAQCciP7+fjU2NpZ9PHLvV8YCVCgU9Prrr8s5p0WLFunQoUNqaGiY7NP6UPT19WnhwoUn1TVLXPfJdN0n4zVLJ991O+fU39+vlpYWJRLlP3makndSiURCCxYsUF9fnySpoaHhpPhDtU7Ga5a47pPJyXjN0sl13e91B/UuGicAAMGiSAEAgjWli1Qmk9HXv/51ZTKZyT6VD83JeM0S130yXffJeM3SyXvd72dKNk4AAE4OU/pOCgAwvVGkAADBokgBAIJFkQIABGvKFqlvfetbWrx4saqrq7V06VL95Cc/mexTqqiOjg6df/75qq+v17x583TFFVfopZdeih3jnNPGjRvV0tKimpoarVq1Svv375+kM668jo4ORVGktra24r7pes2vvfaarrvuOs2ZM0e1tbX62Mc+pr179xYfn27Xnc/n9bWvfU2LFy9WTU2NzjjjDN19990qFArFY6bDNT/11FO6/PLL1dLSoiiK9L3vfS/2+HiuMZvN6rbbbtPcuXNVV1enz3zmM3r11Vc/xKuYZG4KeuSRR1w6nXYPPPCAe/HFF93tt9/u6urq3MGDByf71CrmU5/6lHvwwQfdCy+84Pbt2+cuvfRSt2jRIjcwMFA85t5773X19fXuu9/9rnv++efdVVdd5ebPn+/6+vom8cwr45lnnnGnn366+8hHPuJuv/324v7peM1vv/22O+2009yNN97o/uu//ssdOHDA/eAHP3C//OUvi8dMt+v+8z//czdnzhz3z//8z+7AgQPuH//xH92MGTPcli1bisdMh2t+8skn3YYNG9x3v/tdJ8k9/vjjscfHc40333yzO/XUU11nZ6d79tln3Sc/+Un30Y9+1OXz+Q/5aibHlCxSv/Vbv+Vuvvnm2L4zzzzT3XnnnZN0RhOvu7vbSXK7du1yzjlXKBRcc3Ozu/fee4vHDA8Pu8bGRnf//fdP1mlWRH9/v2ttbXWdnZ1u5cqVxSI1Xa/5jjvucBdeeGHZx6fjdV966aXuS1/6UmzflVde6a677jrn3PS85tIiNZ5rPHLkiEun0+6RRx4pHvPaa6+5RCLh/vVf//VDO/fJNOXivlwup71792r16tWx/atXr9bu3bsn6awmXm9vryRp9uzZkqQDBw6oq6sr9j5kMhmtXLlyyr8Pt9xyiy699FJdfPHFsf3T9ZqfeOIJLVu2TJ/73Oc0b948nXfeeXrggQeKj0/H677wwgv17//+73r55ZclST/72c/09NNP69Of/rSk6XnNpcZzjXv37tXIyEjsmJaWFi1ZsmTavA/vZ8pNMPvWW29pdHRUTU1Nsf1NTU3q6uqapLOaWM45rV27VhdeeKGWLFkiScVrHet9OHjw4Id+jpXyyCOP6Nlnn9WePXuOeWy6XvOvfvUrbdu2TWvXrtWf/Mmf6JlnntGXv/xlZTIZXX/99dPyuu+44w719vbqzDPPVDKZ1OjoqO655x5dffXVkqbvn7U1nmvs6upSVVWVZs2adcwx0/X/u1JTrki9K4qi2PfOuWP2TRe33nqrfv7zn+vpp58+5rHp9D4cOnRIt99+u3bu3Knq6uqyx02na5beWXpm2bJlam9vlySdd9552r9/v7Zt26brr7++eNx0uu5HH31UDz/8sHbs2KFzzjlH+/btU1tbm1paWnTDDTcUj5tO11zOB7nG6fg+lDPl4r65c+cqmUwe81tEd3f3Mb+RTAe33XabnnjiCf3oRz/SggULivubm5slaVq9D3v37lV3d7eWLl2qVCqlVCqlXbt26W/+5m+USqWK1zWdrlmS5s+fr7PPPju276yzztIrr7wiaXr+WX/1q1/VnXfeqS984Qs699xz9cUvflFf+cpX1NHRIWl6XnOp8Vxjc3Ozcrmcenp6yh4z3U25IlVVVaWlS5eqs7Mztr+zs1MrVqyYpLOqPOecbr31Vj322GP64Q9/qMWLF8ceX7x4sZqbm2PvQy6X065du6bs+3DRRRfp+eef1759+4pfy5Yt07XXXqt9+/bpjDPOmHbXLEkXXHDBMcMLXn75ZZ122mmSpuef9eDg4DEL3SWTyWIL+nS85lLjucalS5cqnU7HjnnjjTf0wgsvTJv34X1NWsvGCXi3Bf1v//Zv3Ysvvuja2tpcXV2d+5//+Z/JPrWK+cM//EPX2NjofvzjH7s33nij+DU4OFg85t5773WNjY3usccec88//7y7+uqrp1yL7vux3X3OTc9rfuaZZ1wqlXL33HOP+8UvfuG+853vuNraWvfwww8Xj5lu133DDTe4U089tdiC/thjj7m5c+e6devWFY+ZDtfc39/vnnvuOffcc885SW7z5s3uueeeKw6XGc813nzzzW7BggXuBz/4gXv22Wfd7/7u79KCPhV885vfdKeddpqrqqpyH//4x4ut2dOFpDG/HnzwweIxhULBff3rX3fNzc0uk8m4T3ziE+7555+fvJOeAKVFarpe8z/90z+5JUuWuEwm484880y3ffv22OPT7br7+vrc7bff7hYtWuSqq6vdGWec4TZs2OCy2WzxmOlwzT/60Y/G/Hd8ww03OOfGd41DQ0Pu1ltvdbNnz3Y1NTXusssuc6+88sokXM3kYKkOAECwptxnUgCAkwdFCgAQLIoUACBYFCkAQLAoUgCAYFGkAADBokgBAIJFkQIABIsiBQAIFkUKABAsihQAIFgUKQBAsP4/Hwtwc4vgGNsAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.imshow(lensed_image)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cbb5a938-b38d-4c41-a058-97fc3790cfa2", + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": {