Verifying pyfssa
================

The *pyfssa* package implements algorithmic finite-size scaling
analysis. In order to verify the autoscaling algorithm, the following
computations aim at reproducing the scaling exponents from the
literature of the paradigmatic bond percolation problem on a linear
chain and on a square lattice. First, finite-size percolation studies
utilizing the *pypercolate* package generate the finite-size data. Next,
we plot these raw data for the percolation probability, the order
parameter (largest component) and the susceptibility (second moment of
the cluster-size distribution). For each of these quantities, we scale
the simulation data according to the literature exponents and observe
the data collapse onto the master curve. Finally, for each of these
quantities, *pyfssa* auto-scales the simulation data in order to recover
the literature exponents.

Percolation simulation design
-----------------------------

We utilize the
`pypercolate <http://pypercolate.readthedocs.org/en/stable/>`__ package
to generate the percolation finite-size data, according to the `template
for a HPC bond percolation finite-size scaling
study <http://pypercolate.readthedocs.org/en/stable/pypercolate-hpc.html#template-for-a-hpc-bond-percolation-finite-size-scaling-study>`__.
To carry out the simulations in a GridEngine powered high-performance
computing environment, we employ the
`jug <http://luispedro.org/software/jug/>`__ and
`gridjug <http://gridjug.readthedocs.org/en/stable/>`__ packages. The
*jugfiles* contain the control flow for both studies on the `linear
chain <_static/bond_percolation_1d_jugfile.py>`__ and the `square
lattice <_static/bond_percolation_jugfile.py>`__. The finite system
sizes of the linear chain (number of nodes) are
:math:`2^6, 2^8, \ldots, 2^{20}`, and the linear extensions :math:`L_i`
of the :math:`L_i \times L_i` square lattices (number of nodes) are
:math:`2^3, 2^4, \ldots, 2^{10}`. The total number of runs is
:math:`N_{runs} = 10^4`, respectively. The occupation probabilities are
:math:`\{1 - 10^{-(1 + i/10)} : i = 0, \ldots, 100\}` for the linear
chain and :math:`\{0.5 + \frac{i}{500} : i = -50, \ldots, +50\}` for the
square lattice. The confidence level is at the :math:`1 \sigma` level,
hence :math:`1 - \alpha = \text{erf}\frac{1}{\sqrt{2}} \approx 0.683`.
The master seeds that initialize the Mersenne-Twister RNG to draw the
seeds for each run are :math:`3138595788` for the linear chain and
:math:`3939666632` for the square lattice.

Preamble
--------

.. code-block:: python

    from __future__ import division

.. code-block:: python

    # configure plotting
    %config InlineBackend.rc = {'figure.dpi': 300, 'savefig.dpi': 300, \
                                'figure.figsize': (3, 3 / 1.6), 'font.size': 10, \
                                'figure.facecolor': (1, 1, 1, 0)}
    %matplotlib inline

.. code-block:: python

    from pprint import pprint as print
    from collections import OrderedDict
    import warnings
    
    import h5py
    import matplotlib as mpl
    import numpy as np
    import os
    import scipy.stats
    import seaborn as sns
    
    import matplotlib.pyplot as plt
    
    from cycler import cycler
    
    import fssa


.. parsed-literal::

    /home/sorge/repos/sci/pyfssa/.devenv35/lib/python3.5/site-packages/matplotlib/__init__.py:872: UserWarning: axes.color_cycle is deprecated and replaced with axes.prop_cycle; please use the latter.
      warnings.warn(self.msg_depr % (key, alt_key))


Loading and inspecting the raw data
-----------------------------------

.. code-block:: python

    bond_percolation_data = OrderedDict()
    bond_percolation_data['1D'] = h5py.File(
        os.path.join('_static', 'bond-percolation-hpc-chain-canonical-averages.hdf5')
    )
    bond_percolation_data['2D'] = h5py.File(
        os.path.join('_static', 'bond-percolation-hpc-grid-canonical-averages.hdf5')
    )

.. code-block:: python

    system_sizes = OrderedDict()
    for key, data in bond_percolation_data.items():
        system_sizes[key] = sorted([int(group_key) for group_key in data.keys()])
    print(system_sizes)


.. parsed-literal::

    OrderedDict([('1D', [64, 256, 1024, 4096, 16384, 65536, 262144, 1048576]),
                 ('2D', [8, 16, 32, 64, 128, 256, 512, 1024])])


.. code-block:: python

    print(bond_percolation_data['1D'][str(system_sizes['1D'][0])].dtype.descr)


.. parsed-literal::

    [('number_of_runs', '<u4'),
     ('p', '<f8'),
     ('alpha', '<f8'),
     ('percolation_probability_mean', '<f8'),
     ('percolation_probability_std', '<f8'),
     ('percolation_probability_ci', '<f8', (2,)),
     ('percolation_strength_mean', '<f8'),
     ('percolation_strength_std', '<f8'),
     ('percolation_strength_ci', '<f8', (2,)),
     ('moments_mean', '<f8', (5,)),
     ('moments_std', '<f8', (5,)),
     ('moments_ci', '<f8', (5, 2))]


.. code-block:: python

    def ci_to_std(a, ci):
        return np.max(
            np.abs(
                np.transpose(ci)
                - np.tile(a, (2, 1))
            ),
            axis=0
        )

.. code-block:: python

    ps = dict()
    n = dict()
    alpha = dict()
    for key, my_system_sizes in system_sizes.items():
        ps[key] = bond_percolation_data[key][str(my_system_sizes[0])]['p']
        n[key] = bond_percolation_data[key][str(my_system_sizes[0])]['number_of_runs']
        alpha[key] = bond_percolation_data[key][str(my_system_sizes[0])]['alpha']
        assert np.all(n[key] == n[key][0])
        assert np.all(alpha[key] == alpha[key][0])
        for system_size in my_system_sizes:
            assert np.all(bond_percolation_data[key][str(system_size)]['p'] == ps[key])
            assert np.all(bond_percolation_data[key][str(system_size)]['number_of_runs'] == n[key])
            assert np.all(bond_percolation_data[key][str(system_size)]['alpha'] == alpha[key])
        print((key, ps[key], n[key][0], alpha[key][0]))


.. parsed-literal::

    ('1D',
     array([ 0.9       ,  0.91087491,  0.92056718,  0.92920542,  0.93690427,
            0.94376587,  0.94988128,  0.95533164,  0.96018928,  0.96451866,
            0.96837722,  0.97181617,  0.97488114,  0.97761279,  0.98004738,
            0.98221721,  0.98415107,  0.98587462,  0.98741075,  0.98877982,
            0.99      ,  0.99108749,  0.99205672,  0.99292054,  0.99369043,
            0.99437659,  0.99498813,  0.99553316,  0.99601893,  0.99645187,
            0.99683772,  0.99718162,  0.99748811,  0.99776128,  0.99800474,
            0.99822172,  0.99841511,  0.99858746,  0.99874107,  0.99887798,
            0.999     ,  0.99910875,  0.99920567,  0.99929205,  0.99936904,
            0.99943766,  0.99949881,  0.99955332,  0.99960189,  0.99964519,
            0.99968377,  0.99971816,  0.99974881,  0.99977613,  0.99980047,
            0.99982217,  0.99984151,  0.99985875,  0.99987411,  0.9998878 ,
            0.9999    ,  0.99991087,  0.99992057,  0.99992921,  0.9999369 ,
            0.99994377,  0.99994988,  0.99995533,  0.99996019,  0.99996452,
            0.99996838,  0.99997182,  0.99997488,  0.99997761,  0.99998005,
            0.99998222,  0.99998415,  0.99998587,  0.99998741,  0.99998878,
            0.99999   ,  0.99999109,  0.99999206,  0.99999292,  0.99999369,
            0.99999438,  0.99999499,  0.99999553,  0.99999602,  0.99999645,
            0.99999684,  0.99999718,  0.99999749,  0.99999776,  0.999998  ,
            0.99999822,  0.99999842,  0.99999859,  0.99999874,  0.99999888,
            0.999999  ]),
     10000,
     0.31731050786291415)
    ('2D',
     array([ 0.49  ,  0.4902,  0.4904,  0.4906,  0.4908,  0.491 ,  0.4912,
            0.4914,  0.4916,  0.4918,  0.492 ,  0.4922,  0.4924,  0.4926,
            0.4928,  0.493 ,  0.4932,  0.4934,  0.4936,  0.4938,  0.494 ,
            0.4942,  0.4944,  0.4946,  0.4948,  0.495 ,  0.4952,  0.4954,
            0.4956,  0.4958,  0.496 ,  0.4962,  0.4964,  0.4966,  0.4968,
            0.497 ,  0.4972,  0.4974,  0.4976,  0.4978,  0.498 ,  0.4982,
            0.4984,  0.4986,  0.4988,  0.499 ,  0.4992,  0.4994,  0.4996,
            0.4998,  0.5   ,  0.5002,  0.5004,  0.5006,  0.5008,  0.501 ,
            0.5012,  0.5014,  0.5016,  0.5018,  0.502 ,  0.5022,  0.5024,
            0.5026,  0.5028,  0.503 ,  0.5032,  0.5034,  0.5036,  0.5038,
            0.504 ,  0.5042,  0.5044,  0.5046,  0.5048,  0.505 ,  0.5052,
            0.5054,  0.5056,  0.5058,  0.506 ,  0.5062,  0.5064,  0.5066,
            0.5068,  0.507 ,  0.5072,  0.5074,  0.5076,  0.5078,  0.508 ,
            0.5082,  0.5084,  0.5086,  0.5088,  0.509 ,  0.5092,  0.5094,
            0.5096,  0.5098,  0.51  ]),
     10000,
     0.31731050786291415)


Scaling and visualization routines
----------------------------------

.. code-block:: python

    def autoscale(system_sizes, ps, a, da, rho_c, nu, zeta):
        
        # optimize scaling quality;
        # initialize with exponents from literature
        with warnings.catch_warnings():
            warnings.simplefilter('ignore')
            res = fssa.autoscale(
                l=system_sizes,
                rho=ps,
                a=a,
                da=da,
                rho_c0=rho_c,
                nu0=nu,
                zeta0=zeta,
            )
        
        if not res['success']:
            raise RuntimeError
        
        return {
            res_key: res[res_key]
            for res_key in [
                'rho', 'drho', 'nu', 'dnu', 'zeta', 'dzeta'
            ]
        }

Visualizing the raw data and the algorithmic scaling
----------------------------------------------------

Configuration
~~~~~~~~~~~~~

.. code-block:: python

    config = {
        '1D': {
            # 'percolation_probability': {
            #     'slice': np.s_[:],
            #     'exponents': {
            #         'rho_c': 1.0,
            #         'nu': 1.0,
            #         'zeta': 0.0,
            #     },
            # },
            'quantities': [
                'percolation_strength', 'moments',
            ],
            'percolation_strength': {
                'autoscale_slice': np.s_[4:],
                'slice': np.s_[:],
                'exponents': {
                    'rho_c': 1.0,
                    'nu': 1.0,
                    'zeta': 0.0,
                },
            },
            'moments': {
                'autoscale_slice': np.s_[4:],
                'slice': np.s_[:, 2],
                'exponents': {
                    'rho_c': 1.0,
                    'nu': 1.0,
                    'zeta': 1.0,
                },
                'log_data': True,
            },
        },
        '2D': {
            # 'percolation_probability': {
            #     'slice': np.s_[:],
            #     'exponents': {
            #         'rho_c': 0.5,
            #         'nu': 1.333,
            #         'zeta': 0.0,
            #     },            
            # }, 
            'quantities': [
                'percolation_strength', 'moments',
            ],        
            'percolation_strength': {
                'autoscale_slice': np.s_[4:],
                'slice': np.s_[:],
                'exponents': {
                    'rho_c': 0.5,
                    'nu': 1.333,
                    'zeta': -0.133,
                },
            },
            'moments': {
                'autoscale_slice': np.s_[4:],
                'slice': np.s_[:, 2],
                'exponents': {
                    'rho_c': 0.5,
                    'nu': 1.333,
                    'zeta': 2.4,
                },
            },        
        },
        # 'percolation_probability': {
        #     'ylabel': '$P$',            
        # },
        'percolation_strength': {
            'ylabel': '$C_1/N$',            
        },    
        'moments': {
            'ylabel': '$\chi$',
        }
    }

Preparing the data
~~~~~~~~~~~~~~~~~~

.. code-block:: python

    a = dict()
    da = dict()
    all_exponents = dict()
    
    for key, my_system_sizes in system_sizes.items():
        a[key] = dict()
        da[key] = dict()
        all_exponents[key] = dict()
        
        for quantity in config[key]['quantities']:
            my_config = config[key][quantity]
            all_exponents[key][quantity] = dict()
            my_slice = my_config['slice']
            exponents = my_config['exponents']
            all_exponents[key][quantity]['literature'] = exponents
            mean_field = '{}_mean'.format(quantity)
            ci_field = '{}_ci'.format(quantity)
            
            a[key][quantity] = np.empty((len(my_system_sizes), ps[key].size))
            da[key][quantity] = np.empty((len(my_system_sizes), ps[key].size))
    
            for i, system_size in enumerate(my_system_sizes):
                ds = bond_percolation_data[key][str(system_size)]
                a[key][quantity][i] = ds[mean_field][my_slice]
                da[key][quantity][i] = ci_to_std(a[key][quantity][i], ds[ci_field][my_slice])

Auto-scaling the data
~~~~~~~~~~~~~~~~~~~~~

.. code-block:: python

    for key, my_system_sizes in system_sizes.items():  # '1D', '2D'
        for quantity in config[key]['quantities']:
            my_config = config[key][quantity]
            my_slice = my_config['autoscale_slice']
            all_exponents[key][quantity]['autoscaled'] = autoscale(
                system_sizes=my_system_sizes[my_slice],
                ps=ps[key],
                a=a[key][quantity][my_slice],
                da=da[key][quantity][my_slice],
                **my_config['exponents']
            )

.. code-block:: python

    # scale data with literature and optimized exponents
    scaled_data = dict()
    for key, my_system_sizes in system_sizes.items():  # '1D', '2D'
        scaled_data[key] = dict()
        for quantity in config[key]['quantities']:
            my_config = config[key][quantity]
            scaled_data[key][quantity] = dict()
            
            opt_exponents = {
                'rho_c': all_exponents[key][quantity]['autoscaled']['rho'],
                'nu': all_exponents[key][quantity]['autoscaled']['nu'],
                'zeta': all_exponents[key][quantity]['autoscaled']['zeta'],            
            }
            for exp_key, exponents in zip(
                ['literature', 'optimized'],
                [my_config['exponents'], opt_exponents]
            ):
                with warnings.catch_warnings():
                    warnings.simplefilter('ignore')
    
                    scaled_data[key][quantity][exp_key] = fssa.scaledata(
                        l=my_system_sizes,
                        rho=ps[key],
                        a=a[key][quantity],
                        da=da[key][quantity],
                        **exponents
                    )

Plotting
~~~~~~~~

.. code-block:: python

    assert len(system_sizes['1D']) == len(system_sizes['2D'])
    num_system_sizes = len(system_sizes['1D'])
    
    assert len(config['1D']['quantities']) == len(config['2D']['quantities'])
    num_quantities = len(config['1D']['quantities'])
    
    # Define colors
    palette = sns.cubehelix_palette(
        n_colors=num_system_sizes,
        start=2.0, rot=0.4, gamma=1.0, hue=1.0, light=0.8, dark=0.2,
    )
    
    plot_properties = {
        'ms': 2,
        'fmt': '.',
    }
    
    for key, my_system_sizes in system_sizes.items():  # '1D', '2D'
        fig, axes = plt.subplots(
            nrows=num_quantities, ncols=3,
            figsize=(6, 1 + 2 * num_quantities),
        )
    
        # plot original data
        axes[0, 0].set_title('Simulation data')
        for quantity, my_axes in zip(config[key]['quantities'], axes):
            my_config = config[key][quantity]
            
            ax = my_axes[0]
            ax.set_prop_cycle(cycler('color', palette))
            if my_config.get('log_data', False):
                ax.set_xscale('log')
                ax.set_yscale('log')
                for y, dy, size in zip(
                    a[key][quantity], da[key][quantity], my_system_sizes,
                ):
                    ax.errorbar(
                        1 - ps[key], y, yerr=dy, label=size,
                        **plot_properties
                    )
                ax.invert_xaxis()
                ax.locator_params(axis='x', numticks=3)        
                ax.locator_params(axis='y', numticks=4)        
                ax.set_xlabel(r'$1 - p$')
    
            else:
                # default
                for y, dy in zip(a[key][quantity], da[key][quantity]):
                    ax.errorbar(
                        ps[key], y, yerr=dy,
                        **plot_properties
                    )
                ax.locator_params(axis='x', nbins=3)
                ax.locator_params(axis='y', nbins=4)
                ax.set_xlabel(r'$p$')
    
            try:
                ax.set_ylabel(config[quantity]['ylabel'])
            except KeyError:
                pass
            
        # plot data scaled by literature and optimized exponents
        for exp_index, (exp_key, title) in enumerate(zip(
            ['literature', 'optimized'],
            ['Scaled (lit.)', 'Auto-scaled'],
        )):
            axes[0, exp_index + 1].set_title(title)
            
            for quantity_index, quantity in enumerate(
                config[key]['quantities']
            ):
                my_config = config[key][quantity]
                my_scaled_data = scaled_data[key][quantity][exp_key]
            
                ax = axes[quantity_index, exp_index + 1]
                ax.set_prop_cycle(cycler('color', palette))
                if key == '1D':
                    ax.set_xscale('log')
                    for x, y, dy in zip(
                        my_scaled_data.x,
                        my_scaled_data.y,
                        my_scaled_data.dy,
                    ):
                        ax.errorbar(
                            -x, y, yerr=dy, **plot_properties
                        )
                    ax.invert_xaxis()
                    ax.locator_params(axis='x', numticks=3)        
                    ax.set_xlabel(r'$|x|$')
                else:
                    # default
                    for x, y, dy in zip(
                        my_scaled_data.x,
                        my_scaled_data.y,
                        my_scaled_data.dy,
                    ):
                        ax.errorbar(
                            x, y, yerr=dy, **plot_properties
                        )
                    ax.locator_params(axis='x', nbins=3)
                    ax.set_xlabel(r'$x$')
        
                ax.locator_params(axis='y', nbins=4)
        
        import matplotlib.lines as mlines
        
        fig.legend(
            handles=[
                mlines.Line2D(
                    [], [], ls=' ', marker='.', color=color, ms=3,
                )
                for color in palette
            ],
            labels=[
                r'$2^{{{:n}}}$'.format(log)
                for log in np.log2(my_system_sizes)
            ],
            numpoints=3,
            loc='lower center',
            ncol=num_system_sizes,
            columnspacing=0.75,
            handletextpad=0.25,
            title=r'Linear dimension $L$',
            frameon=True,
            handlelength=1,
        )
        plt.tight_layout(
            rect=[0, 0.1, 1, 0.9]
        )
        plt.show()



.. image:: verification_files/verification_25_0.png



.. image:: verification_files/verification_25_1.png


.. code-block:: python

    for data in bond_percolation_data.values():
        data.close()

Results
-------

.. code-block:: python

    print(all_exponents)


.. parsed-literal::

    {'1D': {'moments': {'autoscaled': {'dnu': 4.7167237931903998e-06,
                                       'drho': 1.4479762551576001e-06,
                                       'dzeta': 0.0076361696412922571,
                                       'nu': 1.0034367452207009,
                                       'rho': 1.0000001610135159,
                                       'zeta': 1.011552722196388},
                        'literature': {'nu': 1.0, 'rho_c': 1.0, 'zeta': 1.0}},
            'percolation_strength': {'autoscaled': {'dnu': 0.0066467778306519773,
                                                    'drho': 1.1767283504220902e-07,
                                                    'dzeta': 0.00018273505453284011,
                                                    'nu': 1.0034333511814242,
                                                    'rho': 1.0000000280887937,
                                                    'zeta': 0.00022191407753831049},
                                     'literature': {'nu': 1.0,
                                                    'rho_c': 1.0,
                                                    'zeta': 0.0}}},
     '2D': {'moments': {'autoscaled': {'dnu': 0.029676064038296639,
                                       'drho': 9.0370392368313164e-05,
                                       'dzeta': 0.048353931195945513,
                                       'nu': 1.3319317711769838,
                                       'rho': 0.50006987574569961,
                                       'zeta': 2.3795003879213787},
                        'literature': {'nu': 1.333, 'rho_c': 0.5, 'zeta': 2.4}},
            'percolation_strength': {'autoscaled': {'dnu': 0.020699500197520947,
                                                    'drho': 0.00014571284024314797,
                                                    'dzeta': 0.011702576022673383,
                                                    'nu': 1.3264821248629133,
                                                    'rho': 0.49998937452554693,
                                                    'zeta': -0.13984359978778532},
                                     'literature': {'nu': 1.333,
                                                    'rho_c': 0.5,
                                                    'zeta': -0.133}}}}


Discussion
----------

For the linear chain, *pyfssa* recovers the critical point both for the
percolation strength and the susceptibility, where the analytical value
:math:`p_c = 1.0` lies within the negligibly small error intervals of
sizes :math:`10^{-6}` and :math:`10^{-7}`. The percolation strength
scaling further recovers the :math:`\nu` scaling exponent within the
small error interval of order of magnitude :math:`10^{-3}`. While the
interval for the scaling exponent :math:`\nu` derived from automatically
scaling the susceptibility data does not contain the literature value
:math:`1.0`, it deviates only by a small fraction of :math:`10^{-3}`.
The intervals for both the order parameter exponent :math:`\beta` and
the susceptibility exponent :math:`\gamma` do not contain the true
values. However, the algorithmically determined values for :math:`\beta`
and :math:`\gamma` deviate from the true values by only a small fraction
of the order of :math:`10^{-4}` and :math:`10^{-2}`, respectively.

For the square lattice, the intervals for all exponents and all
magnitudes contain the literature value and are rather small:
:math:`10^{-2}` for the scaling exponents :math:`\nu`, :math:`\beta` and
:math:`\gamma`; :math:`10^{-4}` for the critical point :math:`p_c`.

In the light of this formidable agreement, we need to keep in mind that
we initialized the auto-scaling algorithm with the literature values.
This certainly helps making the Nelder-Mead optimization results and
error estimates well-defined.

Conclusion
----------

Both visual inspection as well as algorithmic scaling analysis of the
bond percolation data reproduce well-established results from the
literature with formidable agreement, particulary in the two-dimensional
setting. Even though the finite-size scaling analyst needs to take the
usual care of manually (visually) approaching data collapse in order to
initialize the Nelder-Mead search, we conclude the foregoing
computational analysis to confirm that *pyfssa* is correctly
implemented.