Diags

@module This module contains the admitted diagnostics. Diagnostics are stored in a dictionary. Keys are literal descriptions of the diagnostic (e.g. '[NI] 5198/5200') and items are tuples including a label ('N1'), an expression for the line ratio ('L(5198)/L(5200)', and an expression for the uncertainty on the line ratio as a function of the uncertainty on the individual lines ('RMS([E(5200), E(5198)])')

Diagnostics

Bases: object

Diagnostics is the class used to manage the diagnostics and to computed physical conditions (electron temperatures and densities) from them. It is also the class that plots the diagnostic Te-Ne diagrams.

Source code in pyneb/core/diags.py

class Diagnostics(object):
    """
    Diagnostics is the class used to manage the diagnostics and to computed physical conditions 
        (electron temperatures and densities) from them. 
    It is also the class that plots the diagnostic Te-Ne diagrams.

    """    
    def __init__(self, addAll:bool=False, OmegaInterp:str='Cheb', NLevels=None):
        """
        Diagnostics constructor

        Parameters:
            addAll: Switch to include all defined diagnostics
            OmegaInterp: Parameter sent to Atom 

                **Options:**

                * 'Cheb'
                * 'linear'
            NLevels: 

        """
        self.log_ = pn.log_ 
        self.calling = 'Diagnostics'
        ##            
        # @var diags
        # The dictionary containing the diagnostics
        self.diags = {}
        ##
        # @var atomDict
        # The dictionary containing the atoms used for the diagnostics        
        self.atomDict = {}
        self.OmegaInterp = OmegaInterp
        self.NLevels = NLevels
        if addAll:
            self.addAll()
        self.ANN_n_tem=30
        self.ANN_n_den=30
        self.ANN_inst_kwargs = {'RM_type' : 'SK_ANN', 
                                'verbose' : False, 
                                'scaling' : True,
                                'use_log' : True,
                                'random_seed' : None
                                }
        self.ANN_init_kwargs = {'solver' : 'lbfgs', 
                                'activation' : 'tanh', 
                                'hidden_layer_sizes' : (10, 30, 10), 
                                'max_iter' : 20000
                                }


    def getDiagFromLabel(self, label):
        """
        Parameters:
            label(str): a diagnostic label 

        Returns:
            Return the definition of a diagnostic (the 3 or 4 elements tuple)

        **Usage:**

            diags.getDiagFromLabel('[NII] 5755/6548') 
        """
        if label in self.diags:
            return self.diags[label]
        else:
            self.log_.warn('{0} not in diagnostic'.format(label), calling=self.calling)
            return None


    def getDiags(self):
        """Return the definitions (tuples) of all the diagnostics defined in the Object.
        """
        return [self.getDiagFromLabel(label) for label in self.diags]


    def getDiagLabels(self):
        """Return the labels of all the diagnostics defined in the Object"""
        return self.diags.keys()    


    def getAllDiags(self):
        """Return the definitions (tuples) of all the possible diagnostics.

        Returns:
            (tuples): All possible diagnostics.
        """
        return diags_dict


    def getAllDiagLabels(self):
        """Return the labels of all the possible diagnostics.
        """
        return diags_dict.keys()    


    def getUniqueAtoms(self):
        """Return a numpy.ndarray of the ions needed by the diagnostics. Unique is applied to the list before returning.

        Returns:
            (np.ndarray): Ions needed by the diagnostics
        """
        return np.unique([self.diags[d][0] for d in self.diags if self.diags[d] is not None])


    def addDiag(self, label=None, diag_tuple=None, atom=None, wave_range=None):
        """Add diagnostics to the list of available diagnostics.

        It can either be one of the built-in diagnostics,
        a new, user-defined one, or a subset of the built-in diagnostics corresponding to a given atom or wavelength.

        Parameters:
            label  (str) or (list): A string or a list of strings describing the diagnostic

                **Example:**

                   - '[OIII] 4363/5007'

                If it is not a key of diags_dict (a diagnostic define by PyNeb), diag_tuple must also be specified

            diag_tuple   a 3 elements tuple containing:
                           + the atom, e.g. 'Ar5'
                           + the algebraic description of the diagnostic, in terms of line wavelengths or blends or levels, 
                             e.g. '(L(6435)+L(7006))/L(4626)'
                           + the algebraic description of the error, e.g. 
                             'RMS([E(6435)*L(6435)/(L(6435)+L(7006)), E(7006)*L(7006)/(L(6435)+L(7006)), E(4626)])'
            atom (str): The selected atom.

                **Example:**

                - 'O3'

            wave_range: the selected wavelength range


        **Usage:**
            ```python
            diags.addDiag('[OIII] 4363/5007')
            ```

            ```python
            diags.addDiag('[OIII] 5007/51m', ('O3', 'L(5007)/L(51800)', 'RMS([E(51800), E(5007)])'))
            ```

            ```python
            diags.addDiag(atom='O3')
            ```

            ```python
            diags.addDiag(wave_range=[4000, 6000])
            ```

        """
        if type(label) is list:
            for lab in label:
                self.addDiag(lab)
        elif label in self.diags:
            self.log_.warn('{0} already in diagnostic'.format(label), calling=self.calling)
        elif label in diags_dict:
            self.log_.message('Adding diag {0}'.format(label), calling=self.calling)
            self.diags[label] = diags_dict[label]
            atom = diags_dict[label][0]
        elif type(diag_tuple) is tuple:
            if len(diag_tuple) == 3:    
                self.diags[label] = diag_tuple
                atom = diag_tuple[0]
            else:
                self.log_.error('{0} is not in the list of diagnostics. The parameter diag_tuple must be a 3-elements tuple describing the diagnostic'.format(label), calling=self.calling + '.addDiag')
        elif atom is not None:
            for label in diags_dict:
                if diags_dict[label][0] == atom:
                    self.addDiag(label)
#                    if atom not in self.atomDict:
#                        self.atomDict[atom] = pn.Atom(parseAtom(atom)[0], parseAtom(atom)[1])
        elif wave_range is not None:
# To be done: add the possibility to use several independent ranges. The following bit does the trick, but
# only if all the wavelengths of a diagnostic are included in the same range
#            if type(wave_range[0]) is list:
#                for subrange in wave_range:
#                    self.addDiag(label, diag_tuple, atom, subrange)

            for label in diags_dict:
                expr = diags_dict[label][1]
                waves = []
                wave = ''
                in_wave = False
                for i in expr:
                    if i.isdigit():
                        wave = wave + i
                        in_wave = True
                    elif (i == 'm'):
                        if in_wave == True:
                            waves.append(int(wave) * 10000.)
                            in_wave = False
                            wave = ''
                    else:
                        if in_wave == True:
                            waves.append(int(wave))
                            in_wave = False
                            wave = ''
                if (min(waves) > min(wave_range)) and (max(waves) < max(wave_range)):
                    self.addDiag(label)
        else:
            self.log_.error('Bad syntax. You must either give the label of an existing diagnostic, or the label and the tuple of a new one,' + 
                            'or an ion, or a wave range. label={0}, diag_tuple={1}, atom={2}, wave_range={3}'.format(label, diag_tuple, atom, wave_range),
                            calling=self.calling + '.addDiag')
        if atom not in self.atomDict and (type(label) is not list):
            this_atom, this_spec, this_rec = parseAtom2(atom)
            if this_rec == '':
                self.atomDict[atom] = pn.Atom(this_atom, this_spec, NLevels=self.NLevels)
            elif this_rec == 'r':
                self.atomDict[atom] = pn.RecAtom(this_atom, this_spec)


    def addAll(self):
        """Insert all the possible diagnostics in the Object"""
        for label in diags_dict:
            self.addDiag(label)


    def delDiag(self, label=None):
        """
        Remove a diagnostic, based on its label.

        Parameters:
            label (str): The label of the diagnostic ('all': removes all the diagnostics)

        """
        if label == 'all':
            for item in self.diags:
                self.delDiag(item)
        elif label in self.diags:
            del self.diags[label]
        else:
            self.log_.warn('{0} not in diagnostic, cannot be removed'.format(label), calling=self.calling)


    def addDiagsFromObs(self, obs):
        """Add all the possible diagnostics that can be computed from an Observation object

        Parameters:
            obs : An Observation object

        **Usage:**

            diags.addDiagsFromObs(obs)
        """

        if not isinstance(obs, pn.Observation):
            pn.log_.error('The argument must be an Observation object', calling=self.calling + 'addDiagsFromObs')
        old_level = pn.log_.level
        def I(i, j):
            wave = atom.wave_Ang[i - 1, j - 1]
            corrIntens = obs.getLine(sym, spec, wave).corrIntens
            return corrIntens
        def L(wave):
            corrIntens = obs.getLine(sym, spec, wave).corrIntens
            return corrIntens
        def B(label):
            full_label = atom + '_' + label
            corrIntens = obs.getLine(label=full_label).corrIntens
            return corrIntens
        def S(label):
            full_label = atom + '_' + label + 'A'
            corrIntens = obs.getLine(label=full_label).corrIntens
            return corrIntens

        for label in diags_dict:
            atom, diag_expression, error = diags_dict[label]
            sym, spec, rec = parseAtom2(atom)
            if label == '[OII] 3727+/7325+c':
                try:
                    diag_value = eval(diag_expression)
                except Exception as ex:
                    pn.log_.level = old_level
                    pn.log_.debug('Diag not valid {} {}'.format(label, diag_expression))
            try:
                pn.log_.level = 1
                diag_value = eval(diag_expression)
                pn.log_.level = old_level
                if atom not in self.atomDict:
                    if rec == 'r':
                        self.atomDict[atom] = pn.RecAtom(atom=sym+spec)
                    else:
                        self.atomDict[atom] = pn.Atom(atom=atom, NLevels=self.NLevels)
                self.addDiag(label)
            except Exception as ex:
                pn.log_.level = old_level
                pn.log_.debug('Diag not valid {} {}'.format(label, diag_expression))


    def setAtoms(self, atom_dic):
        """Define the dictionary containing the atoms used for the diagnostics.

        A dictionary of atom instantiations refereed by atom strings, for 

        Parameters:

            atom_dic : a dictionary of Atom instances, indexed by atom strings.
                    **Example:**

                        {'O3': pn.Atom('O', 3)}

        """
        if type(atom_dic) != type({}):
            self.log_.error('the parameter must be a dictionary.', calling=self.calling + '.setAtoms')
            return None
        for atom in atom_dic:
            if not isinstance(atom_dic[atom], pn.Atom) and not isinstance(atom_dic[atom], pn.RecAtom):
                self.log_.error('the parameter must be a dictionary of Atom.', calling=self.calling + '.setAtoms')
                return None
        self.atomDict = atom_dic


    def addClabel(self, label, clabel):
        """Add an alternative label to a diagnostic that can be used when plotting diagnostic diagrams.

        Parameters:
            label (str):
            clabel (str):

        """
        if label in self.diags:
            self.diags[label] = (self.diags[label][0], self.diags[label][1], self.diags[label][2], clabel)
        else:
            pn.log_.warn('Try to add clabel in undefined label {0}'.format(label), calling=self.calling)


    def plot(self, emis_grids, obs, quad=True, i_obs=None, alpha=0.3, ax=None, error_band=True,
    return_CS=False, 
             col_dic={'C':'cyan', 'N':'blue', 'O':'green', 'Ne':'magenta',
                      'Ar':'red', 'Cl':'magenta', 'S':'black', 'Fe':'blue'}):
        """Plotting tool to generate Te-Ne diagrams.

        Parameters:
            emis_grids:    A dictionary of EmisGrid objects refereed by their atom strings (e.g. 'O3')
                            This can for example be the output of pn.getEmisGridDict()
            obs:           A pn.Observation object that is supposed to contain the line intensities
                            used for the plot (corrected intensities).
            quad:          If True (default) the sum of the error is quadratic,otherwise is linear.
            i_obs:         reference for the observation to be plotted, in case there is more than one
                            in the obs object
            alpha:         Transparency for the error bands in the plot
            error_band:    Boolean: plot [default] an error band

            return_CS     [False] If True, return a list of the contour plots
            col_dic       Colors for the different ions.            
        **Usage:**
            diags.plot(emisgrids, obs, i_obs=3)
        """
        if not pn.config.INSTALLED['plt']: 
            pn.log_.error('Matplotlib not available, no plot', calling=self.calling + '.plot')
            return None
        else:
            import matplotlib.pyplot as plt
        if type(emis_grids) != type({}):
            self.log_.error('the first parameter must be a dictionary', calling=self.calling + '.plot')
            return None
        for em in emis_grids:
            if not isinstance(emis_grids[em], pn.EmisGrid):
                self.log_.error('the first parameter must be a dictionary of EmisGrid', calling=self.calling + '.plot')
                return None
        if not isinstance(obs, pn.Observation):
            self.log_.error('the second parameter must be an Observation', calling=self.calling + '.plot')
            return None
        if (i_obs is None) and (obs.n_obs != 1):
            self.log_.error('i_obs must be specified when obs is multiple. try i_obs=0', calling=self.calling)
            return None
        if ax is None:
            f, ax = plt.subplots()
        else:
            f = plt.gcf()
        X = np.log10(emis_grids[list(emis_grids.keys())[0]].den2D)
        Y = emis_grids[list(emis_grids.keys())[0]].tem2D
        CSs = []
        for label in self.diags:
            self.log_.message('plotting {0}'.format(label), calling=self.calling)
            diag = self.diags[label]
            atom = diag[0]
            def I(i, j):
                return emis_grids[atom].getGrid(lev_i=i, lev_j=j)
            def L(wave):
                return emis_grids[atom].getGrid(wave=wave)
            def S(label):
                return emis_grids[atom].getGrid(label=label)
            def B(label, I=I, L=L):
                full_label = atom + '_' + label
                if full_label in BLEND_LIST:
                    to_eval = BLEND_LIST[full_label]
                else:
                    self.log_.warn('{0} not in BLEND_LIST'.format(full_label), calling=self.calling)
                    return None
                return eval(to_eval)
            diag_map = eval(diag[1])
            try:
                diag_map = eval(diag[1])
            except:
                diag_map = None
                self.log_.warn('diag {0} {1} not used'.format(diag[0], diag[1]), calling=self.calling)
            if diag_map is not None:
                sym, spec, rec = parseAtom2(atom)
                def I(i, j):
                    wave = emis_grids[atom].atom.wave_Ang[i - 1, j - 1]
                    corrIntens = obs.getLine(sym, spec, wave).corrIntens
                    if i_obs is None:
                        return corrIntens
                    else:
                        return corrIntens[i_obs]
                def L(wave):
                    corrIntens = obs.getLine(sym, spec, wave).corrIntens
                    if i_obs is None:
                        return corrIntens
                    else:
                        return corrIntens[i_obs]
                def S(label):
                    full_label = atom + '_' + label + 'A'
                    corrIntens = obs.getLine(label=full_label).corrIntens
                    if i_obs is None:
                        return corrIntens
                    else:
                        return corrIntens[i_obs]                    
                def B(label):
                    full_label = atom + '_' + label
                    corrIntens = obs.getLine(label=full_label).corrIntens
                    if i_obs is None:
                        return corrIntens
                    else:
                        return corrIntens[i_obs]
                try:
                    ee = eval(diag[1])
                    if type(ee) == np.float64:
                        diag_value = ee
                    else:
                        diag_value = ee[0]
                except:
                    #print(ee, type(ee))
                    diag_value = None
                    self.log_.warn('A line in diagnostic {0} of {1}{2} is missing'.format(diag[1], sym, spec),
                                   calling=self.calling)
                if diag_value is not None and not np.isnan(diag_value) and not np.isinf(diag_value):
                    def E(wave):
                        err = obs.getLine(sym, spec, wave).corrError
                        if i_obs is None:
                            return err
                        else:
                            return err[i_obs]
                    def BE(label):
                        full_label = atom + '_' + label
                        err = obs.getLine(label=full_label).corrError
                        if i_obs is None:
                            return err
                        else:
                            return err[i_obs]
                    def SE(label):
                        full_label = atom + '_' + label + 'A'
                        err = obs.getLine(label=full_label).corrError
                        if i_obs is None:
                            return err
                        else:
                            return err[i_obs]                        
                    if quad is True:
                        RMS = lambda err: np.sqrt((np.asarray(err) ** 2.).sum())
                    else:
                        RMS = lambda err: (np.asarray(err)).sum() 
                    tol_value = eval(diag[2])
                    if sym in col_dic:
                        col = col_dic[sym]
                    else:
                        col = 'black'
                    style_dic = {'1':'-', '2':'--', '3':':', '4':'-.', '5':'-', '6':'--'}
                    style = style_dic[spec]
                    if tol_value > 0. and error_band:
                        levels = [(1 - tol_value) * diag_value, (1 + tol_value) * diag_value]
                        if levels[0] < levels[1]:
                            #pn.log_.debug('{} levels {}'.format(label, levels), calling=self.calling)
                            CS = ax.contourf(X, Y, diag_map, levels=levels, alpha=alpha, colors=col)
                            CSs.append(CS)
                    CS = ax.contour(X, Y, diag_map, levels=[diag_value], colors=col, linestyles=style)
                    CSs.append(CS)
                    try:
                        ax.set_xlabel(r'log(n$_{\rm e}$) [cm$^{-3}$]')
                        ax.set_ylabel(r'T$_{\rm e}$ [K]')
                        if len(diag) >= 4:
                            fmt = diag[3]
                        else:
                            fmt = '[{0}{1}]'.format(sym, int_to_roman(int(spec)))
                        ax.clabel(CS, inline=True, fmt=fmt, fontsize=15, colors=col)
                        if type(diag_value) is np.ndarray:
                            diag_value = diag_value[0]
                        self.log_.message('plotted {0}: {1} = {2:.2} with error of {3:.2} %'.format(fmt, label, diag_value, tol_value * 100),
                                          calling=self.calling)
                    except:
                        self.log_.message('NOT plotted {0} {1}'.format(fmt, label),
                                          calling=self.calling)

        if return_CS:
            return CSs

    def getCrossTemDen(self, diag_tem, diag_den, value_tem=None, value_den=None, obs=None, i_obs=None,
                       guess_tem=10000, tol_tem=1., tol_den=1., max_iter=5, maxError=1e-3,
                       start_tem= -1, end_tem= -1, start_den= -1, end_den= -1, use_ANN=False, 
                       limit_res=False, ANN=None):
        """
        Cross-converge the temperature and density derived from two sensitive line ratios, by inputting the quantity 
        derived with one line ratio into the other and then iterating.
        The temperature- and density-sensitive ratios can be input directly or as an Observation object

        Parameters:

            diag_tem:   temperature-sensitive diagnostic line ratio
            diag_den:   density-sensitive diagnostic line ratio
            value_tem:  value of the temperature-sensitive diagnostic
            value_den:  value of the density-sensitive diagnostic
            obs:        np.Observation object. Values for observed temperature and density diagnostics are
                            taken from it if value_tem and value_den are None
            i_obs:      index of the observations to be used from obs. 
                            If None, all the observations are considered.
                            May be a boolean array
            guess_tem:  temperature assumed in the first iteration, in K
            tol_tem:    tolerance of the temperature result, in %
            tol_den:    tolerance of the density result, in %
            max_iter:   maximum number of iterations to be performed, after which the function will throw a result
            maxError:   maximum error in the calls to getTemDen, in %
            start_tem:, end_tem  lower and upper limit of the explored temperature range 
            start_den:, end_den  lower and upper limit of the explored density range 
            use_ANN:    if True, an Analogic Neural Network will be used for the determination of Te and Ne.
                            manage_RM from mwinai library is used.
                            the hyper_parameters can be set-up with the self.ANN_inst_kwargs and
                            self.ANN_init_kwargs dictionnaries.
                            self.ANN_n_tem=30 and self.ANN_n_den=30 are the number of Te and Ne
                            used to train. May also be changed before calling getCrossTemDen
            limit_res:  in case of using ANN, if limit_res, the tem and den values out of the start_tem, end_tem,
                        start_den, end_den are set to np.nan. Otherwise, extrapolation is allowed.
            ANN:        if string, filename where to read AI4neb ANN. Otherwise, ANN is a manage_RM object.
                        In both casesm the ANN needs to already be trained

        **Example:**

            tem, den = diags.getCrossTemDen('[OIII] 4363/5007', '[SII] 6731/6716', 0.0050, 1.0, 
                        guess_tem=10000, tol_tem = 1., tol_den = 1., max_iter = 5)

        """
        # TODO:
        # Define a lower/upper density and temperature 
        # to be used in case the diag is pointing to values out of the boundaires
        if diag_tem not in self.diags:
            self.addDiag(diag_tem)
        atom_tem = self.diags[diag_tem][0]
        elem_tem, spec_tem, rec = parseAtom2(atom_tem)
        if atom_tem not in self.atomDict:
            self.atomDict[atom_tem] = pn.Atom(elem_tem, spec_tem, self.OmegaInterp, NLevels=self.NLevels)
        atom_tem = self.atomDict[atom_tem]
        if diag_den not in self.diags:
            self.addDiag(diag_den)
        atom_den = self.diags[diag_den][0]
        elem_den, spec_den, rec = parseAtom2(self.diags[diag_den][0])
        if (atom_den) not in self.atomDict:
            self.atomDict[atom_den] = pn.Atom(elem_den, spec_den, self.OmegaInterp, NLevels=self.NLevels)
        atom_den = self.atomDict[atom_den]
        eval_tem = self.diags[diag_tem][1]
        eval_den = self.diags[diag_den][1]
        calling = 'Diag.getCrossTemDen %s %s' % (diag_tem, diag_den)
        if value_tem is None:
            def L(wave):
                if i_obs is None:
                    return obs.getLine(elem_tem, spec_tem, wave).corrIntens
                else:
                    return obs.getLine(elem_tem, spec_tem, wave).corrIntens[i_obs]
            def I(i, j):
                wave = atom_tem.wave_Ang[i - 1, j - 1]
                if i_obs is None:
                    return obs.getLine(elem_tem, spec_tem, wave).corrIntens
                else:
                    return obs.getLine(elem_tem, spec_tem, wave).corrIntens[i_obs]
            def B(label, I=I, L=L):
                full_label = elem_tem + spec_tem + '_' + label
                if i_obs is None:
                    corrIntens = obs.getLine(label=full_label).corrIntens
                else:
                    corrIntens = obs.getLine(label=full_label).corrIntens[i_obs]
                return corrIntens
            #pn.log_.debug('to eval is {0}'.format(eval_tem), calling=calling + 'TEST')
            try:
                value_tem = eval(eval_tem)
                #pn.log_.debug('to eval = {0}'.format(value_tem), calling=calling + 'TEST')
            except:
                pn.log_.warn('No value for {0} {1}: {2} in obs'.format(elem_tem, spec_tem, diag_tem), calling=calling)
                return None
        else:
            if type(value_tem) == type([]): value_tem = np.asarray(value_tem)
        if value_den is None:
            def L(wave): 
                if i_obs is None:
                    return obs.getLine(elem_den, spec_den, wave).corrIntens
                else:
                    return obs.getLine(elem_den, spec_den, wave).corrIntens[i_obs]
            def I(i, j):
                wave = atom_den.wave_Ang[i - 1, j - 1]
                pn.log_.debug('wave is {0}'.format(wave), calling=calling + 'TEST3')
                if i_obs is None:
                    return obs.getLine(elem_den, spec_den, wave).corrIntens
                else:
                    return obs.getLine(elem_den, spec_den, wave).corrIntens[i_obs]
            def B(label, I=I, L=L):
                full_label = elem_den + spec_den + '_' + label
                if i_obs is None:
                    corrIntens = obs.getLine(label=full_label).corrIntens
                else:
                    corrIntens = obs.getLine(label=full_label).corrIntens[i_obs]
                return corrIntens
            #pn.log_.debug('to eval is {0}'.format(eval_den), calling=calling + ' TEST')
            try:
                value_den = eval(eval_den)
                #pn.log_.debug('to eval = {0}'.format(value_den), calling=calling + ' TEST1')
            except:
                pn.log_.warn('No value for {0} {1}: {2} in obs'.format(elem_den, spec_den, diag_den), calling=calling)
                return None
        else:
            if type(value_den) == type([]): value_den = np.asarray(value_den)
        if use_ANN:

            if config.INSTALLED['ai4neb']:
                from ai4neb import manage_RM

            if not config.INSTALLED['ai4neb']:
                self.log_.error('_getPopulations_ANN cannot be used in absence of ai4neb package',
                              calling=self.calling)
                return None
            if start_tem == -1:
                tem_min = 3000.
            else:
                tem_min = start_tem
            if end_tem == -1:
                tem_max = 20000.
            else:
                tem_max = end_tem
            if start_den == -1:
                den_min = 10.
            else:
                den_min = start_den
            if end_den == -1:
                den_max = 1e6
            else:
                den_max = end_den
            if ANN is None:
                # define emisGrid objects to generate Te-Ne emissionmaps
                tem_EG = pn.EmisGrid(atomObj=atom_tem, 
                                     n_tem=self.ANN_n_tem, n_den=self.ANN_n_den, 
                                     tem_min=tem_min, tem_max=tem_max,
                                     den_min=den_min, den_max=den_max)
                den_EG = pn.EmisGrid(atomObj=atom_den, 
                                     n_tem=self.ANN_n_tem, n_den=self.ANN_n_den, 
                                     tem_min=tem_min, tem_max=tem_max,
                                     den_min=den_min, den_max=den_max)
                # compute emission line ratio maps in the Te-Ne space
                tem_2D = tem_EG.getGrid(to_eval = eval_tem)
                den_2D = den_EG.getGrid(to_eval = eval_den)
                # X is a set of line ratio pairs
                X = np.array((tem_2D.ravel(), den_2D.ravel())).T
                # y is a set of corresponding Te-Ne pairs
                y = np.array((tem_EG.tem2D.ravel()/1e4, np.log10(den_EG.den2D.ravel()))).T
                # Instantiate, init and train the ANN
                self.ANN = manage_RM(X_train=X, y_train=y, **self.ANN_inst_kwargs)
                self.ANN.init_RM(**self.ANN_init_kwargs)
                self.ANN.train_RM()
            else:
                if type(ANN) is str:
                    self.ANN = manage_RM(RM_filename=ANN)
                else:
                    self.ANN = ANN
            # set the test values to the one we are looking for
            shape = value_tem.shape
            self.ANN.set_test(np.array((value_tem.ravel(), value_den.ravel())).T)
            # predict the result and denormalize them
            self.ANN.predict()
            if self.ANN.isfin is None:
                tem = self.ANN.pred[:,0]*1e4
                den = 10**self.ANN.pred[:,1]
            else:
                tem = np.zeros_like(value_tem.ravel()) * -10
                tem[self.ANN.isfin] = self.ANN.pred[:,0]*1e4
                den = np.zeros_like(value_tem.ravel()) * -10
                den[self.ANN.isfin] = 10**self.ANN.pred[:,1]
            tem = np.reshape(tem, shape)
            den = np.reshape(den, shape)
            if limit_res:
                mask = (tem<tem_min) | (tem>tem_max)
                tem[mask] = np.nan
                pn.log_.debug('Removing {} points out of Te range'.format(mask.sum()), calling='getCrossTemDen')
                mask = (den<den_min) | (den>den_max)
                den[mask] = np.nan
                pn.log_.debug('Removing {} points out of Ne range'.format(mask.sum()), calling='getCrossTemDen')
        else:
            den = atom_den.getTemDen(value_den, tem=guess_tem, to_eval=eval_den,
                                     maxError=maxError, start_x=start_den, end_x=end_den)

            tem = atom_tem.getTemDen(value_tem, den=den, to_eval=eval_tem,
                                     maxError=maxError, start_x=start_tem, end_x=end_tem)
    #        self.log_.debug('tem: ' + str(tem) + ' den:' + str(den), calling='getCrossTemDen')
            no_conv = np.ones_like(den).astype(bool)
            n_tot = np.asarray(value_tem).size
            for i in np.arange(max_iter):
                if type(tem) == type(1.):
                    tem_old = tem
                else:
                    tem_old = tem.copy()
                if type(den) == type(1.):
                    den_old = den
                else:
                    den_old = den.copy()

                if n_tot > 1:
                    den[no_conv] = atom_den.getTemDen(value_den[no_conv], tem=tem_old[no_conv],
                                                      to_eval=eval_den, start_x=start_den, end_x=end_den)
                    tem[no_conv] = atom_tem.getTemDen(value_tem[no_conv], den=den_old[no_conv],
                                                      to_eval=eval_tem, start_x=start_tem, end_x=end_tem)
                else:
                    den = atom_den.getTemDen(value_den, tem=tem_old, to_eval=eval_den, start_x=start_den, end_x=end_den)
                    tem = atom_tem.getTemDen(value_tem, den=den_old, to_eval=eval_tem, start_x=start_tem, end_x=end_tem)

                no_conv = ((abs(den_old - den) / den * 100) > tol_den) | ((abs(tem_old - tem) / tem * 100) > tol_tem)
                if type(no_conv) == type(True):
                    n_no_conv = int(no_conv)
                else:
                    n_no_conv = no_conv.sum()

                pn.log_.message('{0} (max={1}): not converged {2} of {3}.'.format(i, max_iter, n_no_conv, n_tot),
                                calling=calling)
                if n_no_conv == 0:
                    return tem, den
            if n_tot == 1:
                tem = np.nan
                den = np.nan
            else:
                tem[no_conv] = np.nan
                den[no_conv] = np.nan
        return tem, den

    def getDiagLimits(self, diag):
        """Return the low and high density values for a given diagnostic
        """
        atom = self.atomDict[self.diags[diag][0]]
        to_eval = self.diags[diag][1]
        HDR = atom.getHighDensRatio(to_eval = to_eval)
        LDR = atom.getLowDensRatio(to_eval = to_eval)
        return(np.sort((LDR, HDR)))

    def eval_diag(self, label, obs):
        """
        Parameters
        ----------
            label : diagnostic label(e.g. '[OIII] 4363/5007')
                A string of a key included in the self.diags dictionnary.
            obs : an Observation object 
        Returns
        -------
            (np.array): The evaluation of the diagnostic corresponding to the label.

        """
        if label not in self.diags:
            self.log_.error('Unknown diagnostic: {}'.format(label), calling='eval_diag')
        atom, diag_expression, error = self.diags[label]
        sym, spec, rec = parseAtom2(atom)
        def I(i, j):
            wave = atom.wave_Ang[i - 1, j - 1]
            corrIntens = obs.getLine(sym, spec, wave).corrIntens
            return corrIntens
        def L(wave):
            corrIntens = obs.getLine(sym, spec, wave).corrIntens
            return corrIntens
        def B(label):
            full_label = atom + '_' + label
            corrIntens = obs.getLine(label=full_label).corrIntens
            return corrIntens
        def S(label):
            full_label = atom + '_' + label + 'A'
            corrIntens = obs.getLine(label=full_label).corrIntens
            return corrIntens
        diag_value = eval(diag_expression)
        return diag_value    

__init__(addAll=False, OmegaInterp='Cheb', NLevels=None)

Diagnostics constructor

Parameters:

Name	Type	Description	Default
addAll	bool	Switch to include all defined diagnostics	False
OmegaInterp	str	Parameter sent to Atom Options: 'Cheb' 'linear'	'Cheb'
NLevels			None

Source code in pyneb/core/diags.py

def __init__(self, addAll:bool=False, OmegaInterp:str='Cheb', NLevels=None):
    """
    Diagnostics constructor

    Parameters:
        addAll: Switch to include all defined diagnostics
        OmegaInterp: Parameter sent to Atom 

            **Options:**

            * 'Cheb'
            * 'linear'
        NLevels: 

    """
    self.log_ = pn.log_ 
    self.calling = 'Diagnostics'
    ##            
    # @var diags
    # The dictionary containing the diagnostics
    self.diags = {}
    ##
    # @var atomDict
    # The dictionary containing the atoms used for the diagnostics        
    self.atomDict = {}
    self.OmegaInterp = OmegaInterp
    self.NLevels = NLevels
    if addAll:
        self.addAll()
    self.ANN_n_tem=30
    self.ANN_n_den=30
    self.ANN_inst_kwargs = {'RM_type' : 'SK_ANN', 
                            'verbose' : False, 
                            'scaling' : True,
                            'use_log' : True,
                            'random_seed' : None
                            }
    self.ANN_init_kwargs = {'solver' : 'lbfgs', 
                            'activation' : 'tanh', 
                            'hidden_layer_sizes' : (10, 30, 10), 
                            'max_iter' : 20000
                            }

addAll()

Insert all the possible diagnostics in the Object

Source code in pyneb/core/diags.py

def addAll(self):
    """Insert all the possible diagnostics in the Object"""
    for label in diags_dict:
        self.addDiag(label)

addClabel(label, clabel)

Add an alternative label to a diagnostic that can be used when plotting diagnostic diagrams.

Parameters:

Name	Type	Description	Default
label	str		required
clabel	str		required

Source code in pyneb/core/diags.py

def addClabel(self, label, clabel):
    """Add an alternative label to a diagnostic that can be used when plotting diagnostic diagrams.

    Parameters:
        label (str):
        clabel (str):

    """
    if label in self.diags:
        self.diags[label] = (self.diags[label][0], self.diags[label][1], self.diags[label][2], clabel)
    else:
        pn.log_.warn('Try to add clabel in undefined label {0}'.format(label), calling=self.calling)

addDiag(label=None, diag_tuple=None, atom=None, wave_range=None)

Add diagnostics to the list of available diagnostics.

It can either be one of the built-in diagnostics, a new, user-defined one, or a subset of the built-in diagnostics corresponding to a given atom or wavelength.

Parameters:

Name	Type	Description	Default
label	(str) or (list	A string or a list of strings describing the diagnostic Example: '[OIII] 4363/5007' If it is not a key of diags_dict (a diagnostic define by PyNeb), diag_tuple must also be specified	None
diag_tuple	a 3 elements tuple containing	+ the atom, e.g. 'Ar5' + the algebraic description of the diagnostic, in terms of line wavelengths or blends or levels, e.g. '(L(6435)+L(7006))/L(4626)' + the algebraic description of the error, e.g. 'RMS([E(6435)*L(6435)/(L(6435)+L(7006)), E(7006)*L(7006)/(L(6435)+L(7006)), E(4626)])'	None
atom	str	The selected atom. Example: 'O3'	None
wave_range		the selected wavelength range	None

Usage: python diags.addDiag('[OIII] 4363/5007')

```python
diags.addDiag('[OIII] 5007/51m', ('O3', 'L(5007)/L(51800)', 'RMS([E(51800), E(5007)])'))
```

```python
diags.addDiag(atom='O3')
```

```python
diags.addDiag(wave_range=[4000, 6000])
```

Source code in pyneb/core/diags.py

    def addDiag(self, label=None, diag_tuple=None, atom=None, wave_range=None):
        """Add diagnostics to the list of available diagnostics.

        It can either be one of the built-in diagnostics,
        a new, user-defined one, or a subset of the built-in diagnostics corresponding to a given atom or wavelength.

        Parameters:
            label  (str) or (list): A string or a list of strings describing the diagnostic

                **Example:**

                   - '[OIII] 4363/5007'

                If it is not a key of diags_dict (a diagnostic define by PyNeb), diag_tuple must also be specified

            diag_tuple   a 3 elements tuple containing:
                           + the atom, e.g. 'Ar5'
                           + the algebraic description of the diagnostic, in terms of line wavelengths or blends or levels, 
                             e.g. '(L(6435)+L(7006))/L(4626)'
                           + the algebraic description of the error, e.g. 
                             'RMS([E(6435)*L(6435)/(L(6435)+L(7006)), E(7006)*L(7006)/(L(6435)+L(7006)), E(4626)])'
            atom (str): The selected atom.

                **Example:**

                - 'O3'

            wave_range: the selected wavelength range


        **Usage:**
            ```python
            diags.addDiag('[OIII] 4363/5007')
            ```

            ```python
            diags.addDiag('[OIII] 5007/51m', ('O3', 'L(5007)/L(51800)', 'RMS([E(51800), E(5007)])'))
            ```

            ```python
            diags.addDiag(atom='O3')
            ```

            ```python
            diags.addDiag(wave_range=[4000, 6000])
            ```

        """
        if type(label) is list:
            for lab in label:
                self.addDiag(lab)
        elif label in self.diags:
            self.log_.warn('{0} already in diagnostic'.format(label), calling=self.calling)
        elif label in diags_dict:
            self.log_.message('Adding diag {0}'.format(label), calling=self.calling)
            self.diags[label] = diags_dict[label]
            atom = diags_dict[label][0]
        elif type(diag_tuple) is tuple:
            if len(diag_tuple) == 3:    
                self.diags[label] = diag_tuple
                atom = diag_tuple[0]
            else:
                self.log_.error('{0} is not in the list of diagnostics. The parameter diag_tuple must be a 3-elements tuple describing the diagnostic'.format(label), calling=self.calling + '.addDiag')
        elif atom is not None:
            for label in diags_dict:
                if diags_dict[label][0] == atom:
                    self.addDiag(label)
#                    if atom not in self.atomDict:
#                        self.atomDict[atom] = pn.Atom(parseAtom(atom)[0], parseAtom(atom)[1])
        elif wave_range is not None:
# To be done: add the possibility to use several independent ranges. The following bit does the trick, but
# only if all the wavelengths of a diagnostic are included in the same range
#            if type(wave_range[0]) is list:
#                for subrange in wave_range:
#                    self.addDiag(label, diag_tuple, atom, subrange)

            for label in diags_dict:
                expr = diags_dict[label][1]
                waves = []
                wave = ''
                in_wave = False
                for i in expr:
                    if i.isdigit():
                        wave = wave + i
                        in_wave = True
                    elif (i == 'm'):
                        if in_wave == True:
                            waves.append(int(wave) * 10000.)
                            in_wave = False
                            wave = ''
                    else:
                        if in_wave == True:
                            waves.append(int(wave))
                            in_wave = False
                            wave = ''
                if (min(waves) > min(wave_range)) and (max(waves) < max(wave_range)):
                    self.addDiag(label)
        else:
            self.log_.error('Bad syntax. You must either give the label of an existing diagnostic, or the label and the tuple of a new one,' + 
                            'or an ion, or a wave range. label={0}, diag_tuple={1}, atom={2}, wave_range={3}'.format(label, diag_tuple, atom, wave_range),
                            calling=self.calling + '.addDiag')
        if atom not in self.atomDict and (type(label) is not list):
            this_atom, this_spec, this_rec = parseAtom2(atom)
            if this_rec == '':
                self.atomDict[atom] = pn.Atom(this_atom, this_spec, NLevels=self.NLevels)
            elif this_rec == 'r':
                self.atomDict[atom] = pn.RecAtom(this_atom, this_spec)

addDiagsFromObs(obs)

Add all the possible diagnostics that can be computed from an Observation object

Parameters:

Name	Type	Description	Default
obs		An Observation object	required

Usage:

diags.addDiagsFromObs(obs)

Source code in pyneb/core/diags.py

def addDiagsFromObs(self, obs):
    """Add all the possible diagnostics that can be computed from an Observation object

    Parameters:
        obs : An Observation object

    **Usage:**

        diags.addDiagsFromObs(obs)
    """

    if not isinstance(obs, pn.Observation):
        pn.log_.error('The argument must be an Observation object', calling=self.calling + 'addDiagsFromObs')
    old_level = pn.log_.level
    def I(i, j):
        wave = atom.wave_Ang[i - 1, j - 1]
        corrIntens = obs.getLine(sym, spec, wave).corrIntens
        return corrIntens
    def L(wave):
        corrIntens = obs.getLine(sym, spec, wave).corrIntens
        return corrIntens
    def B(label):
        full_label = atom + '_' + label
        corrIntens = obs.getLine(label=full_label).corrIntens
        return corrIntens
    def S(label):
        full_label = atom + '_' + label + 'A'
        corrIntens = obs.getLine(label=full_label).corrIntens
        return corrIntens

    for label in diags_dict:
        atom, diag_expression, error = diags_dict[label]
        sym, spec, rec = parseAtom2(atom)
        if label == '[OII] 3727+/7325+c':
            try:
                diag_value = eval(diag_expression)
            except Exception as ex:
                pn.log_.level = old_level
                pn.log_.debug('Diag not valid {} {}'.format(label, diag_expression))
        try:
            pn.log_.level = 1
            diag_value = eval(diag_expression)
            pn.log_.level = old_level
            if atom not in self.atomDict:
                if rec == 'r':
                    self.atomDict[atom] = pn.RecAtom(atom=sym+spec)
                else:
                    self.atomDict[atom] = pn.Atom(atom=atom, NLevels=self.NLevels)
            self.addDiag(label)
        except Exception as ex:
            pn.log_.level = old_level
            pn.log_.debug('Diag not valid {} {}'.format(label, diag_expression))

delDiag(label=None)

Remove a diagnostic, based on its label.

Parameters:

Name	Type	Description	Default
label	str	The label of the diagnostic ('all': removes all the diagnostics)	None

Source code in pyneb/core/diags.py

def delDiag(self, label=None):
    """
    Remove a diagnostic, based on its label.

    Parameters:
        label (str): The label of the diagnostic ('all': removes all the diagnostics)

    """
    if label == 'all':
        for item in self.diags:
            self.delDiag(item)
    elif label in self.diags:
        del self.diags[label]
    else:
        self.log_.warn('{0} not in diagnostic, cannot be removed'.format(label), calling=self.calling)

eval_diag(label, obs)

Parameters

label : diagnostic label(e.g. '[OIII] 4363/5007')
    A string of a key included in the self.diags dictionnary.
obs : an Observation object

Returns

(np.array): The evaluation of the diagnostic corresponding to the label.

Source code in pyneb/core/diags.py

def eval_diag(self, label, obs):
    """
    Parameters
    ----------
        label : diagnostic label(e.g. '[OIII] 4363/5007')
            A string of a key included in the self.diags dictionnary.
        obs : an Observation object 
    Returns
    -------
        (np.array): The evaluation of the diagnostic corresponding to the label.

    """
    if label not in self.diags:
        self.log_.error('Unknown diagnostic: {}'.format(label), calling='eval_diag')
    atom, diag_expression, error = self.diags[label]
    sym, spec, rec = parseAtom2(atom)
    def I(i, j):
        wave = atom.wave_Ang[i - 1, j - 1]
        corrIntens = obs.getLine(sym, spec, wave).corrIntens
        return corrIntens
    def L(wave):
        corrIntens = obs.getLine(sym, spec, wave).corrIntens
        return corrIntens
    def B(label):
        full_label = atom + '_' + label
        corrIntens = obs.getLine(label=full_label).corrIntens
        return corrIntens
    def S(label):
        full_label = atom + '_' + label + 'A'
        corrIntens = obs.getLine(label=full_label).corrIntens
        return corrIntens
    diag_value = eval(diag_expression)
    return diag_value    

getAllDiagLabels()

Return the labels of all the possible diagnostics.

Source code in pyneb/core/diags.py

def getAllDiagLabels(self):
    """Return the labels of all the possible diagnostics.
    """
    return diags_dict.keys()    

getAllDiags()

Return the definitions (tuples) of all the possible diagnostics.

Returns:

Type	Description
tuples	All possible diagnostics.

Source code in pyneb/core/diags.py

def getAllDiags(self):
    """Return the definitions (tuples) of all the possible diagnostics.

    Returns:
        (tuples): All possible diagnostics.
    """
    return diags_dict

getCrossTemDen(diag_tem, diag_den, value_tem=None, value_den=None, obs=None, i_obs=None, guess_tem=10000, tol_tem=1.0, tol_den=1.0, max_iter=5, maxError=0.001, start_tem=-1, end_tem=-1, start_den=-1, end_den=-1, use_ANN=False, limit_res=False, ANN=None)

Cross-converge the temperature and density derived from two sensitive line ratios, by inputting the quantity derived with one line ratio into the other and then iterating. The temperature- and density-sensitive ratios can be input directly or as an Observation object

Parameters:

Name	Type	Description	Default
diag_tem		temperature-sensitive diagnostic line ratio	required
diag_den		density-sensitive diagnostic line ratio	required
value_tem		value of the temperature-sensitive diagnostic	None
value_den		value of the density-sensitive diagnostic	None
obs		np.Observation object. Values for observed temperature and density diagnostics are taken from it if value_tem and value_den are None	None
i_obs		index of the observations to be used from obs. If None, all the observations are considered. May be a boolean array	None
guess_tem		temperature assumed in the first iteration, in K	10000
tol_tem		tolerance of the temperature result, in %	1.0
tol_den		tolerance of the density result, in %	1.0
max_iter		maximum number of iterations to be performed, after which the function will throw a result	5
maxError		maximum error in the calls to getTemDen, in %	0.001
start_tem		, end_tem lower and upper limit of the explored temperature range	-1
start_den		, end_den lower and upper limit of the explored density range	-1
use_ANN		if True, an Analogic Neural Network will be used for the determination of Te and Ne. manage_RM from mwinai library is used. the hyper_parameters can be set-up with the self.ANN_inst_kwargs and self.ANN_init_kwargs dictionnaries. self.ANN_n_tem=30 and self.ANN_n_den=30 are the number of Te and Ne used to train. May also be changed before calling getCrossTemDen	False
limit_res		in case of using ANN, if limit_res, the tem and den values out of the start_tem, end_tem, start_den, end_den are set to np.nan. Otherwise, extrapolation is allowed.	False
ANN		if string, filename where to read AI4neb ANN. Otherwise, ANN is a manage_RM object. In both casesm the ANN needs to already be trained	None

Example:

tem, den = diags.getCrossTemDen('[OIII] 4363/5007', '[SII] 6731/6716', 0.0050, 1.0, 
            guess_tem=10000, tol_tem = 1., tol_den = 1., max_iter = 5)

Source code in pyneb/core/diags.py

def getCrossTemDen(self, diag_tem, diag_den, value_tem=None, value_den=None, obs=None, i_obs=None,
                   guess_tem=10000, tol_tem=1., tol_den=1., max_iter=5, maxError=1e-3,
                   start_tem= -1, end_tem= -1, start_den= -1, end_den= -1, use_ANN=False, 
                   limit_res=False, ANN=None):
    """
    Cross-converge the temperature and density derived from two sensitive line ratios, by inputting the quantity 
    derived with one line ratio into the other and then iterating.
    The temperature- and density-sensitive ratios can be input directly or as an Observation object

    Parameters:

        diag_tem:   temperature-sensitive diagnostic line ratio
        diag_den:   density-sensitive diagnostic line ratio
        value_tem:  value of the temperature-sensitive diagnostic
        value_den:  value of the density-sensitive diagnostic
        obs:        np.Observation object. Values for observed temperature and density diagnostics are
                        taken from it if value_tem and value_den are None
        i_obs:      index of the observations to be used from obs. 
                        If None, all the observations are considered.
                        May be a boolean array
        guess_tem:  temperature assumed in the first iteration, in K
        tol_tem:    tolerance of the temperature result, in %
        tol_den:    tolerance of the density result, in %
        max_iter:   maximum number of iterations to be performed, after which the function will throw a result
        maxError:   maximum error in the calls to getTemDen, in %
        start_tem:, end_tem  lower and upper limit of the explored temperature range 
        start_den:, end_den  lower and upper limit of the explored density range 
        use_ANN:    if True, an Analogic Neural Network will be used for the determination of Te and Ne.
                        manage_RM from mwinai library is used.
                        the hyper_parameters can be set-up with the self.ANN_inst_kwargs and
                        self.ANN_init_kwargs dictionnaries.
                        self.ANN_n_tem=30 and self.ANN_n_den=30 are the number of Te and Ne
                        used to train. May also be changed before calling getCrossTemDen
        limit_res:  in case of using ANN, if limit_res, the tem and den values out of the start_tem, end_tem,
                    start_den, end_den are set to np.nan. Otherwise, extrapolation is allowed.
        ANN:        if string, filename where to read AI4neb ANN. Otherwise, ANN is a manage_RM object.
                    In both casesm the ANN needs to already be trained

    **Example:**

        tem, den = diags.getCrossTemDen('[OIII] 4363/5007', '[SII] 6731/6716', 0.0050, 1.0, 
                    guess_tem=10000, tol_tem = 1., tol_den = 1., max_iter = 5)

    """
    # TODO:
    # Define a lower/upper density and temperature 
    # to be used in case the diag is pointing to values out of the boundaires
    if diag_tem not in self.diags:
        self.addDiag(diag_tem)
    atom_tem = self.diags[diag_tem][0]
    elem_tem, spec_tem, rec = parseAtom2(atom_tem)
    if atom_tem not in self.atomDict:
        self.atomDict[atom_tem] = pn.Atom(elem_tem, spec_tem, self.OmegaInterp, NLevels=self.NLevels)
    atom_tem = self.atomDict[atom_tem]
    if diag_den not in self.diags:
        self.addDiag(diag_den)
    atom_den = self.diags[diag_den][0]
    elem_den, spec_den, rec = parseAtom2(self.diags[diag_den][0])
    if (atom_den) not in self.atomDict:
        self.atomDict[atom_den] = pn.Atom(elem_den, spec_den, self.OmegaInterp, NLevels=self.NLevels)
    atom_den = self.atomDict[atom_den]
    eval_tem = self.diags[diag_tem][1]
    eval_den = self.diags[diag_den][1]
    calling = 'Diag.getCrossTemDen %s %s' % (diag_tem, diag_den)
    if value_tem is None:
        def L(wave):
            if i_obs is None:
                return obs.getLine(elem_tem, spec_tem, wave).corrIntens
            else:
                return obs.getLine(elem_tem, spec_tem, wave).corrIntens[i_obs]
        def I(i, j):
            wave = atom_tem.wave_Ang[i - 1, j - 1]
            if i_obs is None:
                return obs.getLine(elem_tem, spec_tem, wave).corrIntens
            else:
                return obs.getLine(elem_tem, spec_tem, wave).corrIntens[i_obs]
        def B(label, I=I, L=L):
            full_label = elem_tem + spec_tem + '_' + label
            if i_obs is None:
                corrIntens = obs.getLine(label=full_label).corrIntens
            else:
                corrIntens = obs.getLine(label=full_label).corrIntens[i_obs]
            return corrIntens
        #pn.log_.debug('to eval is {0}'.format(eval_tem), calling=calling + 'TEST')
        try:
            value_tem = eval(eval_tem)
            #pn.log_.debug('to eval = {0}'.format(value_tem), calling=calling + 'TEST')
        except:
            pn.log_.warn('No value for {0} {1}: {2} in obs'.format(elem_tem, spec_tem, diag_tem), calling=calling)
            return None
    else:
        if type(value_tem) == type([]): value_tem = np.asarray(value_tem)
    if value_den is None:
        def L(wave): 
            if i_obs is None:
                return obs.getLine(elem_den, spec_den, wave).corrIntens
            else:
                return obs.getLine(elem_den, spec_den, wave).corrIntens[i_obs]
        def I(i, j):
            wave = atom_den.wave_Ang[i - 1, j - 1]
            pn.log_.debug('wave is {0}'.format(wave), calling=calling + 'TEST3')
            if i_obs is None:
                return obs.getLine(elem_den, spec_den, wave).corrIntens
            else:
                return obs.getLine(elem_den, spec_den, wave).corrIntens[i_obs]
        def B(label, I=I, L=L):
            full_label = elem_den + spec_den + '_' + label
            if i_obs is None:
                corrIntens = obs.getLine(label=full_label).corrIntens
            else:
                corrIntens = obs.getLine(label=full_label).corrIntens[i_obs]
            return corrIntens
        #pn.log_.debug('to eval is {0}'.format(eval_den), calling=calling + ' TEST')
        try:
            value_den = eval(eval_den)
            #pn.log_.debug('to eval = {0}'.format(value_den), calling=calling + ' TEST1')
        except:
            pn.log_.warn('No value for {0} {1}: {2} in obs'.format(elem_den, spec_den, diag_den), calling=calling)
            return None
    else:
        if type(value_den) == type([]): value_den = np.asarray(value_den)
    if use_ANN:

        if config.INSTALLED['ai4neb']:
            from ai4neb import manage_RM

        if not config.INSTALLED['ai4neb']:
            self.log_.error('_getPopulations_ANN cannot be used in absence of ai4neb package',
                          calling=self.calling)
            return None
        if start_tem == -1:
            tem_min = 3000.
        else:
            tem_min = start_tem
        if end_tem == -1:
            tem_max = 20000.
        else:
            tem_max = end_tem
        if start_den == -1:
            den_min = 10.
        else:
            den_min = start_den
        if end_den == -1:
            den_max = 1e6
        else:
            den_max = end_den
        if ANN is None:
            # define emisGrid objects to generate Te-Ne emissionmaps
            tem_EG = pn.EmisGrid(atomObj=atom_tem, 
                                 n_tem=self.ANN_n_tem, n_den=self.ANN_n_den, 
                                 tem_min=tem_min, tem_max=tem_max,
                                 den_min=den_min, den_max=den_max)
            den_EG = pn.EmisGrid(atomObj=atom_den, 
                                 n_tem=self.ANN_n_tem, n_den=self.ANN_n_den, 
                                 tem_min=tem_min, tem_max=tem_max,
                                 den_min=den_min, den_max=den_max)
            # compute emission line ratio maps in the Te-Ne space
            tem_2D = tem_EG.getGrid(to_eval = eval_tem)
            den_2D = den_EG.getGrid(to_eval = eval_den)
            # X is a set of line ratio pairs
            X = np.array((tem_2D.ravel(), den_2D.ravel())).T
            # y is a set of corresponding Te-Ne pairs
            y = np.array((tem_EG.tem2D.ravel()/1e4, np.log10(den_EG.den2D.ravel()))).T
            # Instantiate, init and train the ANN
            self.ANN = manage_RM(X_train=X, y_train=y, **self.ANN_inst_kwargs)
            self.ANN.init_RM(**self.ANN_init_kwargs)
            self.ANN.train_RM()
        else:
            if type(ANN) is str:
                self.ANN = manage_RM(RM_filename=ANN)
            else:
                self.ANN = ANN
        # set the test values to the one we are looking for
        shape = value_tem.shape
        self.ANN.set_test(np.array((value_tem.ravel(), value_den.ravel())).T)
        # predict the result and denormalize them
        self.ANN.predict()
        if self.ANN.isfin is None:
            tem = self.ANN.pred[:,0]*1e4
            den = 10**self.ANN.pred[:,1]
        else:
            tem = np.zeros_like(value_tem.ravel()) * -10
            tem[self.ANN.isfin] = self.ANN.pred[:,0]*1e4
            den = np.zeros_like(value_tem.ravel()) * -10
            den[self.ANN.isfin] = 10**self.ANN.pred[:,1]
        tem = np.reshape(tem, shape)
        den = np.reshape(den, shape)
        if limit_res:
            mask = (tem<tem_min) | (tem>tem_max)
            tem[mask] = np.nan
            pn.log_.debug('Removing {} points out of Te range'.format(mask.sum()), calling='getCrossTemDen')
            mask = (den<den_min) | (den>den_max)
            den[mask] = np.nan
            pn.log_.debug('Removing {} points out of Ne range'.format(mask.sum()), calling='getCrossTemDen')
    else:
        den = atom_den.getTemDen(value_den, tem=guess_tem, to_eval=eval_den,
                                 maxError=maxError, start_x=start_den, end_x=end_den)

        tem = atom_tem.getTemDen(value_tem, den=den, to_eval=eval_tem,
                                 maxError=maxError, start_x=start_tem, end_x=end_tem)
#        self.log_.debug('tem: ' + str(tem) + ' den:' + str(den), calling='getCrossTemDen')
        no_conv = np.ones_like(den).astype(bool)
        n_tot = np.asarray(value_tem).size
        for i in np.arange(max_iter):
            if type(tem) == type(1.):
                tem_old = tem
            else:
                tem_old = tem.copy()
            if type(den) == type(1.):
                den_old = den
            else:
                den_old = den.copy()

            if n_tot > 1:
                den[no_conv] = atom_den.getTemDen(value_den[no_conv], tem=tem_old[no_conv],
                                                  to_eval=eval_den, start_x=start_den, end_x=end_den)
                tem[no_conv] = atom_tem.getTemDen(value_tem[no_conv], den=den_old[no_conv],
                                                  to_eval=eval_tem, start_x=start_tem, end_x=end_tem)
            else:
                den = atom_den.getTemDen(value_den, tem=tem_old, to_eval=eval_den, start_x=start_den, end_x=end_den)
                tem = atom_tem.getTemDen(value_tem, den=den_old, to_eval=eval_tem, start_x=start_tem, end_x=end_tem)

            no_conv = ((abs(den_old - den) / den * 100) > tol_den) | ((abs(tem_old - tem) / tem * 100) > tol_tem)
            if type(no_conv) == type(True):
                n_no_conv = int(no_conv)
            else:
                n_no_conv = no_conv.sum()

            pn.log_.message('{0} (max={1}): not converged {2} of {3}.'.format(i, max_iter, n_no_conv, n_tot),
                            calling=calling)
            if n_no_conv == 0:
                return tem, den
        if n_tot == 1:
            tem = np.nan
            den = np.nan
        else:
            tem[no_conv] = np.nan
            den[no_conv] = np.nan
    return tem, den

getDiagFromLabel(label)

Parameters:

Name	Type	Description	Default
label(str)		a diagnostic label	required

Returns:

Type	Description
	Return the definition of a diagnostic (the 3 or 4 elements tuple)

Usage:

diags.getDiagFromLabel('[NII] 5755/6548')

Source code in pyneb/core/diags.py

def getDiagFromLabel(self, label):
    """
    Parameters:
        label(str): a diagnostic label 

    Returns:
        Return the definition of a diagnostic (the 3 or 4 elements tuple)

    **Usage:**

        diags.getDiagFromLabel('[NII] 5755/6548') 
    """
    if label in self.diags:
        return self.diags[label]
    else:
        self.log_.warn('{0} not in diagnostic'.format(label), calling=self.calling)
        return None

getDiagLabels()

Return the labels of all the diagnostics defined in the Object

Source code in pyneb/core/diags.py

def getDiagLabels(self):
    """Return the labels of all the diagnostics defined in the Object"""
    return self.diags.keys()    

getDiagLimits(diag)

Return the low and high density values for a given diagnostic

Source code in pyneb/core/diags.py

def getDiagLimits(self, diag):
    """Return the low and high density values for a given diagnostic
    """
    atom = self.atomDict[self.diags[diag][0]]
    to_eval = self.diags[diag][1]
    HDR = atom.getHighDensRatio(to_eval = to_eval)
    LDR = atom.getLowDensRatio(to_eval = to_eval)
    return(np.sort((LDR, HDR)))

getDiags()

Return the definitions (tuples) of all the diagnostics defined in the Object.

Source code in pyneb/core/diags.py

def getDiags(self):
    """Return the definitions (tuples) of all the diagnostics defined in the Object.
    """
    return [self.getDiagFromLabel(label) for label in self.diags]

getUniqueAtoms()

Return a numpy.ndarray of the ions needed by the diagnostics. Unique is applied to the list before returning.

Returns:

Type	Description
np.ndarray	Ions needed by the diagnostics

Source code in pyneb/core/diags.py

def getUniqueAtoms(self):
    """Return a numpy.ndarray of the ions needed by the diagnostics. Unique is applied to the list before returning.

    Returns:
        (np.ndarray): Ions needed by the diagnostics
    """
    return np.unique([self.diags[d][0] for d in self.diags if self.diags[d] is not None])

plot(emis_grids, obs, quad=True, i_obs=None, alpha=0.3, ax=None, error_band=True, return_CS=False, col_dic={'C': 'cyan', 'N': 'blue', 'O': 'green', 'Ne': 'magenta', 'Ar': 'red', 'Cl': 'magenta', 'S': 'black', 'Fe': 'blue'})

Plotting tool to generate Te-Ne diagrams.

Parameters:

Name	Type	Description	Default
emis_grids		A dictionary of EmisGrid objects refereed by their atom strings (e.g. 'O3') This can for example be the output of pn.getEmisGridDict()	required
obs		A pn.Observation object that is supposed to contain the line intensities used for the plot (corrected intensities).	required
quad		If True (default) the sum of the error is quadratic,otherwise is linear.	True
i_obs		reference for the observation to be plotted, in case there is more than one in the obs object	None
alpha		Transparency for the error bands in the plot	0.3
error_band		Boolean: plot [default] an error band	True

Usage: diags.plot(emisgrids, obs, i_obs=3)

Source code in pyneb/core/diags.py

def plot(self, emis_grids, obs, quad=True, i_obs=None, alpha=0.3, ax=None, error_band=True,
return_CS=False, 
         col_dic={'C':'cyan', 'N':'blue', 'O':'green', 'Ne':'magenta',
                  'Ar':'red', 'Cl':'magenta', 'S':'black', 'Fe':'blue'}):
    """Plotting tool to generate Te-Ne diagrams.

    Parameters:
        emis_grids:    A dictionary of EmisGrid objects refereed by their atom strings (e.g. 'O3')
                        This can for example be the output of pn.getEmisGridDict()
        obs:           A pn.Observation object that is supposed to contain the line intensities
                        used for the plot (corrected intensities).
        quad:          If True (default) the sum of the error is quadratic,otherwise is linear.
        i_obs:         reference for the observation to be plotted, in case there is more than one
                        in the obs object
        alpha:         Transparency for the error bands in the plot
        error_band:    Boolean: plot [default] an error band

        return_CS     [False] If True, return a list of the contour plots
        col_dic       Colors for the different ions.            
    **Usage:**
        diags.plot(emisgrids, obs, i_obs=3)
    """
    if not pn.config.INSTALLED['plt']: 
        pn.log_.error('Matplotlib not available, no plot', calling=self.calling + '.plot')
        return None
    else:
        import matplotlib.pyplot as plt
    if type(emis_grids) != type({}):
        self.log_.error('the first parameter must be a dictionary', calling=self.calling + '.plot')
        return None
    for em in emis_grids:
        if not isinstance(emis_grids[em], pn.EmisGrid):
            self.log_.error('the first parameter must be a dictionary of EmisGrid', calling=self.calling + '.plot')
            return None
    if not isinstance(obs, pn.Observation):
        self.log_.error('the second parameter must be an Observation', calling=self.calling + '.plot')
        return None
    if (i_obs is None) and (obs.n_obs != 1):
        self.log_.error('i_obs must be specified when obs is multiple. try i_obs=0', calling=self.calling)
        return None
    if ax is None:
        f, ax = plt.subplots()
    else:
        f = plt.gcf()
    X = np.log10(emis_grids[list(emis_grids.keys())[0]].den2D)
    Y = emis_grids[list(emis_grids.keys())[0]].tem2D
    CSs = []
    for label in self.diags:
        self.log_.message('plotting {0}'.format(label), calling=self.calling)
        diag = self.diags[label]
        atom = diag[0]
        def I(i, j):
            return emis_grids[atom].getGrid(lev_i=i, lev_j=j)
        def L(wave):
            return emis_grids[atom].getGrid(wave=wave)
        def S(label):
            return emis_grids[atom].getGrid(label=label)
        def B(label, I=I, L=L):
            full_label = atom + '_' + label
            if full_label in BLEND_LIST:
                to_eval = BLEND_LIST[full_label]
            else:
                self.log_.warn('{0} not in BLEND_LIST'.format(full_label), calling=self.calling)
                return None
            return eval(to_eval)
        diag_map = eval(diag[1])
        try:
            diag_map = eval(diag[1])
        except:
            diag_map = None
            self.log_.warn('diag {0} {1} not used'.format(diag[0], diag[1]), calling=self.calling)
        if diag_map is not None:
            sym, spec, rec = parseAtom2(atom)
            def I(i, j):
                wave = emis_grids[atom].atom.wave_Ang[i - 1, j - 1]
                corrIntens = obs.getLine(sym, spec, wave).corrIntens
                if i_obs is None:
                    return corrIntens
                else:
                    return corrIntens[i_obs]
            def L(wave):
                corrIntens = obs.getLine(sym, spec, wave).corrIntens
                if i_obs is None:
                    return corrIntens
                else:
                    return corrIntens[i_obs]
            def S(label):
                full_label = atom + '_' + label + 'A'
                corrIntens = obs.getLine(label=full_label).corrIntens
                if i_obs is None:
                    return corrIntens
                else:
                    return corrIntens[i_obs]                    
            def B(label):
                full_label = atom + '_' + label
                corrIntens = obs.getLine(label=full_label).corrIntens
                if i_obs is None:
                    return corrIntens
                else:
                    return corrIntens[i_obs]
            try:
                ee = eval(diag[1])
                if type(ee) == np.float64:
                    diag_value = ee
                else:
                    diag_value = ee[0]
            except:
                #print(ee, type(ee))
                diag_value = None
                self.log_.warn('A line in diagnostic {0} of {1}{2} is missing'.format(diag[1], sym, spec),
                               calling=self.calling)
            if diag_value is not None and not np.isnan(diag_value) and not np.isinf(diag_value):
                def E(wave):
                    err = obs.getLine(sym, spec, wave).corrError
                    if i_obs is None:
                        return err
                    else:
                        return err[i_obs]
                def BE(label):
                    full_label = atom + '_' + label
                    err = obs.getLine(label=full_label).corrError
                    if i_obs is None:
                        return err
                    else:
                        return err[i_obs]
                def SE(label):
                    full_label = atom + '_' + label + 'A'
                    err = obs.getLine(label=full_label).corrError
                    if i_obs is None:
                        return err
                    else:
                        return err[i_obs]                        
                if quad is True:
                    RMS = lambda err: np.sqrt((np.asarray(err) ** 2.).sum())
                else:
                    RMS = lambda err: (np.asarray(err)).sum() 
                tol_value = eval(diag[2])
                if sym in col_dic:
                    col = col_dic[sym]
                else:
                    col = 'black'
                style_dic = {'1':'-', '2':'--', '3':':', '4':'-.', '5':'-', '6':'--'}
                style = style_dic[spec]
                if tol_value > 0. and error_band:
                    levels = [(1 - tol_value) * diag_value, (1 + tol_value) * diag_value]
                    if levels[0] < levels[1]:
                        #pn.log_.debug('{} levels {}'.format(label, levels), calling=self.calling)
                        CS = ax.contourf(X, Y, diag_map, levels=levels, alpha=alpha, colors=col)
                        CSs.append(CS)
                CS = ax.contour(X, Y, diag_map, levels=[diag_value], colors=col, linestyles=style)
                CSs.append(CS)
                try:
                    ax.set_xlabel(r'log(n$_{\rm e}$) [cm$^{-3}$]')
                    ax.set_ylabel(r'T$_{\rm e}$ [K]')
                    if len(diag) >= 4:
                        fmt = diag[3]
                    else:
                        fmt = '[{0}{1}]'.format(sym, int_to_roman(int(spec)))
                    ax.clabel(CS, inline=True, fmt=fmt, fontsize=15, colors=col)
                    if type(diag_value) is np.ndarray:
                        diag_value = diag_value[0]
                    self.log_.message('plotted {0}: {1} = {2:.2} with error of {3:.2} %'.format(fmt, label, diag_value, tol_value * 100),
                                      calling=self.calling)
                except:
                    self.log_.message('NOT plotted {0} {1}'.format(fmt, label),
                                      calling=self.calling)

    if return_CS:
        return CSs

setAtoms(atom_dic)

Define the dictionary containing the atoms used for the diagnostics.

A dictionary of atom instantiations refereed by atom strings, for

Parameters:

Name	Type	Description	Default
atom_dic		a dictionary of Atom instances, indexed by atom strings. Example: {'O3': pn.Atom('O', 3)}	required

Source code in pyneb/core/diags.py

def setAtoms(self, atom_dic):
    """Define the dictionary containing the atoms used for the diagnostics.

    A dictionary of atom instantiations refereed by atom strings, for 

    Parameters:

        atom_dic : a dictionary of Atom instances, indexed by atom strings.
                **Example:**

                    {'O3': pn.Atom('O', 3)}

    """
    if type(atom_dic) != type({}):
        self.log_.error('the parameter must be a dictionary.', calling=self.calling + '.setAtoms')
        return None
    for atom in atom_dic:
        if not isinstance(atom_dic[atom], pn.Atom) and not isinstance(atom_dic[atom], pn.RecAtom):
            self.log_.error('the parameter must be a dictionary of Atom.', calling=self.calling + '.setAtoms')
            return None
    self.atomDict = atom_dic