LSSTDESC/Spectractor | Build 4890731301 | spectractor/fit/fit_multispectra.py | Coveralls

1

import matplotlib.pyplot as plt

2

from matplotlib import cm

3

from matplotlib import colors

4

from matplotlib.ticker import MaxNLocator

5

from astropy.io import ascii

6

7

import copy

8

import os

9

import numpy as np

10

from scipy.interpolate import interp1d

11

from scipy.integrate import quad

12

13

from spectractor import parameters

14

from spectractor.config import set_logger

15

from spectractor.simulation.atmosphere import Atmosphere, AtmosphereGrid

16

from spectractor.fit.fitter import (FitWorkspace, run_minimisation_sigma_clipping, run_minimisation, FitParameters,

17

                                    write_fitparameter_json)

18

from spectractor.tools import from_lambda_to_colormap, fftconvolve_gaussian

19

from spectractor.extractor.spectrum import Spectrum

20

from spectractor.extractor.spectroscopy import HALPHA, HBETA, HGAMMA, HDELTA, O2_1, O2_2, O2B

21

22

23

class MultiSpectraFitWorkspace(FitWorkspace):

24

25

    def __init__(self, output_file_name, file_names, fixed_A1s=True, inject_random_A1s=False, bin_width=-1,

26

                 verbose=False, plot=False, live_fit=False):

27

        """Class to fit jointly multiple spectra extracted with Spectractor.

28

29

        The spectrum is supposed to be the product of the star SED, a common instrumental throughput,

30

        a grey term (clouds) and a common atmospheric transmission, with the second order diffraction removed.

31

        The truth parameters are loaded from the file header if provided.

32

        If provided, the atmospheric grid files are used for the atmospheric transmission simulations and interpolated

33

        with splines, otherwise Libradtran is called at each step (slower). The files should have the same name as

34

        the spectrum files but with the atmsim suffix.

35

36

        Parameters

37

        ----------

38

        output_file_name: str

39

            Generic file name to output results.

40

        file_names: list

41

            List of spectrum file names.

42

        bin_width: float

43

            Size of the wavelength bins in nm. If negative, no binning.

44

        verbose: bool, optional

45

            Verbosity level (default: False).

46

        plot: bool, optional

47

            If True, many plots are produced (default: False).

48

        live_fit: bool, optional

49

            If True, many plots along the fitting procedure are produced to see convergence in live (default: False).

50

51

        Examples

52

        --------

53

        >>> from spectractor.config import load_config

54

        >>> load_config("./config/ctio.ini")

55

        >>> file_names = ["./tests/data/reduc_20170530_134_spectrum.fits"]

56

        >>> w = MultiSpectraFitWorkspace("./outputs/test", file_names, bin_width=5, verbose=True)

57

        >>> w.output_file_name

58

        './outputs/test'

59

        >>> w.spectra  #doctest: +ELLIPSIS

60

        [<spectractor.extractor.spectrum.Spectrum object at ...>]

61

        >>> w.lambdas  #doctest: +ELLIPSIS

62

        array([[ ...

63

"""

64

        # load spectra

65

        for name in file_names:

66

            if "spectrum" not in name:

67

                raise ValueError(f"ALl file names must contain spectrum keyword and be an output from Spectractor. "

68

                                 f"I found {name} in file_names list.")

69

        self.my_logger = set_logger(self.__class__.__name__)

70

        self.spectrum = Spectrum(file_names[0])

71

        self.spectra = []

72

        self.atmospheres = []

73

        self.file_names = file_names

74

        for name in file_names:

75

            spectrum = Spectrum(name, fast_load=True)

76

            self.spectra.append(spectrum)

77

        self.nspectra = len(self.spectra)

78

        # prepare parameters

79

        p = np.array([0.015, 260, 3, -1, *np.ones(self.nspectra)])

80

        self.A1_first_index = 4

81

        fixed = [False] * p.size

82

        # fixed[0] = True

83

        fixed[3] = True

84

        fixed[self.A1_first_index] = True

85

        if fixed_A1s:

86

            for ip in range(self.A1_first_index, len(fixed)):

87

                fixed[ip] = True

88

        input_labels = ["VAOD", "ozone", "PWV", "reso"] + [f"A1_{k}" for k in range(self.nspectra)]

89

        axis_names = ["VAOD", "ozone", "PWV", "reso"] + ["$A_1^{(" + str(k) + ")}$" for k in range(self.nspectra)]

90

        bounds = [(0, 0.01), (100, 700), (0, 10), (0.1, 100)] + [(1e-3, 2)] * self.nspectra

91

        params = FitParameters(p, input_labels=input_labels, axis_names=axis_names, bounds=bounds, fixed=fixed)

92

        FitWorkspace.__init__(self, params, output_file_name, verbose, plot, live_fit)

93

        self.output_file_name = output_file_name

94

        self.bin_widths = bin_width

95

        self.spectrum_lambdas = [self.spectra[k].lambdas for k in range(self.nspectra)]

96

        self.spectrum_data = [self.spectra[k].data for k in range(self.nspectra)]

97

        self.spectrum_err = [self.spectra[k].err for k in range(self.nspectra)]

98

        self.spectrum_data_cov = [self.spectra[k].cov_matrix for k in range(self.nspectra)]

99

        for k, name in enumerate(file_names):

100

            atmgrid_file_name = name.replace("sim", "reduc").replace("spectrum.fits", "atmsim.fits")

101

            if os.path.isfile(atmgrid_file_name):

102

                self.atmospheres.append(AtmosphereGrid(spectrum_filename=name, atmgrid_filename=atmgrid_file_name))

103

            else:

104

                self.my_logger.warning(f"\n\tNo atmosphere grid {atmgrid_file_name}, the fit will be slower...")

105

                self.atmospheres.append(Atmosphere(self.spectra[k].airmass, self.spectra[k].pressure, self.spectra[k].temperature,

106

                                                   lambda_min=np.min(self.spectrum_lambdas),

107

                                                   lambda_max=np.max(self.spectrum_lambdas)+1))

108

        self.lambdas = np.empty(1)

109

        self.lambdas_bin_edges = np.empty(1)

110

        self.ref_spectrum_cube = []

111

        self.random_A1s = None

112

        self._prepare_data()

113

        for atmosphere in self.atmospheres:

114

            if isinstance(atmosphere, AtmosphereGrid):

115

                self.params.bounds[0] = (min(self.atmospheres[0].AER_Points), max(self.atmospheres[0].AER_Points))

116

                self.params.bounds[1] = (min(self.atmospheres[0].OZ_Points), max(self.atmospheres[0].OZ_Points))

117

                self.params.bounds[2] = (min(self.atmospheres[0].PWV_Points), max(self.atmospheres[0].PWV_Points))

118

                break

119

        self.amplitude_truth = None

120

        self.lambdas_truth = None

121

        self.atmosphere = Atmosphere(airmass=1,

122

                                     pressure=float(np.mean([self.spectra[k].header["OUTPRESS"]

123

                                                             for k in range(self.nspectra)])),

124

                                     temperature=float(np.mean([self.spectra[k].header["OUTTEMP"]

125

                                                                for k in range(self.nspectra)])),

126

                                     lambda_min=np.min(self.lambdas_bin_edges),

127

                                     lambda_max=np.max(self.lambdas_bin_edges))

128

        self.true_instrumental_transmission = None

129

        self.true_atmospheric_transmission = None

130

        self.true_A1s = None

131

        self.get_truth()

132

        if inject_random_A1s:

133

            self.inject_random_A1s()

134

135

        # design matrix

136

        self.M = np.zeros((self.nspectra, self.lambdas.size, self.lambdas.size))

137

        self.M_dot_W_dot_M = np.zeros((self.lambdas.size, self.lambdas.size))

138

139

        # prepare results

140

        self.amplitude_params = np.ones(self.lambdas.size)

141

        self.amplitude_params_err = np.zeros(self.lambdas.size)

142

        self.amplitude_cov_matrix = np.zeros((self.lambdas.size, self.lambdas.size))

143

144

        # regularisation

145

        self.amplitude_priors_method = "noprior"

146

        self.reg = parameters.PSF_FIT_REG_PARAM * self.bin_widths

147

        if self.amplitude_priors_method == "spectrum":

148

            self.amplitude_priors = np.copy(self.true_instrumental_transmission)

149

            self.amplitude_priors_cov_matrix = np.eye(self.lambdas[0].size)  # np.diag(np.ones_like(self.lambdas))

150

            self.U = np.diag([1 / np.sqrt(self.amplitude_priors_cov_matrix[i, i]) for i in range(self.lambdas[0].size)])

151

            L = np.diag(-2 * np.ones(self.lambdas[0].size)) + np.diag(np.ones(self.lambdas[0].size), -1)[:-1, :-1] \

152

                + np.diag(np.ones(self.lambdas[0].size), 1)[:-1, :-1]

153

            L[0, 0] = -1

154

            L[-1, -1] = -1

155

            self.L = L.astype(float)

156

            self.Q = L.T @ np.linalg.inv(self.amplitude_priors_cov_matrix) @ L

157

            self.Q_dot_A0 = self.Q @ self.amplitude_priors

158

159

    @property

160

    def A1s(self):

161

        return self.params.p[self.A1_first_index:]

162

163

    def _prepare_data(self):

164

        # rebin wavelengths

165

        if self.bin_widths > 0:

166

            if self.nspectra > 1:

167

                lambdas_bin_edges = np.arange(int(np.min(np.concatenate(list(self.spectrum_lambdas)))),

168

                                              int(np.max(np.concatenate(list(self.spectrum_lambdas)))) + 1,

169

                                              self.bin_widths)

170

            else:

171

                lambdas_bin_edges = np.arange(int(np.min(self.spectrum_lambdas[0])),

172

                                              int(np.max(self.spectrum_lambdas[0])) + 1,

173

                                              self.bin_widths)

174

            lbdas = []

175

            for i in range(1, lambdas_bin_edges.size):

176

                lbdas.append(0.5 * (0*lambdas_bin_edges[i] + 2*lambdas_bin_edges[i - 1]))  # lambda bin value on left

177

            self.lambdas = []

178

            for k in range(self.nspectra):

179

                self.lambdas.append(np.asarray(lbdas))

180

            self.lambdas = np.asarray(self.lambdas)

181

        else:

182

            for k in range(1, len(self.spectrum_lambdas)):

183

                if self.spectrum_lambdas[k].size != self.spectrum_lambdas[0].size or \

184

                        not np.all(np.isclose(self.spectrum_lambdas[k], self.spectrum_lambdas[0])):

185

                    raise ValueError("\nIf you don't rebin your spectra, "

186

                                     "they must share the same wavelength arrays (in length and values).")

187

            self.lambdas = np.copy(self.spectrum_lambdas)

188

            dlbda = self.lambdas[0, -1] - self.lambdas[0, -2]

189

            lambdas_bin_edges = list(self.lambdas[0]) + [self.lambdas[0, -1] + dlbda]

190

        self.lambdas_bin_edges = lambdas_bin_edges

191

        # mask

192

        lambdas_to_mask = [np.arange(300, 355, self.bin_widths)]

193

        for line in [HALPHA, HBETA, HGAMMA, HDELTA, O2_1, O2_2, O2B]:

194

            width = line.width_bounds[1]

195

            lambdas_to_mask += [np.arange(line.wavelength - width, line.wavelength + width, self.bin_widths)]

196

        lambdas_to_mask = np.concatenate(lambdas_to_mask).ravel()

197

        lambdas_to_mask_indices = []

198

        for k in range(self.nspectra):

199

            lambdas_to_mask_indices.append(np.asarray([np.argmin(np.abs(self.lambdas[k] - lambdas_to_mask[i]))

200

                                                       for i in range(lambdas_to_mask.size)]))

201

        # rebin atmosphere

202

        if self.bin_widths > 0 and isinstance(self.atmospheres[0], AtmosphereGrid):

203

            self.atmosphere_lambda_bins = []

204

            for i in range(0, lambdas_bin_edges.size):

205

                self.atmosphere_lambda_bins.append([])

206

                for j in range(0, self.atmospheres[0].lambdas.size):

207

                    if self.atmospheres[0].lambdas[j] >= lambdas_bin_edges[i]:

208

                        self.atmosphere_lambda_bins[-1].append(j)

209

                    if i < lambdas_bin_edges.size - 1 and self.atmospheres[0].lambdas[j] >= lambdas_bin_edges[i + 1]:

210

                        self.atmosphere_lambda_bins[-1] = np.array(self.atmosphere_lambda_bins[-1])

211

                        break

212

            self.atmosphere_lambda_bins = np.array(self.atmosphere_lambda_bins, dtype=object)

213

            self.atmosphere_lambda_step = np.gradient(self.atmospheres[0].lambdas)[0]

214

        # rescale data lambdas

215

        # D2CCD = np.median([self.spectra[k].header["D2CCD"] for k in range(self.nspectra)])

216

        # for k in range(self.nspectra):

217

        #     self.spectra[k].disperser.D = self.spectra[k].header["D2CCD"]

218

        #     dist = self.spectra[k].disperser.grating_lambda_to_pixel(self.spectra[k].lambdas, x0=self.spectra[k].x0)

219

        #     self.spectra[k].disperser.D = D2CCD

220

        #     self.spectra[k].lambdas = self.spectra[k].disperser.grating_pixel_to_lambda(dist, x0=self.spectra[k].x0)

221

        # rebin data

222

        self.data = np.empty(self.nspectra, dtype=object)

223

        if self.bin_widths > 0:

224

            for k in range(self.nspectra):

225

                data_func = interp1d(self.spectra[k].lambdas, self.spectra[k].data,

226

                                     kind="cubic", fill_value="extrapolate", bounds_error=None)

227

                # lambdas_truth = np.fromstring(self.spectra[k].header['LBDAS_T'][1:-1], sep=' ')

228

                # amplitude_truth = np.fromstring(self.spectra[k].header['AMPLIS_T'][1:-1], sep=' ', dtype=float)

229

                # data_func = interp1d(lambdas_truth, amplitude_truth,

230

                #                      kind="cubic", fill_value="extrapolate", bounds_error=None)

231

                data = []

232

                for i in range(1, lambdas_bin_edges.size):

233

                    data.append(quad(data_func, lambdas_bin_edges[i - 1], lambdas_bin_edges[i])[0] / self.bin_widths)

234

                self.data[k] = np.copy(data)

235

                # if parameters.DEBUG:

236

                #     if "LBDAS_T" in self.spectra[k].header:

237

                #         lambdas_truth = np.fromstring(self.spectra[k].header['LBDAS_T'][1:-1], sep=' ')

238

                #         amplitude_truth = np.fromstring(self.spectra[k].header['AMPLIS_T'][1:-1],sep=' ',dtype=float)

239

                #         plt.plot(lambdas_truth, amplitude_truth, label="truth")  # -amplitude_truth)

240

                #     plt.plot(self.lambdas, self.data_cube[-1], label="binned data")  # -amplitude_truth)

241

                #     plt.plot(self.spectra[k].lambdas, self.spectra[k].data, label="raw data")  # -amplitude_truth)

242

                #     # plt.title(self.spectra[k].filename)

243

                #     # plt.xlim(480,700)

244

                #     plt.grid()

245

                #     plt.legend()

246

                #     plt.show()

247

        else:

248

            for k in range(self.nspectra):

249

                self.data[k] = np.copy(self.spectrum_data[k])

250

        # rebin reference star

251

        self.ref_spectrum_cube = []

252

        if self.bin_widths > 0:

253

            for k in range(self.nspectra):

254

                data_func = interp1d(self.spectra[k].target.wavelengths[0], self.spectra[k].target.spectra[0],

255

                                     kind="cubic", fill_value="extrapolate", bounds_error=None)

256

                data = []

257

                for i in range(1, lambdas_bin_edges.size):

258

                    data.append(quad(data_func, lambdas_bin_edges[i - 1], lambdas_bin_edges[i])[0] / self.bin_widths)

259

                self.ref_spectrum_cube.append(np.copy(data))

260

        else:

261

            for k in range(self.nspectra):

262

                ref = interp1d(self.spectra[k].target.wavelengths[0], self.spectra[k].target.spectra[0],

263

                               kind="cubic", fill_value="extrapolate", bounds_error=None)(self.lambdas[k])

264

                self.ref_spectrum_cube.append(np.copy(ref))

265

        self.ref_spectrum_cube = np.asarray(self.ref_spectrum_cube)

266

        # rebin errors

267

        self.err = np.empty(self.nspectra, dtype=object)

268

        if self.bin_widths > 0:

269

            for k in range(self.nspectra):

270

                err_func = interp1d(self.spectra[k].lambdas, self.spectra[k].err ** 2,

271

                                    kind="cubic", fill_value="extrapolate", bounds_error=False)

272

                err = []

273

                for i in range(1, lambdas_bin_edges.size):

274

                    if i in lambdas_to_mask_indices[k]:

275

                        err.append(np.nan)

276

                    else:

277

                        err.append(np.sqrt(np.abs(quad(err_func, lambdas_bin_edges[i - 1], lambdas_bin_edges[i])[0])

278

                                           / self.bin_widths))

279

                self.err[k] = np.copy(err)

280

        else:

281

            for k in range(self.nspectra):

282

                self.err[k] = np.copy(self.spectrum_err[k])

283

        if parameters.DEBUG:

284

            for k in range(self.nspectra):

285

                plt.errorbar(self.lambdas[k], self.data[k], self.err[k], label=f"spectrum {k}")

286

                plt.ylim(0, 1.2 * np.max(self.data[k]))

287

            plt.grid()

288

            # plt.legend()

289

            plt.show()

290

        # rebin W matrices

291

        # import time

292

        # start = time.time()

293

        self.data_cov = np.empty(self.nspectra, dtype=object)

294

        self.W = np.empty(self.nspectra, dtype=object)

295

        if self.bin_widths > 0:

296

            lmins = []

297

            lmaxs = []

298

            for k in range(self.nspectra):

299

                lmins.append([])

300

                lmaxs.append([])

301

                for i in range(self.lambdas[k].size):

302

                    lmins[-1].append(max(0, int(np.argmin(np.abs(self.spectrum_lambdas[k] - lambdas_bin_edges[i])))))

303

                    lmaxs[-1].append(min(self.spectrum_data_cov[k].shape[0] - 1,

304

                                         np.argmin(np.abs(self.spectrum_lambdas[k] - lambdas_bin_edges[i + 1]))))

305

            for k in range(self.nspectra):

306

                cov = np.zeros((self.lambdas[k].size, self.lambdas[k].size))

307

                for i in range(cov.shape[0]):

308

                    # imin = max(0, int(np.argmin(np.abs(self.spectrum_lambdas[k] - lambdas_bin_edges[i]))))

309

                    # imax = min(self.spectrum_data_cov[k].shape[0] - 1,

310

                    #           np.argmin(np.abs(self.spectrum_lambdas[k] - lambdas_bin_edges[i + 1])))

311

                    imin = lmins[k][i]

312

                    imax = lmaxs[k][i]

313

                    if imin == imax:

314

                        cov[i, i] = (i + 1) * 1e10

315

                        continue

316

                    if i in lambdas_to_mask_indices[k]:

317

                        cov[i, i] = (i + 1e10)

318

                        continue

319

                    for j in range(i, cov.shape[1]):

320

                        # jmin = max(0, int(np.argmin(np.abs(self.spectrum_lambdas[k] - lambdas_bin_edges[j]))))

321

                        # jmax = min(self.spectrum_data_cov[k].shape[0] - 1,

322

                        #           np.argmin(np.abs(self.spectrum_lambdas[k] - lambdas_bin_edges[j + 1])))

323

                        jmin = lmins[k][j]

324

                        jmax = lmaxs[k][j]

325

                        # if imin == imax:

326

                        #     cov[i, i] = (i + 1) * 1e10

327

                        # elif jmin == jmax:

328

                        #     cov[j, j] = (j + 1) * 1e10

329

                        # else:

330

                        if jmin == jmax:

331

                            cov[j, j] = (j + 1) * 1e10

332

                        else:

333

                            if j in lambdas_to_mask_indices[k]:

334

                                cov[j, j] = (j + 1e10)

335

                            else:

336

                                mean = np.mean(self.spectrum_data_cov[k][imin:imax, jmin:jmax])

337

                                cov[i, j] = mean

338

                                cov[j, i] = mean

339

                self.data_cov[k] = np.copy(cov)

340

            # self.data_cov = np.zeros(self.nspectra * np.array(self.data_cov_cube[0].shape))

341

            # for k in range(self.nspectra):

342

            #     self.data_cov[k * self.lambdas[k].size:(k + 1) * self.lambdas[k].size,

343

            #     k * self.lambdas[k].size:(k + 1) * self.lambdas[k].size] = \

344

            #         self.data_cov_cube[k]

345

            # self.data_cov = self.data_cov_cube

346

            # print("fill data_cov_cube", time.time() - start)

347

            # start = time.time()

348

            for k in range(self.nspectra):

349

                try:

350

                    L = np.linalg.inv(np.linalg.cholesky(self.data_cov[k]))

351

                    invcov_matrix = L.T @ L

352

                except np.linalg.LinAlgError:

353

                    invcov_matrix = np.linalg.inv(self.data_cov[k])

354

                self.W[k] = invcov_matrix

355

            # self.data_invcov = np.zeros(self.nspectra * np.array(self.data_cov_cube[0].shape))

356

            # for k in range(self.nspectra):

357

            #     self.data_invcov[k * self.lambdas[k].size:(k + 1) * self.lambdas[k].size,

358

            #     k * self.lambdas[k].size:(k + 1) * self.lambdas[k].size] = \

359

            #         self.data_invcov_cube[k]

360

            # self.data_invcov = self.data_invcov_cube

361

            # print("inv data_cov_cube", time.time() - start)

362

            # start = time.time()

363

        else:

364

            self.W = np.empty(self.nspectra, dtype=np.object)

365

            for k in range(self.nspectra):

366

                try:

367

                    L = np.linalg.inv(np.linalg.cholesky(self.spectrum_data_cov[k]))

368

                    invcov_matrix = L.T @ L

369

                except np.linalg.LinAlgError:

370

                    invcov_matrix = np.linalg.inv(self.spectrum_data_cov[k])

371

                invcov_matrix[lambdas_to_mask_indices[k], :] = 0

372

                invcov_matrix[:, lambdas_to_mask_indices[k]] = 0

373

                self.W[k] = invcov_matrix

374

375

    def inject_random_A1s(self):  # pragma: no cover

376

        random_A1s = np.random.uniform(0.5, 1, size=self.nspectra)

377

        for k in range(self.nspectra):

378

            self.data[k] *= random_A1s[k]

379

            self.err[k] *= random_A1s[k]

380

            self.data_cov[k] *= random_A1s[k] ** 2

381

            self.W[k] /= random_A1s[k] ** 2

382

        if self.true_A1s is not None:

383

            self.true_A1s *= random_A1s

384

385

    def get_truth(self):

386

        """Load the truth parameters (if provided) from the file header.

387

388

"""

389

        if 'A1_T' in list(self.spectra[0].header.keys()):

390

            ozone_truth = self.spectrum.header['OZONE_T']

391

            pwv_truth = self.spectrum.header['PWV_T']

392

            aerosols_truth = self.spectrum.header['VAOD_T']

393

            self.truth = (aerosols_truth, ozone_truth, pwv_truth)

394

            self.true_atmospheric_transmission = []

395

            tatm = self.atmosphere.simulate(ozone=ozone_truth, pwv=pwv_truth, aerosols=aerosols_truth)

396

            if self.bin_widths > 0:

397

                for i in range(1, self.lambdas_bin_edges.size):

398

                    self.true_atmospheric_transmission.append(quad(tatm, self.lambdas_bin_edges[i - 1],

399

                                                                   self.lambdas_bin_edges[i])[0] / self.bin_widths)

400

            else:

401

                self.true_atmospheric_transmission = tatm(self.lambdas[0])

402

            self.true_atmospheric_transmission = np.array(self.true_atmospheric_transmission)

403

            self.true_A1s = np.array([self.spectra[k].header["A1_T"] for k in range(self.nspectra)], dtype=float)

404

        else:

405

            self.truth = None

406

        self.true_instrumental_transmission = []

407

        tinst = lambda lbda: self.spectrum.disperser.transmission(lbda) * self.spectrum.throughput.transmission(lbda)

408

        if self.bin_widths > 0:

409

            for i in range(1, self.lambdas_bin_edges.size):

410

                self.true_instrumental_transmission.append(quad(tinst, self.lambdas_bin_edges[i - 1],

411

                                                                self.lambdas_bin_edges[i])[0] / self.bin_widths)

412

        else:

413

            self.true_instrumental_transmission = tinst(self.lambdas[0])

414

        self.true_instrumental_transmission = np.array(self.true_instrumental_transmission)

415

416

    def simulate(self, aerosols, ozone, pwv, reso, *A1s):

417

        """Interface method to simulate multiple spectra with a single atmosphere.

418

419

        Parameters

420

        ----------

421

        aerosols: float

422

            Vertical Aerosols Optical Depth quantity for Libradtran (no units).

423

        ozone: float

424

            Ozone parameter for Libradtran (in db).

425

        pwv: float

426

            Precipitable Water Vapor quantity for Libradtran (in mm).

427

        reso: float

428

            Width of the gaussian kernel to smooth the spectra (if <0: no convolution).

429

430

        Returns

431

        -------

432

        lambdas: array_like

433

            Array of wavelengths (1D).

434

        model: array_like

435

            2D array of the spectrogram simulation.

436

        model_err: array_like

437

            2D array of the spectrogram simulation uncertainty.

438

439

        Examples

440

        --------

441

        >>> from spectractor.config import load_config

442

        >>> load_config("./config/ctio.ini")

443

        >>> file_names = ["./tests/data/reduc_20170530_134_spectrum.fits"]

444

        >>> w = MultiSpectraFitWorkspace("./outputs/test", file_names, bin_width=5, verbose=True)

445

        >>> lambdas, model, model_err = w.simulate(*w.params.p)

446

        >>> assert np.sum(model) > 0

447

        >>> assert np.all(lambdas == w.lambdas)

448

        >>> assert np.sum(w.amplitude_params) > 0

449

450

"""

451

        # linear regression for the instrumental transmission parameters T

452

        # first: force the grey terms to have an average of 1

453

        A1s = np.array(A1s)

454

        if A1s.size > 1:

455

            m = 1

456

            A1s[0] = m * A1s.size - np.sum(A1s[1:])

457

            self.params.p[self.A1_first_index] = A1s[0]

458

        # Matrix M filling: hereafter a fast integration is used

459

        M = []

460

        for k in range(self.nspectra):

461

            atm = []

462

            a = self.atmospheres[k].simulate(aerosols=aerosols, ozone=ozone, pwv=pwv)

463

            lbdas = self.atmospheres[k].lambdas

464

            for i in range(1, self.lambdas_bin_edges.size):

465

                if isinstance(self.atmospheres[k], AtmosphereGrid):

466

                    delta = self.atmosphere_lambda_bins[i][-1] - self.atmosphere_lambda_bins[i][0]

467

                    if delta > 0:

468

                        atm.append(

469

                            np.trapz(a(lbdas[self.atmosphere_lambda_bins[i]]), dx=self.atmosphere_lambda_step) / delta)

470

                    else:

471

                        atm.append(1)

472

                else:

473

                    delta = self.lambdas_bin_edges[i] - self.lambdas_bin_edges[i-1]

474

                    if delta > 0:

475

                        atm.append(quad(a, self.lambdas_bin_edges[i-1], self.lambdas_bin_edges[i])[0] / delta)

476

                    else:

477

                        atm.append(1)

478

            if reso > 0:

479

                M.append(A1s[k] * np.diag(fftconvolve_gaussian(self.ref_spectrum_cube[k] * np.array(atm), reso)))

480

            else:

481

                M.append(A1s[k] * np.diag(self.ref_spectrum_cube[k] * np.array(atm)))

482

        # hereafter: no binning but gives unbiased result on extracted spectra from simulations and truth spectra

483

        # if self.reso > 0:

484

        #     M = np.array([A1s[k] * np.diag(fftconvolve_gaussian(self.ref_spectrum_cube[k] *

485

        #                                    self.atmospheres[k].simulate(aerosols, ozone, pwv)(self.lambdas[k]), reso))

486

        #                   for k in range(self.nspectra)])

487

        # else:

488

        #     M = np.array([A1s[k] * np.diag(self.ref_spectrum_cube[k] *

489

        #                                    self.atmospheres[k].simulate(aerosols, ozone, pwv)(self.lambdas[k]))

490

        #                   for k in range(self.nspectra)])

491

        # print("compute M", time.time() - start)

492

        # start = time.time()

493

        # for k in range(self.nspectra):

494

        #     plt.plot(self.atmospheres[k].lambdas, [M[k][i,i] for i in range(self.atmospheres[k].lambdas.size)])

495

        #     # plt.plot(self.lambdas, self.ref_spectrum_cube[k], linestyle="--")

496

        # plt.grid()

497

        # plt.title(f"reso={reso:.3f}")

498

        # plt.show()

499

        # Matrix W filling: if spectra are not independent, use these lines with einstein summations:

500

        # W = np.zeros((self.nspectra, self.nspectra, self.lambdas.size, self.lambdas.size))

501

        # for k in range(self.nspectra):

502

        #     W[k, k, ...] = self.data_invcov[k]

503

        # W_dot_M = np.einsum('lkji,kjh->lih', W, M)

504

        # M_dot_W_dot_M = np.einsum('lkj,lki->ij', M, W_dot_M)

505

        # M_dot_W_dot_M = np.zeros_like(M_dot_W_dot_M)

506

        # otherwise, this is much faster:

507

        M_dot_W_dot_M = np.sum([M[k].T @ self.W[k] @ M[k] for k in range(self.nspectra)], axis=0)

508

        M_dot_W_dot_D = np.sum([M[k].T @ self.W[k] @ self.data[k] for k in range(self.nspectra)], axis=0)

509

        if self.amplitude_priors_method != "spectrum":

510

            for i in range(self.lambdas[0].size):

511

                if np.sum(M_dot_W_dot_M[i]) == 0:

512

                    M_dot_W_dot_M[i, i] = 1e-10 * np.mean(M_dot_W_dot_M) * np.random.random()

513

            try:

514

                L = np.linalg.inv(np.linalg.cholesky(M_dot_W_dot_M))

515

                cov_matrix = L.T @ L

516

            except np.linalg.LinAlgError:

517

                cov_matrix = np.linalg.inv(M_dot_W_dot_M)

518

            amplitude_params = cov_matrix @ M_dot_W_dot_D

519

        else:

520

            M_dot_W_dot_M_plus_Q = M_dot_W_dot_M + self.reg * self.Q

521

            try:

522

                L = np.linalg.inv(np.linalg.cholesky(M_dot_W_dot_M_plus_Q))

523

                cov_matrix = L.T @ L

524

            except np.linalg.LinAlgError:

525

                cov_matrix = np.linalg.inv(M_dot_W_dot_M_plus_Q)

526

            amplitude_params = cov_matrix @ (M_dot_W_dot_D + self.reg * self.Q_dot_A0)

527

528

        self.M = M

529

        self.M_dot_W_dot_M = M_dot_W_dot_M

530

        self.M_dot_W_dot_D = M_dot_W_dot_D

531

        model_cube = []

532

        model_err_cube = []

533

        for k in range(self.nspectra):

534

            model_cube.append(M[k] @ amplitude_params)

535

            model_err_cube.append(np.zeros_like(model_cube[-1]))

536

        self.model = np.asarray(model_cube)

537

        self.model_err = np.asarray(model_err_cube)

538

        self.amplitude_params = np.copy(amplitude_params)

539

        self.amplitude_params_err = np.array([np.sqrt(cov_matrix[i, i])

540

                                              if cov_matrix[i, i] > 0 else 0 for i in range(self.lambdas[0].size)])

541

        self.amplitude_cov_matrix = np.copy(cov_matrix)

542

        # print("algebra", time.time() - start)

543

        # start = time.time()

544

        return self.lambdas, self.model, self.model_err

545

546

    def plot_fit(self):

547

        """Plot the fit result.

548

549

        Examples

550

        --------

551

        >>> from spectractor.config import load_config

552

        >>> load_config("./config/ctio.ini")

553

        >>> file_names = 3 * ["./tests/data/reduc_20170530_134_spectrum.fits"]

554

        >>> w = MultiSpectraFitWorkspace("./outputs/test", file_names, bin_width=5, verbose=True)

555

        >>> w.simulate(*w.params.p)  #doctest: +ELLIPSIS

556

        (array(...

557

        >>> w.plot_fit()

558

"""

559

        cmap_bwr = copy.copy(cm.get_cmap('bwr'))

560

        cmap_bwr.set_bad(color='lightgrey')

561

        cmap_viridis = copy.copy(cm.get_cmap('viridis'))

562

        cmap_viridis.set_bad(color='lightgrey')

563

564

        data = copy.deepcopy(self.data)

565

        for k in range(self.nspectra):

566

            data[k][np.isnan(data[k]/self.err[k])] = np.nan

567

        if len(self.outliers) > 0:

568

            bad_indices = self.get_bad_indices()

569

            for k in range(self.nspectra):

570

                data[k][bad_indices[k]] = np.nan

571

                data[k] = np.ma.masked_invalid(data[k])

572

        data = np.array([data[k] for k in range(self.nspectra)], dtype=float)

573

        model = np.array([self.model[k] for k in range(self.nspectra)], dtype=float)

574

        err = np.array([self.err[k] for k in range(self.nspectra)], dtype=float)

575

        gs_kw = dict(width_ratios=[3, 0.13], height_ratios=[1, 1, 1])

576

        fig, ax = plt.subplots(nrows=3, ncols=2, figsize=(7, 6), gridspec_kw=gs_kw)

577

        # aerosols, ozone, pwv, reso, *A1s = self.p

578

        # plt.suptitle(f'VAOD={aerosols:.3f}, ozone={ozone:.0f}db, PWV={pwv:.2f}mm, reso={reso:.2f}', y=0.995)

579

        norm = np.nanmax(data)

580

        y = np.arange(0, self.nspectra+1).astype(int) - 0.5

581

        xx, yy = np.meshgrid(self.lambdas[0], y)

582

        ylbda = -0.45 * np.ones_like(self.lambdas[0][1:-1])

583

        # model

584

        im = ax[1, 0].pcolormesh(xx, yy, model / norm, vmin=0, vmax=1, cmap=cmap_viridis, shading="auto")

585

        plt.colorbar(im, cax=ax[1, 1], label='1/max(data)', format="%.1f")

586

        ax[1, 0].set_title("Model", fontsize=12, color='white', x=0.91, y=0.76)

587

        ax[1, 0].grid(color='silver', ls='solid')

588

        ax[1, 0].scatter(self.lambdas[0][1:-1], ylbda, cmap=from_lambda_to_colormap(self.lambdas[0][1:-1]),

589

                         edgecolors='None', c=self.lambdas[0][1:-1], label='', marker='o', s=20)

590

        # data

591

        im = ax[0, 0].pcolormesh(xx, yy, data / norm, vmin=0, vmax=1, cmap=cmap_viridis, shading="auto")

592

        plt.colorbar(im, cax=ax[0, 1], label='1/max(data)', format="%.1f")

593

        ax[0, 0].set_title("Data", fontsize=12, color='white', x=0.91, y=0.76)

594

        ax[0, 0].grid(color='silver', ls='solid')

595

        ax[0, 0].scatter(self.lambdas[0][1:-1], ylbda, cmap=from_lambda_to_colormap(self.lambdas[0][1:-1]),

596

                         edgecolors='None', c=self.lambdas[0][1:-1], label='', marker='o', s=20)

597

        # residuals

598

        residuals = (data - model)

599

        norm = err

600

        residuals /= norm

601

        std = float(np.nanstd(residuals))

602

        im = ax[2, 0].pcolormesh(xx, yy, residuals, vmin=-3 * std, vmax=3 * std, cmap=cmap_bwr, shading="auto")

603

        plt.colorbar(im, cax=ax[2, 1], label='(Data-Model)/Err', format="%.0f")

604

        # ax[2, 0].set_title('(Data-Model)/Err', fontsize=10, color='black', x=0.84, y=0.76)

605

        ax[2, 0].grid(color='silver', ls='solid')

606

        ax[2, 0].scatter(self.lambdas[0][1:-1], ylbda, cmap=from_lambda_to_colormap(self.lambdas[0][1:-1]),

607

                         edgecolors='None', c=self.lambdas[0][1:-1], label='', marker='o', s=10*self.nspectra)

608

        ax[2, 0].text(0.05, 0.8, f'mean={np.nanmean(residuals):.3f}\nstd={np.nanstd(residuals):.3f}',

609

                      horizontalalignment='left', verticalalignment='bottom',

610

                      color='black', transform=ax[2, 0].transAxes)

611

        ax[2, 0].set_xlabel(r"$\lambda$ [nm]")

612

        for i in range(3):

613

            ax[i, 0].set_xlim(self.lambdas[0, 0], self.lambdas[0, -1])

614

            ax[i, 0].set_ylim(-0.5, self.nspectra-0.5)

615

            ax[i, 0].yaxis.set_major_locator(MaxNLocator(integer=True))

616

            ax[i, 0].set_ylabel("Spectrum index")

617

            ax[i, 1].get_yaxis().set_label_coords(2.6, 0.5)

618

            ax[i, 0].get_yaxis().set_label_coords(-0.06, 0.5)

619

        fig.tight_layout()

620

        if parameters.SAVE:

621

            fig.savefig(self.output_file_name + '_bestfit.pdf', dpi=100, bbox_inches='tight')

622

        if self.live_fit:  # pragma: no cover

623

            plt.draw()

624

            plt.pause(1e-8)

625

            plt.close()

626

        else:  # pragma: no cover

627

            if parameters.DISPLAY and self.verbose:

628

                plt.show()

629

630

    def plot_transmissions(self):

631

        """Plot the fit result for transmissions.

632

633

        Examples

634

        --------

635

        >>> from spectractor.config import load_config

636

        >>> load_config("./config/ctio.ini")

637

        >>> file_names = ["./tests/data/sim_20170530_134_spectrum.fits"]

638

        >>> w = MultiSpectraFitWorkspace("./outputs/test", file_names, bin_width=5, verbose=True)

639

        >>> w.plot_transmissions()

640

"""

641

        gs_kw = dict(width_ratios=[1, 1], height_ratios=[1, 0.15])

642

        fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(9, 6), gridspec_kw=gs_kw, sharex="all")

643

        aerosols, ozone, pwv, reso, *A1s = self.params.p

644

        plt.suptitle(f'VAOD={aerosols:.3f}, ozone={ozone:.0f}db, PWV={pwv:.2f}mm', y=1)

645

        masked = self.amplitude_params_err > 1e6

646

        transmission = np.copy(self.amplitude_params)

647

        transmission_err = np.copy(self.amplitude_params_err)

648

        transmission[masked] = np.nan

649

        transmission_err[masked] = np.nan

650

        ax[0, 0].errorbar(self.lambdas[0], transmission, yerr=transmission_err,

651

                          label=r'$T_{\mathrm{inst}} * \left\langle A_1 \right\rangle$', fmt='k.')  # , markersize=0.1)

652

        ax[0, 0].set_ylabel(r'Instrumental transmission')

653

        ax[0, 0].set_xlim(self.lambdas[0][0], self.lambdas[0][-1])

654

        ax[0, 0].set_ylim(0, 1.1 * np.nanmax(transmission))

655

        ax[0, 0].grid(True)

656

        ax[0, 0].set_xlabel(r'$\lambda$ [nm]')

657

        if self.true_instrumental_transmission is not None:

658

            ax[0, 0].plot(self.lambdas[0], self.true_instrumental_transmission, "g-",

659

                          label=r'true $T_{\mathrm{inst}}* \left\langle A_1 \right\rangle$')

660

            ax[1, 0].set_xlabel(r'$\lambda$ [nm]')

661

            ax[1, 0].grid(True)

662

            ax[1, 0].set_ylabel(r'(Data-Truth)/Err')

663

            norm = transmission_err

664

            residuals = (self.amplitude_params - self.true_instrumental_transmission) / norm

665

            residuals[masked] = np.nan

666

            ax[1, 0].errorbar(self.lambdas[0], residuals, yerr=transmission_err / norm,

667

                              label=r'$T_{\mathrm{inst}}$', fmt='k.')  # , markersize=0.1)

668

            ax[1, 0].set_ylim(-1.1 * np.nanmax(np.abs(residuals)), 1.1 * np.nanmax(np.abs(residuals)))

669

        else:

670

            ax[1, 0].remove()

671

        ax[0, 0].legend()

672

673

        tatm = self.atmosphere.simulate(ozone=ozone, pwv=pwv, aerosols=aerosols)

674

        tatm_binned = []

675

        for i in range(1, self.lambdas_bin_edges.size):

676

            tatm_binned.append(quad(tatm, self.lambdas_bin_edges[i - 1], self.lambdas_bin_edges[i])[0] /

677

                               (self.lambdas_bin_edges[i] - self.lambdas_bin_edges[i - 1]))

678

679

        ax[0, 1].errorbar(self.lambdas[0], tatm_binned,

680

                          label=r'$T_{\mathrm{atm}}$', fmt='k.')  # , markersize=0.1)

681

        ax[0, 1].set_ylabel(r'Atmospheric transmission')

682

        ax[0, 1].set_xlabel(r'$\lambda$ [nm]')

683

        ax[0, 1].set_xlim(self.lambdas[0][0], self.lambdas[0][-1])

684

        ax[0, 1].grid(True)

685

        if self.truth is not None:

686

            ax[0, 1].plot(self.lambdas[0], self.true_atmospheric_transmission, "b-", label=r'true $T_{\mathrm{atm}}$')

687

            ax[1, 1].set_xlabel(r'$\lambda$ [nm]')

688

            ax[1, 1].set_ylabel(r'Data-Truth')

689

            ax[1, 1].grid(True)

690

            residuals = np.asarray(tatm_binned) - self.true_atmospheric_transmission

691

            ax[1, 1].errorbar(self.lambdas[0], residuals, label=r'$T_{\mathrm{inst}}$', fmt='k.')  # , markersize=0.1)

692

            ax[1, 1].set_ylim(-1.1 * np.max(np.abs(residuals)), 1.1 * np.max(np.abs(residuals)))

693

        else:

694

            ax[1, 1].remove()

695

        ax[0, 1].legend()

696

        fig.tight_layout()

697

        if parameters.SAVE:

698

            fig.savefig(self.output_file_name + '_Tinst_best_fit.pdf', dpi=100, bbox_inches='tight')

699

        if self.live_fit:  # pragma: no cover

700

            plt.draw()

701

            plt.pause(1e-8)

702

            plt.close()

703

        else:  # pragma: no cover

704

            if parameters.DISPLAY and self.verbose:

705

                plt.show()

706

707

    def plot_A1s(self):

708

"""

709

        Examples

710

        --------

711

        >>> from spectractor.config import load_config

712

        >>> load_config("./config/ctio.ini")

713

        >>> file_names = ["./tests/data/sim_20170530_134_spectrum.fits"]

714

        >>> w = MultiSpectraFitWorkspace("./outputs/test", file_names, bin_width=5, verbose=True)

715

        >>> w.cov = np.eye(3 + w.nspectra - 1)

716

        >>> w.plot_A1s()

717

718

"""

719

        aerosols, ozone, pwv, reso, *A1s = self.params.p

720

        zs = [self.spectra[k].header["AIRMASS"] for k in range(self.nspectra)]

721

        err = np.sqrt([0] + [self.params.cov[ip, ip] for ip in range(self.A1_first_index, self.params.cov.shape[0])])

722

        spectra_index = np.arange(self.nspectra)

723

        sc = plt.scatter(spectra_index, A1s, c=zs, s=0)

724

        plt.colorbar(sc, label="Airmass")

725

726

        # convert time to a color tuple using the colormap used for scatter

727

        norm = colors.Normalize(vmin=np.min(zs), vmax=np.max(zs), clip=True)

728

        mapper = cm.ScalarMappable(norm=norm, cmap='viridis')

729

        z_color = np.array([(mapper.to_rgba(z)) for z in zs])

730

731

        # loop over each data point to plot

732

        for k, A1, e, color in zip(spectra_index, A1s, err, z_color):

733

            plt.plot(k, A1, 'o', color=color)

734

            plt.errorbar(k, A1, e, lw=1, capsize=3, color=color)

735

736

        if self.true_A1s is not None:

737

            plt.plot(spectra_index, self.true_A1s, 'b-', label="true relative $A_1$'s")

738

739

        plt.axhline(1, color="k", linestyle="--")

740

        plt.axhline(np.mean(A1s), color="b", linestyle="--",

741

                    label=rf"$\left\langle A_1\right\rangle = {np.mean(A1s):.3f}$ (std={np.std(A1s):.3f})")

742

        plt.grid()

743

        plt.ylabel("Relative grey transmissions")

744

        plt.xlabel("Spectrum index")

745

        plt.gca().xaxis.set_major_locator(MaxNLocator(integer=True))

746

        plt.legend()

747

        if parameters.SAVE:

748

            plt.gcf().savefig(self.output_file_name + '_A1s.pdf', dpi=100, bbox_inches='tight')

749

        plt.show()

750

751

    def save_transmissions(self):

752

        aerosols, ozone, pwv, reso, *A1s = self.params.p

753

        tatm = self.atmosphere.simulate(ozone=ozone, pwv=pwv, aerosols=aerosols)

754

        tatm_binned = []

755

        for i in range(1, self.lambdas_bin_edges.size):

756

            tatm_binned.append(quad(tatm, self.lambdas_bin_edges[i - 1], self.lambdas_bin_edges[i])[0] /

757

                               (self.lambdas_bin_edges[i] - self.lambdas_bin_edges[i - 1]))

758

759

        throughput = self.amplitude_params / self.spectrum.disperser.transmission(self.lambdas[0])

760

        throughput_err = self.amplitude_params_err / self.spectrum.disperser.transmission(self.lambdas[0])

761

        if "sim" in self.file_names[0]:

762

            file_name = self.output_file_name + f"_sim_transmissions.txt"

763

        else:

764

            file_name = self.output_file_name + f"_transmissions.txt"

765

        ascii.write([self.lambdas[0], self.amplitude_params, self.amplitude_params_err,

766

                     throughput, throughput_err, tatm_binned], file_name,

767

                    names=["wl", "Tinst", "Tinst_err", "Ttel", "Ttel_err", "Tatm"], overwrite=True)

768

769

    def jacobian(self, params, epsilon, fixed_params=None, model_input=None):

770

        """Generic function to compute the Jacobian matrix of a model, with numerical derivatives.

771

772

        Parameters

773

        ----------

774

        params: array_like

775

            The array of model parameters.

776

        epsilon: array_like

777

            The array of small steps to compute the partial derivatives of the model.

778

        fixed_params: array_like

779

            List of boolean values. If True, the parameter is considered fixed and no derivative are computed.

780

        model_input: array_like, optional

781

            A model input as a list with (x, model, model_err) to avoid an additional call to simulate().

782

783

        Returns

784

        -------

785

        J: np.array

786

            The Jacobian matrix.

787

788

"""

789

        if model_input:

790

            x, model, model_err = model_input

791

        else:

792

            x, model, model_err = self.simulate(*params)

793

        if self.W.dtype == object and self.W[0].ndim == 2:

794

            J = [[] for _ in range(params.size)]

795

        else:

796

            model = model.flatten()

797

            J = np.zeros((params.size, model.size))

798

        for ip, p in enumerate(params):

799

            if self.params.fixed[ip]:

800

                continue

801

            tmp_p = np.copy(params)

802

            if tmp_p[ip] + epsilon[ip] < self.params.bounds[ip][0] or tmp_p[ip] + epsilon[ip] > self.params.bounds[ip][1]:

803

                epsilon[ip] = - epsilon[ip]

804

            tmp_p[ip] += epsilon[ip]

805

            # if "A1_" not in self.input_labels[ip]:

806

            tmp_x, tmp_model, tmp_model_err = self.simulate(*tmp_p)

807

            if self.W.dtype == object and self.W[0].ndim == 2:

808

                for k in range(model.shape[0]):

809

                    J[ip].append((tmp_model[k] - model[k]) / epsilon[ip])

810

            else:

811

                J[ip] = (tmp_model.flatten() - model) / epsilon[ip]

812

        return np.asarray(J)

813

814

815

def run_multispectra_minimisation(fit_workspace, method="newton", verbose=False):

816

    """Interface function to fit spectrum simulation parameters to data.

817

818

    Parameters

819

    ----------

820

    fit_workspace: MultiSpectraFitWorkspace

821

        An instance of the SpectrogramFitWorkspace class.

822

    method: str, optional

823

        Fitting method (default: 'newton').

824

825

    Examples

826

    --------

827

    >>> from spectractor.config import load_config

828

    >>> load_config("./config/ctio.ini")

829

    >>> file_names = 4 * ["./tests/data/reduc_20170530_134_spectrum.fits"]

830

    >>> w = MultiSpectraFitWorkspace("./outputs/test", file_names, bin_width=5, verbose=True, fixed_A1s=False)

831

    >>> parameters.VERBOSE = True

832

    >>> run_multispectra_minimisation(w, method="newton")

833

    >>> assert np.all(np.isclose(w.A1s, 1))

834

835

"""

836

    my_logger = set_logger(__name__)

837

    guess = np.asarray(fit_workspace.params.p)

838

    if method != "newton":

839

        run_minimisation(fit_workspace, method=method)

840

    else:

841

        my_logger.info(f"\n\tStart guess: {guess}\n\twith {fit_workspace.params.input_labels}")

842

        epsilon = 1e-2 * guess

843

        epsilon[epsilon == 0] = 1e-2

844

        if isinstance(fit_workspace.atmospheres[0], AtmosphereGrid):

845

            epsilon = np.array([np.gradient(fit_workspace.atmospheres[0].AER_Points)[0],

846

                                np.gradient(fit_workspace.atmospheres[0].OZ_Points)[0],

847

                                np.gradient(fit_workspace.atmospheres[0].PWV_Points)[0], 0.04]) / 2

848

        else:

849

            epsilon = np.array([0.0001, 1, 0.01])

850

        epsilon = np.array(list(epsilon) + [1e-4] * fit_workspace.A1s.size)

851

852

        run_minimisation_sigma_clipping(fit_workspace, method="newton", epsilon=epsilon, xtol=1e-6,

853

                                        ftol=1 / fit_workspace.data.size, sigma_clip=5, niter_clip=3, verbose=verbose)

854

855

        # w_reg = RegFitWorkspace(fit_workspace, opt_reg=parameters.PSF_FIT_REG_PARAM, verbose=parameters.VERBOSE)

856

        # run_minimisation(w_reg, method="minimize", ftol=1e-4, xtol=1e-2, verbose=parameters.VERBOSE, epsilon=[1e-1],

857

        #                  minimizer_method="Nelder-Mead")

858

        # w_reg.opt_reg = 10 ** w_reg.p[0]

859

        # w_reg.my_logger.info(f"\n\tOptimal regularisation parameter: {w_reg.opt_reg}")

860

        # fit_workspace.reg = np.copy(w_reg.opt_reg)

861

        # fit_workspace.opt_reg = w_reg.opt_reg

862

        # Recompute and save params in class attributes

863

        fit_workspace.simulate(*fit_workspace.params.p)

864

865

        # Renormalize A1s and instrumental transmission

866

        aerosols, ozone, pwv, reso, *A1s = fit_workspace.params.p

867

        mean_A1 = np.mean(A1s)

868

        fit_workspace.amplitude_params /= mean_A1

869

        fit_workspace.amplitude_params_err /= mean_A1

870

        if fit_workspace.true_A1s is not None:

871

            fit_workspace.true_instrumental_transmission *= np.mean(fit_workspace.true_A1s)

872

            fit_workspace.true_A1s /= np.mean(fit_workspace.true_A1s)

873

874

        tinst = np.array(fit_workspace.amplitude_params)

875

        for k in range(fit_workspace.nspectra):

876

            plt.plot(fit_workspace.lambdas[k],

877

                     tinst * np.array([fit_workspace.M[k][i, i] for i in range(fit_workspace.lambdas[k].size)]))

878

            plt.ylim(0, 1.2 * np.max(fit_workspace.data[k]))

879

        plt.grid()

880

        plt.title(f"reso={reso:.3f}")

881

        plt.show()

882

        if fit_workspace.filename != "":

883

            parameters.SAVE = True

884

            fit_workspace.params.plot_correlation_matrix()

885

            write_fitparameter_json(fit_workspace.params.json_filename, fit_workspace.params,

886

                                    extra={"chi2": fit_workspace.costs[-1] / fit_workspace.data.size,

887

                                           "date-obs": fit_workspace.spectrum.date_obs})

888

            fit_workspace.plot_fit()

889

            fit_workspace.plot_transmissions()

890

            fit_workspace.plot_A1s()

891

            fit_workspace.save_transmissions()

892

            parameters.SAVE = False

893

894

895

if __name__ == "__main__":

896

    import doctest

897

898

    doctest.testmod()

LSSTDESC / Spectractor / 4890731301

Source File
Press 'n' to go to next uncovered line, 'b' for previous

LSSTDESC / Spectractor / 4890731301

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous