import sys
import warnings
sys.path.append('../..')

from zephyr.backend import SplineGridInterpolator
from zephyr.middleware import FullwvDatastore

"""
Explanation: Test Interpolators
Brendan Smithyman | March, 2016
End of explanation
"""


fds = FullwvDatastore('xhlayr')
systemConfig = fds.systemConfig
scu = {
    'scale': 1.5,
}
systemConfig.update(scu)

"""
Explanation: Setup
Read a configuration from xhlayr.ini and xhlayr.vp, and set up a scale factor of 1.5; i.e., $sdx = 1.5\times dx$.
End of explanation
"""


sgi = SplineGridInterpolator(systemConfig)
print(sgi.shape)

"""
Explanation: Create a SplineGridInterpolator, which behaves like an $8911\times 20000$ matrix.
End of explanation
"""


with warnings.catch_warnings():
    warnings.simplefilter('ignore')
    
    c = systemConfig['c']
    cn = (sgi * c).reshape((sgi.snz, sgi.snx))
    cr = (sgi.T * cn).reshape((sgi.nz, sgi.nx))

"""
Explanation: Get the velocity c from the systemConfig, downsample it to create cn, and then apply the transpose of the operation to get cr.
End of explanation
"""


%pylab inline

with warnings.catch_warnings():
    warnings.simplefilter('ignore')
    

    plotOpts = {
        'vmin': c.min(),
        'vmax': c.max(),
    }

    subplot(1,3,1)
    imshow(c, **plotOpts)
    colorbar(shrink=0.5)
    title('Original')

    subplot(1,3,2)
    imshow(cn, **plotOpts)
    colorbar(shrink=0.5)
    title('Rescaled')

    subplot(1,3,3)
    imshow(cr, **plotOpts)
    colorbar(shrink=0.5)
    title('After adjoint')

    tight_layout()

nLoops = 100

ca = c
sgiT = sgi.T
for i in xrange(nLoops):
    cb = sgi*ca
    ca = sgiT*cb
ca.shape = (sgi.nz,sgi.nx)

plotOpts = {
    'vmin': c.min(),
    'vmax': c.max(),
}

subplot(1,2,1)
imshow(c, **plotOpts)
colorbar(shrink=0.5)
title('Original')

subplot(1,2,2)
imshow(ca, **plotOpts)
colorbar(shrink=0.5)
title('After %d loops'%nLoops)

"""
Explanation: Results
End of explanation
"""


