1+ import json
12import os
23from os import path
34import numpy as np
4- from typing import Optional
5+ from typing import List , Optional
6+
7+ from autonerves .dictable import to_dict
8+
9+ import autofit as af
510
611def plot_profile_1d (
712 xvalues : np .ndarray ,
@@ -51,4 +56,364 @@ def plot_profile_1d(
5156 os .makedirs (output_path )
5257 plt .savefig (output_path / f"{ output_filename } .png" )
5358 plt .clf ()
54- plt .close ()
59+ plt .close ()
60+
61+
62+ """
63+ The functions below simulate the example 1D datasets fitted throughout the PyAutoFit workspaces.
64+
65+ They live in the library rather than beside the simulator scripts because a Jupyter notebook kernel
66+ has no script directory on `sys.path` — a local `import util` resolves for `python
67+ scripts/simulators/simulators.py` but can never resolve in the generated notebook.
68+ """
69+
70+ PIXELS = 100
71+ SIGNAL_TO_NOISE_RATIO = 25.0
72+
73+
74+ def _kernel_1d () -> np .ndarray :
75+ """
76+ Return the normalized 1D Gaussian kernel that the blurred datasets are convolved with.
77+ """
78+ kernel_pixels = 21
79+ kernel_sigma = 5.0
80+ kernel_centre = 10.0
81+
82+ kernel_xvalues = np .subtract (np .arange (kernel_pixels ), kernel_centre )
83+
84+ kernel = np .multiply (
85+ np .divide (1.0 , kernel_sigma * np .sqrt (2.0 * np .pi )),
86+ np .exp (- 0.5 * np .square (np .divide (kernel_xvalues , kernel_sigma ))),
87+ )
88+
89+ return kernel / np .sum (kernel )
90+
91+
92+ def _noise_map_1d () -> np .ndarray :
93+ """
94+ Return the noise-map, which is required when evaluating the chi-squared value of the likelihood.
95+ """
96+ return (1.0 / SIGNAL_TO_NOISE_RATIO ) * np .ones (PIXELS )
97+
98+
99+ def _output_dataset (
100+ dataset_path : str ,
101+ data : np .ndarray ,
102+ noise_map : np .ndarray ,
103+ kernel : Optional [np .ndarray ] = None ,
104+ ):
105+ """
106+ Output the data, noise-map and (optionally) the convolution kernel to the `dataset` folder, so
107+ they can be loaded and used in other example scripts.
108+
109+ Parameters
110+ ----------
111+ dataset_path
112+ The folder the dataset .json files are output to.
113+ data
114+ The simulated data that is fitted.
115+ noise_map
116+ The noise-map of the simulated data.
117+ kernel
118+ The convolution kernel, output only for the blurred datasets.
119+ """
120+ af .util .numpy_array_to_json (
121+ array = data , file_path = path .join (dataset_path , "data.json" ), overwrite = True
122+ )
123+ af .util .numpy_array_to_json (
124+ array = noise_map ,
125+ file_path = path .join (dataset_path , "noise_map.json" ),
126+ overwrite = True ,
127+ )
128+
129+ if kernel is not None :
130+ af .util .numpy_array_to_json (
131+ array = kernel ,
132+ file_path = path .join (dataset_path , "kernel.json" ),
133+ overwrite = True ,
134+ )
135+
136+
137+ def _output_model_json (
138+ dataset_path : str ,
139+ profile_1d_list : List ,
140+ indexed : bool ,
141+ strict : bool = False ,
142+ ):
143+ """
144+ Output each profile to a .json file so its parameters can be referred to in the future.
145+
146+ Parameters
147+ ----------
148+ dataset_path
149+ The folder the model .json files are output to.
150+ profile_1d_list
151+ The profiles the dataset was simulated from.
152+ indexed
153+ If `True` the profiles are written as `model_{i}.json`, else the single profile is written
154+ as `model.json`. Set by the caller rather than inferred from the list length, so a
155+ one-profile list keeps the indexed filenames its dataset folder is read with.
156+ strict
157+ If `False` a profile that cannot be serialized is skipped rather than raising.
158+ """
159+ for i , profile in enumerate (profile_1d_list ):
160+ filename = f"model_{ i } .json" if indexed else "model.json"
161+
162+ with open (path .join (dataset_path , filename ), "w+" ) as f :
163+ if strict :
164+ json .dump (to_dict (profile ), f , indent = 4 )
165+ else :
166+ try :
167+ json .dump (to_dict (profile ), f , indent = 4 )
168+ except (TypeError , ValueError ):
169+ pass
170+
171+
172+ def _output_figure (
173+ dataset_path : str ,
174+ xvalues : np .ndarray ,
175+ data : np .ndarray ,
176+ noise_map : np .ndarray ,
177+ title : str ,
178+ model_data_1d : Optional [np .ndarray ] = None ,
179+ model_data_1d_list : Optional [List [np .ndarray ]] = None ,
180+ ):
181+ """
182+ Output a .png of the simulated dataset, optionally overlaying the model profiles it was
183+ simulated from.
184+
185+ Parameters
186+ ----------
187+ dataset_path
188+ The folder the `image.png` is output to.
189+ xvalues
190+ The x-coordinates the profile is defined on.
191+ data
192+ The simulated data that is fitted.
193+ noise_map
194+ The noise-map of the simulated data, plotted as the error bars.
195+ title
196+ The title of the figure.
197+ model_data_1d
198+ The summed model profile, overlaid in red if input.
199+ model_data_1d_list
200+ The individual model profiles, overlaid as dashed lines if input.
201+ """
202+ import matplotlib .pyplot as plt
203+
204+ plt .errorbar (
205+ x = xvalues ,
206+ y = data ,
207+ yerr = noise_map ,
208+ linestyle = "" ,
209+ color = "k" ,
210+ ecolor = "k" ,
211+ elinewidth = 1 ,
212+ capsize = 2 ,
213+ )
214+
215+ if model_data_1d is not None :
216+ plt .plot (range (data .shape [0 ]), model_data_1d , color = "r" )
217+
218+ if model_data_1d_list is not None :
219+ for model_data_1d_individual in model_data_1d_list :
220+ plt .plot (range (data .shape [0 ]), model_data_1d_individual , "--" )
221+
222+ plt .title (title )
223+ plt .xlabel ("x values of profile" )
224+ plt .ylabel ("Profile normalization" )
225+ plt .savefig (path .join (dataset_path , "image.png" ))
226+ plt .close ()
227+
228+
229+ def simulate_dataset_1d_via_gaussian_from (gaussian , dataset_path : str ):
230+ """
231+ Simulate a 1D dataset from a single `Gaussian` and output it to the `dataset` folder.
232+
233+ The `Gaussian` is evaluated on `PIXELS` xvalues to create its model profile, noise is added at
234+ a fixed signal-to-noise ratio, and the data, noise-map, an `image.png` and the model .json are
235+ written to `dataset_path`.
236+
237+ Parameters
238+ ----------
239+ gaussian
240+ The `Gaussian` profile the dataset is simulated from.
241+ dataset_path
242+ The folder the simulated dataset is output to.
243+ """
244+ xvalues = np .arange (PIXELS )
245+
246+ model_data_1d = gaussian .model_data_from (xvalues = xvalues )
247+
248+ noise = np .random .normal (0.0 , 1.0 / SIGNAL_TO_NOISE_RATIO , PIXELS )
249+
250+ data = model_data_1d + noise
251+ noise_map = _noise_map_1d ()
252+
253+ _output_dataset (dataset_path = dataset_path , data = data , noise_map = noise_map )
254+
255+ _output_figure (
256+ dataset_path = dataset_path ,
257+ xvalues = xvalues ,
258+ data = data ,
259+ noise_map = noise_map ,
260+ title = "1D Gaussian Dataset." ,
261+ )
262+
263+ _output_model_json (
264+ dataset_path = dataset_path , profile_1d_list = [gaussian ], indexed = False
265+ )
266+
267+
268+ def simulate_data_1d_with_kernel_via_gaussian_from (gaussian , dataset_path : str ):
269+ """
270+ Simulate a 1D dataset from a single `Gaussian` convolved with a Gaussian kernel, and output it
271+ to the `dataset` folder.
272+
273+ Identical to `simulate_dataset_1d_via_gaussian_from` except the model profile is convolved
274+ before the noise is added, and the kernel is output alongside the data so the fit can deconvolve
275+ it.
276+
277+ Parameters
278+ ----------
279+ gaussian
280+ The `Gaussian` profile the dataset is simulated from.
281+ dataset_path
282+ The folder the simulated dataset is output to.
283+ """
284+ xvalues = np .arange (PIXELS )
285+
286+ model_data_1d = gaussian .model_data_from (xvalues = xvalues )
287+
288+ kernel = _kernel_1d ()
289+
290+ blurred_model_data_1d = np .convolve (model_data_1d , kernel , mode = "same" )
291+
292+ noise = np .random .normal (0.0 , 1.0 / SIGNAL_TO_NOISE_RATIO , PIXELS )
293+
294+ data = blurred_model_data_1d + noise
295+ noise_map = _noise_map_1d ()
296+
297+ _output_dataset (
298+ dataset_path = dataset_path , data = data , noise_map = noise_map , kernel = kernel
299+ )
300+
301+ _output_figure (
302+ dataset_path = dataset_path ,
303+ xvalues = xvalues ,
304+ data = data ,
305+ noise_map = noise_map ,
306+ title = "1D Gaussian Dataset with Convolver Blurring." ,
307+ )
308+
309+ _output_model_json (
310+ dataset_path = dataset_path ,
311+ profile_1d_list = [gaussian ],
312+ indexed = False ,
313+ strict = True ,
314+ )
315+
316+
317+ def simulate_dataset_1d_via_profile_1d_list_from (profile_1d_list : List , dataset_path : str ):
318+ """
319+ Simulate a 1D dataset from a list of profiles and output it to the `dataset` folder.
320+
321+ Every profile is evaluated on `PIXELS` xvalues and summed to create the overall model profile.
322+ The maximum log likelihood of the dataset (that of the profiles it was simulated from) is
323+ output alongside it as `max_log_likelihood.json`.
324+
325+ Parameters
326+ ----------
327+ profile_1d_list
328+ The profiles which are summed to simulate the dataset.
329+ dataset_path
330+ The folder the simulated dataset is output to.
331+ """
332+ xvalues = np .arange (PIXELS )
333+
334+ model_data_1d_list = [
335+ profile_1d .model_data_from (xvalues = xvalues ) for profile_1d in profile_1d_list
336+ ]
337+
338+ model_data_1d = sum (model_data_1d_list )
339+
340+ noise = np .random .normal (0.0 , 1.0 / SIGNAL_TO_NOISE_RATIO , PIXELS )
341+
342+ data = model_data_1d + noise
343+ noise_map = _noise_map_1d ()
344+
345+ _output_dataset (dataset_path = dataset_path , data = data , noise_map = noise_map )
346+
347+ _output_figure (
348+ dataset_path = dataset_path ,
349+ xvalues = xvalues ,
350+ data = data ,
351+ noise_map = noise_map ,
352+ title = "1D Profiles Dataset." ,
353+ model_data_1d = model_data_1d ,
354+ model_data_1d_list = model_data_1d_list ,
355+ )
356+
357+ _output_model_json (
358+ dataset_path = dataset_path , profile_1d_list = profile_1d_list , indexed = True
359+ )
360+
361+ chi_squared = np .sum (((data - model_data_1d ) / noise_map ) ** 2 )
362+ noise_normalization = np .sum (np .log (2 * np .pi * noise_map ** 2.0 ))
363+ log_likelihood = - 0.5 * (chi_squared + noise_normalization )
364+
365+ with open (path .join (dataset_path , "max_log_likelihood.json" ), "w+" ) as f :
366+ json .dump ({"log_likelihood" : log_likelihood }, f , indent = 4 )
367+
368+
369+ def simulate_data_1d_with_kernel_via_profile_1d_list_from (
370+ profile_1d_list : List , dataset_path : str
371+ ):
372+ """
373+ Simulate a 1D dataset from a list of profiles convolved with a Gaussian kernel, and output it
374+ to the `dataset` folder.
375+
376+ Identical to `simulate_dataset_1d_via_profile_1d_list_from` except the summed model profile is
377+ convolved before the noise is added, and the kernel is output alongside the data.
378+
379+ Parameters
380+ ----------
381+ profile_1d_list
382+ The profiles which are summed to simulate the dataset.
383+ dataset_path
384+ The folder the simulated dataset is output to.
385+ """
386+ xvalues = np .arange (PIXELS )
387+
388+ model_data_1d = np .zeros (shape = PIXELS )
389+
390+ for profile in profile_1d_list :
391+ model_data_1d += profile .model_data_from (xvalues = xvalues )
392+
393+ kernel = _kernel_1d ()
394+
395+ blurred_model_data_1d = np .convolve (model_data_1d , kernel , mode = "same" )
396+
397+ noise = np .random .normal (0.0 , 1.0 / SIGNAL_TO_NOISE_RATIO , PIXELS )
398+
399+ data = blurred_model_data_1d + noise
400+ noise_map = _noise_map_1d ()
401+
402+ _output_dataset (
403+ dataset_path = dataset_path , data = data , noise_map = noise_map , kernel = kernel
404+ )
405+
406+ _output_figure (
407+ dataset_path = dataset_path ,
408+ xvalues = xvalues ,
409+ data = data ,
410+ noise_map = noise_map ,
411+ title = "1D Profiles Dataset with Convolver Blurring." ,
412+ )
413+
414+ _output_model_json (
415+ dataset_path = dataset_path ,
416+ profile_1d_list = profile_1d_list ,
417+ indexed = True ,
418+ strict = True ,
419+ )
0 commit comments