#!/usr/bin/env python
# coding: utf-8

# In[1]:


get_ipython().run_line_magic('matplotlib', 'inline')
# Loading python modules used in this notebook
#import pygments
#from pygments import highlight
#from pygments.lexers import MakefileLexer, CppLexer
#from pygments.formatters import HtmlFormatter
from IPython.core.display import display,HTML

import sys
print(sys.version, sys.path)

import sqlite3
import xml.etree.ElementTree as ET
import ruamel.yaml as yaml
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from IPython.display import Image
from matplotlib.ticker import FormatStrFormatter
import os
import fileinput
import re
import subprocess
from displayNC import displayVariables,temperatureDisplay, getSingleArray
from netCDF4 import Dataset
from scipy import stats
import datetime
from plotly import __version__
from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot
import plotly.graph_objs as go
import itertools
import pandas as pd

init_notebook_mode(connected=True)

plt.style.use('seaborn-whitegrid')

    
def plot_dotty(param,perf,names,figname) :
    num_rows = len(names)
    fig, axes = plt.subplots(nrows = num_rows, ncols = 1, figsize = (10, num_rows * 2.5))
    cnt = 0
    formatter = matplotlib.ticker.ScalarFormatter(useOffset=False)
    for name,aH in zip(names, axes) :
        plt.sca(aH)
        lH, = plt.plot(param[:,cnt],perf,'b.',markersize = 1)
        aH.set_xlabel(name)
        #print name
        #plt.xlim = (-100000,-70000)  #This is hardcoded and is horrible. Will eventually fail
        aH.yaxis.set_major_formatter(formatter)
        aH.xaxis.set_major_formatter(formatter)
        aH.locator_params(axis='y',tight=True,nbins=6)
        cnt = cnt + 1
    fig.tight_layout()
    #plt.savefig(figname)
    plt.close()
    
def setInYaml(cur, list, value):
    if len(list) == 1:
        cur[list[0]] = value
        return
    if not list[0] in cur:
        cur[list[0]] = {}
    setInYaml(cur[list[0]], list[1:], value)


# # Setup of the GOTM/FABM model for Langtjern
# 
# ## Compiling GOTM and FABM
# 
# Please follow the instructions at <https://github.com/gotm-model/code>. Please note that the provided instructions are only available for Linux systems at present. 
# 
# For Windows, a combination of Visual Studio, CMake and the Intel Fortran compiler should work with minimal modifications when CMake is setup to generate Visual Studio project files.
# 
# For the rest of this document, it is assumed that Linux is used in combination with CMake, using GCC-5.4. 
# 
# ### Getting the source code
# 
# Credentials might be required for repositories. The text highlighted in a box should be run in a shell.
# 
# #### GOTM
# 
# Clone the source code to a directory of your choice.
# <section>
#    <div class="warn">
#    git clone https://github.com/gotm-model/code.git
#    <div/>
# </section>
# 
# Switch to lake branch.
# <section>
#    <div class="warn">
#    git checkout lake
#    <div/>
# </section>
# 
# #### FABM
# Clone the source code to a directory of your choice.
# <section>
#    <div class="warn">
#    git clone git@gitlab.au.dk:anbios/fabm-prognos.git
#    <div/>
# </section>
# 
# #### Langtjern setup
# Clone the source code to a directory of your choice.
# <section>
#    <div class="warn">
#    git clone git@gitlab.au.dh:anbios/PROGNOS.git
#    <div/>
# </section>
# 
# ### Compiling
# 
# #### FABM
# 
# From within FABM source directory.
# 
# Out-of-source builds should be preferred. Neater.
# Create a _./build_ directory in the root folder for FABM (or wherever you like but point to the right _/src_ folder in the _cmake_ command below:
# 
# <section>
#    <div class="warn">
#    mkdir build <br>
#    cd build <br>
#    cmake ../src -DFABM_HOST=gotm <br>
#    make <br>
#    make install
#    <div/>
# </section>
# 
# Note that make will create libfabm.a in the _.build_ directory. _make install_ will place the library and some other include files in the _~/local_ directory. The later is convenient to have all the builds for both FABM and GOTM in a single directory (that can be added to the PATH) thus simplifying linking.
# 
# #### GOTM
# 
# From within GOTM source directory.
# 
# Repeating the comments for the FABM compilation. Create a _./build_ directory in the root folder for GOTM. You will need to point _cmake_ to FABM root dir. Modify the cmake command below to point the the right locations if necessary:
# 
# <section>
#    <div class="warn">
#    mkdir build <br>
#    cd build <br>
#    cmake ../src -DFABM_BASE=../../fabm-gotm <br>
#    make <br>
#    make install
#    <div/>
# </section>
# 
# ### Running GOTM from the "PROGNOS/Langtjern" folder.
# 
# 
# GOTM and FABM can take a variety of inputs, depending on data availability. The format required is described in the [FABM repository](https://github.com/fabm-model/fabm/wiki/GOTM). In short, model setup relies on three different types of files:
# 
# 1. **.xml** General setup file. The different settings in _.nml_ files (see next item) will be read from the _langtjern.xml_ file inside the "_Langtjern_" folder.
# 1. **.nml** (nameless) files, to regulate the interaction between FABM and GOTM. And also to specify processes to be used within GOTM.
# 2. **.yaml** to specify which biogeochemical models will be used and their parameterizations.
# 3. **.dat** data files used to force and evaluate the model. They are in the format specified in the above link. Their are two basic flavors of data files: depth-dependent or independent.
# 
# In order to start running simulations with GOTM and FABM, some basic setup needs to be peformed first.
# 
# From the "_Langtjern_" folder in the PROGNOS repository, _editscenario.sh_ &mdash;a wrapper to a python script&mdash; needs to be run. This basically creates the <i>gotm_fabm.nml</i> and other nameless files in the "_Langtjern_" folder, filling with options read from the _Langtjern.xml_ file. 
# 
# <b>Important note:</b> The default option is NOT to use FABM when running GOTM. In order to change this *fabm_calc*  needs to be set to _.true._ in <i>gotm_fabm.nml</i>.
# 
# If you want to modify any setting there are two options:
# 
# 1. Modify the _Langtjern.xml_ file and run _editscenario.sh_.
# 2. Modify the nameless files.
# 
# #### Setting up the meteo data
# 
# In order to run GOTM some basic meteo data are required. In this case, the file is called 
# <i>meteo_file.dat</i> and should be located inside the _Langtjern_ folder. When instantiated through _editscenario.sh_ this file while be a copy of the file specified in  _langtjern.xml_. The file telling GOTM where the meteo data are located is _airsea.nml_.
# 
# In our particular setup, the meteo file was generated by running a [Python jupyter notebook](http://nbviewer.jupyter.org/github/Lecheps/PROGNOS_data/blob/master/PROGNOS_preprocessing.ipynb) that loads data from [AquaMonitor](http://www.aquamonitor.no/Portal/) and puts in in the format required by GOTM/FABM.
# 
# #### Running GOTM
# 
# The _Langtjern.xml_ file was modified to point to file generated by the Jupyter notebook that processes the AquaMonitor data.
# 
# After running "./editscenario.sh" the meteo data should be in place and GOTM should run.
# 
# <section>
#    <div class="warn">
#    ./gotm
#    <div/>
# </section>
# 

# ### Calibrating GOTM
# 
# *ACPY*, a python package developped by B&B should be installed first with the provided wheel file.
# 
# Once ACPY is installed and on the path, it should be run from the ACPY folder in the Langtjern folder. 
# 
# <section>
#    <div class="warn">
#    acpy run config_acpy.xml
#    <div/>
# </section>
# 
# A _temp.obs_ file containing the meteo data should be present in the ACPY folder. Note that the format is slightly different than the one required by FABM/GOTM.
# 
# 
# #### Results
# All the parameters tested during the calibration are saved in a sqlite database, whose name is specified in <i>config_acpy.xml</i>. Let's visualize the results of the calibration.

# In[2]:


basePath  = '/home/jose-luis/Dropbox/NIVA/PROGNOS_Development/PROGNOS/langtjern/acpy_feb01_physical_model/'
configFile = 'config_acpy_fabm.xml'

#Getting name of database from config_acpy.xml
tree = ET.parse(basePath + configFile)
root = tree.getroot()
p = tree.find("transports")
db_element = p.find("transport")
db_element = db_element.attrib
#db_path = basePath + db_element["path"]
db_path = basePath + db_element["path"] 

#Getting parameter names from config_acpy.xml
p = tree.find("parameters")
cnt = 0
variable_list = []
for i in p :
    cnt += 1
    dummy = i.attrib
    if 'variable' in dummy :
        variable_list.append(dummy["variable"])
    else :
        variable_list.append("dummy")
        
print("The variables being calibrated are: ", variable_list)
print(db_path)


# In[3]:


#Open database connection
print(db_path)
conn = sqlite3.connect(db_path)
c = conn.cursor()
#Getting the "best" parameter set according to the log-likelihood
c.execute("SELECT parameters,lnlikelihood FROM results WHERE lnlikelihood NOT NULL AND run=1 ORDER BY lnlikelihood DESC LIMIT 1")
params = c.fetchone()
log_like = float(params[1])
params = [float (x) for x in params[0].split(';')]


print("The best performing parameters are: ")
for a,b in zip(variable_list,params) :
    print (a, "\t->\t ", b)
print("With log-likelihood -> ", log_like)
conn.close()


# In[4]:


#Dotty plots simulation performed during calibration
#Storing the parameters as a numpy array
conn = sqlite3.connect(db_path)
c = conn.cursor()
c.execute("SELECT parameters,lnlikelihood FROM results WHERE lnlikelihood NOT NULL ORDER BY lnlikelihood DESC");
dotty_data = c.fetchall()
param_array = np.zeros((len(dotty_data), len(params)))
perf_array = np.zeros((len(dotty_data)))
cnt = 0
for i in dotty_data :
    #print(i[0])
    param_array[cnt][:] = np.fromstring(i[0],sep=';')
    perf_array[cnt] = float(i[1])
    cnt = cnt + 1
    
num_dots = int(len(perf_array)*0.5) 
plot_dotty(param_array[:num_dots,:],perf_array[:num_dots],variable_list,"dotty.jpg")
del param_array, perf_array
conn.close()

Image(filename='dotty.jpg')


# #### Running GOTM/FABM with the best performing parameter set
# 
# Setting parameter values in the .nml and .yaml files. 
# The list of files is read from an acpy configuration file, e.g. <i>config_acpy.xml</i>.
# The values to be set were previously obtained from the acpy result database, selecting the parameter set with the best _lnlikelihood_.

# In[5]:


#Loading parameter file information from config_acpy.xml   
tree = ET.parse(basePath + configFile)
root = tree.getroot()
p = tree.find("parameters")
parentPath = os.path.realpath(os.path.join(basePath,'..'))

#Replacing paramters in file
#This is laughingly inefficient but should do for smallish nml files
print("Setting parameter value in file:")
for child, value in zip(p,params):
    db_element = child.attrib
    if 'file' in db_element :
        _,ext = os.path.splitext(db_element['file'])
        filename = parentPath + '/' + db_element['file']
        pattern = db_element['variable'] + '\s*=.*$'
        print(filename)
        if ext == ".nml" :        
            replacement = db_element['variable'] + ' = ' +  str(value) + ','
            cnt = 0
            f = fileinput.FileInput(filename, inplace=True)
            for line in f :
                line = re.sub(pattern,replacement,line.rstrip())
                print(line)
            f.close()
        if ext == ".yaml" :
            strMap = db_element["variable"].split('/')
            with open(filename, 'r') as f:
                yamlData = yaml.round_trip_load(f)
            f.close()        
            setInYaml(yamlData,strMap,value)
            with open(filename, "w") as f:
                yaml.round_trip_dump(yamlData,f,default_flow_style=False)
            f.close()


# In[6]:


#And now running gotm
callStr = parentPath + '/' + 'gotm'
print(parentPath)
p = subprocess.Popen('./gotm', cwd = parentPath)
p.wait()


# #### Evaluating results
# 
# Loading the observed data:
# 

# In[7]:


filename = parentPath + '/output.nc'
var_to_show = ['tprof','temp']
fig_name = "langtjern_temperatures.png"
temperatureDisplay(filename,fig_name)
Image(filename=fig_name) 


# And showing some statistics for the residuals:

# In[8]:


nc = Dataset(filename, mode='r')
obs = nc.variables['tprof'][:]
obs = np.flipud(np.transpose(np.squeeze(obs)))
sim = nc.variables['temp'][:]
sim = np.flipud(np.transpose(np.squeeze(sim)))
zz = np.flipud(np.array(nc.variables['z'][0]))
tt = nc.variables['time']
startTime = datetime.datetime.strptime(tt.units,'seconds since %Y-%m-%d %H:%M:%S')
timestamp = np.squeeze(np.array(nc.variables['time'][:]))
timestamp = [(startTime + datetime.timedelta(seconds = x)) for x in timestamp]
nc.close()
residuals = obs-sim;
#fig = plt.plot(residuals.flatten(),'k-', markersize = 1)
fig = plt.hist(residuals.flatten(), 100, normed=1,facecolor='green')
#plt.show()
#plt.tight_layout()

plt.xlabel('Observed - simulated temperature ($^{o}C$)',usetex=True)
plt.savefig('residuals.png')
print(stats.describe(residuals.flatten()))    

#Various statistics
print('MAE',np.mean(np.absolute(residuals.flatten())))
print('NSE', 1 - np.sum((residuals.flatten())**2) / np.sum((obs.flatten()-np.mean(obs.flatten()))**2) ) 
print('RMSE', np.sqrt(np.mean(residuals.flatten()**2)))


# Also showing modelling results by depth:

# In[9]:


data = []
for a,b,dd in zip(obs,sim,zz) :
    dummy = dict(
                 visible = False,
                 line=dict(color='blue', width = 1),
                 x = timestamp,
                 y = a,
                 name = 'Observed temperature, depth = ' + ("%4.2f m" % (dd)) 
            )
    dummy0 = dict(
                 visible = False,
                 line=dict(color='red', width = 1),
                 x = timestamp,
                 y = b,
                 mode = 'markers',
                 name = 'Simulated temperature, depth = ' + ("%4.2f m" % (dd)) 
            )
    data.append([dummy,dummy0])
    
for i in data[0][:] :
    i['visible'] = True

flattenList = lambda x : list(itertools.chain(*x)) 
layout = go.Layout(title='Obs. and Sim. temperature at different depths', 
                   yaxis=dict(range=(-2,35)),
                  legend =dict(orientation = "h"))
fig = go.Figure(data=flattenList(data),layout=layout)

#iplot(fig, filename='dummy.html')

steps = []

#getting tempStr is farcically convoluted
for i in range(len(data)) :
    tempStr = data[i][0]['name']
    m = re.search('[0-9.]+\sm$',tempStr)
    tempStr = m.group(0)
    step = dict(
        label = tempStr,
        method = 'restyle',
        args = ['visible', [False] * len(flattenList(data))]
    )
    step['args'][1][2*i] = True # Toggle i'th trace to "visible"
    step['args'][1][2*i+1] = True
    steps.append(step)

sliders = [dict(
    active = 0,
    currentvalue = {"prefix": " Depth: "},
    pad = {"t": 50},
    steps = steps
)]

#layout.append(dict(sliders=sliders))

layout = go.Layout(title='Obs. and Sim. temperature at different depths', 
                   yaxis=dict(range=(-2,35),title='Temperature (°C)'),
                   legend =dict(orientation = "v"),
                   sliders = sliders)

fig = go.Figure(data=flattenList(data), layout=layout)

iplot(fig, filename='slider')
    
    


# In[10]:


data = []
for a,b,dd in zip(np.transpose(obs),np.transpose(sim),timestamp) :
    dummy = dict(
                 visible = False,
                 line=dict(color='blue', width = 1),
                 x = a,
                 y = np.squeeze(zz),
                 name = 'Observed temperature: ' + datetime.datetime.strftime(dd,"%d-%m-%Y")
            )
    dummy0 = dict(
                 visible = False,
                 line=dict(color='red', width = 1),
                 x = b,
                 y = np.squeeze(zz),
                 mode = 'markers',
                 name = 'Simulated values'
            )
    data.append([dummy,dummy0])

for i in data[0][:] :
    i['visible'] = True
    
#fig = go.Figure(data=flattenList(data))    
#iplot(fig, filename='dummy.html')

steps = []
for i in range(len(data)):
    tempStr = data[i][0]['name']
    m = re.search('[0-9-]+$',tempStr)
    tempStr = m.group(0)
    step = dict(
        label = tempStr,
        method = 'restyle',
        args = ['visible', [False] * len(flattenList(data))]
    )
    step['args'][1][2*i] = True # Toggle i'th trace to "visible"
    step['args'][1][2*i+1] = True
    steps.append(step)

sliders = [dict(
    active = 0,
    currentvalue = {"prefix": " Date: "},
    pad = {"t": 50},
    steps = steps
)]

layout = go.Layout(title='Obs. and Sim. temperature at different times',
                   xaxis=dict(range=(-2,25), title = 'Temperature (°C) '),
                   yaxis =dict(title='Depth (m)'),
                   legend =dict(orientation = "v"),
                   sliders = sliders) 
fig = go.Figure(data=flattenList(data), layout=layout)

iplot(fig, filename='slider')


# ### Adding chemical variables
# 
# The available chemistry data can be seen in the data pre-processing [notebook](http://nbviewer.jupyter.org/github/Lecheps/PROGNOS_data/blob/master/PROGNOS_preprocessing.ipynb).
# 
# The chemical variables are CDOM at the outlet and inlet, Total phosphorus, O2 profile, NH3, and NH4.
# 
# In order to include these data, the fabm

# In[11]:


filename = parentPath + '/output.nc'
var_to_show = ['dom_DOMb','dom_DOMa']
fig_name = "langtjern_dom.png"
tt = displayVariables(filename,var_to_show,fig_name)
Image(filename=fig_name)


# In[12]:


#Getting simulated DOMb
DOMb = getSingleArray(filename,'dom_DOMb')

#Getting DOMb at the outlet, reading from obs file
filename=basePath + 'cdom_out.obs'
DOM_out = pd.read_table(filename,'r',delimiter=" ")
DOM_out.set_index(pd.to_datetime(DOM_out.iloc[:,0] + ' ' + DOM_out.iloc[:,1] ,format = '%Y-%m-%d %H:%M:%S',utc=True),inplace=True)
DOM_out = DOM_out.drop(DOM_out.columns[[0, 1, 2]], axis=1)
DOM_out.plot()

plt.close("all")
plt.figure(num=None, figsize=(10, 4), dpi=80, facecolor='w', edgecolor='k')

plt.plot(tt,DOMb[0,:],color='r')
plt.plot(DOM_out.index,DOM_out.iloc[:,0],color='b')
plt.legend(('Top layer in lake','Observed at outlet'),fancybox=True)



# In[13]:


def css_styling():
    styles = open("./styles/custom.css", "r").read()
    return HTML(styles)
css_styling()


# In[14]:


get_ipython().run_cell_magic('javascript', '', '$(\'<div id="toc"></div>\').css({position: \'fixed\', top: \'120px\', left: 0}).appendTo(document.body);\n$.getScript(\'https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js\');\n')