%matplotlib inline

from spectral_cube import SpectralCube
cube = SpectralCube.read("/Users/eric/Dropbox/AstroStatistics/Design4_21_0_0_flatrho_0021_13co.fits")

from turbustat.statistics import VCS

vcs_clean = VCS(cube).run(verbose=True)

print([(slop, err) for slop, err in zip(vcs_clean.slopes, vcs_clean.slope_errs)])

"""
Explanation: Does the mask effect the VCS break point?
Our thorough and on-going exploration of the structure function for oddly shaped masks has demonstrated that the slope is dominated by the mask on most of the relevant scales. Since the VCS is inherently linked to this slope, it suggests that it too should show this sensitivity. Thus the slopes will definitely be affected, and should return $m\approx1.5$, where $m$ is the structure function slope, but will it also cause the break point to change? This break point is the resolution limit of the data (Lazarian & Pogosyan, 2006). Even without the inclusion of noise in simulation, we see a change in the position of the break.
End of explanation
"""


masked_cube = cube.with_mask(cube > cube.mean())

p.imshow(masked_cube.moment0().value, origin='lower', cmap='afmhot')

vcs_masked = VCS(masked_cube).run(verbose=True)

print(vcs_masked.brk_err)

print([(slop, err) for slop, err in zip(vcs_masked.slopes, vcs_masked.slope_errs)])

"""
Explanation: The fit gives slopes of $-1.91\pm0.02$ and $-7.54\pm0.04$, and a break at ~5.1 pixels. The difference between the slopes is $-4/m$, so we find $m=0.7$.
Now what happens when I start masking?
End of explanation
"""


masked_cube_70 = cube.with_mask(cube > cube.percentile(70))
p.imshow(masked_cube_70.moment0().value, origin='lower', cmap='afmhot')

vcs_masked_70 = VCS(masked_cube_70).run(verbose=True)

print([(slop, err) for slop, err in zip(vcs_masked_70.slopes, vcs_masked_70.slope_errs)])

"""
Explanation: The difference in the slopes has decreased from $-5.62\pm0.04$ to $-5.11\pm0.04$. The break point has not moved significantly. We now have $m=0.78$.
Increasing the mask threshold:
End of explanation
"""


masked_cube_93 = cube.with_mask(cube > cube.percentile(93))
p.imshow(masked_cube_93.moment0().value, origin='lower', cmap='afmhot')

vcs_masked_93 = VCS(masked_cube_93).run(verbose=True, breaks=-0.9)

print([(slop, err) for slop, err in zip(vcs_masked_93.slopes, vcs_masked_93.slope_errs)])

"""
Explanation: The fit has converged remarkably close to the unmasked data. Hmm...
Masking at the 93th percentile? This completely removes some spectra.
End of explanation
"""


from astropy.wcs import WCS
downgraded_wcs = cube.wcs.copy()
downgraded_wcs.wcs.cdelt[2] *= 2.
data = cube.filled_data[:].value
data[np.isnan(data)] = 0.0
downgraded_cube = SpectralCube(data=nd.zoom(data, (0.5, 1., 1.))*cube.unit, wcs=downgraded_wcs)
downgraded_cube

vcs_downgrade = VCS(downgraded_cube).run(verbose=True, breaks=None)

print([(slop, err) for slop, err in zip(vcs_downgrade.slopes, vcs_downgrade.slope_errs)])

"""
Explanation: $m=1.3$. As soon as those masked regions show up, the break point becomes particularly indistinct. This causes the slopes to converge, the $m=3/2$ limit. The increased power in the noise is having the same effect as smoothing over a larger region.
Does this hold true if I actually smooth/downsample the data?
End of explanation
"""


p.loglog(vcs_clean.vel_freqs, vcs_clean.ps1D, 'bD', alpha=0.05)
p.loglog(vcs_masked.vel_freqs, vcs_masked.ps1D, 'go', alpha=0.05)
p.loglog(vcs_masked_93.vel_freqs, vcs_masked_93.ps1D, 'ro', alpha=0.05)
p.loglog(vcs_downgrade.vel_freqs, vcs_downgrade.ps1D, 'cD', alpha=0.05)

"""
Explanation: This would give $m=2.4$, which doesn't make any sense. I think this has reached the unresolved regime and the asymptotic solution no longer holds.
End of explanation
"""


