context('Test food demand model outputs')
library(dplyr)

test_that('Food demand model runs and produces the expected result with old parameters from Edmonds et al.', {
    ps <- 0.1
    pn <- 0.5
    y <- seq(0.1, 10, 0.1)
    samp.params <- vec2param(c(1.28,1.14,-0.19,0.21,-0.33,0.5,0.1,16,5.06,100,20))
    expect_silent(rslt <- food.dmnd(ps, pn, y, samp.params))
    expect_equal_to_reference(rslt, 'food_demand_result.rds', update=FALSE)
})


test_that('New food demand budget shares are not unreasonable', {
    parameter_data <- read.csv("test_outputs/parameter_data.csv")
    params <- vec2param(c(parameter_data$params_vector.par))
    ps <- 0.1
    pn <- 0.2
    y <- seq(0.1, 10, 0.1)
    expect_silent(rslt <- food.dmnd(ps, pn, y, params))

    rslt %>% dplyr:::mutate(alpha.t = alpha.s+alpha.n) %>%
                    filter(alpha.t >0.9)->tmp

    as.numeric(nrow(tmp))->e_rows

    expect(e_rows==0,"Budget shares generated by the food demand model are higher than 90 percent at some income levels.")
})

test_that('Income elasticities calculated are valid',{

    #First get parameters
    samp.params <- vec2param(c(1.28,1.14,-0.19,0.21,-0.33,0.5,0.1,16,5.06,100,20))

    #Get new demand
    ps <- 0.1
    pn <- 0.5
    y <- seq(0.1, 10, 0.1)


    #Calculate income elasticities on demand
    expect_silent(eta.s <- samp.params$yfunc[[1]](Y=y,FALSE))
    expect_silent(eta.n <- samp.params$yfunc[[2]](Y=y,FALSE))

    #Calculate Y terms
    expect_silent(eta.s <- samp.params$yfunc[[1]](Y=y,TRUE))
    expect_silent(eta.n <- samp.params$yfunc[[2]](Y=y,TRUE))

})

test_that('Price elasticities calculated are valid',{

    #First get parameters
    samp.params <- vec2param(c(1.28,1.14,-0.19,0.21,-0.33,0.5,0.1,16,5.06,100,20))

    #Get new demand
    ps <- 0.1
    pn <- 0.5
    y <- seq(0.1, 10, 0.1)
    rslt <- food.dmnd(ps,pn,y,samp.params)

    #Calculate income elasticities on demand
    expect_silent(eta.s <- samp.params$yfunc[[1]](Y=y,FALSE))
    expect_silent(eta.n <- samp.params$yfunc[[2]](Y=y,FALSE))

    #Calculate new price elasticities
    expect_silent(epsilon_matrix_new <- calc1eps(rslt$alpha.s,rslt$alpha.n,eta.s,eta.n,samp.params$xi))

    #Calculate old price elasticities
    Old_Demand <- readRDS("food_demand_result.rds")

    #Calculate  price elasticities
    expect_silent(epsilon_matrix_old <- calc1eps(Old_Demand$alpha.s,Old_Demand$alpha.n,eta.s,eta.n,samp.params$xi))

    expect_equal(epsilon_matrix_old,epsilon_matrix_new,tolerance=0.001,info=paste("New price elasticities are not equal to old price elasticities."))
})

test_that("Actual elasticities calculated are valid ",{

    ps <- 0.1
    pn <- 0.5
    y <- seq(0.1, 10, 0.1)
    samp.params <- vec2param(c(1.28,1.14,-0.19,0.21,-0.33,0.5,0.1,16,5.06,100,20))

    expect_silent(tmp <- calc.elas.actual(ps,pn,y,params = samp.params ))

    columns <- colnames(tmp)

    for (i in columns){
        tmpna <- tmp[is.na(toString(i)),]
        expect((nrow(tmpna))==0,"There are NA values in actual elasticity values calculated.")
    }


})

test_that("Food demand by year is reasonable compared to observations",{

    parameter_data <- read.csv("test_outputs/parameter_data.csv")
    params <- vec2param(c(parameter_data$params_vector.par))

    raw_data <- read.csv("test_outputs/Training_Data.csv") %>% filter(year %in% (2013:2015))

    (food_demand <- food.dmnd.byyear(raw_data,params = params))

    food_demand %>%  group_by(year) %>% mutate(Qs_ratio = mean(Qs/Qs.Obs), Qn_ratio= mean(Qn/Qn.Obs)) %>% ungroup()->food_demand

    food_demand$rgn <- raw_data$iso
    food_demand$income <- raw_data$gdp_pcap_thous/1000

    tmpna <- food_demand %>%
             filter(Qs_ratio>1.5  | Qn_ratio>1.5)


   expect(nrow(tmpna)==0,"The estimated values are unreasonably higher than actual values by a ratio higher than 1.5 globally for a single year. ")

})

test_that("Constant elasticities returned are valid ",{

    expect_silent(elas_func <- eta.constant(0.39))

    expect_equal(elas_func(Y=10),elas_func(Y=20),tolerance=0.01,info=("Constant elasticities functions are valid for non-staples"))

    expect_equal(elas_func(Y=10,calcQ = TRUE),elas_func(Y=11,calcQ=TRUE),tolerance=0.05,info=(" Y terms returned by constant elasticities functions are valid for non-staples"))

    expect_silent(staples_elasticity <- eta.s(0.39,9.7))
    expect_silent(staples_elasticity <- eta.s(0.39,9.7,mc.mode = TRUE))

    expect_equal(staples_elasticity(Y=10),staples_elasticity(Y=11),tolerance=0.01,info=("Constant elasticities functions are valid for staples"))

    expect_equal(staples_elasticity(Y=10,calcQ = TRUE),staples_elasticity(Y=12,calcQ=TRUE),tolerance=0.05,info=(" Y terms returned by constant elasticities functions are valid for staples"))



})