# pip install APMonitor
from APMonitor.apm import *
import numpy as np

s = 'http://byu.apmonitor.com'
a = 'measles_sq_error'

apm(s,a,'clear all')

apm_load(s,a,'measles_biweek.apm')
csv_load(s,a,'measles_biweek.csv')

# known parameters
# number of biweeks in a year
nb = 26
apm_info(s,a,'FV','N')
apm_info(s,a,'FV','MU')

apm_info(s,a,'FV','REP_FRAC')
# beta values (unknown parameters in the model)
for i in range(0,nb):
   apm_info(s,a,'FV','BETA['+str(i+1)+']')
   apm_option(s,a,'BETA['+str(i+1)+'].status',1)
# predicted values
apm_info(s,a,'SV','S')
apm_info(s,a,'SV','I')
# measured values
apm_info(s,a,'CV','CASES')

# options
# imode = 5, dynamic estimation
apm_option(s,a,'nlc.imode',5)
apm_option(s,a,'nlc.time_shift',0)
# ev_type = 1 (L1-norm) or 2 (squared error)
apm_option(s,a,'nlc.ev_type',2)
# meas_gap = deadband that represents level of accuracy / measurement noise
apm_option(s,a,'CASES.meas_gap',100)
# ctrl_units = 5, time units are in years
apm_option(s,a,'nlc.ctrl_units',5)
# turn on feedback status for CASES
apm_option(s,a,'CASES.fstatus',1)
apm_option(s,a,'CASES.wmodel',0)

# solve model and print solver output
npts = 50
gamma = np.linspace(0.05,1.5,npts)
obj = np.empty((npts))
apm_option(s,a,'gamma.lower',0.05)
apm_option(s,a,'gamma.upper',1)
apm_option(s,a,'gamma.status',1)
output = apm(s,a,'solve')
print(output)
apm_option(s,a,'gamma.upper',100)

for i in range(npts):
   apm_info(s,a,'FV','GAMMA')
   apm_option(s,a,'gamma.status',0)
   apm_meas(s,a,'gamma',gamma[i])
   output = apm(s,a,'solve')
   success = int(apm_tag(s,a,'nlc.appstatus'))
   if (success==1):
      obj[i] = apm_tag(s,a,'nlc.objfcnval')
   else:
      obj[i] = np.nan      
   print(output)

import matplotlib.pyplot as plt

plt.figure(1)
plt.plot(gamma,obj,'ro-',linewidth=2.0)
plt.xlabel('Gamma parameter value')
plt.ylabel('Objective Function')

plt.show()

# open web-interface to view results
#apm_web(s,a)
