diff --git a/.gitignore b/.gitignore index 018dd29..8d1117c 100644 --- a/.gitignore +++ b/.gitignore @@ -52,6 +52,7 @@ cover/ # Sphinx documentation docs/_build/ +docs/source/images/ # Jupyter Notebook .ipynb_checkpoints @@ -71,4 +72,4 @@ venv.bak/ ## Misc .DS_Store -.vscode \ No newline at end of file +.vscode diff --git a/README.rst b/README.rst index c4ddcd6..1df3fde 100644 --- a/README.rst +++ b/README.rst @@ -13,8 +13,7 @@ to fit heavy-tailed distributions like power laws. Academics, please cite as: Jeff Alstott, Ed Bullmore, Dietmar Plenz. (2014). powerlaw: a Python package for analysis of heavy-tailed distributions. `PLoS ONE 9(1): e85777 `_ - -Also available at `arXiv:1305.0215 [physics.data-an] `_ +(also available at `arXiv:1305.0215 [physics.data-an] `_) Basic Usage @@ -79,9 +78,8 @@ Alternatively, you can install directly from the source: This library depends on the usual scientific computing libraries that you probably already have installed: ``numpy``, ``scipy``, ``matplotlib``, and -``mpmath``. - -The package ``tqdm`` is used for creating progress bars. +``mpmath``, as well as ``dill`` and ``h5py`` for caching objects and ``tqdm`` +for creating progress bars. The requirement of ``mpmath`` will be dropped if/when the scipy functions ``gamma``, ``gammainc`` and ``gammaincc`` are updated to have sufficient numerical @@ -119,6 +117,7 @@ their code available. Their implementations were a critical starting point for making ``powerlaw``. + Power Laws vs. Lognormals and powerlaw's 'lognormal_positive' option -------------------------------------------------------------------- When fitting a power law to a data set, one should compare the goodness of fit to that of a `lognormal distribution `__. This is done because lognormal distributions are another heavy-tailed distribution, but they can be generated by a very simple process: multiplying random positive variables together. The lognormal is thus much like the normal distribution, which can be created by adding random variables together; in fact, the log of a lognormal distribution is a normal distribution (hence the name), and the exponential of a normal distribution is the lognormal (which maybe would be better called an expnormal). In contrast, creating a power law generally requires fancy or exotic generative mechanisms (this is probably why you're looking for a power law to begin with; they're sexy). So, even though the power law has only one parameter (``alpha``: the slope) and the lognormal has two (``mu``: the mean of the random variables in the underlying normal and ``sigma``: the standard deviation of the underlying normal distribution), we typically consider the lognormal to be a simpler explanation for observed data, as long as the distribution fits the data just as well. For most data sets, a power law is actually a worse fit than a lognormal distribution, or perhaps equally good, but rarely better. This fact was one of the central empirical results of the paper `Clauset et al. 2007 `__, which developed the statistical methods that ``powerlaw`` implements. diff --git a/docs/source/index.rst b/docs/source/index.rst index 34f6c8c..bfe1291 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -20,7 +20,7 @@ package for analysis of heavy-tailed distributions. PLoS ONE 9(1): e85777 Quick links ----------- - Original paper for the library: http://arxiv.org/abs/1305.0215 -- Source code: https://github.com/jeffalstott/powerlaw +- Source code: https://github.com/powerlaw-devs/powerlaw Installation @@ -36,21 +36,20 @@ Alternatively, you can install directly from the source: .. code-block:: console - $ git clone https://github.com/jeffalstott/powerlaw + $ git clone https://github.com/powerlaw-devs/powerlaw $ cd powerlaw $ pip install . This library depends on the usual scientific computing libraries that you probably already have installed: ``numpy``, ``scipy``, ``matplotlib``, and -``mpmath``. - -The package ``tqdm`` is used for creating progress bars. +``mpmath``, as well as ``dill`` and ``h5py`` for caching objects and ``tqdm`` +for creating progress bars. The requirement of ``mpmath`` will be dropped if/when the scipy functions ``gamma``, ``gammainc`` and ``gammaincc`` are updated to have sufficient numerical accuracy for negative numbers. -See the `powerlaw home page `_ for more +See the `powerlaw home page `_ for more information and examples. Basic usage diff --git a/docs/source/tutorials/saving_fits.rst b/docs/source/tutorials/saving_fits.rst new file mode 100644 index 0000000..085b6e4 --- /dev/null +++ b/docs/source/tutorials/saving_fits.rst @@ -0,0 +1,186 @@ +Saving and loading fits +======================= + +After choosing parameters, calculating ``xmin`` values, or fitting various +distributions, you might want to save the ``Fit`` object to a file. This +allows you to easily load it in during a future session, without having +to recalculate everything (particularly the ``xmin`` value, which is usually +somewhat computationally expensive). + +This can be done easily using :meth:`Fit.save` and :meth:`Fit.load`. + +.. code-block:: + + data = [1.1, 5.3, 3.7, ...] + fit = powerlaw.Fit(data, xmin=0.1) + + fit.save('output.h5') + +.. code-block:: + + # In another session + + fit = powerlaw.Fit.load('output.h5') + + fit.plot_pdf() + ... + +The saving and loading functions currently support two different file formats: +`pickle `_ and +`hdf5 `_. +A pickle file is Python's way of serializing an object, +which saves the entire object to a file that can then be loaded later. The +hdf5 format is a more universal format that allows you to save numerical data +alongside various metadata. This means that the hdf5 file doesn't contain the +actual ``Fit`` object like the pickle file does, but instead saves all of +the important information and then reconstructs the ``Fit`` when you load it +back in. Pickling is done using the `dill `_ +library (which improves on the standard library ``pickle``) and hdf5 file +operations are done using the `h5py `_ library. + +You can choose which format to use by either including it in the filename, +or with the ``format`` keyword: + +.. code-block:: + + fit.save('output.h5') # saves in hdf5 format + fit.save('output', format='h5') # saves in hdf5 format + fit.save('output.pkl') # saves in pickle format + fit.save('output', format='pkl') # saves in pickle format + +If you're just working with the ``powerlaw`` library, these two formats are +almost entirely interchangeable, with hdf5 files being slightly smaller than +pickle files. That being said, hdf5 files do have the advantage of being +easily read and interpreted outside of this library, or even outside of +Python altogether. If you're worried about future-proofing your data, or +want to use this data in other programming languages, hdf5 is probably better. + + +Automatic caching +----------------- + +``powerlaw`` offers the option automatically cache *all* fits, if you don't +want to have to manually save files. This is disabled by default, but can +be enabled by setting the cache directory with :meth:`powerlaw.Fit.set_cache_folder()`. + +.. code-block:: + + powerlaw.Fit.set_cache_folder('data/') + +For the rest of the session, all ``Fit`` objects will automatically be +saved in this folder after creation. And if you create a ``Fit`` object that +is identical to a cached one, it will be loaded instead of recalculating things. +This might be useful if you are working on a project where you are consistently +working with several predefined datasets, and you don't want to have to, for +example, recalculate ``xmin`` during each session. + +.. code-block:: + + powerlaw.Fit.set_cache_folder('data/') + + data = np.genfromtxt('data.txt') + + # This will calculate xmin, and then cache the object + fit = powerlaw.Fit(data) + +.. code-block:: + + # In another session + + powerlaw.Fit.set_cache_folder('data/') + + # The same data as before + data = np.genfromtxt('data.txt') + + # This will just load the previously cached file + fit = powerlaw.Fit(data) + +This replacement only happens when the data and all of the parameters of +fitting are exactly the same. + +.. code-block:: + + # In another session + + # The same data as before + data = np.genfromtxt('data.txt') + + # This will *not* load the previously cached file since xmin is different + fit = powerlaw.Fit(data, xmin=1) + + +A note on constraints +--------------------- + +Constraint functions are a little tricky to save since they might have +dependencies on variables, functions or libraries beyond the function itself. +For example, the following constraint could very likely give an error: + +.. code-block:: + + import powerlaw + import numpy as np + + data = np.genfromtxt('data.txt') + + def constraint(dist): + """ + Some constraint that depends on the library numpy + """ + E = np.exp(...) + ... + + constraint_dict = {"type": 'eq', + "fun": constraint} + + fit = powerlaw.Fit(data, parameter_constraints=constraint_dict) + + fit.save('output.h5') + +.. code-block:: + + # In another session + + import powerlaw + # numpy is *not* imported + + fit = powerlaw.Fit.load('output.h5') + + constraint = fit.parameter_constraints[0]["fun"] + + # This will give an error that the function can't find numpy since we + # haven't imported it. + constraint(...) + +The best practice here is to have constraint functions be fully self contained, +including definitions of variables and library imports. + +.. code-block:: + + # Best practice: fully self-contained + def constraint(dist): + import numpy as np + + T = 100 + E = np.exp(-dist.Lambda * T) + ... + +.. code-block:: + + # Not good practice but will still work + T = 100 + def constraint(dist): + import numpy as np + + E = np.exp(-dist.Lambda * T) + ... + +.. code-block:: + + # Will not work! + import numpy as np + T = 100 + def constraint(dist): + + E = np.exp(-dist.Lambda * T) + ... diff --git a/docs/source/tutorials_top.rst b/docs/source/tutorials_top.rst index 6cdd864..bb95e63 100644 --- a/docs/source/tutorials_top.rst +++ b/docs/source/tutorials_top.rst @@ -12,6 +12,7 @@ Tutorials tutorials/discrete_continuous tutorials/ranges_and_constraints tutorials/comparing_distributions + tutorials/saving_fits tutorials/generating_data tutorials/advanced_topics tutorials/warnings diff --git a/powerlaw/distributions.py b/powerlaw/distributions.py index 0f05a60..c7fd390 100644 --- a/powerlaw/distributions.py +++ b/powerlaw/distributions.py @@ -320,7 +320,10 @@ def initialize_parameters(self, initial_parameters=None, **kwargs): # (since we may be given parameters from a fit object that # contains information for other distributions). for k,v in initial_parameters.items(): - if k in self.parameter_names: + # It's possible that we could be passed an empty initial parameter + # value, so we have to check that we actually have a valid + # value. + if k in self.parameter_names and v not in [None, np.nan]: initial_parameters_dict[k] = v elif hasattr(initial_parameters, '__iter__') and len(initial_parameters) == len(self.parameter_names): diff --git a/powerlaw/fitting.py b/powerlaw/fitting.py index 854c663..e076998 100644 --- a/powerlaw/fitting.py +++ b/powerlaw/fitting.py @@ -19,6 +19,22 @@ # So we can ignore this warning while fitting xmin from scipy.optimize import OptimizeWarning +# So we can hash the fit object for caching +import hashlib + +# For saving and loading +import h5py +# For saving function source code when saving/loading +import inspect +import json +import textwrap + +# Dill is a drop-in replacement for pickle that allows you to serialize +# local functions (among other improvements); used for saving and loading +# fit objects, and parallelization. +#import pickle +import dill as pickle + import warnings from tqdm import tqdm @@ -26,6 +42,18 @@ from .statistics import * from .distributions import * +# Try and grab the current version of the package from the _version.py file. +# This is for comparing with saved/loaded Fit objects. If we can't find that +# file, we have to work without it. +try: + from ._version import __version__ + POWERLAW_VERSION = __version__ +except: + # Raise a warning, in case the user isn't aware that the package + # isn't properly installed + warnings.warn('powerlaw version not found, likely because the package isn\'t installed properly. Not critical, but could affect caching of files.') + POWERLAW_VERSION = 'unknown' + # This needs to be a list of the keys in the supported_distributions # attribute of the Fit class. The __getattr__ method needs the list. # If it uses supported_distributions.keys(), then it gets into an @@ -42,20 +70,106 @@ SUPPORTED_DISTRIBUTION_LIST = list(SUPPORTED_DISTRIBUTIONS.keys()) """ -Currently just templated; doesn't work yet. - Whether to enable parallelization for certain heavy calculations, eg. fitting the xmin value. """ -PARALLEL_ENABLE = False +_parallel_enable = False + +""" +This is the number of cores/processes that the library should use. +""" +_parallel_cores = 1 + +""" +By default, the multiprocessing library has many limitations when it comes +to what functions can be used in a Pool or Process object. This is becuase +these require that you serialize the function using pickle, which is then +pass between processes. + +Unfortunately, the standard library `pickle` doesn't have support for many +types of functions, including those not defined at the root level. The +function we want to parallelize for `find_xmin` is not defined at the +root level, so normally we would get an error along the lines of: + Can't pickle local object 'Fit.find_xmin..fit_function' + +We can fix this by manually replacing multiprocessing's use of pickle to +`dill`'s version of it, which can handle these types of functions. +""" +# We imported dill as 'pickle', so all references to 'pickle' are actually +# dill. +pickle.Pickler.dumps, pickle.Pickler.loads = pickle.dumps, pickle.loads +multiprocessing.reduction.ForkingPickler = pickle.Pickler +multiprocessing.reduction.dump = pickle.dump + +""" +The file types and extensions supported for saving and loading files. + +Each key is the name of the format, and the value should be the list of +extensions that correspond to this format. +""" +SUPPORTED_SAVE_FORMATS = {"hdf5": ['h5', 'hdf5'], + "pickle": ['pkl', 'pickle']} + +SUPPORTED_SAVE_FILE_EXTENSIONS = [ext for v in SUPPORTED_SAVE_FORMATS.values() for ext in v] + +DEFAULT_SAVE_FORMAT = 'h5' + +""" +Whether fits are cached automatically or not. + +This variable should not be manually set, but rather is automatically +set to True when a valid path is given using Fit.set_cache_folder(). +""" +_cache_enabled = False + """ -Currently just templated; doesn't work yet. +The path to the cache folder. -This is the number of cores that the library should leave free when doing -certain heavy calculations. For example, if you have 8 cores, and this is -set to 2, then the processing would use (up to) 6 cores. +This variable should not be manually set, but rather set using +Fit.set_cache_folder(). """ -PARALLEL_UNUSED_CORES = 2 +_cache_path = None + +""" +The format to use when automatically caching files. + +This variable should not be manually set, but rather set using +Fit.set_cache_format(). The default is hdf5. +""" +_cache_format = DEFAULT_SAVE_FORMAT + + +def set_parallel_cores(num_cores): + """ + Set the number of cores to use in parallel for computing xmin. + + If a negative number, will leave that many cores open. + """ + # See the os documentation on for this below; note that newer + # versions of python (3.13+) have a function `process_cpu_count()` + # that does this in a cleaner way, but it's probably better + # to be backwards compatible. + # https://docs.python.org/3.11/library/os.html#os.cpu_count + total_cores = len(os.sched_getaffinity(0)) + + global _parallel_enable, _parallel_cores + + if num_cores < 0: + usable_cores = total_cores + num_cores + usable_cores = max(usable_cores, 1) + else: + if num_cores > total_cores: + raise ValueError(f'Attempted to parallelize using {num_cores} cores, but only {total_cores} are available.') + + usable_cores = num_cores + + _parallel_cores = usable_cores + + if _parallel_cores == 0 or _parallel_cores == 1: + _parallel_enable = False + + else: + _parallel_enable = True class Fit(object): @@ -172,6 +286,10 @@ def constraint(dist): xmin (True) or generate a uniform distribution of values that spans the data. + ignore_cache : bool + Whether to ignore cached files, even if automatic cacheing is + enabled. + verbose: {0, 1, 2} or bool, optional Whether to print updates about where we are in the fitting process. @@ -195,6 +313,7 @@ def __init__(self, xmin_distance='D', xmin_distribution='power_law', test_all_xmin=False, + ignore_cache=False, verbose=1): self.verbose = verbose @@ -212,9 +331,29 @@ def __init__(self, self.discrete_normalization = discrete_normalization self.sigma_threshold = sigma_threshold + # For initial parameters and ranges, we need to do some standardization + # such that we can nicely save and load Fit objects to files. This + # is primarily relevant to the hashing, since we need to make sure + # the type of numbers in either case is a float or None. + if hasattr(initial_parameters, '__iter__'): + for k, v in initial_parameters.items(): + initial_parameters[k] = float(v) if v else None + self.initial_parameters = initial_parameters + + if hasattr(parameter_ranges, '__iter__'): + for k, v in parameter_ranges.items(): + parameter_ranges[k][0] = float(v[0]) if v[0] else None + parameter_ranges[k][1] = float(v[1]) if v[1] else None + self.parameter_ranges = parameter_ranges - self.parameter_constraints = parameter_constraints + + # Make sure that we are given a list of constraints; so even if we + # are only given a single one, added it to a list + if type(parameter_constraints) is dict: + self.parameter_constraints = [parameter_constraints] + else: + self.parameter_constraints = parameter_constraints # We keep track of the xmin and xmax values if they are provided. # I don't really see the purpose for this variable, but I'll leave @@ -282,7 +421,37 @@ def __init__(self, self.supported_distributions = SUPPORTED_DISTRIBUTIONS self.xmin_distribution_cls: type[Distribution] = self.supported_distributions[xmin_distribution] - + # We need to save this for hashing; otherwise, we never use it + # again after passing to find_xmin(). + self.test_all_xmin = test_all_xmin + + #################################### + # CHECK CACHE + # Now that we have all of our variables set, we should check to see + # if we need to look for a cached copy of this fit. Notably, this + # happens before we fit xmin, since the point of caching is to + # avoid having to repeat that calculation if possible. + if _cache_enabled and not ignore_cache: + potential_cache_file = os.path.join(_cache_path, f'{hash(self)}.{_cache_format}') + + if os.path.exists(potential_cache_file): + if self.verbose: + print(f'Found cached file: {potential_cache_file}') + + # You can't just overwrite self, so we have to update + # every variable + loaded_fit = Fit.load(potential_cache_file) + print(self.__dict__) + self.__dict__.update(loaded_fit.__dict__) + print(self.__dict__) + return + + # If we don't find a cache file, no problem, we just continue on with + # calculating xmin. + + + #################################### + # FIT XMIN # If we have a fixed xmin, we can directly fit a power law distribution if self.fixed_xmin: self.xmin = float(xmin) @@ -291,13 +460,24 @@ def __init__(self, print(f'Calculating best minimal value for {xmin_distribution.replace("_"," ")} fit') # This function tries to optimize the fit based on the xmin - self.find_xmin(xmin_distance, test_all_xmin) + self.find_xmin(self.xmin_distance, self.test_all_xmin) # Crop the data to the xmin and self.data = self.data[self.data >= self.xmin] self.n = float(len(self.data)) self.n_tail = self.n + n_above_max + #################################### + # CREATE CACHE + if _cache_enabled and not ignore_cache: + new_cache_file = os.path.join(_cache_path, f'{hash(self)}.{_cache_format}') + + if not os.path.exists(new_cache_file): + if verbose: + print(f'Caching file at: {new_cache_file}') + self.save(new_cache_file) + + def __dir__(self): """ @@ -351,6 +531,7 @@ def __getattr__(self, name): else: raise AttributeError(name) + @property def xmin_distribution(self): return getattr(self, self.xmin_distribution_cls.name) @@ -408,10 +589,10 @@ def find_xmin(self, xmin_distance=None, test_all_xmin=False): else: max_bin_value = np.sort(possible_xmin)[-3] - # 10% of the number of datapoints sounds good. And note that we - # only generate bins up into the 3rd to last point so we always - # have enough points to calculate distance metrics. - num_bins = max(100, len(self.data) // 10) + # Use 1% of the datapoints, but always at least 100 if the 1% + # is less than that. This should work well for most cases, but + # might need to be improved in the future. + num_bins = max(100, len(self.data) // 100) # These are logarithmically spaced possible_xmin = np.logspace(np.log10(np.min(self.data)), np.log10(max_bin_value), num_bins) @@ -476,34 +657,39 @@ def fit_function(xmin): warnings.filterwarnings('ignore', category=OptimizeWarning) warnings.filterwarnings('ignore', category=UserWarning) - # TODO parallelize - # This is slightly harder than I thought it would be since I can't - # directly parallelize a function that isn't defined at the top - # level. The alternative is to use a third-party library like - # multiprocess but this is something to ask. - if PARALLEL_ENABLE: - raise NotImplementedError('Parallelization not yet implemented! Use `PARALLEL_ENABLE=False`') - - # See the os documentation on for this below; note that newer - # versions of python (3.13+) have a function `process_cpu_count()` - # that does this in a cleaner way, but it's probably better - # to be backwards compatible. - # https://docs.python.org/3.11/library/os.html#os.cpu_count - usable_cores = len(os.sched_getaffinity(0)) - PARALLEL_UNUSED_CORES - usable_cores = max(usable_cores, 1) + # TODO improve parallelization + # This currently works, but we see less time reduction than + # I would expect. For example, I would expect that using 4 + # processes/cores would give something like a speedup of close + # 4 (though definitely less), but it usually ends up just around + # 2x. This means either the majority of the time is spent + # communicating between processes (eg. writing arrays), or + # something else is weird. + + global _parallel_enable, _parallel_cores + if _parallel_enable: + #raise NotImplementedError('Parallelization not yet implemented! Use `PARALLEL_ENABLE=False`') # Create a pool of workers - with multiprocessing.Pool(usable_cores) as pool: + with multiprocessing.Pool(_parallel_cores) as pool: # We don't need the xmin to be tested in order, so we # use an unordered map - result_mapping = pool.imap_unordered(fit_function, possible_xmin) - for result in tqdm(result_mapping, desc="Fitting xmin") if self.verbose else result_mapping: - distances[i], valid_fits[i], params[:, i] = result + + # chunksize controls how many iterations a process works + # on before communicating back to the main thread. I + # experimented with this a bit, but I think we can get + # better results by choosing this well. + #chunksize = max(int(len(possible_xmin) / PARALLEL_CORES / 100), 1) + chunksize = 1 + + result_mapping = enumerate(pool.imap_unordered(fit_function, possible_xmin, chunksize=chunksize)) + for i, result in tqdm(result_mapping, desc="Fitting xmin") if self.verbose else result_mapping: + distances[i], valid_fits[i], params[:,i] = result else: # For non-parallel case, we just use a simple for loop for i in tqdm(range(num_xmin), desc='Fitting xmin') if self.verbose else range(num_xmin): - distances[i], valid_fits[i], params[:, i] = fit_function(possible_xmin[i]) + distances[i], valid_fits[i], params[:,i] = fit_function(possible_xmin[i]) # The possible xmin values should of course have all parameters @@ -858,4 +1044,540 @@ def plot_pdf(self, original_data=False, linear_bins=False, bins=None, ax=None, * return plot_pdf(data, xmin=xmin, xmax=xmax, linear_bins=linear_bins, bins=bins, ax=ax, **kwargs) + def save(self, filename, format=None): + """ + Save the fit object to a file. + + Note that this doesn't use Python's pickling framework, but instead + saves the relevant data and fitting information in an hdf5 file. + This way, the data can easily be recovered even if you aren't + working with this library or even with Python. The cost of this is that + the saving and loading methods are relatively complex (or at least + just long) since we have to parse all of the important information + into generic formats. As an end-user, this isn't a problem at all, + but might make maintenance slightly more difficult. + + Note that saving and loading imposes some restrictions on the + form of constraint functions. This saving and loading is done by + using the actual source code of the constraint functions, which + means that each function must be fully self-contained. For more + information, see the tutorial page about parameter constraints. + + The current version of ``powerlaw`` will be saved within the file; + if you try to load a file from a different version, you will be shown + a warning. + + Parameters + ---------- + filename : str or Path + The file to which to save the ``Fit`` data. + """ + + # Determine what type of file we have + file_extension = filename.split('.')[-1] + + # If we have no extension, we follow the value of the format kw + if '.' not in filename or file_extension not in SUPPORTED_SAVE_FILE_EXTENSIONS: + + if format is not None and format in SUPPORTED_SAVE_FILE_EXTENSIONS: + full_filename = filename + '.' + format + + elif format is not None: + raise ValueError('Desired format ({format}) is unsupported.') + + else: + # Just use the default format + full_filename = filename + '.' + DEFAULT_SAVE_FORMAT + format = DEFAULT_SAVE_FORMAT + + else: + full_filename = filename + format = file_extension + + if format in ['h5', 'hdf5']: + _write_hdf5_file(full_filename, self) + + elif format in ['pkl', 'pickle']: + _write_pickle_file(full_filename, self) + + else: + raise ValueError('Unsupported save format!') + + + @staticmethod + def load(filename, verbose=1): + """ + Load a saved fit object. + + See also ``powerlaw.Fit.save()``. + + Note that the loading and saving imposes some restrictions on the + form of constraint functions. This saving and loading is done by + using the actual source code of the constraint functions, which + means that each function must be fully self-contained. For more + information, see the tutorial page about parameter constraints. + + This function can load the following types, signified by their + respective extensions: + + HDF5: .h5, .hdf5 + Pickle: .pkl, .pickle + Parameters + ---------- + filename : str or Path + The path to the file that will be loaded. + + Returns + ------- + fit : Fit + The loaded fit object, which should be functionally identical + to the saved object. + """ + # Determine what type of file we have + file_extension = filename.split('.')[-1] + + # HDF5 + if file_extension.lower() in ['h5', 'hdf5']: + return _parse_hdf5_file(filename, verbose) + + # Pickle + elif file_extension.lower() in ['pkl', 'pickle']: + return _parse_pickle_file(filename) + + else: + raise ValueError(f'Unknown filetype passed ({filename}); make sure that your file has an appropriate extension. See documentation for this function for acceptable extensions.') + + + def __eq__(self, other): + """ + Check for equality between two ``Fit`` objects. + + This is done using a custom implementation that checks the actual + fit data and parameters against each other, such that two separate + instances with the exact same information will be deemed equal. + """ + return hash(self) == hash(other) + + + def __hash__(self): + """ + Generate a unique hash for this fit based on the data and other + fitting parameters, such that the hash is the same for any two + identical fits. + + Used for automatic caching of ``Fit`` objects. + + Note that this doesn't do any checks about specific distributions, + eg. `fit.power_law`; this is because these may or may not exist + depending on whether the user has accessed them. + + Returns + ------- + + hash : int + The integer hash of the ``Fit`` object + """ + # We cast the data to a specific type, since it could be passed + # as a float32, float64, int32, int64, etc. + # We choose float32 since it's not too large, but includes enough + # precision. + retyped_data = np.sort(np.array(self.data_original, dtype=np.float32)) + # Numpy arrays are unhashable, so we need a fixed tuple. + retyped_data = tuple(retyped_data) + + # We need to be able to create this hash *before* computing xmin + # (if requested), so we can't use the actual xmin value in the hashing + # if one isn't explicitly specified. + xmin_value = np.float32(self.xmin if self.fixed_xmin else -1) + xmax_value = np.float32(self.xmax if self.xmax else -1) + + # For the parameter ranges and initial values, we need to dump the + # dictionary objects with sorted keys, otherwise we could see + # differences based on the arbitrary order in which keys are added. + if hasattr(self.parameter_ranges, '__iter__'): + parameter_ranges_value = json.dumps(self.parameter_ranges, sort_keys=True) + else: + parameter_ranges_value = None + + if hasattr(self.initial_parameters, '__iter__'): + initial_parameters_value = json.dumps(self.initial_parameters, sort_keys=True) + else: + initial_parameters_value = None + + # TODO: It might be good to include parameter constraint functions + # here too, but finding a consistent hash for a function like that + # would be very difficult... + + # Also, we have to cast the bool variables to actual bools, since + # for some reason they might get transformed into numpy versions + # of True and False (numpy.True_ and numpy.False_) when being read + # from an hdf5 file. + + # The actual object we will hash is a tuple of the important + # identifying information + hash_data = retyped_data + (bool(self.fixed_xmin), + xmin_value, + xmax_value, + bool(self.discrete), + self.fit_method, + bool(self.estimate_discrete) if self.discrete else False, + self.discrete_normalization if self.discrete else '', + self.xmin_distribution_cls.name, + self.xmin_distance, + bool(self.test_all_xmin), + parameter_ranges_value, + initial_parameters_value) + + # Python's basic hash() function isn't consistent across different + # runs, so we can't use it to identify the specific properties of a + # Fit object. In contrast, hashlib's functions are consistent. + + return int(hashlib.sha256(str(hash_data).encode("utf-8")).hexdigest(), 16) + + + @staticmethod + def set_cache_folder(path): + """ + """ + global _cache_enabled, _cache_path + + # First, create the folder if it doesn't exist + os.makedirs(path, exist_ok=True) + + # Save the path and set caching enabled + _cache_enabled = True + _cache_path = os.path.abspath(path) + + + @staticmethod + def set_cache_format(format): + """ + """ + global _cache_format + if format in SUPPORTED_SAVE_FILE_EXTENSIONS: + _cache_format = format + + else: + raise ValueError(f'Invalid save file format provided: {format}. Available formats are: {SUPPORTED_SAVE_FILE_EXTENSIONS}') + + +def _write_hdf5_file(filename, fit): + """ + Save a fit object to a file, using the hdf5 file format. + + Parameters + ---------- + filename : str or Path + The path to the file that will be created. + + fit : Fit + The fit object. + """ + with h5py.File(filename, 'w') as f: + # Create the main dataset + dataset = f.create_dataset('data', data=fit.data_original) + + # h5 files can't save Python's NoneType (since they should + # work with any language) so we have to do a little parsing. + metadata = {} + metadata["xmin"] = fit.xmin + metadata["fixed_xmin"] = fit.fixed_xmin + metadata["xmax"] = fit.xmax if fit.xmax else np.nan + + metadata["discrete"] = fit.discrete + metadata["fit_method"] = fit.fit_method + metadata["estimate_discrete"] = fit.estimate_discrete if fit.discrete else False + metadata["discrete_normalization"] = fit.discrete_normalization + metadata["sigma_threshold"] = fit.sigma_threshold if fit.sigma_threshold else np.nan + metadata["xmin_distribution_name"] = fit.xmin_distribution_cls.name + metadata["xmin_distance"] = fit.xmin_distance + metadata["test_all_xmin"] = fit.test_all_xmin + metadata["noise_flag"] = getattr(fit, 'noise_flag', False) + + # For the initial parameters and parameter ranges, we're going + # to have to unpack each entry in those dictionaries (since we + # can't store a proper dictionary). + if fit.initial_parameters: + for k, v in fit.initial_parameters.items(): + metadata[f"initial_{k}"] = v if v else np.nan + + if fit.parameter_ranges: + for k, v in fit.parameter_ranges.items(): + metadata[f"range_{k}_min"] = v[0] if v[0] else np.nan + metadata[f"range_{k}_max"] = v[1] if v[1] else np.nan + + # The parameter constraints are quite tricky as well, since we + # can't store a Python function object in an h5 file. Instead, we + # just copy the source code for the constraint. Note that this + # does mean you can't use values defined outside of your constraint + # function, for example: + # + # some_value = 5 + # def constraint(dist): + # return dist.value < some_value + # + # This above constraint function cannot be saved properly because + # some_value is defined outside of its scope. + if fit.parameter_constraints: + for constr in fit.parameter_constraints: + function_source = inspect.getsource(constr["fun"]) + function_name = constr["fun"].__name__ + + # This serialized function is the thing we will actually + # use to reconstruct the function. It can also be done + # by just executing the source code, though this will + # lose access to any variable defined outside of the function. + # The serialized code will include those. + serialized_function = pickle.dumps(constr["fun"]) + + # But the serialized function isn't easily readable, so we + # want to save the source code too. + function_data = np.array([function_source, serialized_function], dtype='S') + + # We create a new dataset for each constraint, with the + # main data being the function source code + constraint_dataset = f.create_dataset(f'constraint_{function_name}', data=function_data) + + constraint_dataset.attrs.update({'type': constr["type"], + 'dists': constr.get("dists", ''), + 'name': function_name}) + + # Include the xmin fitting results if available. + if hasattr(fit, 'xmin_fitting_results'): + # We'll make a new dataset folder with entries for each of + # the arrays + + # We also save the minimum index of these arrays in each + # one. This might be a little redundant, but again the idea + # of using h5 in the first place is that it is accessible + # outside of this library and outside of even Python. + distances = f.create_dataset('xmin_fitting_results/distances', data=fit.xmin_fitting_results["distances"]) + distances.attrs['min_index'] = fit.xmin_fitting_results["min_index"] + + xmins = f.create_dataset('xmin_fitting_results/xmins', data=fit.xmin_fitting_results["xmins"]) + xmins.attrs['min_index'] = fit.xmin_fitting_results["min_index"] + + valid_fits = f.create_dataset('xmin_fitting_results/valid_fits', data=fit.xmin_fitting_results["valid_fits"]) + valid_fits.attrs['min_index'] = fit.xmin_fitting_results["min_index"] + + # Now we have to create arrays for the parameters of the specific + # xmin distribution we are using + xmin_fitting_parameters = fit.xmin_distribution_cls.parameter_names + for param in xmin_fitting_parameters: + param_dataset = f.create_dataset(f'xmin_fitting_results/{param}', data=fit.xmin_fitting_results[param]) + param_dataset.attrs['min_index'] = fit.xmin_fitting_results["min_index"] + + # Finally, we include the version, since we want to warn if we + # open a file made in another version. + metadata["powerlaw_version"] = POWERLAW_VERSION if POWERLAW_VERSION else 'Unknown' + + dataset.attrs.update(metadata) + + +def _write_pickle_file(filename, fit): + """ + Save a fit object to a file, using the Python pickle format. + + Parameters + ---------- + filename : str or Path + The path to the file that will be created. + + fit : Fit + The fit object. + """ + with open(filename, 'wb') as f: + pickle.dump((fit, POWERLAW_VERSION), f) + + +def _parse_hdf5_file(filename, verbose=1): + """ + Parse an hdf5 file created using the ``powerlaw.Fit.save()`` function. + + Not intended to be called directly, but rather through + ``powerlaw.Fit.load()``. + + Parameters + ---------- + filename : str or Path + The path to the file that will be loaded. + + Returns + ------- + fit : Fit + The loaded fit object, which should be functionally identical + to the saved object. + """ + with h5py.File(filename) as f: + data = f["data"][:] + # Most of the metadata here we will directly pass to the Fit + # object, but we do have to parse a few things. + metadata = dict(f["data"].attrs) + + # Parse the initial parameter values + initial_parameters = {} + for k, v in metadata.items(): + if "initial_" in k: + variable_name = k.split('_')[1] + initial_parameters[variable_name] = v if (not np.isnan(v)) else None + + if len(initial_parameters) == 0: + initial_parameters = None + + # Parse the parameter ranges + # These will be entries of the form "range_{var}_min" or "range_{var}_max" + parameter_ranges = {} + for k, v in metadata.items(): + if "range_" in k: + if k.split('_')[-1] not in ["max", "min"]: + warnings.warn(f'Malformed h5 file formatting for parameter ranges: {k}:{v}') + continue + + variable_name = k.split('_')[1] + range_index = int(k.split('_')[-1] == "max") + + current_range = parameter_ranges.get(variable_name, [None, None]) + current_range[range_index] = float(v) if (not np.isnan(v)) else None + + parameter_ranges[variable_name] = current_range + + if len(parameter_ranges) == 0: + parameter_ranges = None + + # Parse the parameter constraints + parameter_constraints = [] + for dataset in f.keys(): + if "constraint_" in dataset: + function_source = f[dataset][:][0] + serialized_function = f[dataset][:][1] + + # Execute the source code to load the function in + try: + # Note that we have to "dedent" this function since it + # will maintain any indentation from exactly where + # it was written. This gets rid of extra indentation + # and puts the "def ..." line at zero indent. + #exec(textwrap.dedent(function_source.decode('utf-8'))) + # It is better to use the serialization to reconstruct + # the function since it can save the values of variables + # defined outside of the function. + function = pickle.loads(serialized_function) + + except: + raise ValueError(f'Malformed constraint function {dataset} in file {filename}.') + + constraint_metadata = dict(f[dataset].attrs) + + function_type = constraint_metadata["type"] + dists = constraint_metadata["dists"] + + constraint_dict = {"type": function_type, + "fun": function} + if len(dists) > 0: + constraint_dict["dists"] = list(dists) + + parameter_constraints.append(constraint_dict) + + + if len(parameter_constraints) == 0: + parameter_constraints = None + + # Parse xmin fitting results. These are stored in separate + # datasets within the folder (within the file) 'xmin_fitting_results'. + if not metadata["fixed_xmin"]: + + xmin_fitting_results = {} + xmin_fitting_results["distances"] = f["xmin_fitting_results/distances"][:] + xmin_fitting_results["valid_fits"] = f["xmin_fitting_results/valid_fits"][:] + xmin_fitting_results["xmins"] = f["xmin_fitting_results/xmins"][:] + + xmin_fitting_results["min_index"] = f["xmin_fitting_results/xmins"].attrs["min_index"] + + # Now get the parameter arrays + xmin_fitting_parameters = SUPPORTED_DISTRIBUTIONS[metadata["xmin_distribution_name"]].parameter_names + for param in xmin_fitting_parameters: + xmin_fitting_results[param] = f[f"xmin_fitting_results/{param}"][:] + + + # There are various issues that might arise if you use a different + # version of the package from the one that created the cached + # file. For example, if there was a bug in some calculation + # that was fixed in a newer version, you might not recalculate + # the value to fix the issue. Most of the time you probably + # shouldn't have issues, but we should give a warning just in + # case. + if POWERLAW_VERSION.lower() != "unknown" and metadata["powerlaw_version"].lower() != "unknown": + if POWERLAW_VERSION != metadata["powerlaw_version"]: + warnings.warn(f'Cached file {filename} was saved with a different version of powerlaw {metadata["powerlaw_version"]} than what you are currently using {POWERLAW_VERSION}! This may cause issues...') + + # Create the fit object + fit = Fit(data=data, + discrete=bool(metadata["discrete"]), + xmin=metadata["xmin"], + xmax=metadata["xmax"] if (not np.isnan(metadata["xmax"])) else None, + fit_method=metadata["fit_method"], + estimate_discrete=bool(metadata["estimate_discrete"] if metadata["discrete"] else None), + discrete_normalization=metadata["discrete_normalization"], + sigma_threshold=metadata["sigma_threshold"] if (not np.isnan(metadata["sigma_threshold"])) else None, + initial_parameters=initial_parameters, + parameter_ranges=parameter_ranges, + parameter_constraints=parameter_constraints, + xmin_distance=metadata["xmin_distance"], + xmin_distribution=metadata["xmin_distribution_name"], + test_all_xmin=bool(metadata["test_all_xmin"]), + ignore_cache=True, # We have to ignore cacheing otherwise we'll get an infinite loop. + verbose=verbose) + + # Now we have to adjust the fact that maybe we did actually do + # xmin fitting (but we just cached the result) + if not metadata["fixed_xmin"]: + fit.fixed_xmin = False + fit.noise_flag = metadata["noise_flag"] + + fit.xmin_fitting_results = xmin_fitting_results + + # Set the Fit's xmin to the optimal xmin + setattr(fit, metadata["xmin_distance"], xmin_fitting_results["distances"][xmin_fitting_results["min_index"]]) + + # Update the fitting CDF given the new xmin, in case other objects, like + # Distributions, want to use it for fitting (like if they do KS fitting) + fit.fitting_cdf_bins, fit.fitting_cdf = fit.cdf() + + return fit + + +def _parse_pickle_file(filename): + """ + Parse a pickle file created using the ``powerlaw.Fit.save()`` function. + + Not intended to be called directly, but rather through + ``powerlaw.Fit.load()``. + + Parameters + ---------- + filename : str or Path + The path to the file that will be loaded. + + Returns + ------- + fit : Fit + The loaded fit object, which should be functionally identical + to the saved object. + """ + with open(filename, 'rb') as f: + # Note that we save the version with the fit object + fit, version = pickle.load(f) + + # There are various issues that might arise if you use a different + # version of the package from the one that created the cached + # file. For example, if there was a bug in some calculation + # that was fixed in a newer version, you might not recalculate + # the value to fix the issue. Most of the time you probably + # shouldn't have issues, but we should give a warning just in + # case. + if POWERLAW_VERSION.lower() != "unknown" and version.lower() != "unknown": + if POWERLAW_VERSION != version: + warnings.warn(f'Cached file {filename} was saved with a different version of powerlaw {version} than what you are currently using {POWERLAW_VERSION}! This may cause issues...') + + return fit diff --git a/pyproject.toml b/pyproject.toml index e9c1a7b..7c1fa10 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -33,7 +33,9 @@ dependencies = [ "numpy", "matplotlib", "mpmath", - "tqdm" + "tqdm", + "h5py", + "dill" ] [project.optional-dependencies] diff --git a/setup.py b/setup.py index 4223e2b..265a18c 100644 --- a/setup.py +++ b/setup.py @@ -22,7 +22,7 @@ author='Jeff Alstott', author_email='jeffalstott@gmail.com', url='http://www.github.com/jeffalstott/powerlaw', - install_requires=['scipy', 'numpy', 'matplotlib', 'mpmath', 'tqdm'], + install_requires=['scipy', 'numpy', 'matplotlib', 'mpmath', 'tqdm', 'h5py', 'dill'], license='MIT', classifiers=[ 'License :: OSI Approved :: MIT License', diff --git a/testing/test_saving_loading.py b/testing/test_saving_loading.py new file mode 100644 index 0000000..d562a41 --- /dev/null +++ b/testing/test_saving_loading.py @@ -0,0 +1,277 @@ +""" +This class tests that you can properly pickle fit objects. + +The old version used the reference datasets, but I think this isn't necessary; +it is totally reasonable just to use synthetic data, since this is cleaner +and requires less overhead (making this file easier to read). +""" +import numpy as np +from numpy.testing import assert_equal, assert_allclose + +import os +import unittest +import powerlaw +import pickle + +class TestHashing(unittest.TestCase): + + def test_hash_equal(self): + """ + """ + # The specific dataset is arbitrary + data = powerlaw.load_test_dataset('blackouts') + + # xmin given in different type (but same value) + fit_1 = powerlaw.Fit(data, xmin=1) + fit_2 = powerlaw.Fit(data, xmin=1.0) + + assert hash(fit_1) == hash(fit_2) + + # estimate_discrete for non discrete distribution + fit_1 = powerlaw.Fit(data, xmin=1) + fit_2 = powerlaw.Fit(data, xmin=1, estimate_discrete=True) + + assert hash(fit_1) == hash(fit_2) + + # Fitting xmin + fit_1 = powerlaw.Fit(data) + fit_2 = powerlaw.Fit(data) + + assert hash(fit_1) == hash(fit_2) + + # Data type + fit_1 = powerlaw.Fit(data.astype(np.float32)) + fit_2 = powerlaw.Fit(data.astype(np.float64)) + + assert hash(fit_1) == hash(fit_2) + + # Rounding error + fit_1 = powerlaw.Fit(data) + fit_2 = powerlaw.Fit(data + 1e-20) + + assert hash(fit_1) == hash(fit_2) + + # Accessing a distribution + fit_1 = powerlaw.Fit(data, discrete=True) + fit_2 = powerlaw.Fit(data, discrete=True) + fit_2.power_law + + assert hash(fit_1) == hash(fit_2) + + + def test_hash_not_equal(self): + """ + """ + # The specific dataset is arbitrary + data = powerlaw.load_test_dataset('fires') + + # Changing the data + fit_1 = powerlaw.Fit(data, xmin=1) + fit_2 = powerlaw.Fit(data + 1e-5, xmin=1) + + assert hash(fit_1) != hash(fit_2) + + # Changing the xmin + fit_1 = powerlaw.Fit(data, xmin=1) + fit_2 = powerlaw.Fit(data, xmin=2) + + assert hash(fit_1) != hash(fit_2) + + # Discrete vs continuous + fit_1 = powerlaw.Fit(data, xmin=1, discrete=False) + fit_2 = powerlaw.Fit(data, xmin=1, discrete=True) + + assert hash(fit_1) != hash(fit_2) + + # xmax + fit_1 = powerlaw.Fit(data, xmin=1, xmax=100) + fit_2 = powerlaw.Fit(data, xmin=1, xmax=None) + + assert hash(fit_1) != hash(fit_2) + + +# The save function will determine what format to use based on the filename +# given, so we can test both h5 and pickle saving just by changing the name +# of the output path. +TEST_FILE_H5 = 'unit_test_saving_loading.h5' +TEST_FILE_PKL = 'unit_test_saving_loading.pkl' + +TEST_FILES = [TEST_FILE_H5, TEST_FILE_PKL] + +class TestSavingLoading_H5(unittest.TestCase): + + def test_compare_hash(self): + """ + This test makes sure that saved and loaded files have the same + hash. This covers pretty much all of the variables in the class, + though there are a few outside of this that we test separately. + """ + # The specific dataset is arbitrary + data = powerlaw.load_test_dataset('blackouts') + + for test_f in TEST_FILES: + # With fixed xmin + try: + original_fit = powerlaw.Fit(data, xmin=1) + original_fit.save(test_f) + loaded_fit = powerlaw.Fit.load(test_f) + + assert hash(original_fit) == hash(loaded_fit) + + finally: + os.remove(test_f) + + # Fitting xmin + try: + original_fit = powerlaw.Fit(data) + original_fit.save(test_f) + loaded_fit = powerlaw.Fit.load(test_f) + + assert hash(original_fit) == hash(loaded_fit) + + finally: + os.remove(test_f) + + + def test_compare_xmin_fitting(self): + """ + The xmin fitting results can't be included in the hash because + we might not know the results at the time we need to hash, so + we test them separately. + """ + # The specific dataset is arbitrary + data = powerlaw.load_test_dataset('blackouts') + + for test_f in TEST_FILES: + try: + original_fit = powerlaw.Fit(data) + original_fit.save(test_f) + loaded_fit = powerlaw.Fit.load(test_f) + + fitting_result_keys = ['distances', 'xmins', 'valid_fits'] + for key in fitting_result_keys: + assert_equal(original_fit.xmin_fitting_results[key], loaded_fit.xmin_fitting_results[key]) + + finally: + os.remove(test_f) + + + def test_compare_fit(self): + """ + This test compares the fitted parameters for a specific distribution + between the original and loaded fits. + + This is quite an important test even if it doesn't seem like it. This + is because the actual distribution objects aren't saved with the + cahched file (for h5 format), so they have to be reconstructed based on the parameters + and data. So if the distribution fitting is exactly the same, that + indicates that indeed everything about the fit is functionally identical. + """ + # The specific dataset is arbitrary + data = powerlaw.load_test_dataset('flares') + + for test_f in TEST_FILES: + try: + original_fit = powerlaw.Fit(data) + original_fit.save(test_f) + loaded_fit = powerlaw.Fit.load(test_f) + + for dist in original_fit.supported_distributions.keys(): + original_values = list(getattr(original_fit, dist).parameters.values()) + loaded_values = list(getattr(loaded_fit, dist).parameters.values()) + + assert_allclose(original_values, loaded_values, rtol=0.01, atol=0.01) + + finally: + os.remove(test_f) + + + def test_compare_fit_initial_parameters(self): + """ + This test compares the fitted parameters for a specific distribution + between the original and loaded fits. + """ + # The specific dataset is arbitrary + data = powerlaw.load_test_dataset('flares') + + for test_f in TEST_FILES: + try: + # The parameters are intentionally scattered (and you would almost + # never given a None value for one, but since this is for unit + # testing we want to push things a bit). + original_fit = powerlaw.Fit(data, initial_parameters={"alpha": 1.2, "Lambda": 1e-2, "mu": None}) + original_fit.save(test_f) + loaded_fit = powerlaw.Fit.load(test_f) + + for dist in original_fit.supported_distributions.keys(): + original_values = list(getattr(original_fit, dist).parameters.values()) + loaded_values = list(getattr(loaded_fit, dist).parameters.values()) + + assert_allclose(original_values, loaded_values, rtol=0.01, atol=0.01) + + finally: + os.remove(test_f) + + + def test_compare_fit_parameter_ranges(self): + """ + This test compares the fitted parameters for a specific distribution + between the original and loaded fits. + """ + # The specific dataset is arbitrary + data = powerlaw.load_test_dataset('flares') + + for test_f in TEST_FILES: + try: + # The parameters are intentionally scattered (and you would almost + # never given a None value for one, but since this is for unit + # testing we want to push things a bit). + original_fit = powerlaw.Fit(data, parameter_ranges={"alpha": [1.2, 1.5], "Lambda": [1e-5, None]}) + original_fit.save(test_f) + loaded_fit = powerlaw.Fit.load(test_f) + + for dist in original_fit.supported_distributions.keys(): + original_values = list(getattr(original_fit, dist).parameters.values()) + loaded_values = list(getattr(loaded_fit, dist).parameters.values()) + + assert_allclose(original_values, loaded_values, rtol=0.01, atol=0.01) + + finally: + os.remove(test_f) + + + def test_compare_fit_parameter_constraints(self): + """ + This test compares the fitted parameters for all distributions with + a constraint. + """ + # The specific dataset is arbitrary + data = powerlaw.load_test_dataset('blackouts') + + for test_f in TEST_FILES: + try: + # As discussed in the documentation, this should be entirely + # self contained. + def constr(dist): + N = 100 + return len(dist.data) - N + + constraint_dict = {"type": 'ineq', + "fun": constr, + "dists": ['power_law']} + + # The parameters are intentionally scattered (and you would almost + # never given a None value for one, but since this is for unit + # testing we want to push things a bit). + original_fit = powerlaw.Fit(data, parameter_constraints=[constraint_dict]) + original_fit.save(test_f) + loaded_fit = powerlaw.Fit.load(test_f) + + for dist in original_fit.supported_distributions.keys(): + original_values = list(getattr(original_fit, dist).parameters.values()) + loaded_values = list(getattr(loaded_fit, dist).parameters.values()) + + assert_allclose(original_values, loaded_values, rtol=0.01, atol=0.01) + + finally: + os.remove(test_f)