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

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

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
beta = np.empty((nb))
apm_info(s,a,'FV','N')
apm_info(s,a,'FV','MU')
apm_info(s,a,'FV','GAMMA')
apm_option(s,a,'gamma.status',1)
apm_option(s,a,'gamma.upper',0.5);
apm_option(s,a,'gamma.lower',0.05);

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)
# ev_type = 1 (L1-norm) or 2 (squared error)
apm_option(s,a,'nlc.ev_type',1)
# 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
output = apm(s,a,'solve')
print(output)

est = apm_sol(s,a)

apm_web(s,a)


################################
a2 = 'control'
################################

apm(s,a2,'clear all')
apm_load(s,a2,'measles_biweek.apm')
csv_load(s,a2,'measles_biweek.csv')

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

# options
# imode = 6, dynamic control
apm_option(s,a2,'nlc.imode',6)
# cv_type = 1 (L1-norm) or 2 (squared error)
apm_option(s,a2,'nlc.cv_type',1)
# ctrl_units = 5, time units are in years
apm_option(s,a2,'nlc.ctrl_units',5)

# turn off beta and gamma values
apm_option(s,a2,'gamma.status',0)
# turn off cases feedback status (no measurements)
apm_option(s,a2,'CASES.fstatus',0)
# configure controller objective
apm_option(s,a2,'CASES.status',1)
apm_option(s,a2,'CASES.tr_init',0)
apm_option(s,a2,'CASES.sphi',50)
apm_option(s,a2,'CASES.splo',0)
# no hard constraints - makes it infeasible
# apm_option(s,a2,'CASES.upper',50)
# apm_option(s,a2,'CASES.lower',0)

# insert parameters from estimation
gamma = apm_tag(s,a,'gamma.newval')
apm_meas(s,a2,'gamma',gamma)
for i in range(0,nb):
   apm_option(s,a2,'BETA['+str(i+1)+'].status',0)
   apm_option(s,a2,'BETA['+str(i+1)+'].fstatus',1)
   beta[i] = apm_tag(s,a,'BETA['+str(i+1)+'].newval')
   apm_meas(s,a2,'BETA['+str(i+1)+']',beta[i])
# enable vaccine rate as a decision
apm_option(s,a2,'Vr.status',1)
apm_option(s,a2,'Vr.upper',10000)
apm_option(s,a2,'Vr.lower',0)
apm_option(s,a2,'Vr.cost',1e-5)
# lower bound for vaccine supply
apm_option(s,a2,'V.lower',0)
apm_option(s,a2,'V.sphi',200000)
apm_option(s,a2,'V.splo',0)
apm_option(s,a2,'V.status',0)
apm_option(s,a2,'V.tr_init',0)

# solve model and print solver output
output = apm(s,a2,'solve')
print(output)

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