%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np

import open_cp.retrohotspot
import open_cp.data
import open_cp.predictors
import open_cp.evaluation

n = 1000
times = np.random.random(n) * 365
times = times * (np.timedelta64(1, "D") / np.timedelta64(1, "s")) * np.timedelta64(1, "s")
times = np.sort(np.datetime64("2016-01-01") + times)
points = open_cp.data.TimedPoints(times, np.random.random((2,n)) * 1000)

mask = [[False]*100]*100
grid = open_cp.data.MaskedGrid(10, 10, 0, 0, mask)

predictor = open_cp.retrohotspot.RetroHotSpot()
predictor.data = points
predictor.weight = open_cp.retrohotspot.TruncatedGaussian(100, 10)
cts_pred = predictor.predict(end_time=np.datetime64("2016-10-01"))
cts_pred.samples = 5
pred = open_cp.predictors.GridPredictionArray.from_continuous_prediction_grid(cts_pred, grid)

fig, ax = plt.subplots(ncols=2, figsize=(14,8))

ax[0].pcolor(*pred.mesh_data(), pred.intensity_matrix, cmap="Greys")
p = open_cp.evaluation.top_slice_prediction(pred, 0.1)
ax[1].pcolor(*p.mesh_data(), p.intensity_matrix, cmap="Greys")

"""
Explanation: Problems with Monte Carlo integration
By default, we convert "continuous" predictions into grid based predictions by using monte carlo integration: evaluate the kernel at random points and average.
This seems to introduce a large amount of error...  The following demonstrates this with artificial, random data.  However, I noticed this problem on data from Chicago.
End of explanation
"""


cts_pred = predictor.predict(end_time=np.datetime64("2016-10-01"))
cts_pred.samples = 5
pred = open_cp.predictors.GridPredictionArray.from_continuous_prediction_grid(cts_pred, grid)
p1 = open_cp.evaluation.top_slice_prediction(pred, 0.1)

np.sum(p.intensity_matrix.mask ^ p.intensity_matrix.mask)

np.sum(p1.intensity_matrix.mask ^ p1.intensity_matrix.mask)

np.sum(p.intensity_matrix.mask ^ p1.intensity_matrix.mask)

"""
Explanation: Repeat: Form the prediction again, top slice the top 10% of risk again, and compare the selected cells (same as looking at the "mask" of the top sliced prediction.  They differ by a lot!)
End of explanation
"""


cts_pred = predictor.predict(end_time=np.datetime64("2016-10-01"))
cts_pred.samples = 25
pred = open_cp.predictors.GridPredictionArray.from_continuous_prediction_grid(cts_pred, grid)
p = open_cp.evaluation.top_slice_prediction(pred, 0.1)

cts_pred = predictor.predict(end_time=np.datetime64("2016-10-01"))
cts_pred.samples = 25
pred = open_cp.predictors.GridPredictionArray.from_continuous_prediction_grid(cts_pred, grid)
p1 = open_cp.evaluation.top_slice_prediction(pred, 0.1)

np.sum(p.intensity_matrix.mask ^ p1.intensity_matrix.mask)

"""
Explanation: You might argue that 5 samples is too small...
End of explanation
"""


cts_pred = predictor.predict(end_time=np.datetime64("2016-10-01"))
cts_pred.samples = -5
pred = open_cp.predictors.GridPredictionArray.from_continuous_prediction_grid(cts_pred, grid)
p = open_cp.evaluation.top_slice_prediction(pred, 0.1)

cts_pred = predictor.predict(end_time=np.datetime64("2016-10-01"))
cts_pred.samples = -5
pred = open_cp.predictors.GridPredictionArray.from_continuous_prediction_grid(cts_pred, grid)
p1 = open_cp.evaluation.top_slice_prediction(pred, 0.1)

np.sum(p.intensity_matrix.mask ^ p1.intensity_matrix.mask)

"""
Explanation: We have added a "reproducible" mode, whereby we don't sample at random points, but rather on a regular sub-grid
End of explanation
"""


