library("inline")
library("Rcpp")
library("pomp")

registerPlugin(name="pomp",function(){
       Rcpp<-getPlugin("Rcpp")
       modifyList(Rcpp,list(
               LinkingTo=c('pomp',"Rcpp"),
               includes=paste(Rcpp$includes,"#include <pomp.h>",'\n'),
               Depends = c("pomp","Rcpp")
               ))
})

CppColonized<-cxxfunction(
       sig=signature(x="numeric", t="numeric", params="numeric", dt="numeric"),
       plugin='pomp',
       includes = '
       	#include <Rmath.h>
       	#include <math.h>
       	
       typedef void (*FUNC)(int,double,double*,double,double*);

       	',
       body='
       
       			FUNC p_reulermultinom = (FUNC) R_GetCCallable("pomp", "reulermultinom");
       
               GetRNGstate();
               enum Xind{C, S};  // X indices

               Rcpp::NumericVector X(x);
               double T = as<double>(t);
               Rcpp::List P(params);
               Rcpp::NumericVector D(dt);

               Rcpp::Environment global = Rcpp::Environment::global_env();
               Rcpp::NumericVector Admissions = global["Admissions"];
               Rcpp::NumericVector Discharge = global["Discharge"];

               int N = X[C] + X[S];

               double mu = P["mu"];
               double trans=0;
               double rate = 1/(1+exp(-mu))*X[C]/N;
               
               p_reulermultinom(1 , X[S], &rate, D[0], &trans);

               double nu = P["nu"];
               double import = Rf_rbinom(Admissions[T+1], 1/(1+exp( -nu)));
               double ColDis = floor(X[C]*Discharge[T+1]/N);

               std::map<std::string,double> updated;
               updated["C"] = X[C] + trans + import - ColDis;
               updated["S"] = X[S] - trans + Admissions[T+1] - import - Discharge[T+1] + ColDis;

               PutRNGstate();
               return wrap(updated);
       ',verbose=T)

       
Admissions = rpois(1,10000)
Discharge  = rpois(1,10000)
params=c(mu=log(0.02/0.98), nu=log(0.04/0.96), p=0, S.0=18000, C.0=1500)
CppColonized(x=c(C=100,S=1000),t=1,params,dt=1)