__init__.py 10.5 KB
Newer Older
1
# -*- coding: utf-8 -*-
2
3
# Copyright (C) 2014 Laboratoire d'Astrophysique de Marseille, AMU
# Copyright (C) 2012, 2014 Centre de données Astrophysiques de Marseille
4
# Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
5
# Author: Yannick Roehlly & Denis Burgarella
6
7

from importlib import import_module
8
9
10
import numpy as np
from scipy import stats
from astropy.table import Column
11
12
13
14
15
16


class AnalysisModule(object):
    """Abstract class, the pCigale analysis modules are based on.
    """

17
    # parameter_list is a dictionary containing all the parameters
18
19
20
21
22
    # used by the module. Each parameter name is associate to a tuple
    # (variable type, description [string], default value). Each module must
    # define its parameter list, unless it does not use any parameter. Using
    # None means that there is no description, unit or default value. If None
    # should be the default value, use the 'None' string instead.
23
    parameter_list = dict()
24
25
26
27
28
29
30
31
32
33
34

    def __init__(self, **kwargs):
        """Instantiate a analysis module

        The module parameters values can be passed as keyword parameters.
        """
        # parameters is a dictionary containing the actual values for each
        # module parameter.
        self.parameters = kwargs

    def _process(self, data_file, column_list, creation_modules,
35
                 creation_modules_params, parameters):
36
37
38
39
40
41
42
        """Do the actual analysis

        This method is responsible for the fitting / analysis process
        and must be implemented by each real module.

        Parameters
        ----------
43
        data_file: string
44
            Name of the file containing the observations to be fitted.
45
        column_list: array of strings
46
            Names of the columns from the data file to use in the analysis.
47
        creation_modules: array of strings
48
            Names (in the right order) of the modules to use to build the SED.
49
        creation_modules_params: array of array of dictionaries
50
51
            Array containing all the possible combinations of configurations
            for the creation_modules. Each 'inner' array has the same length as
Médéric Boquien's avatar
Médéric Boquien committed
52
53
            the creation_modules array and contains the configuration
            dictionary for the corresponding module.
54
        parameters: dictionary
55
56
57
58
59
60
61
62
63
64
            Configuration for the module.

        Returns
        -------
        The process results are saved to disk by the analysis module.

        """
        raise NotImplementedError()

    def process(self, data_file, column_list, creation_modules,
65
                creation_modules_params, parameters):
66
67
68
69
70
71
72
73
74
        """Process with the analysis

        This method is responsible for checking the module parameters before
        doing the actual processing (_process method). If a parameter is not
        given but exists in the parameter_list with a default value, this
        value is used.

        Parameters
        ----------
75
        data_file: string
76
            Name of the file containing the observations to be fitted.
77
        column_list: array of strings
78
            Names of the columns from the data file to use in the analysis.
79
        creation_modules: array of strings
80
            Names (in the right order) of the modules to use to build the SED.
81
        creation_modules_params: array of array of dictionaries
82
83
84
85
            Array containing all the possible combinations of configurations
            for the creation_modules. Each 'inner' array has the same length as
            the creation_modules array and contains the configuration
            dictionary for the corresponding module.
86
        parameters: dictionary
87
88
89
90
91
92
93
94
            Configuration for the module.

        Returns
        -------
        The process results are saved to disk by the analysis module

        Raises
        ------
95
        KeyError: when not all the needed parameters are given.
96
97
98
99
100
101

        """
        # For parameters that are present on the parameter_list with a default
        # value and that are not in the parameters dictionary, we add them
        # with their default value.
        for key in self.parameter_list:
Médéric Boquien's avatar
Médéric Boquien committed
102
            if (key not in parameters) and (
103
104
105
106
107
108
109
                    self.parameter_list[key][2] is not None):
                parameters[key] = self.parameter_list[key][2]

        # If the keys of the parameters dictionary are different from the one
        # of the parameter_list dictionary, we raises a KeyError. That means
        # that a parameter is missing (and has no default value) or that an
        # unexpected one was given.
Yannick Roehlly's avatar
Yannick Roehlly committed
110
111
        if not set(parameters) == set(self.parameter_list):
            missing_parameters = (set(self.parameter_list) - set(parameters))
Médéric Boquien's avatar
Médéric Boquien committed
112
113
            unexpected_parameters = (set(parameters) -
                                     set(self.parameter_list))
114
115
116
117
118
119
120
121
122
123
124
125
            message = ""
            if missing_parameters:
                message += ("Missing parameters: " +
                            ", ".join(missing_parameters) +
                            ".")
            if unexpected_parameters:
                message += ("Unexpected parameters: " +
                            ", ".join(unexpected_parameters) +
                            ".")
            raise KeyError("The parameters passed are different from the "
                           "expected one." + message)

Médéric Boquien's avatar
Médéric Boquien committed
126
        # We do the actual processing
127
        self._process(data_file, column_list, creation_modules,
128
                      creation_modules_params, parameters)
129
130
131
132
133
134
135


def get_module(module_name):
    """Return the main class of the module provided

    Parameters
    ----------
136
    module_name: string
137
138
139
140
        The name of the module we want to get the class.

    Returns
    -------
141
    module_class: class
142
143
144
145
146
147
148
149
    """

    try:
        module = import_module('.' + module_name, 'pcigale.analysis_modules')
        return module.Module()
    except ImportError:
        print('Module ' + module_name + ' does not exists!')
        raise
150
151


152
def adjust_errors(flux, error, tolerance, lim_flag, default_error=0.1,
153
154
155
156
157
158
159
160
161
                  systematic_deviation=0.1):
    """Adjust the errors replacing the 0 values by the default error and
    adding the systematic deviation.

    The systematic deviation change the error to:
    sqrt( error² + (flux * deviation)² )

    Parameters
    ----------
162
    flux: array of floats
163
        Fluxes.
164
    error: array of floats
165
        Observational error in the same unit as the fluxes.
166
    tolerance: float
167
        Tolerance threshold under flux error is considered as 0.
168
169
    lim_flag: boolean
        Do we process upper limits (True) or treat them as no-data (False)?
170
    default_error: float
171
172
        Default error factor used when the provided error in under the
        tolerance threshold.
173
    systematic_deviation: float
174
175
176
177
        Systematic deviation added to the error.

    Returns
    -------
178
    error: array of floats
179
180
181
182
183
184
185
186
187
188
189
        The corrected errors.

    """
    # The arrays must have the same lengths.
    if len(flux) != len(error):
        raise ValueError("The flux and error arrays must have the same "
                         "length.")

    # We copy the error array not to modify the original one.
    error = np.copy(error)

190
191
192
193
194
195
196
    # We define:
    # 1) upper limit == (lim_flag==True) and
    #                [(flux > tolerance) and (-9990. < error < tolerance)]
    # 2) no data == (flux < -9990.) and (error < -9990.)
    # Note that the upper-limit test must be performed before the no-data test
    # because if lim_flag is False, we process upper limits as no-data.
    #
197
    # Replace errors below tolerance by the default one.
Médéric Boquien's avatar
Médéric Boquien committed
198
    mask_noerror = np.logical_and(flux > tolerance, error < -9990.)
199
    error[mask_noerror] = (default_error * flux[mask_noerror])
200

201
    mask_limflag = np.logical_and.reduce(
Médéric Boquien's avatar
Médéric Boquien committed
202
                       (flux > tolerance, error >= -9990., error < tolerance))
203
204
205
206
207

    # Replace upper limits by no data if lim_flag==False
    if not lim_flag:
        flux[mask_limflag] = -9999.
        error[mask_limflag] = -9999.
208

Médéric Boquien's avatar
Médéric Boquien committed
209
    mask_ok = np.logical_and(flux > tolerance, error > tolerance)
210
211

    # Add the systematic error.
Médéric Boquien's avatar
Médéric Boquien committed
212
213
    error[mask_ok] = np.sqrt(error[mask_ok]**2 +
                             (flux[mask_ok]*systematic_deviation)**2)
214
215
216
    return error


Médéric Boquien's avatar
Médéric Boquien committed
217
218
def complete_obs_table(obs_table, used_columns, filter_list, tolerance,
                       lim_flag, default_error=0.1, systematic_deviation=0.1):
219
220
221
222
223
224
225
226
227
228
    """Complete the observation table

    For each filter:
    * If the corresponding error is not present in the used column list or in
      the table columns, add (or replace) an error column with the default
      error.
    * Adjust the error value.

    Parameters
    ----------
229
    obs_table: astropy.table.Table
230
        The observation table.
231
    used_columns: list of strings
232
        The list of columns to use in the observation table.
233
    filter_list: list of strings
234
        The list of filters used in the analysis.
235
    tolerance: float
236
        Tolerance threshold under flux error is considered as 0.
237
238
    lim_flag: boolean
        Do we process upper limits (True) or treat them as no-data (False)?
239
    default_error: float
240
241
        Default error factor used when the provided error in under the
        tolerance threshold.
242
    systematic_deviation: float
243
244
245
246
247
248
249
250
251
        Systematic deviation added to the error.

    Returns
    -------
    obs_table = astropy.table.Table
        The completed observation table

    Raises
    ------
252
    Exception: When a filter is not present in the observation table.
253
254
255
256
257
258

    """
    # TODO Print or log a warning when an error column is in the used column
    # list but is not present in the observation table.
    for name in filter_list:
        if name not in obs_table.columns:
259
            raise Exception("The filter <{}> (at least) is not present in "
260
261
262
263
264
265
266
                                "the observation table.".format(name))

        name_err = name + "_err"
        if name_err not in used_columns:
            if name_err not in obs_table.columns:
                obs_table.add_column(Column(
                    name=name_err,
267
                    data=np.ones(len(obs_table), dtype=float))*-9999.,
268
                    index=obs_table.colnames.index(name)+1
269
                )
270
271
272
273
274
275
            else:
                obs_table[name_err] = np.zeros(len(obs_table))

        obs_table[name_err] = adjust_errors(obs_table[name],
                                            obs_table[name_err],
                                            tolerance,
276
                                            lim_flag,
277
278
279
                                            default_error,
                                            systematic_deviation)
    return obs_table