Merge pull request #31 from gaia-dpci/cholesky-addition

Cholesky addition

Author: druzm

.github/workflows/python-package.yml

@@ -39,4 +39,4 @@ jobs:
         flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics
     - name: Test with pytest
       run: |
-        pytest
+        pytest -s

.gitignore

@@ -11,4 +11,8 @@ __pycache__/
 build/
 dist/
 docs/_build/
+htmlcov/
+junit/
 tests_output_files/
+coverage.xml
+query_tests.py

MANIFEST.in

@@ -1,5 +1,5 @@
 include LICENCE
 include requirements.txt
-include gaiaxpy/headers/headers_dict.txt
+include gaiaxpy/output/ecsv_headers/headers_dict.txt
 graft gaiaxpy/config
 global-exclude __pycache__

docs/source/releasenotes.rst

@@ -14,8 +14,8 @@ Version 1.1.3
 -------------
 Released on 2022/06/16.
 
-* Fix legend bug in plotter.
-* Restrict pandas version, >= 1.0.0.
+* Fixed legend bug in plotter.
+* Restricted pandas version, >= 1.0.0.
 
 Version 1.1.2
 -------------

gaiaxpy/__init__.py

@@ -1,9 +1,10 @@
 from .calibrator.calibrator import calibrate
+from .cholesky.cholesky import get_inverse_covariance_matrix, get_chi2
 from .converter.converter import convert
-from .core import pwl_to_wl, wl_to_pwl, pwl_range, wl_range
-from .error_correction import apply_error_correction
+from .core.dispersion_function import pwl_to_wl, wl_to_pwl, pwl_range, wl_range
+from .error_correction.error_correction import apply_error_correction
 from .generator.generator import generate
 from .generator.photometric_system import PhotometricSystem
 from .plotter.plot_spectra import plot_spectra
 
-__version__ = '1.1.4'
+__version__ = '1.2.0'

gaiaxpy/calibrator/__init__.py

@@ -1,5 +0,0 @@
-from . import calibrator
-from .calibrator import _calibrate, _create_merge, _create_spectrum
-
-from . import external_instrument_model
-from .external_instrument_model import ExternalInstrumentModel

gaiaxpy/calibrator/calibrator.py

@@ -8,21 +8,27 @@
 import pandas as pd
 from configparser import ConfigParser
 from pathlib import Path
+from tqdm import tqdm
 from os.path import join
 from .external_instrument_model import ExternalInstrumentModel
-from gaiaxpy.config import config_path
-from gaiaxpy.core.satellite import BANDS
-from gaiaxpy.core import _get_spectra_type, _load_xpmerge_from_csv, \
-                         _load_xpsampling_from_csv, _progress_tracker, \
-                         _validate_arguments, _validate_wl_sampling, satellite
-from gaiaxpy.input_reader import InputReader
-from gaiaxpy.output import SampledSpectraData
-from gaiaxpy.spectrum import _get_covariance_matrix, AbsoluteSampledSpectrum, \
-                             SampledBasisFunctions, XpContinuousSpectrum
+from gaiaxpy.config.paths import config_path
+from gaiaxpy.core.config import _load_xpmerge_from_csv, _load_xpsampling_from_csv
+from gaiaxpy.core.generic_functions import cast_output, _get_spectra_type, \
+                                           _validate_arguments, \
+                                           _validate_wl_sampling
+from gaiaxpy.core.satellite import BANDS, BP_WL, RP_WL
+from gaiaxpy.core.generic_variables import pbar_colour, pbar_units
+from gaiaxpy.input_reader.input_reader import InputReader
+from gaiaxpy.output.sampled_spectra_data import SampledSpectraData
+from gaiaxpy.spectrum.absolute_sampled_spectrum import AbsoluteSampledSpectrum
+from gaiaxpy.spectrum.sampled_basis_functions import SampledBasisFunctions
+from gaiaxpy.spectrum.utils import _get_covariance_matrix
+from gaiaxpy.spectrum.xp_continuous_spectrum import XpContinuousSpectrum
 
 config_parser = ConfigParser()
 config_parser.read(join(config_path, 'config.ini'))
-
+tqdm.pandas(desc='Processing data', unit=pbar_units['calibrator'], leave=False, \
+            colour=pbar_colour) # Activate tqdm for pandas
 
 def calibrate(
         input_object,
@@ -107,18 +113,16 @@ def _calibrate(
     """
     _validate_wl_sampling(sampling)
     _validate_arguments(_calibrate.__defaults__[3], output_file, save_file)
-    parsed_input_data, extension = InputReader(input_object, _calibrate, username, password)._read()
+    parsed_input_data, extension = InputReader(input_object, _calibrate, \
+                                               username, password)._read()
     label = 'calibrator'
-
-    xp_design_matrices, xp_merge = _generate_xp_matrices_and_merge(label, sampling, bp_model, rp_model)
+    xp_design_matrices, xp_merge = _generate_xp_matrices_and_merge(label, \
+                                   sampling, bp_model, rp_model)
     # Create sampled basis functions
-    spectra_list = _create_spectra(parsed_input_data, truncation, xp_design_matrices, xp_merge)
-    # Generate output
-    spectra_df = pd.DataFrame.from_records([spectrum._spectrum_to_dict() for spectrum in spectra_list])
-    spectra_type = _get_spectra_type(spectra_list)
-    spectra_df.attrs['data_type'] = spectra_type
-    positions = spectra_list[0]._get_positions()
+    spectra_df, positions = _create_spectra(parsed_input_data, truncation, \
+                                            xp_design_matrices, xp_merge)
     output_data = SampledSpectraData(spectra_df, positions)
+    output_data.data = cast_output(output_data)
     # Save output
     Path(output_path).mkdir(parents=True, exist_ok=True)
     output_data.save(save_file, output_path, output_file, output_format, extension)
@@ -136,8 +140,8 @@ def _create_merge(xp, sampling):
     Returns:
         dict: A dictionary containing a BP and an RP array with weights.
     """
-    wl_high = satellite.BP_WL.high
-    wl_low = satellite.RP_WL.low
+    wl_high = BP_WL.high
+    wl_low = RP_WL.low
 
     if xp == BANDS.bp:
         weight = np.array([1.0 if wl < wl_low else 0.0 if wl > wl_high else (
@@ -184,22 +188,15 @@ def _get_file_for_xp(xp, key, bp_model=bp_model, rp_model=rp_model):
 
 
 def _create_spectra(parsed_spectrum_file, truncation, design_matrices, merge):
-    """
-    Internal wrapper function. Allows _create_spectrum to use the generic
-    progress tracker.
-    """
-    spectra_list = []
     nrows = len(parsed_spectrum_file)
-
-    @_progress_tracker
-    def create_spectrum(row, *args):
-        truncation, design_matrices, merge = args[:3]
-        spectrum = _create_spectrum(
-            row, truncation, design_matrices, merge)
-        spectra_list.append(spectrum)
-    for index, row in parsed_spectrum_file.iterrows():
-        create_spectrum(row, truncation, design_matrices, merge, index, nrows)
-    return spectra_list
+    spectra_series = parsed_spectrum_file.progress_apply(lambda row: \
+                 _create_spectrum(row, truncation, design_matrices, merge), axis=1)
+    positions = spectra_series.iloc[0]._get_positions()
+    spectra_type = _get_spectra_type(spectra_series.iloc[0])
+    spectra_series = spectra_series.map(lambda x: x._spectrum_to_dict())
+    spectra_df = pd.DataFrame(spectra_series.tolist())
+    spectra_df.attrs['data_type'] = spectra_type
+    return spectra_df, positions
 
 
 def _create_spectrum(row, truncation, design_matrix, merge):
@@ -224,9 +221,6 @@ def _create_spectrum(row, truncation, design_matrix, merge):
     source_id = row['source_id']
     cont_dict = {}
     # Split both bands
-    source_id = row['source_id']
-    cont_dict = {}
-    # Split both bands
     for band in BANDS:
         try:
             covariance_matrix = _get_covariance_matrix(row, band)

gaiaxpy/calibrator/external_instrument_model.py

@@ -7,7 +7,7 @@
 
 import numpy as np
 from scipy import interpolate
-from gaiaxpy.file_parser import GenericParser
+from gaiaxpy.file_parser.parse_generic import GenericParser
 
 
 class ExternalInstrumentModel(object):

gaiaxpy/cholesky/__init__.py

@@ -0,0 +1,75 @@
+import pandas as pd
+from numpy import diag, dot, identity
+from scipy.linalg import cholesky, solve_triangular
+from gaiaxpy.core.satellite import BANDS
+from gaiaxpy.input_reader.input_reader import InputReader
+
+
+def __get_dot_product(L_inv):
+    try:
+        return dot(L_inv.T, L_inv)
+    except AttributeError:
+        return None
+
+
+def __get_inv_cholesky_decomp_lower(xp_errors, xp_correlation_matrix):
+    try:
+        L = cholesky(xp_correlation_matrix, lower=True)
+        # Invert lower triangular matrix.
+        L_inv = solve_triangular(L, identity(len(L)), lower=True)
+        # Matrix of inverse errors.
+        E_inv = diag(1.0 / xp_errors)
+        return dot(L_inv, E_inv)
+    except ValueError:
+        return None
+
+
+def get_inverse_covariance_matrix(input_object, band=None):
+    """
+    Compute the inverse covariance matrix.
+
+    Args:
+        input_object (object): Path to the file containing the mean spectra
+            as downloaded from the archive in their continuous representation,
+            a list of sources ids (string or long), or a pandas DataFrame.
+        band (str): Chosen band: 'bp' or 'rp'. If no band is passed, the function
+            will compute the inverse covariance for both 'bp' and 'rp''
+
+    Returns:
+        DataFrame or ndarray of ndarrays: DataFrame containing the source IDs and
+            the output inverse covariance matrices for the sources in the input
+            object if it contains more than one source or no band is passed to
+            the function.
+            The function will return a ndarray (of shape (55, 55)) if there is
+            only one source ID in the input data and a single band is selected.
+    """
+    parsed_input_data, extension = InputReader(input_object, get_inverse_covariance_matrix)._read()
+    bands_output = []
+    if band is None:
+        bands_to_process = BANDS
+        output_columns = ['source_id', 'bp_inverse_covariance', 'rp_inverse_covariance']
+    else:
+        bands_to_process = [band]
+        output_columns = ['source_id', f'{band}_inverse_covariance']
+    for b in bands_to_process:
+        xp_errors = parsed_input_data[f'{b}_coefficient_errors']
+        xp_correlation_matrix = parsed_input_data[f'{b}_coefficient_correlations']
+        L_inv_iterable = map(__get_inv_cholesky_decomp_lower, xp_errors, xp_correlation_matrix)
+        band_output = map(__get_dot_product, L_inv_iterable)
+        bands_output.append(band_output)
+    output_list = [parsed_input_data['source_id']]
+    for element in bands_output:
+        output_list.append(element)
+    output_df = pd.DataFrame(zip(*output_list), columns=output_columns)
+    if len(bands_to_process) == 1 and len(output_df) == 1:
+        return output_df[f'{band}_inverse_covariance'].iloc[0]
+    else:
+        return output_df
+
+def get_chi2(L_inv, residuals):
+    if L_inv.shape != (55, 55):
+        raise ValueError('Inverse covariance matrix shape must be (55, 55).')
+    if residuals.shape != (55,):
+        raise ValueError('Residuals shape must be (55,).')
+    x = dot(L_inv.T, residuals)
+    return dot(x.T, x)