R/parameters.R

Defines functions generate_parameters lhs_parameters_parallel delete_parameter_sets lhs_parameters fix_parameters FormerFSWtoGPF

Documented in fix_parameters

#' @export
#' @useDynLib cotonou
FormerFSWtoGPF <- function(x) {
  if(!is.null(dim(x))) {
    x[4,] = x[3,]
    x[,4] = x[,3]
  }
  return(x)
}


#' parameters which depend on others, etc
#'
#' Eugene's function to be documented
#' @export
#' @useDynLib cotonou
fix_parameters <- function(y, Ncat, Nage, par_seq, condom_seq, groups_seq, years_seq) {
  useless_parameter = 1
  # why not working
  if(y$ignore_ranges_fc_c == 0) {



    if(y$prep_dropout - y$rate_leave_pro_FSW - y$muF - 1/45 - 0.008 > 0)
    {
      y$PrEP_loss_to_follow_up <- y$prep_dropout - y$rate_leave_pro_FSW - y$muF - 1/45 - 0.008
    } else {
      y$PrEP_loss_to_follow_up = 0
    }


    y$frac_women_exFSW = y$frac_women_ProFSW


    # getting client - low FSW n to be same as for pro fsw
    # y$n_y_comm_1985_LowFSW_Client = y$n_y_comm_1985_ProFSW_Client
    y$n_y_noncomm_2002_LowFSW_Client = y$n_y_noncomm_2002_ProFSW_Client
    y$n_y_noncomm_2015_LowFSW_Client = y$n_y_noncomm_2015_ProFSW_Client

    y$n_y_noncomm_1985_GPF_Client = y$n_y_noncomm_1985_GPF_GPM
    y$n_y_noncomm_1998_GPF_Client = y$n_y_noncomm_1998_GPF_GPM
    y$n_y_noncomm_2011_GPF_Client = y$n_y_noncomm_2011_GPF_GPM



    # getting low fsw condom with client for comm and noncomm partnerships same as pro fsw
    # y$fc_y_comm_1985_LowFSW_Client = y$fc_y_comm_1985_ProFSW_Client
    # y$fc_y_comm_1993_LowFSW_Client = y$fc_y_comm_1993_ProFSW_Client
    # y$fc_y_comm_1998_LowFSW_Client = y$fc_y_comm_1998_ProFSW_Client
    # y$fc_y_comm_2002_LowFSW_Client = y$fc_y_comm_2002_ProFSW_Client

    # y$fc_y_noncomm_1985_LowFSW_Client = y$fc_y_noncomm_1985_ProFSW_Client
    # y$fc_y_noncomm_2002_LowFSW_Client = y$fc_y_noncomm_2002_ProFSW_Client


    # getting client condom with GPF same as GPF-GPM
    y$fc_y_noncomm_1985_GPF_Client = y$fc_y_noncomm_1985_GPF_GPM
    y$fc_y_noncomm_1998_GPF_Client = y$fc_y_noncomm_1998_GPF_GPM
    y$fc_y_noncomm_2011_GPF_Client = y$fc_y_noncomm_2011_GPF_GPM

    # CONDOMS

    what_we_got_condom = c()
    for(year in 1:length(years_seq))
    {for(group in 1:length(groups_seq))
    {for(group2 in 1:length(groups_seq))
    {for(par in 1:length(condom_seq))
    {if(paste0(condom_seq[par], "_", years_seq[year], "_", groups_seq[group], "_", groups_seq[group2]) %in% names(y)){

      # average if the opposite value is there too - note must be same year...
      if(paste0(condom_seq[par], "_", years_seq[year], "_", groups_seq[group2], "_", groups_seq[group]) %in% names(y)) {
        av=(y[[which(paste0(condom_seq[par], "_", years_seq[year], "_", groups_seq[group], "_", groups_seq[group2]) == names(y))]] +
              y[[which(paste0(condom_seq[par], "_", years_seq[year], "_", groups_seq[group2], "_", groups_seq[group]) == names(y))]]) / 2;
        y[[which(paste0(condom_seq[par], "_", years_seq[year], "_", groups_seq[group], "_", groups_seq[group2]) == names(y))]] = av;
        y[[which(paste0(condom_seq[par], "_", years_seq[year], "_", groups_seq[group2], "_", groups_seq[group]) == names(y))]] = av;
      }

      # replacing the non-group specific parameter to include these precisions
      y[paste0(condom_seq[par], "_", years_seq[year])][[paste0(condom_seq[par], "_", years_seq[year])]][group, group2] = y[paste0(condom_seq[par], "_", years_seq[year], "_", groups_seq[group], "_", groups_seq[group2])][[1]]
      y[paste0(condom_seq[par], "_", years_seq[year])][[paste0(condom_seq[par], "_", years_seq[year])]][group2, group] = y[paste0(condom_seq[par], "_", years_seq[year])][[paste0(condom_seq[par], "_", years_seq[year])]][group, group2]
      what_we_got_condom = rbind(what_we_got_condom, c(par, year, group, group2))
    }}}}}

    if (length(what_we_got_condom[,1]) > 0)
    {
      colnames(what_we_got_condom) = c("par", "year", "group", "group2")

      # now we have to fill in the rest of the years... IF 2 OR MORE YEARS FOR SAME GROUP AND PARM
      # all years before earliest one is same as earliest estimate,
      # and all years after last one is same as last...
      # and everything in between is interpolated!
      d <- data.frame(what_we_got_condom)
      par_counts = plyr::count(d, c('par', 'group', 'group2'))

      for(i in 1:length(par_counts[,1]))
      {
        if(par_counts[i,"freq"] > 0)
        {
          the_years = subset(d, c(par == par_counts[i, "par"] & group == par_counts[i, "group"] & group2 == par_counts[i, "group2"]), year)

          #before
          if(min(the_years) > 1)
          {
            for(year in 1:(min(the_years)-1))
            {
              if(paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year]) %in% names(y))
              {y[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year])][[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year])]][par_counts[i, "group"], par_counts[i, "group2"]] =
                y[paste0(condom_seq[par_counts[i,"par"]], "_", years_seq[min(the_years)], "_", groups_seq[par_counts[i, "group"]], "_", groups_seq[par_counts[i, "group2"]])][[1]]
              # equalising the complementary condom use
              y[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year])][[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year])]][par_counts[i, "group2"], par_counts[i, "group"]] =
                y[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year])][[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year])]][par_counts[i, "group"], par_counts[i, "group2"]]

              }
            }
          }

          #after
          if(max(the_years) < length(years_seq))
          {
            for(year in (max(the_years)+1):length(years_seq))
            {
              if(paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year]) %in% names(y))
              {y[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year])][[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year])]][par_counts[i, "group"], par_counts[i, "group2"]] =
                y[paste0(condom_seq[par_counts[i,"par"]], "_", years_seq[max(the_years)], "_", groups_seq[par_counts[i, "group"]], "_", groups_seq[par_counts[i, "group2"]])][[1]]
              # equalising the complementary condom use
              y[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year])][[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year])]][par_counts[i, "group2"], par_counts[i, "group"]] =
                y[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year])][[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[year])]][par_counts[i, "group"], par_counts[i, "group2"]]
              }
            }
          }

          # interpolate between
          # take the years that have values
          # calculate the slopes between them
          # apply the slope and time difference to the pars in between

          if(par_counts[i,"freq"] > 1) {
            for(j in 1:(length(unlist(the_years))-1)) {
              slope = (y[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[unlist(the_years)][j+1])][[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[unlist(the_years)][j+1])]][par_counts[i, "group"], par_counts[i, "group2"]] -
                         y[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[unlist(the_years)][j])][[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[unlist(the_years)][j])]][par_counts[i, "group"], par_counts[i, "group2"]]) /
                (years_seq[unlist(the_years)][j+1] - years_seq[unlist(the_years)][j])

              for(k in (years_seq[unlist(the_years)][j]+1):(years_seq[unlist(the_years)][j+1]-1))
              {
                if(paste0(condom_seq[par_counts[i, "par"]], "_", k) %in% names(y))
                { y[paste0(condom_seq[par_counts[i, "par"]], "_", k)][[paste0(condom_seq[par_counts[i, "par"]], "_", k)]][par_counts[i, "group"], par_counts[i, "group2"]] = slope *
                  (k - years_seq[unlist(the_years)][j]) +
                  y[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[unlist(the_years)][j])][[paste0(condom_seq[par_counts[i, "par"]], "_", years_seq[unlist(the_years)][j])]][par_counts[i, "group"], par_counts[i, "group2"]]

                # equalising complementary condom rate
                y[paste0(condom_seq[par_counts[i, "par"]], "_", k)][[paste0(condom_seq[par_counts[i, "par"]], "_", k)]][par_counts[i, "group2"], par_counts[i, "group"]] =
                  y[paste0(condom_seq[par_counts[i, "par"]], "_", k)][[paste0(condom_seq[par_counts[i, "par"]], "_", k)]][par_counts[i, "group"], par_counts[i, "group2"]]}
              }
            }
          }


        }
      }
    }

    # MAKING FORMER FSW IN COTONOU EQUAL CONDOM USE TO GPF _ NOTE DO I NEED TO DO THIS FOR ALL GPF PARAMETERS???




    y[names(y)[grep("fc_y_", names(y))]] = lapply(invisible(y[names(y)[grep("fc_y_", names(y))]]), FormerFSWtoGPF)

    y[names(y)[grep("n_y_", names(y))]] = lapply(invisible(y[names(y)[grep("n_y_", names(y))]]), FormerFSWtoGPF)


    ##########################################################################

    # c_comm and c_noncomm
    # what_we_got = data.frame(par = c(), year = c(), group = c())
    what_we_got = c()
    for(year in 1:length(years_seq))
    {for(group in 1:length(groups_seq))
    {for(par in 1:length(par_seq))
      # {if(grepl(years_seq[year], names(y)) && grepl(groups_seq[group], names(y)) && grepl(par_seq[par], names(y))){
    {if(paste0(paste0(par_seq[par], "_", years_seq[year], "_", groups_seq[group])) %in% names(y)){

      y[paste0(par_seq[par], "_", years_seq[year])][[paste0(par_seq[par], "_", years_seq[year])]][group] = y[paste0(par_seq[par], "_", years_seq[year], "_", groups_seq[group])][[1]]
      # what_we_got = c(what_we_got, paste0(par_seq[par], "_", years_seq[year], "_", groups_seq[group]))
      what_we_got = rbind(what_we_got, c(par, year, group))
      # print(c(par_seq[par], groups_seq[group], years_seq[year]))
    }}}}
    # print(list(what_we_got, years_seq, groups_seq, par_seq, what_we_got_condom))
    if(length(what_we_got) > 0)
    {colnames(what_we_got) = c("par", "year", "group")

    # now we have to fill in the rest of the years... IF 2 OR MORE YEARS FOR SAME GROUP AND PARM
    # all years before earliest one is same as earliest estimate,
    # and all years after last one is same as last...
    # and everything in between is interpolated!
    d <- data.frame(what_we_got)
    par_counts = plyr::count(d, c('par', 'group'))

    for(i in 1:length(par_counts[,1]))
    {
      if(par_counts[i,"freq"] > 0)
      {
        the_years = subset(d, c(par == par_counts[i, "par"] & group == par_counts[i, "group"]), year)

        #before
        if(min(the_years) > 1)
        {
          for(year in 1:(min(the_years)-1))
          {
            if(paste0(par_seq[par_counts[i, "par"]], "_", years_seq[year]) %in% names(y))
            {y[paste0(par_seq[par_counts[i, "par"]], "_", years_seq[year])][[paste0(par_seq[par_counts[i, "par"]], "_", years_seq[year])]][par_counts[i, "group"]] =
              y[paste0(par_seq[par_counts[i,"par"]], "_", years_seq[min(the_years)], "_", groups_seq[par_counts[i, "group"]])][[1]]}
          }
        }

        #after
        if(max(the_years) < length(years_seq))
        {
          for(year in (max(the_years)+1):length(years_seq))
          {
            if(paste0(par_seq[par_counts[i, "par"]], "_", years_seq[year]) %in% names(y))
            {y[paste0(par_seq[par_counts[i, "par"]], "_", years_seq[year])][[paste0(par_seq[par_counts[i, "par"]], "_", years_seq[year])]][par_counts[i, "group"]] =
              y[paste0(par_seq[par_counts[i,"par"]], "_", years_seq[max(the_years)], "_", groups_seq[par_counts[i, "group"]])][[1]]}
          }
        }



        # interoplate between
        # take the years that have values
        # calculate the slopes between them
        # apply the slope and time difference to the pars in between


        if(par_counts[i,"freq"] > 1) {
          for(j in 1:(length(unlist(the_years))-1)) {
            slope = (y[paste0(par_seq[par_counts[i, "par"]], "_", years_seq[unlist(the_years)][j+1])][[paste0(par_seq[par_counts[i, "par"]], "_", years_seq[unlist(the_years)][j+1])]][par_counts[i, "group"]] -
                       y[paste0(par_seq[par_counts[i, "par"]], "_", years_seq[unlist(the_years)][j])][[paste0(par_seq[par_counts[i, "par"]], "_", years_seq[unlist(the_years)][j])]][par_counts[i, "group"]]) /
              (years_seq[unlist(the_years)][j+1] - years_seq[unlist(the_years)][j])

            for(k in (years_seq[unlist(the_years)][j]+1):(years_seq[unlist(the_years)][j+1]-1))
            {
              if(paste0(par_seq[par_counts[i, "par"]], "_", k) %in% names(y))
                y[paste0(par_seq[par_counts[i, "par"]], "_", k)][[paste0(par_seq[par_counts[i, "par"]], "_", k)]][par_counts[i, "group"]] = slope *
                  (k - years_seq[unlist(the_years)][j]) +
                  y[paste0(par_seq[par_counts[i, "par"]], "_", years_seq[unlist(the_years)][j])][[paste0(par_seq[par_counts[i, "par"]], "_", years_seq[unlist(the_years)][j])]][par_counts[i, "group"]]
            }
          }
        }
      }
    }

    }
  }
  #   print("after")
  #   print(y)


  # number of clients at initial conditions shaped by commercial partner change rates and the N of Pro FSW!
  # partnerships offered by ProFSW / P.C.R of clients / fraction of men * fraction of men that are clients
  if(y$init_clientN_from_PCR == 1)
    y$frac_men_client = y$fraction_F * y$frac_women_ProFSW * y$initial_Ntot * y$c_comm_1985[1] / y$c_comm_1985[5] / (y$initial_Ntot  * (1-y$fraction_F))


  # fractions of groups at initialisation
  y$N_init = y$initial_Ntot*c(y$fraction_F*y$frac_women_ProFSW,
                              y$fraction_F*y$frac_women_ProFSW*y$frac_women_LowFSW,
                              y$fraction_F - (y$fraction_F*y$frac_women_ProFSW + y$fraction_F*y$frac_women_ProFSW*y$frac_women_LowFSW + y$fraction_F*y$frac_women_virgin + y$fraction_F*y$frac_women_exFSW),
                              y$fraction_F*y$frac_women_exFSW,
                              (1-y$fraction_F)*y$frac_men_client,
                              (1 - y$fraction_F - ((1-y$fraction_F)*y$frac_men_client + (1-y$fraction_F) * y$frac_men_virgin)),
                              y$fraction_F*y$frac_women_virgin,
                              (1-y$fraction_F) * y$frac_men_virgin,
                              0
  )





  y$betaMtoF_FSW = y$betaMtoF_baseline * (1 + (y$RR_beta_HSV2_a_FSW - 1) * y$prev_HSV2_FSW)

  y$betaFtoM_client = y$betaMtoF_baseline * y$RR_beta_FtM * y$RR_beta_circum * (1 + (y$RR_beta_HSV2_a_client - 1) * y$prev_HSV2_Client)

  y$betaMtoF_GPF = y$betaMtoF_baseline * (1 + (y$RR_beta_HSV2_a_GPF - 1) * y$prev_HSV2_GPF)

  y$betaFtoM_GPM = y$betaMtoF_baseline * y$RR_beta_FtM * y$RR_beta_circum * (1 + (y$RR_beta_HSV2_a_GPM - 1) * y$prev_HSV2_GPM)


  # y$betaMtoF_noncomm_FSW = y$betaMtoF_baseline * (1 + (y$RR_beta_HSV2_noncomm_a - 1) * y$prev_HSV2_FSW + (y$RR_beta_HSV2_noncomm_t - 1) * y$prev_HSV2_Client)
  # y$betaFtoM_noncomm_client = y$betaMtoF_baseline * y$RR_beta_FtM * y$RR_beta_circum * (1 + (y$RR_beta_HSV2_noncomm_a - 1) * y$prev_HSV2_Client + (y$RR_beta_HSV2_noncomm_t - 1) * y$prev_HSV2_FSW)
  #
  # y$betaMtoF_noncomm = y$betaMtoF_baseline * (1 + (y$RR_beta_HSV2_noncomm_a - 1) * y$prev_HSV2_GPF + (y$RR_beta_HSV2_noncomm_t - 1) * y$prev_HSV2_GPM)
  # y$betaFtoM_noncomm = y$betaMtoF_baseline * y$RR_beta_FtM * y$RR_beta_circum * (1 + (y$RR_beta_HSV2_noncomm_a - 1) * y$prev_HSV2_GPM + (y$RR_beta_HSV2_noncomm_t - 1) * y$prev_HSV2_GPF)
  #
  # y$betaMtoF_comm = y$betaMtoF_baseline * (1 + (y$RR_beta_HSV2_comm_a - 1) * y$prev_HSV2_FSW + (y$RR_beta_HSV2_comm_t - 1) * y$prev_HSV2_Client)
  # y$betaFtoM_comm = y$betaMtoF_baseline * y$RR_beta_FtM * y$RR_beta_circum * (1 + (y$RR_beta_HSV2_comm_a - 1) * y$prev_HSV2_Client + (y$RR_beta_HSV2_comm_t - 1) * y$prev_HSV2_FSW)

  # if any beta becomes > 1, then make them all zero and flag it
  # if(y$betaMtoF_noncomm_FSW * y$infect_acute >= 1 || y$betaFtoM_noncomm_client * y$infect_acute >= 1 || y$betaMtoF_noncomm * y$infect_acute >= 1 || y$betaMtoF_comm * y$infect_acute >= 1 || y$betaFtoM_noncomm * y$infect_acute >= 1 || y$betaFtoM_comm * y$infect_acute >= 1)
  if(y$betaMtoF_FSW * y$infect_acute >= 1 || y$betaFtoM_client * y$infect_acute >= 1 || y$betaMtoF_GPF * y$infect_acute >= 1 || y$betaFtoM_GPM * y$infect_acute >= 1)
    {
    # y$betaMtoF_noncomm_FSW = 0
    # y$betaFtoM_noncomm_client = 0
    # y$betaMtoF_noncomm = 0
    # y$betaFtoM_noncomm = 0
    # y$betaMtoF_comm = 0
    # y$betaFtoM_comm = 0
    y$betaMtoF_FSW = 0
    y$betaFtoM_client = 0
    y$betaMtoF_GPF = 0
    y$betaFtoM_GPM = 0

    y$beta_above_1 = 1
  }

  y$epsilon_y = c(y$epsilon_1985, y$epsilon_1992, y$epsilon_2002, y$epsilon_2013, y$epsilon_2016)


  init_prev = if(Ncat == 9) c(y$prev_init_FSW, rep_len(y$prev_init_rest, 5), 0, 0, 0) else c(y$prev_init_FSW, rep_len(y$prev_init_rest, Ncat-1))
  y$S0_init = (1-init_prev) * y$N_init
  y$I01_init = init_prev * y$N_init

  N = y$S0_init + y$I01_init

  # BIOLOGICAL

  #dur_primary_phase, dur_200_349 input as a DURATION
  #SC_to_death input as a duration

  y$ART_RR_mort = y$ART_RR_prog

  y$gamma01 = 1/y$dur_primary_phase
  y$gamma02 = 2/(y$SC_to_death - y$dur_primary_phase - y$dur_200_349 - y$dur_below_200)
  y$gamma03 = y$gamma02
  y$gamma04 = 1/y$dur_200_349

  y$gamma01 <- rep_len(y$gamma01, Ncat)
  y$gamma02 <- rep_len(y$gamma02, Ncat)
  y$gamma03 <- rep_len(y$gamma03, Ncat)
  y$gamma04 <- rep_len(y$gamma04, Ncat)

  y$gamma11 <- y$gamma01 # loss to follow up same as normal movement?


  y$gamma22 <- y$gamma02
  y$gamma23 <- y$gamma03
  y$gamma24 <- y$gamma04

  y$gamma42 <- y$gamma02
  y$gamma43 <- y$gamma03
  y$gamma44 <- y$gamma04

  # progression is slowed by ART_RR
  y$gamma32_without_supp <- (y$gamma02)/y$ART_RR_prog
  y$gamma33_without_supp <- (y$gamma03)/y$ART_RR_prog
  y$gamma34_without_supp <- (y$gamma04)/y$ART_RR_prog


  # mortality
  y$alpha03 <- rep_len(y$alpha03, Ncat)
  y$alpha04 <- rep_len(y$alpha04, Ncat)

  y$alpha05 <- 1/y$dur_below_200
  y$alpha05 <- rep_len(y$alpha05, Ncat)

  y$alpha23 <- y$alpha03
  y$alpha24 <- y$alpha04
  y$alpha25 <- y$alpha05

  y$alpha43 <- y$alpha03
  y$alpha44 <- y$alpha04
  y$alpha45 <- y$alpha05

  y$alpha33_without_supp <- (y$alpha03)/y$ART_RR_mort
  y$alpha34_without_supp <- (y$alpha04)/y$ART_RR_mort
  y$alpha35_without_supp <- (y$alpha05)/y$ART_RR_mort



  # PCR



  y$c_t_comm = c(1985, 1993, 1995, 1998, 2002, 2005, 2008, 2012, 2015, 2016)
  y$c_t_noncomm = c(1985, 1993, 1995, 1998, 2002, 2005, 2008, 2012, 2015, 2016)

  y$c_y_comm = rbind(y$c_comm_1985, y$c_comm_1993, y$c_comm_1995, y$c_comm_1998, y$c_comm_2002,
                     y$c_comm_2005, y$c_comm_2008, y$c_comm_2012, y$c_comm_2015, y$c_comm_2016)

  y$c_y_noncomm = rbind(y$c_noncomm_1985, y$c_noncomm_1993, y$c_noncomm_1995, y$c_noncomm_1998, y$c_noncomm_2002,
                        y$c_noncomm_2005, y$c_noncomm_2008, y$c_noncomm_2012, y$c_noncomm_2015, y$c_noncomm_2016)






  y$fc_y_comm = array(data = c(y$fc_y_comm_1985, y$fc_y_comm_1993,
                               y$fc_y_comm_1995, y$fc_y_comm_1998,
                               y$fc_y_comm_2002, y$fc_y_comm_2005,
                               y$fc_y_comm_2008, y$fc_y_comm_2012,
                               y$fc_y_comm_2015, y$fc_y_comm_2016), dim=c(Ncat, Ncat, 10))

  y$fc_y_comm = aperm(y$fc_y_comm, c(3, 1, 2))

  y$fc_y_noncomm = array(data = c(y$fc_y_noncomm_1985, y$fc_y_noncomm_1993, y$fc_y_noncomm_1998, y$fc_y_noncomm_2002,
                                  y$fc_y_noncomm_2008, y$fc_y_noncomm_2011, y$fc_y_noncomm_2015,
                                  y$fc_y_noncomm_2016), dim=c(Ncat, Ncat, 8))


  y$fc_y_noncomm = aperm(y$fc_y_noncomm, c(3, 1, 2))



  y$n_y_comm = array(data = c(y$n_y_comm_1985, y$n_y_comm_2002,
                              y$n_y_comm_2015, y$n_y_comm_2016), dim=c(Ncat, Ncat, 4))
  y$n_y_noncomm = array(data = c(y$n_y_noncomm_1985, y$n_y_noncomm_1998, y$n_y_noncomm_2002, y$n_y_noncomm_2011,
                                 y$n_y_noncomm_2015, y$n_y_noncomm_2016), dim=c(Ncat, Ncat, 6))

  y$n_y_comm = aperm(y$n_y_comm, c(3, 1, 2))
  y$n_y_noncomm = aperm(y$n_y_noncomm, c(3, 1, 2))


  y$omega = N/sum(N)


  # when sampling pars, sometimes instead of length Ncat, I only sampled length one
  if(length(y$ec) == 1)
    y$ec = rep_len(y$ec, 9)
  if(length(y$eP0) == 1)
    y$eP0 = c(y$eP0, rep_len(0, 8))
  if(length(y$eP1a) == 1)
    y$eP1a = c(y$eP1a, rep_len(0, 8))
  if(length(y$eP1b) == 1)
    y$eP1b = c(y$eP1b, rep_len(0, 8))
  if(length(y$eP1c) == 1)
    y$eP1c = c(y$eP1c, rep_len(0, 8))
  if(length(y$eP1d) == 1)
    y$eP1d = c(y$eP1d, rep_len(0, 8))

  if(length(y$psia) == 1)
    y$psia = c(y$psia, rep_len(0, 8))
  if(length(y$psib) == 1)
    y$psib = c(y$psib, rep_len(0, 8))


  y$phi2 = c(y$dropout_rate_FSW, rep(y$dropout_rate_not_FSW, (Ncat-1)))
  y$phi3 = y$phi2
  y$phi4 = y$phi2
  y$phi5 = y$phi2



  if (Ncat == 9) {

    y$who_believe_comm = round(y$who_believe_comm)


    y$omega = c(
      y$fraction_F * y$frac_women_ProFSW * y$fraction_FSW_foreign, # some FSW come from outside Cotonou
      y$fraction_F * y$frac_women_LowFSW * y$frac_women_ProFSW * y$fraction_FSW_foreign,
      y$fraction_F - y$fraction_F * (1-y$fraction_sexually_active_15_F) - y$fraction_F * y$frac_women_ProFSW * y$fraction_FSW_foreign - y$fraction_F * y$frac_women_LowFSW * y$frac_women_ProFSW * y$fraction_FSW_foreign,
      0,
      0,
      (1 - y$fraction_F) - (1 - y$fraction_F) * (1 - y$fraction_sexually_active_15_M),
      y$fraction_F * (1-y$fraction_sexually_active_15_F),
      (1 - y$fraction_F) * (1 - y$fraction_sexually_active_15_M),
      0)
    #
    # MIXING
    ###############################################

    y$M_noncomm = matrix(c(0, 0, 0, 0, 1, 0, 0, 0, 0,
                           0, 0, 0, 0, 1, 0, 0, 0, 0,
                           0, 0, 0, 0, 1, 1, 0, 0, 0,
                           0, 0, 0, 0, 1, 1, 0, 0, 0,
                           1, 1, 1, 1, 0, 0, 0, 0, 0,
                           0, 0, 1, 1, 0, 0, 0, 0, 0,
                           0, 0, 0, 0, 0, 0, 0, 0, 0,
                           0, 0, 0, 0, 0, 0, 0, 0, 0,
                           0, 0, 0, 0, 0, 0, 0, 0, 0),
                         nrow = 9, ncol = 9, byrow = T)
    y$M_comm = matrix(c(0, 0, 0, 0, 1, 0, 0, 0, 0,
                        0, 0, 0, 0, 1, 0, 0, 0, 0,
                        0, 0, 0, 0, 0, 0, 0, 0, 0,
                        0, 0, 0, 0, 0, 0, 0, 0, 0,
                        1, 1, 0, 0, 0, 0, 0, 0, 0,
                        0, 0, 0, 0, 0, 0, 0, 0, 0,
                        0, 0, 0, 0, 0, 0, 0, 0, 0,
                        0, 0, 0, 0, 0, 0, 0, 0, 0,
                        0, 0, 0, 0, 0, 0, 0, 0, 0),
                      nrow = 9, ncol = 9, byrow = T)

    # MOVEMENT
    ###############################################

    y$rate_leave_low_FSW =  y$rate_leave_pro_FSW


    # virgin movement
    y$rate_move_out[7] = - y$rate_enter_sexual_pop_F
    y$rate_move_out[8] = - y$rate_enter_sexual_pop_M

    y$rate_move_in[3,7] = y$rate_enter_sexual_pop_F
    y$rate_move_in[6,8] = y$rate_enter_sexual_pop_M

    #this is just to show what happens when you increase movement
    #     y$rate_leave_pro_FSW = y$rate_leave_pro_FSW * 10
    #     y$rate_leave_low_FSW = y$rate_leave_low_FSW * 10
    #     y$rate_leave_client = y$rate_leave_client * 10

    y$prop_client_GPM = y$N_init[5] / y$N_init[6]
    y$prop_pro_FSW_GPF = y$N_init[1] / y$N_init[3]
    y$prop_low_FSW_GPF = y$N_init[2] / y$N_init[3]

    # FEMALE MOVEMENT

    y$rate_move_out[1] = - y$rate_leave_pro_FSW
    y$rate_move_out[2] = - y$rate_leave_low_FSW
    y$rate_move_out[3] = - (y$rate_leave_pro_FSW + y$muF + y$nu) * y$prop_pro_FSW_GPF * (1-y$fraction_FSW_foreign) -
      (y$rate_leave_low_FSW + y$muF + y$nu) * y$prop_low_FSW_GPF* (1-y$fraction_FSW_foreign)

    y$rate_move_in[1,3] = (y$rate_leave_pro_FSW + y$muF + y$nu) * y$prop_pro_FSW_GPF * (1-y$fraction_FSW_foreign) # moving from GPF to pro-FSW
    y$rate_move_in[2,3] = (y$rate_leave_low_FSW + y$muF + y$nu) * y$prop_low_FSW_GPF * (1-y$fraction_FSW_foreign) # moving from GPF to low-FSW

    y$rate_move_in[4,1] = y$rate_leave_pro_FSW * (1 - y$fraction_FSW_foreign) # moving from pro-FSW to former FSW in Cot
    y$rate_move_in[4,2] = y$rate_leave_low_FSW * (1 - y$fraction_FSW_foreign) # moving from low-FSW to former FSW in Cot

    y$rate_move_in[9,1] = y$rate_leave_pro_FSW * y$fraction_FSW_foreign # moving from low-FSW to former FSW NOT in Cot
    y$rate_move_in[9,2] = y$rate_leave_low_FSW * y$fraction_FSW_foreign # moving from low-FSW to former FSW NOT in Cot


    # MALE MOVEMENT

    y$rate_move_out[5] = - y$rate_leave_client
    y$rate_move_out[6] = - y$prop_client_GPM * (y$rate_leave_client + y$muM + y$nu)

    y$rate_move_in[6,5] = y$rate_leave_client  # moving from client to GPM
    y$rate_move_in[5,6] = y$prop_client_GPM * (y$rate_leave_client + y$muM + y$nu) # moving from GPM to client

    y$rate_move_out_PrEP = y$rate_move_out



    y$beta_comm = c(y$betaMtoF_FSW, y$betaMtoF_GPF, y$betaMtoF_GPF, y$betaMtoF_GPF, y$betaFtoM_client, y$betaFtoM_GPM, 0, 0, 0)
    y$beta_noncomm = c(y$betaMtoF_FSW, y$betaMtoF_GPF, y$betaMtoF_GPF, y$betaMtoF_GPF, y$betaFtoM_client, y$betaFtoM_GPM, 0, 0, 0)



    # c_t_comm = c(1985, 1993, 1995, 1998, 2002, 2005, 2008, 2012, 2015, 2016),
    # y$c_y_comm =
    y$mu = c(y$muF, y$muF, y$muF, y$muF, y$muM, y$muM, y$muF, y$muM, y$muF)

  } else {
    # y$omega = y$omega/sum(y$omega)
  }

  if(y$movement == 0) {
    y$rate_move_in = y$rate_move_in * 0
    y$rate_move_out = y$rate_move_out * 0}


  # extending all interpolating parameters beyond 2016 to the max of time

  y$c_y_comm <- rbind(y$c_y_comm, y$c_y_comm[nrow(y$c_y_comm),])
  y$c_t_comm <- c(y$c_t_comm, max(y$time))

  y$c_y_noncomm <- rbind(y$c_y_noncomm, y$c_y_noncomm[nrow(y$c_y_noncomm),])
  y$c_t_noncomm <- c(y$c_t_noncomm, max(y$time))

  y$fc_y_comm <- abind::abind(y$fc_y_comm, y$fc_y_comm[length(y$fc_y_comm[,1,1]),,], along = 1)
  y$fc_t_comm <- c(y$fc_t_comm, max(y$time))

  y$fc_y_noncomm <- abind::abind(y$fc_y_noncomm, y$fc_y_noncomm[length(y$fc_y_noncomm[,1,1]),,], along = 1)
  y$fc_t_noncomm <- c(y$fc_t_noncomm, max(y$time))


  y$n_y_comm <- abind::abind(y$n_y_comm, y$n_y_comm[length(y$n_y_comm[,1,1]),,], along = 1)
  y$n_t_comm <- c(y$n_t_comm, max(y$time))

  y$n_y_noncomm <- abind::abind(y$n_y_noncomm, y$n_y_noncomm[length(y$n_y_noncomm[,1,1]),,], along = 1)
  y$n_t_noncomm <- c(y$n_t_noncomm, max(y$time))



  y$zetaa_t <- c(y$zetaa_t, max(y$time))
  y$zetab_t <- c(y$zetab_t, max(y$time))
  y$zetac_t <- c(y$zetac_t, max(y$time))

  y$zetaa_y <- rbind(y$zetaa_y, y$zetaa_y[nrow(y$zetaa_y),])
  y$zetab_y <- rbind(y$zetab_y, y$zetab_y[nrow(y$zetab_y),])
  y$zetac_y <- rbind(y$zetac_y, y$zetac_y[nrow(y$zetac_y),])




  y$testing_prob_y[which(y$testing_prob_t == 2006), c(2,3,4)] <- y$testing_prob_women_2006
  y$testing_prob_y[which(y$testing_prob_t == 2008), c(2,3,4)] <- y$testing_prob_women_2008
  y$testing_prob_y[which(y$testing_prob_t == 2012), c(2,3,4)] <- y$testing_prob_women_2012
  y$testing_prob_y[which(y$testing_prob_t == 2013), c(2,3,4)] <- y$testing_prob_women_2012
  y$testing_prob_y[which(y$testing_prob_t == 2015), c(2,3,4)] <- y$testing_prob_women_2012
  y$testing_prob_y[which(y$testing_prob_t == 2016), c(2,3,4)] <- y$testing_prob_women_2012

  y$testing_prob_y[which(y$testing_prob_t == 2006), c(5,6)] <- y$testing_prob_men_2006
  y$testing_prob_y[which(y$testing_prob_t == 2008), c(5,6)] <- y$testing_prob_men_2008
  y$testing_prob_y[which(y$testing_prob_t == 2012), c(5,6)] <- y$testing_prob_men_2012
  y$testing_prob_y[which(y$testing_prob_t == 2013), c(5,6)] <- y$testing_prob_men_2012
  y$testing_prob_y[which(y$testing_prob_t == 2015), c(5,6)] <- y$testing_prob_men_2012
  y$testing_prob_y[which(y$testing_prob_t == 2016), c(5,6)] <- y$testing_prob_men_2012

  y$testing_prob_y <- rbind(y$testing_prob_y, y$testing_prob_y[nrow(y$testing_prob_y),])
  y$testing_prob_t <- c(y$testing_prob_t, max(y$time))
  y$ART_prob_y <- rbind(y$ART_prob_y, y$ART_prob_y[nrow(y$ART_prob_y),])
  y$ART_prob_t <- c(y$ART_prob_t, max(y$time))

  y$fPb = 1 - y$fPa - y$fPc



  # ART and efficacy
  y$viral_supp_y[1,] = y$viral_supp_y_1986_rest
  y$viral_supp_y[2,] = c(y$viral_supp_y_2015_ProFSW, rep_len(y$viral_supp_y_1986_rest, 8))
  y$viral_supp_y[3,] = y$viral_supp_y[2,]

  y$infect_ART_y = 1-(y$viral_supp_y * y$ART_eff)
  y$infect_ART_t = y$viral_supp_t




  y$rho = c(y$ART_recruit_rate_FSW, y$ART_recruit_rate_rest,
            y$ART_recruit_rate_rest,
            y$ART_recruit_rate_rest,
            y$ART_recruit_rate_rest * y$ART_init_ratio_MF,
            y$ART_recruit_rate_rest * y$ART_init_ratio_MF, 0, 0, 0)



  y$iota = c(y$ART_reinit_rate_FSW, y$ART_reinit_rate_rest,
             y$ART_reinit_rate_rest,
             y$ART_reinit_rate_rest,
             y$ART_reinit_rate_rest * y$ART_init_ratio_MF,
             y$ART_reinit_rate_rest * y$ART_init_ratio_MF, 0, 0, 0)




  if(y$long_intervention == 0)
  {
    y$tau_intervention_y = matrix(c(rep(0, 9),  y$intervention_testing_increase, c(rep(0, 8)),  y$intervention_testing_increase, c(rep(0, 8)), rep(0, 9)), ncol = 9, nrow = 4, byrow = T)
    y$rho_intervention_y = matrix(c(rep(0, 9), y$intervention_ART_increase, c(rep(0, 8)), y$intervention_ART_increase, c(rep(0, 8)), rep(0, 9)), ncol = 9, nrow = 4, byrow = T)
    y$prep_intervention_y = matrix(c(rep(0, 9), y$prep_offering_rate , rep(0, 9-1), y$prep_offering_rate , rep(0, 9-1), rep(0, 9)), ncol = 9, byrow = T)
  } else {
    y$tau_intervention_y = matrix(c(rep(0, 9),  y$intervention_testing_increase, c(rep(0, 8)),  y$intervention_testing_increase, c(rep(0, 8)), y$intervention_testing_increase, c(rep(0, 8))), ncol = 9, nrow = 4, byrow = T)
    y$rho_intervention_y = matrix(c(rep(0, 9), y$intervention_ART_increase, c(rep(0, 8)), y$intervention_ART_increase, c(rep(0, 8)), y$intervention_ART_increase, c(rep(0, 8))), ncol = 9, nrow = 4, byrow = T)
    y$prep_intervention_y = matrix(c(rep(0, 9), y$prep_offering_rate , rep(0, 9-1), y$prep_offering_rate , rep(0, 9-1), y$prep_offering_rate , rep(0, 9-1)), ncol = 9, byrow = T)

  }

  y$kappaa = c(y$PrEP_loss_to_follow_up, rep_len(0,(8)))

  # y$kappaa = c(y$prep_dropout, rep_len(0,(8)))

  y$kappab = y$kappaa
  y$kappac = y$kappaa
  y$kappa1 = y$kappaa

  y$above_500_by_group = c(y$FSW_eligible, y$GP_eligible, y$GP_eligible,
                           y$GP_eligible, y$GP_eligible, y$GP_eligible,
                           y$GP_eligible, y$GP_eligible, y$GP_eligible)

  # y$W1 = y$W1
  # y$W2 = y$W2
  # y$W3 = y$W3
  # y$W4 = y$W4

  y$pfFSW_y = matrix(c(rep(0, 9), y$prev_non_ben_fsw_1993, c(rep(0, 8)), y$prev_non_ben_fsw_2015, c(rep(0, 8)), y$prev_non_ben_fsw_2015, c(rep(0, 8))), ncol = 9, nrow = 4, byrow = T)





  return(y)
}

# N_Ncat7 = function(y) return S0_init[1] + I01_init[1] + I02_init[1] + I03_init[1] + I04_init[1] + I05_init[1]

#
# the parameters below will be sampled from an LHS and will replace their respective defaults
# unless I put something in the args of the function, eg sample = mu

#' @export
#' @useDynLib cotonou
lhs_parameters <- function(n, sample = NULL, Ncat = 9, Nage = 1, ..., set_pars = list(...), forced_pars = list(...), set_null= list(...), ranges = NULL, par_seq, condom_seq, groups_seq, years_seq) {


  set_pars <- modifyList(set_pars, forced_pars)

  # parameters that I have defined in ranges will be removed from set_pars and fixed pars (fixed pars below)
  if(length(which(names(set_pars) %in% rownames(ranges))) > 0)
    set_pars <- set_pars[-which(names(set_pars) %in% rownames(ranges))]


  #fixed pars list i think for the fix parameters function
  fixed_pars = list(
    # test_rate_prep = c(4, 0, 0, 0, 0, 0, 0, 0, 0),

    fraction_sexually_active_15_F = 0,
    fraction_sexually_active_15_M = 0,

    initial_Ntot = 286114,

    frac_women_ProFSW = 0.0024,
    frac_women_LowFSW = 3,
    frac_women_exFSW = 0.0024,

    frac_men_client = 0.2,
    frac_women_virgin = 0.1,
    frac_men_virgin = 0.1,


    fraction_F = 0.516,

    rate_move_in = matrix(0, nrow = Ncat, ncol = Ncat),
    rate_move_out = rep_len(0, Ncat),
    rate_move_out_PrEP = rep_len(0, Ncat),
    epsilon_y = 0,
    rate_enter_sexual_pop_F = 0.4,
    rate_enter_sexual_pop_M = 0.4,

    fraction_FSW_foreign = 0,
    movement = 1,
    dur_below_200 = rep_len(0.3,Ncat),
    fc_y_comm_1985 = matrix(0, Ncat, Ncat),
    fc_y_comm_1993 = matrix(0, Ncat, Ncat),
    fc_y_comm_1995 = matrix(0, Ncat, Ncat),
    fc_y_comm_1998 = matrix(0, Ncat, Ncat),
    fc_y_comm_2002 = matrix(0, Ncat, Ncat),
    fc_y_comm_2005 = matrix(0, Ncat, Ncat),
    fc_y_comm_2008 = matrix(0, Ncat, Ncat),
    fc_y_comm_2012 = matrix(0, Ncat, Ncat),
    fc_y_comm_2015 = matrix(0, Ncat, Ncat),
    fc_y_comm_2016 = matrix(0, Ncat, Ncat),

    fc_y_noncomm_1985 = matrix(0, Ncat, Ncat),
    fc_y_noncomm_1993 = matrix(0, Ncat, Ncat),
    fc_y_noncomm_1998 = matrix(0, Ncat, Ncat),
    fc_y_noncomm_2002 = matrix(0, Ncat, Ncat),

    fc_y_noncomm_2008 = matrix(0, Ncat, Ncat),
    fc_y_noncomm_2011 = matrix(0, Ncat, Ncat),
    fc_y_noncomm_2015 = matrix(0, Ncat, Ncat),
    fc_y_noncomm_2016 = matrix(0, Ncat, Ncat),

    n_y_comm_1985 = matrix(1.02, Ncat, Ncat),
    n_y_comm_2002 = matrix(1.02, Ncat, Ncat),
    n_y_comm_2015 = matrix(1.02, Ncat, Ncat),

    n_y_comm_2016 = matrix(1.02, Ncat, Ncat),

    n_y_noncomm_1985 = matrix(1.03, Ncat, Ncat),

    n_y_noncomm_1998 = matrix(1.03, Ncat, Ncat),
    n_y_noncomm_2011 = matrix(1.03, Ncat, Ncat),

    n_y_noncomm_2002 = matrix(1.03, Ncat, Ncat),
    n_y_noncomm_2015 = matrix(1.03, Ncat, Ncat),

    n_y_noncomm_2016 = matrix(1.03, Ncat, Ncat),

    c_comm_1985 = rep_len(2, Ncat),
    c_comm_1993 = rep_len(2, Ncat),
    c_comm_1995 = rep_len(2, Ncat),
    c_comm_1998 = rep_len(2, Ncat),
    c_comm_2002 = rep_len(2, Ncat),
    c_comm_2005 = rep_len(2, Ncat),
    c_comm_2008 = rep_len(2, Ncat),
    c_comm_2012 = rep_len(2, Ncat),
    c_comm_2015 = rep_len(2, Ncat),
    c_comm_2016 = rep_len(2, Ncat),

    c_noncomm_1985 = rep_len(1, Ncat),
    c_noncomm_1993 = rep_len(1, Ncat),
    c_noncomm_1995 = rep_len(1, Ncat),
    c_noncomm_1998 = rep_len(1, Ncat),
    c_noncomm_2002 = rep_len(1, Ncat),
    c_noncomm_2005 = rep_len(1, Ncat),
    c_noncomm_2008 = rep_len(1, Ncat),
    c_noncomm_2012 = rep_len(1, Ncat),
    c_noncomm_2015 = rep_len(1, Ncat),
    c_noncomm_2016 = rep_len(1, Ncat),

    betaMtoF_noncomm = 0.00193,
    betaFtoM_noncomm = 0.00867,
    betaMtoF_comm = 0.00193,
    betaFtoM_comm = 0.00867,

    betaMtoF_baseline = 0.00081,

    muF = 0.02597403,
    muM = 0.02739726,
    RR_beta_FtM = 1,
    RR_beta_circum = 0.44,

    prev_HSV2_FSW = 1,
    prev_HSV2_Client = 1,
    prev_HSV2_GPF = 1,
    prev_HSV2_GPM = 1,
    RR_beta_HSV2_comm_t = 1,
    RR_beta_HSV2_noncomm_t = 1,
    RR_beta_HSV2_comm_a = 2,
    RR_beta_HSV2_noncomm_a = 2,

    who_believe_comm = 0,
    init_clientN_from_PCR = 0,
    beta_above_1 = 0,
    ignore_ranges_fc_c = 0,
    dropout_rate_not_FSW = 0.025,
    dropout_rate_FSW = 0.025,

    delete = 0,
    nu = 0.02222222,
    PrEPOnOff = 0,
    fPa = 0.5,
    fPc = 0.4,
    iota  = c(0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0, 0, 0),
    viral_supp_t = c(1986, 2015, 2016),
    viral_supp_y = matrix(0, nrow = 3, ncol = 9),
    viral_supp_y_1986_rest = 0.6,
    viral_supp_y_2015_ProFSW = 0.7,
    ART_eff = 0.98,
    rho = c(0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0, 0, 0),
    ART_recruit_rate = 0.2,
    ART_reinit_rate = 0.1,
    RR_test_CD4200 = 5.4,

    tau_intervention_t = c(1986, 2015, 2016, 2017),
    rho_intervention_t = c(1986, 2015, 2016, 2017),

    tau_intervention_y = matrix(c(rep(0, 9), 1.5, c(rep(0, 8)), 1.5, c(rep(0, 8)), rep(0, 9)), ncol = 9, nrow = 4, byrow = T),
    rho_intervention_y = matrix(c(rep(0, 9), 6, c(rep(0, 8)), 6, c(rep(0, 8)), rep(0, 9)), ncol = 9, nrow = 4, byrow = T),

    intervention_testing_increase = 1.4,
    intervention_ART_increase = 6.5,

    psia = rep_len(0.1,9),
    psib = rep_len(0.1,9),

    prep_offering_rate = 1,
    kappaa = c(0.2, rep_len(0,(9-1))),
    kappab = c(0.2, rep_len(0,(9-1))),
    kappac = c(0.2, rep_len(0,(9-1))),
    kappa1 = c(0.2, rep_len(0,(9-1))),

    prep_dropout = 2,
    eP1a = 0.9,
    dur_FSW = 30,


    cost_ART_initiation_study_FSW = 1,
    cost_1_year_ART_study_FSW = 1,
    cost_1_year_ART_rest = 1,
    W0 = 1,
    W1 = 1,
    W2 = 1,
    W3 = 1,

    cost_Initiation_of_ART_study_FSW = 1,
    cost_Initiation_of_ART_government_FSW = 1,
    cost_1_year_of_ART_study_FSW = 1,
    cost_1_year_of_ART_government_FSW = 1,
    cost_Initiation_ART_rest_of_population = 1,
    cost_1_year_of_ART_rest_of_population = 1,
    cost_FSW_initiation_ART_Patient_costs = 1,
    cost_FSW_1_year_ART_Patient_costs = 1,

    cost_Initiation_of_PrEP_study = 1,
    cost_1_year_PrEP_perfect_adherence_study = 1,
    cost_1_year_PrEP_intermediate_adherence_study = 1,
    cost_1_year_PrEP_non_adherence_study = 1,
    cost_Initiation_of_PrEP_government = 1,
    cost_1_year_PrEP_perfect_adherence_government =1,
    cost_1_year_PrEP_intermediate_adherence_government = 1,
    cost_1_year_PrEP_non_adherence_government = 1,
    cost_PREP_initiation_Patient_costs = 1,
    cost_PREP_1_year_ART_Patient_costs = 1,
    TasP_testing = 1,
    long_intervention = 0,
    above_500_by_group = c(1, 1, 1, 1, 1, 1, 1, 1, 1),
    FSW_eligible = 1,
    GP_eligible = 1,
    prev_non_ben_fsw_1993 = 0.55,
    prev_non_ben_fsw_2015 = 0.2,
    infected_FSW_incoming = 1,
    ART_recruit_rate_rest = 0.1,
    ART_reinit_rate_rest = 0.1,

    ART_init_ratio_MF = 2,

    RR_beta_HSV2_a_FSW = 1,
    RR_beta_HSV2_a_client = 1,
    RR_beta_HSV2_a_GPF = 1,
    RR_beta_HSV2_a_GPM = 1,


    testing_prob_men_2006 = 0.111,
    testing_prob_men_2008 = 0.112,
    testing_prob_men_2012 = 0.113,

    testing_prob_women_2006 = 0.114,
    testing_prob_women_2008 = 0.115,
    testing_prob_women_2012 = 0.116,

    ART_recruit_rate_FSW = 1,
    ART_reinit_rate_FSW = 1


  )

  if(length(which(names(fixed_pars) %in% rownames(ranges))) > 0)
    fixed_pars <- fixed_pars[-which(names(fixed_pars) %in% rownames(ranges))]



  mu <- matrix(rep(c(1/50, 1/42), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("mu", Ncat), NULL))
  omega <- if(Ncat == 9) matrix(c(0.0017, 0.0067, 0, 0, 0, 0, 0, 0, 0.1, 0.2, 0, 0, 0, 0, 0, 0, 0, 0), nrow = Ncat, byrow = TRUE, dimnames = list(rep("omega", Ncat), NULL)) else
    matrix(rep(c(0.4, 0.6), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("omega", Ncat), NULL))

  # c_y_comm <- if(Ncat == 9) matrix(c(300, 1400, 40, 64, 0, 0, 0, 0, 18.67, 37.5, 0, 0, 0, 0, 0, 0, 0, 0), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_y_comm", Ncat), NULL)) else c(1, 3)

  #these parameters need to be here so fix_parameters works? below are fixed...
  S0_init = matrix(rep(c(4000, 4000), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("S0_init", Ncat), NULL))
  I01_init = matrix(rep(c(1000, 1000), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("I01_init", Ncat), NULL))



  N_init = if(Ncat == 9) matrix(c(672, 672, 757, 757, 130895, 130895, 672, 672, 27091, 27091, 100335, 100335, 14544, 14544, 11148, 11148, 0, 0), nrow = Ncat, byrow = TRUE, dimnames = list(rep("N_init", Ncat), NULL)) else c(300000, 300000)
  #   c_comm = if(Ncat == 9) matrix(c(1,1,1,1,1,1,1,1,1,1,1,1,1,1), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_comm", Ncat), NULL)) else
  #     matrix(rep(c(1,3), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_comm", Ncat), NULL))
  #   c_comm = if(Ncat == 9) matrix(c(272, 1439, 40, 64, 0, 0, 0, 0, 18.67, 37.5, 0, 0, 0, 0, 0, 0, 0, 0), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_comm", Ncat), NULL)) else
  #     matrix(rep(c(1,3), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_comm", Ncat), NULL))
  #
  #   c_noncomm = if(Ncat == 9) matrix(c(0.2729358, 0.4682779, 0.2729358, 0.4682779, 0.90, 1.02, 0.90, 1.02, 1.21, 2.5, 1.28, 1.40, 0, 0, 0, 0, 0, 0), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_noncomm", Ncat), NULL)) else
  #     matrix(rep(c(1,3), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_noncomm", Ncat), NULL))
  #

  old_ranges <- rbind(
    # c_y_comm,

    epsilon_1985 = c(0.059346131*1.5, 0.059346131*1.5),
    epsilon_1992 = c(0.053594832*1.5, 0.053594832*1.5),
    epsilon_2002 = c(0.026936907*1.5, 0.026936907*1.5),
    epsilon_2013 = c(0.026936907*1.5, 0.026936907*1.5),
    epsilon_2016 = c(0.026936907*1.5, 0.026936907*1.5),


    rate_leave_pro_FSW = c(0.2, 0.2),
    rate_leave_low_FSW = c(0.1, 0.1),
    rate_leave_client = c(0.05, 0.05),


    betaMtoF = c(0.00086, 0.00433),
    betaFtoM = c(0.00279, 0.02701),

    prev_init_FSW = c(0.01318836, 0.06592892),
    prev_init_rest = c(0.0003134459, 0.0029420363),
    #     c_comm,
    #     c_noncomm,
    mu,

    eff_ART = c(0.96, 1), # infectiousness RR when on ART
    infect_acute = c(4, 18), # RR for acute phase

    # dur_primary_phase = c(1/0.5, 1/0.16), # rate
    dur_primary_phase = c(0.16, 0.5), # from Mathieu's parameters  IN YEARS
    #     dur_primary_phase = c(2, 6.25), # from Mathieu's parameters

    SC_to_death = c(8.7, 12.3),

    dur_200_349 = c(3.9, 5), # from Mathieu's parameters  IN YEARS
    #     dur_200_349 = c(3.9, 5), # from Mathieu's parameters

    ART_RR_prog = c(8.8, 12.1),
    # ART_RR_mort = c(2, 3),

    # omega,

    #below are fixed...
    S0_init,
    I01_init,
    N_init


  )

  if (is.null(ranges)) {
    ranges <- old_ranges
  }

  if (!is.null(sample)) {
    ranges <- ranges[rownames(ranges) %in% sample,  drop=FALSE]
  }
  samples <- tgp::lhs(n, ranges)
  nms <- rownames(ranges)
  i <- split(seq_along(nms), nms)
  f <- function(x) {
    lapply(i, function(j) x[j])
  }
  samples_list <- apply(samples, 1, f)

  #   print("ehre:")
  #   print(samples_list)

  samples_list <- lapply(samples_list, function(x) modifyList(x, fixed_pars)) # btw earlier, fixed pars is replaced by ranges where they overlap

  # HERE!
  # samples_list <- lapply(samples_list, after_LHS, set_pars = set_pars)


  samples_list <- lapply(samples_list, function(x) modifyList(x, set_pars)) # set pars before fixed pars in order to get right fixed

  #   print("ehre2:")
  #   print(samples_list)
  # samples_list_test <<- samples_list

  samples_list <- lapply(samples_list, fix_parameters, Ncat = Ncat, par_seq = par_seq, condom_seq = condom_seq, groups_seq = groups_seq, years_seq = years_seq)

  # samples_list_test_fixed <<- samples_list


  # samples_list <- lapply(samples_list, function(x) modifyList(x, set_pars)) # set pars after fixed pars in order to get right set pars
  samples_list <- lapply(samples_list, function(x) modifyList(x, forced_pars)) # set pars after fixed pars in order to get right set pars

  samples_list <- lapply(samples_list, delete_parameter_sets)


  samples_list_test <<- samples_list


  if(n > 1)
    samples_list <- samples_list[!unlist(lapply(samples_list, function(x) x$delete))] # removing those runs which do not pass initial criteria


  lapply(samples_list, function(x) generate_parameters(parameters = x, Ncat = Ncat, set_null = set_null))
}

#' @export
#' @useDynLib cotonou
delete_parameter_sets <- function(x) {
  x$delete <- ifelse(x$beta_above_1 == 1, 1, # if beta is above 1 due to the various RRs

                     # if believing FSW partner change rates, then remove runs where clients pop size too big or too small ONLY IF client size is fixed at start from PCR balancing
                     ifelse(((x$c_comm_1985[1] * x$N_init[1] + x$c_comm_1985[2] * x$N_init[2])/ (x$c_comm_1985[5])) > (0.4 * sum(c(x$N_init[5], x$N_init[6], x$N_init[8]))) && (x$init_clientN_from_PCR == 1), 1,

                            # if believing client partner change rates, we are now also believing pro FSW and altering low FSW. killing runs where frac of low fsw is too low
                            ifelse(((x$c_comm_1985[5] * x$N_init[5] - x$c_comm_1985[2] * x$N_init[2])/ (x$c_comm_1985[1])) < (0* sum(x$N_init)) && (x$init_clientN_from_PCR == 0), 1, 0)))

  return(x)
}



#' @export
#' @useDynLib cotonou
lhs_parameters_parallel <- function(n, sample = NULL, Ncat = 9, Nage = 1, ..., set_pars = list(...), forced_pars = list(...), set_null= list(...), ranges = NULL, par_seq, condom_seq, groups_seq, years_seq) {


  set_pars <- modifyList(set_pars, forced_pars)

  # parameters that I have defined in ranges will be removed from set_pars and fixed pars (fixed pars below)
  if(length(which(names(set_pars) %in% rownames(ranges))) > 0)
    set_pars <- set_pars[-which(names(set_pars) %in% rownames(ranges))]


  #fixed pars list i think for the fix parameters function
  fixed_pars = list(
    # test_rate_prep = c(4, 0, 0, 0, 0, 0, 0, 0, 0),

    fraction_sexually_active_15_F = 0,
    fraction_sexually_active_15_M = 0,

    initial_Ntot = 286114,

    frac_women_ProFSW = 0.0024,
    frac_women_LowFSW = 3,
    frac_women_exFSW = 0.0024,

    frac_men_client = 0.2,
    frac_women_virgin = 0.1,
    frac_men_virgin = 0.1,


    fraction_F = 0.516,

    rate_move_in = matrix(0, nrow = Ncat, ncol = Ncat),
    rate_move_out = rep_len(0, Ncat),
    rate_move_out_PrEP = rep_len(0, Ncat),


    epsilon_y = 0,
    rate_enter_sexual_pop_F = 0.4,
    rate_enter_sexual_pop_M = 0.4,

    fraction_FSW_foreign = 0,
    movement = 1,
    dur_below_200 = rep_len(0.3,Ncat),
    fc_y_comm_1985 = matrix(0, Ncat, Ncat),
    fc_y_comm_1993 = matrix(0, Ncat, Ncat),
    fc_y_comm_1995 = matrix(0, Ncat, Ncat),
    fc_y_comm_1998 = matrix(0, Ncat, Ncat),
    fc_y_comm_2002 = matrix(0, Ncat, Ncat),
    fc_y_comm_2005 = matrix(0, Ncat, Ncat),
    fc_y_comm_2008 = matrix(0, Ncat, Ncat),
    fc_y_comm_2012 = matrix(0, Ncat, Ncat),
    fc_y_comm_2015 = matrix(0, Ncat, Ncat),
    fc_y_comm_2016 = matrix(0, Ncat, Ncat),

    fc_y_noncomm_1985 = matrix(0, Ncat, Ncat),
    fc_y_noncomm_1993 = matrix(0, Ncat, Ncat),
    fc_y_noncomm_1998 = matrix(0, Ncat, Ncat),
    fc_y_noncomm_2002 = matrix(0, Ncat, Ncat),

    fc_y_noncomm_2008 = matrix(0, Ncat, Ncat),
    fc_y_noncomm_2011 = matrix(0, Ncat, Ncat),
    fc_y_noncomm_2015 = matrix(0, Ncat, Ncat),
    fc_y_noncomm_2016 = matrix(0, Ncat, Ncat),

    n_y_comm_1985 = matrix(1.02, Ncat, Ncat),
    n_y_comm_2002 = matrix(1.02, Ncat, Ncat),
    n_y_comm_2015 = matrix(1.02, Ncat, Ncat),

    n_y_comm_2016 = matrix(1.02, Ncat, Ncat),

    n_y_noncomm_1985 = matrix(1.03, Ncat, Ncat),
    n_y_noncomm_2002 = matrix(1.03, Ncat, Ncat),
    n_y_noncomm_2015 = matrix(1.03, Ncat, Ncat),

    n_y_noncomm_1998 = matrix(1.03, Ncat, Ncat),
    n_y_noncomm_2011 = matrix(1.03, Ncat, Ncat),

    n_y_noncomm_2016 = matrix(1.03, Ncat, Ncat),

    c_comm_1985 = rep_len(2, Ncat),
    c_comm_1993 = rep_len(2, Ncat),
    c_comm_1995 = rep_len(2, Ncat),
    c_comm_1998 = rep_len(2, Ncat),
    c_comm_2002 = rep_len(2, Ncat),
    c_comm_2005 = rep_len(2, Ncat),
    c_comm_2008 = rep_len(2, Ncat),
    c_comm_2012 = rep_len(2, Ncat),
    c_comm_2015 = rep_len(2, Ncat),
    c_comm_2016 = rep_len(2, Ncat),

    c_noncomm_1985 = rep_len(1, Ncat),
    c_noncomm_1993 = rep_len(1, Ncat),
    c_noncomm_1995 = rep_len(1, Ncat),
    c_noncomm_1998 = rep_len(1, Ncat),
    c_noncomm_2002 = rep_len(1, Ncat),
    c_noncomm_2005 = rep_len(1, Ncat),
    c_noncomm_2008 = rep_len(1, Ncat),
    c_noncomm_2012 = rep_len(1, Ncat),
    c_noncomm_2015 = rep_len(1, Ncat),
    c_noncomm_2016 = rep_len(1, Ncat),

    betaMtoF_noncomm = 0.00193,
    betaFtoM_noncomm = 0.00867,
    betaMtoF_comm = 0.00193,
    betaFtoM_comm = 0.00867,

    betaMtoF_baseline = 0.00081,

    muF = 0.02597403,
    muM = 0.02739726,
    RR_beta_FtM = 1,
    RR_beta_circum = 0.44,

    prev_HSV2_FSW = 1,
    prev_HSV2_Client = 1,
    prev_HSV2_GPF = 1,
    prev_HSV2_GPM = 1,
    RR_beta_HSV2_comm_t = 1,
    RR_beta_HSV2_noncomm_t = 1,
    RR_beta_HSV2_comm_a = 2,
    RR_beta_HSV2_noncomm_a = 2,

    who_believe_comm = 0,
    init_clientN_from_PCR = 0,
    beta_above_1 = 0,
    ignore_ranges_fc_c = 0,
    dropout_rate_not_FSW = 0.025,
    dropout_rate_FSW = 0.025,

    delete = 0,
    nu = 0.02222222,
    PrEPOnOff = 0,
    fPa = 0.5,
    fPc = 0.4,
    iota  = c(0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0, 0, 0),
    viral_supp_t = c(1986, 2015, 2016),
    viral_supp_y = matrix(0, nrow = 3, ncol = 9),
    viral_supp_y_1986_rest = 0.6,
    viral_supp_y_2015_ProFSW = 0.7,
    ART_eff = 0.98,
    rho = c(0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0, 0, 0),
    ART_recruit_rate = 0.2,
    ART_reinit_rate = 0.1,
    RR_test_CD4200 = 5.4,

    tau_intervention_t = c(1986, 2015, 2016, 2017),
    rho_intervention_t = c(1986, 2015, 2016, 2017),

    tau_intervention_y = matrix(c(rep(0, 9), 1.5, c(rep(0, 8)), 1.5, c(rep(0, 8)), rep(0, 9)), ncol = 9, nrow = 4, byrow = T),
    rho_intervention_y = matrix(c(rep(0, 9), 6, c(rep(0, 8)), 6, c(rep(0, 8)), rep(0, 9)), ncol = 9, nrow = 4, byrow = T),

    intervention_testing_increase = 1.4,
    intervention_ART_increase = 6.5,

    psia = rep_len(0.1,9),
    psib = rep_len(0.1,9),

    prep_offering_rate = 1,
    kappaa = c(0.2, rep_len(0,(9-1))),
    kappab = c(0.2, rep_len(0,(9-1))),
    kappac = c(0.2, rep_len(0,(9-1))),
    kappa1 = c(0.2, rep_len(0,(9-1))),

    prep_dropout = 2,
    eP1a = 0.9,
    dur_FSW = 30,


    cost_ART_initiation_study_FSW = 1,
    cost_1_year_ART_study_FSW = 1,
    cost_1_year_ART_rest = 1,
    W0 = 1,
    W1 = 1,
    W2 = 1,
    W3 = 1,
    cost_Initiation_of_ART_study_FSW = 1,
    cost_Initiation_of_ART_government_FSW = 1,
    cost_1_year_of_ART_study_FSW = 1,
    cost_1_year_of_ART_government_FSW = 1,
    cost_Initiation_ART_rest_of_population = 1,
    cost_1_year_of_ART_rest_of_population = 1,
    cost_FSW_initiation_ART_Patient_costs = 1,
    cost_FSW_1_year_ART_Patient_costs = 1,

    cost_Initiation_of_PrEP_study = 1,
    cost_1_year_PrEP_perfect_adherence_study = 1,
    cost_1_year_PrEP_intermediate_adherence_study = 1,
    cost_1_year_PrEP_non_adherence_study = 1,
    cost_Initiation_of_PrEP_government = 1,
    cost_1_year_PrEP_perfect_adherence_government =1,
    cost_1_year_PrEP_intermediate_adherence_government = 1,
    cost_1_year_PrEP_non_adherence_government = 1,
    cost_PREP_initiation_Patient_costs = 1,
    cost_PREP_1_year_ART_Patient_costs = 1,
    TasP_testing = 1,
    long_intervention = 0,
    above_500_by_group = c(1, 1, 1, 1, 1, 1, 1, 1, 1),
    FSW_eligible = 1,
    GP_eligible = 1,

    prev_non_ben_fsw_1993 = 0.55,
    prev_non_ben_fsw_2015 = 0.2,
    infected_FSW_incoming = 1,
    ART_recruit_rate_rest = 0.1,
    ART_reinit_rate_rest = 0.1,

    ART_init_ratio_MF = 2,

    RR_beta_HSV2_a_FSW = 1,
    RR_beta_HSV2_a_client = 1,
    RR_beta_HSV2_a_GPF = 1,
    RR_beta_HSV2_a_GPM = 1,


    testing_prob_men_2006 = 0.111,
    testing_prob_men_2008 = 0.112,
    testing_prob_men_2012 = 0.113,

    testing_prob_women_2006 = 0.114,
    testing_prob_women_2008 = 0.115,
    testing_prob_women_2012 = 0.116,

    ART_recruit_rate_FSW = 1,
    ART_reinit_rate_FSW = 1



  )

  if(length(which(names(fixed_pars) %in% rownames(ranges))) > 0)
    fixed_pars <- fixed_pars[-which(names(fixed_pars) %in% rownames(ranges))]



  mu <- matrix(rep(c(1/50, 1/42), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("mu", Ncat), NULL))
  omega <- if(Ncat == 9) matrix(c(0.0017, 0.0067, 0, 0, 0, 0, 0, 0, 0.1, 0.2, 0, 0, 0, 0, 0, 0, 0, 0), nrow = Ncat, byrow = TRUE, dimnames = list(rep("omega", Ncat), NULL)) else
    matrix(rep(c(0.4, 0.6), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("omega", Ncat), NULL))

  # c_y_comm <- if(Ncat == 9) matrix(c(300, 1400, 40, 64, 0, 0, 0, 0, 18.67, 37.5, 0, 0, 0, 0, 0, 0, 0, 0), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_y_comm", Ncat), NULL)) else c(1, 3)

  #these parameters need to be here so fix_parameters works? below are fixed...
  S0_init = matrix(rep(c(4000, 4000), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("S0_init", Ncat), NULL))
  I01_init = matrix(rep(c(1000, 1000), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("I01_init", Ncat), NULL))



  N_init = if(Ncat == 9) matrix(c(672, 672, 757, 757, 130895, 130895, 672, 672, 27091, 27091, 100335, 100335, 14544, 14544, 11148, 11148, 0, 0), nrow = Ncat, byrow = TRUE, dimnames = list(rep("N_init", Ncat), NULL)) else c(300000, 300000)
  #   c_comm = if(Ncat == 9) matrix(c(1,1,1,1,1,1,1,1,1,1,1,1,1,1), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_comm", Ncat), NULL)) else
  #     matrix(rep(c(1,3), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_comm", Ncat), NULL))
  #   c_comm = if(Ncat == 9) matrix(c(272, 1439, 40, 64, 0, 0, 0, 0, 18.67, 37.5, 0, 0, 0, 0, 0, 0, 0, 0), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_comm", Ncat), NULL)) else
  #     matrix(rep(c(1,3), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_comm", Ncat), NULL))
  #
  #   c_noncomm = if(Ncat == 9) matrix(c(0.2729358, 0.4682779, 0.2729358, 0.4682779, 0.90, 1.02, 0.90, 1.02, 1.21, 2.5, 1.28, 1.40, 0, 0, 0, 0, 0, 0), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_noncomm", Ncat), NULL)) else
  #     matrix(rep(c(1,3), Ncat), nrow = Ncat, byrow = TRUE, dimnames = list(rep("c_noncomm", Ncat), NULL))
  #

  old_ranges <- rbind(
    # c_y_comm,

    epsilon_1985 = c(0.059346131*1.5, 0.059346131*1.5),
    epsilon_1992 = c(0.053594832*1.5, 0.053594832*1.5),
    epsilon_2002 = c(0.026936907*1.5, 0.026936907*1.5),
    epsilon_2013 = c(0.026936907*1.5, 0.026936907*1.5),
    epsilon_2016 = c(0.026936907*1.5, 0.026936907*1.5),


    rate_leave_pro_FSW = c(0.2, 0.2),
    rate_leave_low_FSW = c(0.1, 0.1),
    rate_leave_client = c(0.05, 0.05),


    betaMtoF = c(0.00086, 0.00433),
    betaFtoM = c(0.00279, 0.02701),

    prev_init_FSW = c(0.01318836, 0.06592892),
    prev_init_rest = c(0.0003134459, 0.0029420363),
    #     c_comm,
    #     c_noncomm,
    mu,

    eff_ART = c(0.96, 1), # infectiousness RR when on ART
    infect_acute = c(4, 18), # RR for acute phase

    # dur_primary_phase = c(1/0.5, 1/0.16), # rate
    dur_primary_phase = c(0.16, 0.5), # from Mathieu's parameters  IN YEARS
    #     dur_primary_phase = c(2, 6.25), # from Mathieu's parameters

    SC_to_death = c(8.7, 12.3),

    dur_200_349 = c(3.9, 5), # from Mathieu's parameters  IN YEARS
    #     dur_200_349 = c(3.9, 5), # from Mathieu's parameters

    ART_RR_prog = c(8.8, 12.1),
    # ART_RR_mort = c(2, 3),

    # omega,

    #below are fixed...
    S0_init,
    I01_init,
    N_init


  )

  if (is.null(ranges)) {
    ranges <- old_ranges
  }

  if (!is.null(sample)) {
    ranges <- ranges[rownames(ranges) %in% sample,  drop=FALSE]
  }
  samples <- tgp::lhs(n, ranges)
  nms <- rownames(ranges)
  i <- split(seq_along(nms), nms)
  f <- function(x) {
    lapply(i, function(j) x[j])
  }
  samples_list <- apply(samples, 1, f)

  #   print("ehre:")
  #   print(samples_list)

  samples_list <- parallel::parLapply(NULL, samples_list, function(x) modifyList(x, fixed_pars)) # btw earlier, fixed pars is replaced by ranges where they overlap

  # HERE!
  # samples_list <- lapply(samples_list, after_LHS, set_pars = set_pars)


  samples_list <- parallel::parLapply(NULL, samples_list, function(x) modifyList(x, set_pars)) # set pars before fixed pars in order to get right fixed

  #   print("ehre2:")
  #   print(samples_list)
  # samples_list_test <<- samples_list

  samples_list <- parallel::parLapply(NULL, samples_list, fix_parameters, Ncat = Ncat, par_seq = par_seq, condom_seq = condom_seq, groups_seq = groups_seq, years_seq = years_seq)

  # samples_list_test_fixed <<- samples_list


  # samples_list <- lapply(samples_list, function(x) modifyList(x, set_pars)) # set pars after fixed pars in order to get right set pars
  samples_list <- parallel::parLapply(NULL, samples_list, function(x) modifyList(x, forced_pars)) # set pars after fixed pars in order to get right set pars


  samples_list <- parallel::parLapply(NULL, samples_list, delete_parameter_sets)



  if(n > 1)
    samples_list <- samples_list[!unlist(parallel::parLapply(NULL, samples_list, function(x) x$delete))] # removing those runs which do not pass initial criteria

  parallel::parLapply(NULL, samples_list, function(x) generate_parameters(parameters = x, Ncat = Ncat, set_null = set_null))
}








#' @export
#' @useDynLib cotonou
generate_parameters <- function(..., parameters = list(...), set_null = list(...), Ncat = 9, Nage = 1) {
  defaults <- list(Ncat = Ncat,
                   Nage = Nage,

                   N_init = 300000, #ignore that this is weird - I replace with a vector later

                   initial_Ntot = 286114,
                   frac_women_ProFSW = 0.0024,
                   frac_women_LowFSW = 3,
                   frac_women_exFSW = 0.0024,

                   frac_men_client = 0.2,
                   frac_women_virgin = 0.1,
                   frac_men_virgin = 0.1,


                   above_500_by_group = c(1, 1, 1, 1, 1, 1, 1, 1, 1),
                   ART_eligible_CD4_above_500_t = c(1985, 2002, 2012, 2015, 2016),
                   ART_eligible_CD4_350_500_t = c(1985, 2002, 2012, 2015, 2016),
                   ART_eligible_CD4_200_349_t = c(1985, 2002, 2012, 2015, 2016),
                   ART_eligible_CD4_below_200_t = c(1985, 2002, 2012, 2015, 2016),

                   ART_eligible_CD4_above_500_y = c(0, 0, 0, 0, 1),
                   ART_eligible_CD4_350_500_y = c(0, 0, 0, 1, 1),
                   ART_eligible_CD4_200_349_y = c(0, 0.1, 1, 1, 1),
                   ART_eligible_CD4_below_200_y = c(0, 1, 1, 1, 1),


                   #                    epsilon_t_comm = c(1985, 1991.99, 1992, 2001.99, 2002, 2012.99, 2013, 2016),
                   #                    epsilon_y_comm = c(0.059346131, 0.059346131, 0.053594832, 0.053594832, 0.026936907, 0.026936907, 0.026936907, 0.026936907),
                   epsilon_t = c(1985, 1992, 2002, 2013, 2016),
                   epsilon_y = c(0.059346131, 0.053594832, 0.026936907, 0.026936907, 0.026936907),
                   epsilon_1985 = 0.059346131,
                   epsilon_1992 = 0.053594832,
                   epsilon_2002 = 0.026936907,
                   epsilon_2013 = 0.026936907,
                   epsilon_2016 = 0.026936907,

                   S0_init = rep_len(2000, Ncat),
                   S1a_init = rep_len(0, Ncat),
                   S1b_init = rep_len(0, Ncat),
                   S1c_init = rep_len(0, Ncat),
                   S1d_init = rep_len(0, Ncat),

                   I01_init = rep_len(1000, Ncat),
                   I11_init = rep_len(0, Ncat),
                   I02_init = rep_len(0, Ncat),
                   I03_init = rep_len(0, Ncat),
                   I04_init = rep_len(0, Ncat),
                   I05_init = rep_len(0, Ncat),
                   I22_init = rep_len(0, Ncat),
                   I23_init = rep_len(0, Ncat),
                   I24_init = rep_len(0, Ncat),
                   I25_init = rep_len(0, Ncat),
                   I32_init = rep_len(0, Ncat),
                   I33_init = rep_len(0, Ncat),
                   I34_init = rep_len(0, Ncat),
                   I35_init = rep_len(0, Ncat),
                   I42_init = rep_len(0, Ncat),
                   I43_init = rep_len(0, Ncat),
                   I44_init = rep_len(0, Ncat),
                   I45_init = rep_len(0, Ncat),

                   cumuInf_init = rep_len(0, Ncat),


                   mu = rep_len(0.02, Ncat),
                   dur_primary_phase = 0.3,
                   gamma02 = rep_len(0.2, Ncat),
                   gamma03 = rep_len(0.2, Ncat),
                   dur_200_349 = 3,

                   gamma11 = rep_len(0.2, Ncat),

                   gamma22 = rep_len(0.2, Ncat),
                   gamma23 = rep_len(0.2, Ncat),
                   gamma24 = rep_len(0.2, Ncat),

                   gamma32 = rep_len(0.2, Ncat),
                   gamma33 = rep_len(0.2, Ncat),
                   gamma34 = rep_len(0.2, Ncat),

                   gamma42 = rep_len(0.2, Ncat),
                   gamma43 = rep_len(0.2, Ncat),
                   gamma44 = rep_len(0.2, Ncat),



                   ART_prob_t = c(1985, 2005, 2012, 2015),
                   ART_prob_y = matrix(
                     rep(c(0, 0, 0.2, 0.4), Ncat), ncol = Ncat),

                   rho = c(0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0, 0, 0),

                   phi2 = rep_len(0.004,Ncat), # sort out later
                   phi3 = rep_len(0.004,Ncat),
                   phi4 = rep_len(0.004,Ncat),
                   phi5 = rep_len(0.004,Ncat),

                   psia = rep_len(0.1,Ncat),
                   psib = rep_len(0.1,Ncat),

                   testing_prob_t = c(1985, 2001, 2005, 2006, 2008, 2012, 2013, 2015, 2016),
                   testing_prob_y = matrix(
                     rep(c(0, 0.1, 0.2, 0.4, 0.5, 0.7, 0.5, 0.8, 0.7), Ncat), ncol = Ncat),

                   test_rate_prep = c(4, 0, 0, 0, 0, 0, 0, 0, 0),
                   sigma = c(0.82, 0, 0, 0, 0, 0, 0, 0, 0),
                   prep_intervention_t = c(1985, 2015, 2016, 2017),
                   prep_intervention_y = matrix(c(rep(0, Ncat), 1, rep(0, Ncat-1), 1, rep(0, Ncat-1), rep(0, Ncat)), ncol = Ncat, byrow = T), # offering rate
                   PrEPOnOff = 0,


                   RR_test_onPrEP = 2,
                   RR_test_CD4200 = 5.4,
                   RR_ART_CD4200 = 2,

                   tau01 = rep_len(1,Ncat),
                   tau11 = rep_len(1,Ncat),
                   tau2 = rep_len(1,Ncat),
                   tau3 = rep_len(1,Ncat),
                   tau4 = rep_len(1,Ncat),
                   tau5 = rep_len(1,Ncat),

                   # PREP
                   zetaa_t = c(1985, 2013, 2015, 2016),
                   zetaa_y = matrix(c(rep(0, Ncat), 0.0075, rep(0, Ncat-1), rep(0, Ncat), rep(0, Ncat)), ncol = Ncat, byrow = T),
                   zetab_t = c(1985, 2013, 2015, 2016),
                   zetab_y = matrix(c(rep(0, Ncat), 0.0075, rep(0, Ncat-1), rep(0, Ncat), rep(0, Ncat)), ncol = Ncat, byrow = T),
                   zetac_t = c(1985, 2013, 2015, 2016),
                   zetac_y = matrix(c(rep(0, Ncat), 0.0075, rep(0, Ncat-1), rep(0, Ncat), rep(0, Ncat)), ncol = Ncat, byrow = T),

                   eP = rep_len(0.6,Ncat),

                   eP0 = c(0, rep_len(0, (Ncat-1))),
                   eP1a = c(0.9, rep_len(0, (Ncat-1))),
                   eP1b = c(0.45, rep_len(0, (Ncat-1))),
                   eP1c = c(0, rep_len(0, (Ncat-1))),
                   eP1d = c(0, rep_len(0, (Ncat-1))),


                   # partner change rate
                   c_comm_1985 = rep_len(2, Ncat),
                   c_comm_1993 = rep_len(2, Ncat),
                   c_comm_1995 = rep_len(2, Ncat),
                   c_comm_1998 = rep_len(2, Ncat),
                   c_comm_2002 = rep_len(2, Ncat),
                   c_comm_2005 = rep_len(2, Ncat),
                   c_comm_2008 = rep_len(2, Ncat),
                   c_comm_2012 = rep_len(2, Ncat),
                   c_comm_2015 = rep_len(2, Ncat),
                   c_comm_2016 = rep_len(2, Ncat),

                   c_noncomm_1985 = rep_len(1, Ncat),
                   c_noncomm_1993 = rep_len(1, Ncat),
                   c_noncomm_1995 = rep_len(1, Ncat),
                   c_noncomm_1998 = rep_len(1, Ncat),
                   c_noncomm_2002 = rep_len(1, Ncat),
                   c_noncomm_2005 = rep_len(1, Ncat),
                   c_noncomm_2008 = rep_len(1, Ncat),
                   c_noncomm_2012 = rep_len(1, Ncat),
                   c_noncomm_2015 = rep_len(1, Ncat),
                   c_noncomm_2016 = rep_len(1, Ncat),

                   c_t_comm = c(1985, 1993, 1995, 1998, 2002, 2005, 2008, 2012, 2015, 2016),
                   c_t_noncomm = c(1985, 1993, 1995, 1998, 2002, 2005, 2008, 2012, 2015, 2016),

                   c_y_comm = matrix(2, nrow = 10, ncol = Ncat),
                   c_y_noncomm = matrix(1, nrow = 10, ncol = Ncat),

                   # condoms
                   #                    fc_t_noncomm = c(1985, 1990, 1998, 2016),
                   #                    fc_y_noncomm = matrix(
                   #                      rep(c(0, 0, 0.3, 0.5), Ncat), ncol = Ncat),

                   fc_y_comm_1985 = matrix(0.2, Ncat, Ncat),
                   fc_y_comm_1993 = matrix(0.5, Ncat, Ncat),
                   fc_y_comm_1995 = matrix(0.7, Ncat, Ncat),
                   fc_y_comm_1998 = matrix(0.3, Ncat, Ncat),
                   fc_y_comm_2002 = matrix(0.4, Ncat, Ncat),
                   fc_y_comm_2005 = matrix(0.1, Ncat, Ncat),
                   fc_y_comm_2008 = matrix(0.1, Ncat, Ncat),
                   fc_y_comm_2012 = matrix(0.6, Ncat, Ncat),
                   fc_y_comm_2015 = matrix(0, Ncat, Ncat),
                   fc_y_comm_2016 = matrix(0, Ncat, Ncat),

                   fc_y_noncomm_1985 = matrix(0.2, Ncat, Ncat),
                   fc_y_noncomm_1993 = matrix(0.2, Ncat, Ncat),
                   fc_y_noncomm_1998 = matrix(0.4, Ncat, Ncat),
                   fc_y_noncomm_2002 = matrix(0.3, Ncat, Ncat),
                   fc_y_noncomm_2008 = matrix(0.3, Ncat, Ncat),
                   fc_y_noncomm_2011 = matrix(0.3, Ncat, Ncat),
                   fc_y_noncomm_2015 = matrix(0.5, Ncat, Ncat),
                   fc_y_noncomm_2016 = matrix(0.5, Ncat, Ncat),



                   fc_y_comm = array(0.5,dim=c(10,Ncat,Ncat)),
                   fc_y_noncomm = array(0.5,dim=c(7,Ncat,Ncat)),

                   fc_t_comm = c(1985, 1993, 1995, 1998, 2002, 2005, 2008, 2012, 2015, 2016),
                   #                    fc_y_comm = matrix(
                   #                      rep(c(0.5, 0.5, 0.9, 0.99), Ncat), ncol = Ncat),
                   fc_t_noncomm = c(1985, 1993, 1998, 2002, 2008, 2011, 2015, 2016),
                   #                    fc_y_comm = matrix(
                   #                      rep(c(0.5, 0.5, 0.9, 0.99), Ncat), ncol = Ncat),

                   kappaa = c(0.2, rep_len(0,(Ncat-1))),
                   kappab = c(0.2, rep_len(0,(Ncat-1))),
                   kappac = c(0.2, rep_len(0,(Ncat-1))),
                   kappa1 = c(0.2, rep_len(0,(Ncat-1))),


                   alpha01 = rep_len(0,Ncat),
                   alpha02 = rep_len(0,Ncat),
                   alpha03 = 0.03,
                   alpha04 = 0.07,
                   dur_below_200 = rep_len(0.3, Ncat),

                   alpha11 = rep_len(0.01,Ncat),

                   alpha21 = rep_len(0.01,Ncat),
                   alpha22 = rep_len(0.01,Ncat),
                   alpha23 = rep_len(0.01,Ncat),
                   alpha24 = rep_len(0.01,Ncat),
                   alpha25 = rep_len(0.3,Ncat),

                   alpha32 = rep_len(0.01,Ncat),
                   alpha33 = rep_len(0.01,Ncat),
                   alpha34 = rep_len(0.01,Ncat),
                   alpha35 = rep_len(0.3,Ncat),

                   alpha42 = rep_len(0.01,Ncat),
                   alpha43 = rep_len(0.01,Ncat),
                   alpha44 = rep_len(0.01,Ncat),
                   alpha45 = rep_len(0.3,Ncat),


                   beta = rep_len(0.005,Ncat),
                   betaMtoF = 0.00193,
                   betaFtoM = 0.00867,

                   betaMtoF_noncomm = 0.00193,
                   betaFtoM_noncomm = 0.00867,
                   betaMtoF_comm = 0.00193,
                   betaFtoM_comm = 0.00867,

                   #beta = 0,

                   p_comm = matrix(1, ncol = Ncat, nrow = Ncat),
                   p_noncomm = matrix(1, ncol = Ncat, nrow = Ncat),


                   ec = rep_len(0.8,Ncat),
                   # ec = rep_len(1,1),




                   epsilon = 0.001,

                   c_comm = rep_len(2, Ncat),
                   c_noncomm = rep_len(2, Ncat),
                   #                    c_t_comm = c(1985, 1993, 1995, 1998, 2002, 2005, 2008, 2012, 2015, 2016),
                   #                    c_y_comm = matrix(rep(c(1, 2, 3, 1, 3, 4, 2, 1, 2, 4), Ncat), ncol = Ncat),
                   #
                   #                    c_t_noncomm = c(1985, 1990, 1998, 2016),
                   #                    c_y_noncomm = matrix(rep(c(0.5, 1, 0.7, 0.9), Ncat), ncol = Ncat),



                   fP_t_comm = c(1985, 2014, 2015, 2016, 2060),
                   fP_y_comm = matrix(
                     rep(c(1, 1, 1, 1, 1), Ncat), ncol = Ncat),



                   fP_t_noncomm = c(1985, 2014, 2015, 2016, 2060),
                   fP_y_noncomm = matrix(
                     rep(c(1, 1, 1, 1, 1), Ncat), ncol = Ncat),

                   n_y_comm_1985 = matrix(1.02, Ncat, Ncat),
                   n_y_comm_2002 = matrix(1.02, Ncat, Ncat),
                   n_y_comm_2015 = matrix(1.02, Ncat, Ncat),

                   n_y_comm_2016 = matrix(1.02, Ncat, Ncat),

                   n_y_noncomm_1985 = matrix(1.03, Ncat, Ncat),
                   n_y_noncomm_2002 = matrix(1.03, Ncat, Ncat),
                   n_y_noncomm_2015 = matrix(1.03, Ncat, Ncat),

                   n_y_noncomm_1998 = matrix(1.03, Ncat, Ncat),
                   n_y_noncomm_2011 = matrix(1.03, Ncat, Ncat),

                   n_y_noncomm_2016 = matrix(1.03, Ncat, Ncat),

                   #n = 0,

                   R = 1,
                   omega = rep_len(1/Ncat, Ncat),
                   theta = matrix(0.5, ncol = Ncat, nrow = Ncat),

                   M_comm = matrix(1/Ncat, Ncat, Ncat),
                   M_noncomm = matrix(1/Ncat, Ncat, Ncat),

                   # A_F = matrix(1/NAge, NAge, NAge),
                   # A_M = matrix(1/NAge, NAge, NAge),

                   #                    p = matrix(1, NAge, NAge),

                   ART_RR_prog = 10, # survival extension cofactor
                   # ART_RR_mort = 2.5, # survival extension cofactor

                   infect_ART = c(0, rep_len(0, 8)),
                   infect_acute = 9, # RR for acute phase
                   infect_AIDS = 7.27, # RR for AIDS phase

                   dur_FSW = 30,
                   OnPrEP_init = rep_len(0, Ncat),


                   # have to include the follow in the function for it to work just using generate_parameters(), and not lhs_parameters()
                   SC_to_death = 10,
                   prev_init_FSW = 0.04,
                   prev_init_rest = 0.0008,
                   rate_leave_pro_FSW = 0.2,
                   rate_leave_client = 0.05,
                   rate_leave_low_FSW = 0.1,
                   prop_client_GPM = 0.2430057, # 27091/111483
                   prop_pro_FSW_GPF = 0.004620494, # 672 / 145439
                   prop_low_FSW_GPF = 0.005204931, # 757 / 145439
                   rate_move_in = matrix(0, ncol = Ncat, nrow = Ncat),
                   rate_move_out = rep_len(0, Ncat),
                   rate_move_out_PrEP = rep_len(0, Ncat),

                   rate_enter_sexual_pop_F = 1,
                   rate_enter_sexual_pop_M = 1,

                   fraction_F = 0.51,
                   fraction_FSW_foreign = 0.5,
                   replaceDeaths = 0,
                   movement = 1,
                   beta_comm = rep_len(0.002, Ncat),

                   beta_noncomm = rep_len(0.001, Ncat),
                   muF = 0.02597403,
                   muM = 0.02739726,
                   RR_beta_FtM = 1,
                   RR_beta_circum = 0.44,

                   prev_HSV2_FSW = 1,
                   prev_HSV2_Client = 1,
                   prev_HSV2_GPF = 1,
                   prev_HSV2_GPM = 1,
                   RR_beta_HSV2_comm_t = 1,
                   RR_beta_HSV2_noncomm_t = 1,
                   RR_beta_HSV2_comm_a = 2,
                   RR_beta_HSV2_noncomm_a = 2,

                   betaMtoF_baseline = 0.00081,


                   who_believe_comm = 0,
                   init_clientN_from_PCR = 0,
                   fraction_sexually_active_15_F = 0,
                   fraction_sexually_active_15_M = 0,
                   beta_above_1 = 0,
                   ignore_ranges_fc_c = 0,
                   dropout_rate_not_FSW = 0.025,
                   dropout_rate_FSW = 0.025,

                   delete = 0,
                   nu = 0.02222222,
                   fPa = 0.5,
                   fPc = 0.4,
                   iota  = c(0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0, 0, 0),
                   viral_supp_t = c(1985, 2015, 2016),
                   viral_supp_y = matrix(0, nrow = 3, ncol = 9),
                   viral_supp_y_1986_rest = 0.6,
                   viral_supp_y_2015_ProFSW = 0.7,
                   ART_eff = 0.98,
                   ART_recruit_rate = 0.2,
                   ART_reinit_rate = 0.1,

                   tau_intervention_t = c(1986, 2015, 2016, 2017),
                   rho_intervention_t = c(1986, 2015, 2016, 2017),

                   tau_intervention_y = matrix(c(rep(0, 9), 1.5, c(rep(0, 8)), 1.5, c(rep(0, 8)), rep(0, 9)), ncol = 9, nrow = 4, byrow = T),
                   rho_intervention_y = matrix(c(rep(0, 9), 6, c(rep(0, 8)), 6, c(rep(0, 8)), rep(0, 9)), ncol = 9, nrow = 4, byrow = T),

                   intervention_testing_increase = 1.4,
                   intervention_ART_increase = 6.5,

                   prep_offering_rate = 1,

                   FSW_ONLY = c(1, 1, rep(0, 7)),


                   cost_ART_initiation_study_FSW = 1,
                   cost_1_year_ART_study_FSW = 1,
                   cost_1_year_ART_rest = 1,
                   W0 = 1,
                   W1 = 1,
                   W2 = 1,
                   W3 = 1,
                   cost_Initiation_of_ART_study_FSW = 1,
                   cost_Initiation_of_ART_government_FSW = 1,
                   cost_1_year_of_ART_study_FSW = 1,
                   cost_1_year_of_ART_government_FSW = 1,
                   cost_Initiation_ART_rest_of_population = 1,
                   cost_1_year_of_ART_rest_of_population = 1,
                   cost_FSW_initiation_ART_Patient_costs = 1,
                   cost_FSW_1_year_ART_Patient_costs = 1,

                   cost_Initiation_of_PrEP_study = 1,
                   cost_1_year_PrEP_perfect_adherence_study = 1,
                   cost_1_year_PrEP_intermediate_adherence_study = 1,
                   cost_1_year_PrEP_non_adherence_study = 1,
                   cost_Initiation_of_PrEP_government = 1,
                   cost_1_year_PrEP_perfect_adherence_government =1,
                   cost_1_year_PrEP_intermediate_adherence_government = 1,
                   cost_1_year_PrEP_non_adherence_government = 1,
                   cost_PREP_initiation_Patient_costs = 1,
                   cost_PREP_1_year_ART_Patient_costs = 1,
                   TasP_testing = 1,
                   long_intervention = 0,
                   FSW_eligible = 1,
                   GP_eligible = 1,
                   pfFSW_t = c(1985, 1993, 2015, 2060),

                   pfFSW_y = matrix(c(rep(0, 9), 0.55, c(rep(0, 8)), 0.2, c(rep(0, 8)), 0.2, c(rep(0, 8))), ncol = 9, nrow = 4, byrow = T),

                   prev_non_ben_fsw_1993 = 0.55,
                   prev_non_ben_fsw_2015 = 0.2,


                   infected_FSW_incoming = 1,

                   prep_efficacious_y = c(0, 1, 1, 0, 0),
                   prep_efficacious_t = c(1985, 2015, 2017, 2017.01, 2060),

                   PrEP_loss_to_follow_up = 0.1,

                   ART_recruit_rate_rest = 0.1,
                   ART_reinit_rate_rest = 0.1,

                   ART_init_ratio_MF = 2,

                   RR_beta_HSV2_a_FSW = 1,
                   RR_beta_HSV2_a_client = 1,
                   RR_beta_HSV2_a_GPF = 1,
                   RR_beta_HSV2_a_GPM = 1,


                   testing_prob_men_2006 = 0.111,
                   testing_prob_men_2008 = 0.112,
                   testing_prob_men_2012 = 0.113,

                   testing_prob_women_2006 = 0.114,
                   testing_prob_women_2008 = 0.115,
                   testing_prob_women_2012 = 0.116,

                   ART_recruit_rate_FSW = 1,
                   ART_reinit_rate_FSW = 1






  )



  if (length(parameters) == 0L) {
    return(defaults)
  }

  if (is.null(names(parameters)) || !all(nzchar(names(parameters)))) {
    stop("All arguments must be named")
  }
  #   extra <- setdiff(names(parameters), names(defaults))
  #   if (length(extra) > 0L) {
  #     stop("Unknown arguments: ", extra)
  #   }

  # list of parameters that depend on others
  ret <- modifyList(defaults, parameters)

  if (length(set_null) > 0L) {
    ret = modifyList(ret, lapply(ret[match(set_null, names(ret))], function(x) x*0))
  }

  ret
}
Geidelberg/cotonou documentation built on Nov. 15, 2018, 6:49 p.m.