DAISIEsim:

DAISIE_sim = function(
  time,
  M,
  pars,
  replicates,
  mainland_params = NULL,
  divdepmodel = 'CS',
  island_type = 'oceanic', # 'oceanic' = intially 0 species; 'nonoceanic' = requires probability of
  #being sampled and probability of being non-endemic
  nonoceanic = NULL,
  prop_type2_pool = NA,
  replicates_apply_type2 = TRUE,
  sample_freq = 25,
  plot_sims = TRUE,
  island_ontogeny = 'const', # const = no effect; 'linear' = linear decreasing function; 'beta' = beta function;
  sea_level = 'const', # const = no effect; 'linear_pos' = linear increasing function,
  #'linear_neg' = linear decreasing function; 'sine' = symmetric fluctuations of a sin wave
  Apars = NULL,
  Epars = NULL,
  divdep = c('lac', 'gam'), #'lac is cladogenesis, 'mu' is extinction, 'gam' is immigration,and any combination
  keep_final_state = FALSE,
  stored_data = NULL,
  verbose = TRUE,
  ...
) {

  testit::assert(
    "island_ontogeny is not valid input. Specify 'const',\n
    'linear' or  ' beta'", DAISIE::is_island_ontogeny_input(island_ontogeny)
  )

  #TODO: TEST island_replicates INPUT! SANITIZE STORED_DATA INPUT! ASSERT + TEST

  #TODO: write test for oceanic and nonoceanic island type

  #TODO: write test for nonoceanic sampling parameters

  if (!is.null(stored_data)) {
    start_midway <- TRUE
  } else {
    start_midway <- FALSE
  }

  # @richelbilderbeek
  if (!is.null(mainland_params)) {
    return(
      DAISIE_sim_with_mainland(
        time = time,
        M = M,
        pars = pars,
        replicates = replicates,
        mainland_params = mainland_params,
        divdepmodel = divdepmodel,
        prop_type2_pool = prop_type2_pool,
        replicates_apply_type2 = replicates_apply_type2,
        sample_freq = sample_freq
      )
    )
  }

  # Classic behavior
  totaltime <- time
  island_replicates  = list()

  if(divdepmodel =='IW')
  {
    if(length(pars) > 5)
    {
      stop('Island-wide carrying capacity model not yet implemented for two types of mainland species')
    }
    for(rep in 1:replicates)
    {
      island_replicates[[rep]] <-DAISIE_sim_core(
        time = totaltime,
        mainland_n = M,
        pars = pars,
        island_type = island_type,
        nonoceanic = nonoceanic,
        island_ontogeny = island_ontogeny,
        sea_level = sea_level,
        Apars = Apars,
        Epars = Epars,
        keep_final_state = keep_final_state,
        island_spec = NULL
      )

      if (verbose == TRUE) {
        print(paste("Island replicate ", rep, sep = ""))
      }
    }


    if(island_type == "oceanic")
    {
      island_replicates = DAISIE_format_IW_oceanic(island_replicates = island_replicates,
                                                   time = totaltime,
                                                   M = M,
                                                   sample_freq = sample_freq)
    }
    if(island_type == "nonoceanic")
    {
      island_replicates = DAISIE_format_IW_nonoceanic(island_replicates = island_replicates,
                                                      time = totaltime,
                                                      M = M,
                                                      sample_freq = sample_freq)
    }
  }

  if(divdepmodel == 'CS')
  {
    if(length(pars) == 5)
    {
      # Midway simulation

      if (!is.null(stored_data)) {
        for (rep in 1:replicates)
        {
          n_colonized_replicates <- length(stored_data[[rep]]) - 1
          colonized_island_spec <- list()

          for (k in 1:n_colonized_replicates) {
            colonized_island_spec[[k]] <- stored_data[[rep]][[k + 1]]$island_spec
          }

          island_replicates = list()

          # Run each clade seperately

          full_list = list()

          # Run midway clades
          #currently only run for oceanic DAISIE

          for (m_spec in 1:n_colonized_replicates)
          {
            full_list[[m_spec]] <- DAISIE_sim_core(
              time = totaltime,
              mainland_n = 1,
              pars = pars,
              island_ontogeny = island_ontogeny,
              Apars = Apars,
              Epars = Epars,
              keep_final_state = keep_final_state,
              island_spec = colonized_island_spec[[m_spec]]
            )
          }

          # Run empty clades that didn't get colonists

          for (m_spec in (n_colonized_replicates + 1):1000)
          {
            full_list[[m_spec]] <- DAISIE_sim_core(
              time = totaltime,
              mainland_n = 1,
              pars = pars,
              island_ontogeny = island_ontogeny,
              Apars = Apars,
              Epars = Epars,
              keep_final_state = keep_final_state,
              island_spec = NULL
            )
          }

          island_replicates[[rep]] = full_list

          if (verbose == TRUE) {
            print(paste("Island replicate ",rep,sep = ""))
          }
        }

      } else {

        # Simulation from empty island

        for(rep in 1:replicates)
        {
          island_replicates[[rep]] = list()

          # Run each clade seperately

          full_list = list()

          for(m_spec in 1:M)
          {
            full_list[[m_spec]] <- DAISIE_sim_core(
              time = totaltime,
              mainland_n = 1,
              pars = pars,
              island_type = island_type,
              nonoceanic = nonoceanic,
              island_ontogeny = island_ontogeny,
              sea_level = sea_level,
              Apars = Apars,
              Epars = Epars,
              keep_final_state = keep_final_state,
              island_spec = NULL
            )
          }
          island_replicates[[rep]] = full_list

          if (verbose == TRUE) {
            print(paste("Island replicate ",rep,sep = ""))
          }
        }
      }
    }

    if(length(pars) == 10)
    {
      if(is.na(prop_type2_pool))
      {
        stop('prop_type2_pool (fraction of mainland species that belongs to the second subset of species) must be specified when running model with two species types')
      }
      if(island_type == "nonoceanic")
      {
        stop('nonoceanic islands cannot have two type islands')
      }
      if(replicates_apply_type2 == TRUE)
      {
        island_replicates = DAISIE_sim_min_type2(time = totaltime,
                                                 M = M,
                                                 pars = pars,
                                                 replicates = replicates,
                                                 prop_type2_pool = prop_type2_pool)
      } else
      {
        for(rep in 1:replicates)
        {
          pool2 = DDD::roundn(M * prop_type2_pool)
          pool1 = M - pool2
          lac_1 = pars[1]
          mu_1 = pars[2]
          K_1 = pars[3]
          gam_1 = pars[4]
          laa_1 = pars[5]
          lac_2 = pars[6]
          mu_2 = pars[7]
          K_2 = pars[8]
          gam_2 = pars[9]
          laa_2 = pars[10]
          full_list = list()

          #### species of pool1

          for(m_spec in 1:pool1)
          {
            full_list[[m_spec]] = DAISIE_sim_core(time = totaltime,mainland_n = 1,pars = c(lac_1,mu_1,K_1,gam_1,laa_1))
            full_list[[m_spec]]$type1or2  = 1
          }

          #### species of pool2

          for(m_spec in (pool1 + 1):(pool1 + pool2))
          {
            full_list[[m_spec]] = DAISIE_sim_core(time = totaltime,
                                                          mainland_n = 1,
                                                          pars = c(lac_2,mu_2,K_2,gam_2,laa_2))
            full_list[[m_spec]]$type1or2 = 2
          }

          island_replicates[[rep]] = full_list

          if (verbose == TRUE) {
            print(paste("Island replicate ",rep,sep = ""))
          }
        }
      }
    }

    if (start_midway == TRUE)
    {
      island_replicates <- DAISIE_format_CS_oceanic(
        island_replicates = island_replicates,
        time = totaltime,
        M = M,
        sample_freq = sample_freq,
        start_midway = start_midway,
        verbose = verbose
      )
    }


    if(island_type == "oceanic")
    {
      island_replicates = DAISIE_format_CS_oceanic(island_replicates = island_replicates,
                                                   time = totaltime,
                                                   M = M,
                                                   sample_freq = sample_freq)
    }
    if(island_type == "nonoceanic")
    {
      island_replicates = DAISIE_format_CS_nonoceanic(island_replicates = island_replicates,
                                                      time = totaltime,
                                                      M = M,
                                                      sample_freq = sample_freq)
    }

  }

  if(plot_sims == TRUE)
  {
    DAISIE_plot_sims(island_replicates)
  }
  return(island_replicates)
}