Saving and loading for Fit objects

.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
+.vscode

README.rst

@@ -1,8 +1,8 @@
 powerlaw: A Python Package for Analysis of Heavy-Tailed Distributions
 =====================================================================
 
-.. image:: https://github.com/jeffalstott/powerlaw/workflows/Tests/badge.svg
-   :target: https://github.com/jeffalstott/powerlaw/actions
+.. image:: https://github.com/powerlaw-devs/powerlaw/workflows/Tests/badge.svg
+   :target: https://github.com/powerlaw-devs/powerlaw/actions
    :alt: Tests
 
 ``powerlaw`` is a toolbox implementing the statistical methods developed in
@@ -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 <https://doi.org/10.1371/journal.pone.0085777>`_
-
-Also available at `arXiv:1305.0215 [physics.data-an] <http://arxiv.org/abs/1305.0215>`_
+(also available at `arXiv:1305.0215 [physics.data-an] <http://arxiv.org/abs/1305.0215>`_)
 
 
 Basic Usage
@@ -73,15 +72,14 @@ 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
@@ -117,16 +115,3 @@ discussions. Thanks also to Andreas Klaus,
 and `Adam Ginsburg <https://github.com/keflavich/plfit>`_ for making
 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 <https://en.wikipedia.org/wiki/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 <http://arxiv.org/abs/0706.1062>`__, which developed the statistical methods that ``powerlaw`` implements.
-
-However, for many data sets, the superior lognormal fit is only possible if one allows the fitted parameter ``mu`` to go negative. Whether or not this is sensible depends on your theory of what's generating the data. If the data is thought to be generated by multiplying random positive variables, ``mu`` is just the log of the distribution's median; a negative ``mu`` just indicates those variables' products are typically below 1. However, if the data is thought to be generated by exponentiating a normal distribution, then ``mu`` is interpreted as the median of the underlying normal data. In that case, the normal data is likely generated by summing random variables (positive and negative), and ``mu`` is those sums' median (and mean). A negative ``mu``, then, indicates that the random variables are typically negative. For some physical systems, this is perfectly possible. For the data you're studying, though, it may be a weird assumption. For starters, all of the data points you're fitting to are positive by definition, since power laws must have positive values (indeed, ``powerlaw`` throws out 0s or negative values). Why would those data be generated by a process that sums and exponentiates *negative* variables?
-
-If you think that your physical system could be modeled by summing and exponentiating random variables, but you think that those random variables should be positive, one possible hacks is ``powerlaw``'s ``lognormal_positive``. This is just a regular lognormal distribution, except ``mu`` must be positive. Note that this does not force the underlying normal distribution to be the sum of only positive variables; it only forces the sums' *average* to be positive, but it's a start. You can compare a power law to this distribution in the normal way shown above::
-
-    R, p = results.distribution_compare('power_law', 'lognormal_positive')
-
-You may find that a lognormal where ``mu`` must be positive gives a much worse fit to your data, and that leaves the power law looking like the best explanation of the data. Before concluding that the data is in fact power law distributed, consider carefully whether a more likely explanation is that the data was generated by multiplying positive random variables, or even by summing and exponentiating random variables; either one would allow for a lognormal with an intelligible negative value of ``mu``.

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 <https://github.com/jeffalstott/powerlaw>`_ for more
+See the `powerlaw home page <https://github.com/powerlaw-devs/powerlaw>`_ for more
 information and examples.
 
 Basic usage

docs/source/tutorials/saving_fits.rst

@@ -0,0 +1,182 @@
+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.
+
+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)
+        ...

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

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):