<commit_before>import numpy as np
import h5py
from sklearn.decomposition import PCA

from Starfish.grid_tools import HDF5Interface, determine_chunk_log
from Starfish.covariance import sigma, V12, V22
from Starfish import constants as C

class PCAGrid:
    '''
    Create and query eigenspectra.
    '''

    def __init__(self, wl, dv, flux_mean, flux_std, eigenspectra, w, gparams):
        '''

        :param wl: wavelength array
        :type wl: 1D np.array
        :param dv: delta velocity
        :type dv: float
        :param flux_mean: mean flux spectrum
        :type flux_mean: 1D np.array
        :param flux_std: standard deviation flux spectrum
        :type flux_std: 1D np.array
        :param eigenspectra: the principal component eigenspectra
        :type eigenspectra: 2D np.array
        :param w: weights to reproduce any spectrum in the original grid
        :type w: 1D np.array
        :param gparams: The stellar parameters of the synthetic library
        :type gparams: 1D list of dictionaries


        '''
        self.wl = wl
        self.dv = dv
        self.flux_mean = flux_mean
        self.flux_std = flux_std
        self.eigenspectra = eigenspectra
        self.ncomp = len(self.eigenspectra)
        self.w = w
        self.gparams = gparams

        self.npix = len(self.wl)
        self.m = self.w.shape[1] # The number of spectra in the synthetic grid


    @classmethod
    def create(cls, interface, ncomp):
        '''
        Create a PCA grid object from a synthetic spectral library, with
        configuration options specified in a dictionary.

        :param interface: HDF5Interface containing the instrument-processed spectra.
        :type interface: HDF5Interface
        :param ncomp: number of eigenspectra to keep
        :type ncomp: int

        '''

        wl = interface.wl
        dv = interface.dv

        npix = len(wl)
        # number of spectra in the synthetic library
        m = len(interface.list_grid_points)

        fluxes = np.empty((m, npix))

        z = 0
        for i, spec in enumerate(interface.fluxes):
            fluxes[z,:] = spec
            z += 1

        # Normalize all of the fluxes to an average value of 1
        # In order to remove uninteresting correlations
        fluxes = fluxes/np.average(fluxes, axis=1)[np.newaxis].T

        # Subtract the mean from all of the fluxes.
        flux_mean = np.average(fluxes, axis=0)
        fluxes -= flux_mean

        # "Whiten" each spectrum such that the variance for each wavelength is 1
        flux_std = np.std(fluxes, axis=0)
        fluxes /= flux_std

        # Use the scikit-learn PCA module
        pca = PCA()
        pca.fit(fluxes)
        comp = pca.transform(fluxes)
        components = pca.components_
        mean = pca.mean_
        print("Shape of PCA components {}".format(components.shape))

        if not np.allclose(mean, np.zeros_like(mean)):
            import sys
            sys.exit("PCA mean is more than just numerical noise. Something's wrong!")
            # Otherwise, the PCA mean is just numerical noise that we can ignore.

        print("Keeping only the first {} components".format(ncomp))
        eigenspectra = components[0:ncomp]

        gparams = np.empty((m, 3))
        z = 0
        for i, params in enumerate(interface.list_grid_points):
            if i == test_index:
                test_params = np.array([params["temp"], params["logg"], params["Z"]])
                continue
            gparams[z, :] = np.array([params["temp"], params["logg"], params["Z"]])
            z += 1

        #Create w, the weights corresponding to the synthetic grid

        w = np.empty((ncomp, m))
        for i,pcomp in enumerate(eigenspectra):
            for j,spec in enumerate(fluxes):
                w[i,j] = np.sum(pcomp * spec)

        return cls(wl, dv, flux_mean, flux_std, eigenspectra, w, gparams)

    def write(self, filename):
        '''
        Write the PCA decomposition to an HDF5 file.

        :param filename: name of the HDF5 to create.
        :type filename: str

        '''

        hdf5 = h5py.File(filename, "w")

        hdf5.attrs["dv"] = self.dv

        pdset = hdf5.create_dataset("eigenspectra", (self.ncomp + 3, self.npix),
            compression='gzip', dtype="f8", compression_opts=9)<commit_msg>Removed test_index<commit_after>import numpy as np
import h5py
from sklearn.decomposition import PCA

from Starfish.grid_tools import HDF5Interface, determine_chunk_log
from Starfish.covariance import sigma, V12, V22
from Starfish import constants as C

class PCAGrid:
    '''
    Create and query eigenspectra.
    '''

    def __init__(self, wl, dv, flux_mean, flux_std, eigenspectra, w, gparams):
        '''

        :param wl: wavelength array
        :type wl: 1D np.array
        :param dv: delta velocity
        :type dv: float
        :param flux_mean: mean flux spectrum
        :type flux_mean: 1D np.array
        :param flux_std: standard deviation flux spectrum
        :type flux_std: 1D np.array
        :param eigenspectra: the principal component eigenspectra
        :type eigenspectra: 2D np.array
        :param w: weights to reproduce any spectrum in the original grid
        :type w: 1D np.array
        :param gparams: The stellar parameters of the synthetic library
        :type gparams: 1D list of dictionaries


        '''
        self.wl = wl
        self.dv = dv
        self.flux_mean = flux_mean
        self.flux_std = flux_std
        self.eigenspectra = eigenspectra
        self.ncomp = len(self.eigenspectra)
        self.w = w
        self.gparams = gparams

        self.npix = len(self.wl)
        self.m = self.w.shape[1] # The number of spectra in the synthetic grid


    @classmethod
    def create(cls, interface, ncomp):
        '''
        Create a PCA grid object from a synthetic spectral library, with
        configuration options specified in a dictionary.

        :param interface: HDF5Interface containing the instrument-processed spectra.
        :type interface: HDF5Interface
        :param ncomp: number of eigenspectra to keep
        :type ncomp: int

        '''

        wl = interface.wl
        dv = interface.dv

        npix = len(wl)
        # number of spectra in the synthetic library
        m = len(interface.list_grid_points)

        fluxes = np.empty((m, npix))

        z = 0
        for i, spec in enumerate(interface.fluxes):
            fluxes[z,:] = spec
            z += 1

        # Normalize all of the fluxes to an average value of 1
        # In order to remove uninteresting correlations
        fluxes = fluxes/np.average(fluxes, axis=1)[np.newaxis].T

        # Subtract the mean from all of the fluxes.
        flux_mean = np.average(fluxes, axis=0)
        fluxes -= flux_mean

        # "Whiten" each spectrum such that the variance for each wavelength is 1
        flux_std = np.std(fluxes, axis=0)
        fluxes /= flux_std

        # Use the scikit-learn PCA module
        pca = PCA()
        pca.fit(fluxes)
        comp = pca.transform(fluxes)
        components = pca.components_
        mean = pca.mean_
        print("Shape of PCA components {}".format(components.shape))

        if not np.allclose(mean, np.zeros_like(mean)):
            import sys
            sys.exit("PCA mean is more than just numerical noise. Something's wrong!")
            # Otherwise, the PCA mean is just numerical noise that we can ignore.

        print("Keeping only the first {} components".format(ncomp))
        eigenspectra = components[0:ncomp]

        gparams = np.empty((m, 3))
        z = 0
        for i, params in enumerate(interface.list_grid_points):
            # if i == test_index:
                # test_params = np.array([params["temp"], params["logg"], params["Z"]])
                # continue
            gparams[z, :] = np.array([params["temp"], params["logg"], params["Z"]])
            z += 1

        #Create w, the weights corresponding to the synthetic grid

        w = np.empty((ncomp, m))
        for i,pcomp in enumerate(eigenspectra):
            for j,spec in enumerate(fluxes):
                w[i,j] = np.sum(pcomp * spec)

        return cls(wl, dv, flux_mean, flux_std, eigenspectra, w, gparams)

    def write(self, filename):
        '''
        Write the PCA decomposition to an HDF5 file.

        :param filename: name of the HDF5 to create.
        :type filename: str

        '''

        hdf5 = h5py.File(filename, "w")

        hdf5.attrs["dv"] = self.dv

        pdset = hdf5.create_dataset("eigenspectra", (self.ncomp + 3, self.npix),
            compression='gzip', dtype="f8", compression_opts=9)
{'fraction_non_alphanumeric': 0.09118224893100658, 'fraction_numerical': 0.00808967987981047, 'mean_word_length': 2.154939846882975, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6396146', 'n_tokens_mistral': 2680, 'n_tokens_neox': 2484, 'n_words': 818}