Repository: StingraySoftware/stingray

Issue:
Optimize Stingray
I would like to discuss here the bottlenecks of the library, and all techniques that can be used to speed-up the processing.
We should start from numbers. Let us make use of the `% prun` magic on ipython to profile the significant bits of code (I would personally start from `Lightcurve` objects with >10^6 bins), together with visualization tools like `snakeviz`, to get them.
After, we can talk about strategies. I would absolutely avoid talking about generic optimization strategies when we do not have an idea of the real bottlenecks.

Discussion:

Patch:
diff --git a/setup.cfg b/setup.cfg
index cc901abd..3a82e5dc 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -40,6 +40,8 @@ all =
     pyfftw
     h5py
     pint-pulsar
+    zarr
+    numcodecs
 test =
     pytest<6
     pytest-astropy
diff --git a/stingray/__init__.py b/stingray/__init__.py
index 007b6a35..d2c3ecf3 100644
--- a/stingray/__init__.py
+++ b/stingray/__init__.py
@@ -19,3 +19,4 @@ if not _ASTROPY_SETUP_:
     from stingray.stats import *
     from stingray.bispectrum import *
     from stingray.varenergyspectrum import *
+    from stingray.largememory import *
diff --git a/stingray/crossspectrum.py b/stingray/crossspectrum.py
index c8060b06..14e53a01 100644
--- a/stingray/crossspectrum.py
+++ b/stingray/crossspectrum.py
@@ -1,12 +1,20 @@
-
+import copy
 import warnings
 from collections.abc import Iterable, Iterator
+
 import numpy as np
 import scipy
-import scipy.stats
-import scipy.fftpack
 import scipy.optimize
-import copy
+import scipy.stats
+
+from stingray.exceptions import StingrayError
+from stingray.gti import bin_intervals_from_gtis, check_gtis, cross_two_gtis
+from stingray.largememory import createChunkedSpectra, saveData
+from stingray.utils import genDataPath, rebin_data, rebin_data_log, simon
+
+from .events import EventList
+from .lightcurve import Lightcurve
+from .utils import show_progress
 
 # location of factorial moved between scipy versions
 try:
@@ -15,30 +23,15 @@ except ImportError:
     from scipy.special import factorial
 
 try:
-    import pyfftw
     from pyfftw.interfaces.scipy_fftpack import fft, fftfreq
 except ImportError:
     warnings.warn("Using standard scipy fft")
     from scipy.fftpack import fft, fftfreq
 
-
-from stingray.lightcurve import Lightcurve
-from stingray.utils import rebin_data, simon, rebin_data_log
-from stingray.exceptions import StingrayError
-from stingray.gti import cross_two_gtis, bin_intervals_from_gtis, check_gtis
-from .events import EventList
-from .utils import show_progress
-
-
-try:
-    from tqdm import tqdm as show_progress
-except ImportError:
-    def show_progress(a, **kwargs):
-        return a
-
-__all__ = ["Crossspectrum", "AveragedCrossspectrum",
-           "coherence", "time_lag", "cospectra_pvalue",
-            "normalize_crossspectrum"]
+__all__ = [
+    "Crossspectrum", "AveragedCrossspectrum", "coherence", "time_lag",
+    "cospectra_pvalue", "normalize_crossspectrum"
+]
 
 
 def normalize_crossspectrum(unnorm_power, tseg, nbins, nphots1, nphots2, norm="none", power_type="real"):
@@ -107,7 +100,6 @@ def normalize_crossspectrum(unnorm_power, tseg, nbins, nphots1, nphots2, norm="n
     return power
 
 
-
 def _averaged_cospectra_cdf(xcoord, n):
     """
     Function calculating the cumulative distribution function for
@@ -139,8 +131,8 @@ def _averaged_cospectra_cdf(xcoord, n):
                 n - 1 - j) * factorial(j)
             prefac_bottom3 = 2.0 ** (n + j)
 
-            prefac = prefac_top / (
-            prefac_bottom1 * prefac_bottom2 * prefac_bottom3)
+            prefac = prefac_top / (prefac_bottom1 * prefac_bottom2 *
+                                   prefac_bottom3)
 
             gf = -j + n
 
@@ -210,7 +202,7 @@ def cospectra_pvalue(power, nspec):
     if not np.all(np.isfinite(power)):
         raise ValueError("power must be a finite floating point number!")
 
-    #if power < 0:
+    # if power < 0:
     #    raise ValueError("power must be a positive real number!")
 
     if not np.isfinite(nspec):
@@ -690,7 +682,6 @@ class Crossspectrum(object):
             unnorm_power, tseg, self.n, self.nphots1, self.nphots2, self.norm,
             self.power_type)
 
-
     def rebin_log(self, f=0.01):
         """
         Logarithmic rebin of the periodogram.
@@ -825,13 +816,24 @@ class Crossspectrum(object):
         except ImportError:
             raise ImportError("Matplotlib required for plot()")
 
-        fig = plt.figure('crossspectrum')
-        fig = plt.plot(self.freq, np.abs(self.power), marker, color='b',
-                       label='Amplitude')
-        fig = plt.plot(self.freq, np.abs(self.power.real), marker, color='r',
-                       alpha=0.5, label='Real Part')
-        fig = plt.plot(self.freq, np.abs(self.power.imag), marker, color='g',
-                       alpha=0.5, label='Imaginary Part')
+        plt.figure('crossspectrum')
+        plt.plot(self.freq,
+                 np.abs(self.power),
+                 marker,
+                 color='b',
+                 label='Amplitude')
+        plt.plot(self.freq,
+                 np.abs(self.power.real),
+                 marker,
+                 color='r',
+                 alpha=0.5,
+                 label='Real Part')
+        plt.plot(self.freq,
+                 np.abs(self.power.imag),
+                 marker,
+                 color='g',
+                 alpha=0.5,
+                 label='Imaginary Part')
 
         if labels is not None:
             try:
@@ -990,6 +992,9 @@ class AveragedCrossspectrum(Crossspectrum):
         For backwards compatibility only. Like ``data2``, but no
         :class:`stingray.events.EventList` objects allowed
 
+    large_data : bool, default False
+        Use only for data larger than 10**7 data points!! Uses zarr and dask for computation.
+
     Attributes
     ----------
     freq: numpy.ndarray
@@ -1024,12 +1029,11 @@ class AveragedCrossspectrum(Crossspectrum):
         ``[[gti0_0, gti0_1], [gti1_0, gti1_1], ...]`` -- Good Time intervals.
         They are calculated by taking the common GTI between the
         two light curves
-
     """
 
     def __init__(self, data1=None, data2=None, segment_size=None, norm='none',
                  gti=None, power_type="real", silent=False, lc1=None, lc2=None,
-                 dt=None):
+                 dt=None, large_data=False):
 
         if lc1 is not None or lc2 is not None:
             warnings.warn("The lcN keywords are now deprecated. Use dataN "
@@ -1040,13 +1044,46 @@ class AveragedCrossspectrum(Crossspectrum):
         if data2 is None:
             data2 = lc2
 
-        self.type = "crossspectrum"
-
         if segment_size is None and data1 is not None:
             raise ValueError("segment_size must be specified")
         if segment_size is not None and not np.isfinite(segment_size):
             raise ValueError("segment_size must be finite!")
 
+        if large_data and data1 is not None and data2 is not None:
+            if isinstance(data1, EventList):
+                input_data = 'EventList'
+            elif isinstance(data1, Lightcurve):
+                input_data = 'Lightcurve'
+                chunks = int(np.rint(segment_size // data1.dt))
+                segment_size = chunks * data1.dt
+            else:
+                raise ValueError(
+                    f'Invalid input data type: {type(data1).__name__}')
+
+            dir_path1 = saveData(data1, persist=False, chunks=chunks)
+            dir_path2 = saveData(data2, persist=False, chunks=chunks)
+
+            data_path1 = genDataPath(dir_path1)
+            data_path2 = genDataPath(dir_path2)
+
+            spec = createChunkedSpectra(input_data,
+                                        'AveragedCrossspectrum',
+                                        data_path=list(data_path1 +
+                                                       data_path2),
+                                        segment_size=segment_size,
+                                        norm=norm,
+                                        gti=gti,
+                                        power_type=power_type,
+                                        silent=silent,
+                                        dt=dt)
+
+            for key, val in spec.__dict__.items():
+                setattr(self, key, val)
+
+            return
+
+        self.type = "crossspectrum"
+
         self.segment_size = segment_size
         self.power_type = power_type
 
@@ -1155,7 +1192,6 @@ class AveragedCrossspectrum(Crossspectrum):
         nphots1_all = []
         nphots2_all = []
 
-
         start_inds, end_inds = \
             bin_intervals_from_gtis(current_gtis, segment_size, lc1.time,
                                     dt=lc1.dt)
diff --git a/stingray/events.py b/stingray/events.py
index 16e8e3d9..d8ebf954 100644
--- a/stingray/events.py
+++ b/stingray/events.py
@@ -9,13 +9,13 @@ import copy
 import numpy as np
 import numpy.random as ra
 
-from .io import read, write
-from .utils import simon, assign_value_if_none
-from .filters import get_deadtime_mask
-from .gti import cross_gtis, append_gtis, check_separate
+from stingray.simulator.base import simulate_times
 
+from .filters import get_deadtime_mask
+from .gti import append_gtis, check_separate, cross_gtis
+from .io import read, write
 from .lightcurve import Lightcurve
-from stingray.simulator.base import simulate_times
+from .utils import assign_value_if_none, simon
 
 __all__ = ['EventList']
 
@@ -92,7 +92,7 @@ class EventList(object):
         self.notes = notes
         self.dt = dt
         self.mjdref = mjdref
-        self.gti = gti
+        self.gti = np.asarray(gti) if gti is not None else None
         self.pi = pi
         self.ncounts = ncounts
 
@@ -139,7 +139,6 @@ class EventList(object):
                                           gti=self.gti, tseg=tseg,
                                           mjdref=self.mjdref)
 
-
     def to_lc_list(self, dt):
         """Convert event list to a generator of Lightcurves.
 
@@ -258,8 +257,9 @@ class EventList(object):
 
         # Assign energies to events corresponding to the random numbers drawn
         self.energy = \
-            np.asarray([energy[
-                np.argwhere(cum_prob == np.min(cum_prob[(cum_prob - r) > 0]))]
+            np.asarray([
+                energy[np.argwhere(
+                    cum_prob == np.min(cum_prob[(cum_prob - r) > 0]))]
                       for r in R])
 
     def join(self, other):
diff --git a/stingray/io.py b/stingray/io.py
index 6c217c45..5e3774a7 100644
--- a/stingray/io.py
+++ b/stingray/io.py
@@ -1,35 +1,28 @@
-
-
-from collections.abc import Iterable
 import logging
 import math
-import numpy as np
 import os
-import six
+import pickle
 import warnings
+from collections.abc import Iterable
 
+import numpy as np
 from astropy.io import fits
 from astropy.table import Table
 from astropy.logger import AstropyUserWarning
 
 import stingray.utils as utils
-from .utils import order_list_of_arrays, is_string
-from .utils import assign_value_if_none
 
+from .utils import assign_value_if_none, is_string, order_list_of_arrays
 from .gti import get_gti_from_all_extensions, load_gtis
 
-try:
-    # Python 2
-    import cPickle as pickle
-except:
-    # Python 3
-    import pickle
+# Python 3
+import pickle
 
 _H5PY_INSTALLED = True
 
 try:
     import h5py
-except:
+except ImportError:
     _H5PY_INSTALLED = False
 
 
@@ -122,7 +115,7 @@ def high_precision_keyword_read(hdr, keyword):
     try:
         value = np.longdouble(hdr[keyword])
         return value
-    except:
+    except KeyError:
         pass
     try:
         if len(keyword) == 8:
@@ -130,7 +123,7 @@ def high_precision_keyword_read(hdr, keyword):
         value = np.longdouble(hdr[keyword + 'I'])
         value += np.longdouble(hdr[keyword + 'F'])
         return value
-    except:
+    except KeyError:
         return None
 
 
@@ -140,7 +133,7 @@ def _get_additional_data(lctable, additional_columns):
         for a in additional_columns:
             try:
                 additional_data[a] = np.array(lctable.field(a))
-            except:  # pragma: no cover
+            except KeyError:  # pragma: no cover
                 if a == 'PI':
                     logging.warning('Column PI not found. Trying with PHA')
                     additional_data[a] = np.array(lctable.field('PHA'))
@@ -256,8 +249,8 @@ def load_events_and_gtis(
     # Load data table
     try:
         lctable = lchdulist[hduname].data
-    except Exception:  # pragma: no cover
-        warnings.warn("HDU %s not found. Trying first extension" % hduname)
+    except KeyError:  # pragma: no cover
+        logging.warning('HDU %s not found. Trying first extension' % hduname)
         lctable = lchdulist[1].data
         hduname = 1
 
@@ -269,7 +262,7 @@ def load_events_and_gtis(
     # Read TIMEZERO keyword and apply it to events
     try:
         timezero = np.longdouble(header["TIMEZERO"])
-    except Exception:  # pragma: no cover
+    except KeyError:  # pragma: no cover
         timezero = np.longdouble(0.0)
 
     instr = mission = 'unknown'
@@ -284,7 +277,7 @@ def load_events_and_gtis(
     try:
         t_start = np.longdouble(header["TSTART"])
         t_stop = np.longdouble(header["TSTOP"])
-    except Exception:  # pragma: no cover
+    except KeyError:  # pragma: no cover
         warnings.warn(
             "Tstart and Tstop error. using defaults", AstropyUserWarning
         )
@@ -420,7 +413,7 @@ def read_header_key(fits_file, key, hdu=1):
     hdulist = fits.open(fits_file, ignore_missing_end=True)
     try:
         value = hdulist[hdu].header[key]
-    except:  # pragma: no cover
+    except KeyError:  # pragma: no cover
         value = ''
     hdulist.close()
     return value
@@ -727,7 +720,7 @@ def _retrieve_ascii_object(filename, **kwargs):
     data : astropy.Table object
         An astropy.Table object with the data from the file
     """
-    if not isinstance(filename, six.string_types):
+    if not isinstance(filename, str):
         raise TypeError("filename must be string!")
 
     if 'usecols' in list(kwargs.keys()):
@@ -878,7 +871,7 @@ def _retrieve_fits_object(filename, **kwargs):
         The name of file with which object was saved
 
     Other Parameters
-    -----------------------------
+    ----------------
     cols: str iterable
         The names of columns to extract from fits tables.
 
diff --git a/stingray/largememory.py b/stingray/largememory.py
new file mode 100644
index 00000000..51d92f91
--- /dev/null
+++ b/stingray/largememory.py
@@ -0,0 +1,1129 @@
+import os
+import random
+import string
+import warnings
+
+import numpy as np
+from astropy import log
+from astropy.io import fits
+
+import stingray
+
+from .events import EventList
+from .gti import cross_two_gtis
+from .io import high_precision_keyword_read
+from .lightcurve import Lightcurve
+from .utils import genDataPath
+
+HAS_ZARR = False
+try:
+    import zarr
+
+    HAS_ZARR = True
+    from numcodecs import Blosc
+except ImportError:
+    warnings.warn(
+        "Large Datasets may not be processed efficiently due to "
+        "computational constraints"
+    )
+
+__all__ = ["createChunkedSpectra", "saveData", "retrieveData"]
+
+
+def _saveChunkLC(lc, dir_name, chunks):
+    """
+    Save Lightcurve in chunks on disk.
+
+    Parameters
+    ----------
+    lc: :class:`stingray.Lightcurve` object
+        Lightcurve to be saved
+
+    dir_name: string
+        Top Level diretory name where Lightcurve is to be saved
+
+    chunks: int
+        The number of elements per chunk
+    """
+    # Creating a Nested Store and multiple groups for temporary saving
+    store = zarr.NestedDirectoryStore(dir_name)
+    lc_data_group = zarr.group(store=store, overwrite=True)
+    main_data_group = lc_data_group.create_group("main_data", overwrite=True)
+    meta_data_group = lc_data_group.create_group("meta_data", overwrite=True)
+
+    compressor = Blosc(cname="lz4", clevel=1, shuffle=-1)  # Optimal
+
+    main_data_group.create_dataset(
+        name="times",
+        data=lc.time,
+        compressor=compressor,
+        overwrite=True,
+        chunks=(chunks,),
+    )
+
+    main_data_group.create_dataset(
+        name="counts",
+        data=lc.counts,
+        compressor=compressor,
+        overwrite=True,
+        chunks=(chunks,),
+    )
+
+    if lc._counts_err is not None:
+        main_data_group.create_dataset(
+            name="count_err",
+            data=lc.counts_err,
+            compressor=compressor,
+            overwrite=True,
+            chunks=(chunks,),
+        )
+
+    main_data_group.create_dataset(
+        name="gti", data=lc.gti.flatten(), overwrite=True
+    )
+
+    meta_data_group.create_dataset(
+        name="dt", data=lc.dt, compressor=compressor, overwrite=True
+    )
+
+    meta_data_group.create_dataset(
+        name="err_dist",
+        data=lc.err_dist,
+        compressor=compressor,
+        overwrite=True,
+    )
+
+    meta_data_group.create_dataset(
+        name="mjdref", data=lc.mjdref, compressor=compressor, overwrite=True
+    )
+
+
+def _saveChunkEV(ev, dir_name, chunks):
+    """
+    Save EventList in chunks on disk.
+
+    Parameters
+    ----------
+    ev: :class:`stingray.events.EventList` object
+        EventList to be saved
+
+    dir_name: string
+        Top Level diretory name where EventList is to be saved
+
+    chunks: int
+        The number of elements per chunk
+
+    Raises
+    ------
+    ValueError
+        If there is no data being saved
+    """
+    # Creating a Nested Store and multiple groups for temporary saving
+    store = zarr.NestedDirectoryStore(dir_name)
+    ev_data_group = zarr.group(store=store, overwrite=True)
+    main_data_group = ev_data_group.create_group("main_data", overwrite=True)
+    meta_data_group = ev_data_group.create_group("meta_data", overwrite=True)
+
+    compressor = Blosc(cname="lz4", clevel=1, shuffle=-1)
+
+    if ev.time is not None and (ev.time.all() or ev.time.size != 0):
+        main_data_group.create_dataset(
+            name="times",
+            data=ev.time,
+            compressor=compressor,
+            overwrite=True,
+            chunks=(chunks,),
+        )
+
+    if ev.energy is not None and (ev.energy.all() or ev.energy.size != 0):
+        main_data_group.create_dataset(
+            name="energy",
+            data=ev.energy,
+            compressor=compressor,
+            overwrite=True,
+            chunks=(chunks,),
+        )
+
+    if ev.pi is not None and (ev.pi.all() or ev.pi.size != 0):
+        main_data_group.create_dataset(
+            name="pi_channel",
+            data=ev.pi,
+            compressor=compressor,
+            overwrite=True,
+            chunks=(chunks,),
+        )
+
+    if ev.gti is not None and (ev.gti.all() or ev.gti.shape[0] != 0):
+        main_data_group.create_dataset(
+            name="gti", data=ev.gti.flatten(), overwrite=True
+        )
+
+    if ev.dt != 0:
+        meta_data_group.create_dataset(
+            name="dt", data=ev.dt, compressor=compressor, overwrite=True
+        )
+
+    if ev.ncounts:
+        meta_data_group.create_dataset(
+            name="ncounts",
+            data=ev.ncounts,
+            compressor=compressor,
+            overwrite=True,
+        )
+
+    if ev.notes:
+        meta_data_group.create_dataset(
+            name="notes", data=ev.notes, compressor=compressor, overwrite=True
+        )
+
+    meta_data_group.create_dataset(
+        name="mjdref", data=ev.mjdref, compressor=compressor, overwrite=True
+    )
+
+
+def _saveFITSZarr(f_name, dir_name, chunks):
+    """
+    Read a FITS file and save it for further processing.
+
+    Parameters
+    ----------
+    f_name: string
+        The name of file with which object was saved
+
+    dir_name: string
+        The name of the top level directory where the file is to be stored
+
+    chunks: int
+        The number of elements per chunk
+    """
+
+    compressor = Blosc(cname="lz4", clevel=1, shuffle=-1)
+
+    store = zarr.NestedDirectoryStore(dir_name)
+    fits_data_group = zarr.group(store=store, overwrite=True)
+    main_data_group = fits_data_group.create_group("main_data", overwrite=True)
+    meta_data_group = fits_data_group.create_group("meta_data", overwrite=True)
+
+    with fits.open(f_name, memmap=True) as fits_data:
+        for HDUList in fits_data:
+            if HDUList.name == "EVENTS":
+                times = HDUList.data["TIME"]
+                chunks = times.size if times.size < chunks else chunks
+
+                main_data_group.create_dataset(
+                    name="times",
+                    data=times,
+                    compressor=compressor,
+                    overwrite=True,
+                    chunks=(chunks,),
+                )
+
+                for col in ["PI", "PHA"]:
+                    if col in HDUList.data.columns.names:
+                        main_data_group.create_dataset(
+                            name=f"{col.lower()}_channel",
+                            data=HDUList.data[col],
+                            compressor=compressor,
+                            overwrite=True,
+                            chunks=(chunks,),
+                        )
+
+                meta_data_group.create_dataset(
+                    name="tstart",
+                    data=HDUList.header["TSTART"],
+                    compressor=compressor,
+                    overwrite=True,
+                )
+
+                meta_data_group.create_dataset(
+                    name="tstop",
+                    data=HDUList.header["TSTOP"],
+                    compressor=compressor,
+                    overwrite=True,
+                )
+
+                meta_data_group.create_dataset(
+                    name="mjdref",
+                    data=high_precision_keyword_read(HDUList.header, "MJDREF"),
+                    compressor=compressor,
+                    overwrite=True,
+                )
+
+            elif HDUList.name == "GTI":
+                # TODO: Needs to be generalized
+                start, stop = HDUList.data["START"], HDUList.data["STOP"]
+                gti = np.array(list(zip(start, stop)))
+                main_data_group.create_dataset(
+                    name="gti",
+                    data=gti.flatten(),
+                    compressor=compressor,
+                    overwrite=True,
+                )
+
+
+def saveData(data, persist=False, dir_name=None, chunks=None):
+    """
+    Saves Lightcurve/EventList or any such data in chunks to disk.
+
+    Parameters
+    ----------
+    data: :class:`stingray.Lightcurve` or :class:`stingray.events.EventList` object or string
+        Data to be stored on the disk.
+
+    persist: bool
+        If the data is to be stored on the disk permanently.
+
+    dir_name: string, optional
+        Name of top level directory where data is to be stored, by default randomNameGenerate()
+
+    chunks: int, optional
+        Length of data chunks in number of bins. If None, it is calculated
+        based on the system resources
+
+    Returns
+    -------
+    string
+        Name of top level directory where data is to be stored
+
+    Raises
+    ------
+    ValueError
+        If data is not a Lightcurve or EventList
+    """
+    from sys import platform
+
+    HAS_PSUTIL = False
+    try:
+        import psutil
+
+        HAS_PSUTIL = True
+    except ImportError:
+        if not (platform == "linux" or platform == "linux2"):
+            warnings.warn(
+                "The chunk size will not depend on available RAM and will slowdown execution."
+            )
+
+    if not persist:
+        import tempfile
+
+        dir_name = tempfile.mkdtemp()
+    else:
+        dir_name = (
+            "".join(random.choices(string.ascii_letters + string.digits, k=10))
+            if dir_name is None
+            else dir_name
+        )
+
+    if chunks is None:
+        ideal_chunk, safe_chunk = 8388608, 4193404
+
+        if HAS_PSUTIL:
+            free_m = (
+                psutil.virtual_memory().available + psutil.swap_memory().free
+            ) / 10 ** 9
+            chunks = ideal_chunk if free_m >= 10.0 else safe_chunk
+        else:
+            if platform == "linux" or platform == "linux2":
+                free_m = int(
+                    os.popen("free -t -m").readlines()[-1].split()[-1]
+                )
+                chunks = ideal_chunk if free_m >= 10000 else safe_chunk
+            else:
+                chunks = safe_chunk
+
+    if isinstance(data, Lightcurve):
+        if data.time.size < chunks:
+            chunks = data.time.size
+        _saveChunkLC(data, dir_name, chunks)
+
+    elif isinstance(data, EventList):
+        if not (
+            data.time is not None and (data.time.all() or data.time.size != 0)
+        ):
+            raise ValueError(
+                ("The EventList passed is empty and hence cannot be saved")
+            )
+
+        if data.time.size > 0 and data.time.size < chunks:
+            chunks = data.time.size
+
+        _saveChunkEV(data, dir_name, chunks)
+
+    elif os.path.exists(data) and os.stat(data).st_size > 0:
+        _saveFITSZarr(data, dir_name, chunks)
+
+    else:
+        raise ValueError((f"Invalid data: {data} ({type(data).__name__})"))
+
+    return dir_name
+
+
+def _retrieveDataLC(data_path, chunk_size=0, offset=0, raw=False):
+    """
+    Retrieve data from stored Lightcurve on disk.
+
+    Parameters
+    ----------
+    data_path: list
+        Path to datastore
+
+    chunk_size: int, optional
+        Size of data to be retrieved, default 0
+
+    offset: int, optional
+        Offset or start element to read the array from, default 0
+
+    raw: bool, optional
+        Only to be used for if raw memory mapped zarr arrays are to be obtained, default False
+
+    Returns
+    -------
+    :class:`stingray.Lightcurve` object or tuple
+        Lightcurve retrieved from store or data of Lightcurve
+
+    Raises
+    ------
+    ValueError
+        If offset provided is larger than size of array.
+    """
+    times = zarr.open_array(store=data_path[0], mode="r", path="times")
+    counts = zarr.open_array(store=data_path[0], mode="r", path="counts")
+    try:
+        count_err = zarr.open_array(
+            store=data_path[0], mode="r", path="count_err"
+        )
+    except ValueError:
+        count_err = None
+
+    gti = zarr.open_array(store=data_path[0], mode="r", path="gti")[...]
+    gti = gti.reshape((gti.size // 2, 2))
+
+    dt = zarr.open_array(store=data_path[1], mode="r", path="dt")[...]
+    mjdref = zarr.open_array(store=data_path[1], mode="r", path="mjdref")[...]
+    err_dist = zarr.open_array(store=data_path[1], mode="r", path="err_dist")[
+        ...
+    ]
+
+    if raw:
+        return (
+            times,
+            counts,
+            count_err,
+            gti,
+            float(dt),
+            str(err_dist),
+            float(mjdref),
+        )
+    else:
+        if chunk_size == 0 or chunk_size > times.size:
+            chunk_size = times.size
+            log.info(
+                f"Retrieving the whole array (Chunk size > size of times)"
+            )
+
+        if offset > times.size:
+            raise ValueError(
+                (f"Offset cannot be larger than size of array {times.size}")
+            )
+
+        times = times.get_basic_selection(slice(offset, chunk_size))[...]
+        gti_new = cross_two_gtis(
+            gti,
+            np.asarray(
+                [[times[0] - 0.5 * float(dt), times[-1] + 0.5 * float(dt)]]
+            ),
+        )
+
+        return Lightcurve(
+            time=times,
+            counts=counts.get_basic_selection(slice(offset, chunk_size))[...],
+            err=count_err.get_basic_selection(slice(offset, chunk_size))[...]
+            if count_err is not None
+            else None,
+            gti=gti_new,
+            dt=float(dt),
+            err_dist=str(err_dist),
+            mjdref=float(mjdref),
+            skip_checks=True,
+        )
+
+
+def _retrieveDataEV(data_path, chunk_size=0, offset=0, raw=False):
+    """
+    Retrieve data from stored EventList on disk.
+
+    Parameters
+    ----------
+    data_path: list
+        Path to datastore.
+
+    chunk_size: int
+        Size of data to be retrieved
+
+    offset: int
+        Offset or start element to read the array from
+
+    raw: bool
+        Only to be used for if raw memory mapped zarr arrays are to be obtained
+
+    Returns
+    -------
+    :class:`stingray.events.EventList` object or tuple
+        EventList or data of EventList retrieved from store.
+
+    Raises
+    ------
+    ValueError
+        If array does not exist at path
+
+    ValueError
+        If offset provided is larger than size of array.
+
+    ValueError
+        If the file to read is empty
+    """
+    try:
+        times = zarr.open_array(store=data_path[0], mode="r", path="times")
+    except ValueError:
+        raise ValueError("This file contains no events.")
+
+    try:
+        energy = zarr.open_array(store=data_path[0], mode="r", path="energy")
+    except ValueError:
+        energy = None
+
+    try:
+        pi_channel = zarr.open_array(
+            store=data_path[0], mode="r", path="pi_channel"
+        )
+    except ValueError:
+        pi_channel = None
+
+    try:
+        gti = zarr.open_array(store=data_path[0], mode="r", path="gti")[...]
+        gti = gti.reshape((gti.size // 2, 2))
+    except ValueError:
+        gti = None
+
+    try:
+        dt = zarr.open_array(store=data_path[1], mode="r", path="dt")[...]
+    except ValueError:
+        dt = 0
+
+    try:
+        ncounts = zarr.open_array(
+            store=data_path[1], mode="r", path="ncounts"
+        )[...]
+    except ValueError:
+        ncounts = None
+
+    try:
+        mjdref = zarr.open_array(store=data_path[1], mode="r", path="mjdref")[
+            ...
+        ]
+    except ValueError:
+        mjdref = 0
+
+    try:
+        notes = zarr.open_array(store=data_path[1], mode="r", path="notes")[
+            ...
+        ]
+    except ValueError:
+        notes = ""
+
+    if raw:
+        return (
+            times,
+            energy,
+            int(ncounts),
+            float(mjdref),
+            dt,
+            gti,
+            pi_channel,
+            str(notes),
+        )
+    else:
+        if chunk_size == 0 or chunk_size > times.size:
+            chunk_size = times.size
+            log.info(
+                f"Retrieving the whole array (Chunk size > size of times)"
+            )
+
+        if offset > times.size:
+            raise ValueError(
+                "No element read. Offset cannot be larger than size of array"
+            )
+
+        time = times.get_basic_selection(slice(offset, chunk_size))
+        energy = (
+            energy.get_basic_selection(slice(offset, chunk_size))
+            if energy is not None
+            else None
+        )
+        pi = (
+            pi_channel.get_basic_selection(slice(offset, chunk_size))
+            if pi_channel is not None
+            else None
+        )
+
+        gti_new = None
+        if gti is not None:
+            gti_new = cross_two_gtis(
+                gti,
+                np.asarray(
+                    [[times[0] - 0.5 * float(dt), times[-1] + 0.5 * float(dt)]]
+                ),
+            )
+            warnings.warn(
+                "GTI management is tricky when partially loading event lists "
+                "from ZARR arrays. Please use with care.",
+                UserWarning,
+            )
+
+        if pi_channel is not None:
+            pi_channel = pi_channel.get_basic_selection(
+                slice(offset, chunk_size)
+            )
+        else:
+            pi_channel = None
+
+        return EventList(
+            time=time,
+            energy=energy,
+            ncounts=ncounts,
+            mjdref=float(mjdref),
+            dt=float(dt),
+            gti=gti_new,
+            pi=pi_channel,
+            notes=str(notes),
+        )
+
+
+def retrieveData(
+    data_type, dir_path, chunk_data=False, chunk_size=0, offset=0, raw=False
+):
+    """
+    Retrieves Lightcurve/EventList or any such data from disk.
+
+    Parameters
+    ----------
+    data_type: string
+        Type of data to retrieve i.e. Lightcurve, Eventlist, FITS data to retrieve
+
+    dir_path: string
+        Path to zarr datastore + Top level directory name for data
+
+    chunk_data: bool, optional
+        If only a chunk of data is to be retrieved, by default False
+
+    chunk_size: int, optional
+        Number of values to be retrieved, by default 0
+
+    offset: int, optional
+        Start offset from where values are to be retrieved, by default 0
+
+    raw: bool, optional
+        Only to be used for if raw memory mapped zarr arrays are to be obtained, by default False
+
+    Returns
+    -------
+    :class:`stingray.events.EventList` object or :class:`stingray.Lighrcurve` object or tuple
+        EventList or Lightcurve created from store or raw data
+
+    Raises
+    ------
+    ValueError
+        If datatype is not Lightcurve or EventList of FITS
+    """
+    data_path = genDataPath(dir_path)
+
+    if chunk_data is True and offset >= chunk_size:
+        raise ValueError(("offset should be less than chunk_size"))
+
+    if data_type.lower() == "lightcurve":
+        if chunk_data is True and chunk_size > 0:
+            return _retrieveDataLC(
+                data_path, int(chunk_size), int(offset), raw=raw
+            )
+        else:
+            return _retrieveDataLC(data_path, raw=raw)
+
+    elif data_type.lower() == "eventlist" or data_type.lower() == "fits":
+        if chunk_data is True and chunk_size > 0:
+            return _retrieveDataEV(
+                data_path, int(chunk_size), int(offset), raw=False
+            )
+        else:
+            return _retrieveDataEV(data_path, raw=raw)
+
+    else:
+        raise ValueError(
+            (f"Invalid data: {data_type} ({type(data_type).__name__})")
+        )
+
+
+def _combineSpectra(final_spectra):
+    """
+    Create a final spectra that is the mean of all spectra.
+
+    Parameters
+    ----------
+    final_spectra: :class:`stingray.AveragedCrossspectrum/AveragedPowerspectrum' object
+        Summed spectra of all spectra
+
+    Returns
+    -------
+    :class:`stingray.events.EventList` object or :class:`stingray.Lighrcurve` object
+        Final resulting spectra.
+    """
+    final_spectra.power /= final_spectra.m
+    final_spectra.unnorm_power /= final_spectra.m
+    final_spectra.power_err = (
+        np.sqrt(final_spectra.power_err) / final_spectra.m
+    )
+
+    if isinstance(
+        final_spectra, stingray.AveragedCrossspectrum
+    ) and not isinstance(final_spectra, stingray.AveragedPowerspectrum):
+        final_spectra.pds1.power /= final_spectra.m
+        final_spectra.pds2.power /= final_spectra.m
+        final_spectra.nphots1 /= final_spectra.m
+        final_spectra.nphots2 /= final_spectra.m
+    else:
+        final_spectra.nphots /= final_spectra.m
+
+    return final_spectra
+
+
+def _addSpectra(final_spectra, curr_spec):
+    """
+    Add various Spectra(AveragedCrossspectrum/AveragedPowerspectrum) for combination.
+
+    Parameters
+    ----------
+    final_spectra: object
+        Final Combined AveragedCrossspectrum or AveragedPowerspectrum
+
+    curr_spec: object
+        AveragedCrossspectrum/AveragedPowerspectrum to be combined
+
+    first_iter: bool
+        Check for first iteration variable
+
+    Returns
+    -------
+    final_spectrum: :class:`stingray.crossspectrum.AveragedCrossspectrum` object or :class:`stingray.powerspectrum.AveragedPowerspectrum` object
+        Combined AveragedCrossspectrum/AveragedPowerspectrum
+    """
+
+    final_spectra.power += curr_spec.power * curr_spec.m
+    final_spectra.unnorm_power += curr_spec.unnorm_power * curr_spec.m
+    final_spectra.power_err += curr_spec.power_err ** 2 * curr_spec.m
+
+    final_spectra.m += curr_spec.m
+
+    if isinstance(final_spectra, stingray.AveragedPowerspectrum):
+        final_spectra.nphots += curr_spec.nphots * curr_spec.m
+
+    elif isinstance(final_spectra, stingray.AveragedCrossspectrum):
+        final_spectra.pds1.power += curr_spec.pds1.power * curr_spec.m
+        final_spectra.pds1.power += curr_spec.pds1.power * curr_spec.m
+
+        final_spectra.nphots1 += curr_spec.nphots1 * curr_spec.m
+        final_spectra.nphots2 += curr_spec.nphots2 * curr_spec.m
+
+    return final_spectra
+
+
+def _chunkLCSpec(
+    data_path, spec_type, segment_size, norm, gti, power_type, silent
+):
+    """
+    Create a chunked spectra from Lightcurve stored on disk.
+
+    Parameters
+    ----------
+    data_path: string
+        Path to stored Lightcurve or EventList chunks on disk
+
+    spec_type: string
+        Type of spectra to create AveragedCrossspectrum or AveragedPowerspectrum.
+
+    segment_size: float
+        The size of each segment to average in the AveragedCrossspectrum/AveragedPowerspectrum.
+
+    norm: {``frac``, ``abs``, ``leahy``, ``none``}
+        The normalization of the (real part of the) cross spectrum.
+
+    gti: 2-d float array
+        `[[gti0_0, gti0_1], [gti1_0, gti1_1], ...]`` -- Good Time intervals.
+        This choice overrides the GTIs in the single light curves. Use with
+        care!
+
+    power_type: string
+        Parameter to choose among complete, real part and magnitude of
+         the cross spectrum. None for AveragedPowerspectrum
+
+    silent: bool
+        Do not show a progress bar when generating an averaged cross spectrum.
+        Useful for the batch execution of many spectra
+
+    dt1: float
+        The time resolution of the light curve. Only needed when constructing
+        light curves in the case where data1 or data2 are of :class:EventList
+
+    Returns
+    -------
+    :class:`stingray.events.EventList` object or :class:`stingray.Lighrcurve` object
+        Summed computed spectra
+
+    Raises
+    ------
+    ValueError
+        If spectra is not AveragedCrossspectrum or AveragedPowerspectrum
+
+    ValueError
+        If previous and current spectra frequencies are not identical
+    """
+    times, counts, counts_err, gti_new, dt, err_dist, mjdref = _retrieveDataLC(
+        data_path[0:2], raw=True
+    )
+    final_gti = gti_new.reshape((gti_new.size // 2, 2))
+
+    if spec_type.lower() == "averagedpowerspectrum":
+        fin_spec = stingray.AveragedPowerspectrum()
+
+    elif spec_type.lower() == "averagedcrossspectrum":
+        fin_spec = stingray.AveragedCrossspectrum()
+        (
+            times_other,
+            counts_other,
+            counts_err_other,
+            gti_other,
+            dt_other,
+            err_dist_other,
+            mjdref_other,
+        ) = _retrieveDataLC(data_path[2:4], raw=True)
+        gti_new_other = gti_other.reshape((gti_other.size // 2, 2))
+        final_gti = cross_two_gtis(final_gti, gti_new_other)
+    else:
+        raise ValueError((f"Invalid spectra-type {spec_type}"))
+
+    nchunks_per_segment = max(
+        int(np.rint(segment_size / dt / times.chunks[0])), 1
+    )
+
+    new_segment_size = nchunks_per_segment * times.chunks[0] * dt
+    if not np.isclose(new_segment_size, segment_size):
+        warnings.warn("Adapting segment_size to a fixed number of chunks")
+        segment_size = new_segment_size
+
+    step = times.chunks[0] * nchunks_per_segment
+    prev_freq = None
+    for i in range(0, times.size, step):
+        nmax = i + step
+        if nmax > times.size and (nmax - times.size) * dt < segment_size:
+            continue
+        time = times.get_basic_selection(slice(i, i + step))
+        count = counts.get_basic_selection(slice(i, i + step))
+        count_err = None
+        if counts_err is not None:
+            count_err = counts_err.get_basic_selection(slice(i, i + step))
+
+        gti_new = np.asarray([[time[0] - 0.5 * dt, time[-1] + 0.5 * dt]])
+
+        lc1 = Lightcurve(
+            time=time,
+            counts=count,
+            err=count_err,
+            gti=gti_new,
+            err_dist=err_dist,
+            mjdref=mjdref,
+            dt=dt,
+            skip_checks=True,
+        )
+
+        if isinstance(fin_spec, stingray.AveragedPowerspectrum):
+            avg_spec = stingray.AveragedPowerspectrum(
+                data=lc1,
+                segment_size=segment_size,
+                norm=norm,
+                gti=gti_new,
+                silent=True,
+                large_data=False,
+            )
+
+        else:
+            time_other = times_other.get_basic_selection(slice(i, i + step))
+            count_other = counts_other.get_basic_selection(slice(i, i + step))
+            count_err_other = None
+            if counts_err_other is not None:
+                count_err_other = counts_err_other.get_basic_selection(
+                    slice(i, i + step)
+                )
+
+            gti_new_other = np.asarray(
+                [[time_other[0] - 0.5 * dt, time_other[-1] + 0.5 * dt]]
+            )
+
+            lc2 = Lightcurve(
+                time=time_other,
+                counts=count_other,
+                err=count_err_other,
+                gti=gti_new_other,
+                err_dist=err_dist_other,
+                mjdref=mjdref_other,
+                dt=dt_other,
+                skip_checks=True,
+            )
+
+            assert np.allclose(
+                gti_new, gti_new_other
+            ), "The two light curves should cover identical time frames"
+
+            avg_spec = stingray.AveragedCrossspectrum(
+                data1=lc1,
+                data2=lc2,
+                segment_size=segment_size,
+                norm=norm,
+                gti=gti_new,
+                power_type=power_type,
+                silent=True,
+                large_data=False,
+            )
+
+        # First iteration
+        if i == 0:
+            fin_spec = avg_spec
+            prev_freq = avg_spec.freq
+            continue
+
+        # All others
+        if not np.allclose(prev_freq, avg_spec.freq):
+            bad = ~np.isclose(prev_freq, avg_spec.freq)
+            raise ValueError(
+                (
+                    f"Spectra have unequal frequencies: "
+                    f"{avg_spec.freq[bad]} vs {prev_freq[bad]}"
+                )
+            )
+        fin_spec = _addSpectra(fin_spec, avg_spec)
+        prev_freq = avg_spec.freq
+
+    fin_spec.gti = final_gti
+    return fin_spec
+
+
+# TODO: Add support for eventlist soon!
+# def _chunkEVSpec(data_path, spec_type, segment_size, norm, gti, power_type,
+#                  silent, dt1):
+#     """
+#     Create a chunked spectra from EventList stored on disk.
+#
+#     Parameters
+#     ----------
+#     data_path: string
+#         Path to stored Lightcurve or EventList chunks on disk
+#
+#     spec_type: string
+#         Type of spectra to create AveragedCrossspectrum or AveragedPowerspectrum.
+#
+#     segment_size: float
+#         The size of each segment to average in the AveragedCrossspectrum/AveragedPowerspectrum.
+#
+#     norm: {``frac``, ``abs``, ``leahy``, ``none``}
+#         The normalization of the (real part of the) cross spectrum
+#
+#     gti: 2-d float array
+#         `[[gti0_0, gti0_1], [gti1_0, gti1_1], ...]`` -- Good Time intervals.
+#         This choice overrides the GTIs in the single light curves. Use with
+#         care!
+#
+#     power_type: string
+#         Parameter to choose among complete, real part and magnitude of
+#          the cross spectrum. None for AveragedPowerspectrum
+#
+#     silent: bool
+#         Do not show a progress bar when generating an averaged cross spectrum.
+#         Useful for the batch execution of many spectra
+#
+#     dt1: float
+#         The time resolution of the light curve. Only needed when constructing
+#         light curves in the case where data1 or data2 are of :class:EventList
+#
+#     Returns
+#     -------
+#     :class:`stingray.events.EventList` object or :class:`stingray.Lighrcurve` object
+#         Summed computed spectra
+#
+#     Raises
+#     ------
+#     ValueError
+#         If spectra is not AveragedCrossspectrum or AveragedPowerspectrum
+#
+#     ValueError
+#         If previous and current spectra frequencies are not identical
+#     """
+#     times, energy, ncounts, mjdref, dt, gti_new, pi_channel, notes = _retrieveDataEV(
+#         data_path[0:2], raw=True)
+#     gti_new = gti_new.reshape((gti_new.size // 2, 2))
+#
+#     if spec_type.lower() == 'averagedpowerspectrum':
+#         fin_spec = stingray.AveragedPowerspectrum()
+#
+#     elif spec_type.lower() == 'averagedcrossspectrum':
+#         fin_spec = stingray.AveragedCrossspectrum()
+#         times_other, energy_other, ncounts_other, mjdref_other, dt_other, gti_other, pi_channel_other = _retrieveDataEV(
+#             data_path[2:4], raw=True)
+#         gti_new_other = gti_other.reshape((gti_other.size // 2, 2))
+#
+#     else:
+#         raise ValueError((f"Invalid spectra {spec_type}"))
+#
+#     first_iter = True
+#     for i in range(times.chunks[0], times.size, times.chunks[0]):
+#         gti_new = cross_two_gtis(
+#             gti_new,
+#             np.asarray([[
+#                 times.get_basic_selection(i - times.chunks[0]) - 0.5 * dt,
+#                 times.get_basic_selection(i) + 0.5 * dt
+#             ]]))
+#         ev1 = EventList(
+#             time=times.get_basic_selection(slice(i - times.chunks[0], i))
+#             if times is not None else None,
+#             energy=energy.get_basic_selection(slice(i - times.chunks[0], i))
+#             if energy is not None else None,
+#             ncounts=ncounts,
+#             mjdref=mjdref,
+#             dt=dt,
+#             gti=gti_new,
+#             pi=pi_channel.get_basic_selection(slice(i - times.chunks[0], i))
+#             if isinstance(pi_channel, zarr.core.Array) else pi_channel,
+#             notes=notes)
+#
+#         if isinstance(fin_spec, stingray.AveragedPowerspectrum):
+#             avg_spec = stingray.AveragedPowerspectrum(
+#                 data=ev1,
+#                 segment_size=segment_size,
+#                 norm=norm,
+#                 gti=gti,
+#                 silent=silent,
+#                 dt=dt1,
+#                 large_data=False)
+#
+#         elif isinstance(fin_spec, stingray.AveragedCrossspectrum):
+#             gti_new_other = cross_two_gtis(
+#                 gti_other,
+#                 np.asarray([[
+#                     times_other.get_basic_selection(i - times_other.chunks[0])
+#                     - 0.5 * dt_other,
+#                     times_other.get_basic_selection(i) + 0.5 * dt_other
+#                 ]]))
+#
+#             ev2 = EventList(
+#                 time=times_other.get_basic_selection(slice(i - times.chunks[0], i)) if time_other is not None else None,
+#                 energy=energy_other.get_basic_selection(slice(i - times.chunks[0], i)) if energy_other is not None else None,
+#                 ncounts=ncounts_other,
+#                 mjdref=mjdref_other,
+#                 dt=dt_other,
+#                 gti=gti_new_other,
+#                 pi=pi_channel_other.get_basic_selection(slice(i - times.chunks[0], i)) if isinstance(pi_channel_other, zarr.core.Array) else pi_channel_other,
+#                 notes=notes_other)
+#
+#             avg_spec = stingray.AveragedCrossspectrum(
+#                 data1=ev1,
+#                 data2=ev2,
+#                 segment_size=segment_size,
+#                 norm=norm,
+#                 gti=gti,
+#                 power_type=power_type,
+#                 silent=silent,
+#                 dt=dt1,
+#                 large_data=False)
+#
+#         if first_iter:
+#             prev_freq = avg_spec.freq
+#             fin_spec = _addSpectra(fin_spec, avg_spec, first_iter)
+#         else:
+#             if np.array_equal(prev_freq, avg_spec.freq):
+#                 fin_spec = _addSpectra(fin_spec, avg_spec, first_iter)
+#                 prev_freq = avg_spec.freq
+#             else:
+#                 raise ValueError((
+#                     f"Spectra have unequal frequencies {avg_spec.freq.shape}{prev_freq.shape}"
+#                 ))
+#
+#         first_iter = False
+#     fin_spec.gti = cross_two_gtis(gti_new, gti_new_other)
+#
+#     return fin_spec
+
+
+def createChunkedSpectra(
+    data_type,
+    spec_type,
+    data_path,
+    segment_size,
+    norm,
+    gti,
+    power_type,
+    silent,
+    dt=None,
+):
+    """
+    Create a chunked spectra from zarr files stored on disk.
+
+    Parameters
+    ----------
+    data_type: string
+        Data in Lightcurve or EventList
+
+    spec_type: string
+        Type of spectra to create AveragedCrossspectrum or AveragedPowerspectrum
+
+    data_path: list
+        Path to datastore
+
+    segment_size: float
+        The size of each segment to average in the AveragedCrossspectrum/AveragedPowerspectrum
+
+    norm: {``frac``, ``abs``, ``leahy``, ``none``}
+        The normalization of the (real part of the) cross spectrum
+
+    gti: 2-d float array
+        `[[gti0_0, gti0_1], [gti1_0, gti1_1], ...]`` -- Good Time intervals.
+        This choice overrides the GTIs in the single light curves. Use with
+        care!
+
+    power_type: string
+        Parameter to choose among complete, real part and magnitude of
+         the cross spectrum. None for AveragedPowerspectrum
+
+    silent: bool
+        Do not show a progress bar when generating an averaged cross spectrum.
+        Useful for the batch execution of many spectra
+
+    dt: float, optional
+        The time resolution of the light curve. Only needed when constructing
+        light curves in the case where data1 or data2 are of :class:EventList, by default None
+
+    Returns
+    -------
+    :class:`stingray.events.EventList` object or :class:`stingray.Lighrcurve` object
+        Final computed spectra
+    """
+    if data_type.lower() == "lightcurve":
+        fin_spec = _chunkLCSpec(
+            data_path=data_path,
+            spec_type=spec_type,
+            segment_size=segment_size,
+            norm=norm,
+            gti=gti,
+            power_type=power_type,
+            silent=silent,
+        )
+
+    else:
+        raise RuntimeError("Only input light curves are allowed")
+
+    return _combineSpectra(fin_spec)
diff --git a/stingray/lightcurve.py b/stingray/lightcurve.py
index 467be5fa..196e7e2b 100644
--- a/stingray/lightcurve.py
+++ b/stingray/lightcurve.py
@@ -4,16 +4,18 @@ Definition of :class::class:`Lightcurve`.
 :class::class:`Lightcurve` is used to create light curves out of photon counting data
 or to save existing light curves in a class that's easy to use.
 """
-import warnings
 import logging
+import warnings
+
 import numpy as np
+
 import stingray.io as io
 import stingray.utils as utils
 from stingray.exceptions import StingrayError
-from stingray.utils import simon, assign_value_if_none, baseline_als
-from stingray.utils import poisson_symmetrical_errors
-from stingray.gti import cross_two_gtis, join_gtis, gti_border_bins
-from stingray.gti import check_gtis, create_gti_mask, bin_intervals_from_gtis
+from stingray.gti import (bin_intervals_from_gtis, check_gtis, create_gti_mask,
+                          cross_two_gtis, gti_border_bins, join_gtis)
+from stingray.utils import (assign_value_if_none, baseline_als,
+                            poisson_symmetrical_errors, simon)
 
 __all__ = ["Lightcurve"]
 
@@ -1163,11 +1165,11 @@ class Lightcurve(object):
 
         # tolerance for the newly created GTIs: Note that this seems to work
         # with a tolerance of 2, but not if I substitute 10. I don't know why
-        epsilon = np.min(tdiff)/2.0
+        epsilon = np.min(tdiff) / 2.0
 
         # calculate new GTIs
-        gti_start = np.hstack([self.time[0]-epsilon, self.time[gap_idx+1]-epsilon])
-        gti_stop = np.hstack([self.time[gap_idx]+epsilon, self.time[-1]+epsilon])
+        gti_start = np.hstack([self.time[0] - epsilon, self.time[gap_idx + 1] - epsilon])
+        gti_stop = np.hstack([self.time[gap_idx] + epsilon, self.time[-1] + epsilon])
 
         gti = np.vstack([gti_start, gti_stop]).T
         if hasattr(self, 'gti') and self.gti is not None:
@@ -1417,7 +1419,8 @@ class Lightcurve(object):
             infinite or nan points. Use at your own risk.
         """
         return Lightcurve(time=lk.time, counts=lk.flux,
-                          err=lk.flux_err, input_counts=False, skip_checks=skip_checks)
+                          err=lk.flux_err, input_counts=False,
+                          skip_checks=skip_checks)
 
     def plot(self, witherrors=False, labels=None, axis=None, title=None,
              marker='-', save=False, filename=None):
@@ -1617,7 +1620,7 @@ class Lightcurve(object):
             start = start_bins[i]
             stop = stop_bins[i]
 
-            if np.isclose(stop-start, 1):
+            if np.isclose(stop - start, 1):
                 logging.warning("Segment with a single time bin! Ignoring this segment!")
                 continue
             if (stop - start) < min_points:
diff --git a/stingray/powerspectrum.py b/stingray/powerspectrum.py
index 359596cd..56f3a36f 100755
--- a/stingray/powerspectrum.py
+++ b/stingray/powerspectrum.py
@@ -2,17 +2,19 @@ import warnings
 
 import numpy as np
 import scipy
-import scipy.fftpack
 import scipy.optimize
 import scipy.stats
 
-import stingray.lightcurve as lightcurve
 import stingray.utils as utils
-from stingray.crossspectrum import Crossspectrum, AveragedCrossspectrum
+from stingray.crossspectrum import AveragedCrossspectrum, Crossspectrum
 from stingray.gti import bin_intervals_from_gtis, check_gtis
+from stingray.largememory import createChunkedSpectra, saveData
 from stingray.stats import pds_probability
+from stingray.utils import genDataPath
+
 from .events import EventList
 from .gti import cross_two_gtis
+from .lightcurve import Lightcurve
 
 try:
     from tqdm import tqdm as show_progress
@@ -302,7 +304,6 @@ class AveragedPowerspectrum(AveragedCrossspectrum, Powerspectrum):
     norm: {``leahy`` | ``frac`` | ``abs`` | ``none`` }, optional, default ``frac``
         The normaliation of the periodogram to be used.
 
-
     Other Parameters
     ----------------
     gti: 2-d float array
@@ -318,6 +319,8 @@ class AveragedPowerspectrum(AveragedCrossspectrum, Powerspectrum):
         The time resolution of the light curve. Only needed when constructing
         light curves in the case where data is of :class:EventList
 
+    large_data : bool, default False
+        Use only for data larger than 10**7 data points!! Uses zarr and dask for computation.
 
     Attributes
     ----------
@@ -352,9 +355,8 @@ class AveragedPowerspectrum(AveragedCrossspectrum, Powerspectrum):
 
     """
     def __init__(self, data=None, segment_size=None, norm="frac", gti=None,
-                 silent=False, dt=None, lc=None):
+                 silent=False, dt=None, lc=None, large_data=False):
 
-        self.type = "powerspectrum"
         if lc is not None:
             warnings.warn("The lc keyword is now deprecated. Use data "
                           "instead", DeprecationWarning)
@@ -367,6 +369,37 @@ class AveragedPowerspectrum(AveragedCrossspectrum, Powerspectrum):
         if segment_size is not None and not np.isfinite(segment_size):
             raise ValueError("segment_size must be finite!")
 
+        if large_data and data is not None:
+            chunks = None
+
+            if isinstance(data, EventList):
+                input_data = 'EventList'
+            elif isinstance(data, Lightcurve):
+                input_data = 'Lightcurve'
+                chunks = int(np.rint(segment_size // data.dt))
+                segment_size = chunks * data.dt
+            else:
+                raise ValueError(
+                    f'Invalid input data type: {type(data).__name__}')
+
+            dir_path = saveData(data, persist=False, chunks=chunks)
+
+            data_path = genDataPath(dir_path)
+            spec = createChunkedSpectra(input_data,
+                                        'AveragedPowerspectrum',
+                                        data_path=data_path,
+                                        segment_size=segment_size,
+                                        norm=norm,
+                                        gti=gti,
+                                        power_type=None,
+                                        silent=silent,
+                                        dt=dt)
+            for key, val in spec.__dict__.items():
+                setattr(self, key, val)
+
+            return
+
+        self.type = "powerspectrum"
         self.dt = dt
 
         if isinstance(data, EventList):
@@ -401,8 +434,8 @@ class AveragedPowerspectrum(AveragedCrossspectrum, Powerspectrum):
         nphots_all : ``numpy.ndarray``
             List containing the number of photons for all segments calculated from ``lc``
         """
-        if not isinstance(lc, lightcurve.Lightcurve):
-            raise TypeError("lc must be a lightcurve.Lightcurve object")
+        if not isinstance(lc, Lightcurve):
+            raise TypeError("lc must be a Lightcurve object")
 
         current_gtis = lc.gti
 
@@ -435,9 +468,9 @@ class AveragedPowerspectrum(AveragedCrossspectrum, Powerspectrum):
                     "No counts in interval {}--{}s".format(time[0], time[-1]))
                 continue
 
-            lc_seg = lightcurve.Lightcurve(time, counts, err=counts_err,
-                                           err_dist=lc.err_dist.lower(),
-                                           skip_checks=True, dt=lc.dt)
+            lc_seg = Lightcurve(time, counts, err=counts_err,
+                                err_dist=lc.err_dist.lower(),
+                                skip_checks=True, dt=lc.dt)
 
             power_seg = Powerspectrum(lc_seg, norm=self.norm)
             power_all.append(power_seg)
@@ -538,11 +571,10 @@ class DynamicalPowerspectrum(AveragedPowerspectrum):
             bin_intervals_from_gtis(current_gti, self.segment_size, lc.time,
                                     dt=lc.dt)
 
-
         tstart = lc.time[start_inds]
         tend = lc.time[end_inds]
 
-        self.time = tstart + 0.5*(tend - tstart)
+        self.time = tstart + 0.5 * (tend - tstart)
 
         # Assign length of lightcurve as time resolution if only one value
         if len(self.time) > 1:
diff --git a/stingray/utils.py b/stingray/utils.py
index a79789a5..2fa3967b 100644
--- a/stingray/utils.py
+++ b/stingray/utils.py
@@ -1,16 +1,16 @@
-
+import numbers
+import os
+import random
+import string
 import sys
+import warnings
 from collections.abc import Iterable
-import numbers
-from six import string_types
 
-import warnings
 import numpy as np
 import scipy
 
 # If numba is installed, import jit. Otherwise, define an empty decorator with
 # the same name.
-
 HAS_NUMBA = False
 try:
     from numba import jit
@@ -97,7 +97,7 @@ __all__ = ['simon', 'rebin_data', 'rebin_data_log', 'look_for_array_in_array',
            'optimal_bin_time', 'contiguous_regions', 'is_int',
            'get_random_state', 'baseline_als', 'excess_variance',
            'create_window', 'poisson_symmetrical_errors', 'standard_error',
-           'nearest_power_of_two', 'find_nearest']
+           'nearest_power_of_two', 'find_nearest', 'genDataPath']
 
 
 def _root_squared_mean(array):
@@ -387,6 +387,7 @@ def apply_function_if_none(variable, value, func):
     else:
         return variable
 
+
 def assign_value_if_none(value, default):
     """
     Assign a value to a variable if that variable has value ``None`` on input.
@@ -803,7 +804,7 @@ def create_window(N, window_type='uniform'):
     windows = ['uniform', 'parzen', 'hamming', 'hanning', 'triangular',
                'welch', 'blackmann', 'flat-top']
 
-    if not isinstance(window_type, string_types):
+    if not isinstance(window_type, str):
         raise TypeError('type of window must be specified as string!')
 
     window_type = window_type.lower()
@@ -869,9 +870,9 @@ def create_window(N, window_type='uniform'):
         a3 = 0.388
         a4 = 0.028
         window = a0 - a1 * np.cos((2 * np.pi * n) / N_minus_1) + \
-                 a2 * np.cos((4 * np.pi * n) / N_minus_1) - \
-                 a3 * np.cos((6 * np.pi * n) / N_minus_1) + \
-                 a4 * np.cos((8 * np.pi * n) / N_minus_1)
+                a2 * np.cos((4 * np.pi * n) / N_minus_1) - \
+                a3 * np.cos((6 * np.pi * n) / N_minus_1) + \
+                a4 * np.cos((8 * np.pi * n) / N_minus_1)
 
     return window
 
@@ -998,3 +999,37 @@ def find_nearest(array, value):
         return array[idx - 1], idx - 1
     else:
         return array[idx], idx
+
+
+def genDataPath(dir_path):
+    """Generates data path to chunks.
+
+    Parameters
+    ----------
+    dir_path: string
+        Path to zarr datastore + Top level directory name for data
+
+    Returns
+    -------
+    list
+        List of path's to datastore
+
+    Raises
+    ------
+    IOError
+        If directory does not exist
+    """
+    path_list = []
+    if os.path.isdir(dir_path):
+        if not (os.path.isdir(os.path.join(dir_path, 'main_data/'))
+                or os.path.join(dir_path, 'meta_data/')):
+            raise IOError(("Directory does not exist."))
+
+        else:
+            path_list.append(os.path.join(dir_path, 'main_data/'))
+            path_list.append(os.path.join(dir_path, 'meta_data/'))
+
+            return path_list
+
+    else:
+        raise IOError(("Directory does not exist."))


Test patch:
diff --git a/stingray/tests/test_events.py b/stingray/tests/test_events.py
index aef93d8c..846d4135 100644
--- a/stingray/tests/test_events.py
+++ b/stingray/tests/test_events.py
@@ -236,11 +236,11 @@ class TestEvents(object):
         """Join two non-overlapping event lists.
         """
         ev = EventList(time=[1, 1, 10, 6, 5],
-                       energy=[10, 6, 3, 11, 2], gti=[[1, 3],[5, 6]],
+                       energy=[10, 6, 3, 11, 2], gti=[[1, 3], [5, 6]],
                        mjdref=57001)
         ev_other = EventList(time=np.asarray([5, 7, 6, 6, 10]) + 86400,
                              energy=[2, 3, 8, 1, 2],
-                             gti=np.asarray([[5, 7],[8, 10]]) + 86400,
+                             gti=np.asarray([[5, 7], [8, 10]]) + 86400,
                              mjdref=57000)
         ev_new = ev.join(ev_other)
 
@@ -304,4 +304,3 @@ class TestEvents(object):
 
         with pytest.raises(KeyError):
             ev.read('ev.pickle', format_="unsupported")
-
diff --git a/stingray/tests/test_io.py b/stingray/tests/test_io.py
index d2f045bf..be67921d 100644
--- a/stingray/tests/test_io.py
+++ b/stingray/tests/test_io.py
@@ -8,7 +8,7 @@ from ..io import read, write, split_numbers
 from ..io import ref_mjd
 from ..io import high_precision_keyword_read
 from ..io import load_events_and_gtis
-from ..io import read_header_key
+from ..io import read_header_key, _retrieve_ascii_object
 
 import warnings
 
@@ -192,6 +192,10 @@ class TestFileFormats(object):
               filename="ascii_test.txt", format_="ascii",
               fmt=["%s", "%s"])
 
+    def test_retrieve_bad(self):
+        with pytest.raises(TypeError):
+            _retrieve_ascii_object(1)
+
     def test_read_ascii(self):
         time = [1,2,3,4,5]
         counts = [5,7,8,2,3]
diff --git a/stingray/tests/test_largememory.py b/stingray/tests/test_largememory.py
new file mode 100644
index 00000000..f09ec1bc
--- /dev/null
+++ b/stingray/tests/test_largememory.py
@@ -0,0 +1,577 @@
+import copy
+import os
+import sys
+from sys import platform
+import tempfile
+
+import numpy as np
+import pytest
+from astropy.io import fits
+
+from stingray.crossspectrum import AveragedCrossspectrum
+from stingray.events import EventList
+from stingray.io import load_events_and_gtis, ref_mjd
+from stingray.largememory import retrieveData, saveData, genDataPath
+from stingray.largememory import _retrieveDataEV, _retrieveDataLC
+from stingray.lightcurve import Lightcurve
+from stingray.powerspectrum import AveragedPowerspectrum
+
+HAS_ZARR = False
+try:
+    import zarr
+
+    HAS_ZARR = True
+    from numcodecs import Blosc
+except ImportError:
+    pass
+
+curdir = os.path.abspath(os.path.dirname(__file__))
+datadir = os.path.join(curdir, "data")
+
+IS_LINUX = True
+if not (platform == "linux" or platform == "linux2"):
+    IS_LINUX = False
+
+class TestSaveSpec(object):
+    @classmethod
+    def setup_class(cls):
+        time = np.arange(0, 1e7)
+        counts = np.random.poisson(10, time.size)
+        cls.lc = Lightcurve(time, counts, skip_checks=True)
+
+        evtimes = np.sort(np.random.uniform(0, 1e7, 10 ** 7))
+        pi = np.random.randint(0, 100, evtimes.size)
+        energy = pi * 0.04 + 1.6
+        cls.ev = EventList(
+            time=evtimes,
+            pi=pi,
+            energy=energy,
+            gti=[[0, 1e7]],
+            dt=1e-5,
+            notes="Bu",
+        )
+
+    @pytest.mark.skipif("not HAS_ZARR")
+    def test_save_wrong_data(self):
+        with pytest.raises(ValueError) as excinfo:
+            saveData("A string", "bububu")
+        assert "Invalid data: A string (str)" in str(excinfo.value)
+
+    @pytest.mark.skipif("not HAS_ZARR")
+    def test_save_lc_small(self):
+        test_lc = copy.deepcopy(self.lc)
+        # Make sure counts_err exists
+        _ = test_lc.counts_err
+
+        # Save small part of data, < certainly chunk_size
+        _ = saveData(test_lc[:300], persist=False, chunks=100000)
+
+    @pytest.mark.skipif("not HAS_ZARR")
+    def test_save_lc(self):
+        test_lc = copy.deepcopy(self.lc)
+        # Make sure counts_err exists
+        _ = test_lc.counts_err
+
+        dir_name = saveData(test_lc, persist=False)
+
+        main = os.path.join(dir_name, "main_data")
+        meta = os.path.join(dir_name, "meta_data")
+
+        errors = []
+
+        if (
+            len([f for f in os.listdir(main) if not f.startswith(".")])
+            or len([f for f in os.listdir(meta) if not f.startswith(".")])
+        ) == 0:
+            errors.append("Lightcurve is not saved or does not exist")
+        else:
+            times = zarr.open_array(store=main, mode="r", path="times")[...]
+            counts = zarr.open_array(store=main, mode="r", path="counts")[...]
+            count_err = zarr.open_array(
+                store=main, mode="r", path="count_err"
+            )[...]
+            gti = zarr.open_array(store=main, mode="r", path="gti")[...]
+            gti = gti.reshape((gti.size // 2, 2))
+
+            dt = zarr.open_array(store=meta, mode="r", path="dt")[...]
+            mjdref = zarr.open_array(store=meta, mode="r", path="mjdref")[...]
+            err_dist = zarr.open_array(store=meta, mode="r", path="err_dist")[
+                ...
+            ]
+
+            if not np.array_equal(test_lc.time, times):
+                errors.append("lc.time is not saved precisely")
+            if not np.array_equal(test_lc.counts, counts):
+                errors.append("lc.counts is not saved precisely")
+            if not np.array_equal(test_lc.counts_err, count_err):
+                errors.append("lc.counts_err is not saved precisely")
+            if not np.array_equal(test_lc.gti, gti):
+                errors.append("lc.gti is not saved precisely")
+            if not (test_lc.dt == dt):
+                errors.append("lc.dt is not saved precisely")
+            if not (test_lc.mjdref == mjdref):
+                errors.append("lc.mjdref is not saved precisely")
+            if not (test_lc.err_dist == err_dist):
+                errors.append("lc.err_dist is not saved precisely")
+
+        assert not errors, "Errors encountered:\n{}".format("\n".join(errors))
+
+    @pytest.mark.skipif("not HAS_ZARR")
+    def test_save_ev_small(self):
+        # Save small part of data, < certainly chunk_size
+        ev = EventList(time=np.arange(1000))
+        _ = saveData(self.ev, persist=False, chunks=100000)
+
+    @pytest.mark.skipif("not (HAS_ZARR and IS_LINUX)")
+    def test_save_ev_missing_psutil_linux(self, monkeypatch):
+        monkeypatch.setitem(sys.modules, "psutil", None)
+
+        _ = saveData(self.ev, persist=False)
+
+    @pytest.mark.skipif("not HAS_ZARR or IS_LINUX")
+    def test_save_ev_missing_psutil_not_linux(self, monkeypatch):
+        monkeypatch.setitem(sys.modules, "psutil", None)
+
+        with pytest.warns(UserWarning) as record:
+            _ = saveData(self.ev, persist=False)
+        assert np.any(['will not depend on available RAM' in r.message.args[0]
+                       for r in record])
+
+    @pytest.mark.skipif("not HAS_ZARR")
+    def test_save_ev(self):
+        dir_name = saveData(self.ev, persist=False)
+
+        main = os.path.join(dir_name, "main_data")
+        meta = os.path.join(dir_name, "meta_data")
+
+        errors = []
+
+        if (
+            len([f for f in os.listdir(main) if not f.startswith(".")])
+            or len([f for f in os.listdir(meta) if not f.startswith(".")])
+        ) == 0:
+            errors.append("EventList is not saved or does not exist")
+
+        else:
+            times = zarr.open_array(store=main, mode="r", path="times")[...]
+            energy = zarr.open_array(store=main, mode="r", path="energy")[...]
+            pi_channel = zarr.open_array(
+                store=main, mode="r", path="pi_channel"
+            )[...]
+            gti = zarr.open_array(store=main, mode="r", path="gti")[...]
+            gti = gti.reshape((gti.size // 2, 2))
+            dt = zarr.open_array(store=meta, mode="r", path="dt")[...]
+            ncounts = zarr.open_array(store=meta, mode="r", path="ncounts")[
+                ...
+            ]
+            mjdref = zarr.open_array(store=meta, mode="r", path="mjdref")[...]
+            notes = zarr.open_array(store=meta, mode="r", path="notes")[...]
+
+            if not np.array_equal(self.ev.time, times):
+                errors.append("ev.time is not saved precisely")
+            if not np.array_equal(self.ev.energy, energy):
+                errors.append("ev.energy is not saved precisely")
+            if not np.array_equal(self.ev.pi, pi_channel):
+                errors.append("ev.pi is not saved precisely")
+            if not np.array_equal(self.ev.gti, gti):
+                errors.append("ev.gti is not saved precisely")
+            if not np.isclose(self.ev.dt, dt):
+                errors.append("ev.dt is not saved precisely")
+            if not self.ev.ncounts == ncounts:
+                errors.append("ev.ncounts is not saved precisely")
+            if not np.isclose(self.ev.mjdref, mjdref):
+                errors.append("ev.mjdref is not saved precisely")
+            if not self.ev.notes == notes:
+                errors.append("ev.notes is not saved precisely")
+
+        assert not errors, "Errors encountered:\n{}".format("\n".join(errors))
+
+    @pytest.mark.skipif("not HAS_ZARR")
+    def test_save_fits_data(self):
+        fname = os.path.join(datadir, "monol_testA.evt")
+        dir_name = saveData(fname, persist=False)
+
+        evtdata = load_events_and_gtis(fname, additional_columns=["PI"])
+        mjdref_def = ref_mjd(fname, hdu=1)
+        time_def = evtdata.ev_list
+        pi_channel_def = evtdata.additional_data["PI"]
+        gti_def = evtdata.gti_list
+        tstart_def = evtdata.t_start
+        tstop_def = evtdata.t_stop
+
+        main = os.path.join(dir_name, "main_data")
+        meta = os.path.join(dir_name, "meta_data")
+
+        errors = []
+
+        if (
+            len([f for f in os.listdir(main) if not f.startswith(".")])
+            or len([f for f in os.listdir(meta) if not f.startswith(".")])
+        ) == 0:
+            errors.append("EventList is not saved or does not exist")
+        else:
+            times = zarr.open_array(store=main, mode="r", path="times")[...]
+            pi_channel = zarr.open_array(
+                store=main, mode="r", path="pi_channel"
+            )[...]
+            gti = zarr.open_array(store=main, mode="r", path="gti")[...]
+            gti = gti.reshape((gti.size // 2, 2))
+            tstart = zarr.open_array(store=meta, mode="r", path="tstart")[...]
+            tstop = zarr.open_array(store=meta, mode="r", path="tstop")[...]
+            mjdref = zarr.open_array(store=meta, mode="r", path="mjdref")[...]
+
+            order = np.argsort(times)
+            times = times[order]
+            pi_channel = pi_channel[order]
+
+            if not np.allclose(time_def, times):
+                errors.append("fits.events.data.time is not saved precisely")
+            if not np.array_equal(pi_channel_def, pi_channel):
+                errors.append("fits.events.data.pi is not saved precisely")
+            if not np.allclose(gti_def, gti):
+                errors.append("fits.gti.data is not saved precisely")
+            if not (tstart == tstart_def):
+                errors.append(
+                    "fits.events.header.tstart is not saved precisely"
+                )
+            if not (tstop == tstop_def):
+                errors.append(
+                    "fits.events.header.tstop is not saved precisely"
+                )
+            if not (mjdref == mjdref_def):
+                errors.append(
+                    "fits.events.header.mjdref is not saved precisely"
+                )
+
+        assert not errors, "Errors encountered:\n{}".format("\n".join(errors))
+
+
+class TestRetrieveSpec(object):
+    @classmethod
+    def setup_class(cls):
+        time = np.arange(0, 1e7)
+        counts = np.random.poisson(10, time.size)
+        cls.lc = Lightcurve(time, counts, skip_checks=True)
+
+        evtimes = np.sort(np.random.uniform(0, 1e7, 10 ** 7))
+        pi = np.random.randint(0, 100, evtimes.size)
+        energy = pi * 0.04 + 1.6
+        cls.ev = EventList(
+            time=evtimes,
+            pi=pi,
+            energy=energy,
+            gti=[[0, 1e7]],
+            dt=1e-5,
+            notes="Bu",
+        )
+        cls.ev_noattrs = copy.deepcopy(cls.ev)
+        cls.ev_noattrs.energy = None
+        cls.ev_noattrs.pi = None
+        cls.ev_noattrs.mjdref = 0
+        cls.ev_noattrs.gti = None
+        cls.ev_noattrs.dt = 0
+        cls.ev_noattrs.notes = None
+
+        cls.lc_path = saveData(cls.lc, persist=False)
+        cls.ev_path = saveData(cls.ev, persist=False)
+        cls.ev_path_noattrs = saveData(cls.ev_noattrs, persist=False)
+
+    @pytest.mark.skipif("not HAS_ZARR")
+    def test_retrieve_wrong_data(self):
+        with pytest.raises(ValueError) as excinfo:
+            retrieveData("A string", self.lc_path)
+
+        assert "Invalid data: A string (str)" in str(excinfo.value)
+
+    @pytest.mark.skipif("not HAS_ZARR")
+    def test_retrieve_lc_data(self):
+        lc = retrieveData(data_type="Lightcurve", dir_path=self.lc_path)
+
+        assert self.lc.__eq__(lc) is True
+
+    @pytest.mark.skipif("not HAS_ZARR")
+    def test_retrieve_ev_data(self):
+        with pytest.warns(UserWarning):
+            ev = retrieveData(data_type="EventList", dir_path=self.ev_path)
+
+        gti = np.asarray([[self.ev.time[0], self.ev.time[-1]]])
+
+        assert np.allclose(ev.time, self.ev.time)
+        assert np.allclose(ev.pi, self.ev.pi)
+        assert np.allclose(ev.gti, gti, atol=0.001)
+
+    # @pytest.mark.skipif('not HAS_ZARR')
+    # def test_retrieve_fits_data(self):
+    #     fname = os.path.join(datadir, 'monol_testA.evt')
+    #     saveData(fname, os.path.join(self.path, self.dir))
+    #
+    #     with fits.open(fname) as fits_data:
+    #         time_def = fits_data[1].data['TIME']
+    #         gti_def = fits_data[2].data
+    #         for col in ['PI', 'PHA']:
+    #             if col in fits_data[1].data.columns.names:
+    #                 pi_channel_def = fits_data[1].data[col]
+    #
+    #         tstart_def = fits_data[1].header['TSTART']
+    #         tstop_def = fits_data[1].header['TSTOP']
+    #         mjdref_def = fits_data[1].header['MJDREF']
+    #
+    #
+    #     times, pi_channel, gti, tstart, tstop, mjdref = retrieveData(
+    #         data_type='FITS', dir_name=self.dir, path=self.path)
+    #
+    #     if not np.allclose(time_def, times):
+    #         errors.append("fits.events.data.time is not saved precisely")
+    #     if not np.array_equal(pi_channel_def, pi_channel):
+    #         errors.append("fits.events.data.pi is not saved precisely")
+    #     if not np.allclose(gti_def, gti):
+    #         errors.append("fits.gti.data is not saved precisely")
+    #     if not (tstart == tstart_def):
+    #         errors.append("fits.events.header.tstart is not saved precisely")
+    #     if not (tstop == tstop_def):
+    #         errors.append("fits.events.header.tstop is not saved precisely")
+    #     if not (mjdref == mjdref_def):
+    #         errors.append("fits.events.header.mjdref is not saved precisely")
+    #
+    #     assert not errors, "Errors encountered:\n{}".format('\n'.join(errors))
+
+    @pytest.mark.skipif("not HAS_ZARR")
+    def test_retrieve_lc_chunk_data(self):
+        lc = retrieveData(
+            data_type="Lightcurve",
+            dir_path=self.lc_path,
+            chunk_data=True,
+            chunk_size=10 ** 5,
+            offset=0,
+            raw=False,
+        )
+
+        trunc_lc = self.lc.truncate(stop=10 ** 5)
+
+        assert trunc_lc.__eq__(lc) is True
+
+    @pytest.mark.skipif("not HAS_ZARR")
+    def test_retrieve_ev_chunk_data(self):
+        maxidx = 10 ** 5
+        ev = retrieveData(
+            data_type="EventList",
+            dir_path=self.ev_path,
+            chunk_data=True,
+            chunk_size=maxidx,
+            offset=0,
+            raw=False,
+        )
+
+        gti = np.asarray([[self.ev.time[0], self.ev.time[-1]]])
+
+        assert np.allclose(ev.time, self.ev.time[:maxidx])
+        assert np.allclose(ev.pi, self.ev.pi[:maxidx])
+        assert np.allclose(ev.gti, gti, atol=0.0001)
+
+    @pytest.mark.skipif("not HAS_ZARR")
+    def test_retrieve_lc_offset_data(self):
+        lc = retrieveData(
+            data_type="Lightcurve",
+            dir_path=self.lc_path,
+            chunk_data=True,
+            chunk_size=10 ** 5,
+            offset=10 ** 2,
+            raw=False,
+        )
+
+        trunc_lc = self.lc.truncate(start=10 ** 2, stop=10 ** 5)
+
+        assert trunc_lc.__eq__(lc) is True
+
+    @pytest.mark.skipif("not HAS_ZARR")
+    def test_retrieve_ev_offset_data(self):
+        maxidx = 10 ** 5
+        offset = 100
+        ev = retrieveData(
+            data_type="EventList",
+            dir_path=self.ev_path,
+            chunk_data=True,
+            chunk_size=maxidx,
+            offset=offset,
+            raw=False,
+        )
+
+        gti = np.asarray([[self.ev.time[0], self.ev.time[-1]]])
+
+        assert np.allclose(ev.time, self.ev.time[offset:maxidx])
+        assert np.allclose(ev.pi, self.ev.pi[offset:maxidx])
+        assert np.allclose(ev.gti, gti, atol=0.1)
+
+    @pytest.mark.skipif("not HAS_ZARR")
+    def test_retrieve_data_bad_offset(self):
+        with pytest.raises(ValueError):
+            _ = retrieveData(data_type="EventList", dir_path=self.ev_path,
+                             chunk_data=True,
+                             offset=101010101010)
+
+    @pytest.mark.skipif("not HAS_ZARR")
+    def test_retrieve_ev_data_bad_offset(self):
+        with pytest.raises(ValueError) as excinfo:
+            _ = _retrieveDataEV(data_path=genDataPath(self.ev_path),
+                                offset=101010101010)
+        assert "Offset cannot be larger than size of" in str(excinfo.value)
+
+    @pytest.mark.skipif("not HAS_ZARR")
+    def test_retrieve_lc_data_bad_offset(self):
+        with pytest.raises(ValueError) as excinfo:
+            _ = _retrieveDataLC(data_path=genDataPath(self.lc_path),
+                                offset=101010101010)
+        assert "Offset cannot be larger than size of " in str(excinfo.value)
+
+    @pytest.mark.skipif("not HAS_ZARR")
+    def test_retrieve_ev_data_no_attrs(self):
+        ev = _retrieveDataEV(data_path=genDataPath(self.ev_path_noattrs))
+        assert ev.mjdref == 0
+        assert ev.pi is None
+        assert ev.energy is None
+        assert ev.gti is None
+        assert ev.notes == ""
+
+    @pytest.mark.skipif("not HAS_ZARR")
+    def test_retrieve_ev_data_raw(self):
+        res = _retrieveDataEV(data_path=genDataPath(self.ev_path_noattrs),
+                              raw=True)
+        times, energy, ncounts, mjdref, dt, gti, pi, notes = res
+        assert mjdref == 0
+        assert pi is None
+        assert energy is None
+        assert gti is None
+        assert notes == ""
+
+class TestChunkPS(object):
+    @classmethod
+    def setup_class(cls):
+        maxtime = 2 ** 21
+        time = np.arange(maxtime)
+        counts1 = np.random.poisson(10, time.size)
+        cls.lc1 = Lightcurve(
+            time, counts1, skip_checks=True, gti=[[0, maxtime]]
+        )
+
+        counts2 = np.random.poisson(10, time.size)
+        cls.lc2 = Lightcurve(
+            time, counts2, skip_checks=True, gti=[[0, maxtime]]
+        )
+
+    @pytest.mark.skipif("not HAS_ZARR")
+    def test_invalid_data_to_pds(self):
+        with pytest.raises(ValueError) as excinfo:
+            AveragedPowerspectrum(
+                "sdfasfsa", segment_size=2048, large_data=True, silent=True
+            )
+        assert "Invalid input data type: str" in str(excinfo.value)
+
+    @pytest.mark.skipif("not HAS_ZARR")
+    def test_invalid_data_to_cpds(self):
+        with pytest.raises(ValueError) as excinfo:
+            AveragedCrossspectrum(
+                "sdfasfsa",
+                "sdfasfsa",
+                segment_size=4096,
+                large_data=True,
+                silent=True,
+            )
+        assert "Invalid input data type: str" in str(excinfo.value)
+
+    @pytest.mark.skipif("not HAS_ZARR")
+    def test_calc_pds(self):
+        ps_normal = AveragedPowerspectrum(
+            self.lc1, segment_size=8192, silent=True, norm="leahy"
+        )
+        ps_large = AveragedPowerspectrum(
+            self.lc1,
+            segment_size=8192,
+            large_data=True,
+            silent=True,
+            norm="leahy",
+        )
+
+        attrs = [
+            "freq",
+            "power",
+            "power_err",
+            "unnorm_power",
+            "df",
+            "n",
+            "nphots",
+            "gti",
+            "m",
+        ]
+        allgood = True
+        assert ps_normal.freq.size == ps_large.freq.size
+        for attr in attrs:
+            if not np.allclose(
+                getattr(ps_normal, attr),
+                getattr(ps_large, attr),
+                rtol=0.1,
+                atol=0.1,
+            ):
+                allgood = False
+                print(f"Attribute = {attr} ")
+                print(
+                    f"Raw Array: \nOriginal: {getattr(ps_normal, attr)}, "
+                    f"\nLarge: {getattr(ps_large, attr)}"
+                )
+                maxdev = np.amax(
+                    getattr(ps_normal, attr) - getattr(ps_large, attr)
+                )
+                maxdev_percent = np.abs(
+                    np.max(getattr(ps_normal, attr) - getattr(ps_large, attr))
+                    * 100
+                ) / np.max(getattr(ps_normal, attr))
+                print(f"Max Deviation: {maxdev}, as %: {maxdev_percent}")
+                print("\n")
+        assert allgood
+
+    @pytest.mark.skipif("not HAS_ZARR")
+    def test_calc_cpds(self):
+        cs_normal = AveragedCrossspectrum(
+            self.lc1, self.lc2, segment_size=8192, silent=True
+        )
+        cs_large = AveragedCrossspectrum(
+            self.lc1, self.lc2, segment_size=8192, large_data=True, silent=True
+        )
+
+        attrs = [
+            "freq",
+            "power",
+            "power_err",
+            "unnorm_power",
+            "df",
+            "n",
+            "nphots1",
+            "nphots2",
+            "m",
+            "gti",
+        ]
+        assert cs_normal.freq.size == cs_large.freq.size
+
+        allgood = True
+        for attr in attrs:
+            if not np.allclose(
+                getattr(cs_normal, attr),
+                getattr(cs_large, attr),
+                rtol=0.1,
+                atol=0.1,
+            ):
+                print(f"Attribute = {attr} ")
+                print(
+                    f"Raw Array: \nOriginal: {getattr(cs_normal, attr)}, \n"
+                    f"Large: {getattr(cs_large, attr)}"
+                )
+                maxdev = np.amax(
+                    getattr(cs_normal, attr) - getattr(cs_large, attr)
+                )
+                maxdev_percent = np.abs(
+                    np.max(getattr(cs_normal, attr) - getattr(cs_large, attr))
+                    * 100
+                ) / np.max(getattr(cs_normal, attr))
+                print(f"Max Deviation: {maxdev}, as %: {maxdev_percent}")
+                print("\n")
+                allgood = False
+        assert allgood