# Load Biospytial modules and etc.
%matplotlib inline
import sys
sys.path.append('/apps')
sys.path.append('/apps/external_plugins/spystats/')
import django
django.setup()

## Use the ggplot style
plt.style.use('ggplot')


import pandas as pd
from statsmodels.regression import linear_model
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
import numpy as np

import tools
import HEC_runs.fit_fia_logbiomass_logspp_GLS as auxiliary

## File locations, change accordingly
empirical_data_path = "/RawDataCSV/idiv_share/FIA_Plots_Biomass_11092017.csv"
variogram_path = "/apps/external_plugins/spystats/HEC_runs/results/variogram/data_envelope.csv"


"""
Explanation: Linear Model Fitting using empirical variogram
In classical geoestatistics it is frequent to use a sample-based variogram to fit a parametrised theoretical variogram. This method has proved to present biasses for the estimators of range, sill and shape. There are two reasons for this.
1. We are supposing that the observations are not random and the variogram is the true population variogram. 
2. The variogram is built on summary statistics (bin size) and the curve-fitness can vary depending on this size. 
Having said so, it is still a method used for estimating parameters and when the size of the data is big, being a good alternative to the likelihood methodology, a more robust optimization but computationally more complex. 
So first thing first, import the neessary modules. Depending on your machine this may change.
For the moment I work with this.
End of explanation
"""


minx = -130
maxx = -60
miny = 24
maxy = 50

"""
Explanation: Geographic extent
End of explanation
"""


data = auxiliary.initAnalysis(empirical_data_path=empirical_data_path,plotdata_path=plotdata_path, maxx=maxx,minx=minx,maxy=maxy,miny=miny)

"""
Explanation: Auxiliary functions
* Open file
* Convert it to Pandas with Geometry (GeoPandas)
* Reproject it to Lambert Conic Conformal
* Calculate OLS and append residuals
End of explanation
"""


test_data = data.loc[0:49,:]

"""
Explanation: For testing, subselect the first 50 data points
End of explanation
"""


model = tools.WhittleVariogram(sill=0.340246718396,range_a=41188.0234423,nugget=0.329937603763,alpha=1.12143687914)

"""
Explanation: Now lets instantiate the variogram according to the results obtained in the cluster.
Remember:
    < Whittle Variogram : sill 0.340246718396, range 41188.0234423, nugget 0.329937603763, alpha1.12143687914 >
End of explanation
"""


import numpy as np
X = np.linspace(0,1000000,1000)
plt.plot(X,model.corr_f(X))

vg = tools.Variogram(test_data,'residuals',model=model)
MM = vg.calculateCovarianceMatrix()

"""
Explanation: Let's see the correlation
End of explanation
"""


lmod = linear_model.GLS.from_formula(formula='logBiomass ~ logSppN',data=test_data,sigma=MM)

results = lmod.fit()

results.summary()

MM

DD = vg.distance_coordinates

correl_mat = vg.model.corr_f(DD)

correl_mat.tofile('/outputs/correl_matf50.csv',sep=",",format="%s")

MM.tofile('/outputs/covar_matf50.csv',sep=",",format="%s")

ls /outputs

MM.tofile?


lmodco = linear_model.GLS.from_formula(formula='logBiomass ~ logSppN',data=test_data,sigma=correl_mat)

res=lmodco.fit()

res.summary()

len(test_data)

"""
Explanation: Fit a GLS
End of explanation
"""


