R/breeding.diploid.R

Defines functions breeding.diploid

Documented in breeding.diploid

'#
  Authors
Torsten Pook, torsten.pook@uni-goettingen.de

Copyright (C) 2017 -- 2020  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.
'#


#' Breeding function
#'
#' Function to simulate a step in a breeding scheme
#' @param population Population list
#' @param mutation.rate Mutation rate in each marker (default: 10^-8)
#' @param remutation.rate Remutation rate in each marker (default: 10^-8)
#' @param recombination.rate Average number of recombination per 1 length unit (default: 1M)
#' @param selection.m Selection criteria for male individuals (Set to "random" to randomly select individuals - this happens automatically when no the input in selection.criteria has no input ((usually breeding values)))
#' @param selection.f Selection criteria for female individuals (default: selection.m , alt: "random", function")
#' @param selection.function.matrix Manuel generation of a temporary selection function (Use BVs instead!)
#' @param selection.size Number of selected individuals for breeding (default: c(0,0) - alt: positive numbers)
#' @param breeding.size Number of individuals to generate
#' @param breeding.sex Share of female animals (if single value is used for breeding size; default: 0.5)
#' @param breeding.sex.random If TRUE randomly chose sex of new individuals (default: FALSE - use expected values)
#' @param relative.selection Use best.selection.ratio instead!
#' @param recom.f.indicator Use step function for recombination map (transform snp.positions if possible instead)
#' @param add.gen Generation you want to add the new individuals to (default: New generation)
#' @param duplication.rate Share of recombination points with a duplication (default: 0 - DEACTIVATED)
#' @param duplication.length Average length of a duplication (Exponentially distributed)
#' @param duplication.recombination Average number of recombinations per 1 length uit of duplication (default: 1)
#' @param new.class Migration level of newly generated individuals (default: 0)
#' @param bve If TRUE perform a breeding value estimation (default: FALSE)
#' @param bve.gen Generations of individuals to consider in breeding value estimation (default: NULL)
#' @param bve.cohorts Cohorts of individuals to consider in breeding value estimation (default: NULL)
#' @param bve.database Groups of individuals to consider in breeding value estimation (default: NULL)
#' @param bve.insert.gen Generations of individuals to compute breeding values for (default: all groups in bve.database)
#' @param bve.insert.cohorts Cohorts of individuals to compute breeding values for (default: all groups in bve.database)
#' @param bve.insert.database Groups of individuals to compute breeding values for (default: all groups in bve.database)
#' @param bve.pseudo If set to TRUE the breeding value estimation will be simulated with resulting accuracy bve.pseudo.accuracy (default: 1)
#' @param bve.pseudo.accuracy The accuracy to be obtained in the "pseudo" - breeding value estimation
#' @param sigma.e Enviromental variance (default: 100)
#' @param sigma.e.gen Generations to consider when estimating sigma.e when using hertability
#' @param sigma.e.cohorts Cohorts to consider when estimating sigma.e when using hertability
#' @param sigma.e.database Groups to consider when estimating sigma.e when using hertability
#' @param forecast.sigma.g Set FALSE to not estimate sigma.g (Default: TRUE)
#' @param heritability Use sigma.e to obtain a certain heritability (default: NULL)
#' @param repeatability Set this to control the share of the residual variance (sigma.e) that is permanent (there for each observation)
#' @param sigma.g Genetic variance (default: 100 - only used if not computed via estimate.sigma.g^2 in der Zuchtwertschaetzung (Default: 100)
#' @param sigma.g.gen Generations to consider when estimating sigma.g
#' @param sigma.g.cohorts Cohorts to consider when estimating sigma.g
#' @param sigma.g.database Groups to consider when estimating sigma.g
#' @param phenotyping Quick acces to phenotyping for (all: "all", non-phenotyped: "non_obs", non-phenotyped male: "non_obs_m", non-phenotyped female: "non_obs_f")
#' @param phenotyping.child Starting phenotypes of newly generated individuals (default: "mean" of both parents, "obs" - regular observation, "zero" - 0)
#' @param relationship.matrix Method to calculate relationship matrix for the breeding value estimation (Default: "vanRaden", alt: "kinship", "CE", "non_stand", "CE2", "CM")
#' @param relationship.matrix.ogc Method to calculate relationship matrix for OGC (Default: "kinship", alt: "vanRaden", "CE", "non_stand", "CE2", "CM")
#' @param delete.haplotypes Generations for with haplotypes of founders can be deleted (only use if storage problem!)
#' @param delete.individuals Generations for with individuals are completley deleted (only use if storage problem!)
#' @param praeimplantation Only use matings the lead to a specific genotype in a specific marker
#' @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 fixed.breeding Set of targeted matings to perform
#' @param fixed.breeding.best Perform targeted matings in the group of selected individuals
#' @param max.offspring Maximum number of offspring per individual (default: c(Inf,Inf) - (m,w))
#' @param max.litter Maximum number of offspring per individual (default: c(Inf,Inf) - (m,w))
#' @param max.mating.pair Set to the maximum number of matings between two individuals (default: Inf)
#' @param store.breeding.totals If TRUE store information on selected animals in $info$breeding.totals
#' @param multiple.bve Way to handle multiple traits in bv/selection (default: "add", alt: "ranking")
#' @param multiple.bve.weights.m Weighting between traits when using "add" (default: 1)
#' @param multiple.bve.weights.f Weighting between traits when using "add" (default: same as multiple.bve.weights.m)
#' @param store.bve.data If TRUE store information of bve in $info$bve.data
#' @param fixed.assignment Set TRUE for targeted mating of best-best individual till worst-worst (of selected). set to "bestworst" for best-worst mating
#' @param selection.highest If 0 individuals with lowest bve are selected as best individuals (default c(1,1) - (m,w))
#' @param same.sex.activ If TRUE allow matings of individuals of same sex
#' @param same.sex.sex Probability to use female individuals as parents (default: 0.5)
#' @param same.sex.selfing Set to TRUE to allow for selfing when using same.sex matings
#' @param selfing.mating If TRUE generate new individuals via selfing
#' @param selfing.sex Share of female individuals used for selfing (default: 0.5)
#' @param multiple.bve.scale.m Default: "bv_sd"; Set to "pheno_sd" when using gains per phenotypic SD, "unit" when using gains per unit, "bve" when using estimated breeding values
#' @param multiple.bve.scale.f Default: "bv_sd"; Set to "pheno_sd" when using gains per phenotypic SD, "unit" when using gains per unit, "bve" when using estimated breeding values
#' @param class.m Migrationlevels of male individuals to consider for mating process (default: 0)
#' @param class.f Migrationlevels of female individuals to consider for mating process (default: 0)
#' @param save.recombination.history If TRUE store the time point of each recombination event
#' @param martini.selection If TRUE use the group of non-selected individuals as second parent
#' @param BGLR.bve If TRUE use BGLR to perform breeding value estimation
#' @param BGLR.model Select which BGLR model to use (default: "RKHS", alt: "BRR", "BL", "BayesA", "BayesB", "BayesC")
#' @param BGLR.burnin Number of burn-in steps in BGLR (default: 1000)
#' @param BGLR.iteration Number of iterations in BGLR (default: 5000)
#' @param BGLR.print If TRUE set verbose to TRUE in BGLR
#' @param BGLR.save Method to use in BGLR (default: "RKHS" - alt: NON currently)
#' @param BGLR.save.random Add random number to store location of internal BGLR computations (only needed when simulating a lot in parallel!)
#' @param copy.individual If TRUE copy the selected father for a mating
#' @param copy.individual.m If TRUE generate exactly one copy of all selected male in a new cohort (or more by setting breeding.size)
#' @param copy.individual.f If TRUE generate exactly one copy of all selected female in a new cohort (or more by setting breeding.size)
#' @param copy.individual.keep.pheno Set to FALSE to not keep estimated breeding values in case of use of copy.individuals
#' @param copy.individual.keep.bve Set to FALSE to not keep estimated breeding value in case of use of copy.individuals
#' @param dh.mating If TRUE generate a DH-line in mating process
#' @param dh.sex Share of DH-lines generated from selected female individuals
#' @param n.observation Number of phenotypes generated per individuals (influences enviromental variance)
#' @param bve.0isNA Individuals with phenotype 0 are used as NA in breeding value estimation
#' @param phenotype.bv If TRUE use phenotype as estimated breeding value
#' @param delete.same.origin If TRUE delete recombination points when genetic origin of adjacent segments is the same
#' @param estimate.u If TRUE estimate u in breeding value estimation (Y = Xb + Zu + e)
#' @param gwas.u If TRUE estimate u via GWAS (relevant for gene editing)
#' @param approx.residuals If FALSE calculate the variance for each marker separatly instead of using a set variance (doesnt change order - only p-values)
#' @param gwas.gen Generations to consider in GWAS analysis
#' @param gwas.cohorts Cohorts to consider in GWAS analysis
#' @param gwas.database Groups to consider in GWAS analysis
#' @param gwas.group.standard If TRUE standardize phenotypes by group mean
#' @param y.gwas.used What y value to use in GWAS study (Default: "pheno", alt: "bv", "bve")
#' @param remove.effect.position If TRUE remove real QTLs in breeding value estimation
#' @param estimate.add.gen.var If TRUE estimate additive genetic variance and heritability based on parent model
#' @param estimate.pheno.var If TRUE estimate total variance in breeding value estimation
#' @param store.comp.times If TRUE store computation times in $info$comp.times (default: TRUE)
#' @param store.comp.times.bve If TRUE store computation times of breeding value estimation in $info$comp.times.bve (default: TRUE)
#' @param store.comp.times.generation If TRUE store computation times of mating simulations in $info$comp.times.generation (default: TRUE)
#' @param ogc If TRUE use optimal genetic contribution theory to perform selection ( This requires the use of the R-package optiSel)
#' @param ogc.target Target of OGC (default: "min.sKin" - minimize inbreeding; alt: "max.BV" / "min.BV" - maximize genetic gain; both under constrains selected below)
#' @param ogc.uniform This corresponds to the uniform constrain in optiSel
#' @param ogc.lb This corresponds to the lb constrain in optiSel
#' @param ogc.ub This corresponds to the ub constrain in optiSel
#' @param ogc.ub.sKin This corresponds to the ub.sKin constrain in optiSel
#' @param ogc.lb.BV This corresponds to the lb.BV constrain in optiSel
#' @param ogc.ub.BV This corresponds to the ub.BV constrain in optiSel
#' @param ogc.eq.BV This corresponds to the eq.BV constrain in optiSel
#' @param ogc.ub.sKin.increase This corresponds to the upper bound (current sKin + ogc.ub.sKin.increase) as ub.sKin in optiSel
#' @param ogc.lb.BV.increase This corresponds to the lower bound (current BV + ogc.lb.BV.increase) as lb.BV in optiSel
#' @param gene.editing.offspring If TRUE perform gene editing on newly generated individuals
#' @param gene.editing.best If TRUE perform gene editing on selected individuals
#' @param gene.editing.offspring.sex Which sex to perform editing on (Default c(TRUE,TRUE), mw)
#' @param gene.editing.best.sex Which sex to perform editing on (Default c(TRUE,TRUE), mw)
#' @param nr.edits Number of edits to perform per individual
#' @param import.position.calculation Function to calculate recombination point into adjacent/following SNP
#' @param emmreml.bve If TRUE use REML estimator from R-package EMMREML in breeding value estimation
#' @param rrblup.bve If TRUE use REML estimator from R-package rrBLUP in breeding value estimation
#' @param sommer.bve If TRUE use REML estimator from R-package sommer in breeding value estimation
#' @param bve.direct.est If TRUE predict BVEs in direct estimation according to vanRaden 2008 method 2 (default: TRUE)
#' @param sequenceZ Split genomic matric into parts (relevent if high memory usage)
#' @param maxZ Number of SNPs to consider in each part of sequenceZ
#' @param maxZtotal Number of matrix entries to consider jointly (maxZ = maxZtotal/number of animals)
#' @param delete.sex Remove all individuals from these sex from generation delete.individuals (default: 1:2 ; note:delete individuals=NULL)
#' @param gen.architecture.m Genetic architecture for male animal (default: 0 - no transformation)
#' @param gen.architecture.f Genetic architecture for female animal (default: gen.architecture.m - no transformation)
#' @param add.architecture List with two vectors containing (A: length of chromosomes, B: position in cM of SNPs)
#' @param ncore Cores used for parallel computing in compute.snps
#' @param Z.integer If TRUE save Z as a integer in parallel computing
#' @param store.effect.freq If TRUE store the allele frequency of effect markers per generation
#' @param backend Chose the used backend (default: "doParallel", alt: "doMPI")
#' @param Rprof Store computation times of each function
#' @param miraculix If TRUE use miraculix to perform computations (ideally already generate population in creating.diploid with this; default: automatic detection from population list)
#' @param miraculix.mult If TRUE use miraculix for matrix multiplications even if miraculix is not used for storage
#' @param miraculix.chol Set to FALSE to deactive miraculix based Cholesky-decomposition (default: TRUE)
#' @param miraculix.cores Number of cores used in miraculix applications (default: 1)
#' @param miraculix.destroyA If FALSE A will not be destroyed in the process of inversion (less computing / more memory)
#' @param best.selection.ratio.m Ratio of the frequency of the selection of the best best animal and the worst best animal (default=1)
#' @param best.selection.ratio.f Ratio of the frequency of the selection of the best best animal and the worst best animal (default=1)
#' @param best.selection.criteria.m Criteria to calculate this ratio (default: "bv", alt: "bve", "pheno")
#' @param best.selection.criteria.f Criteria to calculate this ratio (default: "bv", alt: "bve", "pheno")
#' @param best.selection.manual.ratio.m vector containing probability to draw from for every individual (e.g. c(0.1,0.2,0.7))
#' @param best.selection.manual.ratio.f vector containing probability to draw from for every individual (e.g. c(0.1,0.2,0.7))
#' @param best.selection.manual.reorder Set to FALSE to not use the order from best to worst selected individual but plain order based on database-order
#' @param selection.criteria What to use in the selection proces (default: "bve", alt: "bv", "pheno")
#' @param bve.class Consider only animals of those class classes in breeding value estimation (default: NULL - use all)
#' @param parallel.generation Set TRUE to active parallel computing in animal generation
#' @param ncore.generation Number of cores to use in parallel generation
#' @param randomSeed Set random seed of the process
#' @param randomSeed.generation Set random seed for parallel generation process
#' @param new.selection.calculation If TRUE recalculate breeding values obtained by selection.function.matrix
#' @param name.cohort Name of the newly added cohort
#' @param add.class.cohorts Migration levels of all cohorts selected for reproduction are automatically added to class.m/class.f (default: TRUE)
#' @param display.progress Set FALSE to not display progress bars. Setting verbose to FALSE will automatically deactive progress bars
#' @param ignore.best Not consider the top individuals of the selected individuals (e.g. to use 2-10 best individuals)
#' @param combine Copy existing individuals (e.g. to merge individuals from different groups in a joined cohort). Individuals to use are used as the first parent
#' @param repeat.mating Generate multiple mating from the same dam/sire combination (first column: number of offspring; second column: probability)
#' @param repeat.mating.copy Generate multiple copies from a copy action (combine / copy.individuals.m/f) (first column: number of offspring; second column: probability)
#' @param repeat.mating.fixed Vector containing number of times each mating is repeated. This will overwrite sampling from repeat.mating / repeat.mating.copy (default: NULL)
#' @param repeat.mating.overwrite Set to FALSE to not use the current repeat.mating / repeat.mating.copy input as the new standard values (default: TRUE)
#' @param time.point Time point at which the new individuals are generated
#' @param creating.type Technique to generate new individuals (usage in web-based application)
#' @param multiple.observation Set TRUE to allow for more than one phenotype observation per individual (this will decrease enviromental variance!)
#' @param phenotyping.gen Vector of generation from which to generate additional phenotypes
#' @param phenotyping.cohorts Vector of cohorts from which to generate additional phenotype
#' @param phenotyping.database Matrix of groups from which to generate additional phenotypes
#' @param share.phenotyped Share of the individuals to phenotype
#' @param reduced.selection.panel.m Use only a subset of individuals of the potential selected ones ("Split in user-interface")
#' @param reduced.selection.panel.f Use only a subset of individuals of the potential selected ones ("Split in user-interface")
#' @param breeding.all.combination Set to TRUE to automatically perform each mating combination possible exactly ones.
#' @param depth.pedigree Depth of the pedigree in generations (default: 7)
#' @param depth.pedigree.ogc Depth of the pedigree in generations (default: 7)
#' @param bve.avoid.duplicates If set to FALSE multiple generatations of the same individual can be used in the bve (only possible by using copy.individual to generate individuals)
#' @param report.accuracy Report the accuracy of the breeding value estimation
#' @param share.genotyped Share of individuals newly generated individuals that are genotyped
#' @param added.genotyped Share of individuals that is additionally genotyped (only for copy.individuals)
#' @param singlestep.active Set TRUE to use single step in breeding value estimation (only implemented for vanRaden- G matrix and without use sequenceZ) (Legarra 2014)
#' @param remove.non.genotyped Set to FALSE to manually include non-genotyped individuals in genetic BVE, single-step will deactive this as well
#' @param fast.uhat Set to FALSE to  derive inverse of A in rrBLUP
#' @param offspring.bve.parents.gen Generations to consider to derive phenotype from offspring phenotypes
#' @param offspring.bve.parents.database Groups to consider to derive phenotype from offspring phenotypes
#' @param offspring.bve.parents.cohorts Cohorts to consider to derive phenotype from offspring phenotypes
#' @param offspring.bve.offspring.gen Active generations for import of offspring phenotypes
#' @param offspring.bve.offspring.database Active groups for import of offspring phenotypes
#' @param offspring.bve.offspring.cohorts Active cohorts for import of offspring phenotypes
#' @param sommer.multi.bve Set TRUE to use a mulit-trait model in the R-package sommer for BVE
#' @param calculate.reliability Set TRUE to calculate a reliability when performing Direct-Mixed-Model BVE
#' @param selection.m.gen Generations available for selection of paternal parent
#' @param selection.f.gen Generations available for selection of maternal parent
#' @param selection.m.database Groups available for selection of paternal parent
#' @param selection.f.database Groups available for selection of maternal parent
#' @param selection.m.cohorts Cohorts available for selection of paternal parent
#' @param selection.f.cohorts Cohorts available for selection of maternal parent
#' @param selection.m.miesenberger Use Weighted selection index according to Miesenberger 1997 for paternal selection
#' @param selection.f.miesenberger Use Weighted selection index according to Miesenberger 1997 for maternal selection
#' @param miesenberger.trafo Ignore all eigenvalues below this threshold and apply dimension reduction (default: 0 - use all)
#' @param selection.miesenberger.reliability.est If available reliability estimated are used. If not use default:"estimated" (SD BVE / SD Pheno), alt: "heritability", "derived" (cor(BVE,BV)^2) as replacement
#' @param culling.gen Generations to consider to culling
#' @param culling.database Groups to consider to culling
#' @param culling.cohort Cohort to consider to culling
#' @param culling.time Age of the individuals at culling
#' @param culling.name Name of the culling action (user-interface stuff)
#' @param culling.bv1 Reference Breeding value
#' @param culling.share1 Probability of death for individuals with bv1
#' @param culling.bv2 Alternative breeding value (linear extended for other bvs)
#' @param culling.share2 Probability of death for individuals with bv2
#' @param culling.index Genomic index (default:0 - no genomic impact, use: "lastindex" to use the last selection index applied in selection)
#' @param culling.single Set to FALSE to not apply the culling module on all individuals of the cohort
#' @param culling.all.copy Set to FALSE to not kill copies of the same individual in the culling module
#' @param verbose Set to FALSE to not display any prints
#' @param bve.parent.mean Set to TRUE to use the average parental performance as the breeding value estimate
#' @param bve.grandparent.mean Set to TRUE to use the average grandparental performance as the breeding value estimate
#' @param bve.mean.between Select if you want to use the "bve", "bv", "pheno" or "bvepheno" to form the mean (default: "bvepheno" - if available bve, else pheno)
#' @param mas.bve If TRUE use marker assisted selection in the breeding value estimation
#' @param mas.markers Vector containing markers to be used in marker assisted selection
#' @param mas.number If no markers are provided this nr of markers is selected (if single marker QTL are present highest effect markers are prioritized)
#' @param mas.effects Effects assigned to the MAS markers (Default: estimated via lm())
#' @param threshold.selection Minimum value in the selection index selected individuals have to have
#' @param threshold.sign Pick all individuals above (">") the threshold. Alt: ("<", "=", "<=", ">=")
#' @param input.phenotype Select what to use in BVE (default: own phenotype ("own"), offspring phenotype ("off"), their average ("mean") or a weighted average ("weighted"))
#' @param bve.ignore.traits Vector of traits to ignore in the breeding value estimation (default: NULL, use: "zero" to not consider traits with 0 index weight in multiple.bve.weights.m/.w)
#' @param bv.ignore.traits Vector of traits to ignore in the calculation of the genomic value (default: NULL; Only recommended for high number of traits and experienced users!)
#' @param genotyped.gen Generations to generate genotype data (that can be used in a BVE)
#' @param genotyped.database Groups to generate genotype data (that can be used in a BVE)
#' @param genotyped.cohorts Cohorts to generate genotype data (that can be used in a BVE)
#' @param genotyped.share Share of individuals in genotyped.gen/database/cohort to generate genotype data from (default: 1)
#' @param genotyped.array Genotyping array used
#' @param bve.imputation Set to FALSE to not perform imputation up to the highest marker density of genotyping data that is available
#' @param bve.imputation.errorrate Share of errors in the imputation procedure (default: 0)
#' @param sex.s Specify which newly added individuals are male (1) or female (2)
#' @param new.bv.observation.gen (OLD! use phenotyping.gen) Vector of generation from which to generate additional phenotypes
#' @param new.bv.observation.cohorts (OLD! use phenotyping.cohorts)Vector of cohorts from which to generate additional phenotype
#' @param new.bv.observation.database  (OLD! use phenotyping.database) Matrix of groups from which to generate additional phenotypes
#' @param best1.from.group (OLD!- use selection.m.database) Groups of individuals to consider as First Parent / Father (also female individuals are possible)
#' @param best2.from.group (OLD!- use selection.f.database) Groups of individuals to consider as Second Parent / Mother (also male individuals are possible)
#' @param best1.from.cohort (OLD!- use selection.m.cohorts) Groups of individuals to consider as First Parent / Father (also female individuals are possible)
#' @param best2.from.cohort (OLD! - use selection.f.cohorts) Groups of individuals to consider as Second Parent / Mother (also male individuals are possible)
#' @param new.bv.observation (OLD! - use phenotyping) Quick acces to phenotyping for (all: "all", non-phenotyped: "non_obs", non-phenotyped male: "non_obs_m", non-phenotyped female: "non_obs_f")
#' @param reduce.group (OLD! - use culling modules) Groups of animals for reduce to a new size (by changing class to -1)
#' @param reduce.group.selection (OLD! - use culling modules) Selection criteria for reduction of groups (cf. selection.m / selection.f - default: "random")
#' @param new.bv.child (OLD! - use phenotyping.child) Starting phenotypes of newly generated individuals (default: "mean" of both parents, "obs" - regular observation, "zero" - 0)
#' @param computation.A (OLD! - use relationship.matrix) Method to calculate relationship matrix for the breeding value estimation (Default: "vanRaden", alt: "kinship", "CE", "non_stand", "CE2", "CM")
#' @param computation.A.ogc (OLD! use relationship.matrix.ogc) Method to calculate pedigree matrix in OGC (Default: "kinship", alt: "vanRaden", "CE", "non_stand", "CE2", "CM")
#' @param new.phenotype.correlation (OLD! - use new.residual.correlation!) Correlation of the simulated enviromental variance
#' @param avoid.mating.fullsib Set to TRUE to not generate offspring of full siblings
#' @param avoid.mating.halfsib Set to TRUE to not generate offspring from half or full siblings
#' @param bve.per.sample.sigma.e Set to FALSE to deactivate the use of a heritablity based on the number of observations generated per sample
#' @param bve.solve Provide solver to be used in BVE (default: "exact" solution via inversion, alt: "pcg", function with inputs A, b and output y_hat)
#' @examples
#' population <- creating.diploid(nsnp=1000, nindi=100)
#' population <- breeding.diploid(population, breeding.size=100, selection.size=c(25,25))
#' @return Population-list
#' @export


breeding.diploid <- function(population,
            mutation.rate = 10^-8,
            remutation.rate = 10^-8,
            recombination.rate = 1,
            selection.m = NULL,
            selection.f = NULL,
            new.selection.calculation = TRUE,
            selection.function.matrix = NULL,
            selection.size = 0,
            ignore.best = 0,
            breeding.size = 0,
            breeding.sex = NULL,
            breeding.sex.random = FALSE,
            relative.selection = FALSE,
            class.m = 0,
            class.f = 0,
            add.gen = 0,
            recom.f.indicator = NULL,
            duplication.rate = 0,
            duplication.length = 0.01,
            duplication.recombination = 1,
            new.class = 0L,
            bve = FALSE,
            sigma.e = NULL,
            sigma.g = 100,
            new.bv.child = NULL,
            phenotyping.child = NULL,
            relationship.matrix = "vanRaden",
            relationship.matrix.ogc = "kinship",
            computation.A = NULL,
            computation.A.ogc = NULL,
            delete.haplotypes = NULL,
            delete.individuals = NULL,
            fixed.breeding = NULL,
            fixed.breeding.best = NULL,
            max.offspring = Inf,
            max.litter = Inf,
            store.breeding.totals = FALSE,
            forecast.sigma.g = TRUE,
            multiple.bve = "add",
            store.bve.data = FALSE,
            fixed.assignment = FALSE,
            reduce.group = NULL,
            reduce.group.selection = "random",
            selection.highest = c(TRUE,TRUE),
            selection.criteria = NULL,
            same.sex.activ = FALSE,
            same.sex.sex = 0.5,
            same.sex.selfing = FALSE,
            selfing.mating = FALSE,
            selfing.sex = 0.5,
            praeimplantation = NULL,
            heritability = NULL,
            repeatability = NULL,
            save.recombination.history = FALSE,
            martini.selection = FALSE,
            BGLR.bve = FALSE,
            BGLR.model = "RKHS",
            BGLR.burnin = 500,
            BGLR.iteration = 5000,
            BGLR.print = FALSE,
            copy.individual = FALSE,
            copy.individual.m = FALSE,
            copy.individual.f = FALSE,
            dh.mating = FALSE,
            dh.sex = 0.5,
            n.observation = NULL,
            bve.0isNA = FALSE,
            phenotype.bv = FALSE,
            delete.same.origin = FALSE,
            remove.effect.position = FALSE,
            estimate.u = FALSE,
            new.phenotype.correlation = NULL,
            new.residual.correlation = NULL,
            new.breeding.correlation = NULL,
            estimate.add.gen.var = FALSE,
            estimate.pheno.var = FALSE,
            best1.from.group = NULL,
            best2.from.group = NULL,
            best1.from.cohort = NULL,
            best2.from.cohort = NULL,
            add.class.cohorts = TRUE,
            store.comp.times = TRUE,
            store.comp.times.bve = TRUE,
            store.comp.times.generation = TRUE,
            import.position.calculation = NULL,
            BGLR.save = "RKHS",
            BGLR.save.random = FALSE,
            ogc = FALSE,
            ogc.target = "min.sKin",
            ogc.uniform=NULL,
            ogc.ub = NULL,
            ogc.lb = NULL,
            ogc.ub.sKin = NULL,
            ogc.lb.BV = NULL,
            ogc.ub.BV = NULL,
            ogc.eq.BV = NULL,
            ogc.ub.sKin.increase = NULL,
            ogc.lb.BV.increase = NULL,
            emmreml.bve = FALSE,
            rrblup.bve = FALSE,
            sommer.bve = FALSE,
            sommer.multi.bve=FALSE,
            nr.edits = 0,
            gene.editing.offspring = FALSE,
            gene.editing.best = FALSE,
            gene.editing.offspring.sex = c(TRUE,TRUE),
            gene.editing.best.sex = c(TRUE,TRUE),
            gwas.u = FALSE,
            approx.residuals = TRUE,
            sequenceZ = FALSE,
            maxZ = 5000,
            maxZtotal = 0,
            delete.sex = 1:2,
            gwas.group.standard = FALSE,
            y.gwas.used = "pheno",
            gen.architecture.m = 0,
            gen.architecture.f = NULL,
            add.architecture = NULL,
            ncore = 1,
            ncore.generation = 1,
            Z.integer = FALSE,
            store.effect.freq = FALSE,
            backend = "doParallel",
            randomSeed = NULL,
            randomSeed.generation = NULL,
            Rprof = FALSE,
            miraculix = NULL,
            miraculix.cores = 1,
            miraculix.mult = NULL,
            miraculix.chol = TRUE,
            best.selection.ratio.m = 1,
            best.selection.ratio.f = NULL,
            best.selection.criteria.m = "bv",
            best.selection.criteria.f = NULL,
            best.selection.manual.ratio.m = NULL,
            best.selection.manual.ratio.f = NULL,
            best.selection.manual.reorder = TRUE,
            bve.class = NULL,
            parallel.generation = FALSE,
            name.cohort = NULL,
            display.progress = TRUE,
            combine = FALSE,
            repeat.mating = NULL,
            repeat.mating.copy = NULL,
            repeat.mating.fixed = NULL,
            repeat.mating.overwrite = TRUE,
            time.point = 0,
            creating.type = 0,
            multiple.observation = FALSE,
            new.bv.observation = NULL,
            new.bv.observation.gen = NULL,
            new.bv.observation.cohorts = NULL,
            new.bv.observation.database = NULL,
            phenotyping = NULL,
            phenotyping.gen = NULL,
            phenotyping.cohorts = NULL,
            phenotyping.database = NULL,
            bve.gen = NULL,
            bve.cohorts = NULL,
            bve.database = NULL,
            sigma.e.gen = NULL,
            sigma.e.cohorts = NULL,
            sigma.e.database = NULL,
            sigma.g.gen = NULL,
            sigma.g.cohorts = NULL,
            sigma.g.database = NULL,
            gwas.gen = NULL,
            gwas.cohorts = NULL,
            gwas.database = NULL,
            bve.insert.gen = NULL,
            bve.insert.cohorts = NULL,
            bve.insert.database = NULL,
            reduced.selection.panel.m = NULL,
            reduced.selection.panel.f = NULL,
            breeding.all.combination = FALSE,
            depth.pedigree = 7,
            depth.pedigree.ogc = 7,
            copy.individual.keep.bve = TRUE,
            copy.individual.keep.pheno = TRUE,
            bve.avoid.duplicates = TRUE,
            report.accuracy = TRUE,
            share.genotyped = 1,
            singlestep.active = FALSE,
            remove.non.genotyped = TRUE,
            added.genotyped = 0,
            fast.uhat = TRUE,
            offspring.bve.parents.gen = NULL,
            offspring.bve.parents.database = NULL,
            offspring.bve.parents.cohorts = NULL,
            offspring.bve.offspring.gen = NULL,
            offspring.bve.offspring.database = NULL,
            offspring.bve.offspring.cohorts = NULL,
            culling.gen=NULL,
            culling.database=NULL,
            culling.cohort=NULL,
            culling.time = Inf,
            culling.name = "Not_named",
            culling.bv1 = 0,
            culling.share1 = 0,
            culling.bv2 = NULL,
            culling.share2 = NULL,
            culling.index = 0,
            culling.single = TRUE,
            culling.all.copy = TRUE,
            calculate.reliability=FALSE,
            selection.m.gen = NULL,
            selection.f.gen = NULL,
            selection.m.database = NULL,
            selection.f.database = NULL,
            selection.m.cohorts=NULL,
            selection.f.cohorts=NULL,
            selection.m.miesenberger=FALSE,
            selection.f.miesenberger=NULL,
            selection.miesenberger.reliability.est="estimated",
            miesenberger.trafo = 0,
            multiple.bve.weights.m = 1,
            multiple.bve.weights.f = NULL,
            multiple.bve.scale.m = "bv_sd",
            multiple.bve.scale.f = NULL,
            verbose=TRUE,
            bve.parent.mean=FALSE,
            bve.grandparent.mean=FALSE,
            bve.mean.between="bvepheno",
            bve.direct.est=TRUE,
            bve.pseudo=FALSE,
            bve.pseudo.accuracy=1,
            miraculix.destroyA=TRUE,
            mas.bve=FALSE,
            mas.markers=NULL,
            mas.number=5,
            mas.effects=NULL,
            threshold.selection=NULL,
            threshold.sign=">",
            input.phenotype="own",
            bve.ignore.traits=NULL,
            bv.ignore.traits=NULL,
            genotyped.database = NULL,
            genotyped.gen = NULL,
            genotyped.cohorts = NULL,
            genotyped.share = 1,
            genotyped.array = 1,
            sex.s = NULL,
            bve.imputation = TRUE,
            bve.imputation.errorrate = 0,
            share.phenotyped=1,
            avoid.mating.fullsib=FALSE,
            avoid.mating.halfsib=FALSE,
            max.mating.pair = Inf,
            bve.per.sample.sigma.e=TRUE,
            bve.solve = "exact"
            ){


  #######################################################################
  ## Pre-work to allow for flexiblity when inputing parameter values ####
  # Initialisize parameters that were not initialized in early versions #
  #######################################################################
{

  if(is.function(bve.solve) || bve.solve != "exact"){

    if(!is.function(bve.solve) && bve.solve =="pcg"){
      if (requireNamespace("cPCG", quietly = TRUE)) {
        bve.solve <- cPCG::pcgsolve
      } else{
        stop("Selected solver not available!")
      }

    }

  }

  if(length(population$info$default.parameter.name)>0){

    for(index in 1:length(population$info$default.parameter.name)){

      assign(population$info$default.parameter.name[index], value = population$info$default.parameter.value[[index]])

    }


  }
  if(parallel.generation){
    stop("Parallel generation has been enabled for now. Will come back in a later version!")
  }

  if(length(bv.ignore.traits)>0){
    temp123 <- setdiff(population$info$bv.random.activ , bv.ignore.traits)
  } else{
    temp123 <- population$info$bv.random.activ
  }

  {

  if(avoid.mating.halfsib){
    max_rel = 0
  } else if(avoid.mating.fullsib){
    max_rel = 1
  } else{
    max_rel = 2
  }

  if(repeat.mating.overwrite){
    repeat.mating.store <- population$info$repeat.mating
    repeat.mating.copy.store <- population$info$repeat.mating.copy
  }
  if(length(repeat.mating)>0){
    if(length(repeat.mating)==1){
      population$info$repeat.mating <- cbind(repeat.mating, 1)
    } else{
      population$info$repeat.mating <- repeat.mating
    }
    if(verbose & repeat.mating.overwrite) warning("New standard for litter size / repeat.mating set. This will be the new default for all downstream generation of offspring via breeding")
  }
  if(length(population$info$repeat.mating)==0){
    population$info$repeat.mating <- cbind(1,1)
  }

  if(length(repeat.mating.copy)>0){
    if(length(repeat.mating.copy)==1){
      population$info$repeat.mating.copy <- cbind(repeat.mating.copy, 1)
    } else{
      population$info$repeat.mating.copy <- repeat.mating.copy
    }
    if(verbose & repeat.mating.overwrite) warning("New standard for litter size / repeat.mating.copy set. This will be the new default for all downstream generation of offspring via copy/combine.")
  }
  if(length(population$info$repeat.mating.copy)==0){
    population$info$repeat.mating.copy <- cbind(1,1)
  }
  }


  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]])
      }
    }
  }



  if(Rprof){
    Rprof()
  }

  {
    # foreach variables
    indexb <- NULL
    keeps <- NULL
    pheno <- NULL
    if (requireNamespace("foreach", quietly = TRUE)) {
      `%dopar%` <- foreach::`%dopar%`
    }


  }


  if(sum(population$info$is.combi)>0){
    for(index in 1:length(population$info$combi.weights)){
      population$info$combi.weights[[index]] <- c(population$info$combi.weights[[index]], rep(0, population$info$bv.nr - length(population$info$combi.weights[[index]])))
    }
  }
  if(length(population$info$array.name)==0){
    population$info$array.name = "Full_Array"
    population$info$array.markers = list(rep(TRUE,sum(population$info$snp)))
    population$info$array.is_subset = FALSE
  }

  if(dh.mating){
    selfing.mating = TRUE
  }

  reduced.selection.panel <- list(reduced.selection.panel.m, reduced.selection.panel.f)

  if(length(selection.criteria)==0){
    selection.criteria <- c("bve", "bve")
    if(length(selection.m)==0){
      if(population$info$bv.nr>0){
        selection.m <- "function"
      } else{
        selection.m <- "random"
      }

    }
  } else{
    if(length(selection.m)==0){
      selection.m <- "function"
    }
  }
  if(culling.share1>0 || (length(culling.share2)>0 && culling.share2>0)){
    culling <- TRUE
  } else{
    culling <- FALSE
  }
  if(length(new.bv.child)>0){
    phenotyping.child <- new.bv.child
  }

  if(length(computation.A)>0){
    relationship.matrix <- computation.A
  }
  if(length(computation.A.ogc)){
    relationship.matrix.ogc <- computation.A.ogc
  }


  if(length(randomSeed)>0){
    set.seed(randomSeed)
  }

  if(length(bve.database)==0 && length(bve.gen)==0 && length(bve.cohorts)==0){
    bve.gen <- nrow(population$info$size)
  }

  if(length(population$info$phenotypic.transform)==0){
    population$info$phenotypic.transform <- rep(FALSE, population$info$bv.nr)
  }
  activ.trafo <- which(population$info$phenotypic.transform)


  # Fill databases

  genotyped.database <- get.database(population, genotyped.gen, genotyped.database, genotyped.cohorts)

  bve.gen.input <- bve.gen
  bve.database.input <- bve.database
  bve.cohorts.input <- bve.cohorts

  bve.database <- get.database(population, bve.gen, bve.database, bve.cohorts) # NOT DONE

  if(length(bve.insert.gen)==0 && length(bve.insert.cohorts)==0 && length(bve.insert.database)==0){
    bve.insert.gen <- bve.gen.input
    bve.insert.cohorts <- bve.cohorts.input
    bve.insert.database <- bve.database.input
  }
  bve.insert.database <- get.database(population, bve.insert.gen, bve.insert.database, bve.insert.cohorts) # NOT DONE

  # Add Individuals in bve.insert that are not in bve.database
  for(index in 1:nrow(bve.insert.database)){
    included <- (bve.insert.database[index,1] == bve.database[,1]) & (bve.insert.database[index,2] == bve.database[,2]) & (bve.insert.database[index,3] >= bve.database[,3]) & (bve.insert.database[index,4] <= bve.database[,4])
    if(sum(included)==0){
      if(verbose) cat("Not all individuals with breeding values to insert included in breeding value estimation.\n")
      if(verbose) cat("Missing individuals are automatically added in BVE.\n")
      bve.database <- rbind(bve.database, bve.insert.database[index,])
    }
  }
  if(length(sigma.e.gen)==0 && length(sigma.e.cohorts)==0 && length(sigma.e.database)==0){
    sigma.e.gen <- bve.gen.input
    sigma.e.cohorts <- bve.cohorts.input
    sigma.e.database <- bve.database.input
  }
  sigma.e.database <- get.database(population, sigma.e.gen, sigma.e.database, sigma.e.cohorts) # NOT DONE

  if(length(sigma.g.gen)==0 && length(sigma.g.cohorts)==0 && length(sigma.g.database)==0){
    sigma.g.gen <- bve.gen.input
    sigma.g.cohorts <- bve.cohorts.input
    sigma.g.database <- bve.database.input
  }
  sigma.g.database <- get.database(population, sigma.g.gen, sigma.g.database, sigma.g.cohorts) # NOT DONE

  if(length(gwas.gen)==0 && length(gwas.cohorts)==0 && length(gwas.database)==0){
    gwas.gen <- bve.gen.input
    gwas.cohorts <- bve.cohorts.input
    gwas.database <- bve.database.input
  }
  gwas.database <- get.database(population, gwas.gen, gwas.database, gwas.cohorts) # NOT DONE

  if(length(new.bv.observation)>0){
    phenotyping <- new.bv.observation
  }
  if(length(new.bv.observation.gen)>0){
    phenotyping.gen <- new.bv.observation.gen
  }
  if(length(new.bv.observation.cohorts)>0){
    phenotyping.cohorts <- new.bv.observation.cohorts
  }
  if(length(new.bv.observation.database)>0){
    phenotyping.database <- new.bv.observation.database
  }
  if(length(phenotyping)==1 && (phenotyping=="all" || phenotyping=="non_obs" )){
    phenotyping.gen <- 1:length(population$breeding)
    if(phenotyping=="non_obs"){
      multiple.observation <- FALSE
    }
  }
  if(length(phenotyping)==1 && phenotyping=="non_obs_m"){
    phenotyping.database <- cbind(1:length(population$breeding),1)
  }
  if(length(phenotyping)==1 && phenotyping=="non_obs_f"){
    phenotyping.database <- cbind(1:length(population$breeding),2)
  }

  phenotyping.database <- get.database(population, phenotyping.gen, phenotyping.database, phenotyping.cohorts)

  if((copy.individual.m + copy.individual.f + combine)>1){
    stop("Use of multiple copy parameter at the same time is forbidden!")
  }

  if(length(selection.size)==1){
    selection.size <- rep(selection.size,2)
  }


  if(copy.individual.m){
    copy.individual <- TRUE
    selfing.mating <- TRUE
    selfing.sex <- 0
    if(sum(breeding.size)==0){
      breeding.size <- c(selection.size[1],0)
    }
    if(selection.size[1]>=breeding.size[1]){
      max.offspring <- c(1,1)
    }
  }

  if(copy.individual.f){
    copy.individual <- TRUE
    selfing.mating <- TRUE
    selfing.sex <- 1
    if(sum(breeding.size)==0){
      breeding.size <- c(0,selection.size[2])
    }
    if(selection.size[2]>=breeding.size[2]){
      max.offspring <- c(1,1)
    }
  }


  if(combine==TRUE){
    # combine is modelled via cloning with no recombination
    copy.individual <- TRUE
    selfing.mating <- TRUE
    if(selection.size[2]>0 & selection.size[1]==0){
      selfing.sex <- 1
    } else{
      selfing.sex <- 0
    }

    class.m <- unique(c(class.m, class.f))
    best1.from.cohort <- c(best1.from.cohort, best2.from.cohort)
    best2.from.cohort <- NULL
    best1.from.group <- c(best1.from.group, best2.from.group)
    best2.from.group <- NULL
    max.offspring = c(1,1)
  }


  if(length(n.observation)==0){
    if(copy.individual){
      n.observation <- 0L
    } else{
      n.observation <- 1L
    }

  }
  if(length(n.observation)>0){
    n.observation <- as.integer(n.observation)
  }
  if(length(n.observation)<population$info$bv.nr){
    n.observation <- rep(n.observation, length.out=population$info$bv.nr)

  }

  if(length(phenotyping.child)==0){
    if(copy.individual){
      phenotyping.child <- "addobs"
    } else{
      phenotyping.child <- "zero"
    }

  }


  if(length(best1.from.cohort)){
    if(length(selection.m.cohorts)==0){
      selection.m.cohorts <- best1.from.cohort
    } else{
      selection.m.cohorts <- c(selection.m.cohorts,best1.from.cohort)
    }
  }
  if(length(best2.from.cohort)){
    if(length(selection.m.cohorts)==0){
      selection.f.cohorts <- best2.from.cohort
    } else{
      selection.f.cohorts <- c(selection.f.cohorts,best2.from.cohort)
    }
  }
  if(length(best1.from.group)){
    if(length(selection.m.database)==0){
      selection.m.database <- best1.from.group
    } else{
      selection.m.database <- c(selection.m.database,best1.from.group)
    }
  }
  if(length(best2.from.group)){
    if(length(selection.f.database)==0){
      selection.f.database <- best2.from.group
    } else{
      selection.f.database <- c(selection.f.database,best2.from.group)
    }
  }

  if(population$info$bve && (!population$info$bv.calculated || length(  population$info$bv.random.activ)==0 || sum(population$info$is.combi)>0)){
    population$info$bv.random.activ <- which(population$info$bv.random[1:population$info$bv.calc]==FALSE & population$info$is.combi[1:population$info$bv.calc]==FALSE)
  }

  add.selection <- sum(breeding.size)>0 & sum(selection.size)==0
  if(length(selection.m.gen)==0 && length(selection.m.database)==0 && length(selection.m.cohorts)==0 && (selection.size[1]>0 || add.selection) &&
     length(fixed.breeding)==0 ){

    if(sum(population$info$size[,1])>0){
      if(verbose) cat("No individuals for selection provided (male side). Use last available.\n")
      selection.m.database <- cbind(max(which(population$info$size[,1]>0)),1)
    } else{
      if(verbose) cat("No individuals for selection provided (male side). Non available.\n")
    }

  }
  if(length(selection.f.gen)==0 && length(selection.f.database)==0 && length(selection.f.cohorts)==0 && (selection.size[2]>0 ||add.selection)&&
     length(fixed.breeding)==0){
    if(sum(population$info$size[,2])>0){
      if(verbose) cat("No individuals for selection provided (female side). Use last available.\n")
      selection.f.database <- cbind(max(which(population$info$size[,2]>0)),2)
    } else{
      if(verbose) cat("No individuals for selection provided (female side). Non available.\n")
    }

  }

  selection.m.database <- get.database(population, selection.m.gen, selection.m.database, selection.m.cohorts)
  selection.f.database <- get.database(population, selection.f.gen, selection.f.database, selection.f.cohorts)

  if(add.gen==0){
    current.gen <- length(population$breeding)
    add.gen <- current.gen + 1
  } else{
    current.gen <- add.gen - 1
  }

  if(!copy.individual && (length(selection.m.database)>0 && selection.m.database[,1]>= add.gen && (sum(breeding.size)+ sum(selection.size))>0)){
    stop("Parental individuals must be in an earlier generation than offspring! Check add.gen / selection.m.gen/database/cohort!")
  }
  if(!copy.individual && (length(selection.f.database)>0 && selection.f.database[,1]>= add.gen && (sum(breeding.size)+ sum(selection.size))>0)){
    stop("Parental individuals must be in an earlier generation than offspring! Check add.gen / selection.m.gen/database/cohort!")
  }
  if(length(population$info$cumsnp)==0){
    population$info$cumsnp <- cumsum(population$info$snp)
  }
  if(length(gen.architecture.f)==0){
    gen.architecture.f=gen.architecture.m
  }
  if(length(selection.criteria)==1){
    selection.criteria <- rep(selection.criteria,2)
  }
  if(selection.criteria[1]=="random"){
    selection.m = "random"
  }
  if(selection.criteria[2]=="random"){
    selection.f = "random"
  }
  if(length(max.offspring)==1){
    max.offspring <- rep(max.offspring,2)
  }
  if(length(max.litter)==1){
    max.litter <- rep(max.litter,2)
  }
  best.selection.ratio <- list()
  best.selection.ratio[[1]] <- best.selection.ratio.m
  if(length(best.selection.ratio.f)==0){
    best.selection.ratio.f <- best.selection.ratio.m
  }
  best.selection.ratio[[2]] <- best.selection.ratio.f

  best.selection.criteria <- list()
  best.selection.criteria[[1]] <- best.selection.criteria.m
  if(length(best.selection.criteria.f)==0){
    best.selection.criteria.f <- best.selection.criteria.m
  }
  best.selection.criteria[[2]] <- best.selection.criteria.f

  best.selection.manual.ratio <- list()
  best.selection.manual.ratio[[1]] <- best.selection.manual.ratio.m
  if(length(best.selection.manual.ratio.f)==0 && length(best.selection.manual.ratio.f)!=0 ){
    best.selection.manual.ratio.f <- best.selection.manual.ratio.m
  }
  best.selection.manual.ratio[[2]] <- best.selection.manual.ratio.f
  best.selection.manual.ratio[[3]] <- "placeholder"

  if(length(population$info$origin.gen)>0){
    population$info$origin.gen <- base::as.integer(population$info$origin.gen)
  } else if(population$info$miraculix){
    population$info$origin.gen <- 1:32L
  } else{
    population$info$origin.gen <- 1:64L
  }

  if(bve.pseudo && length(bve.pseudo.accuracy) < population$info$bv.nr){
    bve.pseudo.accuracy <- rep(bve.pseudo.accuracy, length.out=population$info$bv.nr)
  }



  if(length(population$info$origin)==0){
    population$info$origin <- 1:64
  }


  if(sum(population$info$snp)<=maxZ){
    sequenceZ <- FALSE
  }
  if(length(miraculix)==0){
    miraculix <- population$info$miraculix
    if(length(miraculix)==0){
      stop("No input on use of miraculix provided!")
    }
  }
  if(miraculix && length(miraculix.mult)==0){
    miraculix.mult <- TRUE
  } else if(length(miraculix.mult)==0){
    miraculix.mult <- FALSE
  }
  if(store.comp.times){
    comp.times <- numeric(7)
    comp.times[1] <- as.numeric(Sys.time())
  }

  if(gene.editing.offspring || gene.editing.best){
    gene.editing <- TRUE
  } else{
    gene.editing <- FALSE
  }
  if(length(population$info$bitstoring)>0){
    bit.storing <- TRUE
    nbits <- population$info$bitstoring
    if(maxZ%%nbits!=0){
      maxZ <- maxZ + nbits - maxZ%%nbits
    }
  } else{
    nbits <- 30
    bit.storing <- FALSE
  }
  # Increase maxZ to improve run-time based on bitwise computing.
  if(miraculix && maxZ%%512!=0){
    maxZ <- maxZ + 512 - maxZ%%512
  }

  if(length(add.architecture)>0){
    population$info$gen.architecture[[length(population$info$gen.architecture)+1]] <- list()
    population$info$gen.architecture[[length(population$info$gen.architecture)]]$length.total <- cumsum(c(0,add.architecture[[1]]))
    population$info$gen.architecture[[length(population$info$gen.architecture)]]$snp.position <- add.architecture[[2]]

  }
  if (requireNamespace("miraculix", quietly = TRUE)) {
    if(miraculix.cores>1){
      RandomFieldsUtils::RFoptions(cores = miraculix.cores)
    }
    codeOriginsU <- miraculix::codeOrigins
    decodeOriginsU <- miraculix::decodeOrigins
  } else{
    codeOriginsU <- codeOriginsR
    decodeOriginsU <- decodeOriginsR
  }

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



  if(length(new.residual.correlation)>0){
    population$info$pheno.correlation <- t(chol(new.residual.correlation))
  }
  if(length(new.breeding.correlation)>0){
    population$info$bv.correlation <- new.breeding.correlation
  }


  class <- list()
  class[[1]] <- class.m
  class[[2]] <- class.f
  if(length(sigma.e)==0 && length(population$info$last.sigma.e.value)>0){
    sigma.e <- population$info$last.sigma.e.value
  } else{
    sigma.e <- 100
  }

  if(length(sigma.e)==1){
    sigma.e <- rep(sigma.e, population$info$bv.nr)
  }
  if(length(sigma.g)==1){
    sigma.g <- rep(sigma.g, population$info$bv.nr)
  }

  if(store.comp.times){
    comp.times[2] <- as.numeric(Sys.time())
  }




  if(length(selection.function.matrix)!=0 && is.matrix(selection.function.matrix)==FALSE){
    selection.function.matrix <- matrix(selection.function.matrix,nrow=1)
  }
  if(length(selection.f)==0){
    selection.f <- selection.m
  }
  if(length(multiple.bve.weights.m)< population$info$bv.nr){
    multiple.bve.weights.m <- rep(multiple.bve.weights.m, length.out = population$info$bv.nr)
  }
  if(length(multiple.bve.weights.f)==0){
    multiple.bve.weights.f <- multiple.bve.weights.m
  }
  if(length(bve.ignore.traits)==1 && bve.ignore.traits=="zero"){
    bve.ignore.traits <- which(multiple.bve.weights.m==0 & multiple.bve.weights.f==0)
  }
  if(length(bve.ignore.traits)==0){
    bve.keeps <- 1:population$info$bv.nr
  } else{
    bve.keeps <- (1:population$info$bv.nr)[-bve.ignore.traits]
  }
  if(length(multiple.bve.weights.f)< population$info$bv.nr){
    multiple.bve.weights.f <- rep(multiple.bve.weights.f, length.out = population$info$bv.nr)
  }
  if(length(multiple.bve.scale.f)==0){
    multiple.bve.scale.f <- multiple.bve.scale.m
  }

  if(length(selection.f.miesenberger)==0){
    selection.f.miesenberger <- selection.m.miesenberger
  }


  if(length(ignore.best)==1){
    ignore.best <- rep(ignore.best,2)
  }
  if(length(breeding.sex)==0 && length(breeding.size)==2){
    breeding.sex <- breeding.size[2] / sum(breeding.size)
  } else if(length(breeding.sex)==0){
    breeding.sex <- 0.5
  }
  if(length(breeding.size)==1){
    breeding.size.total <- breeding.size
    breeding.size.m <- round(breeding.size*(1-breeding.sex))
    if(breeding.sex.random){
      breeding.size.m <- sum(stats::rbinom(breeding.size,1,(1-breeding.sex)))
    }
    breeding.size <- c(breeding.size.m, breeding.size - breeding.size.m)
  } else{
    breeding.size.total <- sum(breeding.size)
  }


  if(length(sex.s)==0){
    sex.animal <- rep(2, breeding.size.total)
    if(breeding.size[1] > 0){
      sex.animalr <- sample(1:breeding.size.total, breeding.size[1])
      sex.animal[sex.animalr] <- 1
    }
  } else{
    sex.animal <- rep(sex.s, length.out = breeding.size.total)
    if(sum(sex.animal==2) != breeding.size[2]){
      stop("Chosen sex ratio in sex.s does not match with breeding.size")
    }
  }


  if(breeding.all.combination){
    if(sum(selection.size)==0){
      selection.size[1] <- sum(selection.m.database[,4]-selection.m.database[,3]+1)
      selection.size[2] <- sum(selection.f.database[,4]-selection.f.database[,3]+1)
    }
    if(selection.size[1]>0 && selection.size[2]>0){
      fixed.breeding.best <- matrix(0, nrow= selection.size[1] * selection.size[2], ncol=4)
      for(index in 1:selection.size[1]){
        fixed.breeding.best[1:selection.size[2] + (index-1) * selection.size[2],] <- cbind(1,index, 2, 1:selection.size[2])
      }
    } else{
      nindi <- max(selection.size[1:2])
      sex.fixed <- 1 + as.numeric(selection.size[2]>0)
      fixed.breeding.best <- matrix(0, nrow= nindi * (nindi-1) / 2, ncol=4)
      prior <- 0
      for(index in 1:(nindi-1)){
        fixed.breeding.best[1:(nindi-index)+ prior,] <- cbind(sex.fixed,index, sex.fixed, (index+1):nindi)
        prior <- prior +nindi-index
      }
    }
  }

  if(length(fixed.breeding) >0 && ncol(fixed.breeding)==6){
    fixed.breeding <- cbind(fixed.breeding,breeding.sex)
  }
  if(length(fixed.breeding.best) >0 && ncol(fixed.breeding.best)==4){
    fixed.breeding.best <- cbind(fixed.breeding.best,breeding.sex)
  }

  if(copy.individual){
    repeat.mating.activ <- population$info$repeat.mating.copy
  } else{
    repeat.mating.activ <- population$info$repeat.mating
  }



  if(nrow(repeat.mating.activ)==1 && repeat.mating.activ[1,1]>1 && length(fixed.breeding)>0){
    fixed.breeding <- matrix(rep(t(fixed.breeding), repeat.mating.activ[1,1]), ncol=7, byrow=TRUE)
  }
  if(nrow(repeat.mating.activ)==1 && repeat.mating.activ[1,1]>1 && length(fixed.breeding.best)>0){
    fixed.breeding.best <- matrix(rep(t(fixed.breeding.best), repeat.mating.activ[1,1]), ncol=5, byrow=TRUE)
  }



}

  #######################################################################
  ############## Start of the actual Simulation #########################
  #######################################################################
{

  # Check if underlying true genomic values need to be calculated ((usually only on first use of breeding.diploid() or when traits are added))

  if(population$info$bv.calculated==FALSE || length(population$info$effect.p)==0 || (length(population$info$effect.p.add)==0 && length(population$info$effect.p.mult1)==0 && length(population$info$effect.p.dice)==0)){
    excludes <- NULL
    snp.before <- cumsum(c(0,population$info$snp))
    population$info$effect.p.add <- list()
    if(population$info$real.bv.length[1]>0){
      for(index in 1:population$info$real.bv.length[1]){
        if(length(population$info$real.bv.add[[index]])>0 && nrow(population$info$real.bv.add[[index]])>0){
          population$info$effect.p.add[[index]] <- as.integer(population$info$real.bv.add[[index]][,1]+ snp.before[population$info$real.bv.add[[index]][,2]])
          excludes <- unique(c(excludes, population$info$effect.p.add))

        }
      }
    }

    if(population$info$real.bv.length[2]>0){
      population$info$effect.p.mult1 <- list()
      population$info$effect.p.mult2 <- list()
      for(index in 1:population$info$real.bv.length[2]){
        if(length(population$info$real.bv.mult[[index]])>0 &&nrow(population$info$real.bv.mult[[index]])>0){
          population$info$effect.p.mult1[[index]] <- as.integer(population$info$real.bv.mult[[index]][,1]+ snp.before[population$info$real.bv.mult[[index]][,2]])
          population$info$effect.p.mult2[[index]] <- as.integer(population$info$real.bv.mult[[index]][,3]+ snp.before[population$info$real.bv.mult[[index]][,4]])
          excludes <- unique(c(excludes, population$info$effect.p.mult1,population$info$effect.p.mult2))
        }
      }
    }
    if(population$info$real.bv.length[3]>0){

      for(index in 1:population$info$real.bv.length[3]){
        if(length(population$info$real.bv.dice[[index]])>0){
          for(index2 in 1:length(population$info$real.bv.dice[[index]][[1]])){
            population$info$effect.p.dice <- as.integer(c(population$info$effect.p.dice, population$info$real.bv.dice[[index]][[1]][[index2]][,1]+ snp.before[population$info$real.bv.dice[[index]][[1]][[index2]][,2]]))
          }
          excludes <- unique(c(excludes, population$info$effect.p.dice))
        }
      }
    }

    population$info$effect.p <- as.integer(unlist(excludes))

    population$info$effect.p.add.same <- rep(FALSE, population$info$bv.nr)
    if(population$info$real.bv.length[1]>1){
      for(index in 2:population$info$real.bv.length[1]){
        if(length(population$info$effect.p.add)>=index && (length(population$info$effect.p.add[[index]]) == length(population$info$effect.p.add[[index-1]]) &&
           length(population$info$effect.p.add[[index]]) > 0 &&
           prod(population$info$effect.p.add[[index]] == population$info$effect.p.add[[index-1]]) == 1)){
          population$info$effect.p.add.same[index] <- TRUE
        }
      }
    }

  }

  temp1 <- population$info$bv.calculated
  if(population$info$bve && population$info$bv.calculated==FALSE){
    if(verbose) cat("Derive genomic values of founders. \n")


    for(index in 1:length(population$breeding)){
      for(sex in 1:2){
        nanimals <- length(population$breeding[[index]][[sex]])
        if(nanimals >0){
          for(nr.animal in 1:nanimals){
            activ_bv <- population$info$bv.random.activ
            if(length(activ_bv)>0){
              temp_out <- calculate.bv(population, index, sex, nr.animal,
                                       activ_bv, import.position.calculation=import.position.calculation,
                                       decodeOriginsU=decodeOriginsU,
                                       store.effect.freq=store.effect.freq,
                                       bit.storing=bit.storing, nbits=nbits, output_compressed=FALSE,
                                       bv.ignore.traits=bv.ignore.traits)
              population$breeding[[index]][[6+sex]][activ_bv,nr.animal] <- temp_out[[1]]
              population$breeding[[index]][[sex]][[nr.animal]][[25]] <- length(bv.ignore.traits)==0

              if(length(temp123)>0){
                population$breeding[[index]][[sex]][[nr.animal]][[26]] <- temp123
              }

              if(store.effect.freq){
                if(length(population$info$store.effect.freq) < index || length(population$info$store.effect.freq[[index]])==0){
                  colnames(temp_out[[2]]) <- c("Homo0", "Hetero", "Homo1")
                  rownames(temp_out[[2]]) <- population$info$snp.name[population$info$effect.p]
                  population$info$store.effect.freq[[index]] <- temp_out[[2]]
                } else{
                  population$info$store.effect.freq[[index]] <- population$info$store.effect.freq[[index]] + temp_out[[2]]
                }
              }
            }
          }
        }
      }
    }
    population$info$bv.calculated <- TRUE
  }

  # Calculate breeding values for non-QTL-based traits
  if(population$info$bv.calc>0 && population$info$bv.random[population$info$bv.calc]){
    mu1 <- numeric((population$info$bv.calc-1))
    # Estimate sigma.g for QTL-based traits ((correlation))
    if(population$info$bv.calc>1){

      n.animals <- 0
      for(index in 1:nrow(sigma.g.database)){
        n.animals <- n.animals + diff(sigma.g.database[index,3:4]) + 1
      }
      y_real <- array(0, dim=c(n.animals,(population$info$bv.nr))) # schaetzung sigma.g
      cindex <- 1
      temp2 <- 1:(population$info$bv.nr)
      temp3 <- which(population$info$bv.random==FALSE)
      for(index in 1:nrow(sigma.g.database)){
        k.database <- sigma.g.database[index,]
        if(diff(k.database[3:4])>=0){
          for(kindex in k.database[3]:k.database[4]){
            y_real[cindex,temp2] <- population$breeding[[k.database[[1]]]][[6+k.database[[2]]]][temp2, kindex]
            cindex <- cindex +1
          }
        }
      }
      for(bven in temp3){
        population$info$bv.random.variance[bven] <- stats::var(y_real[,bven])
        mu1[bven] <- mean(y_real[,bven])
      }

    }
    population$info$current.bv.correlation <- population$info$bv.correlation
    if(sum(is.na(population$info$bv.correlation))>0){
      emp_cor <- which(is.na(population$info$bv.correlation), arr.ind=TRUE)
      for(index in 1:nrow(emp_cor)){
        population$info$current.bv.correlation[emp_cor[index,1], emp_cor[index,2]] <-
          stats::cor(y_real[,emp_cor[index,1]], y_real[,emp_cor[index,2]])
      }
    }
    if(temp1==FALSE){
      population$info$current.bv.random.variance  <- population$info$bv.random.variance
      if(population$info$bv.calc==1){
        bv.var <- diag.mobps(sqrt(population$info$current.bv.random.variance)) %*%population$info$current.bv.correlation %*% diag.mobps(sqrt(population$info$current.bv.random.variance))
      } else{
        AA <- diag.mobps(sqrt(population$info$current.bv.random.variance)[1:(population$info$bv.calc-1)]) %*% population$info$current.bv.correlation[1:(population$info$bv.calc-1), 1:(population$info$bv.calc-1)]%*% diag.mobps(sqrt(population$info$current.bv.random.variance)[(1:(population$info$bv.calc-1))])
        BB <- diag.mobps(sqrt(population$info$current.bv.random.variance)[1:(population$info$bv.calc-1)]) %*%population$info$current.bv.correlation[1:(population$info$bv.calc-1), -(1:(population$info$bv.calc-1))]%*% diag.mobps(sqrt(population$info$current.bv.random.variance)[-(1:(population$info$bv.calc-1))])
        CC <- diag.mobps(sqrt(population$info$current.bv.random.variance)[-(1:(population$info$bv.calc-1))]) %*%population$info$current.bv.correlation[-(1:(population$info$bv.calc-1)), -(1:(population$info$bv.calc-1))] %*% diag.mobps(sqrt(population$info$current.bv.random.variance)[-(1:(population$info$bv.calc-1))])
        if (requireNamespace("MASS", quietly = TRUE)) {
          bv.var <- CC - t(BB) %*% MASS::ginv(AA) %*% BB
          mu_mult <- t(BB) %*% MASS::ginv(AA)
        } else{
          bv.var <- CC - t(BB) %*% solve(AA) %*% BB
          mu_mult <- t(BB) %*% solve(AA)
        }
      }
      bv.var_chol <- t(chol(bv.var))
      population$info$bv.correlation_col <- bv.var_chol
      for(index in 1:length(population$breeding)){
        for(sex in 1:2){
          nanimals <- length(population$breeding[[index]][[sex]])
          if(nanimals >0){
            for(nr.animal in 1:nanimals){
              if(population$info$bv.calc==1){
                bv.mean <- population$info$base.bv
              } else{
                bv.mean <- population$info$base.bv[-(1:(population$info$bv.calc-1))] +  mu_mult %*% (population$breeding[[index]][[sex+6]][1:(population$info$bv.calc-1), nr.animal] - mu1)
              }
              population$breeding[[index]][[6+sex]][population$info$bv.calc:population$info$bv.nr, nr.animal] <- bv.mean + bv.var_chol %*% stats::rnorm(population$info$bv.nr-population$info$bv.calc+1,0,1)
            }
          }
        }
      }
    }

  }

  if(temp1==FALSE && sum(population$info$is.combi)>0){
    combis <- which(population$info$is.combi)
    for(index in 1:length(population$breeding)){
      for(sex in 1:2){
        nanimals <- length(population$breeding[[index]][[sex]])
        if(nanimals > 0){
          for(combi in combis){
            population$breeding[[index]][[6+sex]][combi,] <- colSums(population$info$combi.weights[[combi]] * population$breeding[[index]][[6+sex]])
          }
        }
      }
    }

  }
  if(store.comp.times){
    comp.times[3] <- as.numeric(Sys.time())
  }


  # Compute sigma.e to fullfil a target heritablity in the reference population
  if(length(heritability)>0){

    if(length(population$info$last.sigma.e.database)==0 || !(nrow(population$info$last.sigma.e.database)==nrow(sigma.e.database) && prod(population$info$last.sigma.e.database==sigma.e.database)==1) || length(heritability)>0 && (length(population$info$last.sigma.e.heritability)==0 || prod(population$info$last.sigma.e.heritability==heritability)==0)){
      if(verbose) cat("Start deriving enviromental variance (according to given heritability).\n")
      if(length(heritability)!=population$info$bv.nr){
        heritability <- rep(heritability, length.out = population$info$bv.nr)
      }
      n.animals <- 0
      for(index in 1:nrow(sigma.e.database)){
        n.animals <- n.animals + diff(sigma.e.database[index,3:4]) + 1
      }
      y_real <- array(0, dim=c(n.animals,(population$info$bv.nr))) # schaetzung sigma.g

      cindex <- 1
      temp1 <- 1:(population$info$bv.nr)
      for(index in 1:nrow(sigma.e.database)){
        k.database <- sigma.e.database[index,]
        if(diff(k.database[3:4])>=0){
          for(kindex in k.database[3]:k.database[4]){
            y_real[cindex,temp1] <- population$breeding[[k.database[1]]][[6+k.database[2]]][temp1,kindex]
            cindex <- cindex +1
          }
        }
      }

      for(bven in 1:population$info$bv.nr){
        if(forecast.sigma.g){
          sigma.g.temp <- stats::var(y_real[,bven])
        }
        sigma.e[bven] <- ((1- heritability[bven]) * sigma.g.temp)/ heritability[bven]
      }
    } else{
      sigma.e <- population$info$last.sigma.e.value
    }

  } else{
    if(length(population$info$last.sigma.e.value)>0){
      sigma.e <- population$info$last.sigma.e.value
    }
  }

  if(length(heritability)==0){
    heritability <- population$info$last.sigma.e.heritability
  }

  if(length(repeatability)>0){
    if(length(repeatability)!= population$info$bv.nr){
      repeatability <- rep(repeatability, length.out = population$info$bv.nr)
    }

    repeatability[repeatability<heritability] <-  heritability[repeatability<heritability]
    population$info$repeatability <- repeatability
  } else{

    if(length(population$info$repeatability)==0){
      population$info$repeatability <- heritability
    }

    repeatability <- population$info$repeatability

  }

  if(length(heritability)>0 & length(repeatability)>0){
    sigma.g.temp1 <- (sigma.e * heritability) /  (1 - heritability)
    sigma.total.temp1 <-  sigma.g.temp1 + sigma.e
    sigma.e.perm <- repeatability * sigma.total.temp1 - sigma.g.temp1
    sigma.e.rest <- sigma.e - sigma.e.perm

    sigma.e.perm[sigma.e.perm<0] <- 0
    sigma.e.rest[sigma.e.rest<0] <- 0
  } else{
    sigma.e.perm <- rep(0, population$info$bv.nr)
    sigma.e.rest <- sigma.e
  }

  # Genotyping

  if(length(genotyped.database)>0){
    for(index in 1:nrow(genotyped.database)){
      if((genotyped.database[index,4]-genotyped.database[index,3])>=0){
        for(index2 in genotyped.database[index,3]:genotyped.database[index,4]){
          temp1 <- stats::rbinom(1, 1, genotyped.share)
          population$breeding[[genotyped.database[index,1]]][[genotyped.database[index,2]]][[index2]][[16]] <- max(population$breeding[[genotyped.database[index,1]]][[genotyped.database[index,2]]][[index2]][[16]], stats::rbinom(1, 1, genotyped.share))
          if(temp1==1){
            population$breeding[[genotyped.database[index,1]]][[genotyped.database[index,2]]][[index2]][[22]] <-
              c(population$breeding[[genotyped.database[index,1]]][[genotyped.database[index,2]]][[index2]][[22]], genotyped.array)
          }

        }
      }
    }
  }

  # Phenotyping

  if(length(phenotyping.database)>0 && population$info$bve && sum(n.observation)>0){

    if(verbose) cat("Start simulating phenotypes.\n")

    for(index in 1:nrow(phenotyping.database)){
      gen <- phenotyping.database[index,1]
      sex <- phenotyping.database[index,2]
      if(diff(phenotyping.database[index,3:4])>=0){

        if(share.phenotyped<1){
          to_phenotype <- sample(phenotyping.database[index,3]:phenotyping.database[index,4], share.phenotyped *
                                   length(phenotyping.database[index,3]:phenotyping.database[index,4]))
        } else{
          to_phenotype <- phenotyping.database[index,3]:phenotyping.database[index,4]
        }
        for(nr.animal in to_phenotype){
          if( length(population$breeding[[gen]][[sex]][[nr.animal]][[23]]) == 0){
            population$breeding[[gen]][[sex]][[nr.animal]][[23]] <- stats::rnorm(population$info$bv.nr,0,1)
          }
          multi_check <- (population$breeding[[gen]][[sex]][[nr.animal]][[15]]==0)

          n.observation_temp <- n.observation * (multi_check | multiple.observation)
          population$breeding[[gen]][[sex]][[nr.animal]][[15]] <- population$breeding[[gen]][[sex]][[nr.animal]][[15]] + n.observation_temp

          obsmax <- max(population$breeding[[gen]][[sex]][[nr.animal]][[15]])
          if(length(population$breeding[[gen]][[sex]][[nr.animal]][[24]])==0 ||
             obsmax > ncol(population$breeding[[gen]][[sex]][[nr.animal]][[24]]) &&
             obsmax > 0){

            population$breeding[[gen]][[sex]][[nr.animal]][[24]] <- cbind(population$breeding[[gen]][[sex]][[nr.animal]][[24]],
                                                                          matrix(stats::rnorm(obsmax * population$info$bv.nr - length(population$breeding[[gen]][[sex]][[nr.animal]][[24]]),0,1),
                                                                                 nrow = population$info$bv.nr))
          }

          if(!population$breeding[[gen]][[sex]][[nr.animal]][[25]]){

            if(length(population$breeding[[gen]][[sex]][[nr.animal]][[26]])==0 || sum(n.observation_temp[-population$breeding[[gen]][[sex]][[nr.animal]][[26]]]>0)>0){

              to_ignore_temp <- which(n.observation_temp==0)

              activ_bv <- population$info$bv.random.activ


              temp1234 <- setdiff(population$info$bv.random.activ , to_ignore_temp)

              temp_out <- calculate.bv(population, gen, sex, nr.animal,
                                       activ_bv, import.position.calculation=import.position.calculation,
                                       decodeOriginsU=decodeOriginsU,
                                       store.effect.freq=store.effect.freq,
                                       bit.storing=bit.storing, nbits=nbits, output_compressed=FALSE,
                                       bv.ignore.traits=to_ignore_temp)

              if(length(population$breeding[[gen]][[sex]][[nr.animal]][[26]])>0){
                activ_replace <- !duplicated(c(population$breeding[[gen]][[sex]][[nr.animal]][[26]], activ_bv))[-(1:length(population$breeding[[gen]][[sex]][[nr.animal]][[26]]))]
              } else{
                activ_replace <- rep(TRUE, length(activ_bv))
              }

              population$breeding[[gen]][[6+sex]][activ_bv[activ_replace],nr.animal] <- temp_out[[1]][activ_replace]

              if(length(c(population$breeding[[gen]][[sex]][[nr.animal]][[26]], temp1234))>0){
                population$breeding[[gen]][[sex]][[nr.animal]][[26]] <- c(population$breeding[[gen]][[sex]][[nr.animal]][[26]], temp1234)
              }

              if(length(population$breeding[[gen]][[sex]][[nr.animal]][[26]])==population$info$bv.nr){
                population$breeding[[gen]][[sex]][[nr.animal]][[25]] <- TRUE
              }


            }



          }
          for(bven in intersect(which(n.observation_temp>0), setdiff(1:population$info$bv.nr, activ.trafo))){
            if(population$breeding[[gen]][[sex]][[nr.animal]][[15]][bven]>=1){
              population$breeding[[gen]][[8+sex]][bven, nr.animal] <- (sqrt(sigma.e.rest) * rowMeans(population$info$pheno.correlation %*% population$breeding[[gen]][[sex]][[nr.animal]][[24]][,1:population$breeding[[gen]][[sex]][[nr.animal]][[15]][bven]]) +
                                                                         sqrt(sigma.e.perm) * population$info$pheno.correlation %*% population$breeding[[gen]][[sex]][[nr.animal]][[23]] +
                                                                         population$breeding[[gen]][[6+sex]][, nr.animal])[bven]
            }

          }
          for(bven in intersect( which(n.observation_temp>0), intersect(1:population$info$bv.nr, activ.trafo))){
            if(population$breeding[[gen]][[sex]][[nr.animal]][[15]][bven]>=1){
              new_pheno <- (sqrt(sigma.e.rest) * population$info$pheno.correlation %*% population$breeding[[gen]][[sex]][[nr.animal]][[24]][,1:population$breeding[[gen]][[sex]][[nr.animal]][[15]][bven]])[bven,] +
                              (sqrt(sigma.e.perm) * population$info$pheno.correlation %*% population$breeding[[gen]][[sex]][[nr.animal]][[23]])[bven] +
                              population$breeding[[gen]][[6+sex]][bven, nr.animal]

              population$breeding[[gen]][[8+sex]][bven, nr.animal] <- mean(population$info$phenotypic.transform.function[[bven]](new_pheno))
            }
          }

        }
      }
    }
  }


  # Import offspring phenotypes for parents
  offspring.bve.parents.database <- get.database(population, offspring.bve.parents.gen, offspring.bve.parents.database, offspring.bve.parents.cohorts)
  if(length(offspring.bve.parents.database)>0){
    offspring.bve <- TRUE
  } else{
    offspring.bve <- FALSE
  }

  if(length(offspring.bve.offspring.gen)>0 || length(offspring.bve.offspring.database)>0 || length(offspring.bve.offspring.cohorts)>0){
    offspring.bve.offspring.database <- get.database(population, offspring.bve.offspring.gen, offspring.bve.offspring.database, offspring.bve.offspring.cohorts)

    first_gen <- min(offspring.bve.parents.database[,1])
    list_of_copy <- list()
    for(gen_check in 1:nrow(offspring.bve.parents.database)){
      for(gen_check2 in offspring.bve.parents.database[gen_check,3]:offspring.bve.parents.database[gen_check,4]){
        first_gen <- min(first_gen, population$breeding[[offspring.bve.parents.database[gen_check,1]]][[offspring.bve.parents.database[gen_check,2]]][[gen_check2]][[21]][,1])
        list_of_copy[[length(list_of_copy)+1]] <- rbind(offspring.bve.parents.database[gen_check,1], offspring.bve.parents.database[gen_check,2], gen_check2 ,
                                                        t(population$breeding[[offspring.bve.parents.database[gen_check,1]]][[offspring.bve.parents.database[gen_check,2]]][[gen_check2]][[21]]), deparse.level = 0)
      }
    }
    list_of_copy <- matrix(unlist(list_of_copy), ncol=6, byrow=TRUE)
    list_of_copy <- list_of_copy[!(list_of_copy[,1]==list_of_copy[,4] & list_of_copy[,2]==list_of_copy[,5] & list_of_copy[,3]==list_of_copy[,6]),]

  } else if(offspring.bve){
    if(verbose) cat("No potential offspring for phenotype import given. Consider all potential individuals.\n")
    first_gen <- min(offspring.bve.parents.database[,1])
    list_of_copy <- list()
    for(gen_check in 1:nrow(offspring.bve.parents.database)){
      for(gen_check2 in offspring.bve.parents.database[gen_check,3]:offspring.bve.parents.database[gen_check,4]){
        first_gen <- min(first_gen, population$breeding[[offspring.bve.parents.database[gen_check,1]]][[offspring.bve.parents.database[gen_check,2]]][[gen_check2]][[21]][,1])
        list_of_copy[[length(list_of_copy)+1]] <- rbind(offspring.bve.parents.database[gen_check,1], offspring.bve.parents.database[gen_check,2], gen_check2 ,
                                                         t(population$breeding[[offspring.bve.parents.database[gen_check,1]]][[offspring.bve.parents.database[gen_check,2]]][[gen_check2]][[21]]), deparse.level = 0)
      }
    }
    list_of_copy <- matrix(unlist(list_of_copy), ncol=6, byrow=TRUE)
    list_of_copy <- list_of_copy[!(list_of_copy[,1]==list_of_copy[,4] & list_of_copy[,2]==list_of_copy[,5] & list_of_copy[,3]==list_of_copy[,6]),]
    ped_off <- get.pedigree(population, gen = first_gen:nrow(population$info$size), raw=TRUE)
    candidates <- rep(TRUE,nrow(ped_off))

    # Remove Individuals themselves
    for(index in 1:nrow(offspring.bve.parents.database)){
      act <- offspring.bve.parents.database[index,]
      candidates[ped_off[,1]==act[1] & ped_off[,2] == act[2] & ped_off[,3]>= act[3] & ped_off[,6]<=act[4]] <- FALSE
    }

    # Remove founder
    candidates[ped_off[,1]==ped_off[,4] & ped_off[,1]==ped_off[,7] &
                 ped_off[,2]==ped_off[,5] & ped_off[,2]==ped_off[,8] &
                 ped_off[,3]==ped_off[,6] & ped_off[,3]==ped_off[,9]] <- FALSE

    ped_off <- ped_off[candidates,]

    #candidates <- rep(FALSE,nrow(ped_off))
    #for(index in 1:nrow(offspring.bve.parents.database)){
    #  act <- offspring.bve.parents.database[index,]
    #  candidates[ped_off[,4]==act[1] & ped_off[,5] == act[2] & ped_off[,6]>= act[3] & ped_off[,6]<=act[4]] <- TRUE
    #  candidates[ped_off[,7]==act[1] & ped_off[,8] == act[2] & ped_off[,9]>= act[3] & ped_off[,9]<=act[4]] <- TRUE
    #}
    candidates <- TRUE
    male_candidates <- (ped_off[candidates & ped_off[,2]==1 ,])
    female_candidates <- (ped_off[candidates & ped_off[,2]==2 ,])

    male_gen <- unique(male_candidates[,1])
    female_gen  <- unique(female_candidates[,1])

    offspring.bve.offspring.database <- NULL
    if(length(male_gen)>0){
      for(index in male_gen){
        offspring.bve.offspring.database <- rbind(offspring.bve.offspring.database,
                                                  c(index,1,min(male_candidates[male_candidates[,1]==index,3]),
                                                    max(male_candidates[male_candidates[,1]==index,3])))
      }
    }
    if(length(female_gen)>0){
      for(index in female_gen){
        offspring.bve.offspring.database <- rbind(offspring.bve.offspring.database,
                                                  c(index,2,min(female_candidates[female_candidates[,1]==index,3]),
                                                    max(female_candidates[female_candidates[,1]==index,3])))
      }
    }

  }

  if(offspring.bve){
    offspring_id <- get.id(population, offspring.bve.offspring.database)
    is_double <- duplicated(offspring_id[length(offspring_id):1])[length(offspring_id):1]
    if(verbose) cat("Import phenotypes of offspring.\n")
    for(index in 1:nrow(offspring.bve.parents.database)){
      activ.parents <- offspring.bve.parents.database[index,]
      new.bv <- counter <- matrix(0, nrow=population$info$bv.nr, ncol=activ.parents[4]-activ.parents[3]+1)

      next_indi <- 1
      for(index2 in 1:nrow(offspring.bve.offspring.database)){
        activ.offspring <- offspring.bve.offspring.database[index2,]
        n.animals <- activ.offspring[4]-activ.offspring[3] +1
        if(n.animals>0){
          for(index3 in (activ.offspring[3]:activ.offspring[4])[!is_double[next_indi:(next_indi+n.animals-1)]]){
            parent1 <- population$breeding[[activ.offspring[1]]][[activ.offspring[2]]][[index3]][[7]]
            parent2 <- population$breeding[[activ.offspring[1]]][[activ.offspring[2]]][[index3]][[8]]
            if(length(list_of_copy)>0){
              replace1 <- which(parent1[1]==list_of_copy[,4] & parent1[2]==list_of_copy[,5] & parent1[3]==list_of_copy[,6])
              replace2 <- which(parent2[1]==list_of_copy[,4] & parent2[2]==list_of_copy[,5] & parent2[3]==list_of_copy[,6])
              if(length(replace1)>=1){
                parent1 <- list_of_copy[replace1[1],1:3]
              }
              if(length(replace2)>=1){
                parent2 <- list_of_copy[replace2[1],1:3]
              }
            }
            if((nrow(population$breeding[[activ.offspring[1]]][[activ.offspring[2]]][[index3]][[21]]) ==
               nrow(population$breeding[[parent1[1]]][[parent1[2]]][[parent1[3]]][[21]])) &&
               (prod(population$breeding[[activ.offspring[1]]][[activ.offspring[2]]][[index3]][[21]][1,] ==
               population$breeding[[parent1[1]]][[parent1[2]]][[parent1[3]]][[21]][1,])==1)){
              parent1 <- c(-1,-1,-1)
            }
            if((nrow(population$breeding[[activ.offspring[1]]][[activ.offspring[2]]][[index3]][[21]]) ==
                nrow(population$breeding[[parent2[1]]][[parent2[2]]][[parent2[3]]][[21]])) &&
               (prod(population$breeding[[activ.offspring[1]]][[activ.offspring[2]]][[index3]][[21]][1,] ==
                     population$breeding[[parent2[1]]][[parent2[2]]][[parent2[3]]][[21]][1,])==1)){
              parent2 <- c(-1,-1,-1)
            }
            if(parent1[1]==activ.parents[1] && parent1[2]==activ.parents[2] && parent1[3]>= activ.parents[3] && parent1[3]<= activ.parents[4]){
              activ.take <- population$breeding[[activ.offspring[1]]][[activ.offspring[2]]][[activ.offspring[3]]][[15]]>0 & !is.na(population$breeding[[activ.offspring[1]]][[activ.offspring[2]+8]][,index3])
              new.bv[activ.take,parent1[3] - activ.parents[3]+1] <- (new.bv[,parent1[3]- activ.parents[3]+1] + population$breeding[[activ.offspring[1]]][[activ.offspring[2]+8]][,index3] * population$breeding[[activ.offspring[1]]][[activ.offspring[2]]][[activ.offspring[3]]][[15]])[activ.take]
              counter[activ.take,parent1[3]- activ.parents[3]+1] <- (counter[,parent1[3]- activ.parents[3]+1] + population$breeding[[activ.offspring[1]]][[activ.offspring[2]]][[activ.offspring[3]]][[15]])[activ.take]
            }
            if(parent2[1]==activ.parents[1] && parent2[2]==activ.parents[2] && parent2[3]>= activ.parents[3] && parent2[3]<= activ.parents[4]){
              activ.take <-  population$breeding[[activ.offspring[1]]][[activ.offspring[2]]][[activ.offspring[3]]][[15]]>0 & !is.na(population$breeding[[activ.offspring[1]]][[activ.offspring[2]+8]][,index3])
              new.bv[activ.take,parent2[3]- activ.parents[3]+1] <- (new.bv[,parent2[3]- activ.parents[3]+1] + population$breeding[[activ.offspring[1]]][[activ.offspring[2]+8]][,index3] * population$breeding[[activ.offspring[1]]][[activ.offspring[2]]][[activ.offspring[3]]][[15]])[activ.take]
              counter[activ.take,parent2[3]- activ.parents[3]+1] <- (counter[,parent2[3]- activ.parents[3]+1] + population$breeding[[activ.offspring[1]]][[activ.offspring[2]]][[activ.offspring[3]]][[15]])[activ.take]
            }
            if(sum(is.na(new.bv))>0){
              stop()
            }
          }
        }
        next_indi <- next_indi + n.animals


      }
      population$breeding[[activ.parents[1]]][[activ.parents[2]+26]][,activ.parents[3]:activ.parents[4]] <- new.bv / counter
      population$breeding[[activ.parents[1]]][[activ.parents[2]+28]][,activ.parents[3]:activ.parents[4]] <- counter
      if(sum(counter==0)>0){
        if(verbose) cat(paste0(sum(counter==0), " phenotype entries without valid offspring for phenotype import from offspring! Set offspring phenotype to NA"))
        population$breeding[[activ.parents[1]]][[activ.parents[2]+26]][,activ.parents[3]:activ.parents[4]][counter==0] <- NA
      }

    }

    if(sum(counter)==0){
      if(input.phenotype!="own"){
        input.phenotype <- "own"
        if(verbose) cat("No phenotypes to import. Automatically set input.phenotype to own")
      }
    }
  }


  ## Culling Module
  if(culling){
    culling.database <- get.database(population, gen=culling.gen, database=culling.database, cohorts=culling.cohort)

    if(length(culling.index)==1 && culling.index=="lastindex"){
      culling.index <- get.selectionindex(population, cohorts=culling.database)
    }
    if(population$info$bv.nr>0){
      culling.bv <- colSums(get.bv(population, database = culling.database) * culling.index)
    } else{
      culling.bv <- numeric(length(get.id(population, database= culling.database)))
    }
    n_animals <- length(culling.bv)

    if(sum(abs(culling.index))==0){
      culling.prob <- rep(culling.share1, n_animals)
    } else{
      culling.prob <- (culling.bv - culling.bv1) * (culling.share1 - culling.share2)/ (culling.bv1 - culling.bv2) + culling.share1
      culling.prob[culling.prob>1] <- 1
      culling.prob[culling.prob<0] <- 0
    }
    if(length(culling.prob)>length(culling.single)){
      culling.single <- rep(culling.single, length.out=length(culling.prob))
    }
    culling.prob <- culling.prob * culling.single
    culling.action <- stats::rbinom(n_animals,1,culling.prob)==1


    store <- population$breeding[[culling.database[1]]][[culling.database[2]+4]][culling.database[3]:culling.database[4]]
    population$breeding[[culling.database[1]]][[culling.database[2]+4]][culling.database[3]:culling.database[4]][culling.action] <- (-1)
    population$breeding[[culling.database[1]]][[culling.database[2]+24]][culling.database[3]:culling.database[4]][culling.action] <-
      population$breeding[[culling.database[1]]][[culling.database[2]+22]][culling.database[3]:culling.database[4]][culling.action] + culling.time
    new_death <- population$breeding[[culling.database[1]]][[culling.database[2]+4]][culling.database[3]:culling.database[4]] != store

    if(culling.all.copy){
      for(kindex in (culling.database[3]:culling.database[4])[new_death]){
        active_indi <- c(culling.database[1:2], kindex)
        if(nrow(population$breeding[[active_indi[1]]][[active_indi[2]]][[kindex]][[21]])){
          for(switch in 1:nrow(population$breeding[[active_indi[1]]][[active_indi[2]]][[kindex]][[21]])){
            active_copy <- population$breeding[[active_indi[1]]][[active_indi[2]]][[kindex]][[21]][switch,]

            population$breeding[[active_copy[1]]][[active_copy[2]+4]][active_copy[3]] <-
              population$breeding[[active_indi[1]]][[active_indi[2]+4]][active_indi[3]]
            population$breeding[[active_copy[1]]][[active_copy[2]+24]][active_copy[3]] <-
              population$breeding[[active_indi[1]]][[active_indi[2]+24]][active_indi[3]]

          }
        }
      }
    }
    n_death <- sum(new_death)

    if(n_death>0){
      population$breeding[[culling.database[1]]][[culling.database[2]+16]][culling.database[3]:culling.database[4]][which(new_death)] <- population$breeding[[culling.database[1]]][[culling.database[2]+10]][culling.database[3]:culling.database[4]][which(new_death)]
      active_cohort <- which(population$info$cohorts[,1]==culling.cohort)
      if(length(active_cohort)>0){
        if(length(population$info$culling.stats)<active_cohort){
          population$info$culling.stats[[active_cohort]] <- cbind(culling.name, n_death, deparse.level = 0)
        } else{
          population$info$culling.stats[[active_cohort]] <- rbind(population$info$culling.stats[[active_cohort]],
                                                                  c(culling.name, n_death), deparse.level = 0)

        }
      }


    }
  }
}

  #######################################################################
  ################ Breeding value estimation ############################
  #######################################################################
{
  ## Track computing times
  if(store.comp.times){
    comp.times[4] <- as.numeric(Sys.time())
  }

  if(store.comp.times.bve){
    comp.times.bve <- numeric(5)
    z_chol <- 0 # Cholesky decomposition
    z_uhat <- 0 # rrblup
    zcalc <- 0 # calculation of Z
    z_ped <- 0 # calculation of pedigree relationship matrix
    z_h <- 0 # calculation single-step H matrix
    comp.times.bve[1] <- as.numeric(Sys.time())
  }


  u_hat <- NULL
  gwas_hat <- NULL
  u_hat_possible <- FALSE

  ## BVE according to average performance of parents / grandparents
  if(bve.parent.mean || bve.grandparent.mean){
    addpheno <- FALSE
    if(bve.mean.between=="pheno"){
      addmean <- 8
    } else if(bve.mean.between=="bv"){
      addmean <- 6
    } else if(bve.mean.between=="bve"){
      addmean <- 2
    } else{
      addmean <- 2
      addpheno <- TRUE
    }
    for(index in 1:nrow(bve.insert.database)){
      activ.base <- bve.insert.database[index,]
      if(diff(activ.base[3:4])>=0){
        for(index2 in activ.base[3]:activ.base[4]){
          reference <- rbind(population$breeding[[activ.base[1]]][[activ.base[2]]][[index2]][[7]],population$breeding[[activ.base[1]]][[activ.base[2]]][[index2]][[8]])
          if(bve.grandparent.mean){
            reference <- rbind(population$breeding[[reference[1,1]]][[reference[1,2]]][[reference[1,3]]][[7]],
                              population$breeding[[reference[1,1]]][[reference[1,2]]][[reference[1,3]]][[8]],
                              population$breeding[[reference[2,1]]][[reference[2,2]]][[reference[2,3]]][[7]],
                              population$breeding[[reference[2,1]]][[reference[2,2]]][[reference[2,3]]][[8]])
          }
          bve_import <- matrix(0, nrow=nrow(reference), ncol=population$info$bv.nr)
          for(index3 in 1:nrow(reference)){
            bve_import[index3,] <- population$breeding[[reference[index3,1]]][[reference[index3,2]+addmean]][,reference[index3,3]]
          }
          if(addpheno){
            for(index3 in 1:nrow(reference)){
              bve_import[index3,bve_import[index3,]==0] <- population$breeding[[reference[index3,1]]][[reference[index3,2]+8]][bve_import[index3,]==0,reference[index3,3]]
            }
          }
          bve_import[bve_import==0] <- NA
          population$breeding[[activ.base[1]]][[activ.base[2]+2]][,index2] <- colMeans(bve_import, na.rm=TRUE)
          if(sum(is.na(population$breeding[[activ.base[1]]][[activ.base[2]+2]][,index2]))>0){
            population$breeding[[activ.base[1]]][[activ.base[2]+2]][is.na(population$breeding[[activ.base[1]]][[activ.base[2]+2]][,index2]),index2] <- 0
          }

        }
      }
      if(report.accuracy){
        y_real_report <- NULL
        y_hat_report <- NULL
        for(index in 1:nrow(bve.insert.database)){
          activ.base <- bve.insert.database[index,]
          if(diff(activ.base[3:4])>=0){
            y_real_report <- cbind(y_real_report, population$breeding[[activ.base[1]]][[activ.base[2]+6]][,activ.base[3]:activ.base[4], drop=FALSE], deparse.level = 0)
            y_hat_report <- cbind(y_hat_report, population$breeding[[activ.base[1]]][[activ.base[2]+2]][,activ.base[3]:activ.base[4], drop=FALSE], deparse.level = 0)
          }
        }

        y_hat_report[y_hat_report==0] <- NA
        if(verbose) cat("Correlation between genetic values and BVE (parent-mean BVE):\n")
        acc <- suppressWarnings(stats::cor(t(y_real_report), t(y_hat_report), use="pairwise.complete.obs"))


        if(sum(is.na(acc))>0){
          acc[is.na(acc)] <- 0
        }
        if(verbose) cat(diag(acc))
        if(verbose) cat("\n")
      }

    }
  } else if(phenotype.bv || bve && sum(sigma.e)==0 && BGLR.bve==FALSE ){
    ## Use Phenotypes as breeding values
    ## You can also use phenotypes in the selection procedure - so potentially removable at some point.
    if(store.bve.data){
      n.animals <- sum(bve.database[,4]-bve.database[,3]+1)
      y <- y_real <- y_hat <- array(0, dim=c(n.animals,population$info$bv.nr)) # schaetzung sigma.g
      cindex <- 1
      for(index in 1:nrow(bve.database)){
        k.database <- bve.database[index,]
        temp1 <- 1:population$info$bv.nr
        if(diff(k.database[3:4])>=0){
          for(kindex in k.database[3]:k.database[4]){
            # t() not needed just to be safe when using multiple individuals at once later
            y[cindex,temp1] <- y_hat[cindex,temp1] <- t(population$breeding[[k.database[[1]]]][[8+k.database[[2]]]][temp1,kindex])
            y_real[cindex,temp1] <- t(population$breeding[[k.database[[1]]]][[6+k.database[[2]]]][temp1,kindex])
            cindex <- cindex +1
          }
        }

      }
    }
    for(index in 1:nrow(bve.insert.database)){
      activ.base <- bve.insert.database[index,]
      if(diff(activ.base[3:4])>=0){
        population$breeding[[activ.base[1]]][[activ.base[2]+2]][,activ.base[3]:activ.base[4]] <- population$breeding[[activ.base[1]]][[activ.base[2]+8]][,activ.base[3]:activ.base[4]]
        population$breeding[[activ.base[1]]][[activ.base[2]+2]][,activ.base[3]:activ.base[4]][is.na(population$breeding[[activ.base[1]]][[activ.base[2]+2]][,activ.base[3]:activ.base[4]])] <- 0 # Breeding values cant be missing
      }
      if(report.accuracy){
        y_real_report <- NULL
        y_hat_report <- NULL
        for(index in 1:nrow(bve.insert.database)){
          activ.base <- bve.insert.database[index,]
          if(diff(activ.base[3:4])>=0){
            y_real_report <- cbind(y_real_report, population$breeding[[activ.base[1]]][[activ.base[2]+6]][,activ.base[3]:activ.base[4], drop=FALSE])
            y_hat_report <- cbind(y_hat_report, population$breeding[[activ.base[1]]][[activ.base[2]+8]][,activ.base[3]:activ.base[4], drop=FALSE])
          }
        }

        y_hat_report[y_hat_report==0] <- NA
        if(verbose) cat("Correlation between genetic values and BVE (phenotypic BVE):\n")
        acc <- suppressWarnings(stats::cor(t(y_real_report), t(y_hat_report), use="pairwise.complete.obs"))


        if(sum(is.na(acc))>0){
          acc[is.na(acc)] <- 0
        }
        if(verbose) cat(diag(acc))
        if(verbose) cat("\n")
      }

    }
  } else if(bve.pseudo){

    ## Simulate a breeding value estimation
    ## This lacks structure (e.g. related individuals are usally all overestimated or all underestimated)
    ## but its super fast!

    n.animals <- sum(bve.insert.database[,4]-bve.insert.database[,3]+1)
    y_real <- y_hat <- array(0, dim=c(n.animals,population$info$bv.nr)) # schaetzung sigma.g
    cindex <- 1
    for(index in 1:nrow(bve.insert.database)){
      k.database <- bve.insert.database[index,]
      if(diff(k.database[3:4])>=0){
        # t() not needed just to be safe when using multiple individuals at once later
        y_real[cindex:(cindex+k.database[4]-k.database[3]),] <- t(population$breeding[[k.database[[1]]]][[6+k.database[[2]]]][,k.database[3]:k.database[4]])
        cindex <- cindex + k.database[4] - k.database[3] +1
      }
    }
    for(index in 1:population$info$bv.nr){
      var_trait <- stats::var(y_real[,index])
      if(bve.pseudo.accuracy[index]==0){
        var_residual <- (var_trait - (bve.pseudo.accuracy[index])^2 * var_trait) / (bve.pseudo.accuracy[index])^2
        y_real[,index] <- y_hat[,index] <- 0
      } else{
        var_residual <- (var_trait - (bve.pseudo.accuracy[index])^2 * var_trait) / (bve.pseudo.accuracy[index])^2
        y_hat[,index] <- (y_real[,index] + stats::rnorm(n.animals, sd= sqrt(var_residual)))

        if(bve.pseudo.accuracy[index]<0){
          temp1 <-  -y_hat[,index]
          temp1 <- temp1 + max(y_hat[,index]) -max(temp1)
          y_hat[,index] <- temp1
        }
      }
      if(!is.matrix(y_hat)){
        y_hat <- matrix(y_hat, ncol=1)
        y_real <- matrix(y_real, ncol=1)
      }

    }

    # Enter BVE
    cindex <- 1
    for(index in 1:nrow(bve.insert.database)){
      k.database <- bve.insert.database[index,]
      if(diff(k.database[3:4])>=0){
        # t() not needed just to be safe when using multiple individuals at once later
        population$breeding[[k.database[[1]]]][[2+k.database[[2]]]][,k.database[3]:k.database[4]] <- t(y_hat[cindex:(cindex+k.database[4]-k.database[3]),])
        cindex <- cindex + k.database[4] - k.database[3] +1
      }
    }

    if(report.accuracy){
      y_hat_report <- y_hat
      y_real_report <- y_real
      y_hat_report[y_hat_report==0] <- NA
      if(verbose) cat("Correlation between genetic values and simulated BVE:\n")
      acc <- suppressWarnings(stats::cor(y_real_report, y_hat_report, use="pairwise.complete.obs"))


      if(sum(is.na(acc))>0){
        acc[is.na(acc)] <- 0
      }
      if(verbose) cat(diag(acc))
      if(verbose) cat("\n")
    }

  } else if(bve){
    ## "proper" breeding value estimation

    if(verbose && (relationship.matrix!="kinship" && relationship.matrix !="pedigree")) cat("Start genomic BVE.\n")
    if(verbose && (relationship.matrix=="kinship" || relationship.matrix =="pedigree")) cat("Start pedigree BVE.\n")

    ## Derive if individuals are using multiple times in the BVE
    ## Import Phenotypes from all copies of an individual (potentially use individual entry [[21]] for efficiency?)
    loop_elements_list <- derive.loop.elements(population=population, bve.database=bve.database,
                                          bve.class=bve.class, bve.avoid.duplicates=bve.avoid.duplicates,
                                          store.which.adding = TRUE, list.of.copys = TRUE)
    loop_elements <- loop_elements_list[[1]]
    loop_elements_copy <- loop_elements_list[[3]]

    n.animals <- nrow(loop_elements)
    genotyped <- numeric(n.animals)

    stay.loop.elements <- NULL
    n.rep <- 0
    if(length(loop_elements_copy)>0){
      n.rep <- nrow(loop_elements_copy)
    }
    # remove non-genotyped samples in case no pedigree-based estimation // single-step
    if(remove.non.genotyped && singlestep.active==FALSE && (relationship.matrix!="kinship" && relationship.matrix !="pedigree")){
      for(index in 1:n.animals){

        k.database <- bve.database[loop_elements[index,3],]
        kindex <- loop_elements[index,2]
        genotyped[index] <- population$breeding[[k.database[[1]]]][[k.database[[2]]]][[kindex]][[16]]
      }

      if(n.rep>0){
        genotyped.copy <- numeric(n.rep)
        for(index in 1:n.rep){
          kindex <- loop_elements_copy[index,2]
          k.database <- bve.database[loop_elements_copy[index,3],]
          if(population$breeding[[k.database[[1]]]][[k.database[[2]]]][[kindex]][[16]]==1){
            genotyped[loop_elements_copy[index,6]] <-  genotyped.copy[index] <- population$breeding[[k.database[[1]]]][[k.database[[2]]]][[kindex]][[16]]

          }
        }
        remove.loop.elements <- which(genotyped.copy==0)

        if(length(remove.loop.elements)>0 && FALSE){
          loop_elements_list[[3]] <- loop_elements_list[[3]][-remove.loop.elements,]
        }
      }
      stay.loop.elements <- which(genotyped==1)
      remove.loop.elements <- which(genotyped==0)
      if(length(remove.loop.elements)>0){
        loop_elements_list[[1]] <- loop_elements_list[[1]][-remove.loop.elements,]
        loop_elements_list[[2]] <- loop_elements_list[[2]][-remove.loop.elements]
      }



      loop_elements <- loop_elements_list[[1]]
      n.animals <- nrow(loop_elements)
      loop_elements_list[[1]][,1] <- 1:n.animals
      genotyped <- numeric(n.animals)

    }


    # sequenceZ is to not process all SNPs at the same time. Especially relevant for large scale datasets!
    if((relationship.matrix!="kinship" && relationship.matrix !="pedigree")){
      if(sequenceZ){
        if(maxZtotal>0){
          maxZ <- floor(maxZtotal / n.animals)
        }
        Zt <- array(0L,dim=c(maxZ,n.animals))
      } else{
        Zt <- array(0L,dim=c(sum(population$info$snp), n.animals))
      }
    }
    y <- y_real <- y_real2 <- y_hat <- y_reli <- y_parent <- array(0,dim=c(n.animals,population$info$bv.nr))
    X <- matrix(1, nrow=n.animals,ncol=1)

    grid.position <- numeric(n.animals)
    cindex <- 1
    size <- cumsum(c(0,as.vector(t(population$info$size))))
    y_obs <- matrix(0, nrow=n.animals, ncol=population$info$bv.nr)
    genotyped_array <- matrix(FALSE, ncol=length(population$info$array.name), nrow=n.animals)


    # Import phenotypes / genomic values
    for(index in 1:n.animals){

      k.database <- bve.database[loop_elements[index,3],]
      kindex <- loop_elements[index,2]
      y[index,] <- population$breeding[[k.database[[1]]]][[8+k.database[[2]]]][,kindex]
      y_real[index,] <- y_real2[index,] <-population$breeding[[k.database[[1]]]][[6+k.database[[2]]]][,kindex]

      if(population$breeding[[k.database[1]]][[k.database[[2]]]][[kindex]][[25]]==FALSE){
        if(length(population$breeding[[k.database[1]]][[k.database[[2]]]][[kindex]][[26]])==0){
          y_real2[index,] <- NA
        } else{
          y_real2[index,-population$breeding[[k.database[1]]][[k.database[[2]]]][[kindex]][[26]]] <- NA
        }

      }
      y_obs[index,] <- population$breeding[[k.database[[1]]]][[k.database[[2]]]][[kindex]][[15]]
      grid.position[index] <- kindex + size[sum(k.database[1:2]*c(2,1))-2] # how many individuals are in earlier generations
      genotyped[index] <- population$breeding[[k.database[[1]]]][[k.database[[2]]]][[kindex]][[16]]
      genotyped_array[index, population$breeding[[k.database[[1]]]][[k.database[[2]]]][[kindex]][[22]]] <- TRUE
      if(!remove.non.genotyped && singlestep.active==FALSE){
        genotyped[index] <- 1
        genotyped_array[index, ] <- TRUE
      }

      if(estimate.add.gen.var){
        father <- population$breeding[[k.database[1]]][[k.database[2]]][[kindex]][[7]]
        mother <- population$breeding[[k.database[1]]][[k.database[2]]][[kindex]][[8]]
        if(sum(father[1:3]==c(k.database[1:2], kindex))==3 ||sum(mother[1:3]==c(k.database[1:2], kindex))==3){
          if(verbose) cat("Schaetzung der additiv genetischen Varianz extrem problematisch. Kein Elterntier fuer jedes Tier vorhanden!\n")
        }
        y_parent[index,] <- mean(population$breeding[[father[1]]][[8+father[2]]][,father[3]],population$breeding[[mother[1]]][[8+mother[2]]][,mother[3]])
      }
    }
    if(n.rep>0){
      for(index in 1:n.rep){
        kindex <- loop_elements_copy[index,2]
        k.database <- bve.database[loop_elements_copy[index,3],]
        if(length(stay.loop.elements)>0){
          non_copy <- which(stay.loop.elements==loop_elements_copy[index,6])
        } else{
          non_copy <- loop_elements_copy[index,6]
        }

        if(length(non_copy)==1 && sum(population$breeding[[k.database[[1]]]][[k.database[[2]]]][[kindex]][[15]]> y_obs[non_copy,])>=1){
          switches <- population$breeding[[k.database[[1]]]][[k.database[[2]]]][[kindex]][[15]]> y_obs[non_copy,]
          y[non_copy,switches] <- population$breeding[[k.database[[1]]]][[8+k.database[[2]]]][switches,kindex]
          y_obs[non_copy,switches] <- population$breeding[[k.database[[1]]]][[k.database[[2]]]][[kindex]][[15]][switches]
        }
      }
    }

    # Check if some copy is genotyped
    if(nrow(loop_elements_list[[3]])>0){
      for(index in 1:nrow(loop_elements_list[[3]])){
        kindex <- loop_elements_list[[3]][index,2]
        k.database <- bve.database[loop_elements_list[[3]][index,3],]
        if(length(stay.loop.elements)>0){
          non_copy <- which(stay.loop.elements==loop_elements_copy[index,6])
        } else{
          non_copy <- loop_elements_copy[index,6]
        }
        genotyped[non_copy] <- max(genotyped[non_copy],population$breeding[[k.database[[1]]]][[k.database[[2]]]][[kindex]][[16]])
        genotyped_array[non_copy, population$breeding[[k.database[[1]]]][[k.database[[2]]]][[kindex]][[22]]] <- TRUE
      }
    }
    genotype.included <- which(genotyped==1)

    batch <- ceiling(nrow(loop_elements)/ncore)
    batche <- list()
    for(index in 1:ncore){
      batche[[index]] <- (batch*(index-1)+1):min(batch*index, nrow(loop_elements))
    }

    # in case of copied individuals use the copy that has been phenotyped the most
    if(input.phenotype!="own"){
      for(index2 in 1:nrow(offspring.bve.parents.database)){
        activ.database <- offspring.bve.parents.database[index2,]
        for(index in 1:nrow(loop_elements)){
          kindex <- loop_elements[index,2]
          k.database <- bve.database[loop_elements[index,3],]
          activ.indi <- population$breeding[[k.database[[1]]]][[k.database[[2]]]][[kindex]][[21]]
          import <- activ.indi[,1]==activ.database[1] & activ.indi[,2]==activ.database[2] & activ.indi[,3]>=activ.database[3] & activ.indi[,3] <= activ.database[4]
          if(sum(import)>0){
            own_pheno <- y[index,]
            n_obs <- y_obs[index,]
            off_pheno <- population$breeding[[activ.database[1]]][[activ.database[2]+26]][,kindex]
            n_off <- population$breeding[[activ.database[1]]][[activ.database[2]+28]][,kindex]
            if(input.phenotype=="off"){
              y[index,] <- off_pheno
              y_obs[index,] <- n_off
            } else if(input.phenotype=="mean"){
              take_both <- !is.na(own_pheno) & !is.na(off_pheno)
              take_own <- !is.na(own_pheno) & is.na(off_pheno)
              take_off <- is.na(own_pheno) & !is.na(off_pheno)
              y[index,take_both] <- ((own_pheno + off_pheno)/((own_pheno!=0) + (off_pheno!=0)))[take_both]
              y[index,take_own] <- (own_pheno /own_pheno!=0)[take_own]
              y[index,take_off] <- (off_pheno/off_pheno!=0)[take_off]
              y_obs[index,] <- n_obs + n_off
            } else if(input.phenotype=="weighted"){
              take_both <- !is.na(own_pheno) & !is.na(off_pheno)
              take_own <- !is.na(own_pheno) & is.na(off_pheno)
              take_off <- is.na(own_pheno) & !is.na(off_pheno)
              y[index,take_both] <- ((own_pheno*n_obs*2 + off_pheno*n_off)/(n_obs*2 + n_off))[take_both]
              y[index,take_own] <- (own_pheno /own_pheno!=0)[take_own]
              y[index,take_off] <- (off_pheno/off_pheno!=0)[take_off]
              y_obs[index,] <- n_obs + n_off / 2
            }
          }
        }
      }

    }
    # Import Z
    if(sequenceZ==FALSE && (relationship.matrix != "pedigree" && relationship.matrix != "kinship")){
      if(miraculix){

        if(store.comp.times.bve){
          before <- as.numeric(Sys.time())
        }

        if (requireNamespace("miraculix", quietly = TRUE)) {
          Z.code <- miraculix::computeSNPS(population, loop_elements[,4], loop_elements[,5], loop_elements[,2], what="geno", output_compressed=TRUE)
          if(relationship.matrix!="vanRaden"){
            Zt <- miraculix::computeSNPS(population, loop_elements[,4], loop_elements[,5], loop_elements[,2], what="geno", output_compressed=FALSE)
          }
        }
        if(store.comp.times.bve){
          after <- as.numeric(Sys.time())
          zcalc <- zcalc + after - before
        }

      } else if(ncore>1){
        if(store.comp.times.bve){
          before <- as.numeric(Sys.time())
        }

        if(backend=="doParallel"){
          if (requireNamespace("doParallel", quietly = TRUE)) {
            doParallel::registerDoParallel(cores=ncore)
          } else{
            stop("Usage of doParallel without being installed")
          }
        } else if(backend=="doMPI"){
          if (requireNamespace("doMPI", quietly = TRUE)) {
            cl <- doMPI::startMPIcluster(count=ncore)
            doMPI::registerDoMPI(cl)
          } else{
            stop("Usage of doMPI without being installed!")
          }
        } else{
          if(verbose) cat("No valid backend specified.\n")
        }

        if (requireNamespace("foreach", quietly = TRUE)) {
        } else{
          stop("Usage of foreach without being installed!")
        }
        Zt <- foreach::foreach(indexb=1:ncore, .combine = "cbind", .multicombine = TRUE,.maxcombine = 1000,
                     .packages="MoBPS") %dopar% {
          Ztpar <- array(0,dim=c(sum(population$info$snp), length(batche[[indexb]])))
          col_index <- 1
          for(index in batche[[indexb]]){
            k.database <- bve.database[loop_elements[index,3],]
            cindex <- col_index
            kindex <- loop_elements[index,2]
            Ztpar[,cindex] <- base::as.integer(colSums(compute.snps(population, k.database[1],k.database[2],kindex, import.position.calculation=import.position.calculation,decodeOriginsU=decodeOriginsU, bit.storing=bit.storing, nbits=nbits, output_compressed=FALSE)))
            col_index <- col_index +1
          }
          if(Z.integer){
            storage.mode(Zpar) <- "integer"
          }
          Ztpar
        }
        if(backend=="doParallel"){
          doParallel::stopImplicitCluster()
        } else if(backend=="doMPI"){
          if (requireNamespace("doMPI", quietly = TRUE)) {
            doMPI::closeCluster(cl)
          } else{
            stop("Usage of doMPI without being installed!")
          }


        }
        if(store.comp.times.bve){
          after <- as.numeric(Sys.time())
          zcalc <- zcalc + after - before
        }
      } else{
        if(store.comp.times.bve){
          before <- as.numeric(Sys.time())
        }

        for(index in 1:n.animals){
          k.database <- bve.database[loop_elements[index,3],]
          cindex <- index
          kindex <- loop_elements[index,2]
          Zt[,cindex] <- base::as.integer(colSums(compute.snps(population, k.database[1],k.database[2],kindex, import.position.calculation=import.position.calculation,decodeOriginsU=decodeOriginsU, bit.storing=bit.storing, nbits=nbits, output_compressed=FALSE)))
        }
        if(store.comp.times.bve){
          after <- as.numeric(Sys.time())
          zcalc <- zcalc + after - before
        }
      }

    }

    to_remove <- NULL
    genotyped_array_variants <- unique(genotyped_array)
    used_arrays <- colSums(genotyped_array)
    if(sum(used_arrays == nrow(genotyped_array)  & !population$info$array.is_subset)>0){
      to_remove <- NULL
    } else{

      if(bve.imputation){
        keeps <- rep(FALSE, sum(population$info$snp))
      } else{
        keeps <- rep(TRUE, sum(population$info$snp))
      }

      for(index in 1:nrow(genotyped_array_variants)){
        activ_keep <- rep(FALSE, sum(population$info$snp))
        for(index2 in which(genotyped_array_variants[index,])){
          activ_keep <- activ_keep | population$info$array.markers[[index2]]
        }

        if(bve.imputation){
          keeps <- keeps | activ_keep
        } else{
          keeps <- keeps & activ_keep
        }

      }
      to_remove <- c(to_remove, which(!keeps))

    }
    if(remove.effect.position && sequenceZ==FALSE){
      to_remove <- c(to_remove, population$info$effect.p)
    }

    if(length(to_remove)==sum(population$info$snp)){

      if(relationship.matrix != "pedigree" && relationship.matrix != "kinship"){
        warning("No genotyped individuals!")
      }
      to_remove <- NULL
    }

    if(length(to_remove)>0 && (relationship.matrix != "pedigree" && relationship.matrix != "kinship")){


      if(miraculix && exists("Z.code")){
        if (requireNamespace("miraculix", quietly = TRUE)) {
          Z.code <- as.matrix(Z.code)
          Z.code <- miraculix::genomicmatrix(Z.code[-to_remove,])
          #Z.code <- miraculix::zeroNthGeno(Z.code, to_remove) ### currently not implemented in miraculix
        }
      } else{
        Zt <- Zt[-to_remove,]
      }
      if(verbose) cat(paste0(sum(population$info$snp) - length(to_remove), " markers survived filtering for BVE.\n"))


    }

    if(bve.imputation.errorrate>0 && (relationship.matrix != "pedigree" && relationship.matrix != "kinship")){

      if(length(to_remove)>0){
        is_imputed <- get.genotyped.snp(population, database = cbind(loop_elements[,4], loop_elements[,5], loop_elements[,2]))[-to_remove,]==0
      } else{
        is_imputed <- get.genotyped.snp(population, database = cbind(loop_elements[,4], loop_elements[,5], loop_elements[,2]))==0
      }


      if(miraculix && exists("Z.code")){
        if (requireNamespace("miraculix", quietly = TRUE)) {
          temp1 <- as.matrix(Z.code)[is_imputed]
        }
      } else{
        temp1 <- Zt[is_imputed]
      }


      temp2 <- temp1
      zeros <- which(temp1==0)
      ones <- which(temp1==1)
      twos <- which(temp1==2)

      if(length(zeros)>0){
        temp1[zeros][stats::rbinom(length(zeros),1,bve.imputation.errorrate^2)==1] <- 2
        temp1[zeros][stats::rbinom(length(zeros),1,2*bve.imputation.errorrate * (1-bve.imputation.errorrate))==1 & temp1[zeros]==0] <- 1
      }


      if(length(ones)>0){
        temp1[ones][stats::rbinom(length(ones),1,bve.imputation.errorrate*(1-bve.imputation.errorrate))==1] <- 0
        temp1[ones][stats::rbinom(length(ones),1,bve.imputation.errorrate*(1-bve.imputation.errorrate))==1 & temp1[ones]==1] <- 2
      }


      if(length(twos)>0){
        temp1[twos][stats::rbinom(length(twos),1,bve.imputation.errorrate^2)==1] <- 0
        temp1[twos][stats::rbinom(length(twos),1,2*bve.imputation.errorrate * (1-bve.imputation.errorrate))==1 & temp1[twos]==2] <- 1
      }


      if(miraculix && exists("Z.code")){
        if (requireNamespace("miraculix", quietly = TRUE)) {

          Z.code <- as.matrix(Z.code)
          Z.code[is_imputed] <- temp1
          Z.code <- miraculix::genomicmatrix(Z.code)
        }
      } else{
        Zt[is_imputed] <- temp1
      }

      if(verbose) cat(paste0("A total of ", length(temp1), " entries in the genotype matrix were imputed.\n"))
      if(verbose) cat(paste0("Accuracy of the imputation: ", mean(temp1==temp2), ".\n"))

    }


    if(bve.0isNA){
      y[y==0] <- NA
    }

    if(store.comp.times.bve){
      comp.times.bve[2] <- as.numeric(Sys.time())
    }

    # Derive relationship matrix sequenceZ and single-step only for vanRaden based genomic relationship
    if(relationship.matrix=="kinship" || relationship.matrix == "pedigree"){
      z_ped <- z_ped - as.numeric(Sys.time())
      A <- kinship.exp(population, database=bve.database, depth.pedigree=depth.pedigree, elements = loop_elements_list[[2]], mult=2, verbose=verbose)
      z_ped <- z_ped + as.numeric(Sys.time())
    } else if(relationship.matrix=="vanRaden"){
      if(sequenceZ){
        total_n <- sum(population$info$snp)
        p_i <- numeric(total_n)
        first <- 1
        last <- min(maxZ, total_n)
        A <- matrix(0, n.animals, n.animals)
        for(index3 in 1:ceiling(total_n/maxZ)){
          if(miraculix){
            if(store.comp.times.bve){
              before <- as.numeric(Sys.time())
            }

            if (requireNamespace("miraculix", quietly = TRUE)) {
              Z.code <- miraculix::computeSNPS(population, loop_elements[,4], loop_elements[,5], loop_elements[,2], from_p=first,
                                  to_p=last, what="geno", output_compressed=TRUE)
            }
            if(store.comp.times.bve){
              after <- as.numeric(Sys.time())
              zcalc <- zcalc + after - before
            }



          } else if(ncore>1){
            if(store.comp.times.bve){
              before <- as.numeric(Sys.time())
            }

            if(backend=="doParallel"){
              if (requireNamespace("doParallel", quietly = TRUE)) {
                doParallel::registerDoParallel(cores=ncore)
              } else{
                stop("Usage of doParallel without being installed!")
              }
            } else if(backend=="doMPI"){
              if (requireNamespace("doMPI", quietly = TRUE)) {
                cl <- doMPI::startMPIcluster(count=ncore)
                doMPI::registerDoMPI(cl)
              } else{
                stop("Usage of doMPI without being installed!")
              }
            } else{
              if(verbose) cat("No valid backend specified.\n")
            }
            if (requireNamespace("foreach", quietly = TRUE)) {
            } else{
              stop("Usage of foreach without being installed!")
            }
            Zt <- foreach::foreach(indexb=1:ncore, .combine = "rbind", .multicombine = TRUE,.maxcombine = 1000,
                         .packages="MoBPS") %dopar% {
              Ztpar <- array(0,dim=c(last-first+1, length(batche[[indexb]])))
              sub <- min(batche[[indexb]]) -1
              for(index in batche[[indexb]]){
                k.database <- bve.database[loop_elements[index,3],]
                cindex <- loop_elements[index,1] - sub
                kindex <- loop_elements[index,2]
                Zpar[,cindex] <- base::as.integer(colSums(compute.snps(population, k.database[1],k.database[2],kindex, import.position.calculation=import.position.calculation, from_p=first, to_p=last, decodeOriginsU=decodeOriginsU, bit.storing=bit.storing, nbits=nbits, output_compressed=FALSE)))
              }
              if(Z.integer){
                storage.mode(Ztpar) <- "integer"
              }
              Ztpar
            }
            if(backend=="doParallel"){
              doParallel::stopImplicitCluster()
            } else if(backend=="doMPI"){
              if (requireNamespace("doMPI", quietly = TRUE)) {
                doMPI::closeCluster(cl)
              } else{
                stop("Usage of doMPI without being installed!")
              }

            }
            if(store.comp.times.bve){
              after <- as.numeric(Sys.time())
              zcalc <- zcalc + after - before
            }
          } else{
            if(store.comp.times.bve){
              before <- as.numeric(Sys.time())
            }

            for(index in 1:n.animals){
              k.database <- bve.database[loop_elements[index,3],]
              cindex <- loop_elements[index,1]
              kindex <- loop_elements[index,2]
              Zt[1:(last-first+1), cindex] <- base::as.integer(colSums(compute.snps(population, k.database[1],k.database[2],kindex, import.position.calculation=import.position.calculation, from_p=first, to_p=last,decodeOriginsU=decodeOriginsU, bit.storing=bit.storing, nbits=nbits, output_compressed=FALSE)))
            }
            if(store.comp.times.bve){
              after <- as.numeric(Sys.time())
              zcalc <- zcalc + after - before
            }
          }

          activ_effect <- population$info$effect.p - first + 1
          activ_effect <- activ_effect[activ_effect>0]
          activ_effect <- activ_effect[activ_effect<=last]
          if(remove.effect.position==TRUE){
            if(miraculix && exists("Z.code")){
              if (requireNamespace("miraculix", quietly = TRUE)) {
                Z.code <- miraculix::zeroNthGeno(Z.code, activ_effect)
              }
            } else{
              Zt <- Zt[-activ_effect,]
            }

          } else if(remove.effect.position=="only_effect"){
            if(miraculix && exists("Z.code")){
              if (requireNamespace("miraculix", quietly = TRUE)) {
                Z.code <- miraculix::zeroNthGeno(Z.code, (1:(last-first+1))[-activ_effect])
              }
            } else{
              Zt <- Zt[activ_effect,]
            }
          }

          if(miraculix){
            if (requireNamespace("miraculix", quietly = TRUE)) {
              p_i[first:last] <- miraculix::allele_freq(Z.code)
              A <- A + miraculix::relationshipMatrix(Z.code, centered=TRUE, normalized=FALSE)
            }
          } else if(miraculix.mult){
            storage.mode(Zt) <- "integer"
            p_i[first:last] <- rowSums(Zt[1:(last-first+1),])/ncol(Zt)/2
            if (requireNamespace("miraculix", quietly = TRUE)) {
              Zt_miraculix <- miraculix::genomicmatrix(Zt[1:(last-first+1),])
              A <- A + miraculix::relationshipMatrix(Zt_miraculix, centered=TRUE, normalized=FALSE)
            }
          } else{
            p_i[first:last] <- rowSums(Zt[1:(last-first+1),])/ncol(Zt)/2
            Ztm <- Zt[1:(last-first+1),] - p_i[first:last] * 2
            A <- A + crossprod(Ztm[1:(last-first+1),])
          }

          first <- first + maxZ
          last <- min(maxZ*(index3+1), total_n)
        }
        A <- A / (2 * sum(p_i*(1-p_i)))


      } else{

        if(singlestep.active){
          n.geno <- length(genotype.included)
          n.ped <- n.animals - n.geno
          if(n.ped==0){
            if(verbose) cat("No non genotyped individuals included. Automatically switch to regular GBLUP instead of ssGBLUP\n")
            singlestep.active <- FALSE
          }
          if(n.geno==0){
            z_ped <- z_ped - as.numeric(Sys.time())
            if(verbose) cat("No genotyped individuals included. Automatically switch to pedigree BLUP instead of ssGBLUP\n")
            singlestep.active <- FALSE
            relationship.matrix <- "kinship"
            A <- kinship.exp(population, database=bve.database, depth.pedigree=depth.pedigree, elements = loop_elements_list[[2]], mult = 2, verbose=verbose)
            z_ped <- z_ped + as.numeric(Sys.time())
          }
        }
        if(singlestep.active){
          if(verbose) cat("Start derive Single-step relationship matrix \n")
          if(verbose) cat(paste0("Construct pedigree matrix for ", length(loop_elements_list[[2]]), " individuals.\n"))
          z_ped <- z_ped - as.numeric(Sys.time())

          A_pedigree <-  kinship.exp(population, database=bve.database, depth.pedigree=depth.pedigree, elements = loop_elements_list[[2]], mult = 2, verbose=verbose)
          z_ped <- z_ped + as.numeric(Sys.time())
          if(verbose) cat(paste0("Derived pedigree matrix in  ", round(z_ped, digits=2), " seconds.\n"))
          if(miraculix){
            if (requireNamespace("miraculix", quietly = TRUE)) {

              # Avoid having the full genotype matrix in memory without bitwise storage!
              if(bve.imputation.errorrate>0){
                Z.code.small <- miraculix::genomicmatrix(as.matrix(Z.code)[,genotype.included])
                if(length(to_remove)>0){
                  Z.code.small <- miraculix::genomicmatrix(as.matrix(Z.code.small)[-to_remove,])
                }

              } else{
                Z.code.small <- miraculix::computeSNPS(population, loop_elements[genotype.included,4], loop_elements[genotype.included,5], loop_elements[genotype.included,2], what="geno", output_compressed=TRUE)
                if(length(to_remove)>0){
                  Z.code.small <- miraculix::genomicmatrix(as.matrix(Z.code.small)[-to_remove,])
                }
              }


              p_i <- miraculix::allele_freq(Z.code.small)
              A_geno <- miraculix::relationshipMatrix(Z.code.small, centered=TRUE, normalized=TRUE)
              rm(Z.code.small)
            }
          } else if(miraculix.mult){
            if (requireNamespace("miraculix", quietly = TRUE)) {
              p_i <- rowSums(Zt[,genotype.included])/ncol(Zt[,genotype.included])/2
              Zt_miraculix <- miraculix::genomicmatrix(Zt[,genotype.included])
              A_geno <- miraculix::relationshipMatrix(Zt_miraculix, centered=TRUE, normalized=TRUE)
              rm(Zt_miraculix)
            }
          } else{
            p_i <- rowSums(Zt[,genotype.included])/ncol(Zt[,genotype.included])/2
            Ztm <- Zt[,genotype.included] - p_i * 2
            A_geno <- crossprod(Ztm)/ (2 * sum(p_i*(1-p_i)))
          }

          z_h <- z_h - as.numeric(Sys.time())

          reorder <- c((1:n.animals)[-genotype.included], genotype.included)
          if(verbose) cat("Start deriving of H matrix for", length(genotype.included), "genotyped and", nrow(A_pedigree)-length(genotype.included), "non-genotyped individuals.\n")

          '#
          H_order <- A <- ssGBLUP(A11= A_pedigree[-genotype.included, -genotype.included],
                                  A12 = A_pedigree[-genotype.included, genotype.included],
                                  A22 = A_pedigree[genotype.included, genotype.included], G = A_geno)

          A[genotype.included, genotype.included] <- H_order[(ncol(H_order)-length(genotype.included)+1):ncol(H_order), (ncol(H_order)-length(genotype.included)+1):ncol(H_order)]
          A[-genotype.included, -genotype.included] <- H_order[1:(ncol(H_order)-length(genotype.included)), 1:(ncol(H_order)-length(genotype.included))]
          A[-genotype.included, genotype.included] <- H_order[1:(ncol(H_order)-length(genotype.included)), (ncol(H_order)-length(genotype.included)+1):ncol(H_order)]
          A[genotype.included, -genotype.included] <- H_order[(ncol(H_order)-length(genotype.included)+1):ncol(H_order), 1:(ncol(H_order)-length(genotype.included))]

          rm(A_geno)
          rm(H_order)
          rm(A_pedigree)
          '#

          test1 <- as.numeric(A_geno)
          test2 <- as.numeric(A_pedigree[genotype.included, genotype.included])
          a_step <- mean(test2) - mean(test1)
          A_geno <- A_geno * (1-a_step/2) + a_step # Modification according to Vitezica 2011

          A <- ssGBLUP(A11= A_pedigree[-genotype.included, -genotype.included],
                                  A12 = A_pedigree[-genotype.included, genotype.included],
                                  A22 = A_pedigree[genotype.included, genotype.included], G = A_geno)

          #rm(A_geno)
          #rm(A_pedigree)

          rest <- (1:n.animals)[-genotype.included]
          A[c(genotype.included, rest), c(genotype.included, rest)] <- A[c((ncol(A)-length(genotype.included)+1):ncol(A),1:(ncol(A)-length(genotype.included)) ), c((ncol(A)-length(genotype.included)+1):ncol(A),1:(ncol(A)-length(genotype.included)))]

          z_h <- z_h + as.numeric(Sys.time())
          if(verbose) cat(paste0("Derived H matrix in  ", round(z_h, digits=2), " seconds.\n"))

        } else if(relationship.matrix=="vanRaden"){
          if(miraculix){
            if (requireNamespace("miraculix", quietly = TRUE)) {
              p_i <- miraculix::allele_freq(Z.code) # Noch nicht implementiert?
              A <- miraculix::relationshipMatrix(Z.code, centered=TRUE, normalized=TRUE)
            }
          } else if(miraculix.mult){
            if (requireNamespace("miraculix", quietly = TRUE)) {
              p_i <- rowSums(Zt)/ncol(Zt)/2
              Zt_miraculix <- miraculix::genomicmatrix(Zt)
              A <- miraculix::relationshipMatrix(Zt_miraculix, centered=TRUE, normalized=TRUE)
            }
          } else{
            p_i <- rowSums(Zt)/ncol(Zt)/2
            Ztm <- Zt - p_i * 2
            A <- crossprod(Ztm)/ (2 * sum(p_i*(1-p_i)))
          }
        }

      }


    } else if(relationship.matrix=="CM"){
      #CM SCHAETZER
      Ztm <- rbind(Zt==0, Zt==1, Zt==2)
      A <- crossprod(Ztm) / nrow(Zt)
    } else if(relationship.matrix=="CE"){

      Ztm <- rbind(Zt==0, Zt==1, Zt==2)
      A <- crossprod(Ztm)
      A <- (A^2 - 0.5*A)/(nrow(Zt)^2)
    } else if(relationship.matrix=="CE2"){
      A2 <- crossprod(Zt)
      A <- 0.5 * A2 * A2 + 0.5 * crossprod(Zt*Zt)
      A <- A/mean(diag(A))
    } else if(relationship.matrix=="CE3"){
      A2 <- crossprod(Zt)
      A <- 0.5 * A2 * A2 - 0.5 * crossprod(Zt*Zt)
      A <- A/mean(diag(A))
    } else if(relationship.matrix=="non_stand"){
      A <- crossprod(Zt) / nrow(Zt)
    } else if(relationship.matrix=="vanRaden2"){
      p_i <- rowSums(Zt)/ncol(Zt)/2
      Ztm <- Zt - p_i * 2
      A <- crossprod(Ztm)/ (2 * sum(p_i*(1-p_i)))
    }

    if(store.comp.times.bve){
      comp.times.bve[3] <- as.numeric(Sys.time())
    }

    sigma.a.hat <- numeric(length(sigma.g))
    sigma.e.hat <- sigma.e

    prev_rest_take <- NULL

    n.rep <- 0
    if(length(bve.database)==length(bve.insert.database) && prod(bve.database==bve.insert.database)==1 && nrow(loop_elements_list[[3]])==0){
      bve.insert <- rep(TRUE, n.animals)
      bve.insert.copy <- NULL
    } else{
      bve.insert <- rep(FALSE, n.animals)
      before <- 0
      for(index in 1:nrow(bve.insert.database)){
        add_insert <- loop_elements[,2] <= bve.insert.database[index,4] & loop_elements[,2] >= bve.insert.database[index,3] & loop_elements[,4] == bve.insert.database[index,1] & loop_elements[,5] == bve.insert.database[index,2]
        bve.insert[add_insert] <- TRUE
      }

      loop_elements_copy <- loop_elements_list[[3]]
      n.rep <- nrow(loop_elements_copy)
      bve.insert.copy <- rep(FALSE, n.rep)
      if(n.rep>0){
        for(index in 1:nrow(bve.insert.database)){
          add_insert <- loop_elements_copy[,2] <= bve.insert.database[index,4] & loop_elements_copy[,2] >= bve.insert.database[index,3] & loop_elements_copy[,4] == bve.insert.database[index,1] & loop_elements_copy[,5] == bve.insert.database[index,2]
          bve.insert.copy[add_insert] <- TRUE
        }
      }

    }
    # Breeding value estimation (Solving of the mixed-model / bayesian generalized linear regression) - all single trait models except multi-sommer

    if(length(bve.ignore.traits)>0){
      text <- population$info$trait.name[bve.ignore.traits][1]
      for(write in population$info$trait.name[bve.ignore.traits][-1]){
        text <- paste0(text, ", ", write)
      }
      if(verbose) cat(paste0("BVE for: ", text, " has been skipped.\n"))
    }


    combi_trait <- NULL
    skip_trait <- NULL
    replace_trait <- NULL

    if(length(population$info$trait.combine.name)>0){

      activ_bves <- (1:population$info$bv.nr)[bve.keeps]

      for(check_combi in 1:length(population$info$trait.combine.included)){

        to_combine <- intersect(population$info$trait.combine.included[[check_combi]], activ_bves)

        if(length(to_combine)>1){
          if(verbose){
            cat(paste0("Perform a joint BVE for trait ", population$info$trait.combine.name[[check_combi]], ".\n"))

            text <- "This includes traits: "
            for(combi2 in to_combine){
              text <- paste0(text, population$info$trait.name[combi2], " ")
            }
          }

          y_temp <- y[,to_combine]
          ny_temp <- y_obs[,to_combine]
          y_temp[is.na(y_temp)] <- 0

          new_y <- rowSums(y_temp *ny_temp)
          new_ny <- rowSums(ny_temp)
          y[,to_combine] <- new_y/new_ny
          y_obs[,to_combine] <- new_ny

          combi_trait <- c(combi_trait, to_combine[1])
          skip_trait <- c(skip_trait, to_combine[-1])
          replace_trait <- c(replace_trait, rep(to_combine[1], length(to_combine)-1))
        }
      }

    }

    beta_hat <-  colMeans(y, na.rm=TRUE) # Rest faellt weg! (X' R^-1 X)^-1 X' R^-1 y
    if(sum(is.na(beta_hat))>0){
      if(verbose) cat("No phenotypes available for traits:", population$info$trait.name[which(is.na(beta_hat))],"\n")
      if(verbose) cat("Set all breeding value estimates for these trait(s) to 0. \n")
      beta_hat[is.na(beta_hat)] <- 0
    }

    for(bven in (1:population$info$bv.nr)[bve.keeps]){

      if(sum(bven==skip_trait)==0){

        if(forecast.sigma.g){
          sigma.g[bven] <- stats::var(y_real2[,bven], na.rm = TRUE)
        }
        if(estimate.pheno.var){
          sigma.e.hat[bven] <- max(0, stats::var(y[,bven], na.rm=TRUE) - sigma.g[bven])
        }
        sigma.a.hat[bven] <- sigma.g[bven]
        if(estimate.add.gen.var){
          sigma.a.hat[bven] <- max(min(stats::lm(y[,bven]~y_parent[,bven])$coefficients[2],1),0.001) * sigma.e.hat[bven]
        }
        if(mas.bve){
          if(length(mas.markers)==0){
            if(length(population$info$real.bv.add[[bven]])>24){
              if(verbose) cat("No markers for MAS provided. Use 5 true effect markers with highest effect.\n")
              temp1 <- sort(rowSums(abs(population$info$real.bv.add[[bven]][,3:5])), index.return=TRUE, decreasing=TRUE)$ix[1:mas.number]
              active.snp <- population$info$real.bv.add[[bven]][temp1,]
              mas.markers.temp <- c(0,population$info$cumsnp)[active.snp[,2]]+ active.snp[,1]
            } else{
              mas.markers.temp <- sample(population$info$effect.p, min(mas.number, length(population$info$effect.p)))
            }

          } else{
            mas.markers.temp <- mas.markers
          }

          mas_geno <- as.matrix(Z.code)[mas.markers.temp,]
          if(length(mas.effects)==0){
            model <- stats::lm(y~t(mas_geno))
            y_hat[,bven] <- model$fitted.values
          } else{
            mas.effects <- rep(mas.effects, length.out=length(mas.markers.temp))
            y_hat[,bven] <- mas.effects %*% mas_geno
          }

        } else if(BGLR.bve){

          if(miraculix){
            Zt <- as.matrix(Z.code)
          }
          Zt <- t(scale(t(Zt), center=TRUE, scale=FALSE))
          fixed <- which(is.na(Zt[,1]))
          if(length(fixed)>0){
            Zt <- Zt[-fixed,]
          }
          if(BGLR.model=="RKHS"){
            ETA <- list(list(K=A, model='RKHS'))
          } else if(BGLR.model=="BayesA"){
            ETA <- list(list(X=t(Zt), model='BayesA'))
          } else if(BGLR.model=="BayesB"){
            ETA <- list(list(X=t(Zt), model='BayesB'))
          } else if(BGLR.model=="BayesC"){
            ETA <- list(list(X=t(Zt), model='BayesC'))
          } else if(BGLR.model=="BL"){
            ETA <- list(list(X=t(Zt), model='BL'))
          } else if(BGLR.model=="BRR"){
            ETA <- list(list(X=t(Zt), model='BRR'))
          }

          if(min(y[!is.na(y[,bven]),bven])==max(y[!is.na(y[,bven]),bven])){
            y_hat[,bven] <- rep(max(y[!is.na(y[,bven]),bven]), length(y[,bven]))
          } else{
            if(BGLR.save.random){
              BGLR.save <- paste0(BGLR.save, sample(1:10000,1))
            }
            if(requireNamespace("BGLR", quietly = TRUE)){
              fm <- BGLR::BGLR(y=y[,bven], ETA=ETA, nIter=BGLR.iteration, burnIn=BGLR.burnin, saveAt=BGLR.save, verbose=BGLR.print)
            } else{
              stop("Use of BGLR without being installed!")
            }
            y_hat[,bven] <- fm$yHat
          }
        } else if(rrblup.bve){

          check <- sum(is.na(y[,bven]))
          if(verbose) cat(paste0(length(y[,bven]) - check, " phenotyped individuals in BVE (Trait: ", population$info$trait.name[bven],").\n"))
          if(verbose) cat(paste0(length(y[,bven]), " individuals considered in BVE.\n"))
          if(check >= (length(y[,bven])-1)){
            if(verbose) cat(paste0("No phenotyped individuals for trait ", population$info$trait.name[bven], "\n"))

            if(verbose) cat(paste0("Skip this BVE.\n"))
            next
          }

          t1 <- as.numeric(Sys.time())
          if (requireNamespace("rrBLUP", quietly = TRUE)) {
            test <- rrBLUP::mixed.solve(y[,bven], K = A, method="REML", bounds = c(1e-9,1e9))
          } else{
            stop("Use of rrBLUP without being installed!")
          }
          t2 <- as.numeric(Sys.time())

          if(verbose) cat(paste0(round(t2-t1, digits=2), " seconds for BVE.\n"))
          y_hat[,bven] <- as.numeric(test$beta) + test$u

        } else if(emmreml.bve){
          check <- sum(is.na(y[,bven]))
          if(verbose) cat(paste0(length(y[,bven]) - check, " phenotyped individuals in BVE (Trait: ", population$info$trait.name[bven],").\n"))
          if(verbose) cat(paste0(length(y[,bven]), " individuals considered in BVE.\n"))
          if(check >= (length(y[,bven])-1)){
            if(verbose) cat(paste0("No phenotyped individuals for trait ", population$info$trait.name[bven], "\n"))
            if(verbose) cat(paste0("Skip this BVE."))
            next
          }
          if(check>0){
            if(verbose) cat(paste0("Breeding value estimation with ", check, " NA phenotypes! EMMREML does not support this!\n"))
            if(verbose) cat(paste0("No estimation is performed to NA individuals. \n"))
            take <- which(!is.na(y[,bven]))
          } else{
            take <- 1:length(y[,bven])
          }
          n <- length(y[take,bven])
          p <- nrow(Zt)
          stopifnot(n == ncol(Zt[,take]))

          if(requireNamespace("EMMREML", quietly = TRUE)){
            fm <- EMMREML::emmreml(
              y[take,bven],
              matrix(1,nrow=n),
              diag(n),
              A[take,take])
          } else{
            stop("Usage of EMMREML without being installed!")
          }


          y_hat[take,bven] <- as.numeric(fm$uhat) + as.numeric(fm$betahat)
          if(estimate.u){
            u_hat <- cbind(u_hat, alpha_to_beta(drop(fm$uhat),A[take,take],t(Zt[,take])), deparse.level = 0)
          }

        } else if(sommer.bve){

          check <- sum(is.na(y[,bven]))
          if(verbose) cat(paste0(length(y[,bven]) - check, " phenotyped individuals in BVE (Trait: ", population$info$trait.name[bven],").\n"))
          if(verbose) cat(paste0(length(y[,bven]), " individuals considered in BVE.\n"))
          if(check >= (length(y[,bven])-1)){
            if(verbose) cat(paste0("No phenotyped individuals for trait ", population$info$trait.name[bven], "\n"))
            if(verbose) cat(paste0("Skip this BVE.\n"))
            next
          }

          traitnames <- (paste0("name", 1:ncol(y)))
          traitnames[bven] <- "name"
          traitnames <- as.factor(traitnames)
          colnames(y) <- traitnames
          id <- as.factor(paste0("P", 1:nrow(y)))
          y_som <- data.frame(y, id)
          rownames(y_som) <- id
          colnames(A) <- rownames(A) <- id

          if (requireNamespace("sommer", quietly = TRUE)) {
            test <- sommer::mmer(name ~1, random=~sommer::vs(id, Gu=A), rcov = ~units, data=y_som)
          } else{
            stop("Use of sommer without being installed!")
          }

          y_hat[sort(as.character(id), index.return=TRUE)$ix,bven] <- test$U[[1]][[1]] + as.numeric(test$Beta[3])

        } else if(sommer.multi.bve){

          check <- sum(is.na(y))

          if(verbose) cat(paste0(length(y[,bven]) - check, " phenotyped individuals in BVE (Trait: ", population$info$trait.name[bven],").\n"))
          if(verbose) cat(paste0(length(y[,bven]), " individuals considered in BVE.\n"))
          if(check >= (length(y[,bven])-1)){
            if(verbose) cat(paste0("No phenotyped individuals for multi-trait mixed model\n"))
            if(verbose) cat(paste0("Skip this BVE.\n"))
            next
          }


          if(bven==max((1:population$info$bv.nr)[bve.keeps])){

            traitnames <- paste0("name", 1:ncol(y))
            colnames(y) <- as.factor(traitnames)
            id <- as.factor(paste0("P", 1:nrow(y)))
            y_som <- data.frame(y, id)
            rownames(y_som) <- id
            colnames(A) <- rownames(A) <- id
            text <- "cbind("
            for(index in 1:length(traitnames)){
              if(index==length(traitnames)){
                text <- paste0(text, traitnames[index], ")")
              } else{
                text <- paste0(text, traitnames[index], ",")
              }

            }
            text <- paste0("sommer::mmer(",text,"~1, random=~sommer::vs(id, Gu=A, Gtc=sommer::unsm(bven)), rcov = ~sommer::vs(units, Gtc=diag(bven)), data=y_som)")
            if (requireNamespace("sommer", quietly = TRUE)) {
              test <- eval(parse(text=text))
            } else{
              stop("Use of sommer without being installed!")
            }

            for(bven1 in 1:population$info$bv.nr){
              y_hat[sort(as.character(id), index.return=TRUE)$ix,bven1] <- test$U[[1]][[bven1]]
            }


            if(estimate.u){
              if(verbose) cat("U estimation not available in sommer")
            }
          }
        } else if(sigma.e[bven]>0 || sum(is.na(y[,bven]>0))){
          # sigma.a.hat / sigma.e.hat are modified to avoid h=0, h=1 cases (numeric instability)
          if(sigma.a.hat[bven]==0){
            if(sigma.e.hat[bven]>0){
              sigma.a.hat[bven] <- sigma.e.hat[bven] * 0.001
            } else{
              sigma.e.hat[bven] <- 1
              sigma.a.hat[bven] <- 1000
            }
          }
          if(sigma.e.hat[bven]==0){
            sigma.e.hat[bven] <- sigma.a.hat[bven] * 0.001
          }
          bve.direct.est.now <- bve.direct.est
          check <- sum(is.na(y[,bven]))

          u_hat_possible <- TRUE
          # This will be more complicated in case other fixed-effects are added to the model !
          multi <- y[,bven] - X %*% beta_hat[bven]

          rrblup.required <- FALSE

          if(verbose) cat(paste0(length(y[,bven]) - check, " phenotyped individuals in BVE (Trait: ", population$info$trait.name[bven],").\n"))
          if(verbose) cat(paste0(length(y[,bven]), " individuals considered in BVE.\n"))
          if(check >= (length(y[,bven])-1)){
            if(verbose) cat(paste0("No phenotyped individuals for trait ", population$info$trait.name[bven], "\n"))
            if(verbose) cat(paste0("Skip this BVE.\n"))
            next
          }

          # take Individuals used to trait the mixed model on
          # take2 individuals for which to enter a breeding value
          # take3 individuals to used for estimation of allele effects (rrBLUP)
          if(check>0){
            if(sum(genotyped==1 & is.na(y[,bven]))==sum(genotyped==1)){
              if(bve.direct.est.now==FALSE){
                if(verbose) cat("No genotyped and phenotyped individuals. Application of rrBLUP not possible!\n")
                if(verbose) cat("Assume non phenotyped individuals to have average phenotype.\n")
                if(bve.0isNA && sum(is.na(multi))>0){
                  multi[is.na(multi)] <- 0
                }
                take <- take2 <- 1:length(y[,bven])
              } else{
                take <- which(!is.na(y[,bven]))
                take2 <- 1:length(y[,bven])
              }


            } else if(sum(genotyped==1 & is.na(y[,bven]))>0){
              if(verbose) cat("Some genotyped individuals without phenotype.\n")
              take <- which(!is.na(y[,bven])) # individuals for GBLUP
              if(!bve.direct.est.now){
                rrblup.required <- TRUE
                take2 <- which(!is.na(y[,bven]) & genotyped==1) # individuals for rrBLUP
              } else if(estimate.u){
                take2 <- which(!is.na(y[,bven]) & genotyped==1) # individuals for rrBLUP
              } else{
                take2 <- 1:length(y[,bven])
              }


              bve.insert.full <- bve.insert
              if(length(bve.insert.copy)>0){

                bve.insert.full[loop_elements_copy[bve.insert.copy,6]] <- TRUE
              }
              take3 <- which(is.na(y[,bven]) & genotyped==1 & bve.insert.full) # individuals to estimate via rrBLUP
              if(length(take3)>0 && !bve.direct.est.now){
                if(verbose) cat("Use rrBLUP to estimate breeding value for those individuals!\n")
              }
            } else{
              take <- take2 <- which(!is.na(y[,bven]))
              if(bve.direct.est.now){
                take2 <- 1:length(y[,bven])
              }
            }
          } else{
            take <- take2 <- 1:length(y[,bven])
          }

          if(length(take)==length(take2)){
            bve.direct.est.now <- FALSE
          }



          if(store.comp.times.bve){
            before <- as.numeric(Sys.time())
          }
          # Solve mixed model which assumed known heritability
          # This is not 100% accurate but massively reduces computing time and should be ok for large scale breeding programs
          t1 <- as.numeric(Sys.time())

          if(length(take)==nrow(A)){
            skip.copy = TRUE
          } else{
            skip.copy = FALSE
          }

          if(bve.per.sample.sigma.e & max(y_obs, na.rm = TRUE)>1){
            if(length(repeatability) < bven || length(heritability) <bven || repeatability[bven]==heritability[bven]){
              heri_factor <- 1 / y_obs[take,bven]
            } else{
              heri_factor <- (sigma.e.perm + sigma.e.rest / y_obs[take,bven]) / (sigma.e.perm + sigma.e.rest)
              if(heri_factor< 0.01 || heri_factor>1){
                warning("Breeding value estimation with multiple observations with extrem residual variance scaling!")
              }
            }


          } else{
            heri_factor <- 1
          }

          if(prod(heri_factor==1)!=1 || is.function(bve.solve)){
            miraculix_possible <- FALSE
          } else{
            miraculix_possible <- TRUE
          }


          if(estimate.u || rrblup.required){
            if(calculate.reliability){
              if(verbose) cat("Reliabilities are currently only calculated for phenotyped individuals. Extension according to vanRaden 2008 planned.\n")
              if(verbose) cat("Use of Z2 Z instead of ZZ (C instead of G). \n")
              if(skip.copy){
                GR1 <- chol2inv(chol(add.diag(A,sigma.e.hat[bven] / sigma.a.hat[bven] * heri_factor)))
              } else{
                GR1 <- chol2inv(chol(add.diag(A[take,take],sigma.e.hat[bven] / sigma.a.hat[bven] * heri_factor)))
              }

              Rest_term <- GR1 %*% multi[take]
              y_hat[take,bven] <- A[take,take] %*% Rest_term  + beta_hat[bven]
              y_reli[take,bven] <- diag( A[take,take] %*% GR1 %*% A[take,take])
            } else if(miraculix && miraculix.chol && miraculix_possible){
              if (requireNamespace("miraculix", quietly = TRUE) ) {

                if(skip.copy){
                  temp1 <- miraculix::solveRelMat(A, sigma.e.hat[bven] / sigma.a.hat[bven], multi[take], beta_hat[bven], destroy_A = miraculix.destroyA)
                } else{
                  temp1 <- miraculix::solveRelMat(A[take,take], sigma.e.hat[bven] / sigma.a.hat[bven], multi[take], beta_hat[bven], destroy_A = miraculix.destroyA)
                }
              }
              Rest_term <- temp1[[1]]
              y_hat[take,bven] <- temp1[[2]]
            } else{
              if(skip.copy){
                Rest_term <- (chol2inv(chol(add.diag(A, sigma.e.hat[bven] / sigma.a.hat[bven] * heri_factor))) %*% multi[take])
                y_hat[take,bven] <- A %*% Rest_term  + beta_hat[bven]
              } else{
                Rest_term <- (chol2inv(chol(add.diag(A[take,take], sigma.e.hat[bven] / sigma.a.hat[bven] * heri_factor))) %*% multi[take])
                y_hat[take,bven] <- A[take,take] %*% Rest_term  + beta_hat[bven]
              }

            }
          } else{
            if(calculate.reliability){

              if(skip.copy){
                GR1 <- chol2inv(chol(add.diag(A, sigma.e.hat[bven] / sigma.a.hat[bven] * heri_factor)))
              } else{
                GR1 <- chol2inv(chol(add.diag(A[take,take], sigma.e.hat[bven] / sigma.a.hat[bven] * heri_factor)))
              }

              if(bve.direct.est.now){
                y_hat[take2,bven] <- A[take2,take] %*% (GR1 %*% multi[take])  + beta_hat[bven]
                y_reli[take2,bven] <- diag( A[take2,take] %*% GR1 %*% A[take,take2])
              } else{
                y_hat[take,bven] <- A[take,take] %*% (GR1 %*% multi[take])  + beta_hat[bven]
                y_reli[take,bven] <- diag( A[take,take] %*% GR1 %*% A[take,take])
              }

            } else if(miraculix && miraculix.chol && miraculix_possible){
              if (requireNamespace("miraculix", quietly = TRUE)) {
                if(bve.direct.est.now){
                  y_hat[take2,bven] <-  A[take2,take] %*% miraculix::solveRelMat(A[take,take], sigma.e.hat[bven] / sigma.a.hat[bven], multi[take],betahat = NULL, destroy_A = miraculix.destroyA) + beta_hat[bven]
                } else{
                  if(skip.copy){
                    y_hat[,bven] <- miraculix::solveRelMat(A, sigma.e.hat[bven] / sigma.a.hat[bven], multi[take],beta_hat[bven], destroy_A = miraculix.destroyA)[[2]]
                  } else{
                    y_hat[take,bven] <- miraculix::solveRelMat(A[take,take], sigma.e.hat[bven] / sigma.a.hat[bven], multi[take],beta_hat[bven], destroy_A = miraculix.destroyA)[[2]]
                  }
                }

              }
            } else{
              if(bve.direct.est.now){

                if(is.function(bve.solve)){
                  y_hat[take2,bven] <- A[take2,take] %*% bve.solve(add.diag(A[take,take], sigma.e.hat[bven] / sigma.a.hat[bven] * heri_factor), b= multi[take]) + beta_hat[bven]
                } else{
                  y_hat[take2,bven] <- A[take2,take] %*% (chol2inv(chol(add.diag(A[take,take], sigma.e.hat[bven] / sigma.a.hat[bven] * heri_factor))) %*% multi[take]) + beta_hat[bven]
                }

              } else{
                if(skip.copy){
                  if(is.function(bve.solve)){
                    y_hat[,bven] <- A %*% bve.solve(A=add.diag(A,sigma.e.hat[bven] / sigma.a.hat[bven] * heri_factor), b=multi[take]) + beta_hat[bven]
                  } else{
                    y_hat[,bven] <- A %*% (chol2inv(chol(add.diag(A,sigma.e.hat[bven] / sigma.a.hat[bven] * heri_factor))) %*% multi[take]) + beta_hat[bven]

                  }

                } else{

                  if(is.function(bve.solve)){
                    y_hat[take,bven] <- A[take,take] %*% bve.solve(A=add.diag(A[take,take],sigma.e.hat[bven] / sigma.a.hat[bven] * heri_factor), b=multi[take]) + beta_hat[bven]
                  } else{
                    y_hat[take,bven] <- A[take,take] %*% (chol2inv(chol(add.diag(A[take,take],sigma.e.hat[bven] / sigma.a.hat[bven] * heri_factor))) %*% multi[take]) + beta_hat[bven]
                  }


                }
              }
            }
          }

          t2 <- as.numeric(Sys.time())
          if(verbose) cat(paste0(round(t2-t1, digits=2), " seconds for BVE.\n"))

          if(store.comp.times.bve){
            after <- as.numeric(Sys.time())
            z_chol <- after - before
          }

          # Estimation of marker effects via rrBLUP

          if(estimate.u || rrblup.required){

            while(bven>1 && (length(u_hat)==0 || ncol(u_hat)<(bven-1))){
              u_hat <- cbind(u_hat,rep(0, sum(population$info$snp)), deparse.level = 0)
            }
            rest_take <- which(duplicated(c(take,take2))[-(1:length(take))])

            if(sequenceZ){
              total_n <- sum(population$info$snp)
              u_hat_new <- numeric(total_n)
              first <- 1
              last <- min(maxZ, total_n)
              for(index3 in 1:ceiling(total_n/maxZ)){
                if(miraculix){
                  if(store.comp.times.bve){
                    before <- as.numeric(Sys.time())
                  }
                  if (requireNamespace("miraculix", quietly = TRUE)) {
                    Z.code2 <- miraculix::computeSNPS(population, loop_elements[take2,4], loop_elements[take2,5], loop_elements[take2,2],
                                                      from_p=first, to_p=last, what="geno", output_compressed=TRUE
                    )
                  }

                  if(store.comp.times.bve){
                    after <- as.numeric(Sys.time())
                    zcalc <- zcalc + after - before
                  }
                } else if(ncore>1){
                  if(store.comp.times.bve){
                    before <- as.numeric(Sys.time())
                  }

                  if(backend=="doParallel"){
                    if (requireNamespace("doParallel", quietly = TRUE)) {
                      doParallel::registerDoParallel(cores=ncore)
                    } else{
                      stop("Use of doParallel without being installed!")
                    }
                  } else if(backend=="doMPI"){
                    if (requireNamespace("doMPI", quietly = TRUE)) {
                      cl <- doMPI::startMPIcluster(count=ncore)
                      doMPI::registerDoMPI(cl)
                    } else{
                      stop("Usage of doMPI without being installed!")
                    }
                  } else{
                    if(verbose) cat("No valid backend specified.\n")
                  }
                  if (requireNamespace("foreach", quietly = TRUE)) {
                  } else{
                    stop("Usage of foreach without being installed!")
                  }
                  Zt <- foreach::foreach(indexb=1:ncore, .combine = "cbind", .multicombine = TRUE,.maxcombine = 1000,
                                         .packages="MoBPS") %dopar% {
                                           Ztpar <- array(0,dim=c(sum(population$info$snp), length(batche[[indexb]])))
                                           sub <- min(batche[[indexb]]) -1
                                           for(index in batche[[indexb]]){
                                             k.database <- bve.database[loop_elements[index,3],]
                                             cindex <- loop_elements[index,1] - sub
                                             kindex <- loop_elements[index,2]
                                             Ztpar[,cindex] <- base::as.integer(colSums(compute.snps(population, k.database[1],k.database[2],kindex, import.position.calculation=import.position.calculation, from_p=first, to_p=last, decodeOriginsU=decodeOriginsU, bit.storing=bit.storing, nbits=nbits, output_compressed=FALSE)))
                                           }
                                           if(Z.integer){
                                             storage.mode(Ztpar) <- "integer"
                                           }
                                           Ztpar
                                         }
                  if(backend=="doParallel"){
                    doParallel::stopImplicitCluster()
                  } else if(backend=="doMPI"){
                    if (requireNamespace("doMPI", quietly = TRUE)) {
                      doMPI::closeCluster(cl)
                    } else{
                      stop("Usage of doMPI without being installed!")
                    }

                  }
                  if(store.comp.times.bve){
                    after <- as.numeric(Sys.time())
                    zcalc <- zcalc + after - before
                  }
                } else{
                  if(store.comp.times.bve){
                    before <- as.numeric(Sys.time())
                  }
                  for(index in 1:n.animals){
                    k.database <- bve.database[loop_elements[index,3],]
                    cindex <- loop_elements[index,1]
                    kindex <- loop_elements[index,2]
                    Zt[1:(last-first+1), cindex] <- base::as.integer(colSums(compute.snps(population, k.database[1],k.database[2],kindex, import.position.calculation=import.position.calculation, from_p=first, to_p=last,decodeOriginsU=decodeOriginsU, bit.storing=bit.storing, nbits=nbits, output_compressed=FALSE)))
                  }
                  if(store.comp.times.bve){
                    after <- as.numeric(Sys.time())
                    zcalc <- zcalc + after - before
                  }
                }

                if(store.comp.times.bve){
                  before <- as.numeric(Sys.time())
                }

                if(!(length(rest_take)==length(prev_rest_take) && prod(rest_take==prev_rest_take)==1) && !fast.uhat){
                  if (requireNamespace("MASS", quietly = TRUE)) {
                    A1 <- MASS::ginv(A[take2[rest_take], take2[rest_take]])
                  } else{
                    stop("Use of MASS without being installed!")
                  }
                }
                if(relationship.matrix!="vanRaden"){
                  if(miraculix){
                    if (requireNamespace("miraculix", quietly = TRUE)) {
                      p_i <- miraculix::allele_freq(Z.code2)
                    }
                  } else{
                    p_i <- rowSums(Zt[,take2[rest_take]])/2
                  }
                }
                if(miraculix){

                  if (requireNamespace("miraculix", quietly = TRUE)) {
                    if(fast.uhat){
                      u_hat_new[first:last] <- 1/ 2 / sum(p_i*(1-p_i))* miraculix::genoVector(Z.code2, Rest_term[rest_take])
                    } else{
                      u_hat_new[first:last] <- 1/ 2 / sum(p_i*(1-p_i))* miraculix::genoVector(Z.code2, (A1 %*% (y_hat[take2[rest_take],bven] - beta_hat[bven])))
                    }
                  }
                } else if(miraculix.mult){
                  if (requireNamespace("miraculix", quietly = TRUE)) {
                    if(fast.uhat){
                      u_hat_new[first:last] <- 1/ 2 / sum(p_i*(1-p_i))* miraculix::genoVector(miraculix::genomicmatrix(Zt[,take2[rest_take]]), Rest_term[rest_take])
                    } else{
                      u_hat_new[first:last] <- 1/ 2 / sum(p_i*(1-p_i))* miraculix::genoVector(miraculix::genomicmatrix(Zt[,take2[rest_take]]), (A1 %*% (y_hat[take2[rest_take],bven] - beta_hat[bven])))
                    }
                  }
                } else{
                  if(fast.uhat){
                    u_hat_new[first:last] <- 1/ 2 / sum(p_i*(1-p_i))*(Zt[,rest_take] %*% Rest_term[rest_take])
                  } else{
                    u_hat_new[first:last] <- 1/ 2 / sum(p_i*(1-p_i))*(Zt[,rest_take] %*% (A1 %*% (y_hat[rest_take,bven] - beta_hat[bven])))
                  }
                }

                if(store.comp.times.bve){
                  after <- as.numeric(Sys.time())
                  z_uhat <- z_uhat + after - before
                }


                first <- first + maxZ
                last <- min(maxZ*(index3+1), total_n)
              }
              u_hat <- cbind(u_hat, u_hat_new, deparse.level = 0)

            } else{

              if(store.comp.times.bve){
                before <- as.numeric(Sys.time())
              }
              rest_take <- which(duplicated(c(take,take2))[-(1:length(take))])
              if(!(length(rest_take)==length(prev_rest_take) && prod(rest_take==prev_rest_take)==1)){
                if(miraculix && length(take2)!=nrow(loop_elements)){
                  if (requireNamespace("miraculix", quietly = TRUE)) {
                    Z.code2 <- miraculix::computeSNPS(population, loop_elements[take2,4], loop_elements[take2,5], loop_elements[take2,2], what="geno",
                                                      output_compressed = TRUE)
                  }
                } else if(length(take2)!=nrow(loop_elements)){
                  Zt2 <- Zt[,take2[rest_take]]
                } else if(miraculix){
                  Z.code2 <- Z.code
                } else {
                  Zt2 <- Zt
                }
                if(!fast.uhat){
                  if (requireNamespace("MASS", quietly = TRUE)) {
                    A1 <- MASS::ginv(A[take2[rest_take], take2[rest_take]])
                  } else{
                    stop("Use of MASS without being installed!")
                  }
                }

                prev_rest_take <- rest_take
              } else if(length(prev_rest_take)==0){
                if(miraculix){
                  Z.code2 <- Z.code
                } else{
                  Zt2 <- Zt
                }

              }

              if(relationship.matrix!="vanRaden"){
                if(miraculix){
                  if (requireNamespace("miraculix", quietly = TRUE)) {
                    p_i <- miraculix::allele_freq(Z.code2)
                  }
                } else{
                  p_i <- rowSums(Zt[,take2[rest_take]])/2
                }
              }

              if(miraculix){
                if (requireNamespace("miraculix", quietly = TRUE)) {
                  if(fast.uhat){
                    u_hat <- cbind(u_hat, 1/ 2 / sum(p_i*(1-p_i))* miraculix::genoVector(Z.code2, Rest_term[rest_take]), deparse.level = 0)
                  } else{
                    u_hat <- cbind(u_hat, 1/ 2 / sum(p_i*(1-p_i)) * miraculix::genoVector(Z.code2, A1 %*% (y_hat[take2[rest_take],bven] - beta_hat[bven])), deparse.level = 0)
                  }
                }
              } else if(miraculix.mult){
                if (requireNamespace("miraculix", quietly = TRUE)) {
                  if(fast.uhat){
                    u_hat <- cbind(u_hat, 1/ 2 / sum(p_i*(1-p_i))* miraculix::genoVector(miraculix::genomicmatrix(Zt[,take2[rest_take]]), Rest_term[rest_take]), deparse.level = 0)
                  } else{
                    u_hat <- cbind(u_hat, 1/ 2 / sum(p_i*(1-p_i))* miraculix::genoVector(miraculix::genomicmatrix(Zt[,take2[rest_take]]), A1 %*% (y_hat[take2[rest_take],bven] - beta_hat[bven])), deparse.level = 0)
                  }
                }
              } else{
                if(fast.uhat){
                  u_hat <- cbind(u_hat, 1/ 2 / sum(p_i*(1-p_i))*(Zt[,take2[rest_take]] %*% Rest_term[rest_take]), deparse.level = 0)
                } else{
                  u_hat <- cbind(u_hat, 1/ 2 / sum(p_i*(1-p_i))*(Zt[,take2[rest_take]] %*% (A1 %*% (y_hat[take2[rest_take],bven] - beta_hat[bven]))), deparse.level = 0)
                }
              }

              if(store.comp.times.bve){
                after <- as.numeric(Sys.time())
                z_uhat <- z_uhat + after - before
              }

            }

            if(rrblup.required){
              if(sequenceZ){
                stop("Not implemented!")
              } else{
                if(miraculix){
                  if (requireNamespace("miraculix", quietly = TRUE)) {
                    y_hat[take3,bven] <- u_hat[,bven] %*% (as.matrix(Z.code)[,take3]-2*p_i) + beta_hat[bven]
                  }
                } else {
                  y_hat[take3,bven] <- u_hat[,bven] %*% Ztm[,take3] + beta_hat[bven]
                }
              }

            }
          }
        } else{
          y_hat[,bven] <- y[,bven]
        }

      } else{

        which_skip <- which(skip_trait==bven)

        y_hat[,bven] <- y_hat[, replace_trait[which_skip]]
      }


      if(store.comp.times.bve){
        comp.times.bve[4] <- as.numeric(Sys.time())
      }





      for(index in (1:nrow(loop_elements))[bve.insert]){
        population$breeding[[loop_elements[index,4]]][[loop_elements[index,5]+2]][bve.keeps, loop_elements[index,2]] <- y_hat[index,bve.keeps]
      }
      if(calculate.reliability){
        for(index in (1:nrow(loop_elements))[bve.insert]){
          population$breeding[[loop_elements[index,4]]][[loop_elements[index,5]+18]][bve.keeps, loop_elements[index,2]] <- y_reli[index,bve.keeps]
        }
      }

      if(n.rep>0){
        for(index in (1:nrow(loop_elements_copy))[bve.insert.copy]){
          if(length(stay.loop.elements)>0){
            non_copy <- which(stay.loop.elements==loop_elements_copy[index,6])
          } else{
            non_copy <- loop_elements_copy[index,6]
          }
          if(length(non_copy)==1){
            population$breeding[[loop_elements_copy[index,4]]][[loop_elements_copy[index,5]+2]][bve.keeps, loop_elements_copy[index,2]] <- y_hat[non_copy,bve.keeps]

          }
        }
      }



      #  GWAS CODE NOT WRITEN FOR PARALLEL COMPUTING
      if(gwas.u){

        if(y.gwas.used=="pheno"){
          y_gwas <- y
        } else if(y.gwas.used=="bv"){
          y_gwas <- y_real
        } else if(y.gwas.used=="bve"){
          y_gwas <- y_hat
        }

        if(nrow(gwas.database)!=nrow(bve.database) || prod(gwas.database==bve.database)==0){

          loop_elements_gwas_list <- derive.loop.elements(population=population, bve.database=bve.database,
                                                bve.class=bve.class, bve.avoid.duplicates=bve.avoid.duplicates,
                                                store.adding=TRUE)
          loop_elements_gwas <- loop_elements_gwas_list[[1]]
          n.animals.gwas <- nrow(loop_elements_gwas)
        }
        if(gwas.group.standard){
          gwas_start <- loop_elements_gwas_list[[2]]
        }
        if(sequenceZ){
          if(nrow(gwas.database)!=nrow(bve.database) || prod(gwas.database==bve.database)==0){
            if(maxZtotal>0){
              maxZ <- floor(maxZtotal / n.animals.gwas)
            }
            if(ncore<=1 && miraculix==FALSE){
              Zt <- array(0L,dim=c(maxZ,n.animals.gwas))
            }
            y_gwas <- array(0, dim=c(n.animals.gwas, population$info$bv.nr))

          }
          total_n <- sum(population$info$snp)
          x_mean <- x2_mean <- xy_mean <- numeric(total_n)
          first <- 1
          last <- min(maxZ, total_n)
          for(index3 in 1:ceiling(total_n/maxZ)){

            cindex <- 1
            for(index in 1:n.animals.gwas){
              k.database <- gwas.database[loop_elements_gwas[index,3],]
              if(miraculix){
                if (requireNamespace("miraculix", quietly = TRUE)) {
                  Zt[1:(last-first+1), cindex] <- miraculix::computeSNPS(population, k.database[1],k.database[2],kindex, from_p=first, to_p=last, what="geno", output_compressed=FALSE)
                }
              } else{
                Zt[1:(last-first+1), cindex] <- base::as.integer(colSums(compute.snps(population, k.database[1],k.database[2],kindex, import.position.calculation=import.position.calculation, from_p=first, to_p=last, decodeOriginsU=decodeOriginsU, bit.storing=bit.storing, nbits=nbits, output_compressed=FALSE)))
              }
              population$breeding[[k.database[1]]][[k.database[2]]][[kindex]][[16]] <- 1
              for(bven in 1:population$info$bv.nr){
                # HIER EVENTUELL SNPs auslesen!
                if(y.gwas.used=="pheno"){
                  y_gwas[cindex,bven] <- population$breeding[[k.database[1]]][[8+k.database[2]]][bven,kindex]
                } else if(y.gwas.used=="bv"){
                  y_gwas[cindex,bven] <- population$breeding[[k.database[1]]][[6+k.database[2]]][bven,kindex]
                } else if(y.gwas.used=="bve"){
                  y_gwas[cindex,bven] <- population$breeding[[k.database[1]]][[2+k.database[2]]][bven,kindex]
                }

              }
              cindex <- cindex +1
            }
          }
          if(gwas.group.standard){
            for(indexg in 1:(length(gwas_start)-1)){
              if(gwas_start[indexg]<=(gwas_start[indexg+1]-1)){
                y_gwas[gwas_start[indexg]:(gwas_start[indexg+1]-1), bven] <- y_gwas[gwas_start[indexg]:(gwas_start[indexg+1]-1), bven] - mean(y_gwas[gwas_start[indexg]:(gwas_start[indexg+1]-1), bven])
              }
            }
          }

          x_mean[first:last] <- rowMeans(Zt[1:(last-first+1),])
          x2_mean[first:last] <- rowMeans(Zt[1:(last-first+1),]^2)
          xy_mean[first:last] <- rowMeans(Zt[1:(last-first+1),]*y_gwas[,bven])
          first <- first + maxZ
          last <- min(maxZ*(index3+1), total_n)
        } else{
          if(nrow(gwas.database)!=nrow(bve.database) || prod(gwas.database==bve.database)==0){
            Zt <- array(0L,dim=c(sum(population$info$snp), n.animals.gwas))
            y_gwas <- array(0, dim=c(n.animals.gwas, population$info$bv.nr))
            cindex <- 1
            for(index in 1:n.animals.gwas){
              k.database <- gwas.database[loop_elements_gwas[index,3],]
              if(miraculix){
                if (requireNamespace("miraculix", quietly = TRUE)) {
                  Zt[,cindex] <- miraculix::computeSNPS(population, k.database[1],k.database[2],kindex, what="geno", output_compressed=FALSE)
                }
              } else{
                Zt[,cindex] <- base::as.integer(colSums(compute.snps(population, k.database[1],k.database[2],kindex, import.position.calculation=import.position.calculation, decodeOriginsU=decodeOriginsU, bit.storing=bit.storing, nbits=nbits, output_compressed=FALSE)))
              }
              for(bven in 1:population$info$bv.nr){
                # HIER EVENTUELL SNPs auslesen!
                if(y.gwas.used=="pheno"){
                  y_gwas[cindex,bven] <- population$breeding[[k.database[1]]][[8+k.database[2]]][bven,kindex]
                } else if(y.gwas.used=="bv"){
                  y_gwas[cindex,bven] <- population$breeding[[k.database[1]]][[6+k.database[2]]][bven,kindex]
                } else if(y.gwas.used=="bve"){
                  y_gwas[cindex,bven] <- population$breeding[[k.database[1]]][[2+k.database[2]]][bven,kindex]
                }
              }
              cindex <- cindex +1
            }
          }

          if(gwas.group.standard){

            for(indexg in 1:(length(gwas_start)-1)){
              if(gwas_start[indexg]>=(gwas_start[indexg+1]-1)){
                y_gwas[gwas_start[indexg]:(gwas_start[indexg+1]-1), bven] <- y_gwas[gwas_start[indexg]:(gwas_start[indexg+1]-1), bven] - mean(y_gwas[gwas_start[indexg]:(gwas_start[indexg+1]-1), bven])
              }
            }
          }
          x_mean <- rowMeans(Zt)
          x2_mean <- rowMeans(Zt^2)
#          xy_mean <- colMeans(Z*y_gwas[,bven])
          xy_mean <- colMeans(t(Zt)*y_gwas[,bven])
        }

        n <- length(y_gwas[,bven])
        y_mean <- mean(y_gwas[,bven])

        b1 <- (n *xy_mean - n *x_mean * y_mean) / (x2_mean*n - n *x_mean^2)


        if(approx.residuals==FALSE){
          sigma1 <- 1/(n * (x2_mean-(x_mean)^2))
          b0 <- y_mean - b1 * x_mean
          var1 <- numeric(length(sigma1))
          for(index in 1:length(sigma1)){
            var1[index] <- sigma1[index] * stats::var(y_gwas[,bven] - b1[index] * Zt[index,] - b0[index]) * (n-1)/(n-2)
          }

        } else{
          var1 <- 1/(n * (x2_mean-(x_mean)^2)) * stats::var(y_gwas[,bven])
        }
        test <- b1/sqrt(var1)
        gwas_hat <- cbind(gwas_hat, test, deparse.level = 0)
        #sorted <- sort(abs(test), index.return=TRUE)
      }

    }

    if(report.accuracy){
      if(verbose) cat("Correlation between genetic values and BVE:\n")
      if(n.rep==0){
        y_hat_temp <- y_hat
        y_hat_temp[y_hat_temp==0] <- NA
        acc <- suppressWarnings(stats::cor(y_real2[bve.insert,], y_hat_temp[bve.insert,], use="pairwise.complete.obs"))
      } else{
        insert.temp <- numeric(length(bve.insert.copy))

        if(length(stay.loop.elements)>0){
          for(index in (1:nrow(loop_elements_copy))[bve.insert.copy]){
            inserter <- which(stay.loop.elements==loop_elements_copy[index,6])
            insert.temp[index] <- if(length(inserter)==1){ inserter} else{NA}
          }
        } else{
          for(index in (1:nrow(loop_elements_copy))[bve.insert.copy]){
            insert.temp[index] <- loop_elements_copy[index,6]
          }

        }
        y_hat_temp <- rbind(y_hat[bve.insert,,drop=FALSE], y_hat[insert.temp,,drop=FALSE])
        y_hat_temp[y_hat_temp==0] <- NA
        acc <- suppressWarnings(stats::cor(rbind(y_real2[bve.insert,,drop=FALSE], y_real2[insert.temp,, drop=FALSE]),
                                           y_hat_temp, use="pairwise.complete.obs"))
      }
      if(length(acc)==1){
        acc <- matrix(acc,nrow=1)
      }

      if(sum(is.na(acc))>0){
        acc[is.na(acc)] <- 0
      }
      if(verbose) cat(diag(acc))
      if(verbose) cat("\n")
    }

  }
  if(u_hat_possible && bve && estimate.u && relationship.matrix=="vanRaden"){
    population$info$u_hat[[length(population$info$u_hat)+1]] <- u_hat
    population$info$u_hat_single[[length(population$info$u_hat)]] <- list()
    for(bven in 1:ncol(u_hat)){
      population$info$u_hat_single[[length(population$info$u_hat)]][[bven]] <- cbind((-2*p_i) *u_hat[,bven],(-2*p_i+1) *u_hat[,bven],(-2*p_i+2) *u_hat[,bven], deparse.level = 0)
    }
  } else if(u_hat_possible && bve && estimate.u && relationship.matrix=="CM"){
    population$info$u_hat[[length(population$info$u_hat)+1]] <- u_hat
    population$info$u_hat_single[[length(population$info$u_hat)]] <- list()
    for(bven in 1:ncol(u_hat)){
      population$info$u_hat_single[[length(population$info$u_hat)]][[bven]] <- cbind(u_hat[1:nrow(Zt),bven],u_hat[1:nrow(Zt)+ nrow(Zt),bven],u_hat[1:nrow(Zt)+2*nrow(Zt),bven], deparse.level = 0)
    }
  }
  if(gwas.u){
    if(sum(is.na(gwas_hat)>0)){
      gwas_hat[is.na(gwas_hat)] <- 0
    }
    population$info$gwas_hat[[length(population$info$gwas_hat)+1]] <- gwas_hat
  }

  if(store.comp.times.bve){
    comp.times.bve[5] <- as.numeric(Sys.time())
  }

  if(store.comp.times){
    comp.times[5] <- as.numeric(Sys.time())
  }
}

  #######################################################################
  ############## Selection of top individuals (sire/dams) ###############
  #######################################################################
{

  if(length(selection.m.database)>0 & selection.size[1]==0){
    selection.size[1] <- sum(selection.m.database[,4] - selection.m.database[,3] + 1)
    if(verbose) cat("No selection.size provided. Use all available selected individuals.\n")
    if(verbose) cat(paste0(selection.size[1], " male individuals selected.\n"))
  }
  if(length(selection.f.database)>0 & selection.size[2]==0){
    selection.size[2] <- sum(selection.f.database[,4] - selection.f.database[,3] + 1)
    if(verbose) cat("No selection.size provided. Use all available selected individuals.\n")
    if(verbose) cat(paste0(selection.size[2], " female individuals selected.\n"))
  }



  if(length(threshold.selection)==1 && is.na(threshold.selection)){
    threshold.selection <- NULL
  }

  if(length(threshold.selection)>0){
    if(threshold.sign=="<"){
      selection.highest <- c(FALSE, FALSE)
    }
  }

    {
    best <- list(matrix(0, nrow=selection.size[1],ncol=5),
                 matrix(0, nrow=selection.size[2],ncol=5))


    addsel <- c(2,2)
    if(selection.criteria[1]=="bv"){
      addsel[1] = 6
    }
    if(selection.criteria[2]=="bv"){
      addsel[2] = 6
    }
    if(selection.criteria[1]=="pheno"){
      addsel[1] = 8
    }
    if(selection.criteria[2]=="pheno"){
      addsel[2] = 8
    }

    if(selection.criteria[1]=="offpheno"){
      addsel[1] = 26
    }
    if(selection.criteria[2]=="offpheno"){
      addsel[2] = 26
    }

    if(add.class.cohorts && length(selection.m.cohorts)>0){
      for(index in 1:length(selection.m.cohorts)){
        take <- which(population$info$cohorts[,1]==selection.m.cohorts[index])
        class.m <- c(class.m, as.numeric(population$info$cohorts[take,5]))
      }
      class[[1]] <- class.m
    }
    if(add.class.cohorts && length(selection.f.cohorts)>0){
      for(index in 1:length(selection.f.cohorts)){
        take <- which(population$info$cohorts[,1]==selection.f.cohorts[index])
        class.f <- c(class.f, as.numeric(population$info$cohorts[take,5]))
      }
      class[[2]] <- class.f
    }




  }

  activ_groups_list <- list(selection.m.database, selection.f.database)
  chosen.animals.list <- list()
  selection.sex <- c(selection.m, selection.f)
  selection.miesenberger <- c(selection.m.miesenberger, selection.f.miesenberger)
  multiple.bve.weights <- list(multiple.bve.weights.m, multiple.bve.weights.f)
  multiple.bve.scale <- c(multiple.bve.scale.m, multiple.bve.scale.f)

  multiple.bve.scale[multiple.bve.scale=="bve"] <- "bve_sd"
  multiple.bve.scale[multiple.bve.scale=="pheno"] <- "pheno_sd"
  multiple.bve.scale[multiple.bve.scale=="bv"] <- "bv_sd"


  sd_scaling <- rep(1, population$info$bv.nr)
  if(length(fixed.breeding)==0 || length(fixed.breeding.best)>0){
    if(sum(selection.size)>0){
      if(verbose) cat("Start selection procedure.\n")
      for(sex in (1:2)[selection.size>0]){
        possible_animals <- NULL
        activ_groups <- activ_groups_list[[sex]]
        for(index5 in 1:nrow(activ_groups)){
          possible_animals <- rbind(possible_animals, cbind(activ_groups[index5,1], activ_groups[index5,2], activ_groups[index5,3]:activ_groups[index5,4], population$breeding[[activ_groups[index5,1]]][[4+activ_groups[index5,2]]][activ_groups[index5,3]:activ_groups[index5,4]]))
        }

        if(length(reduced.selection.panel[[sex]])>0 && length(possible_animals)>0){
          possible_animals <- possible_animals[reduced.selection.panel[[sex]],,drop=FALSE]
        }
        relevant.animals <- rep(possible_animals[,4], length(class[[sex]])) == rep(class[[sex]],each=nrow(possible_animals))
        relevant.animals <- unique(c(0,relevant.animals * 1:nrow(possible_animals)))[-1]
        possible_animals <- rbind(NULL, possible_animals[unique(relevant.animals),])
        n.animals <- nrow(possible_animals)


        if(n.animals==0){
          stop("No available individuals for selection provided - check gen/database/cohorts and classes of your input!")
        }else if(n.animals<selection.size[sex]){
          warnings(paste0("Less individuals available than to select. Reduce number of selected individuals to ", n.animals))
          selection.size[sex] <- n.animals
          best[[sex]] <- matrix(0, nrow=selection.size[sex],ncol=5)
        }

        if(selection.sex[sex]=="random"){
          if(population$info$bv.nr==0){
            import.bv <- rep(0, n.animals)
          } else if(population$info$bv.nr==1){
            import.bv <- rep(0, n.animals)
            for(index5 in 1:nrow(possible_animals)){
              import.bv[index5] <- population$breeding[[possible_animals[index5,1]]][[possible_animals[index5,2]+addsel[possible_animals[index5,2]]]][,possible_animals[index5,3]]
            }
          } else{
            if(multiple.bve=="add"){
              breeding.values <- matrix(0, ncol=nrow(possible_animals), nrow= population$info$bv.nr)
              for(index5 in 1:nrow(possible_animals)){
                breeding.values[,index5] <- population$breeding[[possible_animals[index5,1]]][[possible_animals[index5,2]+addsel[possible_animals[index5,2]]]][,possible_animals[index5,3]]
              }


              if(multiple.bve.scale[sex]=="bve_sd" || multiple.bve.scale[sex]=="pheno_sd" || multiple.bve.scale[sex]=="bv_sd"){

                sd_store <- numeric(population$info$bv.nr)
                if(ncol(breeding.values)>1){
                  P <- stats::var(t(breeding.values))
                } else{
                  P <- diag(nrow(breeding.values))
                }

                if(multiple.bve.scale[sex]=="pheno_sd" || multiple.bve.scale[sex] =="bv_sd"){
                  pheno.values <- genomic.values <-  matrix(0, ncol=nrow(possible_animals), nrow= population$info$bv.nr)
                  for(index5 in 1:nrow(possible_animals)){
                    pheno.values[,index5] <- population$breeding[[possible_animals[index5,1]]][[possible_animals[index5,2]+8]][,possible_animals[index5,3]]
                    genomic.values[,index5] <- population$breeding[[possible_animals[index5,1]]][[possible_animals[index5,2]+6]][,possible_animals[index5,3]]
                  }
                }

                for(bven in 1:population$info$bv.nr){
                  breeding.values[bven,] <- breeding.values[bven,] - mean(breeding.values[bven,])
                  if(ncol(breeding.values)!=1){
                    if(multiple.bve.scale[sex]=="bve_sd"){
                      sd_store[bven] <- stats::sd(breeding.values[bven,])
                    } else if(multiple.bve.scale[sex]=="pheno_sd"){
                      sd_store[bven] <- stats::sd(pheno.values[bven,])
                      if(is.na(sd_store[bven])){
                        sd_store[bven] <- 0
                      }
                      if(sd_store[bven]==0){
                        if(verbose) cat("No observed phenotypes in the group of selected individuals. Sure you want to scale according to phenotypes?\n")
                        if(verbose) cat("Use residual variance for scaling!\n")
                        sd_store[bven] <- sqrt(sigma.e[bven])
                        if(is.na(sd_store[bven])){
                          sd_store[bven] <- 0
                        }
                      }
                    } else {
                      sd_store[bven] <- stats::sd(genomic.values[bven,])
                      if(sd_store[bven]==0){
                        if(verbose & multiple.bve.weights[[sex]][bven]!=0) cat("No variation in true genomic values. Sure you want to scale according to true genomic values?\n")
                        if(verbose & multiple.bve.weights[[sex]][bven]!=0) cat("Use residual variance for scaling!\n")
                        sd_store[bven] <- sqrt(sigma.e[bven])
                        if(is.na(sd_store[bven])){
                          sd_store[bven] <- 0
                        }
                      }
                    }

                  } else{
                    sd_store[bven] <- 1
                  }

                  if(sd_store[bven]>0){
                    breeding.values[bven,] <- breeding.values[bven,] / sd_store[bven]
                  }
                }
              }

              bve.sum <- colSums(breeding.values * multiple.bve.weights[[sex]])

            } else if(multiple.bve=="ranking"){
              breeding.values <- matrix(0, ncol=nrow(possible_animals), nrow= population$info$bv.nr)
              for(index5 in 1:nrow(possible_animals)){
                breeding.values[,index5] <- population$breeding[[possible_animals[index5,1]]][[possible_animals[index5,2]+addsel[possible_animals[index5,2]]]][,possible_animals[index5,3]]
              }
              ranking <- matrix(0, ncol=nrow(possible_animals), nrow= population$info$bv.nr)
              for(bven in 1:population$info$bv.nr){
                order <- sort(breeding.values[bven,], index.return=TRUE, decreasing=selection.highest[sex])$ix
                ranking[bven,order] <- length(order):1
              }
              bve.sum <- colSums(ranking*multiple.bve.weights[[sex]])
            }
            import.bv <- bve.sum
          }

          if(nrow(possible_animals)>0 & sum(abs(import.bv))>0){
            for(entry in 1:nrow(possible_animals)){
              population$breeding[[possible_animals[entry,1]]][[possible_animals[entry,2]+30]][possible_animals[entry,3]] <- import.bv[entry]
            }
          }


          chosen.animals <- sample(1:n.animals, selection.size[sex])
          best[[sex]][,1:3] <- possible_animals[chosen.animals,1:3]
          best[[sex]][,4] <- import.bv[chosen.animals]
        } else if(selection.sex[sex]=="function"){
          if(population$info$bv.nr==1){
            import.bv <- numeric(nrow(possible_animals))
            for(index5 in 1:nrow(possible_animals)){
              import.bv[index5] <- population$breeding[[possible_animals[index5,1]]][[possible_animals[index5,2]+addsel[possible_animals[index5,2]]]][,possible_animals[index5,3]]
            }
            import.bv[is.na(import.bv)] <- -Inf
            chosen.animals <- sort(import.bv, index.return=TRUE, decreasing=selection.highest[sex])$ix[1:sum(selection.size[[sex]])] # Diese 3er werden zu 4 in Weiblich
            if(sum(import.bv[chosen.animals[length(chosen.animals)]]==import.bv)>1){
              cutoff1 <- import.bv[chosen.animals[length(chosen.animals)]]
              chosen.animals <- NULL
              if(selection.highest[sex]){
                chosen.animals <- which(import.bv>cutoff1)
              } else{
                chosen.animals <- which(import.bv<cutoff1)
              }
              chosen.animals <- c(chosen.animals, sample(which(import.bv==cutoff1),  sum(selection.size[[sex]]) - length(chosen.animals)))
            }
            best[[sex]][,1:4] <- cbind(possible_animals[chosen.animals, 1], possible_animals[chosen.animals, 2], possible_animals[chosen.animals, 3], import.bv[chosen.animals])
          } else{
            if(multiple.bve=="add"){
              breeding.values <- matrix(0, ncol=nrow(possible_animals), nrow= population$info$bv.nr)
              for(index5 in 1:nrow(possible_animals)){
                breeding.values[,index5] <- population$breeding[[possible_animals[index5,1]]][[possible_animals[index5,2]+addsel[possible_animals[index5,2]]]][,possible_animals[index5,3]]
              }
              breeding.values[is.na(breeding.values)] <- -Inf
              if((multiple.bve.scale[sex]=="bve_sd" || multiple.bve.scale[sex]=="pheno_sd" || multiple.bve.scale[sex] =="bv_sd") && selection.miesenberger[sex]==FALSE){

                sd_scaling <- numeric(population$info$bv.nr)
                if(ncol(breeding.values)>1){
                  P <- stats::var(t(breeding.values), use="pairwise.complete.obs")
                } else{
                  P <- diag(nrow(breeding.values))
                }

                if(multiple.bve.scale[sex]=="pheno_sd" || multiple.bve.scale[sex] =="bv_sd"){
                  pheno.values <- genomic.values <-  matrix(0, ncol=nrow(possible_animals), nrow= population$info$bv.nr)
                  for(index5 in 1:nrow(possible_animals)){
                    pheno.values[,index5] <- population$breeding[[possible_animals[index5,1]]][[possible_animals[index5,2]+8]][,possible_animals[index5,3]]
                    genomic.values[,index5] <- population$breeding[[possible_animals[index5,1]]][[possible_animals[index5,2]+6]][,possible_animals[index5,3]]
                  }
                }

                if(sum(abs(breeding.values))==0){
                  if(verbose) cat(paste0("No selection criteria or breeding values for selection provided. Continue with random selection (", if(sex==1){"male side"}else{"female side"} ,").\n"))
                  sd_scaling <- rep(1, population$info$bv.nr)
                } else{
                  for(bven in 1:population$info$bv.nr){
                    breeding.values[bven,] <- breeding.values[bven,] - mean(breeding.values[bven,])
                    if(ncol(breeding.values)!=1){
                      if(multiple.bve.scale[sex]=="bve_sd"){
                        sd_scaling[bven] <- stats::sd(breeding.values[bven,], na.rm=TRUE)
                        if(sd_scaling[bven]==0 || is.na(sd_scaling[bven])){
                          if(verbose & multiple.bve.weights[[sex]][bven]!=0) cat(paste0("No estimated breeding values entered in the group of selected individuals (Trait: ", population$info$trait.name[bven], " - ", if(sex==1){"male side"}else{"female side"}, "). Please check your input variables!?\n"))
                          if(verbose & multiple.bve.weights[[sex]][bven]!=0) cat("Use residual variance!\n")
                          sd_scaling[bven] <- sqrt(sigma.e[bven])
                        }
                      } else if(multiple.bve.scale[sex]=="pheno_sd") {
                        sd_scaling[bven] <- stats::sd(pheno.values[bven,], na.rm=TRUE)

                        if(is.na(sd_scaling[bven])){
                          sd_scaling[bven] <- 0
                        }
                        if(sd_scaling[bven]==0){
                          if(verbose & multiple.bve.weights[[sex]][bven]!=0) cat(paste0("No observed phenotypes in the group of selected individuals (Trait: ", population$info$trait.name[bven], " - ", if(sex==1){"male side"}else{"female side"},"). Sure you want to scale according to phenotypes?\n"))
                          if(verbose & multiple.bve.weights[[sex]][bven]!=0) cat("Expected phenotypic sd based on one observation was used!\n")
                          sd_scaling[bven] <- sqrt(stats::var(genomic.values[bven,]) + sigma.e[bven])

                          if(is.na(sd_scaling[bven])){
                            sd_scaling[bven] <- 0
                          }

                        }
                      } else{
                        sd_scaling[bven] <- stats::sd(genomic.values[bven,], na.rm=TRUE)
                        if(sd_scaling[bven]==0){
                          if(verbose & multiple.bve.weights[[sex]][bven]!=0) cat(paste0("No variation in true genomic values in the group of selected individuals (Trait: ", population$info$trait.name[bven], " - ", if(sex==1){"male side"}else{"female side"},"). Sure you want to scale according to true genomic values?\n"))
                          if(verbose & multiple.bve.weights[[sex]][bven]!=0) cat("Expected phenotypic sd based on one observation was used!\n")
                          sd_scaling[bven] <- sqrt(stats::var(genomic.values[bven,]) + sigma.e[bven])
                          if(is.na(sd_scaling[bven])){
                            sd_scaling[bven] <- 0
                          }
                        }
                      }

                    } else{
                      sd_scaling[bven] <- 1
                    }

                    if(sd_scaling[bven]>0){
                      breeding.values[bven,] <- breeding.values[bven,] / sd_scaling[bven]
                    }
                  }
                }

              }

              if(selection.miesenberger[sex]){
                genomic.values <- matrix(0, ncol=nrow(possible_animals), nrow= population$info$bv.nr)
                bve.values <- matrix(0, ncol=nrow(possible_animals), nrow= population$info$bv.nr)
                pheno.values <- matrix(0, ncol=nrow(possible_animals), nrow= population$info$bv.nr)
                reliability.values <- matrix(0, ncol=nrow(possible_animals), nrow= population$info$bv.nr)
                bve.index <- matrix(0, ncol=nrow(possible_animals), nrow= population$info$bv.nr)

                for(index5 in 1:nrow(possible_animals)){
                  genomic.values[,index5] <- population$breeding[[possible_animals[index5,1]]][[possible_animals[index5,2]+6]][,possible_animals[index5,3]]
                  bve.values[,index5] <- population$breeding[[possible_animals[index5,1]]][[possible_animals[index5,2]+2]][,possible_animals[index5,3]]
                  pheno.values[,index5] <- population$breeding[[possible_animals[index5,1]]][[possible_animals[index5,2]+8]][,possible_animals[index5,3]]
                  reliability.values[,index5] <- population$breeding[[possible_animals[index5,1]]][[possible_animals[index5,2]+18]][,possible_animals[index5,3]]
                }
                if(sum(bve.values==0)==length(bve.values)){
                  if(verbose) cat("No breeding values for selected individuals! Miesenberger selection was skip.\n")

                  index.weights <- multiple.bve.weights[[sex]]
                  for(index5 in 1:nrow(possible_animals)){
                    bve.index[,index5] <- index.weights
                    population$breeding[[possible_animals[index5,1]]][[possible_animals[index5,2]+20]][,possible_animals[index5,3]] <- bve.index[,index5]
                  }
                } else{
                  active_traits <- which(rowSums(bve.values==0)<nrow(bve.values))
                  if(length(active_traits)<population$info$bv.nr){
                    if(verbose) cat("Only use traits", population$info$trait.name[active_traits], "for Miesenberger.\n")
                    if(verbose) cat("Traits", population$info$trait.name[-active_traits], "with no BVE.\n")
                  }
                  genomic.cov <- stats::var(t(genomic.values[active_traits,,drop=FALSE]))
                  bve.cov <- stats::var(t(bve.values[active_traits,,drop=FALSE]))
                  pheno.cov <- stats::var(t(pheno.values[active_traits,colSums(is.na(pheno.values))==0,drop=FALSE]))

#                  diag(genomic.cov) <- diag(genomic.cov) + 10^-8
#                  diag(bve.cov) <- diag(bve.cov) + 10^-8
#                  diag(pheno.cov) <- diag(pheno.cov) + 10^-8

                  heri.emp <- diag(genomic.cov) / diag(pheno.cov)
                  derived.reliability <- diag(stats::cor(t(genomic.values[active_traits,,drop=FALSE]), t(bve.values[active_traits,,drop=FALSE])))^2
                  estimated.reliability <- sqrt(diag(bve.cov) / diag(pheno.cov))
                  estimated.reliability[diag(pheno.cov)==0] <- derived.reliability[diag(pheno.cov)==0]
                  estimated.reliability[estimated.reliability>1] <- 1


                  for(bven in 1:length(active_traits)){
                    trait_index <- active_traits[bven]
                    if(selection.miesenberger.reliability.est=="estimated"){
                      reliability.values[trait_index, reliability.values[trait_index,]==0] <- estimated.reliability[bven]
                    } else if(selection.miesenberger.reliability.est=="heritability"){
                      reliability.values[trait_index, reliability.values[trait_index,]==0] <- heri.emp[bven]
                    } else if(selection.miesenberger.reliability.est=="derived"){
                      reliability.values[trait_index, reliability.values[trait_index,]==0] <- derived.reliability[bven]
                    }

                  }
                  reliability.values[is.na(reliability.values)] <- 0

                  V <- bve.cov
                  V1 <- MASS::ginv(V)

                  # V1 <- chol2inv(chol(V))

                  G_cov <- genomic.cov
                  RG <- sqrt(diag(1/diag(G_cov))) %*% G_cov %*% sqrt(diag(1/diag(G_cov)))

                  if(multiple.bve.scale[sex]=="pheno_sd"){
                    index.weights <- multiple.bve.weights[[sex]][active_traits] / sqrt(diag(pheno.cov))
                    if(sum(is.na(diag(pheno.cov))) >0 || sum(diag(pheno.cov)==0)>0){
                      if(verbose) cat("Are you sure you want to scale Miesenberger w by phenotypic sd? No phenotypes for some traits!\n")
                      if(verbose) cat("Expected phenotypic sd was used!\n")
                      index.weights <- multiple.bve.weights[[sex]][active_traits] / sqrt(diag(genomic.cov) + sigma.e[active_traits])
                    }
                  } else if(multiple.bve.scale[sex]=="unit"){
                    index.weights <- multiple.bve.weights[[sex]][active_traits]
                  } else if(multiple.bve.scale[sex]=="bve_sd"){
                    index.weights <- multiple.bve.weights[[sex]][active_traits] / sqrt(diag(bve.cov))
                    index.weights[diag(bve.cov)==0] <- 0
                  } else if(multiple.bve.scale[sex]=="bv_sd"){
                    index.weights <- multiple.bve.weights[[sex]][active_traits] / sqrt(diag(genomic.cov))
                    index.weights[diag(genomic.cov)==0] <- 0
                  }
                  if(miesenberger.trafo>0){

                    eigenV <- eigen(V)
                    keeps <- which(eigenV$values > (miesenberger.trafo*sum(eigenV$values)))

                    P <- eigenV$vectors[,keeps]

                    V_trafo <- t(P) %*% V %*% P
                    G_trafo <- t(P) %*% G_cov %*% P

                    V1_trafo <- MASS::ginv(V_trafo)
                    RG_trafo <- sqrt(diag(1/diag(G_trafo))) %*% G_trafo %*% sqrt(diag(1/diag(G_trafo)))

                    w_trafo <- t(P) %*% index.weights

                    breeding.values1 <- t(P) %*% breeding.values[active_traits, ]
                    bve.index1 <- matrix(0, nrow=nrow(breeding.values1), ncol=ncol(breeding.values1))
                    for(index5 in 1:nrow(possible_animals)){
                      r_trafo <- sqrt(colSums(reliability.values[active_traits,index5] * (eigenV$vectors[,keeps]^2)))
                      bve.index1[,index5] <- miesenberger.index(V1 = V1_trafo, V = V_trafo, G_trafo, RG = RG_trafo, r = r_trafo, w = w_trafo)
                    }


                  } else{
                    for(index5 in 1:nrow(possible_animals)){
                      bve.index[active_traits,index5] <- miesenberger.index(V1=V1, V= V, G = G_cov, RG = RG, r = sqrt(reliability.values[active_traits,index5]), w = index.weights)
                      population$breeding[[possible_animals[index5,1]]][[possible_animals[index5,2]+20]][,possible_animals[index5,3]] <- bve.index[,index5]
                    }

                  }

                  if(verbose) cat("Average index used for in selection according to Miesenberger:\n")

                  if(miesenberger.trafo>0){
                    bve.index2 <-  P %*% bve.index1
                  } else{
                    bve.index2 <- bve.index
                  }
                  if(multiple.bve.scale[sex]=="pheno_sd"){
                    if(sum(is.na(diag(pheno.cov)))>0 || sum(diag(pheno.cov)==0)>0){
                      if(verbose) cat(paste0(rowMeans(bve.index2) * sqrt(diag(genomic.cov) + sigma.e[active_traits]), " avg. index weightings.\n"))
                    } else{
                      if(verbose) cat(paste0(rowMeans(bve.index2) * sqrt(diag(pheno.cov)), " avg. index weightings.\n"))
                    }
                  } else if(multiple.bve.scale[sex]=="unit"){
                    if(verbose) cat(paste0(rowMeans(bve.index2), " avg. index weightings.\n"))
                  } else if(multiple.bve.scale[sex]=="bve_sd"){
                    if(verbose) cat(paste0(rowMeans(bve.index2) * sqrt(diag(bve.cov)), " avg. index weightings.\n"))
                  } else if(multiple.bve.scale[sex]=="bv_sd"){
                    if(verbose) cat(paste0(rowMeans(bve.index2) * sqrt(diag(genomic.cov)), " avg. index weightings.\n"))
                  }


                }

                if(miesenberger.trafo>0){
                  bve.sum <- colSums(bve.index1 * (breeding.values1-rowMeans(breeding.values1)), na.rm=TRUE)
                } else{
                  bve.sum <- colSums(bve.index * (breeding.values-rowMeans(breeding.values)), na.rm=TRUE)
                }


              } else{
                bve.sum <- colSums(breeding.values * multiple.bve.weights[[sex]], na.rm=TRUE)
              }
            } else if(multiple.bve=="ranking"){
              breeding.values <- matrix(0, ncol=nrow(possible_animals), nrow= population$info$bv.nr)
              for(index5 in 1:nrow(possible_animals)){
                breeding.values[,index5] <- population$breeding[[possible_animals[index5,1]]][[possible_animals[index5,2]+addsel[possible_animals[index5,2]]]][,possible_animals[index5,3]]
              }
              breeding.values[is.na(breeding.values)] <- -Inf
              ranking <- matrix(0, ncol=nrow(possible_animals), nrow= population$info$bv.nr)
              for(bven in 1:population$info$bv.nr){
                order <- sort(breeding.values[bven,], index.return=TRUE, decreasing=selection.highest[sex])$ix
                ranking[bven,order] <- length(order):1
              }
              bve.sum <- colSums(ranking*multiple.bve.weights[[sex]], na.rm=TRUE)
            }

            if(nrow(possible_animals)>0 & sum(abs(bve.sum))>0){
              for(entry in 1:nrow(possible_animals)){
                population$breeding[[possible_animals[entry,1]]][[possible_animals[entry,2]+30]][possible_animals[entry,3]] <- bve.sum[entry]
              }
            }
            chosen.animals <- sort(bve.sum, index.return=TRUE, decreasing=selection.highest[sex])$ix[1:sum(selection.size[sex])] # Diese 3er werden zu 4 in Weiblich

            if(sum(bve.sum[chosen.animals[length(chosen.animals)]]==bve.sum)>1){
              cutoff1 <- bve.sum[chosen.animals[length(chosen.animals)]]
              chosen.animals <- NULL
              if(selection.highest[sex]){
                chosen.animals <- which(bve.sum>cutoff1)
              } else{
                chosen.animals <- which(bve.sum<cutoff1)
              }
              chosen.animals <- c(chosen.animals, sample(which(bve.sum==cutoff1),  sum(selection.size[[sex]]) - length(chosen.animals)))
            }
            best[[sex]][,1:4] <- cbind(possible_animals[chosen.animals,1:3, drop=FALSE], bve.sum[chosen.animals])
          }

        }

        if(selection.miesenberger[sex]==FALSE){
          for(index5 in 1:nrow(possible_animals)){
            population$breeding[[possible_animals[index5,1]]][[possible_animals[index5,2]+20]][,possible_animals[index5,3]] <- multiple.bve.weights[[sex]]
          }
        }

        ## Just some follow up computation in case individuasl are not picked with equal share
        if(population$info$bv.nr==1){
          breeding.values <- numeric(length(chosen.animals))
          for(running in 1:length(chosen.animals)){
            if(best.selection.criteria[[sex]]=="bv"){
              breeding.values[running] <- population$breeding[[possible_animals[chosen.animals[running],1]]][[possible_animals[chosen.animals[running],2]+6]][1,possible_animals[chosen.animals[running],3]]
            } else if(best.selection.criteria[[sex]]=="pheno"){
              breeding.values[running] <- population$breeding[[possible_animals[chosen.animals[running],1]]][[possible_animals[chosen.animals[running],2]+8]][1,possible_animals[chosen.animals[running],3]]
            } else{
              breeding.values[running] <- population$breeding[[possible_animals[chosen.animals[running],1]]][[possible_animals[chosen.animals[running],2]+2]][1,possible_animals[chosen.animals[running],3]]
            }
          }
          breeding.values[is.na(breeding.values)] <- -Inf
          best[[sex]][,5] <- breeding.values
        } else if(population$info$bv.nr>1){
          if(multiple.bve=="add"){
            breeding.values <- matrix(0, nrow=population$info$bv.nr, ncol=length(chosen.animals))
            breeding.values.scaling <- matrix(0, nrow=population$info$bv.nr, ncol=length(chosen.animals))
            for(running in 1:length(chosen.animals)){
              if(best.selection.criteria[[sex]]=="bv"){
                breeding.values[,running] <- population$breeding[[possible_animals[chosen.animals[running],1]]][[possible_animals[chosen.animals[running],2]+6]][,possible_animals[chosen.animals[running],3]]
              } else if(best.selection.criteria[[sex]]=="pheno"){
                breeding.values[,running] <- population$breeding[[possible_animals[chosen.animals[running],1]]][[possible_animals[chosen.animals[running],2]+8]][,possible_animals[chosen.animals[running],3]]
              } else{
                breeding.values[,running] <- population$breeding[[possible_animals[chosen.animals[running],1]]][[possible_animals[chosen.animals[running],2]+2]][,possible_animals[chosen.animals[running],3]]
              }
              breeding.values.scaling[,running] <- population$breeding[[possible_animals[chosen.animals[running],1]]][[possible_animals[chosen.animals[running],2]+20]][,possible_animals[chosen.animals[running],3]]

            }
            breeding.values[is.na(breeding.values)] <- -Inf

            if((multiple.bve.scale[sex]=="bve_sd" || multiple.bve.scale[sex]=="pheno_sd" || multiple.bve.scale[sex]=="bv_sd") && selection.miesenberger[sex]==FALSE){

              sd_scaling <- numeric(population$info$bv.nr)
              if(ncol(breeding.values)>1){
                P <- stats::var(t(breeding.values))
              } else{
                P <- diag(nrow(breeding.values))
              }

              if(multiple.bve.scale[sex]=="pheno_sd" || multiple.bve.scale[sex]=="bv_sd"){
                pheno.values <- genomic.values <-  matrix(0, ncol=nrow(possible_animals), nrow= population$info$bv.nr)
                for(index5 in 1:nrow(possible_animals)){
                  pheno.values[,index5] <- population$breeding[[possible_animals[index5,1]]][[possible_animals[index5,2]+8]][,possible_animals[index5,3]]
                  genomic.values[,index5] <- population$breeding[[possible_animals[index5,1]]][[possible_animals[index5,2]+6]][,possible_animals[index5,3]]
                }
              }
              for(bven in 1:population$info$bv.nr){
                breeding.values[bven,] <- breeding.values[bven,] - mean(breeding.values[bven,])
                if(ncol(breeding.values)!=1){
                  if(multiple.bve.scale[sex]=="bve_sd"){
                    sd_scaling[bven] <- stats::sd(breeding.values[bven,])
                    if(sd_scaling[bven]==0 || is.na(sd_scaling[bven])){
                      if(verbose & multiple.bve.weights[[sex]][bven]!=0 & selection.sex[sex]!= "random") cat("No estimated breeding values entered in the group of selected individuals. Please check your input variables!?\n")
                      if(verbose & multiple.bve.weights[[sex]][bven]!=0 & selection.sex[sex]!= "random") cat("Use residual variance!\n")
                      sd_scaling[bven] <- sqrt(sigma.e[bven])
                    }
                  } else if(multiple.bve.scale[sex] =="pheno_sd"){
                    sd_scaling[bven] <- stats::sd(pheno.values[bven,])
                    if(is.na(sd_scaling[bven]) || sd_scaling[bven]==0){
                      if(verbose & multiple.bve.weights[[sex]][bven]!=0 & selection.sex[sex]!= "random") cat("No observed phenotypes in the group of selected individuals. Sure you want to scale according to phenotypes?\n")
                      if(verbose & multiple.bve.weights[[sex]][bven]!=0 & selection.sex[sex]!= "random") cat("Expected phenotypic sd based on one observation was used!\n")
                      sd_scaling[bven] <- sqrt(stats::var(genomic.values[bven,]) + sigma.e[bven])
                    }
                  } else{
                    sd_scaling[bven] <- stats::sd(genomic.values[bven,])
                    if(sd_scaling[bven]==0){
                      if(verbose & multiple.bve.weights[[sex]][bven]!=0 & selection.sex[sex]!= "random") cat("No variation in true genomic values in the group of selected individuals. Sure you want to scale according to true genomic values?\n")
                      if(verbose & multiple.bve.weights[[sex]][bven]!=0 & selection.sex[sex]!= "random") cat("Expected phenotypic sd based on one observation was used!\n")
                      sd_scaling[bven] <- sqrt(stats::var(genomic.values[bven,]) + sigma.e[bven])
                      if(is.na(sd_scaling[bven])){
                        sd_scaling[bven] <- 0
                      }
                    }
                  }

                } else{
                  sd_scaling[bven] <- 1
                }

                if(sd_scaling[bven]>0){
                  breeding.values[bven,] <- breeding.values[bven,] / sd_scaling[bven]
                }
              }
            }
            bve.sum <- colSums(breeding.values * breeding.values.scaling, na.rm=TRUE)
          } else if(multiple.bve=="ranking"){
            ranking <- matrix(0, nrow=population$info$bv.nr, ncol=length(chosen.animals))
            for(running in 1:length(chosen.animals)){
              if(best.selection.criteria[[sex]]=="bv"){
                ranking[,running] <- population$breeding[[possible_animals[chosen.animals[running],1]]][[possible_animals[chosen.animals[running],2]+6]][,possible_animals[chosen.animals[running],3]]
              } else if(best.selection.criteria[[sex]]=="pheno"){
                ranking[,running] <- population$breeding[[possible_animals[chosen.animals[running],1]]][[possible_animals[chosen.animals[running],2]+8]][,possible_animals[chosen.animals[running],3]]
              } else{
                ranking[,running] <- population$breeding[[possible_animals[chosen.animals[running],1]]][[possible_animals[chosen.animals[running],2]+2]][,possible_animals[chosen.animals[running],3]]
              }
            }
            ranking[is.na(ranking)] <- -Inf
            for(bven in 1:population$info$bv.nr){
              order <- sort(ranking[bven,], index.return=TRUE, decreasing=selection.highest[sex])$ix
              ranking[bven,order] <- length(order):1
            }
            bve.sum <- colSums(ranking*breeding.values.scaling, na.rm=TRUE)
          }
          best[[sex]][,5] <- bve.sum
        }

        reorder <- sort(best[[sex]][,4], index.return=TRUE, decreasing = selection.highest[sex])$ix
        best[[sex]] <- best[[sex]][reorder,,drop=FALSE]
        chosen.animals <- chosen.animals[reorder]
        if(ignore.best[sex]>0){
          best[[sex]] <- best[[sex]][-(1:ignore.best[sex]),,drop=FALSE]
          chosen.animals <- chosen.animals[-(1:ignore.best[sex])]
          selection.size[sex] <- selection.size[sex] - ignore.best[sex]

        }

        chosen.animals.list[[sex]] <- chosen.animals

        if(!best.selection.manual.reorder){
          new_order <- numeric(nrow(best[[sex]]))
          temp1 <- t(best[[sex]][,1:3])
          for(index in 1:nrow(best[[sex]])){
            new_order[index] <- which(colSums(temp1==possible_animals[index,1:3])==3)
          }
          best.selection.manual.ratio[[sex]][new_order] <- best.selection.manual.ratio[[sex]]
        }
      }

      if(fixed.assignment!=FALSE){
        nm <- nrow(best[[1]])
        nw <- nrow(best[[2]])
        male.rounds <- breeding.size.total / nm
        female.rounds <- breeding.size.total / nw
        full.m <- nm*(male.rounds - round(male.rounds - 0.5))
        full.f <- nw*(female.rounds - round(female.rounds -0.5))
        fixed.assignment.m <- c(rep(1:full.m,each=male.rounds-full.m/nm+1), rep((full.m+1):nm, each =male.rounds-full.m/nm))
        fixed.assignment.f <- c(rep(1:full.f,each=female.rounds - full.f/nw +1), rep((full.f+1):nw, each =female.rounds-full.f/nw))
        fixed.assignment.m <- fixed.assignment.m[1:breeding.size.total]
        fixed.assignment.f <- fixed.assignment.f[1:breeding.size.total]
        if(fixed.assignment=="bestworst" || fixed.assignment=="worstbest"){
          fixed.assignment.f <- fixed.assignment.f[breeding.size.total:1]
        }
        fixed.breeding <- cbind(best[[1]][fixed.assignment.m,1:3], best[[2]][fixed.assignment.f,1:3], sex.animal)
      }
      if(ogc){
        animallist <- rbind(cbind(best[[1]],1), cbind(best[[2]],2))
        n.animals <- nrow(animallist)

        if(sum(abs(animallist[,4]))==0){
          if(verbose) cat("No breeding values stored for OGC. Use breeding value of first trait as u!\n")
          for(index in 1:nrow(animallist)){
            animallist[index,4] <- population$breeding[[animallist[index,1]]][[2+animallist[index,2]]][1,animallist[index,3]]
          }
        }

        if(sum(abs(animallist[,4]))==0){
          if(verbose) cat("No breeding value estimated available! Use genomic value of first trait as u!\n")
          for(index in 1:nrow(animallist)){
            animallist[index,4] <- population$breeding[[animallist[index,1]]][[6+animallist[index,2]]][1,animallist[index,3]]
          }
        }

        BV <- animallist[,4]
        Sex <- rep("male", length(BV))
        Sex[animallist[,6]==2] <- "female"

        if(relationship.matrix.ogc != "kinship" || relationship.matrix.ogc != "pedigree"){
          if(miraculix){
            Z.code <- miraculix::computeSNPS(population, animallist[,1], animallist[,2], animallist[,3], what="geno", output_compressed = TRUE)
          } else{
            Zt <- array(0,dim=c(sum(population$info$snp), n.animals))
            for(index in 1:n.animals){
              Zt[,index] <- colSums(compute.snps(population, animallist[index,1], animallist[index,2], animallist[index,3], import.position.calculation=import.position.calculation, decodeOriginsU=decodeOriginsU, bit.storing=bit.storing, nbits=nbits, output_compressed=FALSE))
            }
          }
        }


        # Verwandtschaftsmatrix:
        if(relationship.matrix.ogc=="kinship" || relationship.matrix.ogc == "pedigree"){
          A <- kinship.exp(population, database = animallist[,c(1,2,3,3)], depth.pedigree = depth.pedigree.ogc, verbose=verbose)
        } else if(relationship.matrix.ogc=="vanRaden"){
          if(miraculix){
            p_i <- miraculix::allele_freq(Z.code) # Noch nicht implementiert?
            A <- miraculix::relationshipMatrix(Z.code, centered=TRUE, normalized=TRUE)

          } else if(miraculix.mult){
            p_i <- rowSums(Zt)/ncol(Zt)/2
            Zt_miraculix <- miraculix::genomicmatrix(Zt)
            A <- miraculix::relationshipMatrix(Zt_miraculix, centered=TRUE, normalized=TRUE)
          } else{
            p_i <- rowSums(Zt)/ncol(Zt)/2
            Ztm <- Zt - p_i * 2
            A <- crossprod(Ztm)/ (2 * sum(p_i*(1-p_i)))
          }

        } else if(relationship.matrix.ogc=="CM"){
          #CM SCHAETZER
          Ztm <- rbind(Zt==0, Zt==1, Zt==2)
          A <- crossprod(Ztm) / ncol(Zt)
        } else if(relationship.matrix.ogc=="CE"){
          Ztm <- rbind(Zt==0, Zt==1, Zt==2)
          A <- crossprod(Ztm)
          A <- (A^2 - 0.5*A)/(nrow(Zt)^2)

        } else if(relationship.matrix.ogc=="non_stand"){
          A <- crossprod(Zt) / nrow(Zt)
        }

        Indiv <- paste0("Indi", 1:length(BV))
        colnames(A) <- rownames(A) <- names(BV) <- Indiv
        cont <- cbind(1,0.5,0.5)
        colnames(cont) <- c("age", "male", "female")
        Born <- rep(1, length(BV))
        Breed <- rep("Breed1", length(BV))
        herd <- rep(NA, length(BV))
        isCandidate <- rep(TRUE, length(BV))
        phen <- data.frame(Indiv, Born, Breed, BV, Sex, herd, isCandidate)

        sKin <- A/2
        colnames(sKin) <- rownames(sKin) <- Indiv
        cand <- optiSel::candes(phen = phen, sKin = sKin, cont = cont)
        con <- list()
        if(length(ogc.uniform)>0){
          con$uniform = ogc.uniform
        }
        if(length(ogc.lb)>0){
          con$lb = ogc.lb
        }
        if(length(ogc.ub)>0){
          con$ub = ogc.ub
        }
        if(length(ogc.ub.sKin)>0){
          con$ub.sKin = ogc.ub.sKin
        }
        if(length(ogc.lb.BV)>0){
          con$lb.BV = ogc.lb.BV
        }
        if(length(ogc.ub.BV)>0){
          con$ub.BV = ogc.ub.BV
        }
        if(length(ogc.eq.BV)>0){
          con$eq.BV = ogc.eq.BV
        }



        if(length(ogc.ub.sKin.increase)>0){
          con$ub.sKin = ogc.ub.sKin.increase + cand$current[2,4]
        }
        if(length(ogc.lb.BV.increase)>0){
          con$lb.BV = ogc.lb.BV.increase + cand$current[1,4]
        }

        Offspring <- optiSel::opticont(ogc.target, cand, con, quiet = !verbose)

        contribution <- list(Offspring$parent$oc[Sex=="male"], Offspring$parent$oc[Sex=="female"])

        if(verbose){
          cat(paste0(sum(contribution[[1]]>0), " male individuals with positive contribution ((", sum(contribution[[1]]>(0.001 * max(contribution[[1]]))), " with major contribution).\n"))
          cat(paste0(sum(contribution[[2]]>0), " female individuals with positive contribution ((", sum(contribution[[2]]>(0.001 * max(contribution[[2]]))), " with major contribution)\n."))
        }
      }
    }

  } else{
    breeding.size.total <- nrow(fixed.breeding)
    sex.animal <- fixed.breeding[,7] <- stats::rbinom(breeding.size.total, 1, fixed.breeding[,7]) +1
    breeding.size <- c(sum(fixed.breeding[,7]==1), sum(fixed.breeding[,7]==2))

  }

  if(length(threshold.selection)>0){
    if(length(best[[1]])>0){
      eval(parse(text=paste0("keeps <- which(best[[1]][,4]",threshold.sign,"threshold.selection)")))

      best[[1]] <- best[[1]][keeps,]

      if(verbose) cat(paste0(length(keeps), " sires fullfil threshold selection!\n"))
      selection.size[1] <- nrow(best[[1]])
      if(max.offspring[1]*selection.size[1] < breeding.size[1]){
        breeding.size[1] <- max.offspring[1]*selection.size[1]
        if(verbose) cat("Breeding size reduced based on threshold selection.\n")
      }
    }
    if(length(best[[2]])>0){
      eval(parse(text=paste0("keeps <- which(best[[2]][,4]",threshold.sign,"threshold.selection)")))
      best[[2]] <- best[[2]][keeps,]
      if(verbose) cat(paste0(length(keeps), " dams fullfil threshold selection!\n"))
      selection.size[2] <- nrow(best[[2]])
      if(max.offspring[2]*selection.size[2] < breeding.size[2]){
        breeding.size[2] <- max.offspring[2]*selection.size[2]
        if(verbose) cat("Breeding size reduced based on threshold selection")
      }
    }

    breeding.size.total <- sum(breeding.size)
  }

  sample_prob <- list()
  for(sex in 1:2){
    if(length(best[[sex]])>0){
      if(ogc){
        sample_prob[[sex]] <- contribution[[sex]]
      } else if(length(best.selection.manual.ratio[[sex]])==0){
        individual_bv <- best[[sex]][,5]
        min_bv <- min(individual_bv)
        max_bv <- max(individual_bv)
        if(min_bv==max_bv){
          sample_prob[[sex]] <- rep(1, length(individual_bv))
        } else{
          sample_prob[[sex]] <- (individual_bv- min_bv) / (max_bv-min_bv) * (best.selection.ratio[[sex]]-1) + 1
        }
      } else{
        sample_prob[[sex]] <- best.selection.manual.ratio[[sex]]
      }

      if(length(unique(sample_prob[[sex]]))==1){
        sample_prob[[sex]] <- NULL
      }
    }
  }
  sample_prob[[3]] <- "placeholder"



  if(length(fixed.breeding.best)>0){
    fixed.breeding <- matrix(0, nrow=nrow(fixed.breeding.best), ncol=7)

    for(index in 1:nrow(fixed.breeding.best)){
      fixed.breeding[index,1:3] <- best[[fixed.breeding.best[index,1]]][fixed.breeding.best[index,2], 1:3]
      fixed.breeding[index,4:6] <- best[[fixed.breeding.best[index,3]]][fixed.breeding.best[index,4], 1:3]
    }

    fixed.breeding[,7] <- fixed.breeding.best[,5]
    if(breeding.size.total==0){
      breeding.size.total <- nrow(fixed.breeding)
    }
    sex.animal <- fixed.breeding[,7] <- stats::rbinom(breeding.size.total, 1, fixed.breeding[,7]) +1
    breeding.size <- c(sum(fixed.breeding[,7]==1), sum(fixed.breeding[,7]==2))
  }



  if(store.comp.times){
    comp.times[6] <- as.numeric(Sys.time())
  }
}

  #######################################################################
  ############# Gene editing based on Simianer et al. 2018 ##############
  #######################################################################
{
    if(gene.editing){
      if(estimate.u){
        effect_order<- sort(abs(population$info$u_hat[[length(population$info$u_hat)]]), index.return=TRUE, decreasing=TRUE)$ix
        direction <- (population$info$u_hat[[length(population$info$u_hat)]] < 0)
        distance_actual <- integer(length(effect_order))
        for(index in 1:length(effect_order)){
          distance_actual[index] <- min(abs(effect_order[index]-population$info$effect.p))
        }
        if(length(population$info$editing_info)==0){
          population$info$editing_info <- list()
        }
        population$info$editing_info[[length(population$info$editing_info)+1]] <- cbind(effect_order, direction[effect_order], distance_actual, deparse.level = 0)
      }
      if(gwas.u){
        effect_order <- sort(abs(population$info$gwas_hat[[length(population$info$gwas_hat)]]), index.return=TRUE, decreasing=TRUE)$ix
        direction <- (population$info$gwas_hat[[length(population$info$gwas_hat)]] <0 )
        distance_actual <- integer(length(effect_order))
        for(index in 1:length(effect_order)){
          distance_actual[index] <- min(abs(effect_order[index]-population$info$effect.p))
        }
        if(length(population$info$editing_info)==0){
          population$info$editing_info <- list()
        }
        population$info$editing_info[[length(population$info$editing_info)+1]] <- cbind(effect_order, direction[effect_order], distance_actual, deparse.level = 0)
      }
    }



  if(gene.editing.best){
    for(sex in (1:2)[gene.editing.best.sex]){
      if(length(best[[sex]])>0){
        for(index in 1:nrow(best[[sex]])){
          activ <- best[[sex]][index,]
          population$breeding[[activ[1]]][[activ[2]]][[activ[3]]] <- edit_animal(population, activ[1], activ[2], activ[3], nr.edits, decodeOriginsU=decodeOriginsU, bit.storing=bit.storing, nbits=nbits)
        }
        for(index in 1:nrow(best[[sex]])){
          activ <- best[[sex]][index,]
          if(length(population$breeding[[activ[1]]][[activ[2]]][[activ[3]]])>=17){
            population$breeding[[activ[1]]][[activ[2]]][[activ[3]]][[17]] <- c(population$breeding[[activ[1]]][[activ[2]]][[activ[3]]][[17]],population$breeding[[activ[1]]][[6+activ[2]]][,activ[3]])
          } else{
            population$breeding[[activ[1]]][[activ[2]]][[activ[3]]][[17]] <- population$breeding[[activ[1]]][[6+activ[2]]][,activ[3]]
          }
          population$breeding[[activ[1]]][[activ[2]]][[activ[3]]][[15]] <- rep(0L, population$info$bv.nr)
          activ_bv <- population$info$bv.random.activ
          if(length(activ_bv)>0){
            temp_out <- calculate.bv(population, activ[1], activ[2], activ[3], activ_bv, import.position.calculation=import.position.calculation, decodeOriginsU=decodeOriginsU, bit.storing=bit.storing, nbits=nbits, output_compressed=FALSE, bv.ignore.traits=bv.ignore.traits)
            population$breeding[[activ[1]]][[6+activ[2]]][activ_bv,activ[3]] <- temp_out[[1]]
          }
          population$breeding[[activ[1]]][[activ[2]]][[activ[3]]][[17]] <- c(population$breeding[[activ[1]]][[activ[2]]][[activ[3]]][[17]] ,population$breeding[[activ[1]]][[6+activ[2]]][,activ[3]])
        }
      }
    }
  }
  }

  #######################################################################
  ############## Old Culling model - make ready for delete ##############
  #######################################################################
{
  if(is.matrix(reduce.group)==0 && length(reduce.group)>0){
    reduce.group <- t(reduce.group)
  }
  if(length(reduce.group)>0){
    reduce.group <- reduce.group[reduce.group[,1]!=0,]
  }
  if(length(reduce.group)>0){
    for(index in 1:nrow(reduce.group)){
      activ.reduce <- reduce.group[index,]
      group.animals <- (population$breeding[[activ.reduce[1]]][[activ.reduce[2]+4]] == activ.reduce[4])
      n.animals <- sum(group.animals)
      to.kill <- n.animals - activ.reduce[3]
      animal.position <- unique(c(0,group.animals * 1:length(group.animals)))[-1]
      if(to.kill>0){
        if(reduce.group.selection=="random"){
          death <- sample(animal.position, to.kill)

        } else if(reduce.group.selection=="selection"){
          if(multiple.bve=="add"){
            if(verbose) cat(population$breeding[[activ.reduce[1]]][[activ.reduce[2]+2]])
            bve.sum <- colSums(rbind(population$breeding[[activ.reduce[1]]][[activ.reduce[2]+2]][,animal.position]*multiple.bve.weights[[sex]],0))
          } else if(multiple.bve=="ranking"){
            ranking <- population$breeding[[index]][[sex+2]][,relevant.animals]
            for(bven in 1:population$info$bv.nr){
              order <- sort(ranking[bven,], index.return=TRUE, decreasing=selection.highest[activ.reduce[2]])$ix
              ranking[bven,order] <- length(order):1
            }
            bve.sum <- colSums(ranking*multiple.bve.weights[[sex]])
          }
          chosen.animals <- sort(bve.sum, index.return=TRUE, decreasing=selection.highest[activ.reduce[2]])$ix
          death <- chosen.animals[(length(chosen.animals)-to.kill+1):length(chosen.animals)]
        }
        for(modanimal in death){
          population$breeding[[activ.reduce[1]]][[activ.reduce[2]]][[modanimal]][[18]] <- c(population$breeding[[activ.reduce[1]]][[activ.reduce[2]+4]][modanimal], current.gen)
          population$breeding[[activ.reduce[1]]][[activ.reduce[2]+4]][modanimal] <- -1
        }
      }


    }
  }

}

  #######################################################################
  #################### Generation of new individuals ####################
  #######################################################################
{

  # preparation of the population - list // counting for max.offspring etc.
{

  if(length(population$breeding)<=current.gen && sum(breeding.size)>0 ){
    population$breeding[[current.gen+1]] <- list()
    population$info$size <- rbind(population$info$size, 0)
  }

  selection.rate <- list(numeric(selection.size[1]), numeric(selection.size[2]))
  selection.rate.litter <- list(numeric(selection.size[1]), numeric(selection.size[2]))
  activ.selection.size <- selection.size
  availables.m <- 1:selection.size[1]
  availables.f <- 1:selection.size[2]
  availables <- list(availables.m, availables.f)

  if(relative.selection){
    if(sum(best[[1]][,4]<0)>0){
      best[[1]][(best[[1]][,4]<0),4] <- 0
    }
    if(sum(best[[2]][,4]<0)>0){
      best[[2]][(best[[2]][,4]<0),4] <- 0
    }

    cum.sum.m <- cumsum(best[[1]][,4])
    cum.sum.f <- cumsum(best[[2]][,4])
    sum.m <- sum(best[[1]][,4])
    sum.f <- sum(best[[2]][,4])
    cum.sum <- list(cum.sum.m, cum.sum.f)
    sum <- list(sum.m,sum.f)
  }
  current.size <- c(1,1)
  for(sex in 1:2){
    if(sum(breeding.size)>0){
      if(length(population$breeding[[current.gen+1]])<=2 || length(population$breeding[[current.gen+1]][[sex]])==0){
        population$breeding[[current.gen+1]][[2+sex]] <- matrix(0, nrow=population$info$bv.nr, ncol=breeding.size[sex])
        population$breeding[[current.gen+1]][[4+sex]] <- rep(new.class, breeding.size[sex])
        population$breeding[[current.gen+1]][[6+sex]] <- matrix(0, nrow=population$info$bv.nr, ncol=breeding.size[sex])
        population$breeding[[current.gen+1]][[8+sex]] <- matrix(NA, nrow=population$info$bv.nr, ncol=breeding.size[sex])
        population$breeding[[current.gen+1]][[10+sex]] <- rep(time.point, breeding.size[sex])
        population$breeding[[current.gen+1]][[12+sex]] <- rep(creating.type, breeding.size[sex])
        if(copy.individual){
          population$breeding[[current.gen+1]][[14+sex]] <- rep(0, breeding.size[sex])
        } else{
          population$breeding[[current.gen+1]][[14+sex]] <- seq(population$info$next.animal, population$info$next.animal + breeding.size[sex] -1, length.out= breeding.size[sex])
          population$info$next.animal <- population$info$next.animal + breeding.size[sex]
        }
        population$breeding[[current.gen+1]][[16+sex]] <- rep(NA, breeding.size[sex])
        population$breeding[[current.gen+1]][[18+sex]] <- matrix(0, nrow=population$info$bv.nr, ncol=breeding.size[sex])
        population$breeding[[current.gen+1]][[20+sex]] <- matrix(0, nrow=population$info$bv.nr, ncol=breeding.size[sex])
        if(copy.individual){
          #placeholder
          population$breeding[[current.gen+1]][[22+sex]] <- rep(-1, breeding.size[sex])
        } else{
          population$breeding[[current.gen+1]][[22+sex]] <- rep(time.point, breeding.size[sex])
        }
        population$breeding[[current.gen+1]][[24+sex]] <- rep(NA, breeding.size[sex])
        population$breeding[[current.gen+1]][[26+sex]] <- matrix(0, nrow=population$info$bv.nr, ncol=breeding.size[sex])
        population$breeding[[current.gen+1]][[28+sex]] <- matrix(0, nrow=population$info$bv.nr, ncol=breeding.size[sex])

        population$breeding[[current.gen+1]][[30+sex]] <- rep(0, ncol=breeding.size[sex])

      } else{
        current.size[sex] <- length(population$breeding[[current.gen+1]][[4+sex]]) + 1
        population$breeding[[current.gen+1]][[2+sex]] <- cbind(population$breeding[[current.gen+1]][[sex+2]], matrix(0, nrow= population$info$bv.nr, ncol=breeding.size[sex]))
        population$breeding[[current.gen+1]][[4+sex]] <- c(population$breeding[[current.gen+1]][[sex+4]], rep(new.class, breeding.size[sex]))
        population$breeding[[current.gen+1]][[6+sex]] <- cbind(population$breeding[[current.gen+1]][[6+sex]], matrix(0, nrow= population$info$bv.nr, ncol=breeding.size[sex]))
        population$breeding[[current.gen+1]][[8+sex]] <- cbind(population$breeding[[current.gen+1]][[8+sex]], matrix(NA, nrow= population$info$bv.nr, ncol=breeding.size[sex]))
        population$breeding[[current.gen+1]][[10+sex]] <- c(population$breeding[[current.gen+1]][[sex+10]], rep(time.point, breeding.size[sex]))
        population$breeding[[current.gen+1]][[12+sex]] <- c(population$breeding[[current.gen+1]][[sex+12]], rep(creating.type, breeding.size[sex]))

        if(copy.individual){
          population$breeding[[current.gen+1]][[14+sex]] <- c(population$breeding[[current.gen+1]][[14+sex]], rep(0,breeding.size[sex]))
        } else{
          population$breeding[[current.gen+1]][[14+sex]] <- c(population$breeding[[current.gen+1]][[14+sex]], seq(population$info$next.animal, population$info$next.animal + breeding.size[sex] -1, length.out= breeding.size[sex]))
          population$info$next.animal <- population$info$next.animal + breeding.size[sex]
        }
        population$breeding[[current.gen+1]][[16+sex]] <- c(population$breeding[[current.gen+1]][[sex+16]], rep(NA, breeding.size[sex]))
        population$breeding[[current.gen+1]][[18+sex]] <- cbind(population$breeding[[current.gen+1]][[18+sex]], matrix(0, nrow= population$info$bv.nr, ncol=breeding.size[sex]))
        population$breeding[[current.gen+1]][[20+sex]] <- cbind(population$breeding[[current.gen+1]][[20+sex]], matrix(0, nrow= population$info$bv.nr, ncol=breeding.size[sex]))
        if(copy.individual){
          #placeholder
          population$breeding[[current.gen+1]][[22+sex]] <- c(population$breeding[[current.gen+1]][[sex+22]], rep(-1, breeding.size[sex]))

        } else{
          population$breeding[[current.gen+1]][[22+sex]] <- c(population$breeding[[current.gen+1]][[sex+22]], rep(time.point, breeding.size[sex]))

        }#
        population$breeding[[current.gen+1]][[24+sex]] <- c(population$breeding[[current.gen+1]][[sex+24]], rep(NA, breeding.size[sex]))
        population$breeding[[current.gen+1]][[26+sex]] <- cbind(population$breeding[[current.gen+1]][[26+sex]], matrix(0, nrow= population$info$bv.nr, ncol=breeding.size[sex]))
        population$breeding[[current.gen+1]][[28+sex]] <- cbind(population$breeding[[current.gen+1]][[28+sex]], matrix(0, nrow= population$info$bv.nr, ncol=breeding.size[sex]))
        population$breeding[[current.gen+1]][[30+sex]] <- c(population$breeding[[current.gen+1]][[30+sex]], rep(0, breeding.size[sex]))

      }
      if(length(population$breeding[[current.gen+1]][[sex]])==0){
        population$breeding[[current.gen+1]][[sex]] <- list()
      }

    }
  }
  if(length(praeimplantation)>0){
    nsnps <- c(0,population$info$cumsnp)
    praeimplantation <- cbind(praeimplantation, nsnps[praeimplantation[,1]] + praeimplantation[,2])
    praeimplantation.max <- list()
    for(sex in 1:2){
      praeimplantation.max[[sex]] <- numeric(nrow(best[[sex]]))
      for(index in 1:length(best[[sex]])){
        activ <- best[[sex]][index,]
        hap1 <- population$breeding[[activ[1]]][[activ[2]]][[activ[3]]][[9]][praeimplantation[,4]]
        hap2 <- population$breeding[[activ[1]]][[activ[2]]][[activ[3]]][[10]][praeimplantation[,4]]
        pos <- ((hap1==praeimplantation[,3]) + (hap2==praeimplantation[,3]))>0
        praeimplantation.max[[sex]][index] <- sum(pos)
      }
    }
    if(verbose) cat("praeimplantation bisher nicht in Selektion enthalten!")
  }


##  store.effect.freq and multiple correlated bvs deactivated
  if(store.comp.times.generation){
    pre_stuff <- 0
    generation_stuff <- 0
    bv_stuff <- 0
  }

  if(length(name.cohort)==0 && breeding.size.total>0){
    name.cohort <- paste0("Cohort_", population$info$cohort.index)
    population$info$cohort.index <- population$info$cohort.index + 1
  }

  if(copy.individual==FALSE && breeding.size.total>0){
    if(nrow(best[[1]])==0){
      if(same.sex.activ==FALSE || same.sex.sex > 0){
        if(verbose){
          if(!dh.mating && !selfing.mating) {cat("No male individuals provided for reproduction. Automatically allow female X female matings.\n")}
        }
      }
      same.sex.activ <- TRUE
      same.sex.sex <- selfing.sex <- 1
    } else if(nrow(best[[2]])==0){
      if(same.sex.activ==FALSE || same.sex.sex < 1){
        if(verbose){
          if(!dh.mating && !selfing.mating) cat("No female individuals provided for reproduction. Automatically allow male X male matings.\n")
        }
      }
      same.sex.activ <- TRUE
      same.sex.sex <- selfing.sex <-  0
    }
  }

}




  # Actually generation ((need to revisit parallel generation - currently basically no computational gains))
  if(parallel.generation && breeding.size.total>0){

    if(requireNamespace("doParallel", quietly = TRUE)) {

      if(length(name.cohort)>0){
        if(verbose) cat(paste0("Start generation of new individuals (cohort: ", name.cohort,").\n"))
      } else{
        if(verbose) cat("Start generation of new individuals.\n")
      }


      if(store.comp.times.generation){
        tick <- as.numeric(Sys.time())
      }
      info_father_list <- info_mother_list <- matrix(0, nrow=breeding.size.total, ncol=5)

      if(length(repeat.mating.fixed)>0){
        repeat.mating.fixed <- rep(repeat.mating.fixed, length.out = breeding.size.total)
      }
      rmf <- 1

      runs <- 0
      for(animal.nr in 1:breeding.size.total){

        sex <- sex.animal[animal.nr]
        if(length(fixed.breeding)>0){
          info.father <- fixed.breeding[animal.nr,1:3]
          info.mother <- fixed.breeding[animal.nr,4:6]
        } else{
          if(runs>0){
            runs <- runs - 1
          } else{

            if(copy.individual){
              mating.temp <- sample(1:nrow(population$info$repeat.mating.copy), prob = population$info$repeat.mating.copy[,2], 1)
              repeat.mating.temp <- population$info$repeat.mating.copy[mating.temp,1]
            } else{
              mating.temp <- sample(1:nrow(population$info$repeat.mating), prob = population$info$repeat.mating[,2], 1)
              repeat.mating.temp <- population$info$repeat.mating[mating.temp,1]
            }
            if(length(repeat.mating.fixed)>0){
              repeat.mating.temp <- repeat.mating.fixed[rmf]
              rmf <- rmf +1
            }

            runs <- repeat.mating.temp - 1

            accepted <- FALSE
            while(accepted==FALSE){

              if(selfing.mating==FALSE){
                sex1 <- 1
                sex2 <- 2
                if(same.sex.activ==FALSE && relative.selection==FALSE){

                  number1 <- availables[[sex1]][sample(1:activ.selection.size[sex1],1, prob=sample_prob[[sex1]][availables[[sex1]]])]
                  number2 <- availables[[sex2]][sample(1:activ.selection.size[sex2],1, prob=sample_prob[[sex2]][availables[[sex2]]])]

                  info.father <- best[[sex1]][number1,]
                  info.mother <- best[[sex2]][number2,]
                } else if(same.sex.activ==FALSE && relative.selection){
                  info.father <- best[[sex1]][sum(cum.sum[[sex1]] <stats::runif(1,0,sum[[1]])) +1 ,1:3]
                  info.mother <- best[[sex2]][sum(cum.sum[[sex2]] <stats::runif(1,0,sum[[2]])) +1 ,1:3]
                } else{
                  sex1 <- stats::rbinom(1,1,same.sex.sex) + 1 # ungleichviele tiere erhoeht

                  sex2 <- stats::rbinom(1,1,same.sex.sex) + 1

                  number1 <- availables[[sex1]][sample(1:activ.selection.size[sex1],1, prob=sample_prob[[sex1]][availables[[sex1]]])]
                  number2 <- availables[[sex2]][sample(1:activ.selection.size[sex2],1, prob=sample_prob[[sex2]][availables[[sex2]]])]
                  test <- 1
                  while(same.sex.selfing==FALSE && number1==number2 && test < 100){
                    number2 <- availables[[sex2]][sample(1:activ.selection.size[sex2],1, prob=sample_prob[[sex2]][availables[[sex2]]])]
                    test <- test+1
                    if(test==100 && number1==number2){
                      warning("Only one remaining individual in the selected cohorts. // proceed with selfing.")
                    }
                  }

                  if(relative.selection==FALSE){
                    info.father <- best[[sex1]][number1,]
                    # waehle fuers bveite Tier ein Tier aus der Gruppe der Nicht-besten Tiere
                    if(martini.selection){
                      options <- (1:population$info$size[best[[sex2]][1,1],sex2])[-best[[sex2]][,3]]
                      number2 <- sample(options,1)
                      info.mother <- c(best[[sex2]][1,1:2], number2)
                    } else{
                      info.mother <- best[[sex2]][number2,]
                    }
                  } else{
                    info.father <- best[[sex1]][sum(cum.sum[[sex1]] <stats::runif(1,0,cum.sum[[sex1]])) +1 ,1:3]
                    info.mother <- best[[sex2]][sum(cum.sum[[sex2]] <stats::runif(1,0,cum.sum[[sex2]])) +1 ,1:3]
                  }

                }
              } else{
                sex1 <- sex2 <- stats::rbinom(1,1,selfing.sex)+1
                if(length(availables[[sex1]])==0){
                  sex1 <- sex2 <- 3 - sex1
                }
                number1 <- number2 <- availables[[sex1]][sample(1:activ.selection.size[sex1],1, prob=sample_prob[[sex1]][availables[[sex1]]])]
                info.father <- best[[sex1]][number1,]
                info.mother <- best[[sex2]][number2,]

              }

              accepted <- TRUE
              if(max.mating.pair < Inf){

                if(info.father[2]==info.mother[2]){
                  if(info.father[1]>info.mother[1] || (info.father[1]==info.mother[1] & info.father[3]> info.mother[3])){
                    info_mother_temp <- info.mother
                    info.mother <- info.father
                    info.father <- info_mother_temp

                  }
                }

                if(  sum((colSums(t(info_father_list) == info.father) + colSums(t(info_mother_list) == info.mother))==10)  >= max.mating.pair){
                  accepted <- FALSE
                }
              }
            }

            if(martini.selection==FALSE){

              selection.rate[[sex1]][number1] <- selection.rate[[sex1]][number1] + repeat.mating.temp
              selection.rate[[sex2]][number2] <- selection.rate[[sex2]][number2] + repeat.mating.temp

              selection.rate.litter[[sex1]][number1] <- selection.rate.litter[[sex1]][number1] + 1
              selection.rate.litter[[sex2]][number2] <- selection.rate.litter[[sex2]][number2] + 1

              if(selection.rate[[sex1]][number1] >= max.offspring[sex1]){
                activ.selection.size[sex1] <-  activ.selection.size[sex1] - 1
                availables[[sex1]] <- availables[[sex1]][availables[[sex1]]!=number1]
              } else if(selection.rate.litter[[sex1]][number1] >= max.litter[sex1]){
                activ.selection.size[sex1] <-  activ.selection.size[sex1] - 1
                availables[[sex1]] <- availables[[sex1]][availables[[sex1]]!=number1]
              }

              if(selection.rate[[sex2]][number2] >= max.offspring[sex2]){
                if(sex1!= sex2 || number1 != number2){
                  activ.selection.size[sex2] <-  activ.selection.size[sex2] -1
                }
                availables[[sex2]] <- availables[[sex2]][availables[[sex2]]!=number2]
              } else if(selection.rate.litter[[sex2]][number2] >= max.litter[sex2]){
                if(sex1!= sex2 || number1 != number2){
                  activ.selection.size[sex2] <-  activ.selection.size[sex2] -1
                }
                availables[[sex2]] <- availables[[sex2]][availables[[sex2]]!=number2]
              }

            }
          }


        }


        info_father_list[animal.nr,] <- info.father
        info_mother_list[animal.nr,] <- info.mother
      }

      if(store.comp.times.generation){
        tock <- as.numeric(Sys.time())
        pre_stuff <- tock-tick
      }

      if(Sys.info()[['sysname']]=="Windows"){
        if (requireNamespace("doParallel", quietly = TRUE)) {
          doParallel::registerDoParallel(cores=ncore.generation)
        } else{
          stop("Usage of doParallel without being installed!")
        }
        if(length(randomSeed.generation)>0){
          if (requireNamespace("doRNG", quietly = TRUE)) {
            doRNG::registerDoRNG(seed=randomSeed.generation)
          } else{
            stop("Usage of doRNG without being installed!")
          }
        }
      }


      for(sex_running in 1:2){

        if(store.comp.times.generation){
          tick <- as.numeric(Sys.time())
        }

        if(Sys.info()[['sysname']]=="Windows"){
          if (requireNamespace("foreach", quietly = TRUE)) {
          } else{
            stop("Usage of foreach without being installed!")
          }
          new_animal <- foreach::foreach(indexb=(1:breeding.size.total)[sex.animal==sex_running],
                                         .packages="MoBPS") %dopar% {

                                           child <- generation.individual(indexb,
                                                                          population, info_father_list, info_mother_list,
                                                                          copy.individual, mutation.rate, remutation.rate, recombination.rate,
                                                                          recom.f.indicator, duplication.rate, duplication.length,
                                                                          duplication.recombination, delete.same.origin,
                                                                          (gene.editing.offspring * gene.editing.offspring.sex[1]), nr.edits,
                                                                          gen.architecture.m, gen.architecture.f, decodeOriginsU,
                                                                          current.gen, save.recombination.history, phenotyping.child,
                                                                          dh.mating, share.genotyped, added.genotyped, genotyped.array,
                                                                          dh.sex, n.observation)
                                           child

          }
        } else {
          activ_stuff <- (1:breeding.size.total)[sex.animal==sex_running]

          if (requireNamespace("parallel", quietly = TRUE)) {
          } else{
            stop("Use of parallel without being installed!")
          }
          new_animal <- parallel::mclapply(activ_stuff, function(x) generation.individual(x,
                                                                                population, info_father_list, info_mother_list,
                                                                                copy.individual, mutation.rate, remutation.rate, recombination.rate,
                                                                                recom.f.indicator, duplication.rate, duplication.length,
                                                                                duplication.recombination, delete.same.origin,
                                                                                (gene.editing.offspring * gene.editing.offspring.sex[1]), nr.edits,
                                                                                gen.architecture.m, gen.architecture.f, decodeOriginsU,
                                                                                current.gen, save.recombination.history, phenotyping.child,
                                                                                dh.mating, share.genotyped, added.genotyped, genotyped.array,
                                                                                dh.sex, n.observation),
                                 mc.cores=ncore.generation)
        }


        prev_ani <- length(population$breeding[[current.gen+1]][[sex_running]])
        present_before <- population$info$size[current.gen+1,sex_running]
        population$breeding[[current.gen+1]][[sex_running]] <-  c(population$breeding[[current.gen+1]][[sex_running]], new_animal)
        population$info$size[current.gen+1,sex_running] <- length(population$breeding[[current.gen+1]][[sex_running]])
        current.size[sex_running] <- length(population$breeding[[current.gen+1]][[sex_running]]) +1

        if(length(new_animal)>0){
          for(index6 in 1:length(new_animal) + prev_ani){
            if(copy.individual){
              first_copy <- population$breeding[[current.gen+1]][[sex_running]][[index6]][[21]]
              new_copy <- rbind(population$breeding[[first_copy[1,1]]][[first_copy[1,2]]][[first_copy[1,3]]][[21]],
                                c(current.gen+1, sex_running, index6))
              storage.mode(new_copy) <- "integer"
              for(index7 in 1:nrow(new_copy)){
                population$breeding[[new_copy[index7,1]]][[new_copy[index7,2]]][[new_copy[index7,3]]][[21]] <- new_copy
              }
              population$breeding[[current.gen+1]][[sex_running+22]][index6] <- population$breeding[[first_copy[1,1]]][[first_copy[1,2]+22]][first_copy[1,3]]

            } else{
              population$breeding[[current.gen+1]][[sex_running]][[index6]][[21]] <- cbind(current.gen+1, sex_running, index6, deparse.level = 0)
              storage.mode(population$breeding[[current.gen+1]][[sex_running]][[index6]][[21]]) <- "integer"
            }
          }
        }

        if(store.comp.times.generation){
          tack <- as.numeric(Sys.time())
          generation_stuff <- tack-tick + generation_stuff
        }
        if(store.effect.freq){
          store.effect.freq <- FALSE
          if(verbose) cat("Effect-Freq not available in parallel computing.\n")
        }

        index_loop <- NULL
        if(length(new_animal)>0) index_loop <- 1:length(new_animal)

        if(Sys.info()[['sysname']]!="Windows"){
          doParallel::registerDoParallel(cores=ncore.generation)
        }

        new.bv.list <- foreach::foreach(indexb=index_loop,
                              .packages=c("MoBPS", "miraculix")) %dopar% {

                              info.father <- info_father_list[indexb,]
                              info.mother <- info_mother_list[indexb,]
                              new.bv <-  new.bve <- numeric(population$info$bv.nr)
                              new.bv_approx <- rep(NA, population$info$bv.nr)
                              activ_bv <- population$info$bv.random.activ

                              if(length(activ_bv)>0){
                                temp_out <- calculate.bv(population, current.gen+1, sex_running, indexb + present_before, activ_bv, import.position.calculation=import.position.calculation, decodeOriginsU=decodeOriginsU, store.effect.freq=store.effect.freq, bit.storing=bit.storing, nbits=nbits, output_compressed=FALSE,bv.ignore.traits=bv.ignore.traits)
                                new.bv[activ_bv] <- temp_out[[1]]
                                if(store.effect.freq){
                                  if(length(population$info$store.effect.freq) < (current.gen+1) || length(population$info$store.effect.freq[[current.gen+1]])==0){
                                    colnames(temp_out[[2]]) <- c("Homo0", "Hetero", "Homo1")
                                    rownames(temp_out[[2]]) <- population$info$snp.name[population$info$effect.p]
                                    population$info$store.effect.freq[[current.gen+1]] <- temp_out[[2]]
                                  } else{
                                    population$info$store.effect.freq[[current.gen+1]] <- population$info$store.effect.freq[[current.gen+1]] + temp_out[[2]]
                                  }
                                }
                              }

                              if(population$info$bv.calc > 0  && population$info$bv.random[population$info$bv.calc]){
                                means <- 0.5*(population$breeding[[info.father[1]]][[6+info.father[2]]][[population$info$bv.calc:population$info$bv.nr,info.father[3]]] + population$breeding[[info.mother[1]]][[6+info.mother[2]]][[population$info$bv.calc:population$info$bv.nr,info.mother[3]]])
                                varp <- kinship.emp(list(population$breeding[[info.father[1]]][[info.father[2]]][[info.father[3]]], population$breeding[[info.mother[1]]][[info.mother[2]]][[info.mother[3]]]))
                                varp <- (2 - (2* (varp[1,1]-0.5) + 2 * (varp[2,2]- 0.5)))/4

                                if(FALSE){
                                  new.bv[population$info$bv.calc:population$info$bv.nr] <- stats::rnorm(1, mean=means, sd= sqrt(var*population$info$bv.random.variance[bven]))
                                } else{
                                  if(population$info$bv.calc==1){
                                    population$info$current.bv.random.variance <- varp * population$info$bv.random.variance
                                    bv.var <- diag.mobps(sqrt(population$info$current.bv.random.variance)) %*%population$info$current.bv.correlation %*% diag.mobps(sqrt(population$info$current.bv.random.variance))
                                    single.mean <- means
                                  } else{
                                    population$info$current.bv.random.variance <- c(population$info$bv.random.variance[1:(population$info$bv.calc-1)],varp * population$info$bv.random.variance[population$info$bv.calc:population$info$bv.nr])
                                    AA <- diag.mobps(sqrt(population$info$current.bv.random.variance)[1:(population$info$bv.calc-1)]) %*% population$info$current.bv.correlation[1:(population$info$bv.calc-1), 1:(population$info$bv.calc-1)]%*% diag.mobps(sqrt(population$info$current.bv.random.variance)[(1:(population$info$bv.calc-1))])
                                    BB <- diag.mobps(sqrt(population$info$current.bv.random.variance)[1:(population$info$bv.calc-1)]) %*%population$info$current.bv.correlation[1:(population$info$bv.calc-1), -(1:(population$info$bv.calc-1))]%*% diag.mobps(sqrt(population$info$current.bv.random.variance)[-(1:(population$info$bv.calc-1))])
                                    CC <- diag.mobps(sqrt(population$info$current.bv.random.variance)[-(1:(population$info$bv.calc-1))]) %*%population$info$current.bv.correlation[-(1:(population$info$bv.calc-1)), -(1:(population$info$bv.calc-1))] %*% diag.mobps(sqrt(population$info$current.bv.random.variance)[-(1:(population$info$bv.calc-1))])
                                    if (requireNamespace("MASS", quietly = TRUE)) {
                                      bv.var <- CC - t(BB) %*% MASS::ginv(AA) %*% BB
                                      single.mean <- means + t(BB) %*% MASS::ginv(AA) %*% ( new.bv[1:(population$info$bv.calc-1)]-mu1[1:(population$info$bv.calc-1)])
                                    } else{
                                      bv.var <- CC - t(BB) %*% solve(AA) %*% BB
                                      single.mean <- means + t(BB) %*% solve(AA) %*% ( new.bv[1:(population$info$bv.calc-1)]-mu1[1:(population$info$bv.calc-1)])
                                    }
                                  }

                                  bv.var_chol <- t(chol(bv.var))
                                  population$info$bv.correlation_col <- bv.var_chol

                                  random_part <- bv.var_chol %*% stats::rnorm(population$info$bv.nr - population$info$bv.calc +1, 0,1 )
                                  means_part <- single.mean
                                  new.bv[population$info$bv.calc:population$info$bv.nr] <-random_part + means_part
                                }

                              }


                              if(sum(population$info$is.combi)>0){
                                combis <- which(population$info$is.combi)
                                for(combi in combis){
                                  new.bv[combi] <- sum(new.bv * population$info$combi.weights[[combi]])
                                }
                              }

                              if(phenotyping.child=="mean" || phenotyping.child=="addobs"){
                                activ_father <- population$breeding[[current.gen+1]][[sex_running]][[indexb + present_before]][[7]]
                                activ_mother <- population$breeding[[current.gen+1]][[sex_running]][[indexb + present_before]][[8]]
                                if(copy.individual){
                                  new.bv_approx <- population$breeding[[info.father[1]]][[8+info.father[2]]][,info.father[3]]
                                } else{
                                  for(bven in 1:population$info$bv.nr){
                                    new.bv_approx[bven] <- mean(c(population$breeding[[activ_father[1]]][[8+activ_father[2]]][bven,activ_father[3]],population$breeding[[activ_mother[1]]][[8+activ_mother[2]]][bven,activ_mother[3]]))
                                  }
                                }

                              }
                              if(phenotyping.child=="obs" || phenotyping.child=="addobs"){
                                if(sum(n.observation)>0){
                                  if(phenotyping.child=="addobs"){
                                    prior_obs <- population$breeding[[info.father[1]]][[info.father[2]]][[info.father[3]]][[15]]
                                    total_obs <- prior_obs + n.observation
                                    new.bv_approx <- (new.bv_approx - population$breeding[[info.father[1]]][[info.father[2]+6]][,info.father[3]]) *  prior_obs / total_obs
                                  } else{
                                    total_obs <- n.observation
                                  }
                                  observation_reps <- sort(unique(c(0L,n.observation)))
                                  for(observation_rep in 2:length(observation_reps)){
                                    new.obs <- observation_reps[observation_rep] - observation_reps[observation_rep-1]
                                    temp_random <- matrix(stats::rnorm(population$info$bv.nr*new.obs,0,1), ncol=new.obs)
                                    active.traits <- (n.observation >=observation_reps[observation_rep])
                                    active.traits <- active.traits*(1:length(active.traits))
                                    for(bven in (1:population$info$bv.nr)[setdiff(active.traits, activ.trafo)]){
                                      if(is.na(new.bv_approx[bven]) & new.obs > 0){
                                        new.bv_approx[bven] <-  new.obs/(total_obs[bven]) * rowMeans(population$info$pheno.correlation %*% temp_random)[bven] * sqrt(sigma.e[bven])
                                      } else{
                                        new.bv_approx[bven] <-  new.bv_approx[bven] + new.obs/(total_obs[bven]) * rowMeans(population$info$pheno.correlation %*% temp_random)[bven] * sqrt(sigma.e[bven])
                                      }
                                    }
                                    for(bven in (1:population$info$bv.nr)[intersect(active.traits, activ.trafo)]){
                                      new_pheno <- (population$info$pheno.correlation %*% temp_random)[bven,] * sqrt(sigma.e[bven]) + new.bv[bven]
                                      new_pheno <- population$info$phenotypic.transform.function[[bven]](new_pheno) - new.bv[bven]
                                      if(is.na(new.bv_approx[bven]) & new.obs > 0){
                                        new.bv_approx[bven] <-  new.obs/(total_obs[bven]) * mean(new_pheno)
                                      } else{
                                        new.bv_approx[bven] <-  new.bv_approx[bven] + new.obs/(total_obs[bven]) * mean(new_pheno)
                                      }

                                    }


                                  }
                                  new.bv_approx <- new.bv + new.bv_approx
                                  new.bv_approx[n.observation==0] <- NA
                                }

                              }

                              new.bve <- population$breeding[[info.father[1]]][[2+info.father[2]]][,info.father[3]]
                              new.reli <- population$breeding[[info.father[1]]][[18+info.father[2]]][,info.father[3]]

                              temp1 <- c(new.bv, new.bv_approx, new.bve, new.reli)
                              temp1
                              }


        if(length(index_loop)>0){
          new.bv.list <- matrix(unlist(new.bv.list), ncol=length(new_animal))
          n_bv <- nrow(new.bv.list) / 4
          population$breeding[[current.gen+1]][[sex_running+6]][,(present_before+1):(present_before+length(new_animal))] <- new.bv.list[1:n_bv,,drop=FALSE]
          population$breeding[[current.gen+1]][[sex_running+8]][,(present_before+1):(present_before+length(new_animal))] <- new.bv.list[1:n_bv+n_bv,,drop=FALSE]
          if(copy.individual && copy.individual.keep.bve){
            population$breeding[[current.gen+1]][[sex_running+2]][,(present_before+1):(present_before+length(new_animal))] <- new.bv.list[1:n_bv+2*n_bv,,drop=FALSE]
            population$breeding[[current.gen+1]][[sex_running+18]][,(present_before+1):(present_before+length(new_animal))] <- new.bv.list[1:n_bv+3*n_bv,,drop=FALSE]
          }
        }

        if(store.comp.times.generation){
          tock <- as.numeric(Sys.time())
          bv_stuff <- tock - tack + bv_stuff
        }
      }

      if(Sys.info()[['sysname']]=="Windows" || TRUE){
        doParallel::stopImplicitCluster()
      }

      if(copy.individual){
        activ_prior <- length(population$breeding[[current.gen+1]][[sex_running]]) - length(new_animal) +1
        for(index in 1:nrow(info_father_list)){
          info.father <- info_father_list[index,]
          population$breeding[[current.gen+1]][[14+sex_running]][activ_prior] <- population$breeding[[info.father[1]]][[info.father[2]+14]][[info.father[3]]]
        }
      }
    } else{
      stop("Usage of doParallel without being installed!")
    }

  } else if(breeding.size.total>0){
    if(length(name.cohort)>0){
      if(verbose) cat(paste0("Start generation of new individuals (cohort: ", name.cohort,").\n"))
    } else{
      if(verbose) cat("Start generation of new individuals.\n")
    }

    pb_temp <- round(breeding.size.total/100)

    if(display.progress & verbose){
      pb <- utils::txtProgressBar(min = 0, max = breeding.size.total, style = 3)
    }

    runs <- 0

    if(max.mating.pair < Inf){
      info_father_list <- info_mother_list <- matrix(0, nrow=breeding.size.total, ncol=5)
    }

    if(length(repeat.mating.fixed)>0){
      repeat.mating.fixed <- rep(repeat.mating.fixed, length.out = breeding.size.total)
    }
    rmf <- 1 # repeat.mating temp-variable

    for(animal.nr in 1:breeding.size.total){
      if(store.comp.times.generation){
        tick <- as.numeric(Sys.time())
      }
      sex <- sex.animal[animal.nr]
      new.bv <-  new.bve <- new.reli <- individual.id <- numeric(population$info$bv.nr)
      new.bv_approx <- rep(NA, population$info$bv.nr)
      if(length(fixed.breeding)>0){
        info.father <- fixed.breeding[animal.nr,1:3]
        info.mother <- fixed.breeding[animal.nr,4:6]
      } else{
        if(runs>0){
          runs <- runs - 1
        } else{
          if(copy.individual){
            mating.temp <- sample(1:nrow(population$info$repeat.mating.copy), prob = population$info$repeat.mating.copy[,2], 1)
            repeat.mating.temp <- population$info$repeat.mating.copy[mating.temp,1]
          } else{
            mating.temp <- sample(1:nrow(population$info$repeat.mating), prob = population$info$repeat.mating[,2], 1)
            repeat.mating.temp <- population$info$repeat.mating[mating.temp,1]
          }

          if(length(repeat.mating.fixed)>0){
            repeat.mating.temp <- repeat.mating.fixed[rmf]
            rmf <- rmf +1
          }
          runs <- repeat.mating.temp - 1




          # Determine Sire/dam combination
          check_rel <- FALSE
          accepted <- FALSE
          max_counter <- 0
          max_rel_temp <- max_rel
          while(!check_rel || accepted==FALSE){
            if(selfing.mating==FALSE){
              sex1 <- 1
              sex2 <- 2
              if(same.sex.activ==FALSE && relative.selection==FALSE){
                number1 <- availables[[sex1]][sample(1:activ.selection.size[sex1],1, prob=sample_prob[[sex1]][availables[[sex1]]])]
                number2 <- availables[[sex2]][sample(1:activ.selection.size[sex2],1, prob=sample_prob[[sex2]][availables[[sex2]]])]
                info.father <- best[[sex1]][number1,]
                info.mother <- best[[sex2]][number2,]
              } else if(same.sex.activ==FALSE && relative.selection){
                info.father <- best[[sex1]][sum(cum.sum[[sex1]] <stats::runif(1,0,sum[[1]])) +1 ,1:3]
                info.mother <- best[[sex2]][sum(cum.sum[[sex2]] <stats::runif(1,0,sum[[2]])) +1 ,1:3]
              } else{
                sex1 <- stats::rbinom(1,1,same.sex.sex) + 1 # ungleichviele tiere erhoeht
                sex2 <- stats::rbinom(1,1,same.sex.sex) + 1
                number1 <- availables[[sex1]][sample(1:activ.selection.size[sex1],1, prob=sample_prob[[sex1]][availables[[sex1]]])]
                number2 <- availables[[sex2]][sample(1:activ.selection.size[sex2],1, prob=sample_prob[[sex2]][availables[[sex2]]])]
                test <- 1
                while(same.sex.selfing==FALSE && number1==number2 && test < 100){
                  number2 <- availables[[sex2]][sample(1:activ.selection.size[sex2],1, prob=sample_prob[[sex2]][availables[[sex2]]])]
                  test <- test+1
                  if(test==100 && number1==number2){
                    warning("Only one remaining individual in the selected cohorts.")
                  }
                }

                if(relative.selection==FALSE){
                  info.father <- best[[sex1]][number1,]
                  # waehle fuers bveite Tier ein Tier aus der Gruppe der Nicht-besten Tiere
                  if(martini.selection){
                    options <- (1:population$info$size[best[[sex2]][1,1],sex2])[-best[[sex2]][,3]]
                    number2 <- sample(options,1)
                    info.mother <- c(best[[sex2]][1,1:2], number2)
                  } else{
                    info.mother <- best[[sex2]][number2,]
                  }
                } else{
                  info.father <- best[[sex1]][sum(cum.sum[[sex1]] <stats::runif(1,0,cum.sum[[sex1]])) +1 ,1:3]
                  info.mother <- best[[sex2]][sum(cum.sum[[sex2]] <stats::runif(1,0,cum.sum[[sex2]])) +1 ,1:3]
                }

              }
            } else{
              sex1 <- sex2 <- stats::rbinom(1,1,selfing.sex)+1
              if(length(availables[[sex1]])==0){
                if(verbose) cat(paste0("No ", if(sex1==1){"male"} else{"female"}, " individuals left. Use other sex individuals\n"))
                sex1 <- sex2 <- 3 - sex1
              }
              number1 <- number2 <- availables[[sex1]][sample(1:activ.selection.size[sex1],1, prob=sample_prob[[sex1]][availables[[sex1]]])]
              info.father <- best[[sex1]][number1,]
              info.mother <- best[[sex2]][number2,]

            }


            accepted <- TRUE
            if(max.mating.pair < Inf){

              if(info.father[2]==info.mother[2]){
                if(info.father[1]>info.mother[1] || (info.father[1]==info.mother[1] & info.father[3]> info.mother[3])){
                  info_mother_temp <- info.mother
                  info.mother <- info.father
                  info.father <- info_mother_temp

                }
              }

              if( sum((colSums(t(info_father_list) == info.father) + colSums(t(info_mother_list) == info.mother))==10) >= max.mating.pair){
                accepted <- FALSE
              }
            }

            check_rel = check.parents(population, info.father, info.mother, max.rel=max_rel_temp)
            max_counter <- max_counter + 1
            if(max_counter>=25){
              if(max_rel<2){
                warning("No remaining possible mating via avoid.mating. Proceed with silbing-mating.\n")
                max_rel_temp <- max_rel + 1
                if(max_counter >= 50){
                  check_rel <- TRUE
                  max_rel_temp <- max_rel + 2
                }
                if(max_counter >= 250){
                  accepted <- TRUE
                  warning("No remaining possible mating via max.mating.pair. Proceed without limitation\n")
                }
              } else{
                if(max_counter >= 1000){
                  accepted <- TRUE
                  warning("A mating was executed more times than allowed in max.mating.pair (check if the total number of allowed mating is sufficient!)\n")
                }
              }

            }





            if(accepted && max.mating.pair< Inf){
              info_father_list[animal.nr,] <- info.father
              info_mother_list[animal.nr,] <- info.mother
            }



          }
          if(martini.selection==FALSE){

            selection.rate[[sex1]][number1] <- selection.rate[[sex1]][number1] + repeat.mating.temp
            selection.rate[[sex2]][number2] <- selection.rate[[sex2]][number2] + repeat.mating.temp

            selection.rate.litter[[sex1]][number1] <- selection.rate.litter[[sex1]][number1] + 1
            selection.rate.litter[[sex2]][number2] <- selection.rate.litter[[sex2]][number2] + 1

            if(selection.rate[[sex1]][number1] >= max.offspring[sex1]){
              activ.selection.size[sex1] <-  activ.selection.size[sex1] -1
              availables[[sex1]] <- availables[[sex1]][availables[[sex1]]!=number1]
            } else if(selection.rate.litter[[sex1]][number1] >= max.litter[sex1]){
              activ.selection.size[sex1] <-  activ.selection.size[sex1] -1
              availables[[sex1]] <- availables[[sex1]][availables[[sex1]]!=number1]
            }
            if(selection.rate[[sex2]][number2] >= max.offspring[sex2]){
              if(sex1!= sex2 || number1 != number2){
                activ.selection.size[sex2] <-  activ.selection.size[sex2] -1
              }
              availables[[sex2]] <- availables[[sex2]][availables[[sex2]]!=number2]
            } else if(selection.rate.litter[[sex2]][number2] >= max.litter[sex2]){
              if(sex1!= sex2 || number1 != number2){
                activ.selection.size[sex2] <-  activ.selection.size[sex2] -1
              }
              availables[[sex2]] <- availables[[sex2]][availables[[sex2]]!=number2]
            }
          }

        }


      }
      if(store.comp.times.generation){
        tack <- as.numeric(Sys.time())
        pre_stuff <- pre_stuff + tack -tick
      }

      father <- population$breeding[[info.father[1]]][[info.father[2]]][[info.father[3]]]
      mother <- population$breeding[[info.mother[1]]][[info.mother[2]]][[info.mother[3]]]

      # Simulation of meiosis
      if(copy.individual){
        info.mother <- info.father
        child1 <- list(father[[1]], father[[3]], father[[5]], father[[7]], father[[11]], 0, if(length(father)>19){father[[19]]} else{0})
        child2 <- list(father[[2]], father[[4]], father[[6]], father[[8]], father[[12]], 0, if(length(father)>19){father[[20]]} else{0})
      } else{
        child1 <- breeding.intern(info.father, father, population,
                                  mutation.rate, remutation.rate, recombination.rate,
                                  recom.f.indicator, duplication.rate, duplication.length,
                                  duplication.recombination, delete.same.origin=delete.same.origin,
                                  gene.editing=(gene.editing.offspring*gene.editing.offspring.sex[1]), nr.edits= nr.edits,
                                  gen.architecture=gen.architecture.m,
                                  decodeOriginsU=decodeOriginsU)

        child2 <- breeding.intern(info.mother, mother, population,
                                  mutation.rate, remutation.rate, recombination.rate,
                                  recom.f.indicator, duplication.rate, duplication.length,
                                  duplication.recombination, delete.same.origin=delete.same.origin,
                                  gene.editing=(gene.editing.offspring * gene.editing.offspring.sex[1]) , nr.edits= nr.edits,
                                  gen.architecture=gen.architecture.f,
                                  decodeOriginsU=decodeOriginsU)
      }
      if(dh.mating){
        if(stats::rbinom(1,1,dh.sex)==0){
          child2 <- child1
        } else{
          child1 <- child2
        }
      }

      # praeimplantation needs to be revisited!!

      if(length(praeimplantation)>0){
        counter <- 1
        good1 <- 0
        n.snps <- sum(population$info$snp)
        while(good1==0){
          hap1 <- rep(0,n.snps)
          temp1 <- 0
          current.animal <- child1
          for(index2 in 1:(length(current.animal[[1]])-1)){
            relevant.snp <- (population$info$snp.position < current.animal[[1]][index2+1])*(population$info$snp.position >= current.animal[[1]][index2])*(1:n.snps)
            ursprung <- decodeOriginsU(current.animal[[3]],index2)
            ursprung[1] <- population$info$origin.gen[ursprung[1]]
            hap1[relevant.snp] <-population$breeding[[ursprung[1]]][[ursprung[2]]][[ursprung[3]]][[ursprung[4]+8]][relevant.snp]
          }
          if(length(current.animal[[2]])>0){
            for(index2 in 1:length(current.animal[[2]])){
              position <- which(population$info$snp.position==current.animal[[2]][index2])
              hap1[position] <- 1-population$breeding[[current.gen+1]][[sex]][[current.size[sex]]][[9+temp1]][position]
            }
          }
          if(hap1[pos] == praeimplantation.max[[sex1]][[number1]] || counter==25){
            good1 <- 1
            if(counter==25){warning("praeimplantation failed.")}
          } else{
            child1 <- breeding.intern(info.father, father, population,
                                      mutation.rate, remutation.rate, recombination.rate,
                                      recom.f.indicator, duplication.rate, duplication.length,
                                      duplication.recombination, delete.same.origin=delete.same.origin,
                                      gene.editing=gene.editing, nr.edits= nr.edits,gen.architecture=gen.architecture.m,
                                      decodeOriginsU=decodeOriginsU)
            counter <- counter +1
          }
        }

        good1 <- 0
        n.snps <- sum(population$info$snp)
        while(good1==0){
          hap1 <- rep(0,n.snps)
          temp1 <- 0
          current.animal <- child2
          for(index2 in 1:(length(current.animal[[1]])-1)){
            relevant.snp <- (population$info$snp.position < current.animal[[1]][index2+1])*(population$info$snp.position >= current.animal[[1]][index2])*(1:n.snps)
            ursprung <-  decodeOriginsU(current.animal[[3]],index2)
            ursprung[1] <- population$info$origin.gen[ursprung[1]]
            hap1[relevant.snp] <-population$breeding[[ursprung[1]]][[ursprung[2]]][[ursprung[3]]][[ursprung[4]+8]][relevant.snp]
          }
          if(length(current.animal[[2]])>0){
            for(index2 in 1:length(current.animal[[2]])){
              position <- which(population$info$snp.position==current.animal[[2]][index2])
              hap1[position] <- 1-population$breeding[[current.gen+1]][[sex]][[current.size[sex]]][[9+temp1]][position]
            }
          }
          if(hap1[pos] == praeimplantation.max[[sex2]][[number2]] || counter==25){
            good1 <- 1
            if(counter==25){if(verbose) cat("praeimplantation failed!")}
          } else{
            child2 <- breeding.intern(info.mother, mother, population,
                                      mutation.rate, remutation.rate, recombination.rate,
                                      recom.f.indicator, duplication.rate, duplication.length,
                                      duplication.recombination, delete.same.origin=delete.same.origin,
                                      gene.editing=gene.editing, nr.edits= nr.edits,
                                      gen.architecture=gen.architecture.f, decodeOriginsU=decodeOriginsU)
            counter <- counter +1
          }
        }
      }

      # Put together offspring

      child_temp <- list(child1[[1]], child2[[1]], child1[[2]], child2[[2]], child1[[3]], child2[[3]], child1[[4]], child2[[4]])


      if(copy.individual){
        population$breeding[[current.gen+1]][[sex+22]][current.size[sex]] <- population$breeding[[info.father[1]]][[info.father[2]+22]][info.father[3]]
      }
      population$info$size[current.gen+1 ,sex] <- population$info$size[current.gen+1,sex] + 1

      if(is.vector(child1[[5]])){
        child_temp[[11]] <- t(as.matrix(child1[[5]]))
      } else{
        child_temp[[11]] <- child1[[5]]
      }
      if(is.vector(child2[[5]])){
        child_temp[[12]] <- t(as.matrix(child2[[5]]))
      } else{
        child_temp[[12]] <- child2[[5]]
      }
      if(save.recombination.history && current.gen==1){
        if(length(child1[[6]][-c(1,length(child1[[6]]))])>0){
          child_temp[[13]] <- cbind(current.gen, child1[[6]][-c(1,length(child1[[6]]))], deparse.level = 0)
        } else{
          child_temp[[13]] <- cbind(0,0, deparse.level = 0)
        }
        if(length( child2[[6]][-c(1,length(child2[[6]]))])>0){
          child_temp[[14]] <- cbind(current.gen, child2[[6]][-c(1,length(child2[[6]]))], deparse.level = 0)
        } else{
          child_temp[[14]] <- cbind(0,0, deparse.level = 0)
        }

      } else if(save.recombination.history && current.gen>1){
        if(length(child1[[6]][-c(1,length(child1[[6]]))])>0){
          child_temp[[13]] <- rbind(population$breeding[[info.father[1]]][[info.father[2]]][[info.father[3]]][[13]], population$breeding[[info.father[1]]][[info.father[2]]][[info.father[3]]][[14]], cbind(current.gen, child1[[6]][-c(1,length(child1[[6]]))], deparse.level = 0), deparse.level = 0)
        } else{
          child_temp[[13]] <- rbind(population$breeding[[info.father[1]]][[info.father[2]]][[info.father[3]]][[13]], population$breeding[[info.father[1]]][[info.father[2]]][[info.father[3]]][[14]], deparse.level = 0)

        }
        if(length( child2[[6]][-c(1,length(child2[[6]]))])>0){
          child_temp[[14]] <- rbind(population$breeding[[info.mother[1]]][[info.mother[2]]][[info.mother[3]]][[13]], population$breeding[[info.mother[1]]][[info.mother[2]]][[info.mother[3]]][[14]], cbind(current.gen, child2[[6]][-c(1,length(child2[[6]]))]))
        } else{
          child_temp[[14]] <- rbind(population$breeding[[info.mother[1]]][[info.mother[2]]][[info.mother[3]]][[13]], population$breeding[[info.mother[1]]][[info.mother[2]]][[info.mother[3]]][[14]], deparse.level = 0)

        }

      } else{
        #population$breeding[[current.gen+1]][[sex]][[current.size[sex]]][[13]] <- "test"
      }

      if(share.phenotyped==1){
        is.obs <- TRUE
      } else if(share.phenotyped==0){
        is.obs <- FALSE
      } else{
        is.obs <- stats::rbinom(1,1, share.phenotyped)==1
      }

      if(phenotyping.child=="obs"){
        child_temp[[15]] <- n.observation * is.obs
      } else if(phenotyping.child=="addobs"){
        child_temp[[15]] <-  colMeans(rbind(population$breeding[[info.mother[1]]][[info.mother[2]]][[info.mother[3]]][[15]],
                                                                                                  population$breeding[[info.father[1]]][[info.father[2]]][[info.father[3]]][[15]]))
        switch <- (child_temp[[15]]< (n.observation*is.obs))
        child_temp[[15]][switch] <- (n.observation*is.obs)[switch]

      } else{
        child_temp[[15]] <- rep(0L, population$info$bv.nr)
      }

      if(copy.individual && !copy.individual.keep.pheno){
        child_temp[[15]] <- rep(0L, population$info$bv.nr)
      }

      child_temp[[27]] <- "placeholder"

      if(copy.individual){
        child_temp[[16]] <- population$breeding[[info.father[1]]][[info.father[2]]][[info.father[3]]][[16]]
        if(length(population$breeding[[info.father[1]]][[info.father[2]]][[info.father[3]]][[22]])>0){
          child_temp[[22]] <- population$breeding[[info.father[1]]][[info.father[2]]][[info.father[3]]][[22]]
        }

        if(added.genotyped>0 && child_temp[[16]]==0){
          if(stats::rbinom(1,1,added.genotyped)==1){
            child_temp[[16]] <- 1
            child_temp[[22]] <- c(child_temp[[22]], genotyped.array)
          }
        }

        first_copy <- population$breeding[[info.father[1]]][[info.father[2]]][[info.father[3]]][[21]]
        new_copy <- rbind(population$breeding[[first_copy[1,1]]][[first_copy[1,2]]][[first_copy[1,3]]][[21]],
                          c(current.gen+1, sex, current.size[sex]), deparse.level = 0)
        storage.mode(new_copy) <- "integer"
        if(nrow(new_copy)>1){
          for(index7 in 1:(nrow(new_copy)-1)){
            population$breeding[[new_copy[index7,1]]][[new_copy[index7,2]]][[new_copy[index7,3]]][[21]] <- new_copy
          }
        }
        child_temp[[21]] <- new_copy

      } else{
        if(share.genotyped==1 || (share.genotyped!=0 && stats::rbinom(1,1,share.genotyped)==1)){
          child_temp[[16]] <- 1
          child_temp[[22]] <- genotyped.array
        } else{
          child_temp[[16]] <- 0
        }

        child_temp[[21]] <- cbind(current.gen+1, sex, current.size[sex], deparse.level = 0)
        storage.mode(child_temp[[21]]) <- "integer"
      }
      if(length(child1[[7]])>0){
        child_temp[[19]] <- child1[[7]]
      }
      if(length(child2[[7]])>0){
        child_temp[[20]] <- child2[[7]]
      }

      if(copy.individual){
        child_temp[[25]] <- population$breeding[[info.father[1]]][[info.father[2]]][[info.father[3]]][[25]]

        if(length(population$breeding[[info.father[1]]][[info.father[2]]][[info.father[3]]][[26]])>0){
          child_temp[[26]] <- population$breeding[[info.father[1]]][[info.father[2]]][[info.father[3]]][[26]]
        }

      } else{
        child_temp[[25]] <- length(bv.ignore.traits)==0
        if(length(temp123)>0){
          child_temp[[26]] <- temp123
        }

      }



      population$breeding[[current.gen+1]][[sex]][[current.size[sex]]] <- child_temp
      if(store.comp.times.generation){
        tock <- as.numeric(Sys.time())
        generation_stuff <- generation_stuff + tock -tack
      }
      if(copy.individual){
        individual.id <- population$breeding[[info.father[1]]][[14+info.father[2]]][info.father[3]]
      }

      # calculate underlying true genomic value
      ## revisit computations for non-QTL based effects for efficiently
      if(population$info$bve){
        activ_bv <- population$info$bv.random.activ
        if(length(activ_bv)>0){
          if(!copy.individual || store.effect.freq){
            temp_out <- calculate.bv(population, current.gen+1, sex, current.size[sex], activ_bv, import.position.calculation=import.position.calculation, decodeOriginsU=decodeOriginsU, store.effect.freq=store.effect.freq, bit.storing=bit.storing, nbits=nbits, output_compressed=FALSE, bv.ignore.traits=bv.ignore.traits)
            new.bv[activ_bv] <- temp_out[[1]]

            if(store.effect.freq){
              if(length(population$info$store.effect.freq) < (current.gen+1) || length(population$info$store.effect.freq[[current.gen+1]])==0){
                colnames(temp_out[[2]]) <- c("Homo0", "Hetero", "Homo1")
                rownames(temp_out[[2]]) <- population$info$snp.name[population$info$effect.p]
                population$info$store.effect.freq[[current.gen+1]] <- temp_out[[2]]
              } else{
                population$info$store.effect.freq[[current.gen+1]] <- population$info$store.effect.freq[[current.gen+1]] + temp_out[[2]]
              }
            }

          } else{
            activ_indi <- info.father
            new.bv[activ_bv] <- population$breeding[[activ_indi[1]]][[activ_indi[2]+6]][activ_bv, activ_indi[3]]
          }


        }

        if(population$info$bv.calc > 0  && population$info$bv.random[population$info$bv.calc] && prod(population$info$is.combi[population$info$bv.calc:population$info$bv.nr])==0){

          #Means passt (Korrelation exakt wie gewuenscht)
          means <- 0.5*(population$breeding[[info.father[1]]][[6+info.father[2]]][population$info$bv.calc:population$info$bv.nr,info.father[3]] + population$breeding[[info.mother[1]]][[6+info.mother[2]]][population$info$bv.calc:population$info$bv.nr,info.mother[3]])

          # Berechnung i (0-0.5)
          varp <- kinship.emp(list(population$breeding[[info.father[1]]][[info.father[2]]][[info.father[3]]], population$breeding[[info.mother[1]]][[info.mother[2]]][[info.mother[3]]]))
          varp <- (2 - (2* (varp[1,1]-0.5) + 2 * (varp[2,2]- 0.5)))/4

          if(FALSE){
            new.bv[population$info$bv.calc:population$info$bv.nr] <- stats::rnorm(1, mean=means, sd= sqrt(var*population$info$bv.random.variance[bven]))
          } else{
            if(population$info$bv.calc==1){
              population$info$current.bv.random.variance <- varp * population$info$bv.random.variance
              bv.var <- diag.mobps(sqrt(population$info$current.bv.random.variance)) %*%population$info$current.bv.correlation %*% diag.mobps(sqrt(population$info$current.bv.random.variance))
              single.mean <- means
            } else{
              #population$info$current.bv.random.variance <- c(population$info$bv.random.variance[1:(population$info$bv.calc-1)], population$info$bv.random.variance[population$info$bv.calc:population$info$bv.nr])
              population$info$current.bv.random.variance <- c(population$info$bv.random.variance[1:(population$info$bv.calc-1)],varp * population$info$bv.random.variance[population$info$bv.calc:population$info$bv.nr])

              AA <- diag.mobps(sqrt(population$info$current.bv.random.variance)[1:(population$info$bv.calc-1)]) %*% population$info$current.bv.correlation[1:(population$info$bv.calc-1), 1:(population$info$bv.calc-1)]%*% diag.mobps(sqrt(population$info$current.bv.random.variance)[(1:(population$info$bv.calc-1))])
              BB <- diag.mobps(sqrt(population$info$current.bv.random.variance)[1:(population$info$bv.calc-1)]) %*%population$info$current.bv.correlation[1:(population$info$bv.calc-1), -(1:(population$info$bv.calc-1))]%*% diag.mobps(sqrt(population$info$current.bv.random.variance)[-(1:(population$info$bv.calc-1))])
              CC <- diag.mobps(sqrt(population$info$current.bv.random.variance)[-(1:(population$info$bv.calc-1))]) %*%population$info$current.bv.correlation[-(1:(population$info$bv.calc-1)), -(1:(population$info$bv.calc-1))] %*% diag.mobps(sqrt(population$info$current.bv.random.variance)[-(1:(population$info$bv.calc-1))])
              if (requireNamespace("MASS", quietly = TRUE)) {
                bv.var <- CC - t(BB) %*% MASS::ginv(AA) %*% BB
                single.mean <- means + t(BB) %*% MASS::ginv(AA) %*% ( new.bv[1:(population$info$bv.calc-1)]-mu1[1:(population$info$bv.calc-1)])
              } else{
                bv.var <- CC - t(BB) %*% solve(AA) %*% BB
                single.mean <- means + t(BB) %*% solve(AA) %*% ( new.bv[1:(population$info$bv.calc-1)]-mu1[1:(population$info$bv.calc-1)])
              }


            }

            bv.var_chol <- t(chol(bv.var))
            population$info$bv.correlation_col <- bv.var_chol

#            random_part <- varp * bv.var_chol %*% stats:rnorm(population$info$bv.nr - population$info$bv.calc +1, 0,1 )
            random_part <- bv.var_chol %*% stats::rnorm(population$info$bv.nr - population$info$bv.calc +1, 0,1 )
            means_part <- single.mean


            new.bv[population$info$bv.calc:population$info$bv.nr] <-   single.mean  + random_part
            #new.bv[population$info$bv.calc:population$info$bv.nr] <- bv.var_chol %*% stats::rnorm(population$info$bv.nr - population$info$bv.calc +1, 0,1 ) + single.mean
          }


        }
        if(sum(population$info$is.combi)>0){
          combis <- which(population$info$is.combi)
          for(combi in combis){
            new.bv[combi] <- sum(new.bv * population$info$combi.weights[[combi]])
          }
        }

        if(copy.individual && copy.individual.keep.bve){
          new.bve <- population$breeding[[info.father[1]]][[2+info.father[2]]][,info.father[3]]
          new.reli <- population$breeding[[info.father[1]]][[18+info.father[2]]][,info.father[3]]
        }

        if(copy.individual && length(population$breeding[[info.father[1]]][[info.father[2]]][[info.father[3]]][[23]])>0){
          population$breeding[[current.gen+1]][[sex]][[current.size[sex]]][[23]] <-
            population$breeding[[info.father[1]]][[info.father[2]]][[info.father[3]]][[23]]
        }
        if(phenotyping.child=="mean" || phenotyping.child=="addobs"){
          if(copy.individual){
            if(length(population$breeding[[info.father[1]]][[info.father[2]]][[info.father[3]]][[24]])>0){
              population$breeding[[current.gen+1]][[sex]][[current.size[sex]]][[24]] <-
                population$breeding[[info.father[1]]][[info.father[2]]][[info.father[3]]][[24]]
            }
          } else{
            new.bv_approx <- 0.5 * (population$breeding[[info.father[1]]][[8+info.father[2]]][,info.father[3]] + population$breeding[[info.mother[1]]][[8+info.mother[2]]][,info.mother[3]])
          }

        }

        # phenotyping

        if(phenotyping.child=="obs" || phenotyping.child=="addobs" || copy.individual){
          if(sum(n.observation)>0 || (copy.individual && copy.individual.keep.pheno)){

            if( length(population$breeding[[current.gen+1]][[sex]][[current.size[sex]]][[23]]) == 0){
              population$breeding[[current.gen+1]][[sex]][[current.size[sex]]][[23]] <- stats::rnorm(population$info$bv.nr,0,1)
            }

            obsmax <- max(population$breeding[[current.gen+1]][[sex]][[current.size[sex]]][[15]])

            if(length(population$breeding[[current.gen+1]][[sex]][[current.size[sex]]][[24]])==0 ||
               obsmax > ncol(population$breeding[[current.gen+1]][[sex]][[current.size[sex]]][[24]]) &&
               obsmax > 0){

              population$breeding[[current.gen+1]][[sex]][[current.size[sex]]][[24]] <- cbind(population$breeding[[current.gen+1]][[sex]][[current.size[sex]]][[24]],
                                                                            matrix(stats::rnorm(obsmax * population$info$bv.nr - length(population$breeding[[current.gen+1]][[sex]][[current.size[sex]]][[24]]),0,1),
                                                                                   nrow = population$info$bv.nr))
            }

            for(bven in setdiff(1:population$info$bv.nr, activ.trafo)){
              if(population$breeding[[current.gen+1]][[sex]][[current.size[sex]]][[15]][bven]>=1){

                if(copy.individual && n.observation[bven]==0){
                  new.bv_approx[bven] <- get.pheno(population, database = matrix(info.father[c(1,2,3,3)], nrow=1))[bven]
                } else{
                  new.bv_approx[bven] <- (sqrt(sigma.e.rest) * rowMeans(population$info$pheno.correlation %*% population$breeding[[current.gen+1]][[sex]][[current.size[sex]]][[24]][,1:population$breeding[[current.gen+1]][[sex]][[current.size[sex]]][[15]][bven]]) +
                                            sqrt(sigma.e.perm) * population$info$pheno.correlation %*% population$breeding[[current.gen+1]][[sex]][[current.size[sex]]][[23]] +
                                            new.bv)[bven]
                }


              }

            }

            for(bven in intersect(1:population$info$bv.nr, activ.trafo)){
              if(population$breeding[[current.gen+1]][[sex]][[current.size[sex]]][[15]][bven]>=1){
                new_pheno <- (sqrt(sigma.e.rest) * population$info$pheno.correlation %*% population$breeding[[current.gen+1]][[sex]][[current.size[sex]]][[24]][,1:population$breeding[[current.gen+1]][[sex]][[current.size[sex]]][[15]][bven]])[bven,] +
                  (sqrt(sigma.e.perm) * population$info$pheno.correlation %*% population$breeding[[current.gen+1]][[sex]][[current.size[sex]]][[23]])[bven] +
                  new.bv[bven]

                if(copy.individual && n.observation[bven]==0){
                  new.bv_approx[bven] <- get.pheno(population, database = matrix(info.father[c(1,2,3,3)], nrow=1))[bven]
                } else{
                  new.bv_approx[bven] <- mean(population$info$phenotypic.transform.function[[bven]](new_pheno))
                }

              }
            }

          }


        }
      }

      population$breeding[[current.gen+1]][[2+sex]][,current.size[sex]] <- new.bve
      population$breeding[[current.gen+1]][[6+sex]][,current.size[sex]] <- new.bv
      population$breeding[[current.gen+1]][[8+sex]][,current.size[sex]] <- new.bv_approx
      population$breeding[[current.gen+1]][[18+sex]][,current.size[sex]] <- new.reli


      if(copy.individual){
        population$breeding[[current.gen+1]][[14+sex]][current.size[sex]] <- individual.id
      }
      #    } else if(length(population$breeding[[current.gen+1]][[sex+2]])==0){
      #      population$breeding[[current.gen+1]][[2+sex]] <- rep(0, breeding.size[sex])
      #      population$breeding[[current.gen+1]][[4+sex]] <- rep(new.class, breeding.size[sex])
      #     population$breeding[[current.gen+1]][[6+sex]] <- new.bv[sex,]
      #      population$breeding[[current.gen+1]][[8+sex]] <- new.bv.approx[sex,]

      if(store.comp.times.generation){
        tock2 <- as.numeric(Sys.time())
        bv_stuff <- bv_stuff+tock2-tock
      }
      if(display.progress & verbose & (breeding.size.total < 100 || animal.nr%%pb_temp==0)){
        utils::setTxtProgressBar(pb, animal.nr)
      }

      current.size[sex] <- current.size[sex] +1
    }
    if(display.progress & verbose){
      close(pb)
    }

  }

}

  #######################################################################
  ########################## Final housekeeping #########################
  #######################################################################
{
  delete.haplotypes <- delete.haplotypes[delete.haplotypes>1]
  if(length(delete.haplotypes)>0){
    for(gen_r in delete.haplotypes){
      for(sex in 1:2){
        n.animals <- length(population$breeding[[gen_r]][[sex]])
        if(n.animals>0){
          for(index2 in 1:n.animals){
            population$breeding[[gen_r]][[sex]][[index2]][[9]] <- "removed"
            population$breeding[[gen_r]][[sex]][[index2]][[10]] <- "removed"
          }
        }
      }
    }
  }
  delete.individuals <- delete.individuals[delete.individuals>1]
  if(length(delete.individuals)>0){
    for(gen_r in delete.individuals){
      for(sex in delete.sex){
        n.animals <- length(population$breeding[[gen_r]][[sex]])
        if(n.animals>0){
          population$breeding[[gen_r]][[sex]]<- "removed"

        }
      }
    }
  }
  if(store.breeding.totals){
    cur <- length(population$info$breeding.totals) + 1L
    population$info$breeding.totals[[cur]] <- list()
    population$info$breeding.totals[[cur]][[1]] <- current.gen +1L
    population$info$breeding.totals[[cur]][[2]] <- breeding.size
    population$info$breeding.totals[[cur]][[3]] <- best[[1]]
    population$info$breeding.totals[[cur]][[4]] <- best[[2]]
    population$info$breeding.totals[[cur]][[5]] <- selection.rate[[1]]
    population$info$breeding.totals[[cur]][[6]] <- selection.rate[[2]]
    population$info$breeding.totals[[cur]][[7]] <- chosen.animals.list
  }

  if(store.bve.data && bve){
    cur <- length(population$info$bve.data) + 1
    population$info$bve.data[[cur]] <- list()
    population$info$bve.data[[cur]][[1]] <- current.gen
    population$info$bve.data[[cur]][[2]] <- sigma.e
    population$info$bve.data[[cur]][[3]] <- sigma.e.hat
    population$info$bve.data[[cur]][[4]] <- sigma.g
    population$info$bve.data[[cur]][[5]] <- sigma.a.hat
    population$info$bve.data[[cur]][[6]] <- bve.database
    population$info$bve.data[[cur]][[7]] <- phenotyping
    population$info$bve.data[[cur]][[8]] <- y_real
    population$info$bve.data[[cur]][[9]] <- y_hat
    population$info$bve.data[[cur]][[10]] <- y
  }

  population$info$last.sigma.e.database <- sigma.e.database
  population$info$last.sigma.e.value <- sigma.e
  if(length(heritability)>0){
    population$info$last.sigma.e.heritability <- heritability
  }


  if(repeat.mating.overwrite){
    population$info$repeat.mating <- repeat.mating.store
    population$info$repeat.mating.copy <- repeat.mating.copy.store
  }

  if(store.comp.times){
    comp.times[7] <- as.numeric(Sys.time())
    comp.times <- c(comp.times[-1] - comp.times[-length(comp.times)], comp.times[length(comp.times)]-comp.times[1])

    comp.times[comp.times<0] <- 0
    comp.times[comp.times>10e6] <- 0

    population$info$comp.times <- round(rbind(population$info$comp.times, comp.times, deparse.level = 0), digits=4)
    if(nrow(population$info$comp.times)==1){
      colnames(population$info$comp.times) <- c("preparation", "new real BV", "phenotypes", "BVE","selection","generate new individuals","total")
    }
  }
  if(store.comp.times.bve){
    comp.times.bve <- c(comp.times.bve[-1] - comp.times.bve[-length(comp.times.bve)], zcalc, z_chol, z_uhat, z_ped, z_h, comp.times.bve[length(comp.times.bve)]-comp.times.bve[1])
    comp.times.bve[comp.times.bve<0] <- 0
    comp.times.bve[comp.times.bve>10e6] <- 0

    population$info$comp.times.bve <- round(rbind(population$info$comp.times.bve, comp.times.bve, deparse.level = 0), digits=4)
    if(nrow(population$info$comp.times.bve)==1){
      colnames(population$info$comp.times.bve) <- c("y_z_import", "A genomic", "solveMixed","Gwas_stuff", "Derive Z", "A inversion", "rrBlup", "A-Pedigree","SingleStep H", "Total")
    }
  }
  if(store.comp.times.generation){
    comp.times.generation <- c(pre_stuff, generation_stuff, bv_stuff, sum(pre_stuff, generation_stuff, bv_stuff))
    population$info$comp.times.generation <- round(rbind(population$info$comp.times.generation, comp.times.generation, deparse.level = 0), digits=4)
    if(nrow(population$info$comp.times.generation)==1){
      colnames(population$info$comp.times.generation) <- c("Preparation", "Generation", "BV-Calculation", "Total")
    }
  }

  if(length(name.cohort)>0){
    if(breeding.size[1] > 0 && breeding.size[2]>0){
      population$info$cohorts <- rbind(population$info$cohorts, c(paste0(name.cohort, "_M"), current.gen+1, breeding.size[1],0, new.class, (current.size-breeding.size)[1], 0,
                                                                  time.point, creating.type),
                                       c(paste0(name.cohort, "_F"), current.gen+1, 0, breeding.size[2], new.class, 0, (current.size-breeding.size)[2],time.point, creating.type))
      if(verbose) cat("Added _M, _F to cohort names!\n")
      rownames(population$info$cohorts)[(nrow(population$info$cohorts)-1):nrow(population$info$cohorts)] <- paste0(name.cohort, c("_M", c("_F")))

      if(verbose){
        posi <- get.database(population, cohorts = paste0(name.cohort, "_M"))
        cat(paste0("Successfully generated cohort: ",  paste0(name.cohort, "_M"), "\n",
                   "Database position: ", posi[1], " (gen), ", posi[2], " (sex), ", posi[3], " (first), ", posi[4], " (last).\n" ))
        posi <- get.database(population, cohorts = paste0(name.cohort, "_F"))
        cat(paste0("Successfully generated cohort: ",  paste0(name.cohort, "_F"), "\n",
                   "Database position: ", posi[1], " (gen), ", posi[2], " (sex), ", posi[3], " (first), ", posi[4], " (last).\n" ))
      }

    } else{
      population$info$cohorts <- rbind(population$info$cohorts, c(name.cohort, current.gen+1, breeding.size[1:2], new.class, current.size-breeding.size,
                                                                  time.point, creating.type))
      rownames(population$info$cohorts)[nrow(population$info$cohorts)] <- paste0(name.cohort)


      if(verbose){
        posi <- get.database(population, cohorts = name.cohort)
        cat(paste0("Successfully generated cohort: ", name.cohort, "\n",
                   "Database position: ", posi[1], " (gen), ", posi[2], " (sex), ", posi[3], " (first), ", posi[4], " (last).\n" ))
      }
    }

    if(nrow(population$info$cohorts)<=2){
      colnames(population$info$cohorts) <- c("name","generation", "male individuals", "female individuals", "class", "position first male", "position first female",
                                             "time point", "creating.type")    }
  }
  if(Rprof){
    Rprof(NULL)
    population$info$Rprof[[length(population$info$Rprof)+1]] <- utils::summaryRprof()
  }
  }

  return(population)
}

Try the MoBPS package in your browser

Any scripts or data that you put into this service are public.

MoBPS documentation built on Nov. 9, 2021, 5:08 p.m.