MPSlib: Variable template size in mps_snesim_tree and mps_snesim_list

mps_snesim_tree and mps_snesim_list support a template size that changes for each multiple grid level. The template size is specified for the coarsest multiple grid and for the final dense simulation grid; linear interpolation is used at intermediate levels.

For example, using a template size of 8x7x4 on the coarsest grid and 4x3x3 on the finest:

Search template size X # 8 4
Search template size Y # 7 3
Search template size Z # 4 3

This can be set in Python using:

O = mps.mpslib(method='mps_snesim_tree')
O.par['template_size'] = np.array([[8, 7, 4], [4, 3, 3]]).T

A simple constant template of size [8, 7, 4] can be set using:

O.par['template_size'] = np.array([8, 7, 4])

The main reason for using a variable template size is that a considerable amount of CPU is spent in the finer simulation grids pruning (removing) conditional data. Reducing the template at the finest grid can yield significant speedups.

See more at https://mpslib.readthedocs.io/en/latest/Examples/ex_varying_template.html

import mpslib as mps
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec

Setup

#TI1, TI_filename1 = mps.trainingimages.strebelle(2, coarse3d=1)
TI1, TI_filename1 = mps.trainingimages.strebelle(1, coarse3d=1)

O1 = mps.mpslib(method='mps_snesim_tree')
O1.ti = TI1
O1.par['n_multiple_grids'] = 4
O1.par['n_cond'] = 81
O1.par['n_real'] = 1
O1.par['rseed'] = 1
O1.par['debug_level'] = -1
O1.par['simulation_grid_size'][0] = 135
O1.par['simulation_grid_size'][1] = 100
O1.par['simulation_grid_size'][2] = 1

Sweep over template sizes at the finest grid

Fix the coarse-grid template size at 11, and vary the fine-grid template size from 11 down to 1.

r1 = 11  # Template size on the coarsest grid
r2 = [11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1]  # Template size on the finest grid

t = []
R = []
for ir in range(len(r2)):

    O1.delete_local_files()

    template = np.array([[r1, r2[ir]], [r1, r2[ir]], [1, 1]])
    O1.par['template_size'] = template
    name = '%s_%d_%d' % (O1.method, r1, r2[ir])
    print(name)
    O1.parameter_filename = name + '.par'
    O1.mps_snesim_par_write()
    O1.run()
    R.append(O1.sim[0])
    t.append(O1.time)

mps_snesim_tree_11_11
mps_snesim_tree_11_10
mps_snesim_tree_11_9
mps_snesim_tree_11_8
mps_snesim_tree_11_7
mps_snesim_tree_11_6
mps_snesim_tree_11_5
mps_snesim_tree_11_4
mps_snesim_tree_11_3
mps_snesim_tree_11_2
mps_snesim_tree_11_1

Plot realizations and computation time

fig = plt.figure(figsize=(15, 15))
outer = gridspec.GridSpec(4, 3, wspace=0.2, hspace=0.2)

for ir in range(len(r2)):
    ax1 = plt.Subplot(fig, outer[ir])
    fig.add_subplot(ax1)
    plt.imshow(np.transpose(np.squeeze(R[ir])))
    plt.title('template=[%d,%d], t=%g s' % (r1, r2[ir], t[ir]))

ax1 = plt.Subplot(fig, outer[-1])
fig.add_subplot(ax1)
plt.bar(r2, t)
plt.xlabel('Template size at fine grid')
plt.ylabel('Computation time (s)')
plt.savefig('varying_template', dpi=600, facecolor='w', edgecolor='w',
            orientation='portrait', transparent=True)
plt.show()

Speedup relative to full template

fig = plt.figure(figsize=(5, 5))
plt.bar(r2, t[0] / np.array(t))
plt.xlabel('Template size at fine grid')
plt.ylabel('Speedup compared to using full template')
plt.grid()
plt.savefig('varying_template_speedup', dpi=600, facecolor='w', edgecolor='w',
            orientation='portrait', transparent=True)
plt.show()