%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from ttim import *

"""
Explanation: Confined Aquifer Test
This test is taken from AQTESOLV examples.
End of explanation
"""


b = -5.4846 #aquifer thickness in m
Q = 1199.218 #constant discharge in m^3/d
r = 251.1552 #distance between observation well to test well in m
rw = 0.1524 #screen radius of test well in m

"""
Explanation: Set basic parameters:
End of explanation
"""


data1 = np.loadtxt('data/gridley_well_1.txt')
t1 = data1[:, 0]
h1 = data1[:, 1]
data2 = np.loadtxt('data/gridley_well_3.txt')
t2 = data2[:, 0]
h2 = data2[:, 1]

"""
Explanation: Load dataset:
End of explanation
"""


ml = ModelMaq(kaq=10, z=[0, b], Saq=0.001, tmin=0.001, tmax=1, topboundary='conf')
w = Well(ml, xw=0, yw=0, rw=rw, tsandQ=[(0, Q)], layers=0)
ml.solve()

"""
Explanation: Create conceptual model:
End of explanation
"""


#unknown parameters: kaq, Saq
ca_0 = Calibrate(ml)
ca_0.set_parameter(name='kaq0', initial=10)
ca_0.set_parameter(name='Saq0', initial=1e-4)
ca_0.series(name='obs1', x=r, y=0, t=t1, h=h1, layer=0)
ca_0.fit(report=True)

display(ca_0.parameters)
print('rmse:', ca_0.rmse())

hm_0 = ml.head(r, 0, t1)
plt.figure(figsize = (8, 5))
plt.semilogx(t1, hm_0[0], label = 'ttim')
plt.semilogx(t1, h1, '.', label = 'obs1')
plt.xlabel('time(d)')
plt.ylabel('drawdown(m)')
plt.legend();

#unknown parameters: kaq, Saq
ca_1 = Calibrate(ml)
ca_1.set_parameter(name='kaq0', initial=10)
ca_1.set_parameter(name='Saq0', initial=1e-4)
ca_1.series(name='obs2', x=0, y=0, t=t2, h=h2, layer=0)
ca_1.fit(report=True)

display(ca_1.parameters)
print('rmse:', ca_1.rmse())

hm_1 = ml.head(0, 0, t2)
plt.figure(figsize = (8, 5))
plt.semilogx(t2, hm_1[0], label = 'ttim')
plt.semilogx(t2, h2, '.', label = 'obs2')
plt.xlabel('time(d)')
plt.ylabel('drawdown(m)')
plt.legend();

"""
Explanation: Calibrate with two datasets simultaneously:
End of explanation
"""


ml_1 = ModelMaq(kaq=10, z=[0, b], Saq=0.001, tmin=0.001, tmax=1, topboundary='conf')
w_1 = Well(ml_1, xw=0, yw=0, rw=rw, tsandQ=[(0, Q)], layers=0)
ml_1.solve()

ca_2 = Calibrate(ml_1)
ca_2.set_parameter(name='kaq0', initial=10)
ca_2.set_parameter(name='Saq0', initial=1e-4, pmin=0)
ca_2.series(name='obs1', x=r, y=0, t=t1, h=h1, layer=0)
ca_2.series(name='obs2', x=0, y=0, t=t2, h=h2, layer=0)
ca_2.fit(report=True)

display(ca_2.parameters)
print('rmse:', ca_2.rmse())

hm1_2 = ml.head(r, 0, t1)
hm2_2 = ml.head(0, 0, t2)
plt.figure(figsize = (8, 5))
plt.semilogx(t1, hm1_2[0], label = 'ttim1')
plt.semilogx(t1, h1, '.', label = 'obs1')
plt.semilogx(t2, hm2_2[0], label = 'ttim3')
plt.semilogx(t2, h2, '.', label = 'obs3')
plt.xlabel('time(d)')
plt.ylabel('drawdown(m)')
plt.legend();

"""
Explanation: Calibrate with two datasets simultaneously:
End of explanation
"""


ml_2 = ModelMaq(kaq=10, z=[0, b], Saq=0.001, tmin=0.001, tmax=1, topboundary='conf')
w_2 = Well(ml_2, xw=0, yw=0, rw=rw, rc=0.2, res=0.2, tsandQ=[(0, Q)], layers=0)
ml_2.solve()

"""
Explanation: Ty adding well skin resistance and wellbore storage:
End of explanation
"""


ca_3 = Calibrate(ml_2)
ca_3.set_parameter(name = 'kaq0', initial = 10)
ca_3.set_parameter(name = 'Saq0', initial = 1e-4, pmin=0)
ca_3.set_parameter_by_reference(name='res', parameter=w_2.res, initial =0.2)
ca_3.series(name='obs1', x=r, y=0, t=t1, h=h1, layer=0)
ca_3.series(name='obs3', x=0, y=0, t=t2, h=h2, layer=0)
ca_3.fit(report=True)

display(ca_3.parameters)
print('rmse:', ca_3.rmse())

hw1 = ml_2.head(r, 0, t1)
hw2 = ml_2.head(0, 0, t2)
plt.figure(figsize = (8, 5))
plt.semilogx(t1, hw1[0], label = 'ttim1')
plt.semilogx(t1, h1, '.', label = 'obs1')
plt.semilogx(t2, hw2[0], label = 'ttim3')
plt.semilogx(t2, h2, '.', label = 'obs3')
plt.xlabel('time(d)')
plt.ylabel('drawdown(m)')
plt.legend();

"""
Explanation: If adding wellbore sotrage to the parameters to be optimized, the fit gives extremely large values of each parameter which is imposiible. However, when remove rc from well function, the fit cannot be completed with uncertainties. Thus, the rc value is determined as 0.2 by trial-and-error procedure.
End of explanation
"""


t = pd.DataFrame(columns=['k [m/d]', 'Ss [1/m]', 'res'], \
                index=['MLU', 'AQTESOLV', 'ttim', 'ttim-res&rc'])
t.loc['MLU'] = [38.094, 1.193E-06, '-']
t.loc['AQTESOLV'] = [37.803, 1.356E-06, '-']
t.loc['ttim'] = np.append(ca_2.parameters['optimal'].values, '-')
t.loc['ttim-res&rc'] = ca_3.parameters['optimal'].values 
t['rc'] = ['-', '-', '-', 0.2]
t['RMSE'] = [0.259, 0.270, 0.272, 0.192]
t

"""
Explanation: Summary of values simulated by AQTESOLV and MLU
The results simulated by different methods with two datasets simultaneously are presented below. In the example of AQTESOLV, result simulated with only observation well is presented. The comparision of results when only observation well is included can be found in the report related to this test.
End of explanation
"""


