# Load Biospytial modules and etc.
%matplotlib inline
import sys
sys.path.append('/apps')
import django
django.setup()
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
## Use the ggplot style
plt.style.use('ggplot')

from external_plugins.spystats import tools
%run ../testvariogram.py

section.shape

"""
Explanation: Here I'm process by chunks the entire region.
End of explanation
"""


minx,maxx,miny,maxy = getExtent(new_data)

maxy

## Let's build the partition
N = 6
xp,dx = np.linspace(minx,maxx,N,retstep=True)
yp,dy = np.linspace(miny,maxy,N,retstep=True)


dx

xx,yy = np.meshgrid(xp,yp)

coordinates_list = [ (xx[i][j],yy[i][j]) for i in range(N) for j in range(N)]

from functools import partial
tuples = map(lambda (x,y) : partial(getExtentFromPoint,x,y,step_sizex=dx,step_sizey=dy)(),coordinates_list)

len(tuples)

chunks = map(lambda (mx,Mx,my,My) : subselectDataFrameByCoordinates(new_data,'newLon','newLat',mx,Mx,my,My),tuples)
chunks_sizes = map(lambda df : df.shape[0],chunks)

chunk_w_size = zip(chunks,chunks_sizes)

## Here we can filter based on a threshold
threshold = 10
nonempty_chunks_w_size = filter(lambda (df,n) : df.shape[0] > threshold ,chunk_w_size)
chunks_non_empty, chunks_sizes = zip(*nonempty_chunks_w_size)

lengths = pd.Series(map(lambda ch : ch.shape[0],chunks_non_empty))

lengths.plot.hist()

cs = chunks_non_empty
variograms =map(lambda chunk : tools.Variogram(chunk,'residuals1',using_distance_threshold=600000),cs)

%time vars = map(lambda v : v.calculateEmpirical(),variograms)
%time vars = map(lambda v : v.calculateEnvelope(num_iterations=1),variograms)

%time lags = map(lambda v : v.lags,variograms)

lags = pd.DataFrame(lags).transpose()

lags = lags[[0]]

"""
Explanation: Algorithm for processing Chunks

Make a partition given the extent
Produce a tuple (minx ,maxx,miny,maxy) for each element on the partition
Calculate the semivariogram for each chunk and save it in a dataframe
Plot Everything
Do the same with a mMatern Kernel
End of explanation
"""


envslow = pd.concat(map(lambda df : df[['envlow']],vars),axis=1)
envhigh = pd.concat(map(lambda df : df[['envhigh']],vars),axis=1)
variogram = pd.concat(map(lambda df : df[['variogram']],vars),axis=1)
n_points = pd.DataFrame(map(lambda v : v.n_points,variograms))

points = n_points.transpose()

print(variogram.shape)
print(points.shape)

ejem1 = pd.DataFrame(variogram.values * points.values)

# Chunks (variograms) columns
# lag rows
vempchunk = ejem1.sum(axis=1) / points.sum(axis=1)
plt.plot(lags,vempchunk,'--',color='blue',lw=2.0)

## Cut some values
vchunk = pd.concat([lags,vempchunk],axis=1)
vchunk.columns = ['lags','semivariance']
v = vchunk[vchunk['lags'] < 500000]
plt.plot(v.lags,v.semivariance,'--',color='blue',lw=2.0)
#vemp2



"""
Explanation: Take an average of the empirical variograms also with the envelope.
We will use the group by directive on the field lags
End of explanation
"""


thrs_dist = 1000000
nt = 30 # num iterations
filename = "../HEC_runs/results/low_q/data_envelope.csv"
envelope_data = pd.read_csv(filename)
gvg = tools.Variogram(new_data,'logBiomass',using_distance_threshold=thrs_dist)
gvg.envelope = envelope_data
gvg.empirical = gvg.envelope.variogram
gvg.lags = gvg.envelope.lags
vdata = gvg.envelope.dropna()
gvg.plot(refresh=False,legend=False,percentage_trunked=20)
plt.plot(v.lags,v.semivariance,'--',color='blue',lw=2.0)

"""
Explanation: Let's bring the whole empirical variogram (calculated in HEC)
For comparison purposes
End of explanation
"""


cs = chunks_non_empty
variograms2 =map(lambda chunk : tools.Variogram(chunk,'residuals2',using_distance_threshold=600000),cs)
%time vars = map(lambda v : v.calculateEmpirical(),variograms2)
%time vars = map(lambda v : v.calculateEnvelope(num_iterations=1),variograms2)
%time lags = map(lambda v : v.lags,variograms2)
lags = pd.DataFrame(lags).transpose()
envslow = pd.concat(map(lambda df : df[['envlow']],vars),axis=1)
envhigh = pd.concat(map(lambda df : df[['envhigh']],vars),axis=1)
variogram2 = pd.concat(map(lambda df : df[['variogram']],vars),axis=1)


n_points = pd.DataFrame(map(lambda v : v.n_points,variograms2))
points = n_points.transpose()
ejem2 = pd.DataFrame(variogram2.values * points.values)
# Chunks (variograms) columns
# lag rows
vempchunk2 = ejem2.sum(axis=1) / points.sum(axis=1)
#plt.plot(lags,vempchunk2,'--',color='blue',lw=2.0)
thrs_dist = 1000000
nt = 30 # num iterations
filename = "../HEC_runs/results/low_q/data_envelope.csv"
envelope_data = pd.read_csv(filename)
gvg = tools.Variogram(new_data,'logBiomass',using_distance_threshold=thrs_dist)
gvg.envelope = envelope_data
gvg.empirical = vempchunk2
gvg.lags = gvg.envelope.lags
vdata = gvg.envelope.dropna()

gvg.plot(refresh=False,legend=False,percentage_trunked=20)
plt.plot(variograms2[0].lags,vempchunk2)
#plt.plot(v.lags,v.semivariance,'--',color='blue',lw=2.0)

## Save the emp variogram
#variograms2[0].lags
outputvar = pd.DataFrame([variograms2[0].lags,vempchunk2]).transpose()
outputvar.columns = ['lags','semivariance']
outputvar.to_csv("/outputs/logbiomas_logsppn_residuals_variogram.csv")


## Here I´ll plot the thing
plt.plot(vdata.lags,vdata.variogram)


"""
Explanation: Ok, now the thing that I was supposed to do since the begining
The log of the species instead of only species number
End of explanation
"""


## test duplicates in new_data
new_data[new_data.duplicated(['newLon'], keep=False)]

## Remove duplications
withoutrep = new_data.drop_duplicates(subset=['newLon','newLat'])
print(new_data).shape
print(withoutrep).shape
new_data = withoutrep


def gaussianVariogram(h,sill=0,range_a=0,nugget=0):
    if isinstance(h,np.ndarray):
        Ih = np.array([1.0 if hx >= 0.0 else 0.0 for hx in h])
    else:
        Ih = 1.0 if h >= 0 else 0.0
    #Ih = 1.0 if h >= 0 else 0.0    
    g_h = ((sill - nugget)*(1 - np.exp(-(h**2 / range_a**2)))) + nugget*Ih
    return g_h
## Fitting model.
### Optimizing the empirical values
def theoreticalVariogram(model_function,sill,range_a,nugget):
    return lambda x : model_function(x,sill,range_a,nugget)

dt = pd.concat([lags[[0]],vempchunk2],axis=1)
dt.columns = ['lags','semivariance']
vdata = dt.dropna()

import matplotlib.pylab as plt
hx = np.linspace(0,600000,100)
from scipy.optimize import curve_fit
s = 0.345
r = 100000.0
nugget = 0.33
init_vals = [0.34, 50000, 0.33]     # for [amp, cen, wid]
bg, covar_gaussian = curve_fit(gaussianVariogram, xdata=vdata.lags.values, ydata=vdata.semivariance.values, p0=init_vals)
#vdata = gvg.envelope.dropna()
## The best parameters asre:
gau_var = gaussianVariogram(hx,bg[0],bg[1],bg[2])
## The best is the gaussian
#model_var = theoreticalVariogram(gaussianVariogram,bg[0],bg[1],bg[2])

#gvg.empirical = vdata.
#empirical variogram logbiomas ~ sppn
gvg.plot(refresh=False,legend=False)

#empirical variogram by chunks logbiomass ~ logsppn 
#plt.plot(dt.lags,dt.semivariance,'--',color='blue',lw=2.0)
#model
plt.plot(hx,gau_var,color='green',lw=2.3)
#plt.plot(hx,exp_var)


"""
Explanation: The variograms for logbiomas and log sppn are different from logbiomass and log sppn.
Being the logSppn lower spatial variance. Indeed the sill reaches lower variance. 
The chunk method provides a faster way of calculating the empirical variogram albeit a slightly more uncertainty (noise)
when compared with the global variogram.
To do now.
1 .Revisit the model for generating the covariance matrix
2. The problem was that the matrix was not positive semidefinite. There could be repetitions in the original data
End of explanation
"""


print("Sill %s, Range %s, Nugget %s"%(bg[0],bg[1],bg[2]))

"""
Explanation: The best fitted parameters are:
Sill, range , nugget (respectively)
End of explanation
"""


import statsmodels.regression.linear_model as lm
%time Mdist = gvg.distance_coordinates.flatten()
#%time vars = np.array(map(model_var,Mdist))
#CovMat = vars.reshape(len(section),len(section))
#X = section.logSppN.values
#Y = section.logBiomass.values

# But the data is massive and it consumes all my RAM (32 GB) I need to  do something clever.
covar = []
for i in range(len(Mdist)):
    x = Mdist.pop(i)
    g = gaussianVariogram(x,bg[0],bg[1],bg[2])
    covar.append(g)

### Apply the distance function, rememeber that now it receives a list (array) as X
%time vars = gaussianVariogram(Mdist,bg[0],bg[1],bg[2])

values = lags[[0]].values

for i in range(len(values)):
    print(values[i] - values[i+1])

"""
Explanation: Now for finishing I need to perform the lucky step,
Re fit the $\beta$ 
Oke, first calculate the distances
End of explanation
"""


