R/creating.trait.R

In MoBPS: Modular Breeding Program Simulator

'#
  Authors
Torsten Pook, torsten.pook@wur.nl

Copyright (C) 2017 -- 2025  Torsten Pook

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 3
of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
'#

#' Generation of genomic traits
#'
#' Generation of the trait in a starting population
#' @param population Population list
#' @param trait.cor Target correlation between QTL-based traits (underlying true genomic values)
#' @param trait.cor.include Vector of traits to be included in the modelling of correlated traits (default: all - needs to match with trait.cor)
#' @param qtl.position.shared Set to TRUE to put QTL effects on the same markers for different traits
#' @param n.traits Number of traits (If more than traits via real.bv.X use traits with no directly underlying QTL)
#' @param polygenic.variance Genetic variance of traits with no underlying QTL
#' @param randomSeed Set random seed of the process
#' @param n.additive Number of additive QTL with effect size drawn from a gaussian distribution
#' @param n.dominant Number of dominant QTL with effect size drawn from a gaussian distribution
#' @param n.overdominant Number of overdominant QTL with effect size drawn from absolute value of a gaussian distribution
#' @param n.equal.additive Number of additive QTL with equal effect size (effect.size)
#' @param n.equal.dominant Number of dominant QTL with equal effect size
#' @param n.equal.overdominant Number of overdominant QTL with equal effect size
#' @param n.qualitative Number of qualitative epistatic QTL
#' @param n.quantitative Number of quantitative epistatic QTL
#' @param effect.distribution Set to "gamma" for gamma distribution effects with gamma.shape1, gamma.shape2 instead of gaussian (default: "gauss")
#' @param gamma.shape1 Default: 1
#' @param gamma.shape2 Default: 1
#' @param dominant.only.positive Set to TRUE to always assign the heterozygous variant with the higher of the two homozygous effects (e.g. hybrid breeding); default: FALSE
#' @param var.additive.l Variance of additive QTL
#' @param var.dominant.l Variance of dominante QTL
#' @param var.overdominant.l Variance of overdominante QTL
#' @param var.qualitative.l Variance of qualitative epistatic QTL
#' @param var.quantitative.l Variance of quantitative epistatic QTL
#' @param effect.size.equal.add Effect size of the QTLs in n.equal.additive
#' @param effect.size.equal.dom Effect size of the QTLs in n.equal.dominant
#' @param effect.size.equal.over Effect size of the QTLs in n.equal.overdominant
#' @param exclude.snps Marker were no QTL are simulated on
#' @param replace.traits If TRUE delete the simulated traits added before
#' @param shuffle.cor OLD! Use trait.cor - Target Correlation between traits
#' @param shuffle.traits OLD! Use trait.cor.include - Vector of traits to be included for modelling of correlated traits (default: all - needs to match with shuffle.cor)
#' @param real.bv.add Single Marker effects
#' @param real.bv.mult Two Marker effects
#' @param real.bv.dice Multi-marker effects
#' @param bve.mult.factor Multiplicate trait value times this
#' @param bve.poly.factor Potency trait value over this
#' @param base.bv Average genetic value of a trait
#' @param new.phenotype.correlation (OLD! - use new.residual.correlation) Correlation of the simulated enviromental variance
#' @param new.residual.correlation Correlation of the simulated enviromental variance
#' @param new.breeding.correlation Correlation of the simulated genetic variance (child share! heritage is not influenced!
#' @param trait.name Name of the trait generated
#' @param remove.invalid.qtl Set to FALSE to deactive the automatic removal of QTLs on markers that do not exist
#' @param bv.standard Set TRUE to standardize trait mean and variance via bv.standardization()
#' @param mean.target Target mean
#' @param var.target Target variance
#' @param verbose Set to FALSE to not display any prints
#' @param is.maternal Vector coding if a trait is caused by a maternal effect (Default: all FALSE)
#' @param is.paternal Vector coding if a trait is caused by a paternal effect (Default: all FALSE)
#' @param fixed.effects Matrix containing fixed effects (p x k -matrix with p being the number of traits and k being number of fixed effects; default: p x 1 matrix with 0s (additional intercept))
#' @param trait.pool Vector providing information for which pools QTLs of this trait are activ (default: 0 - all pools)
#' @param set.zero Set to TRUE to have no effect on the 0 genotype (or 00 for QTLs with 2 underlying SNPs)
#' @param gxe.correlation Correlation matrix between locations / environments (default: only one location, sampled from gxe.max / gxe.min)
#' @param gxe.max Maximum correlation between locations / environments when generating correlation matrix via sampling (default: 0.85)
#' @param gxe.min Minimum correlation between locations / environments when generating correlation matrix via sampling (default: 0.70)
#' @param n.locations Number of locations / environments to consider for the GxE model
#' @param gxe.combine Set to FALSE to not view the same trait from different locations / environments as the sample trait in the prediction model (default: TRUE)
#' @param location.name Same of the different locations / environments used
#' @param use.recalculate.manual Set to TRUE to use recalculate.manual to calculate genomic values (all individuals and traits jointly, default: FALSE)
#' @param bv.total OLD! Use n.traits instead. Number of traits (If more than traits via real.bv.X use traits with no directly underlying QTL)
#' @examples
#' population <- creating.diploid(nsnp=1000, nindi=100)
#' population <- creating.trait(population, n.additive=100)
#' @return Population-list with one or more additional new traits
#' @export


creating.trait <- function(population,
                           ###### Traits
                           trait.name=NULL,
                           mean.target=NULL,
                           var.target=NULL,
                           qtl.position.shared = FALSE,
                           trait.cor = NULL,
                           trait.cor.include = NULL,
                           n.additive=0,
                           n.equal.additive=0,
                           n.dominant=0,
                           n.equal.dominant=0,
                           n.overdominant=0,
                           n.equal.overdominant=0,
                           n.qualitative=0,
                           n.quantitative=0,
                           effect.distribution = "gauss",
                           gamma.shape1 = 1,
                           gamma.shape2 = 1,
                           real.bv.add=NULL,
                           real.bv.mult=NULL,
                           real.bv.dice=NULL,
                           n.traits=0,
                           base.bv=NULL,
                           new.residual.correlation=NULL,
                           new.breeding.correlation=NULL,
                           is.maternal=NULL,
                           is.paternal=NULL,
                           fixed.effects=NULL,
                           trait.pool = 0,
                           gxe.correlation = NULL,
                           n.locations = NULL,
                           gxe.max = 0.85,
                           gxe.min = 0.7,
                           location.name = NULL,
                           gxe.combine = TRUE,
                           dominant.only.positive = FALSE,
                           exclude.snps=NULL,
                           var.additive.l=NULL,
                           var.dominant.l=NULL,
                           var.overdominant.l=NULL,
                           var.qualitative.l=NULL,
                           var.quantitative.l=NULL,
                           effect.size.equal.add = 1,
                           effect.size.equal.dom = 1,
                           effect.size.equal.over = 1,

                           polygenic.variance=100,
                           bve.mult.factor=NULL,
                           bve.poly.factor=NULL,
                           set.zero = FALSE,
                           bv.standard=FALSE,
                           replace.traits=FALSE,
                           remove.invalid.qtl=TRUE,
                           #### Other
                           randomSeed=NULL,
                           verbose=TRUE,
                           use.recalculate.manual = NULL,
                           new.phenotype.correlation=NULL,
                           shuffle.traits=NULL,
                           shuffle.cor= NULL,
                           bv.total = 0){


  {
  if(n.traits>0){
    if(bv.total>0){
      warning("bv.total has been overwritten with value from n.traits")
    }
    bv.total = n.traits
  }

  if(length(trait.cor)>0){
    shuffle.cor = trait.cor
  }
  if(length(trait.cor.include)>0){
    shuffle.traits = trait.cor.include
  }
  }
  # GxE Trait generation module
  {
    if(length(n.locations)>0 && n.locations > 1 && length(gxe.correlation)==0){
      gxe.correlation = matrix(stats::runif(n.locations^2, gxe.min, gxe.max), ncol=n.locations)
      for(i in 2:nrow(gxe.correlation)) {
        for(j in 1:(i-1)) {
          gxe.correlation[i,j]=gxe.correlation[j,i]
        }
      }
      diag(gxe.correlation) = 1
      gxe.correlation = matrix.posdef(A = gxe.correlation)

      if(verbose){
        cat("Generated GxE matrix")
        print(round(gxe.correlation, digits = 3))
      }

    }

    if(length(gxe.correlation)>1 && length(n.locations)==0){
      n.locations = ncol(gxe.correlation)
    }

    if(length(location.name)==0 && length(gxe.correlation)>0){
      location.name = paste0("Location ", 1:ncol(gxe.correlation))
    }


    trait_location = NULL
    trait_nr = NULL

    if(length(gxe.correlation)>0){

      if(length(population)>0 && population$info$bv.nr >0){
        stop("GxE module is only intended for the use when no traits where previously generated")
      }

      if(length(real.bv.add)>0 || length(real.bv.mult)>0 || length(real.bv.dice)>0){
        stop("GxE module is only intended for the use with predefined MoBPS trait architectures")
      }

      # Determine total number of traits

      trait_sum <- n.additive + n.dominant + n.qualitative + n.quantitative + n.equal.additive + n.equal.dominant + n.overdominant + n.equal.overdominant
      n_traits <- length(trait_sum)

      n.additive <- rep(c(n.additive, rep(0, length.out=n_traits-length(n.additive))), n.locations)
      n.dominant <- rep(c(n.dominant, rep(0, length.out=n_traits-length(n.dominant))), n.locations)
      n.equal.additive <- rep(c(n.equal.additive, rep(0, length.out=n_traits-length(n.equal.additive))), n.locations)
      n.equal.dominant <- rep(c(n.equal.dominant, rep(0, length.out=n_traits-length(n.equal.dominant))), n.locations)
      n.overdominant <- rep(c(n.overdominant, rep(0, length.out=n_traits-length(n.overdominant))), n.locations)
      n.equal.overdominant <- rep(c(n.equal.overdominant, rep(0, length.out=n_traits-length(n.equal.overdominant))), n.locations)

      n.qualitative <- rep(c(n.qualitative, rep(0, length.out=n_traits-length(n.qualitative))), n.locations)
      n.quantitative <- rep(c(n.quantitative, rep(0, length.out=n_traits-length(n.quantitative))), n.locations)

      if(length(trait.name) < n_traits){
        trait.name = c(trait.name, paste0("Trait ", (length(trait.name)+1):n_traits))
      }

      # GxE will always result in a multi-trait model

      if(length(shuffle.cor)==0){
        shuffle.cor = diag(1, n_traits)
        shuffle.traits = 1:n_traits
      }
      n.locations = ncol(gxe.correlation)
      if(length(shuffle.cor)>0){

        shuffle.cor =   gxe.correlation  %x% shuffle.cor
        if(length(shuffle.traits)>0){
          shuffle.traits = rep(shuffle.traits, n.locations) + sort(rep(1:n.locations*n_traits - n_traits, length(shuffle.traits)))
        }
      }

      if(length(trait.name) < (n_traits * n.locations)){
        trait.name = paste0(rep(trait.name, n.locations) ," x ", rep(location.name, each = n_traits))
      }


      trait_location = rep(1:n.locations, each = n_traits)
      trait_nr =  rep(1:n_traits, n.locations)

      colnames(shuffle.cor) = rownames(shuffle.cor) = trait.name
      if(verbose && n.locations > 1 && n_traits > 1){
        cat("Used genetic correlation matrix:\n")
        print(shuffle.cor)
      }

    }
  }

  {
  if(length(randomSeed)>0){
    set.seed(randomSeed)
  }
  if(length(mean.target)>0){
    bv.standard <- TRUE
  } else{
    mean.target <- 100
  }
  if(length(var.target)>0){
    bv.standard <- TRUE
  } else{
    var.target <- 10
  }

  if(sum(set.zero)>0){
    bv.standard <- TRUE
  }

  if(replace.traits==FALSE & length(population$info$bv.nr)>0){
    prior_traits = population$info$bv.nr
  } else{
    prior_traits = 0
  }

  preserve.bve <- length(population)==0

  if(length(new.phenotype.correlation)>0){
    new.residual.correlation <- new.phenotype.correlation
  }

  if(!is.list(var.additive.l) ){
    var.additive.l <- list(var.additive.l)
  }
  if(!is.list(var.dominant.l)){
    var.dominant.l <- list(var.dominant.l)
  }
  if(!is.list(var.overdominant.l)){
    var.overdominant.l <- list(var.overdominant.l)
  }
  if(!is.list(var.qualitative.l)){
    var.qualitative.l <- list(var.qualitative.l)
  }
  if(!is.list(var.quantitative.l)){
    var.quantitative.l <- list(var.quantitative.l)
  }

  trait_sum <- n.additive + n.dominant + n.qualitative + n.quantitative + n.equal.additive + n.equal.dominant + n.overdominant + n.equal.overdominant
  test <- list(NULL)

  if(length(var.additive.l) < length(trait_sum)){
    var.additive.l <- c(var.additive.l, rep(test,length.out=length(trait_sum)-length(var.additive.l)))
  }
  if(length(var.dominant.l) < length(trait_sum)){
    var.dominant.l <- c(var.dominant.l, rep(test,length.out=length(trait_sum)-length(var.dominant.l)))
  }

  if(length(var.overdominant.l) < length(trait_sum)){
    var.overdominant.l <- c(var.overdominant.l, rep(test,length.out=length(trait_sum)-length(var.overdominant.l)))
  }

  if(length(var.qualitative.l) < length(trait_sum)){
    var.qualitative.l <- c(var.qualitative.l, rep(test,length.out=length(trait_sum)-length(var.qualitative.l)))
  }
  if(length(var.quantitative.l) < length(trait_sum)){
    var.quantitative.l <- c(var.quantitative.l, rep(test,length.out=length(trait_sum)-length(var.quantitative.l)))
  }

  ntraits <- length(trait_sum)
  n.additive <- c(n.additive, rep(0, length.out=ntraits-length(n.additive)))
  n.dominant <- c(n.dominant, rep(0, length.out=ntraits-length(n.dominant)))
  n.equal.additive <- c(n.equal.additive, rep(0, length.out=ntraits-length(n.equal.additive)))
  n.equal.dominant <- c(n.equal.dominant, rep(0, length.out=ntraits-length(n.equal.dominant)))
  n.overdominant <- c(n.overdominant, rep(0, length.out=ntraits-length(n.overdominant)))
  n.equal.overdominant <- c(n.equal.overdominant, rep(0, length.out=ntraits-length(n.equal.overdominant)))

  n.qualitative <- c(n.qualitative, rep(0, length.out=ntraits-length(n.qualitative)))
  n.quantitative <- c(n.quantitative, rep(0, length.out=ntraits-length(n.quantitative)))

  if(length(unlist(c(var.qualitative.l, var.quantitative.l, var.additive.l, var.dominant.l, var.overdominant.l)))>0){
    ntraits <- max(length(trait_sum), length(var.additive.l),length(var.dominant.l), length(var.overdominant.l), length(var.qualitative.l), length(var.quantitative.l) )
    n.additive <- c(n.additive, rep(0, length.out=ntraits-length(n.additive)))
    n.dominant <- c(n.dominant, rep(0, length.out=ntraits-length(n.dominant)))
    n.equal.additive <- c(n.equal.additive, rep(0, length.out=ntraits-length(n.equal.additive)))
    n.equal.dominant <- c(n.equal.dominant, rep(0, length.out=ntraits-length(n.equal.dominant)))
    n.overdominant <- c(n.overdominant, rep(0, length.out=ntraits-length(n.overdominant)))
    n.equal.overdominant <- c(n.equal.overdominant, rep(0, length.out=ntraits-length(n.equal.overdominant)))

    n.qualitative <- c(n.qualitative, rep(0, length.out=ntraits-length(n.qualitative)))
    n.quantitative <- c(n.quantitative, rep(0, length.out=ntraits-length(n.quantitative)))
    trait_sum <- n.additive + n.dominant + n.qualitative + n.quantitative + n.equal.additive + n.equal.dominant + n.overdominant + n.equal.overdominant
    if(length(var.additive.l) < length(trait_sum)){
      var.additive.l <- rep(var.additive.l, length.out=length(trait_sum))
    }
    if(length(var.dominant.l) < length(trait_sum)){
      var.dominant.l <- rep(var.dominant.l, length.out=length(trait_sum))
    }
    if(length(var.overdominant.l) < length(trait_sum)){
      var.overdominant.l <- rep(var.overdominant.l, length.out=length(trait_sum))
    }
    if(length(var.qualitative.l) < length(trait_sum)){
      var.qualitative.l <- rep(var.qualitative.l, length.out=length(trait_sum))
    }
    if(length(var.quantitative.l) < length(trait_sum)){
      var.quantitative.l <- rep(var.quantitative.l, length.out=length(trait_sum))
    }
  }

  if(is.list(real.bv.dice) && length(unlist(real.bv.dice))==0){
    real.bv.dice = NULL
  }

  if(length(real.bv.dice)>0){
    if(is.list(real.bv.dice)){
      mdepth <- 0
      for(index in 1:length(real.bv.dice)){
        if(is.data.frame(real.bv.dice[[index]][[1]]) || is.matrix(real.bv.dice[[index]][[1]])){
          mdepth <- 1
        }
      }
      if(mdepth==0){
        for(index in 1:length(real.bv.dice)){
          if(length(real.bv.dice[[index]])>0){
            for(index2 in 1:length(real.bv.dice[[index]])){
              if(is.data.frame(real.bv.dice[[index]][[1]][[index2]]) || is.matrix(real.bv.dice[[index]][[1]][[index2]]) ){
                mdepth <- 2
              }
            }
          }
        }
      }
    }
    if(mdepth == 1){
      real.bv.dice <- list(real.bv.dice)
    }
    if(mdepth == 0){
      stop("Illegal input for real.bv.dice")
    }
    for(index in 1:length(real.bv.dice)){
      if(length(real.bv.dice[[index]])>0){
        for(index2 in 1:length(real.bv.dice[[index]][[1]])){
          if(length(real.bv.dice[[index]][[2]][[index2]]) != (3 ^nrow(real.bv.dice[[index]][[1]][[index2]]))){
            warning("")
            stop("Length of effects does not match with involved effect SNPs - should be 3^(effect SNPs) (0..0, 0..01, ..., 2..2)")
          }
        }
      }
    }
  }

  if(length(real.bv.add)>0){
    if(!is.list(real.bv.add)){
      real.bv.add <- list(real.bv.add)
    }

    if(length(trait.pool)< max(length(real.bv.add),length(trait_sum))){
      trait.pool_temp = rep(trait.pool, length.out = length(real.bv.add))
    } else{
      trait.pool_temp = trait.pool
    }

    for(index3 in 1:length(real.bv.add)){

      if(length(real.bv.add[[index3]])==1 && real.bv.add[[index3]]=="placeholder"){
        real.bv.add[[index3]] = NULL
        next
      }
      if(ncol(real.bv.add[[index3]])==5){
        real.bv.add[[index3]] = cbind(real.bv.add[[index3]], NA, trait.pool_temp[index3])
      }

      if(ncol(real.bv.add[[index3]])==6){
        real.bv.add[[index3]] = cbind(real.bv.add[[index3]], trait.pool_temp[index3])
      }

      if(ncol(real.bv.add[[index3]])==7){
        real.bv.add[[index3]] = cbind(real.bv.add[[index3]], FALSE)
      }
    }
  }

}

  if(length(population)>0){
    if(length(real.bv.add)==0 && replace.traits==FALSE){
      real.bv.add <- population$info$real.bv.add
      if(length(population$info$real.bv.add)>0){
        real.bv.add[[length(population$info$real.bv.add)]] <- NULL
      }
    } else if(replace.traits==FALSE){
      if(!is.list(real.bv.add)){
        real.bv.add <- list(real.bv.add)
      }
      real.bv.add <- c(population$info$real.bv.add, real.bv.add)
      if(length(population$info$real.bv.add)>0){
        real.bv.add[[length(population$info$real.bv.add)]] <- NULL
      }
    }
    if(length(real.bv.mult)==0 && replace.traits==FALSE){
      real.bv.mult <- population$info$real.bv.mult
      if(length(population$info$real.bv.mult)>0){
        real.bv.mult[[length(population$info$real.bv.mult)]] <- NULL
      }
    } else if(replace.traits==FALSE){
      if(!is.list(real.bv.mult)){
        real.bv.mult <- list(real.bv.mult)
      }
      real.bv.mult <- c(population$info$real.bv.mult, real.bv.mult)
      if(length(population$info$real.bv.mult)>0){
        real.bv.mult[[length(population$info$real.bv.mult)]] <- NULL
      }
    }
    if(length(real.bv.dice)==0 && replace.traits==FALSE){
      real.bv.dice <- population$info$real.bv.dice
      if(length(population$info$real.bv.dice)>0){
        real.bv.dice[[length(population$info$real.bv.dice)]] <- NULL
      }
    } else if(replace.traits==FALSE){
      if(!is.list(real.bv.dice)){
        real.bv.dice <- list(real.bv.dice)
      }
      real.bv.dice <- c(population$info$real.bv.dice, real.bv.dice)
      if(length(population$info$real.bv.dice)>0){
        real.bv.dice[[length(population$info$real.bv.dice)]] <- NULL
      }
    }

  }
  if(length(real.bv.add)>0 && !is.list(real.bv.add)){
    real.bv.add <- list(real.bv.add)
  }
  if(length(real.bv.mult)>0 && !is.list(real.bv.mult)){
    real.bv.mult <- list(real.bv.mult)
  }
  if(length(real.bv.dice)>0 && !is.list(real.bv.dice)){
    real.bv.dice <- list(real.bv.dice)
  }

  {
    # Check for missingness in real.bvs and replace with reasonable inputs
    cum_snp <- cumsum(population$info$snp)
    snpdata <- snpdata <- population$info$snp

    if(qtl.position.shared){

      n_qtls = max(n.additive + n.dominant + n.equal.additive+ n.equal.dominant+ n.quantitative*2+ n.qualitative*2+ n.overdominant+ n.equal.overdominant)

      so_far <- max(length(real.bv.dice), length(real.bv.add), length(real.bv.mult))
      n_qtls_sofar = numeric(so_far)
      if(length(real.bv.add)>0){
        for(index in 1:length(real.bv.add)){
          if(length(real.bv.add[[index]])>0){
            n_qtls_sofar[index] = n_qtls_sofar[index] + nrow(real.bv.add[[index]])

          }
        }
      }
      if(length(real.bv.mult)>0){
        for(index in 1:length(real.bv.mult)){
          if(length(real.bv.mult[[index]])>0){
            n_qtls_sofar[index] = n_qtls_sofar[index] + nrow(real.bv.mult[[index]]) * 2
          }
        }
      }

      effect_marker <- (1:sum(snpdata))
      if(length(exclude.snps)>0){
        effect_marker <- effect_marker[-exclude.snps]
      }

      max_qtls = max(c(n_qtls, n_qtls_sofar))
      effect_marker = sample(effect_marker, max_qtls)

    } else{
      so_far <- max(length(real.bv.dice), length(real.bv.add), length(real.bv.mult))
      effect_marker <- (1:sum(snpdata))
      if(length(exclude.snps)>0){
        effect_marker <- effect_marker[-exclude.snps]
      }
    }

    if(length(real.bv.add)>0){
      for(index in 1:length(real.bv.add)){
        while(sum(is.na(real.bv.add[[index]][,c(1:2,6)]))>0){

          add_marker <- sample(effect_marker, nrow(real.bv.add[[index]]), replace=if(nrow(real.bv.add[[index]])>length(effect_marker)){TRUE} else{FALSE})
          add_snp <- real.bv.add[[index]][,1]
          add_chromo <- real.bv.add[[index]][,2]

          for(index2 in (1:nrow(real.bv.add[[index]]))[is.na(add_snp) | is.na(add_chromo)]){
            add_chromo[index2] <- sum(add_marker[index2] > cum_snp) + 1
            add_snp[index2] <- add_marker[index2] - c(0,cum_snp)[add_chromo[index2]]
          }


          if(sum(is.na(real.bv.add[[index]][,6]))>0){
            add_marker = add_snp + c(0,cum_snp)[add_chromo]
          }

          enter <- add_chromo==real.bv.add[[index]][,2] | is.na(real.bv.add[[index]][,2])| is.na(real.bv.add[[index]][,1]) | is.na(real.bv.add[[index]][,6])

          real.bv.add[[index]][enter,c(1:2, 6)] <- cbind(add_snp, add_chromo, add_marker)[enter,]
        }
      }
    }

    if(length(real.bv.mult)>0){
      for(index in 1:length(real.bv.mult)){
        for(columns in c(0,2)){
          while(sum(is.na(real.bv.mult[[index]][,1:2+columns]))>0){

            add_marker <- sample(effect_marker, nrow(real.bv.mult[[index]]), replace=if(nrow(real.bv.mult[[index]])>length(effect_marker)){TRUE} else{FALSE})
            add_snp <- real.bv.mult[[index]][,1]
            add_chromo <- real.bv.mult[[index]][,2]

            for(index2 in (1:nrow(real.bv.mult[[index]]))[is.na(add_snp) | is.na(add_chromo)]){
              add_chromo[index2] <- sum(add_marker[index2] > cum_snp) + 1
              add_snp[index2] <- add_marker[index2] - c(0,cum_snp)[add_chromo[index2]]
            }

            enter <- add_chromo==real.bv.mult[[index]][,2+columns] | is.na(real.bv.mult[[index]][,2])

            real.bv.mult[[index]][enter,1:2+columns] <- cbind(add_snp, add_chromo)[enter,]
          }
        }

      }
    }
  }

  if(qtl.position.shared){
    sofar_positions = NULL

    if(length(real.bv.add)>0){
      for(index in 1:length(real.bv.add)){
        if(length(real.bv.add[[index]])>0){
          sofar_positions = c(sofar_positions, real.bv.add[[index]][,6])
        }
      }
    }

    sofar_positions = sofar_positions[!is.na(sofar_positions)]
    sofar_positions = setdiff(sofar_positions, effect_marker)
    if(length(sofar_positions)>0 && length(sofar_positions)<= length(effect_marker)){
      effect_marker[1:length(sofar_positions)] = sofar_positions
    }
  }


  if(length(dominant.only.positive)<length(trait_sum)){
    dominant.only.positive <- rep(dominant.only.positive, length.out = length(trait_sum))
  }
  so_far <- max(length(real.bv.dice), length(real.bv.add), length(real.bv.mult))

  if(length(trait.pool)< length(trait_sum)){
    trait.pool = rep(trait.pool, length.out = length(trait_sum))
  }

  if(length(trait_sum)>0){
    for(index_trait in 1:length(trait_sum)){
      var_additive <- var.additive.l[[index_trait]]
      var_dominant <- var.dominant.l[[index_trait]]
      var_overdominant <- var.overdominant.l[[index_trait]]

      var_qualitative <- var.qualitative.l[[index_trait]]
      var_quantitative <- var.quantitative.l[[index_trait]]
      if(n.additive[index_trait]>0 && length(var_additive)<n.additive[index_trait]){
        if(length(var_additive)==0){
          var_additive <- 1
        }
        var_additive <- rep(var_additive, length.out=n.additive[index_trait])
        var.additive.l[[index_trait]] <- var_additive
      }
      if(n.dominant[index_trait]>0 && length(var_dominant)<n.dominant[index_trait]){
        if(length(var_dominant)==0){
          var_dominant <- 1
        }
        var_dominant <- rep(var_dominant, length.out=n.dominant[index_trait])
        var.dominant.l[[index_trait]] <- var_dominant
      }

      if(n.overdominant[index_trait]>0 && length(var_overdominant)<n.overdominant[index_trait]){
        if(length(var_overdominant)==0){
          var_overdominant <- 1
        }
        var_overdominant <- rep(var_overdominant, length.out=n.overdominant[index_trait])
        var.overdominant.l[[index_trait]] <- var_overdominant
      }

      if(n.qualitative[index_trait]>0 && length(var_qualitative)<n.qualitative[index_trait]){
        if(length(var_qualitative)==0){
          var_qualitative <- 1
        }
        var_qualitative <- rep(var_qualitative, length.out=n.qualitative[index_trait])
        var.qualitative.l[[index_trait]] <- var_qualitative
      }
      if(n.quantitative[index_trait]>0 && length(var_quantitative)<n.quantitative[index_trait]){
        if(length(var_quantitative)==0){
          var_quantitative <- 1
        }
        var_quantitative <- rep(var_quantitative, length.out=n.quantitative[index_trait])
        var.quantitative.l[[index_trait]] <- var_quantitative

      }

      if(length(var_additive)!= n.additive[index_trait]){
        n.additive[index_trait] <- length(var_additive)
      }
      if(length(var_dominant)!= n.dominant[index_trait]){
        n.dominant[index_trait] <- length(var_dominant)
      }
      if(length(var_overdominant)!= n.overdominant[index_trait]){
        n.overdominant[index_trait] <- length(var_overdominant)
      }
      if(length(var_qualitative)!= n.qualitative[index_trait]){
        n.qualitative[index_trait] <- length(var_qualitative)
      }
      if(length(var_quantitative)!= n.quantitative[index_trait]){
        n.quantitative[index_trait] <- length(var_quantitative)
      }






      #stop()

      # Generating additive

      add_marker <- sample(effect_marker, n.additive[index_trait], replace=if(n.additive[index_trait]>length(effect_marker)){TRUE} else{FALSE})
      dom_marker <- sample(effect_marker, n.dominant[index_trait], replace=if(n.dominant[index_trait]>length(effect_marker)){TRUE} else{FALSE})
      over_marker <- sample(effect_marker, n.overdominant[index_trait], replace=if(n.overdominant[index_trait]>length(effect_marker)){TRUE} else{FALSE})

      add_marker1 <- sample(effect_marker, n.equal.additive[index_trait], replace=if(n.equal.additive[index_trait]>length(effect_marker)){TRUE} else{FALSE})
      dom_marker1 <- sample(effect_marker, n.equal.dominant[index_trait], replace=if(n.equal.dominant[index_trait]>length(effect_marker)){TRUE} else{FALSE})
      over_marker1 <- sample(effect_marker, n.equal.overdominant[index_trait], replace=if(n.equal.overdominant[index_trait]>length(effect_marker)){TRUE} else{FALSE})

      epi1_marker <- sample(effect_marker, n.quantitative[index_trait]*2, replace=if(n.quantitative[index_trait]*2>length(effect_marker)){TRUE} else{FALSE})
      epi2_marker <- sample(effect_marker, n.qualitative[index_trait]*2, replace=if(n.qualitative[index_trait]*2>length(effect_marker)){TRUE} else{FALSE})



      cum_snp <- cumsum(snpdata)
      real.bv.add.new <- NULL
      real.bv.mult.new <- NULL
      if(n.additive[index_trait]>0){
        add_snp <- add_chromo <- numeric(n.additive[index_trait])
        for(index in 1:n.additive[index_trait]){
          add_chromo[index] <- sum(add_marker[index] > cum_snp) + 1
          add_snp[index] <- add_marker[index] - c(0,cum_snp)[add_chromo[index]]
        }
        if(effect.distribution == "gauss"){
          add_effect <- stats::rnorm(n.additive[index_trait], 0, var_additive)
        } else{
          add_effect <- stats::rgamma(n.additive[index_trait], gamma.shape1, gamma.shape2) * sample( c(-1,1), n.additive[index_trait], replace = TRUE)
        }

        real.bv.add.new <- cbind(add_snp, add_chromo, add_effect,0,-add_effect, add_marker, trait.pool[index_trait], FALSE)
      }

      if(n.equal.additive[index_trait]>0){
        add_snp1 <- add_chromo1 <- numeric(n.equal.additive[index_trait])
        for(index in 1:n.equal.additive[index_trait]){
          add_chromo1[index] <- sum(add_marker1[index] > cum_snp) + 1
          add_snp1[index] <- add_marker1[index] - c(0,cum_snp)[add_chromo1[index]]
        }
        add_effect1 <- effect.size.equal.add
        real.bv.add.new <- rbind(real.bv.add.new, cbind(add_snp1, add_chromo1,  -add_effect1, 0, add_effect1, add_marker1, trait.pool[index_trait], FALSE))

      }

      if(n.dominant[index_trait]>0){
        dom_snp <- dom_chromo <- numeric(n.dominant[index_trait])
        for(index in 1:n.dominant[index_trait]){
          dom_chromo[index] <- sum(dom_marker[index] > cum_snp) + 1
          dom_snp[index] <- dom_marker[index] - c(0,cum_snp)[dom_chromo[index]]
        }

        if(effect.distribution == "gauss"){
          dom_effect <- stats::rnorm(n.dominant[index_trait], 0, var_dominant)
        } else{
          dom_effect <- stats::rgamma(n.dominant[index_trait], gamma.shape1, gamma.shape2) * sample( c(-1,1), n.dominant[index_trait], replace = TRUE)
        }

        if(dominant.only.positive[index_trait]){
          temp1 <- dom_effect
          temp1[temp1<0] <- 0
        } else{
          temp1 <- dom_effect
        }
        real.bv.add.new <- rbind(real.bv.add.new, cbind(dom_snp, dom_chromo, 0 ,temp1,dom_effect, dom_marker, trait.pool[index_trait], FALSE))

      }

      if(n.equal.dominant[index_trait]>0){
        dom_snp1 <- dom_chromo1 <- numeric(n.equal.dominant[index_trait])
        for(index in 1:n.equal.dominant[index_trait]){
          dom_chromo1[index] <- sum(dom_marker1[index] > cum_snp) + 1
          dom_snp1[index] <- dom_marker1[index] - c(0,cum_snp)[dom_chromo1[index]]
        }
        dom_effect1 <- effect.size.equal.dom
        real.bv.add.new <- rbind(real.bv.add.new, cbind(dom_snp1, dom_chromo1, 0 ,dom_effect1, dom_effect1, dom_marker1, trait.pool[index_trait], FALSE))

      }

      if(n.overdominant[index_trait]>0){
        over_snp <- over_chromo <- numeric(n.overdominant[index_trait])
        for(index in 1:n.overdominant[index_trait]){
          over_chromo[index] <- sum(over_marker[index] > cum_snp) + 1
          over_snp[index] <- over_marker[index] - c(0,cum_snp)[over_chromo[index]]
        }

        if(effect.distribution == "gauss"){
          over_effect <- abs(stats::rnorm(n.overdominant[index_trait], 0, var_overdominant))
        } else{
          over_effect <- abs(stats::rgamma(n.overdominant[index_trait], gamma.shape1, gamma.shape2) * sample( c(-1,1), n.overdominant[index_trait], replace = TRUE))
        }

        temp1 <- over_effect

        real.bv.add.new <- rbind(real.bv.add.new, cbind(over_snp, over_chromo, 0 ,temp1,0, over_marker, trait.pool[index_trait], FALSE))

      }

      if(n.equal.overdominant[index_trait]>0){
        over_snp1 <- over_chromo1 <- numeric(n.equal.overdominant[index_trait])
        for(index in 1:n.equal.overdominant[index_trait]){
          over_chromo1[index] <- sum(over_marker1[index] > cum_snp) + 1
          over_snp1[index] <- over_marker1[index] - c(0,cum_snp)[over_chromo1[index]]
        }
        over_effect1 <- effect.size.equal.over
        real.bv.add.new <- rbind(real.bv.add.new, cbind(dom_snp1, dom_chromo1, 0 ,over_effect1, 0, over_marker1, trait.pool[index_trait], FALSE))

      }

      if(n.quantitative[index_trait]){
        epi1_snp <- epi1_chromo <- numeric(n.quantitative[index_trait]*2)
        for(index in 1:(n.quantitative[index_trait]*2)){
          epi1_chromo[index] <- sum(epi1_marker[index] > cum_snp) + 1
          epi1_snp[index] <- epi1_marker[index] - c(0,cum_snp)[epi1_chromo[index]]
        }

        effect_matrix <- matrix(0,nrow=n.quantitative[index_trait], ncol=9)
        for(index in 1:n.quantitative[index_trait]){


          if(effect.distribution == "gauss"){
            d1 <- sort(abs(stats::rnorm(3, 0, var_quantitative[index])))
            d2 <- sort(abs(stats::rnorm(3, 0, var_quantitative[index])))
          } else{
            d1 <- sort(stats::rgamma(3, gamma.shape1, gamma.shape2))
            d2 <- sort(stats::rgamma(3, gamma.shape1, gamma.shape2))
          }

          effect_matrix[index,] <- c(d1*d2[1], d1*d2[2], d1*d2[3])
        }
        real.bv.mult.new <- cbind(epi1_snp[1:n.quantitative[index_trait]], epi1_chromo[1:n.quantitative[index_trait]],
                                  epi1_snp[-(1:n.quantitative[index_trait])], epi1_chromo[-(1:n.quantitative[index_trait])],
                                  effect_matrix)
      }

      if(n.qualitative[index_trait]>0){
        epi2_snp <- epi2_chromo <- numeric(n.qualitative[index_trait]*2)
        for(index in 1:(n.qualitative[index_trait]*2)){
          epi2_chromo[index] <- sum(epi2_marker[index] > cum_snp) + 1
          epi2_snp[index] <- epi2_marker[index] - c(0,cum_snp)[epi2_chromo[index]]
        }

        effect_matrix <- matrix(0,nrow=n.qualitative[index_trait], ncol=9)
        for(index in 1:n.qualitative[index_trait]){

          if(effect.distribution == "gauss"){
            d1 <- -abs(stats::rnorm(9, 0, var_qualitative[index]))
          } else{
            d1 <- - stats::rgamma(9, gamma.shape1, gamma.shape2)
          }

          d1[c(3,7)] <- -d1[c(3,7)]
          effect_matrix[index,] <- d1
        }
        real.bv.mult.new <- rbind(real.bv.mult.new, cbind(epi2_snp[1:n.qualitative[index_trait]], epi2_chromo[1:n.qualitative[index_trait]],
                                                          epi2_snp[-(1:n.qualitative[index_trait])], epi2_chromo[-(1:n.qualitative[index_trait])],
                                                          effect_matrix))
      }

      real.bv.add[[index_trait+so_far]] <- real.bv.add.new
      real.bv.mult[[index_trait+so_far]] <- real.bv.mult.new

    }
  }





  if(length(real.bv.add)>0 && !is.list(real.bv.add)){
    real.bv.add <- list(real.bv.add)
  }
  if(length(real.bv.mult)>0 && !is.list(real.bv.mult)){
    real.bv.mult <- list(real.bv.mult)
  }
  if(length(real.bv.dice)>0 && !is.list(real.bv.dice)){
    real.bv.dice <- list(real.bv.dice)
  }

  nbv <- max(length(real.bv.add), length(real.bv.mult), length(real.bv.dice), if(length(trait_sum)>1){length(trait_sum)} else{0})
  if(nbv >= bv.total){
    bv.total <- nbv
    bv.calc <- nbv
    bv.random <- rep(FALSE, bv.total)
    bv.random.variance <- c(rep(0, nbv))
  }
  if(bv.total > nbv){
    if(length(polygenic.variance)< (bv.total - nbv)){
      polygenic.variance <- rep(polygenic.variance, bv.total - nbv)
    }
    bv.random <- c(rep(FALSE, nbv), rep(TRUE, bv.total - nbv))

    bv.random.variance <- c(rep(0, nbv), polygenic.variance)
    bv.calc <- nbv +1
  }

  population$info$bve <- FALSE
  population$info$bv.calculated <- FALSE
  if(!replace.traits && prior_traits>0){
    population$info$bv.calculated.partly <- 1:prior_traits
  } else{
    population$info$bv.calculated.partly <- NULL

  }

  population$info$breeding.totals <- list()
  population$info$bve.data <- list()
  population$info$bv.nr <- 1 # default um fallunterscheidung zu vermeiden
  population$info$bv.random <- bv.random
  population$info$bv.random.variance <- bv.random.variance

  population$info$phenotypic.transform <- rep(FALSE, bv.total)
  population$info$phenotypic.transform.function <- list()



  store1 <- population$info$is.maternal
  store2 <- population$info$is.paternal
  store3 <- population$info$is.combi
  store4 <- population$info$phenotypic.transform
  store5 <- population$info$phenotypic.transform.function
  store6 <- population$info$bv.random
  store7 <- population$info$bv.random.variance
  store8 <- population$info$bve.mult.factor
  store9 <- population$info$bve.poly.factor
  store10 <- population$info$base.bv

  if(length(is.maternal)==0){
    population$info$is.maternal <- rep(FALSE, bv.total)
  } else{
    if(length(is.maternal)==bv.total){
      population$info$is.maternal <- is.maternal
    } else {
      population$info$is.maternal <- c(population$info$is.materal, rep(is.maternal, length.out = bv.total - length(population$info$is.materal)))
    }

  }
  if(length(is.paternal)==0){
    population$info$is.paternal <- rep(FALSE, bv.total)
  } else{
    if(length(is.paternal)==bv.total){
      population$info$is.paternal <- is.paternal
    } else {
      population$info$is.paternal <- c(population$info$is.paternal, rep(is.paternal, length.out = bv.total - length(population$info$is.paternal)))
    }
  }



  population$info$is.combi <- rep(FALSE, bv.total)


  if(length(bve.mult.factor)==0){
    population$info$bve.mult.factor <- rep(1L, bv.total)
  } else{

    if(length(bve.mult.factor)==bv.total){
      population$info$bve.mult.factor <- bve.mult.factor
    } else {
      population$info$bve.mult.factor <- c(population$info$bve.mult.factor, rep(bve.mult.factor, length.out = bv.total - length(population$info$bve.mult.factor)))
    }
  }

  if(length(bve.poly.factor)==0){
    population$info$bve.poly.factor <- rep(1L, bv.total)
  } else{

    if(length(bve.poly.factor)==bv.total){
      population$info$bve.poly.factor <- bve.poly.factor
    } else {
      population$info$bve.poly.factor <- c(population$info$bve.poly.factor, rep(bve.poly.factor, length.out = bv.total - length(population$info$bve.poly.factor)))
    }

  }
  if(length(base.bv)==0){
    population$info$base.bv <- rep(100L, bv.total)
  } else{

    if(length(base.bv)==bv.total){
      population$info$base.bv <- base.bv
    } else {
      population$info$base.bv <- c(population$info$base.bv, rep(base.bv, length.out = bv.total - length(population$info$base.bv)))
    }
  }


  if(bv.total>0 ||length(real.bv.add)>0  || length(real.bv.mult) >0 || length(real.bv.dice)>0){
    population$info$bve <- TRUE
    if(is.list(real.bv.add)){
      population$info$real.bv.add <- real.bv.add
    } else{
      population$info$real.bv.add <- list(real.bv.add)
    }
    if(is.list(real.bv.mult)){
      population$info$real.bv.mult <- real.bv.mult
    } else{
      population$info$real.bv.mult <- list(real.bv.mult)
    }
    if(is.list(real.bv.dice)){
      population$info$real.bv.dice <- real.bv.dice
    } else{
      if(length(real.bv.dice)>0){
        warning("Invalid input for real.bv.dice!")
      }
      population$info$real.bv.dice <- list(real.bv.dice)
    }

    population$info$bv.nr <- bv.total
    population$info$bv.calc <- bv.calc

    population$info$real.bv.length <- c(length(population$info$real.bv.add),
                                         length(population$info$real.bv.mult),
                                         length(population$info$real.bv.dice))

    population$info$real.bv.add[[nbv+1]] <- "placeholder"
    population$info$real.bv.mult[[nbv+1]] <- "placeholder"
    population$info$real.bv.dice[[nbv+1]] <- "placeholder"



  } else if(preserve.bve){
    population$info$bve <- FALSE
    population$info$bv.nr <- 0
    population$info$bv.calc <- 0
    population$info$real.bv.length <- c(0,0,0)
  }



  if(length(new.residual.correlation)==0 &&
     length(population$info$pheno.correlation)>0 &&
     sum(population$info$pheno.correlation)>sum(diag(population$info$pheno.correlation))){
    if(verbose) cat("Residual correlation has been set to zero since new traits were added ")
  }

  if(length(new.breeding.correlation)==0 &&
     length(population$info$bv.correlation)>0 &&
     sum(abs(population$info$bv.correlation))>sum(diag(population$info$bv.correlation))&&
     sum(population$info$is.combi | !population$info$bv.random) < population$info$bv.nr){
    if(verbose) cat("Genetic correlation between non-QTL traits has been set to zero since new traits were added ")
  }

  store11 <- population$info$pheno.correlation
  store12 <- population$info$bv.correlation

  if(bv.total>0 && (length(population$info$pheno.correlation)==0 || nrow(population$info$pheno.correlation)<bv.total)){
    population$info$pheno.correlation <- diag(1L, bv.total)
  }
  if(length(new.residual.correlation)>0){

    if(sum(new.residual.correlation==1) > ncol(new.residual.correlation)){

      if (requireNamespace("Matrix", quietly = TRUE)) {
        if(verbose) cat("Residual correlation matrix is only semi-definit. Modify slightly to ensure chol() working.\n")
        new.residual.correlation = as.matrix(Matrix::nearPD(new.residual.correlation)$mat)
      }

    }

    population$info$pheno.correlation <- t(chol(new.residual.correlation))
  }
  if(bv.total>0 && (length(population$info$bv.correlation)==0 || nrow(population$info$bv.correlation)<bv.total)){
    population$info$bv.correlation <- diag(1L, bv.total)
  }
  if(length(new.breeding.correlation)>0){
    population$info$bv.correlation <- new.breeding.correlation
  }


  if(replace.traits==FALSE){
    if(length(store1)>0){
      population$info$is.maternal[1:length(store1)] <- store1
    }
    if(length(store2)>0){
      population$info$is.paternal[1:length(store2)] <- store2
    }
    if(length(store3)>0){
      population$info$is.combi[1:length(store3)] <- store3
    }
    if(length(store4)>0){
      population$info$phenotypic.transform[1:length(store4)] <- store4
    }
    if(length(store5)>0){
      population$info$phenotypic.transform[1:length(store4)] <- store5
    }
    if(length(store6)>0){
      population$info$bv.random[1:length(store6)] <- store6
    }
    if(length(store7)>0){
      population$info$bv.random.variance[1:length(store7)] <- store7
    }
    if(length(store8)>0){
      population$info$bve.mult.factor[1:length(store8)] <- store8
    }
    if(length(store9)>0){
      population$info$bve.poly.factor[1:length(store9)] <- store9
    }
    if(length(store10)>0){
      population$info$base.bv[1:length(store10)] <- store10
    }
    if(length(store11)>0){
      population$info$pheno.correlation[1:nrow(store11), 1:nrow(store11)] <- store11
    }
    if(length(store12)>0){
      population$info$bv.correlation[1:nrow(store12), 1:nrow(store12)] <- store12
    }
  }


  for(generation in 1:nrow(population$info$size)){
    counter <- population$info$size[generation,] + 1

    if(replace.traits){

      to_add = population$info$bv.nr

      population$info$last.sigma.e.heritability <- NULL
      population$info$last.sigma.e.database <- NULL
      population$info$last.sigma.e.value <- NULL

      for(index in c(3,4,7:10, 19:22, 27:30)){

        sex_temp = 2-index%%2
        population$breeding[[generation]][[index]] <-  matrix(if(index==9 || index==10) {NA} else{0L}, nrow= to_add, ncol=counter[sex_temp]-1) # estimated breeding value
        # estimated breeding value 3,4
        # real genomic value 7,8
        # phenotype 9,10
        # Reliabilities 19,20
        # Last applied selection index 21,22
        # offspring phenotype 27,28
        # number of offspring used 29,30
      }
    } else{

      to_add = population$info$bv.nr - if(length(population$breeding[[generation]][[3]])>0){nrow(population$breeding[[generation]][[3]])} else {0}
      for(index in c(3,4,7:10, 19:22, 27:30)){
        sex_temp = 2-index%%2

        population$breeding[[generation]][[index]] <- rbind( population$breeding[[generation]][[index]],
                                                             matrix(if(index==9 || index==10) {NA} else{0L}, nrow= to_add, ncol=counter[sex_temp]-1)) # estimated breeding value
        # estimated breeding value 3,4
        # real genomic value 7,8
        # phenotype 9,10
        # Reliabilities 19,20
        # Last applied selection index 21,22
        # offspring phenotype 27,28
        # number of offspring used 29,30
      }
    }

  }

  if(bv.total>0){
    if(length(population$info$trait.name)>0 & replace.traits==FALSE){
      trait.name <- c(population$info$trait.name, trait.name)
    }
    population$info$trait.name <- trait.name
    if(length(trait.name)<bv.total){
      population$info$trait.name <- c(population$info$trait.name, paste0("Trait ", (length(trait.name)+1):bv.total))
    }
  }

  if(length(shuffle.traits)==0){
    if(length(shuffle.cor)>0){

      if(ncol(shuffle.cor)==population$info$bv.calc){
        shuffle.traits <- 1:population$info$bv.calc
      } else{
        shuffle.traits <- 1:ncol(shuffle.cor)
        warning(paste0("shuffle.traits not specified! use the first ", ncol(shuffle.cor), " traits"))
      }
    }
  }

  if(length(shuffle.traits)>0){
    if(length(shuffle.traits)==1){
      shuffle.traits <- which(population$info$bv.random==FALSE)
    }

    if(population$info$founder_multi){
      population$info$founder_multi_calc = FALSE
      for(index in 1:population$info$bv.nr){
        if(length(population$info$real.bv.add[[index]])>1 && sum(population$info$real.bv.add[[index]][,7]!=0)>0){
          population$info$founder_multi_calc = TRUE            }
      }
    }


    population$info$founder_multi_calc = population$info$founder_multi

    #stop()
    ## only needed for creating.trait()
    population$info$bv.calculated = FALSE
    population <- breeding.diploid(population, verbose = FALSE)


    bvs <- get.bv(population, gen=1)
    scalings <- sqrt(diag(stats::var(t(bvs))))
    for(bvnr in shuffle.traits){
      if(length(population$info$real.bv.add[[bvnr]])>0){
        population$info$real.bv.add[[bvnr]][,3:5] <- population$info$real.bv.add[[bvnr]][,3:5] / scalings[bvnr] * scalings[1]
      }


      if(length(population$info$real.bv.mult[[bvnr]])>0){
        population$info$real.bv.mult[[bvnr]][,5:13] <- population$info$real.bv.mult[[bvnr]][,5:13] / scalings[bvnr] * scalings[1]
      }

      if(length(population$info$real.bv.dice[[bvnr]])>0){
        population$info$real.bv.dice[[bvnr]][[2]] <- population$info$real.bv.dice[[bvnr]][[2]] / scalings[bvnr] * scalings[1]
      }


    }
    population$info$bv.calculated <- FALSE
    population$info$bv.calculated.partly <- NULL

    eigen_gen <- eigen(shuffle.cor)
    if(sum(eigen_gen$values<0)>0){
      if(verbose){
        warning("Genetic covariance matrix is not positive definit.")
        cat("Genetic covariance matrix is not positive definit.\n")
      }
      if(verbose) cat("Generate projection on the set of positive definit matrices:")

      test <- eigen_gen

      test$values[test$values<0] <- 0
      M <- diag(test$values)

      S <- test$vectors

      newA <- S %*% M %*% solve(S)

      diag(newA) <- diag(newA) + 0.001 # Avoid numerical issues with inversion
      newA <- newA * matrix(1/sqrt(diag(newA)), nrow=nrow(newA), ncol=nrow(newA), byrow=TRUE) * matrix(1/sqrt(diag(newA)), nrow=nrow(newA), ncol=nrow(newA), byrow=FALSE)
      if(verbose) cat("new suggested genetic correlation matrix:\n")
      shuffle.cor <- newA
      if(verbose) print(round(shuffle.cor, digits=3))
    }

    if(sum(shuffle.cor==1) > ncol(shuffle.cor)){

      if (requireNamespace("Matrix", quietly = TRUE)) {
        if(verbose) cat("Genetic correlation matrix is only semi-definit. Modify slightly to ensure chol() working.\n")
        shuffle.cor = as.matrix(Matrix::nearPD(shuffle.cor)$mat)
      }

    }

    LT <- chol(shuffle.cor)
    if(nrow(LT)!=length(shuffle.traits)){
      stop("Dimension of shuffle correlation matrix doesnt work with traits to shuffle")
    } else{

      population$info$bv.correlation[shuffle.traits,shuffle.traits] <- t(LT) %*% LT
      if(sum(abs(population$info$bv.correlation[shuffle.traits,shuffle.traits]- shuffle.cor))>0.0001){
        warning("No covariance matrix for genetic correlation given! Values above diagonal used.")
      }

      store.add <- population$info$real.bv.add
      store.mult <- population$info$real.bv.mult
      store.dice <- population$info$real.bv.dice


      col <- 1
      for(index in shuffle.traits){
        new.add <- new.mult <- new.dice1 <- new.dice2 <- NULL
        row <- 1
        add.list = mult.list = list()
        for(index2 in shuffle.traits){
          if(length(store.add[[index2]])>0){

            if(LT[row,col]!=0){
              temp1 <- store.add[[index2]] %*% diag(c(1,1,rep(LT[row,col],3),1,1,1))
            } else{
              temp1 <- NULL
            }
            if(length(temp1)>0){
              zeros <- rowSums(abs(temp1[,3:5, drop=FALSE]))
              temp1 <- temp1[zeros>0,,drop=FALSE]
            }
            add.list[[index2]] = temp1

          }
          if(length(store.mult[[index2]])>0){

            if(LT[row,col]!=0){
              temp1 <- store.mult[[index2]] %*% diag(c(1,1,1,1,rep(LT[row,col],9)))
            } else{
              temp1 <- NULL
            }
            if(length(temp1)>0){
              zeros <- rowSums(abs(temp1[,5:13, drop=FALSE]))
              temp1 <- temp1[zeros>0,,drop=FALSE]
            }
            mult.list[[index2]] = temp1
          }
          if(length(store.dice[[index2]])>0){
            before <- length(new.dice2)
            new.dice1 <- c(new.dice1,store.dice[[index2]][[1]])
            new.dice2 <- c(new.dice2,store.dice[[index2]][[2]])
            for(index3 in (before+1):length(new.dice2)){
              new.dice2[[index3]] <- new.dice2[[index3]] * LT[row,col]
            }
          }
          row <- row +1
        }

        new.add = do.call(rbind, add.list)
        new.mult = do.call(rbind, mult.list)

        # DONT REMOVE NULL - MORE WORK NEEDED HERE!
        if(length(new.add)==0){

        } else{
          population$info$real.bv.add[[index]] <- new.add
        }
        if(length(new.mult)==0){

        } else{
          population$info$real.bv.mult[[index]] <- new.mult
        }
        if(length(new.add)==0){

        } else{
          population$info$real.bv.dice[[index]] <- list(new.dice1,new.dice2)
        }
        col <- col +1
      }

    }

    for(index in shuffle.traits){
      population$info$real.bv.length[1] <- max(population$info$real.bv.length[1], if(length(population$info$real.bv.add[[index]])>0){index} else{0})
      population$info$real.bv.length[2] <- max(population$info$real.bv.length[2], if(length(population$info$real.bv.mult[[index]])>0){index} else{0})
      population$info$real.bv.length[3] <- max(population$info$real.bv.length[3], if(length(population$info$real.bv.dice[[index]][[1]])>0){index} else{0})
    }


  }

  if(population$info$bv.nr>0){
    for(index in 1:population$info$bv.nr){
      if(length(population$info$real.bv.add)>= index){
        if(length(population$info$real.bv.add[[index]])>1){
          t <- population$info$real.bv.add[[index]]
          take <- sort(t[,1]+ cumsum(c(0,population$info$snp))[t[,2]], index.return=TRUE)
          t <- t[take$ix,,drop=FALSE]
          take <- sort(t[,1]+ t[,2] * 10^10 + t[,7]*10^8+ t[,8]*10^7)
          keep <- c(0,which(diff(take)!=0), length(take))
          if(length(keep) < (nrow(t)+1)){

            for(index2 in (2:(length(keep)))[diff(keep)>1]){
              t[keep[index2],3:5] <- colSums(t[(keep[index2-1]+1):keep[index2],3:5, drop=FALSE])
            }
          }
          population$info$real.bv.add[[index]] <- t[keep,,drop=FALSE]
        }
      }

    }
  }



  if(length(fixed.effects)==0){
    fixed.effects <- matrix(0, nrow= population$info$bv.nr, ncol=0)
    temp1 <- TRUE
  } else{
    temp1 <- FALSE
  }

  if(length(population$info$fixed.effects)>0 & temp1){
    fixed.effects[1:nrow(population$info$fixed.effects), 1:ncol(population$info$fixed.effects)] <- population$info$fixed.effects
  }
  population$info$fixed.effects <- fixed.effects


  # Add traits with no generated phenotypes / litter or pen effects.
  temp1 <- rep(0, population$info$bv.nr)
  for(gen in 1:length(population$breeding)){
    for(sex in 1:2){
      if(length(population$breeding[[gen]][[sex]])>0){
        for(index in 1:length(population$breeding[[gen]][[sex]])){
          population$breeding[[gen]][[sex]][[index]][[15]] <- temp1
          population$breeding[[gen]][[sex]][[index]][[28]] <- c(population$breeding[[gen]][[sex]][[index]][[28]], rep(0, ncol(population$info$fixed.effects) - length(population$breeding[[gen]][[sex]][[index]][[28]])))
          population$breeding[[gen]][[sex]][[index]][[29]] <- temp1
          population$breeding[[gen]][[sex]][[index]][[30]] <- temp1
        }
      }
    }
  }


  if(remove.invalid.qtl && length(population$info$real.bv.add)>1){
    for(index in 1:(length(population$info$real.bv.add)-1)){
      if(length(population$info$real.bv.add[[index]])>1){
        removes <- which(population$info$real.bv.add[[index]][,1] > population$info$snp[population$info$real.bv.add[[index]][,2]] |  population$info$real.bv.add[[index]][,1] < 1 | round(population$info$real.bv.add[[index]][,1])!= population$info$real.bv.add[[index]][,1]|  population$info$real.bv.add[[index]][,2] < 1 | round(population$info$real.bv.add[[index]][,2])!= population$info$real.bv.add[[index]][,2])
        if(length(removes)>0){
          population$info$real.bv.add[[index]] <- population$info$real.bv.add[[index]][-removes,,drop=FALSE]
          if(verbose) cat(paste0(length(removes), " QTL-effects entered on markers that do not exist for ", population$info$trait.name[index], ".\n"))
          if(verbose) cat(paste0(nrow(population$info$real.bv.add[[index]]), " QTL-effects remain.\n"))
        }
      }

    }
    for(index in 1:(length(population$info$real.bv.mult)-1)){
      if(length(population$info$real.bv.mult[[index]])>1){
        removes <- which(population$info$real.bv.mult[[index]][,1] > population$info$snp[population$info$real.bv.mult[[index]][,2]] |  population$info$real.bv.mult[[index]][,1] < 1 | round(population$info$real.bv.mult[[index]][,1])!= population$info$real.bv.mult[[index]][,1] |  population$info$real.bv.mult[[index]][,2] < 1 | round(population$info$real.bv.mult[[index]][,2])!= population$info$real.bv.mult[[index]][,2])
        if(length(removes)>0){
          population$info$real.bv.mult[[index]] <- population$info$real.bv.mult[[index]][-removes,,drop=FALSE]
          if(verbose) cat(paste0(length(removes), " QTL-effects entered on markers that do not exist for ", population$info$trait.name[index], ".\n"))
          if(verbose) cat(paste0(nrow(population$info$real.bv.mult[[index]]), " QTL-effects remain.\n"))
        }
      }

    }
    for(index in 1:(length(population$info$real.bv.mult)-1)){
      if(length(population$info$real.bv.mult[[index]])>1){
        removes <- which(population$info$real.bv.mult[[index]][,3] > population$info$snp[population$info$real.bv.mult[[index]][,4]] |  population$info$real.bv.mult[[index]][,3] < 1 | round(population$info$real.bv.mult[[index]][,3])!= population$info$real.bv.mult[[index]][,3]|  population$info$real.bv.mult[[index]][,4] < 1 | round(population$info$real.bv.mult[[index]][,4])!= population$info$real.bv.mult[[index]][,4])
        if(length(removes)>0){
          population$info$real.bv.mult[[index]] <- population$info$real.bv.mult[[index]][-removes,,drop=FALSE]
          if(verbose) cat(paste0(length(removes), " QTL-effects entered on markers that do not exist for ", population$info$trait.name[index], ".\n"))
          if(verbose) cat(paste0(nrow(population$info$real.bv.mult[[index]]), " QTL-effects remain.\n"))
        }
      }
    }

  }

  population$info$neff <- list()
  if(length(population$info$real.bv.add)>1){
    for(index in 1:(length(population$info$real.bv.add)-1)){
      if(length(population$info$real.bv.add[[index]])>0){
        population$info$neff[[index]] <- 1:nrow(population$info$real.bv.add[[index]])
      }
    }
  }

  # recalculate - manual check:
  {
    if(length(unlist(population$info$real.bv.mult)) > 1){
      population$info$recalculate.possible = FALSE
    }

    if(length(unlist(population$info$real.bv.dice)) > 1){
      population$info$recalculate.possible = FALSE
    }

    if(length(population$info$real.bv.add) > 1){
      for(index in 1:(length(population$info$real.bv.add)-1)){
        if(is.matrix(population$info$real.bv.add[[index]]) && sum(population$info$real.bv.add[[index]][,7:8] != 0)>0){
          population$info$recalculate.possible = FALSE
        }
      }
    }

    if(max(population$info$real.bv.length) < population$info$bv.nr){
      population$info$recalculate.possible = FALSE
    }

    population$info$e0_activ = NULL
    population$info$e1_activ = NULL
    population$info$e2_activ = NULL
    population$info$e0_mat = NULL
    population$info$e1_mat = NULL
    population$info$e2_mat = NULL
  }
  #stop()


  population$info$bv.calculated = FALSE
  population = breeding.diploid(population, verbose=FALSE, use.recalculate.manual = use.recalculate.manual)


  population$info$pool_effects_calc = FALSE
  population$info$pool_list = NULL
  population$info$bypool_list = NULL

  if(bv.standard){
    population <- bv.standardization(population, mean.target = mean.target, var.target = var.target, set.zero = set.zero)

      population$info$bv.calculated = FALSE
      population = breeding.diploid(population, verbose=FALSE, use.recalculate.manual = use.recalculate.manual)

  }

  if(length(trait_location)>0){
    population$info$trait.location = trait_location
    population$info$trait.nr = trait_nr
  }

  if(gxe.combine & length(trait_nr)>0){

    traits = unique(trait_nr)
    for(index in traits){
      population <- combine.traits(population, combine.traits = which(trait_nr == index))
    }
  }

  return(population)
}