<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Head data is generated for a pumping test in a two-aquifer model. The well starts pumping at time $t=0$ with a discharge $Q=800$ m$^3$/d. The head is measured in an observation well 10 m from the pumping well. The thickness of the aquifer is 20 m. Questions
Step1:  Model as semi-confined
<ASSISTANT_TASK:>
Python Code:
def generate_data():
    # 2 layer model with some random error
    ml = ModelMaq(kaq=[10, 20], z=[0, -20, -22, -42], c=[1000], 
                  Saq=[0.0002, 0.0001], tmin=0.001, tmax=100)
    w = Well(ml, 0, 0, rw=0.3, tsandQ=[(0, 800)])
    ml.solve()
    t = np.logspace(-2, 1, 100)
    h = ml.head(10, 0, t)
    plt.figure()
    r = 0.01 * np.random.randn(100)
    n = np.zeros_like(r)
    alpha = 0.8
    for i in range(1, len(n)):
        n[i] = 0.8 * n[i - 1] + r[i]
    ho =  h[0] + n
    plt.plot(t, ho, '.')
    data = np.zeros((len(ho), 2))
    data[:, 0] = t
    data[:, 1] = ho
    #np.savetxt('pumpingtestdata.txt', data, fmt='%2.3f', header='time (d), head (m)')
    return data
np.random.seed(11)
data = generate_data()
to = data[:, 0]
ho = data[:, 1]
def func(p, to=to, ho=ho, returnmodel=False):
    k = p[0]
    S = p[1]
    ml = ModelMaq(kaq=k, z=[0, -20], Saq=S, tmin=0.001, tmax=100)
    w = Well(ml, 0, 0, rw=0.3, tsandQ=[(0, 800)])
    ml.solve(silent=True)
    if returnmodel:
        return ml
    h = ml.head(10, 0, to)
    return np.sum((h[0] - ho) ** 2)
from scipy.optimize import fmin
lsopt = fmin(func, [10, 1e-4])
print('optimal parameters:', lsopt)
print('rmse:', np.sqrt(func(lsopt) / len(ho)))
ml = func(lsopt, returnmodel=True)
plt.figure()
plt.plot(data[:, 0], data[:, 1], '.', label='observed')
hm = ml.head(10, 0, to)
plt.plot(to, hm[0], 'r', label='modeled')
plt.legend()
plt.xlabel('time (d)')
plt.ylabel('head (m)');
cal = Calibrate(ml)
cal.set_parameter(name='kaq0', initial=10, pmin=0.1, pmax=1000)
cal.set_parameter(name='Saq0', initial=1e-4, pmin=1e-5, pmax=1e-3)
cal.series(name='obs1', x=10, y=0, layer=0, t=to, h=ho)
cal.fit(report=False)
print('rmse:', cal.rmse())
cal.parameters.style.set_precision(3)
Explanation: Head data is generated for a pumping test in a two-aquifer model. The well starts pumping at time $t=0$ with a discharge $Q=800$ m$^3$/d. The head is measured in an observation well 10 m from the pumping well. The thickness of the aquifer is 20 m. Questions:
Determine the optimal values of the hydraulic conductivity and specific storage coefficient of the aquifer when the aquifer is approximated as confined. Use a least squares approach and make use of the fmin function of scipy.optimize to find the optimal values. Plot the data with dots and the best-fit model in one graph. Print the optimal values of $k$ and $S_s$ to the screen as well as the root mean squared error of the residuals. 
Repeat Question 1 but now approximate the aquifer as semi-confined. Plot the data with dots and the best-fit model in one graph. Print to the screen the optimal values of $k$, $S_s$ and $c$  to the screen as well as the root mean squared error of the residuals. Is the semi-cofined model a better fit than the confined model?
End of explanation
def func2(p, to=to, ho=ho, returnmodel=False):
    k = p[0]
    S = p[1]
    c = p[2]
    ml = ModelMaq(kaq=k, z=[2, 0, -20], Saq=S, c=c, topboundary='semi', 
                  tmin=0.001, tmax=100)
    w = Well(ml, 0, 0, rw=0.3, tsandQ=[(0, 800)])
    ml.solve(silent=True)
    if returnmodel:
        return ml
    h = ml.head(10, 0, to)
    return np.sum((h[0] - ho) ** 2)
lsopt2 = fmin(func2, [10, 1e-4, 1000])
print('optimal parameters:', lsopt2)
print('rmse:', np.sqrt(func2(lsopt2) / len(ho)))
ml = func2(lsopt2, returnmodel=True)
plt.figure()
plt.plot(data[:, 0], data[:, 1], '.', label='observed')
hm = ml.head(10, 0, to)
plt.plot(to, hm[0], 'r', label='modeled')
plt.legend()
plt.xlabel('time (d)')
plt.ylabel('head (m)');
ml = ModelMaq(kaq=10, z=[2, 0, -20], Saq=1e-4, c=1000, topboundary='semi', tmin=0.001, tmax=100)
w = Well(ml, 0, 0, rw=0.3, tsandQ=[(0, 800)])
ml.solve(silent=True)
cal = Calibrate(ml)
cal.set_parameter(name='kaq0', initial=10)
cal.set_parameter(name='Saq0', initial=1e-4)
cal.set_parameter(name='c0', initial=1000)
cal.series(name='obs1', x=10, y=0, layer=0, t=to, h=ho)
cal.fit(report=False)
cal.parameters.style.set_precision(5)
cal.rmse(), ml.aq.kaq
plt.figure()
plt.plot(data[:, 0], data[:, 1], '.', label='observed')
hm = ml.head(10, 0, to)
plt.plot(to, hm[0], 'r', label='modeled')
plt.legend()
plt.xlabel('time (d)')
plt.ylabel('head (m)');
Explanation: Model as semi-confined
End of explanation
<END_TASK>