Basic sums as summary characteristics
=====================================

Example 1: Visualization of dynamic changes of the system
---------------------------------------------------------

.. code:: ipython3

    import numpy as np

Let us import example data.

.. code:: ipython3

    from basicsums.examples import load_example03
    A_hex, r_hex, A_sqr, r_sqr = load_example03()

``A_hex`` and ``A_sqr`` contain 128 samples of 256 equal disks of
concentration 0.7 each. ``A_hex`` and ``r_hex`` represents arrays of the
corresponding radii.

.. code:: ipython3

    np.pi*sum(r_hex[0]**2), len(A_hex), len(A_hex[0])




.. parsed-literal::

    (0.7000000000000015, 128, 256)



First samples of each set represent regular arrays of disks: hexagonal
and square:

.. code:: ipython3

    from basicsums.show_disks import show_disks

.. code:: ipython3

    w1_hex = np.sqrt(2.) / (3**(1/4.))
    w2_hex = w1_hex * np.exp(np.pi/3*1j)
    show_disks(A_hex[0], r_hex[0], w1_hex, w2_hex)



.. image:: tutorial_geometric_parameters_files/tutorial_geometric_parameters_9_0.png


.. code:: ipython3

    w1_sqr = 1
    w2_sqr = 1j
    show_disks(A_sqr[0], r_sqr[0], w1_sqr, w2_sqr)



.. image:: tutorial_geometric_parameters_files/tutorial_geometric_parameters_10_0.png


The remaining samples were created as an consequence of the random walk
of the initial poisitons. Let us visualize the system after 10 steps:

.. code:: ipython3

    t = 10
    show_disks(A_hex[t], r_hex[t], w1_hex, w2_hex)
    show_disks(A_sqr[t], r_sqr[t], w1_sqr, w2_sqr)



.. image:: tutorial_geometric_parameters_files/tutorial_geometric_parameters_12_0.png



.. image:: tutorial_geometric_parameters_files/tutorial_geometric_parameters_12_1.png


Hence, we have dynamicaly changing systems. They could model the stirr casting process [#KURT2013]_. Intuitively, as time of random walk aproaches infinity, both systems should converge to the same distribution (or very similar ones). Since basic sums describes geometry of such systems, the convergence should reflects in the values of certain sums. It is known that  in both initial regular cases :math:`e_{4, 4} = S_4^2`, where :math:`S_4` is a :ref:`Eisenstein-Rayleigh lattice sum <sect-theory_lattice_sums>`. Morover, :math:`S_4` takes different values in both cases. Let us compute :math:`e_{4,4}` for each sample:

.. code:: ipython3

    from basicsums.basic_sums import Cell, BasicSums

.. code:: ipython3

    qmax = 4
    cell_hex = Cell(w1_hex, w2_hex, qmax)
    cell_sqr = Cell(w1_sqr, w2_sqr, qmax)

You can compute the results yourself (it can take some time):

.. code:: ipython3

    res_hex = [BasicSums(A, cell_hex, eisenstein_indexes=[4]).esum(4, 4)
               for A in A_hex]

.. code:: ipython3

    res_sqr = [BasicSums(A, cell_sqr, eisenstein_indexes=[4]).esum(4, 4)
               for A in A_sqr]

or load them directly:

.. code:: ipython3

    from basicsums.examples import load_example03e44
    res_hex, res_sqr = load_example03e44()

.. hint::
   One can use the ``tqdm`` package to track the progress of calculations as follows:

.. code:: ipython3

    from tqdm import tqdm
    res_hex = [BasicSums(A, cell_hex, eisenstein_indexes=[4]).esum(4, 4)
               for A in tqdm(A_hex)]   

Let us create a plot of the results. One can see that the proces of changes of the random structures reflects in the real values of :math:`e_{4,4}`:

.. code:: ipython3

    import matplotlib.pyplot as plt
    plt.plot(np.real(res_hex), '.', np.real(res_sqr), '.')
    plt.grid(True)
    plt.legend(['hexagonal', 'square'])
    plt.show()



.. image:: tutorial_geometric_parameters_files/tutorial_geometric_parameters_24_0.png


Example 2: Summary characteristics of distributions of disks
------------------------------------------------------------

Let us load another example

.. code:: ipython3

    from basicsums.examples import load_example04
    A_norm, r_norm, A_unif, r_unif = load_example04()

``A_norm`` and ``A_unif`` contain 10 samples of 256 equal disks of
concentration 0.5 each. Arrays ``r_norm`` and ``r_unif`` represents
arrays of the corresponding radii.

.. code:: ipython3

    np.pi*sum(r_norm[0]**2), len(A_norm), len(A_norm[0])




.. parsed-literal::

    (0.5000000000000002, 10, 256)



.. code:: ipython3

    w1, w2 = 1, 1j
    show_disks(A_norm[0], r_norm[0], w1, w2)
    show_disks(A_unif[0], r_unif[0], w1, w2)



.. image:: tutorial_geometric_parameters_files/tutorial_geometric_parameters_30_0.png



.. image:: tutorial_geometric_parameters_files/tutorial_geometric_parameters_30_1.png


One can see, that the former samples represent disks of normally
distributed radii, whereas the latter have radii drawn from the uniform
distribution:

.. code:: ipython3

    import matplotlib.pyplot as plt
    
    fig, (ax1, ax2) = plt.subplots(2, sharex=True)
    ax1.hist(r_norm.flatten(), 20)
    ax2.hist(r_unif.flatten(), 20)
    
    ax1.set_ylim([0, 350])
    ax2.set_ylim([0, 350])
    ax1.grid(True)
    ax2.grid(True)
    
    plt.show()



.. image:: tutorial_geometric_parameters_files/tutorial_geometric_parameters_32_0.png


hence, we face two different distributions of disks. In order to characterize them, we need some kind of summary characteristics. Let us apply the following set of basic sums :math:`\{(2), (2, 2), (3, 3), (4, 4), \ldots, (9, 9)\}` as a set of summary characteristics of distributions of disks.

.. code:: ipython3

    qmax = 9
    cell = Cell(w1, w2, qmax)
    sums_indexes = [(2,)] + [(p, p) for p in range(2, qmax+1)]

Let us compute considered basic sums for each sample

.. code:: ipython3

    res_norm = np.array([BasicSums(A, cell, r).esums(sums_indexes)
                         for A, r in zip(A_norm, r_norm)])

.. code:: ipython3

    res_unif = np.array([BasicSums(A, cell, r).esums(sums_indexes)
                         for A, r in zip(A_unif, r_unif)])

The following code creates plots of real parts of each sum:

.. code:: ipython3

    fig, ax = plt.subplots(figsize=(10, 20))
    for k, tup in enumerate(sums_indexes):
        plt.subplot(5, 2, k+1)
        plt.plot(np.real(res_norm[:,k]), 'o',
                 np.real(res_unif[:,k]), 'o')
        plt.title('$e_{' + str(tup)[1:-1] + '}$')
        plt.grid(True)
    plt.show()



.. image:: tutorial_geometric_parameters_files/tutorial_geometric_parameters_39_0.png


One can observe, that the larger order of the sum is, the greater difference between distributions it reveals. It is also known from the theory of composites that :math:`e_2 =\pi` for macroscopically isotropic composites. Hence, one can see that both our distributions are rather isotropic.

It is known from theory of composites that :math:`e_2 =\pi` for macroscopically isotropic composites. Hence, one can see that both our distributions are isotropic.

.. _sect-characteristics_points:

Example 3: Summary characteristics point patterns and regularization
--------------------------------------------------------------------

One can consider centres of disks from the above example as point patterns, neglecting their radii.

.. code:: ipython3

    k = 2
    fig, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize=(9, 4))
    ax1.plot(np.real(A_norm[k]), np.imag(A_norm[k]), '.')
    ax2.plot(np.real(A_unif[k]), np.imag(A_unif[k]), '.')
    ax1.set_aspect('equal')
    ax2.set_aspect('equal')
    plt.show()



.. image:: tutorial_geometric_parameters_files/tutorial_geometric_parameters_45_0.png


In order to apply basic sums as summary characteristics of point patterns, let us first present the simplest approach, i.e. to treat the points as centres of identical disks:

.. code:: ipython3

    res_norm_points = np.array([BasicSums(A, cell).esums(sums_indexes)
                                  for A in A_norm])

.. code:: ipython3

    res_unif_points = np.array([BasicSums(A, cell).esums(sums_indexes)
                                  for A in A_unif])

.. code:: ipython3

    fig, ax = plt.subplots(figsize=(10, 20))
    for k, tup in enumerate(sums_indexes):
        plt.subplot(5, 2, k+1)
        plt.plot(np.real(res_norm_points[:,k]), 'o',
                 np.real(res_unif_points[:,k]), 'o')
        plt.title('$e_{' + str(tup)[1:-1] + '}$')
        plt.grid(True)
    plt.show()



.. image:: tutorial_geometric_parameters_files/tutorial_geometric_parameters_49_0.png


One can see, that :math:`e_{2,2}` distinguishes both class of patterns very well. One can oslo observe large values of sums and deviation of the values, especiali for ``A_unif``. This situation is common in pattern, where a very close points appear. Such a pair of points yields the difference close to zero, hence the values of :math:`E_n` can be very large, since each of the functions :math:`E_n` has a pole of order :math:`n` at :math:`z=0`.

.. code:: ipython3

    np.std(res_norm_points, axis=0)




.. parsed-literal::

    array([9.54853920e-02, 3.06705218e-01, 1.13712103e+00, 3.96363087e+00,
           2.69872549e+01, 8.88881247e+01, 4.15380244e+02, 2.67842762e+03,
           9.78122244e+03])



.. code:: ipython3

    np.std(res_unif_points, axis=0)




.. parsed-literal::

    array([1.60139980e-01, 1.17750500e+00, 4.50692448e+00, 3.32374173e+01,
           3.44128094e+02, 3.35543656e+03, 2.97933172e+04, 2.98269574e+05,
           2.96228716e+06])



In order to obtain more stable characteristics, one can follow another schemme, where the pattern is considered as a marked one. in order to to this, we apply the *regularization* process, in which a set of radii (marks) is generated.
For each point :math:`a_k\in A` we calculate :math:`r_k=\frac12\min_{j\neq k}d(a_k,a_j)` where :math:`d` is a distance in the torus topology. One can consider :math:`r_k` as the *nearest neighbour marks*.

Module :py:mod:`basicsums.preparation` contains the function :py:func:`regularized_radii` computing the marks:

.. code:: ipython3

    from basicsums.preparation import regularized_radii

Let us compute regularized radii for each sample:

.. code:: ipython3

    radii_norm = [regularized_radii(A, w1, w2) for A in A_norm]
    radii_unif = [regularized_radii(A, w1, w2) for A in A_unif]

In this approach, the pattern is treated as a distribution of disks of different radii:

.. code:: ipython3

    show_disks(A_unif[0], radii_unif[0], w1, w2)



.. image:: tutorial_geometric_parameters_files/tutorial_geometric_parameters_59_0.png


Hence, we can calculate basic sums for each pattern:

.. code:: ipython3

    res_norm_points_reg = np.array([BasicSums(A, cell, r).esums(sums_indexes)
                             for A, r in zip(A_norm, radii_norm)])
    res_unif_points_reg = np.array([BasicSums(A, cell, r).esums(sums_indexes)
                             for A, r in zip(A_unif, radii_unif)])

This time, the deviations of the results are much smaller:

.. code:: ipython3

    np.std(res_norm_points_reg, axis=0)




.. parsed-literal::

    array([0.04815164, 0.22962557, 0.27581293, 0.40329808, 0.66573284,
           0.89955487, 0.87183017, 1.64641344, 2.9807516 ])



.. code:: ipython3

    np.std(res_unif_points_reg, axis=0)




.. parsed-literal::

    array([0.06233756, 0.32575137, 0.28875951, 0.60124032, 0.72054763,
           1.4864379 , 2.69567961, 4.80540731, 8.41615877])



.. code:: ipython3

    fig, ax = plt.subplots(figsize=(10, 20))
    for k, tup in enumerate(sums_indexes):
        plt.subplot(5, 2, k+1)
        plt.plot(np.real(res_norm_points_reg[:,k]), 'o',
                 np.real(res_unif_points_reg[:,k]), 'o')
        plt.title('$e_{' + str(tup)[1:-1] + '}$')
        plt.grid(True)
    plt.show()



.. image:: tutorial_geometric_parameters_files/tutorial_geometric_parameters_65_0.png


Poisson
~~~~~~~

.. code:: ipython3

    from basicsums.show_disks import show_disks

.. code:: ipython3

    np.random.seed(10)
    patterns = [np.random.rand(256, 2).dot([1, 1j]) - (0.5 + 0.5j)
                for i in range(10)]

.. code:: ipython3

    show_disks(patterns[0], 0.005*np.ones(256), w1, w2)



.. image:: tutorial_geometric_parameters_files/tutorial_geometric_parameters_69_0.png


.. code:: ipython3

    res_poisson = [BasicSums(A, cell).esums(sums_indexes)
               for A in patterns]

.. code:: ipython3

    res_poisson = np.array(res_poisson)

.. code:: ipython3

    res_poisson.mean(axis=0)




.. parsed-literal::

    array([ 7.43803510e+01-3.32786887e+01j, -7.63239815e+11-7.80213066e-05j,
            7.57274786e+16+1.16441814e+00j, -7.51395614e+21-3.48134672e+03j,
            7.45563025e+26+2.41893146e+10j, -7.39775737e+31-1.44660313e+14j,
            7.34033372e+36-7.23848718e+19j, -7.28335582e+41-2.67437777e+22j])



.. code:: ipython3

    res_poisson.std(axis=0)




.. parsed-literal::

    array([2.37197330e+02, 2.28962306e+12, 2.27182197e+17, 2.25418678e+22,
           2.23668907e+27, 2.21932721e+32, 2.20210012e+37, 2.18500674e+42])



.. code:: ipython3

    import matplotlib.pyplot as plt
    k = 1
    plt.plot(np.real(res_poisson[:,k]), 'o')
    plt.grid(True)
    plt.show()



.. image:: tutorial_geometric_parameters_files/tutorial_geometric_parameters_74_0.png


.. code:: ipython3

    from basicsums.preparation import regularized_radii

.. code:: ipython3

    radii = [regularized_radii(A, w1, w2) for A in patterns]

.. code:: ipython3

    show_disks(patterns[0], radii[0], w1, w2)



.. image:: tutorial_geometric_parameters_files/tutorial_geometric_parameters_77_0.png



.. code:: ipython3

    res_poisson_reg = [BasicSums(A, cell, r).esums(sums_indexes)
               for A, r in zip(patterns, radii)]

.. code:: ipython3

    res_poisson_reg = np.array(res_poisson_reg)

.. code:: ipython3

    res_poisson_reg.mean(axis=0)




.. parsed-literal::

    array([ 3.09909679e+00+6.97105417e-02j, -2.24978693e+01-3.77635105e-17j,
            6.68322118e+01-4.52428521e-16j, -2.01918248e+02-3.22721685e-16j,
            6.07481859e+02+1.44567732e-14j, -1.81014006e+03+2.01664255e-14j,
            5.63754484e+03+1.53328129e-14j, -1.70439305e+04-4.07862527e-13j])



.. code:: ipython3

    res_poisson_reg.std(axis=0)




.. parsed-literal::

    array([3.21408956e-01, 4.24840338e+00, 1.50153967e+01, 5.80584374e+01,
           2.29771452e+02, 8.68856282e+02, 2.99968303e+03, 1.05832972e+04])



.. code:: ipython3

    import matplotlib.pyplot as plt
    k = -1
    plt.plot(np.real(res_poisson_reg[:,k]), 'o')
    plt.grid(True)
    plt.show()



.. image:: tutorial_geometric_parameters_files/tutorial_geometric_parameters_83_0.png


.. rubric:: References


.. [#KURT2013] P. Kurtyka, N. Rylko, *Structure analysis of the modified cast metal matrix composites by use of the RVE theory*. Arch Metall. Mater 58(2), 357-360, 2013. doi: 10.2478/v10172-012-0198-x.