diff --git a/DESCRIPTION b/DESCRIPTION
index b2f4c434..757b7cfe 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: AlphaSimR
 Type: Package
 Title: Breeding Program Simulations
-Version: 1.0.4
-Date: 2021-09-08
+Version: 1.1.2
+Date: 2022-2-17
 Authors@R: c(person("Chris", "Gaynor", email = "gaynor.robert@hotmail.com",
   role = c("aut", "cre"), comment = c(ORCID = "0000-0003-0558-6656")),
   person("Gregor", "Gorjanc", role = "ctb", 
@@ -30,7 +30,7 @@ Encoding: UTF-8
 Depends: R (>= 3.3.0), methods, R6
 Imports: Rcpp (>= 0.12.7)
 LinkingTo: Rcpp, RcppArmadillo (>= 0.7.500.0.0), BH
-RoxygenNote: 7.1.1
+RoxygenNote: 7.1.2
 Suggests: knitr, rmarkdown, testthat
 VignetteBuilder: knitr
 NeedsCompilation: true
diff --git a/NAMESPACE b/NAMESPACE
index 8f0a8a40..9ad4614c 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -26,6 +26,7 @@ export(genicVarA)
 export(genicVarAA)
 export(genicVarD)
 export(genicVarG)
+export(getNumThreads)
 export(getQtlMap)
 export(getSnpMap)
 export(gv)
@@ -73,6 +74,15 @@ export(setPheno)
 export(setPhenoGCA)
 export(setPhenoProgTest)
 export(smithHazel)
+export(solveMKM)
+export(solveMVM)
+export(solveRRBLUP)
+export(solveRRBLUPMK)
+export(solveRRBLUPMV)
+export(solveRRBLUP_EM)
+export(solveRRBLUP_EM2)
+export(solveRRBLUP_EM3)
+export(solveUVM)
 export(usefulness)
 export(varA)
 export(varAA)
diff --git a/NEWS b/NEWS
index 13a4ef0e..a0c83d6e 100644
--- a/NEWS
+++ b/NEWS
@@ -1,3 +1,26 @@
+Changes in version 1.1.2
+    
+  Bug fixes
+    -added missing #ifdef _OPENMP to OCS.cpp
+
+Changes in version 1.1.1
+    
+  Bug fixes
+    -removed use of PI variable in C++ code which is compilier specific
+
+Changes in version 1.1.0
+  
+  New features
+    -added snpChip argument to pullIbdHaplo for backwards compatibility
+    -exposed internal mixed model solvers
+    -all selection functions now return a warning when there are not enough individuals
+    
+  Bug fixes
+    -fixed error in pullIbdHaplo when chr isn't NULL
+    -fixed an error with assigning 1 QTL and/or SNP
+    -changed geno slot from matrix to list to support future RcppArmadillo changes
+    -doubleGenome and reduceGenome now work with IBD tracking
+
 Changes in version 1.0.4
 
   Bug fixes
diff --git a/R/Class-HybridPop.R b/R/Class-HybridPop.R
index 04613ee0..43d87b17 100644
--- a/R/Class-HybridPop.R
+++ b/R/Class-HybridPop.R
@@ -122,3 +122,8 @@ setMethod("c",
             return(x)
           }
 )
+
+isHybridPop = function(x) {
+  ret = is(x, class2 = "HybridPop")
+  return(ret)
+}
\ No newline at end of file
diff --git a/R/Class-LociMap.R b/R/Class-LociMap.R
index 9ab69d52..bbd3cd81 100644
--- a/R/Class-LociMap.R
+++ b/R/Class-LociMap.R
@@ -31,6 +31,11 @@ setValidity("LociMap",function(object){
   }
 })
 
+isLociMap = function(x) {
+  ret = is(x, class2 = "LociMap")
+  return(ret)
+}
+
 #TraitA----
 #' @title Additive trait
 #' 
@@ -58,6 +63,11 @@ setValidity("TraitA",function(object){
   }
 })
 
+isTraitA = function(x) {
+  ret = is(x, class2 = "TraitA")
+  return(ret)
+}
+
 #TraitA2----
 #' @title Sex specific additive trait
 #' 
@@ -84,6 +94,11 @@ setValidity("TraitA2",function(object){
   }
 })
 
+isTraitA2 = function(x) {
+  ret = is(x, class2 = "TraitA2")
+  return(ret)
+}
+
 #TraitAE----
 #' @title Additive and epistatic trait
 #' 
@@ -112,6 +127,11 @@ setValidity("TraitAE",function(object){
   }
 })
 
+isTraitAE = function(x) {
+  ret = is(x, class2 = "TraitAE")
+  return(ret)
+}
+
 #TraitAD----
 #' @title Additive and dominance trait
 #' 
@@ -137,6 +157,11 @@ setValidity("TraitAD",function(object){
   }
 })
 
+isTraitAD = function(x) {
+  ret = is(x, class2 = "TraitAD")
+  return(ret)
+}
+
 #TraitA2D----
 #' @title Sex specific additive and dominance trait
 #' 
@@ -162,6 +187,11 @@ setValidity("TraitA2D",function(object){
   }
 })
 
+isTraitA2D = function(x) {
+  ret = is(x, class2 = "TraitA2D")
+  return(ret)
+}
+
 #TraitADE----
 #' @title Additive, dominance, and epistatic trait
 #' 
@@ -190,6 +220,11 @@ setValidity("TraitADE",function(object){
   }
 })
 
+isTraitADE = function(x) {
+  ret = is(x, class2 = "TraitADE")
+  return(ret)
+}
+
 #TraitAG----
 #' @title Additive and GxE trait
 #' 
@@ -219,6 +254,11 @@ setValidity("TraitAG",function(object){
   }
 })
 
+isTraitAG = function(x) {
+  ret = is(x, class2 = "TraitAG")
+  return(ret)
+}
+
 #TraitAEG----
 #' @title Additive, epistasis and GxE trait
 #' 
@@ -248,6 +288,11 @@ setValidity("TraitAEG",function(object){
   }
 })
 
+isTraitAEG = function(x) {
+  ret = is(x, class2 = "TraitAEG")
+  return(ret)
+}
+
 #TraitADG----
 #' @title Additive, dominance and GxE trait
 #' 
@@ -277,6 +322,11 @@ setValidity("TraitADG",function(object){
   }
 })
 
+isTraitADG = function(x) {
+  ret = is(x, class2 = "TraitADG")
+  return(ret)
+}
+
 #TraitADEG----
 #' @title Additive, dominance, epistasis, and GxE trait
 #' 
@@ -306,3 +356,7 @@ setValidity("TraitADEG",function(object){
   }
 })
 
+isTraitADEG = function(x) {
+  ret = is(x, class2 = "TraitADEG")
+  return(ret)
+}
diff --git a/R/Class-Pop.R b/R/Class-Pop.R
index 8097cd01..8b28f102 100644
--- a/R/Class-Pop.R
+++ b/R/Class-Pop.R
@@ -14,10 +14,9 @@
 #' @slot nChr number of chromosomes
 #' @slot ploidy level of ploidy
 #' @slot nLoci number of loci per chromosome
-#' @slot geno "matrix" containing chromosome genotypes. The "matrix" 
-#' has dimensions nChr by 1 and each element is a three dimensional
-#' array of raw values. The array dimensions are nLoci by ploidy by nInd.
-#' 
+#' @slot geno list of nChr length containing chromosome genotypes.
+#' Each element is a three dimensional array of raw values. 
+#' The array dimensions are nLoci by ploidy by nInd.
 #' 
 #' @export
 setClass("RawPop",
@@ -25,7 +24,7 @@ setClass("RawPop",
                  nChr="integer",
                  ploidy="integer",
                  nLoci="integer",
-                 geno="matrix"))
+                 geno="list"))
 
 setValidity("RawPop",function(object){
   errors = character()
@@ -112,6 +111,11 @@ setMethod("show",
           }
 )
 
+isRawPop = function(x) {
+  ret = is(x, class2 = "RawPop")
+  return(ret)
+}
+
 # MapPop ------------------------------------------------------------------
 
 #' @title Raw population with genetic map
@@ -126,7 +130,7 @@ setMethod("show",
 #' @param i index of individuals
 #' @param ... additional 'MapPop' objects
 #' 
-#' @slot genMap "matrix" of chromosome genetic maps
+#' @slot genMap list of chromosome genetic maps
 #' @slot centromere vector of centromere positions
 #' @slot inbred indicates whether the individuals are fully inbred
 #' 
@@ -192,6 +196,11 @@ setMethod("c",
           }
 )
 
+isMapPop = function(x) {
+  ret = is(x, class2 = "MapPop")
+  return(ret)
+}
+
 # NamedMapPop ------------------------------------------------------------------
 
 #' @title Raw population with genetic map and id
@@ -315,7 +324,7 @@ cChr = function(...){
         stopifnot(x@nInd==y@nInd,
                   x@ploidy==y@ploidy)
         x@nChr = x@nChr+y@nChr
-        x@geno = rbind(x@geno,y@geno)
+        x@geno = c(x@geno,y@geno)
         x@genMap = c(x@genMap,y@genMap)
         x@centromere = c(x@centromere,y@centromere)
         x@nLoci = c(x@nLoci,y@nLoci)
@@ -326,6 +335,11 @@ cChr = function(...){
   return(x)
 }
 
+isNamedMapPop = function(x) {
+  ret = is(x, class2 = "NamedMapPop")
+  return(ret)
+}
+
 # Pop ---------------------------------------------------------------------
 
 #' @title Population
@@ -726,6 +740,11 @@ resetPop = function(pop,simParam=NULL){
   return(pop)
 }
 
+isPop = function(x) {
+  ret = is(x, class2 = "Pop")
+  return(ret)
+}
+
 # MegaPop ------------------------------------------------------------------
 
 #' @title Mega-Population
@@ -833,3 +852,7 @@ newMegaPop = function(...){
   return(output)
 }
 
+isMegaPop = function(x) {
+  ret = is(x, class2 = "MegaPop")
+  return(ret)
+}
diff --git a/R/Class-RRsol.R b/R/Class-RRsol.R
index 46ddb793..c8855674 100644
--- a/R/Class-RRsol.R
+++ b/R/Class-RRsol.R
@@ -20,3 +20,7 @@ setClass("RRsol",
                  Vu="matrix",
                  Ve="matrix"))
 
+isRRsol = function(x) {
+  ret = is(x, class2 = "RRsol")
+  return(ret)
+}
diff --git a/R/Class-SimParam.R b/R/Class-SimParam.R
index 1eb17814..704a978f 100644
--- a/R/Class-SimParam.R
+++ b/R/Class-SimParam.R
@@ -1388,6 +1388,7 @@ SimParam = R6Class(
     #' positions. If NULL, the centromere are assumed to 
     #' be metacentric.
     switchGenMap = function(genMap, centromere=NULL){
+      genMap = lapply(genMap, function(x) x-x[1]) # Set position 1 to 0
       if(is.null(centromere)){
         centromere=sapply(genMap,max)/2
       }
@@ -1413,6 +1414,7 @@ SimParam = R6Class(
     #' positions. If NULL, the centromere are assumed to 
     #' be metacentric.
     switchFemaleMap = function(genMap, centromere=NULL){
+      genMap = lapply(genMap, function(x) x-x[1]) # Set position 1 to 0
       if(is.null(centromere)){
         centromere=sapply(genMap,max)/2
       }
@@ -1443,6 +1445,7 @@ SimParam = R6Class(
     #' positions. If NULL, the centromere are assumed to 
     #' be metacentric.
     switchMaleMap = function(genMap, centromere=NULL){
+      genMap = lapply(genMap, function(x) x-x[1]) # Set position 1 to 0
       if(is.null(centromere)){
         centromere=sapply(genMap,max)/2
       }
@@ -1682,7 +1685,11 @@ SimParam = R6Class(
         if(nSitesPerChr[x]==0){
           return(NULL)
         }else{
-          tmp = sort(sample(pot[[x]],nSitesPerChr[x]))
+          if(length(pot[[x]])==1){
+            tmp = pot[[x]]
+          }else{
+            tmp = sort(sample(pot[[x]],nSitesPerChr[x]))
+          }
           # Add site restrictions
           if(private$.restrSites){
             if(QTL){
@@ -1978,3 +1985,8 @@ sampEpiEff = function(qtlLoci,nTraits,corr,gamma,shape,relVar){
   epiEff = sweep(epiEff,2,sqrt(relVar),"*")
   return(epiEff)
 }
+
+isSimParam = function(x) {
+  ret = is(x, class2 = "SimParam")
+  return(ret)
+}
diff --git a/R/OCS.R b/R/OCS.R
new file mode 100644
index 00000000..fa77d4fb
--- /dev/null
+++ b/R/OCS.R
@@ -0,0 +1,4 @@
+
+# OCS_plant = function(pop, gsModel, simParam=NULL){
+#   
+# }
diff --git a/R/RcppExports.R b/R/RcppExports.R
index bdff8df9..2c36e63d 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -1,6 +1,120 @@
 # Generated by using Rcpp::compileAttributes() -> do not edit by hand
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
+#' @title Solve RR-BLUP
+#'
+#' @description
+#' Solves a univariate mixed model of form \eqn{y=X\beta+Mu+e}
+#'
+#' @param y a matrix with n rows and 1 column
+#' @param X a matrix with n rows and x columns
+#' @param M a matrix with n rows and m columns
+#'
+#' @export
+solveRRBLUP <- function(y, X, M) {
+    .Call(`_AlphaSimR_solveRRBLUP`, y, X, M)
+}
+
+#' @title Solve Multivariate RR-BLUP
+#'
+#' @description
+#' Solves a multivariate mixed model of form \eqn{Y=X\beta+Mu+e}
+#'
+#' @param Y a matrix with n rows and q columns
+#' @param X a matrix with n rows and x columns
+#' @param M a matrix with n rows and m columns
+#' @param tol tolerance for convergence
+#' @param maxIter maximum number of iteration
+#'
+#' @export
+solveRRBLUPMV <- function(Y, X, M, maxIter = 1000L, tol = 1e-6) {
+    .Call(`_AlphaSimR_solveRRBLUPMV`, Y, X, M, maxIter, tol)
+}
+
+#' @title Solve Multikernel RR-BLUP
+#'
+#' @description
+#' Solves a univariate mixed model with multiple random effects.
+#'
+#' @param y a matrix with n rows and 1 column
+#' @param X a matrix with n rows and x columns
+#' @param Mlist a list of M matrices
+#' @param maxIter maximum number of iteration
+#'
+#' @export
+solveRRBLUPMK <- function(y, X, Mlist, maxIter = 40L) {
+    .Call(`_AlphaSimR_solveRRBLUPMK`, y, X, Mlist, maxIter)
+}
+
+#' @title Solve RR-BLUP with EM
+#'
+#' @description
+#' Solves a univariate mixed model of form \eqn{y=X\beta+Mu+e} using
+#' the Expectation-Maximization algorithm.
+#'
+#' @param Y a matrix with n rows and 1 column
+#' @param X a matrix with n rows and x columns
+#' @param M a matrix with n rows and m columns
+#' @param Vu initial guess for variance of marker effects
+#' @param Ve initial guess for error variance
+#' @param tol tolerance for declaring convergence
+#' @param maxIter maximum iteration for attempting convergence
+#' @param useEM should EM algorithm be used. If false, no estimation of
+#' variance components is performed. The initial values are treated as true.
+#'
+#' @export
+solveRRBLUP_EM <- function(Y, X, M, Vu, Ve, tol, maxIter, useEM) {
+    .Call(`_AlphaSimR_solveRRBLUP_EM`, Y, X, M, Vu, Ve, tol, maxIter, useEM)
+}
+
+#' @title Solve RR-BLUP with EM and 2 random effects
+#'
+#' @description
+#' Solves a univariate mixed model of form \eqn{y=X\beta+M_1u_1+M_2u_2+e} using
+#' the Expectation-Maximization algorithm.
+#'
+#' @param Y a matrix with n rows and 1 column
+#' @param X a matrix with n rows and x columns
+#' @param M1 a matrix with n rows and m1 columns
+#' @param M2 a matrix with n rows and m2 columns
+#' @param Vu1 initial guess for variance of the first marker effects
+#' @param Vu2 initial guess for variance of the second marker effects
+#' @param Ve initial guess for error variance
+#' @param tol tolerance for declaring convergence
+#' @param maxIter maximum iteration for attempting convergence
+#' @param useEM should EM algorithm be used. If false, no estimation of
+#' variance components is performed. The initial values are treated as true.
+#'
+#' @export
+solveRRBLUP_EM2 <- function(Y, X, M1, M2, Vu1, Vu2, Ve, tol, maxIter, useEM) {
+    .Call(`_AlphaSimR_solveRRBLUP_EM2`, Y, X, M1, M2, Vu1, Vu2, Ve, tol, maxIter, useEM)
+}
+
+#' @title Solve RR-BLUP with EM and 3 random effects
+#'
+#' @description
+#' Solves a univariate mixed model of form \eqn{y=X\beta+M_1u_1+M_2u_2+M_3u_3+e} using
+#' the Expectation-Maximization algorithm.
+#'
+#' @param Y a matrix with n rows and 1 column
+#' @param X a matrix with n rows and x columns
+#' @param M1 a matrix with n rows and m1 columns
+#' @param M2 a matrix with n rows and m2 columns
+#' @param M3 a matrix with n rows and m3 columns
+#' @param Vu1 initial guess for variance of the first marker effects
+#' @param Vu2 initial guess for variance of the second marker effects
+#' @param Vu3 initial guess for variance of the second marker effects
+#' @param Ve initial guess for error variance
+#' @param tol tolerance for declaring convergence
+#' @param maxIter maximum iteration for attempting convergence
+#' @param useEM should EM algorithm be used. If false, no estimation of
+#' variance components is performed. The initial values are treated as true.
+#'
+#' @export
+solveRRBLUP_EM3 <- function(Y, X, M1, M2, M3, Vu1, Vu2, Vu3, Ve, tol, maxIter, useEM) {
+    .Call(`_AlphaSimR_solveRRBLUP_EM3`, Y, X, M1, M2, M3, Vu1, Vu2, Vu3, Ve, tol, maxIter, useEM)
+}
+
 callFastRRBLUP <- function(y, geno, lociPerChr, lociLoc, Vu, Ve, maxIter, nThreads) {
     .Call(`_AlphaSimR_callFastRRBLUP`, y, geno, lociPerChr, lociLoc, Vu, Ve, maxIter, nThreads)
 }
@@ -41,6 +155,55 @@ callRRBLUP_SCA2 <- function(y, x, reps, geno, lociPerChr, lociLoc, maxIter, Vu1,
     .Call(`_AlphaSimR_callRRBLUP_SCA2`, y, x, reps, geno, lociPerChr, lociLoc, maxIter, Vu1, Vu2, Vu3, Ve, tol, useEM, useReps, nThreads)
 }
 
+#' @title Solve Univariate Model
+#'
+#' @description
+#' Solves a univariate mixed model of form \eqn{y=X\beta+Zu+e}
+#'
+#' @param y a matrix with n rows and 1 column
+#' @param X a matrix with n rows and x columns
+#' @param Z a matrix with n rows and m columns
+#' @param K a matrix with m rows and m columns
+#'
+#' @export
+solveUVM <- function(y, X, Z, K) {
+    .Call(`_AlphaSimR_solveUVM`, y, X, Z, K)
+}
+
+#' @title Solve Multivariate Model
+#'
+#' @description
+#' Solves a multivariate mixed model of form \eqn{Y=X\beta+Zu+e}
+#'
+#' @param Y a matrix with n rows and q columns
+#' @param X a matrix with n rows and x columns
+#' @param Z a matrix with n rows and m columns
+#' @param K a matrix with m rows and m columns
+#' @param tol tolerance for convergence
+#' @param maxIter maximum number of iteration
+#'
+#' @export
+solveMVM <- function(Y, X, Z, K, tol = 1e-6, maxIter = 1000L) {
+    .Call(`_AlphaSimR_solveMVM`, Y, X, Z, K, tol, maxIter)
+}
+
+#' @title Solve Multikernel Model
+#'
+#' @description
+#' Solves a univariate mixed model with multiple random effects.
+#'
+#' @param y a matrix with n rows and 1 column
+#' @param X a matrix with n rows and x columns
+#' @param Zlist a list of Z matrices
+#' @param Klist a list of K matrices
+#' @param maxIter maximum number of iteration
+#' @param tol tolerance for convergence
+#'
+#' @export
+solveMKM <- function(y, X, Zlist, Klist, maxIter = 40L, tol = 1e-4) {
+    .Call(`_AlphaSimR_solveMKM`, y, X, Zlist, Klist, maxIter, tol)
+}
+
 writeASGenotypes <- function(g, locations, allLocations, snpchips, names, missing, fname) {
     invisible(.Call(`_AlphaSimR_writeASGenotypes`, g, locations, allLocations, snpchips, names, missing, fname))
 }
@@ -166,6 +329,18 @@ calcCoef <- function(X, Y) {
     .Call(`_AlphaSimR_calcCoef`, X, Y)
 }
 
+#' @title Number of available threads
+#'
+#' @description
+#' Gets the number of available threads by calling the OpenMP function
+#' \code{omp_get_max_threads()}
+#'
+#' @return integer
+#'
+#' @examples
+#' getNumThreads()
+#'
+#' @export
 getNumThreads <- function() {
     .Call(`_AlphaSimR_getNumThreads`)
 }
diff --git a/R/crossing.R b/R/crossing.R
index de580264..f37552ca 100644
--- a/R/crossing.R
+++ b/R/crossing.R
@@ -66,6 +66,7 @@ makeCross = function(pop,crossPlan,nProgeny=1,
               simParam$maleCentromere,
               simParam$quadProb,
               simParam$nThreads)
+  dim(tmp$geno) = NULL # Account for matrix bug in RcppArmadillo
   rPop = new("RawPop",
              nInd=nrow(crossPlan),
              nChr=pop@nChr,
@@ -329,6 +330,7 @@ makeCross2 = function(females,males,crossPlan,nProgeny=1,simParam=NULL){
             simParam$maleCentromere,
             simParam$quadProb,
             simParam$nThreads)
+  dim(tmp$geno) = NULL # Account for matrix bug in RcppArmadillo
   rPop = new("RawPop",
              nInd=nrow(crossPlan),
              nChr=females@nChr,
@@ -512,6 +514,7 @@ self = function(pop,nProgeny=1,parents=NULL,keepParents=TRUE,
               simParam$maleCentromere,
               simParam$quadProb,
               simParam$nThreads)
+  dim(tmp$geno) = NULL # Account for matrix bug in RcppArmadillo
   rPop = new("RawPop",
              nInd=nrow(crossPlan),
              nChr=pop@nChr,
@@ -603,6 +606,7 @@ makeDH = function(pop,nDH=1,useFemale=TRUE,keepParents=TRUE,
                     simParam$isTrackRec,
                     simParam$nThreads)
   }
+  dim(tmp$geno) = NULL # Account for matrix bug in RcppArmadillo
   rPop = new("RawPop",
              nInd=as.integer(pop@nInd*nDH),
              nChr=pop@nChr,
diff --git a/R/founderPop.R b/R/founderPop.R
index 9fa5c88f..0485f6bf 100644
--- a/R/founderPop.R
+++ b/R/founderPop.R
@@ -73,7 +73,7 @@ newMapPop = function(genMap, haplotypes, inbred=FALSE,
                         nChr=1L,
                         ploidy=ploidy,
                         nLoci=as.integer(segSites[chr]),
-                        geno=as.matrix(list(geno)),
+                        geno=list(geno),
                         genMap=genMap[chr],
                         centromere=max(genMap[[chr]])/2,
                         inbred=inbred)
@@ -204,6 +204,7 @@ runMacs = function(nInd,nChr=1, segSites=NULL, inbred=FALSE, species="GENERIC",
   }
   macsOut = MaCS(command, segSites, inbred, ploidy, 
                  nThreads, seed)
+  dim(macsOut$geno) = NULL # Account for matrix bug in RcppArmadillo
   nLoci = sapply(macsOut$genMap,length)
   genMap = vector("list",nChr)
   for(i in 1:nChr){
@@ -367,7 +368,7 @@ sampleHaplo = function(mapPop,nInd,inbred=FALSE,ploidy=NULL,replace=TRUE){
                         nChr=1L,
                         ploidy=as.integer(ploidy),
                         nLoci=mapPop@nLoci[chr],
-                        geno=as.matrix(list(geno)),
+                        geno=list(geno),
                         genMap=mapPop@genMap[chr],
                         centromere=mapPop@centromere[chr],
                         inbred=inbred)
@@ -431,7 +432,7 @@ quickHaplo = function(nInd,nChr,segSites,genLen=1,ploidy=2L,inbred=FALSE){
              nChr=nChr,
              ploidy=ploidy,
              nLoci=segSites,
-             geno=as.matrix(geno),
+             geno=geno,
              genMap=genMap,
              centromere=centromere,
              inbred=inbred))
diff --git a/R/mergePops.R b/R/mergePops.R
index 4e4e8832..ce4e7fab 100644
--- a/R/mergePops.R
+++ b/R/mergePops.R
@@ -129,6 +129,7 @@ mergePops = function(popList){
   #geno
   nBin = as.integer(nLoci%/%8L + (nLoci%%8L > 0L))
   geno = mergeMultGeno(popList,nInd=nInd,nBin=nBin,ploidy=ploidy)
+  dim(geno) = NULL # Account for matrix bug in RcppArmadillo
   nInd = sum(nInd)
   return(new("Pop",
              nInd=nInd,
diff --git a/R/polyploids.R b/R/polyploids.R
index 678e8b72..5ac79825 100644
--- a/R/polyploids.R
+++ b/R/polyploids.R
@@ -62,6 +62,7 @@ reduceGenome = function(pop,nProgeny=1,useFemale=TRUE,keepParents=TRUE,
                             simParam$femaleCentromere,
                             simParam$quadProb,
                             simParam$nThreads)
+  dim(tmp$geno) = NULL 
   rPop = new("RawPop",
              nInd=as.integer(pop@nInd*nProgeny),
              nChr=pop@nChr,
@@ -69,22 +70,32 @@ reduceGenome = function(pop,nProgeny=1,useFemale=TRUE,keepParents=TRUE,
              nLoci=pop@nLoci,
              geno=tmp$geno)
   if(simParam$isTrackRec){
-    simParam$addToRec(tmp$recHist)
+    hist = tmp$recHist
+  }else{
+    hist = NULL
   }
   if(keepParents){
     return(newPop(rawPop=rPop,
-                  mother=rep(pop@id,each=nProgeny),
-                  father=rep(pop@id,each=nProgeny),
-                  origM=rep(pop@mother,each=nProgeny),
-                  origF=rep(pop@father,each=nProgeny),
-                  isDH=FALSE,
-                  simParam=simParam))
+                  mother=rep(pop@mother,each=nProgeny),
+                  father=rep(pop@father,each=nProgeny),
+                  simParam=simParam,
+                  iMother=rep(pop@iid,each=nProgeny),
+                  iFather=rep(pop@iid,each=nProgeny),
+                  femaleParentPop=pop,
+                  maleParentPop=pop,
+                  hist=hist
+    ))
   }else{
     return(newPop(rawPop=rPop,
                   mother=rep(pop@id,each=nProgeny),
                   father=rep(pop@id,each=nProgeny),
-                  isDH=FALSE,
-                  simParam=simParam))
+                  simParam=simParam,
+                  iMother=rep(pop@iid,each=nProgeny),
+                  iFather=rep(pop@iid,each=nProgeny),
+                  femaleParentPop=pop,
+                  maleParentPop=pop,
+                  hist=hist
+    ))
   }
 }
 
@@ -167,15 +178,19 @@ doubleGenome = function(pop, keepParents=TRUE,
         }
       }
     }
-    simParam$addToRec(newHist)
+  }else{
+    newHist = NULL
   }
   return(newPop(rawPop=rPop,
-                mother=mother,
-                father=father,
-                origM=origM,
-                origF=origF,
+                mother=origM,
+                father=origF,
+                simParam=simParam,
                 isDH=TRUE,
-                simParam=simParam))
+                iMother=pop@iid,
+                iFather=pop@iid,
+                femaleParentPop=pop,
+                maleParentPop=pop,
+                hist=newHist))
 }
 
 #' @title Combine genomes of individuals
@@ -248,7 +263,7 @@ mergeGenome = function(females,males,crossPlan,simParam=NULL){
              nChr=females@nChr,
              ploidy=females@ploidy+males@ploidy,
              nLoci=females@nLoci,
-             geno=as.matrix(geno))
+             geno=geno)
   
   if(simParam$isTrackRec){
     # Duplicate recombination histories
diff --git a/R/popSummary.R b/R/popSummary.R
index 74a52cba..6e9540d7 100644
--- a/R/popSummary.R
+++ b/R/popSummary.R
@@ -161,6 +161,7 @@ genParam = function(pop,simParam=NULL){
     simParam = get("SP",envir=.GlobalEnv)
   }
   stopifnot(class(pop)=="Pop")
+  nInd=nInd(pop)
   gv=NULL
   bv=NULL
   dd=NULL
@@ -215,9 +216,15 @@ genParam = function(pop,simParam=NULL){
       aa = cbind(aa,rep(0,pop@nInd))
       gv_aa = cbind(gv_aa,rep(0,pop@nInd))
     }
-    covAD_L = c(covAD_L,popVar(cbind(bv[,i],dd[,i]))[1,2])
-    covAAA_L = c(covAAA_L,popVar(cbind(bv[,i],aa[,i]))[1,2])
-    covDAA_L = c(covDAA_L,popVar(cbind(dd[,i],aa[,i]))[1,2])
+    if(nInd==1){
+      covAD_L = c(covAD_L,0)
+      covAAA_L = c(covAAA_L,0)
+      covDAA_L = c(covDAA_L,0)
+    } else {
+      covAD_L = c(covAD_L,popVar(cbind(bv[,i],dd[,i]))[1,2])
+      covAAA_L = c(covAAA_L,popVar(cbind(bv[,i],aa[,i]))[1,2])
+      covDAA_L = c(covDAA_L,popVar(cbind(dd[,i],aa[,i]))[1,2])
+    }
   }
   varA = popVar(bv)
   varD = popVar(dd)
diff --git a/R/pullGeno.R b/R/pullGeno.R
index 4429e981..b81a6cea 100644
--- a/R/pullGeno.R
+++ b/R/pullGeno.R
@@ -565,7 +565,9 @@ pullSegSiteHaplo = function(pop, haplo="all",
 #'
 #' @param pop an object of \code{\link{Pop-class}}
 #' @param chr a vector of chromosomes to retrieve. If NULL,
-#' all chromosome are retrieved.
+#' all chromosomes are retrieved.
+#' @param snpChip an integer indicating which SNP array loci 
+#' are to be retrieved. If NULL, all sites are retrieved.
 #' @param simParam an object of \code{\link{SimParam}}
 #'
 #' @return Returns a matrix of SNP haplotypes.
@@ -585,7 +587,7 @@ pullSegSiteHaplo = function(pop, haplo="all",
 #' pullIbdHaplo(pop, simParam=SP)
 #' 
 #' @export
-pullIbdHaplo = function(pop, chr=NULL, simParam=NULL){
+pullIbdHaplo = function(pop, chr=NULL, snpChip=NULL, simParam=NULL){
   if(is.null(simParam)){
     simParam = get("SP",envir=.GlobalEnv)
   }
@@ -614,5 +616,19 @@ pullIbdHaplo = function(pop, chr=NULL, simParam=NULL){
   
   colnames(output) = unlist(lapply(simParam$genMap[chr], names))
   
+  if(!is.null(snpChip)){
+    nLoci = pop@nLoci[chr]
+    tmp = getSnpMap(snpChip=snpChip,simParam=simParam)
+    tmp = tmp[tmp$chr%in%chr,]
+    if(length(chr)>1){
+      for(i in 2:length(chr)){
+        j = chr[i]
+        tmp[tmp$chr==j,"site"] =
+          tmp[tmp$chr==j,"site"] + sum(nLoci[1:(i-1)])
+      }
+    }
+    output = output[,tmp$site,drop=FALSE]
+  }
+  
   return(output)
 }
diff --git a/R/selection.R b/R/selection.R
index c8bd42d7..e2fafe91 100644
--- a/R/selection.R
+++ b/R/selection.R
@@ -130,6 +130,7 @@ getFam = function(pop,famType){
 selectInd = function(pop,nInd,trait=1,use="pheno",sex="B",
                      selectTop=TRUE,returnPop=TRUE,
                      candidates=NULL,simParam=NULL,...){
+  stopifnot(nInd>=0)
   if(is.null(simParam)){
     simParam = get("SP",envir=.GlobalEnv)
   }
@@ -146,7 +147,8 @@ selectInd = function(pop,nInd,trait=1,use="pheno",sex="B",
     eligible = eligible[eligible%in%candidates]
   }
   if(length(eligible)<nInd){
-    stop("Not enough suitable candidates, check request value and sex")
+    nInd = length(eligible)
+    warning("Suitable candidates smaller than nInd, returning ",nInd," individuals")
   }
   response = getResponse(pop=pop,trait=trait,use=use,
                          simParam=simParam,...)
@@ -216,6 +218,7 @@ selectInd = function(pop,nInd,trait=1,use="pheno",sex="B",
 selectFam = function(pop,nFam,trait=1,use="pheno",sex="B",
                      famType="B",selectTop=TRUE,returnPop=TRUE,
                      candidates=NULL,simParam=NULL,...){
+  stopifnot(nFam>=0)
   if(is.null(simParam)){
     simParam = get("SP",envir=.GlobalEnv)
   }
@@ -234,8 +237,8 @@ selectFam = function(pop,nFam,trait=1,use="pheno",sex="B",
   allFam = getFam(pop=pop,famType=famType)
   availFam = allFam[eligible]
   if(nFam>length(unique(availFam))){
-    stop(paste(nFam,"families requested but only",length(unique(availFam)),
-               "families are available"))
+    nFam = length(unique(availFam))
+    warning("Suitable families smaller than nFam, returning ", nFam, " families")
   }
   response = getResponse(pop=pop,trait=trait,use=use,
                          simParam=simParam,...)
@@ -313,6 +316,7 @@ selectFam = function(pop,nFam,trait=1,use="pheno",sex="B",
 selectWithinFam = function(pop,nInd,trait=1,use="pheno",sex="B",
                            famType="B",selectTop=TRUE,returnPop=TRUE,
                            candidates=NULL,simParam=NULL,...){
+  stopifnot(nInd>=0)
   if(is.null(simParam)){
     simParam = get("SP",envir=.GlobalEnv)
   }
@@ -409,6 +413,7 @@ selectOP = function(pop,nInd,nSeeds,probSelf=0,
                     pollenControl=FALSE,trait=1,
                     use="pheno",selectTop=TRUE,
                     candidates=NULL,simParam=NULL,...){
+  stopifnot(nInd>=0)
   if(is.null(simParam)){
     simParam = get("SP",envir=.GlobalEnv)
   }
diff --git a/alphasimr.Rproj b/alphasimr.Rproj
index e4dfdb8d..5a7090ad 100644
--- a/alphasimr.Rproj
+++ b/alphasimr.Rproj
@@ -15,9 +15,9 @@ LaTeX: pdfLaTeX
 BuildType: Package
 PackageUseDevtools: Yes
 PackageInstallArgs: --no-multiarch --with-keep.source
-PackageCheckArgs: --as-cran
+PackageCheckArgs: --as-cran --no-multiarch
 PackageRoxygenize: rd,collate,namespace
 
 PythonType: conda
-PythonVersion: 3.6.10
+PythonVersion: 3.6.13
 PythonPath: C:\Users\gayno\AppData\Local\r-miniconda\envs\r-reticulate\python.exe
diff --git a/man/MapPop-class.Rd b/man/MapPop-class.Rd
index 2310ba1a..ca0033d5 100644
--- a/man/MapPop-class.Rd
+++ b/man/MapPop-class.Rd
@@ -34,7 +34,7 @@ for creating initial populations and setting traits in the
 \section{Slots}{
 
 \describe{
-\item{\code{genMap}}{"matrix" of chromosome genetic maps}
+\item{\code{genMap}}{list of chromosome genetic maps}
 
 \item{\code{centromere}}{vector of centromere positions}
 
diff --git a/man/RawPop-class.Rd b/man/RawPop-class.Rd
index f518a790..890a7697 100644
--- a/man/RawPop-class.Rd
+++ b/man/RawPop-class.Rd
@@ -46,8 +46,8 @@ The raw population class contains only genotype data.
 
 \item{\code{nLoci}}{number of loci per chromosome}
 
-\item{\code{geno}}{"matrix" containing chromosome genotypes. The "matrix" 
-has dimensions nChr by 1 and each element is a three dimensional
-array of raw values. The array dimensions are nLoci by ploidy by nInd.}
+\item{\code{geno}}{list of nChr length containing chromosome genotypes.
+Each element is a three dimensional array of raw values. 
+The array dimensions are nLoci by ploidy by nInd.}
 }}
 
diff --git a/man/getNumThreads.Rd b/man/getNumThreads.Rd
new file mode 100644
index 00000000..cca31e7d
--- /dev/null
+++ b/man/getNumThreads.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/RcppExports.R
+\name{getNumThreads}
+\alias{getNumThreads}
+\title{Number of available threads}
+\usage{
+getNumThreads()
+}
+\value{
+integer
+}
+\description{
+Gets the number of available threads by calling the OpenMP function
+\code{omp_get_max_threads()}
+}
+\examples{
+getNumThreads()
+
+}
diff --git a/man/pullIbdHaplo.Rd b/man/pullIbdHaplo.Rd
index ab7a2674..4f52020d 100644
--- a/man/pullIbdHaplo.Rd
+++ b/man/pullIbdHaplo.Rd
@@ -4,13 +4,16 @@
 \alias{pullIbdHaplo}
 \title{Pull IBD haplotypes}
 \usage{
-pullIbdHaplo(pop, chr = NULL, simParam = NULL)
+pullIbdHaplo(pop, chr = NULL, snpChip = NULL, simParam = NULL)
 }
 \arguments{
 \item{pop}{an object of \code{\link{Pop-class}}}
 
 \item{chr}{a vector of chromosomes to retrieve. If NULL,
-all chromosome are retrieved.}
+all chromosomes are retrieved.}
+
+\item{snpChip}{an integer indicating which SNP array loci 
+are to be retrieved. If NULL, all sites are retrieved.}
 
 \item{simParam}{an object of \code{\link{SimParam}}}
 }
diff --git a/man/solveMKM.Rd b/man/solveMKM.Rd
new file mode 100644
index 00000000..91d7c930
--- /dev/null
+++ b/man/solveMKM.Rd
@@ -0,0 +1,24 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/RcppExports.R
+\name{solveMKM}
+\alias{solveMKM}
+\title{Solve Multikernel Model}
+\usage{
+solveMKM(y, X, Zlist, Klist, maxIter = 40L, tol = 1e-04)
+}
+\arguments{
+\item{y}{a matrix with n rows and 1 column}
+
+\item{X}{a matrix with n rows and x columns}
+
+\item{Zlist}{a list of Z matrices}
+
+\item{Klist}{a list of K matrices}
+
+\item{maxIter}{maximum number of iteration}
+
+\item{tol}{tolerance for convergence}
+}
+\description{
+Solves a univariate mixed model with multiple random effects.
+}
diff --git a/man/solveMVM.Rd b/man/solveMVM.Rd
new file mode 100644
index 00000000..ab9c8582
--- /dev/null
+++ b/man/solveMVM.Rd
@@ -0,0 +1,24 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/RcppExports.R
+\name{solveMVM}
+\alias{solveMVM}
+\title{Solve Multivariate Model}
+\usage{
+solveMVM(Y, X, Z, K, tol = 1e-06, maxIter = 1000L)
+}
+\arguments{
+\item{Y}{a matrix with n rows and q columns}
+
+\item{X}{a matrix with n rows and x columns}
+
+\item{Z}{a matrix with n rows and m columns}
+
+\item{K}{a matrix with m rows and m columns}
+
+\item{tol}{tolerance for convergence}
+
+\item{maxIter}{maximum number of iteration}
+}
+\description{
+Solves a multivariate mixed model of form \eqn{Y=X\beta+Zu+e}
+}
diff --git a/man/solveRRBLUP.Rd b/man/solveRRBLUP.Rd
new file mode 100644
index 00000000..a921eec8
--- /dev/null
+++ b/man/solveRRBLUP.Rd
@@ -0,0 +1,18 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/RcppExports.R
+\name{solveRRBLUP}
+\alias{solveRRBLUP}
+\title{Solve RR-BLUP}
+\usage{
+solveRRBLUP(y, X, M)
+}
+\arguments{
+\item{y}{a matrix with n rows and 1 column}
+
+\item{X}{a matrix with n rows and x columns}
+
+\item{M}{a matrix with n rows and m columns}
+}
+\description{
+Solves a univariate mixed model of form \eqn{y=X\beta+Mu+e}
+}
diff --git a/man/solveRRBLUPMK.Rd b/man/solveRRBLUPMK.Rd
new file mode 100644
index 00000000..393f1cf5
--- /dev/null
+++ b/man/solveRRBLUPMK.Rd
@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/RcppExports.R
+\name{solveRRBLUPMK}
+\alias{solveRRBLUPMK}
+\title{Solve Multikernel RR-BLUP}
+\usage{
+solveRRBLUPMK(y, X, Mlist, maxIter = 40L)
+}
+\arguments{
+\item{y}{a matrix with n rows and 1 column}
+
+\item{X}{a matrix with n rows and x columns}
+
+\item{Mlist}{a list of M matrices}
+
+\item{maxIter}{maximum number of iteration}
+}
+\description{
+Solves a univariate mixed model with multiple random effects.
+}
diff --git a/man/solveRRBLUPMV.Rd b/man/solveRRBLUPMV.Rd
new file mode 100644
index 00000000..40ca24d3
--- /dev/null
+++ b/man/solveRRBLUPMV.Rd
@@ -0,0 +1,22 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/RcppExports.R
+\name{solveRRBLUPMV}
+\alias{solveRRBLUPMV}
+\title{Solve Multivariate RR-BLUP}
+\usage{
+solveRRBLUPMV(Y, X, M, maxIter = 1000L, tol = 1e-06)
+}
+\arguments{
+\item{Y}{a matrix with n rows and q columns}
+
+\item{X}{a matrix with n rows and x columns}
+
+\item{M}{a matrix with n rows and m columns}
+
+\item{maxIter}{maximum number of iteration}
+
+\item{tol}{tolerance for convergence}
+}
+\description{
+Solves a multivariate mixed model of form \eqn{Y=X\beta+Mu+e}
+}
diff --git a/man/solveRRBLUP_EM.Rd b/man/solveRRBLUP_EM.Rd
new file mode 100644
index 00000000..6126c5a4
--- /dev/null
+++ b/man/solveRRBLUP_EM.Rd
@@ -0,0 +1,30 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/RcppExports.R
+\name{solveRRBLUP_EM}
+\alias{solveRRBLUP_EM}
+\title{Solve RR-BLUP with EM}
+\usage{
+solveRRBLUP_EM(Y, X, M, Vu, Ve, tol, maxIter, useEM)
+}
+\arguments{
+\item{Y}{a matrix with n rows and 1 column}
+
+\item{X}{a matrix with n rows and x columns}
+
+\item{M}{a matrix with n rows and m columns}
+
+\item{Vu}{initial guess for variance of marker effects}
+
+\item{Ve}{initial guess for error variance}
+
+\item{tol}{tolerance for declaring convergence}
+
+\item{maxIter}{maximum iteration for attempting convergence}
+
+\item{useEM}{should EM algorithm be used. If false, no estimation of
+variance components is performed. The initial values are treated as true.}
+}
+\description{
+Solves a univariate mixed model of form \eqn{y=X\beta+Mu+e} using
+the Expectation-Maximization algorithm.
+}
diff --git a/man/solveRRBLUP_EM2.Rd b/man/solveRRBLUP_EM2.Rd
new file mode 100644
index 00000000..7e95c6bf
--- /dev/null
+++ b/man/solveRRBLUP_EM2.Rd
@@ -0,0 +1,34 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/RcppExports.R
+\name{solveRRBLUP_EM2}
+\alias{solveRRBLUP_EM2}
+\title{Solve RR-BLUP with EM and 2 random effects}
+\usage{
+solveRRBLUP_EM2(Y, X, M1, M2, Vu1, Vu2, Ve, tol, maxIter, useEM)
+}
+\arguments{
+\item{Y}{a matrix with n rows and 1 column}
+
+\item{X}{a matrix with n rows and x columns}
+
+\item{M1}{a matrix with n rows and m1 columns}
+
+\item{M2}{a matrix with n rows and m2 columns}
+
+\item{Vu1}{initial guess for variance of the first marker effects}
+
+\item{Vu2}{initial guess for variance of the second marker effects}
+
+\item{Ve}{initial guess for error variance}
+
+\item{tol}{tolerance for declaring convergence}
+
+\item{maxIter}{maximum iteration for attempting convergence}
+
+\item{useEM}{should EM algorithm be used. If false, no estimation of
+variance components is performed. The initial values are treated as true.}
+}
+\description{
+Solves a univariate mixed model of form \eqn{y=X\beta+M_1u_1+M_2u_2+e} using
+the Expectation-Maximization algorithm.
+}
diff --git a/man/solveRRBLUP_EM3.Rd b/man/solveRRBLUP_EM3.Rd
new file mode 100644
index 00000000..e9b1fb3f
--- /dev/null
+++ b/man/solveRRBLUP_EM3.Rd
@@ -0,0 +1,38 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/RcppExports.R
+\name{solveRRBLUP_EM3}
+\alias{solveRRBLUP_EM3}
+\title{Solve RR-BLUP with EM and 3 random effects}
+\usage{
+solveRRBLUP_EM3(Y, X, M1, M2, M3, Vu1, Vu2, Vu3, Ve, tol, maxIter, useEM)
+}
+\arguments{
+\item{Y}{a matrix with n rows and 1 column}
+
+\item{X}{a matrix with n rows and x columns}
+
+\item{M1}{a matrix with n rows and m1 columns}
+
+\item{M2}{a matrix with n rows and m2 columns}
+
+\item{M3}{a matrix with n rows and m3 columns}
+
+\item{Vu1}{initial guess for variance of the first marker effects}
+
+\item{Vu2}{initial guess for variance of the second marker effects}
+
+\item{Vu3}{initial guess for variance of the second marker effects}
+
+\item{Ve}{initial guess for error variance}
+
+\item{tol}{tolerance for declaring convergence}
+
+\item{maxIter}{maximum iteration for attempting convergence}
+
+\item{useEM}{should EM algorithm be used. If false, no estimation of
+variance components is performed. The initial values are treated as true.}
+}
+\description{
+Solves a univariate mixed model of form \eqn{y=X\beta+M_1u_1+M_2u_2+M_3u_3+e} using
+the Expectation-Maximization algorithm.
+}
diff --git a/man/solveUVM.Rd b/man/solveUVM.Rd
new file mode 100644
index 00000000..358a0b4c
--- /dev/null
+++ b/man/solveUVM.Rd
@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/RcppExports.R
+\name{solveUVM}
+\alias{solveUVM}
+\title{Solve Univariate Model}
+\usage{
+solveUVM(y, X, Z, K)
+}
+\arguments{
+\item{y}{a matrix with n rows and 1 column}
+
+\item{X}{a matrix with n rows and x columns}
+
+\item{Z}{a matrix with n rows and m columns}
+
+\item{K}{a matrix with m rows and m columns}
+}
+\description{
+Solves a univariate mixed model of form \eqn{y=X\beta+Zu+e}
+}
diff --git a/src/MME.cpp b/src/MME.cpp
index 389277da..e11004e2 100644
--- a/src/MME.cpp
+++ b/src/MME.cpp
@@ -17,6 +17,7 @@ void arma_fortran(arma_dsyevr)(char* JOBZ, char* RANGE, char* UPLO, long long in
                   long long int* LIWORK, long long int* INFO);
 #endif
 
+const double pi = 3.14159265358979323846;
 
 // Replacement for Armadillo's eig_sym
 // Fixes an error with decompisition of large matrices
@@ -110,6 +111,8 @@ arma::mat makeZ(arma::uvec& z, arma::uword nGeno){
   return Z;
 }
 
+
+
 // Generates weighted matrix
 // Allows for heterogenous variance due to unequal replication
 void sweepReps(arma::mat& X, arma::vec reps){
@@ -119,6 +122,17 @@ void sweepReps(arma::mat& X, arma::vec reps){
   }
 }
 
+//' @title Solve RR-BLUP
+//'
+//' @description
+//' Solves a univariate mixed model of form \eqn{y=X\beta+Mu+e}
+//'
+//' @param y a matrix with n rows and 1 column
+//' @param X a matrix with n rows and x columns
+//' @param M a matrix with n rows and m columns
+//'
+//' @export
+// [[Rcpp::export]]
 Rcpp::List solveRRBLUP(const arma::mat& y, const arma::mat& X,
                        const arma::mat& M){
   arma::uword n = y.n_rows;
@@ -165,6 +179,19 @@ Rcpp::List solveRRBLUP(const arma::mat& y, const arma::mat& X,
                             Rcpp::Named("u")=u);
 }
 
+//' @title Solve Multivariate RR-BLUP
+//'
+//' @description
+//' Solves a multivariate mixed model of form \eqn{Y=X\beta+Mu+e}
+//'
+//' @param Y a matrix with n rows and q columns
+//' @param X a matrix with n rows and x columns
+//' @param M a matrix with n rows and m columns
+//' @param tol tolerance for convergence
+//' @param maxIter maximum number of iteration
+//'
+//' @export
+// [[Rcpp::export]]
 Rcpp::List solveRRBLUPMV(const arma::mat& Y, const arma::mat& X,
                          const arma::mat& M, int maxIter=1000, 
                          double tol=1e-6){
@@ -227,6 +254,18 @@ Rcpp::List solveRRBLUPMV(const arma::mat& Y, const arma::mat& X,
                             Rcpp::Named("iter")=iter);
 }
 
+//' @title Solve Multikernel RR-BLUP
+//'
+//' @description
+//' Solves a univariate mixed model with multiple random effects.
+//'
+//' @param y a matrix with n rows and 1 column
+//' @param X a matrix with n rows and x columns
+//' @param Mlist a list of M matrices
+//' @param maxIter maximum number of iteration
+//'
+//' @export
+// [[Rcpp::export]]
 Rcpp::List solveRRBLUPMK(arma::mat& y, arma::mat& X,
                          arma::field<arma::mat>& Mlist,
                          int maxIter=40){
@@ -326,7 +365,24 @@ Rcpp::List solveRRBLUPMK(arma::mat& y, arma::mat& X,
                             Rcpp::Named("iter")=iter);
 }
 
-//Uses EM algorithm to solve a mixed model with 1 random effect
+//' @title Solve RR-BLUP with EM
+//'
+//' @description
+//' Solves a univariate mixed model of form \eqn{y=X\beta+Mu+e} using
+//' the Expectation-Maximization algorithm.
+//'
+//' @param Y a matrix with n rows and 1 column
+//' @param X a matrix with n rows and x columns
+//' @param M a matrix with n rows and m columns
+//' @param Vu initial guess for variance of marker effects
+//' @param Ve initial guess for error variance
+//' @param tol tolerance for declaring convergence
+//' @param maxIter maximum iteration for attempting convergence
+//' @param useEM should EM algorithm be used. If false, no estimation of
+//' variance components is performed. The initial values are treated as true.
+//'
+//' @export
+// [[Rcpp::export]]
 Rcpp::List solveRRBLUP_EM(arma::mat& Y, arma::mat& X,
                           arma::mat& M, double Vu, double Ve, 
                           double tol, int maxIter,
@@ -389,7 +445,26 @@ Rcpp::List solveRRBLUP_EM(arma::mat& Y, arma::mat& X,
                             Rcpp::Named("iter")=iter);
 }
 
-//Uses EM algorithm to solve a mixed model with 2 random effects
+//' @title Solve RR-BLUP with EM and 2 random effects
+//'
+//' @description
+//' Solves a univariate mixed model of form \eqn{y=X\beta+M_1u_1+M_2u_2+e} using
+//' the Expectation-Maximization algorithm.
+//'
+//' @param Y a matrix with n rows and 1 column
+//' @param X a matrix with n rows and x columns
+//' @param M1 a matrix with n rows and m1 columns
+//' @param M2 a matrix with n rows and m2 columns
+//' @param Vu1 initial guess for variance of the first marker effects
+//' @param Vu2 initial guess for variance of the second marker effects
+//' @param Ve initial guess for error variance
+//' @param tol tolerance for declaring convergence
+//' @param maxIter maximum iteration for attempting convergence
+//' @param useEM should EM algorithm be used. If false, no estimation of
+//' variance components is performed. The initial values are treated as true.
+//'
+//' @export
+// [[Rcpp::export]]
 Rcpp::List solveRRBLUP_EM2(const arma::mat& Y, const arma::mat& X,
                            const arma::mat& M1, const arma::mat& M2, 
                            double Vu1, double Vu2, double Ve, 
@@ -488,7 +563,28 @@ Rcpp::List solveRRBLUP_EM2(const arma::mat& Y, const arma::mat& X,
                             Rcpp::Named("iter")=iter);
 }
 
-//Uses EM algorithm to solve a mixed model with 3 random effects
+//' @title Solve RR-BLUP with EM and 3 random effects
+//'
+//' @description
+//' Solves a univariate mixed model of form \eqn{y=X\beta+M_1u_1+M_2u_2+M_3u_3+e} using
+//' the Expectation-Maximization algorithm.
+//'
+//' @param Y a matrix with n rows and 1 column
+//' @param X a matrix with n rows and x columns
+//' @param M1 a matrix with n rows and m1 columns
+//' @param M2 a matrix with n rows and m2 columns
+//' @param M3 a matrix with n rows and m3 columns
+//' @param Vu1 initial guess for variance of the first marker effects
+//' @param Vu2 initial guess for variance of the second marker effects
+//' @param Vu3 initial guess for variance of the second marker effects
+//' @param Ve initial guess for error variance
+//' @param tol tolerance for declaring convergence
+//' @param maxIter maximum iteration for attempting convergence
+//' @param useEM should EM algorithm be used. If false, no estimation of
+//' variance components is performed. The initial values are treated as true.
+//'
+//' @export
+// [[Rcpp::export]]
 Rcpp::List solveRRBLUP_EM3(const arma::mat& Y, const arma::mat& X,
                            const arma::mat& M1, const arma::mat& M2, 
                            const arma::mat& M3, double Vu1, double Vu2, 
@@ -1256,3 +1352,266 @@ Rcpp::List callRRBLUP_SCA2(arma::mat y, arma::uvec x, arma::vec reps,
     Rcpp::Named("Ve")=ans["Ve"]
   );
 }
+
+//' @title Solve Univariate Model
+//'
+//' @description
+//' Solves a univariate mixed model of form \eqn{y=X\beta+Zu+e}
+//'
+//' @param y a matrix with n rows and 1 column
+//' @param X a matrix with n rows and x columns
+//' @param Z a matrix with n rows and m columns
+//' @param K a matrix with m rows and m columns
+//'
+//' @export
+// [[Rcpp::export]]
+Rcpp::List solveUVM(const arma::mat& y, const arma::mat& X,
+                    const arma::mat& Z, const arma::mat& K){
+  arma::uword n = y.n_rows;
+  arma::uword q = X.n_cols;
+  double df = double(n)-double(q);
+  double offset = log(double(n));
+  
+  // Construct system of equations for eigendecomposition
+  arma::mat S = -(X*inv_sympd(X.t()*X)*X.t());
+  S.diag() += 1;
+  arma::mat ZK = Z*K;
+  arma::mat H = ZK*Z.t(); // Used later
+  H.diag() += offset;
+  S = S*H*S;
+  
+  // Compute eigendecomposition
+  arma::vec eigval(n);
+  arma::mat eigvec(n,n);
+  eigen2(eigval, eigvec, S);
+  
+  // Drop eigenvalues
+  eigval = eigval(arma::span(q,eigvec.n_cols-1)) - offset;
+  eigvec = eigvec(arma::span(0,eigvec.n_rows-1),
+                  arma::span(q,eigvec.n_cols-1));
+  
+  // Estimate variances and solve equations
+  arma::vec eta = eigvec.t()*y;
+  Rcpp::List optRes = optimize(*objREML,
+                               Rcpp::List::create(
+                                 Rcpp::Named("df")=df,
+                                 Rcpp::Named("eta")=eta,
+                                 Rcpp::Named("lambda")=eigval),
+                                 1.0e-10, 1.0e10);
+  double delta = optRes["parameter"];
+  H.diag() += (delta-offset);
+  H = inv_sympd(H);
+  arma::mat XH = X.t()*H;
+  arma::mat beta = solve(XH*X,XH*y);
+  arma::mat u = ZK.t()*(H*(y-X*beta));
+  double Vu = sum(eta%eta/(eigval+delta))/df;
+  double Ve = delta*Vu;
+  double ll = -0.5*(double(optRes["objective"])+df+df*log(2*pi/df));
+  return Rcpp::List::create(Rcpp::Named("Vu")=Vu,
+                            Rcpp::Named("Ve")=Ve,
+                            Rcpp::Named("beta")=beta,
+                            Rcpp::Named("u")=u,
+                            Rcpp::Named("LL")=ll);
+}
+
+//' @title Solve Multivariate Model
+//'
+//' @description
+//' Solves a multivariate mixed model of form \eqn{Y=X\beta+Zu+e}
+//'
+//' @param Y a matrix with n rows and q columns
+//' @param X a matrix with n rows and x columns
+//' @param Z a matrix with n rows and m columns
+//' @param K a matrix with m rows and m columns
+//' @param tol tolerance for convergence
+//' @param maxIter maximum number of iteration
+//'
+//' @export
+// [[Rcpp::export]]
+Rcpp::List solveMVM(const arma::mat& Y, const arma::mat& X,
+                    const arma::mat& Z, const arma::mat& K,
+                    double tol=1e-6, int maxIter=1000){
+  arma::uword n = Y.n_rows;
+  arma::uword m = Y.n_cols;
+  arma::mat ZK = Z*K;
+  arma::mat ZKZ = ZK*Z.t();
+  arma::vec eigval(n);
+  arma::mat eigvec(n,n);
+  eigen2(eigval, eigvec, ZKZ);
+  arma::mat Yt = Y.t()*eigvec;
+  arma::mat Xt = X.t()*eigvec;
+  arma::mat Vu = cov(Y)/2;
+  arma::mat Ve = Vu;
+  arma::mat W = Xt.t()*inv_sympd(Xt*Xt.t());
+  arma::mat B = Yt*W; //BLUEs
+  arma::mat Gt(m,n), sigma(m,m), BNew,
+  VeNew(m,m), VuNew(m,m);
+  double denom, numer;
+  bool converging=true;
+  int iter=0;
+  while(converging){
+    ++iter;
+    VeNew.fill(0.0);
+    VuNew.fill(0.0);
+    for(arma::uword i=0; i<n; ++i){
+      Gt.col(i) = eigval(i)*Vu*inv_sympd(eigval(i)*Vu+
+        Ve+tol*arma::eye(m,m))*(Yt.col(i)-B*Xt.col(i));
+    }
+    BNew = (Yt - Gt)*W;
+    for(arma::uword i=0; i<n; ++i){
+      sigma = eigval(i)*Vu-(eigval(i)*Vu)*inv_sympd(eigval(i)*Vu+
+        Ve+tol*arma::eye(m,m))*(eigval(i)*Vu);
+      VuNew += 1.0/(double(n)*eigval(i))*(Gt.col(i)*Gt.col(i).t()+sigma);
+      VeNew += 1.0/double(n)*((Yt.col(i)-BNew*Xt.col(i)-Gt.col(i))*
+        (Yt.col(i)-BNew*Xt.col(i)-Gt.col(i)).t()+sigma);
+    }
+    denom = fabs(sum(Ve.diag()));
+    if(denom>0.0){
+      numer = fabs(sum(VeNew.diag()-Ve.diag()));
+      if((numer/denom)<tol) converging=false;
+    }
+    Ve = VeNew;
+    Vu = VuNew;
+    B = BNew;
+    if(iter>=maxIter){
+      Rf_warning("Reached maxIter without converging");
+      break;
+    }
+  }
+  arma::mat HI = inv_sympd(kron(ZKZ, Vu)+
+    kron(arma::eye(n,n), Ve)+
+    tol*arma::eye(n*m,n*m));
+  arma::mat E = Y.t() - B*X.t();
+  arma::mat U = kron(K, Vu)*kron(Z.t(),
+                     arma::eye(m,m))*(HI*vectorise(E)); //BLUPs
+  U.reshape(m,U.n_elem/m);
+  //Log Likelihood calculation
+  arma::mat ll = -0.5*arma::vectorise(E).t()*HI*vectorise(E);
+  ll -= double(n*m)/2.0*log(2*pi);
+  double value;
+  double sign;
+  log_det(value, sign, kron(ZKZ, Vu)+kron(arma::eye(n,n), Ve));
+  ll -= 0.5*value*sign;
+  return Rcpp::List::create(Rcpp::Named("Vu")=Vu,
+                            Rcpp::Named("Ve")=Ve,
+                            Rcpp::Named("beta")=B.t(),
+                            Rcpp::Named("u")=U.t(),
+                            Rcpp::Named("LL")=arma::as_scalar(ll),
+                            Rcpp::Named("iter")=iter);
+}
+
+//' @title Solve Multikernel Model
+//'
+//' @description
+//' Solves a univariate mixed model with multiple random effects.
+//'
+//' @param y a matrix with n rows and 1 column
+//' @param X a matrix with n rows and x columns
+//' @param Zlist a list of Z matrices
+//' @param Klist a list of K matrices
+//' @param maxIter maximum number of iteration
+//' @param tol tolerance for convergence
+//'
+//' @export
+// [[Rcpp::export]]
+Rcpp::List solveMKM(arma::mat& y, arma::mat& X,
+                    arma::field<arma::mat>& Zlist,
+                    arma::field<arma::mat>& Klist,
+                    int maxIter=40, double tol=1e-4){
+  arma::uword k = Klist.n_elem;
+  arma::uword n = y.n_rows;
+  arma::uword q = X.n_cols;
+  double df = double(n)-double(q);
+  arma::field<arma::mat> V(k);
+  for(arma::uword i=0; i<k; ++i){
+    V(i) = Zlist(i)*Klist(i)*Zlist(i).t();
+  }
+  arma::mat A(k+1,k+1), W0(n,n), W(n,n), WX(n,q), WQX(n,n);
+  arma::vec qvec(k+1), sigma(k+1);
+  double rss, ldet, llik, llik0=0, deltaLlik, taper,
+  value, sign;
+  bool invPass;
+  arma::field<arma::mat> T(k);
+  sigma.fill(var(y.col(0)));
+  int iter = 0;
+  while(true){
+    ++iter;
+    W0 = V(0)*sigma(0);
+    W0.diag() += sigma(k);
+    for(arma::uword i=1; i<k; ++i){
+      W0 += V(i)*sigma(i);
+    }
+    invPass = inv_sympd(W,W0);
+    if(!invPass){
+      W = pinv(W0);
+    }
+    WX = W*X;
+    WQX = W - WX*solve(X.t()*WX, WX.t());
+    rss = as_scalar(y.t()*WQX*y);
+    sigma = sigma*(rss/df);
+    WQX = WQX*(df/rss);
+    log_det(value, sign, WQX);
+    ldet = value*sign;
+    llik = ldet/2 - df/2;
+    if(iter == 1) llik0 = llik;
+    deltaLlik = llik - llik0;
+    llik0 = llik;
+    for(arma::uword i=0; i<k; ++i){
+      T(i) = WQX*V(i);
+    }
+    for(arma::uword i=0; i<k; ++i){
+      qvec(i) = as_scalar(y.t()*T(i)*WQX*y - sum(T(i).diag()));
+      for(arma::uword j=0; j<k; ++j){
+        A(i,j) = accu(T(i)%T(j).t());
+      }
+      A(i,k) = accu(T(i)%WQX.t());
+    }
+    for(arma::uword j=0; j<k; ++j){
+      A(k,j) = accu(WQX%T(j).t());
+    }
+    A(k,k) = accu(WQX%WQX.t());
+    qvec(k) = as_scalar(y.t()*WQX*WQX*y - sum(WQX.diag()));
+    A = pinv(A);
+    qvec = A*qvec;
+    if(iter == 1){
+      taper = 0.5;
+    }else if(iter == 2){
+      taper = 0.7;
+    }else{
+      taper = 0.9;
+    }
+    sigma += taper*qvec;
+    while(sigma.min() < -(1e-6)){
+      sigma(sigma.index_min()) = -(1e-6);
+    }
+    if((iter > 1) & (fabs(deltaLlik) < tol*10)){
+      break;
+    }
+    if(max(abs(qvec)) < tol){
+      break;
+    }
+    if(iter >= maxIter){
+      Rf_warning("Reached maxIter without converging");
+      break;
+    }
+  }
+  while(sigma.min() < 0.0){
+    sigma(sigma.index_min()) = 0.0;
+  }
+  arma::mat beta(q,1), ee(n,1);
+  arma::field<arma::mat> u(k);
+  beta = solve(X.t()*W*X,X.t()*W*y);
+  ee = y - X*beta;
+  for(arma::uword i=0; i<k; ++i){
+    u(i) = (Klist(i)*sigma(i))*Zlist(i).t()*W*ee;
+  }
+  arma::vec Vu(k), Ve(1);
+  Vu = sigma(arma::span(0,k-1));
+  Ve = sigma(k);
+  return Rcpp::List::create(Rcpp::Named("Vu")=Vu,
+                            Rcpp::Named("Ve")=Ve,
+                            Rcpp::Named("beta")=beta,
+                            Rcpp::Named("u")=u,
+                            Rcpp::Named("LL")=llik,
+                            Rcpp::Named("iter")=iter);
+}
diff --git a/src/OCS.cpp b/src/OCS.cpp
new file mode 100644
index 00000000..fba28a1f
--- /dev/null
+++ b/src/OCS.cpp
@@ -0,0 +1,362 @@
+#include "alphasimr.h"
+
+/*
+ * OCS
+ *  1. Crossing
+ *  2. Selection
+ *  
+ * Usefulness
+ *  1. DH
+ *  2. Fx
+ *  3. Gametic variance
+ *  EBV, GV, Pheno, index
+ *  Steps
+ *    1. Pull relevant markers
+ *    2. Create CO progeny mask population
+ *    3. Set parental genotypes to estimated effects
+ *    3. Drop parental effects through mask
+ *    4. Calculate usefulness
+ *  
+ *  Pedigree relationship functions from pedigreeTools
+ * 
+ */
+
+arma::vec calcUsefulness(const arma::imat& parents, //Haplotypes (0-1)
+                         const arma::imat& progeny, //Haplotypes (0-3)
+                         const arma::vec& a,
+                         int nBest=50,
+                         int nThreads=4){
+        arma::uword nParents = parents.n_rows/2;
+        arma::uword nProgeny = progeny.n_rows/2;
+        arma::uword nLoci = parents.n_cols;
+        arma::vec usefulness(nParents*(nParents-1)/2);
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static) num_threads(nThreads)
+#endif
+        for(arma::uword n=0; n<usefulness.n_elem; ++n){
+                arma::uword i = mapRow(n,nParents);
+                arma::uword j = mapCol(i,n,nParents);
+                arma::vec blup(nProgeny, arma::fill::zeros);
+                for(arma::uword k=0; k<nLoci; ++k){
+                        arma::vec x(4);
+                        x(0) = a(k)*parents(2*i,k);
+                        x(1) = a(k)*parents(2*i+1,k);
+                        x(2) = a(k)*parents(2*j,k);
+                        x(3) = a(k)*parents(2*j+1,k);
+                        for(arma::uword l=0; l<(2*nProgeny); ++l){
+                                blup(l/2) += x(progeny(l,k));
+                        }
+                }
+                blup = sort(blup,"descend");
+                usefulness(n) = mean(blup(arma::span(0,nBest-1)));
+        }
+        return usefulness;
+}
+
+// Finds parental contributions for a crossing plan
+// crosses = the crossing plan
+// nInd = number of potential parents
+arma::vec calcContr(arma::uvec crosses, 
+                    const arma::uword& nInd,
+                    arma::uvec& rowPos, 
+                    arma::uvec& colPos){
+        double val = 1/(2*double(crosses.n_elem));
+        arma::vec x(nInd, arma::fill::zeros);
+        for(arma::uvec::iterator i=crosses.begin(); i!=crosses.end(); ++i){
+                x(rowPos(*i)) += val;
+                x(colPos(*i)) += val;
+        }
+        return x;
+};
+
+// Calculate the angle and length of a solution vector
+// angle = angle of solution vector (output variable)
+// length = length of solution vector (output variable)
+// u = usefulness of solution
+// sim = similarity of solution
+// uMax = usefulness of maximum usefulness solution
+// simMax = similarity of maximum usefulness solution
+// uMin = usefulness of minimum similarity solution
+// simMin = similarity of minimum similarity solution
+void calcVec(double& angle, 
+             double& length,
+             double u, 
+             double sim,
+             const double& uMax, 
+             const double& simMax,
+             const double& uMin, 
+             const double& simMin){
+        u = (u-uMin)/(uMax-uMin);
+        sim = (simMax-sim)/(simMax-simMin);
+        length = sqrt(u*u+sim*sim);
+        angle = acos(u/length);
+        if(u<0){
+                length = -length;
+        }
+}
+
+// Performs mating between two crossing plan
+// uses random selection of unique crosses
+// a = a crossing plan
+// b = a crossing plan (can have different length than a)
+// returns crossing plan with same length as a
+// Note: This implementation is slow and could be improved
+arma::uvec mate(arma::uvec a, arma::uvec b){
+        arma::uvec c = shuffle(unique(join_cols(a,b)));
+        c.resize(a.n_elem);
+        return c;
+}
+
+// Peforms mutation of a plan
+// mutation is modeled as mating with a random crossing plan
+// crosses = crossing plan to mutate (input/output variable)
+// nMutate = length of random crossing plan
+// potCross = number of potential crosses
+void mutate(arma::uvec crosses, const arma::uword& nMutate, 
+            const arma::uword& potCross){
+        arma::uvec mutations = sampleInt(nMutate,potCross);
+        crosses = mate(crosses,mutations);
+}
+
+// Quick and dirty implementation of a genetic algorithm (GA) 
+// Uses usefulness criteria for crossing plan optimization
+// Crosses are selected from a half diallele without selfs
+Rcpp::List selectCrosses(arma::uword nCross, //Number of crosses to make
+                         double targetAngle, //Number of degrees to target from maximum usefulness (in radians)
+                         arma::vec& u, //Vector of usefulness criterion for crosses
+                         arma::mat& G, //Relationship matrix for individuals 
+                         double probMut=0.01, //Mutation probability for progeny in GA
+                         arma::uword nMutate=2, //Number of potential mutations in mutated progeny
+                         arma::uword nSel=500, //Number of parents in GA
+                         arma::uword nPop=10000, //Number of progeny in GA
+                         arma::uword maxGen=1000, //Maximum number of generations
+                         arma::uword maxRun=100, //Stopping criteria for maximum number of runs without change
+                         double anglePenalty=0.5, //Penalty to vector length for off angle, higher value emphasises angle more
+                         int nThreads=10){ //Number of cores for OpenMP
+        arma::umat crossPlan(2,nPop); //Crossing plan for GA generations
+        arma::umat outCrossPlan(nCross,2); //Final crossing plan
+        arma::uword nInd = G.n_cols; //Number of parents under consideration
+        arma::uword potCross = nInd*(nInd-1)/2; //Number of potential crosses
+        arma::umat Progeny(nCross,nPop), Parents(nCross,nSel); //Solutions
+        arma::uvec Best(nCross); //Best solution
+        arma::vec uProgeny(nPop), uParents(nSel); //Solution usefulness
+        arma::vec simProgeny(nPop), simParents(nSel); //Solution similarity
+        arma::vec angleProgeny(nPop), angleParents(nSel); //Solution angle
+        arma::vec lenProgeny(nPop), lenParents(nSel); //Solution length
+        arma::vec valProgeny(nPop), valParents(nSel); //Solution value
+        arma::uvec rankProgeny(nPop); //Ranked progreny and select progeny
+        double uBest, simBest, valBest, angleBest, lenBest; //Best solution usefulness, diversity, value, and angle
+        double uMax, uMin, simMax, simMin; //Parameters for min simularity and max usefulness solutions
+        arma::uword currentRun;
+        arma::uvec rowPos(potCross), colPos(potCross); //Map cross number to row and column
+        arma::uword tmpI=0;
+        for(arma::uword i=0; i<(nInd-1); i++){
+                for(arma::uword j=(i+1); j<nInd; j++){
+                        rowPos(tmpI) = i;
+                        colPos(tmpI) = j;
+                        ++tmpI;
+                }
+        }
+        
+        // Calculate maximum usefulness (doesn't require GA)
+        arma::uvec uBestIndex = sort_index(u,"descend");
+        uBestIndex.resize(nCross);
+        arma::vec x = calcContr(uBestIndex.unsafe_col(0), nInd, 
+                                rowPos, colPos);
+        uMax = mean(u(uBestIndex));
+        simMax = as_scalar(x.t()*G*x);
+        
+        if(targetAngle<1e-6){
+                // No need to optimize crossing plan, just return best
+                for(arma::uword i=0; i<nCross; ++i){
+                        outCrossPlan(i,0) = rowPos(uBestIndex(i));
+                        outCrossPlan(i,1) = colPos(uBestIndex(i));
+                }
+                return Rcpp::List::create(Rcpp::Named("crossPlan")=outCrossPlan+1, //C++ to R
+                                          Rcpp::Named("uMax")=uMax,
+                                          Rcpp::Named("simMax")=simMax);
+        }
+        
+        // Calculate minimum simularity solution
+        Rcpp::Rcout<<"Optimize for Simularity"<<std::endl<<std::endl;
+        
+        // Initialize progeny
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static) num_threads(nThreads)
+#endif
+        for(arma::uword i=0; i<nPop; i++){
+                Progeny.col(i) = sampleInt(nCross, potCross);
+                arma::vec x = calcContr(Progeny.unsafe_col(i), nInd, 
+                                        rowPos, colPos);
+                simProgeny(i) = as_scalar(x.t()*G*x);
+        }
+        
+        // Select parents and best solution
+        rankProgeny = sort_index(simProgeny,"ascend");
+        for(arma::uword i=0; i<nSel; i++){
+                Parents.col(i) = Progeny.col(rankProgeny(i));
+                simParents(i) = simProgeny(rankProgeny(i));
+        }
+        simBest = simParents(0);
+        Best = Parents.col(0);
+        
+        // Run GA for simularity
+        Rcpp::Rcout<<"Gen  Simularity"<<std::endl;
+        currentRun = 0;
+        for(arma::uword gen=0; gen<maxGen; gen++){
+                // Mate parents
+                crossPlan = sampHalfDialComb(nSel, nPop);
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static) num_threads(nThreads)
+#endif
+                for(arma::uword i=0; i<nPop; i++){
+                        Progeny.col(i) = 
+                                mate(Parents.unsafe_col(crossPlan(0,i)),
+                                     Parents.unsafe_col(crossPlan(1,i)));
+                        arma::vec p(1,arma::fill::randu);
+                        if(as_scalar(p)<probMut){
+                                mutate(Progeny.unsafe_col(i),nMutate,potCross);
+                        }
+                        arma::vec x = calcContr(Progeny.unsafe_col(i), nInd, 
+                                                rowPos, colPos);
+                        simProgeny(i) = as_scalar(x.t()*G*x);
+                }
+                // Select parents and best solution
+                rankProgeny = sort_index(simProgeny,"ascend");
+                for(arma::uword i=0; i<nSel; i++){
+                        Parents.col(i) = Progeny.col(rankProgeny(i));
+                        simParents(i) = simProgeny(rankProgeny(i));
+                }
+                if(simParents(0)<simBest){
+                        simBest = simParents(0);
+                        Best = Parents.col(0);
+                        currentRun = 0;
+                }else{
+                        ++currentRun;
+                }
+                
+                //Report status
+                if(gen%10 == 0){
+                        Rcpp::Rcout<<gen<<"  "<<simBest<<std::endl;
+                }
+                
+                //Terminate if solution isn't changing
+                if(currentRun>=maxRun){
+                        break;
+                }
+        }
+        simMin = simBest;
+        uMin = mean(u(Best));
+        
+        // Perform optimization for best solution
+        Rcpp::Rcout<<std::endl<<std::endl<<"Optimize for Crossing Plan"<<std::endl<<std::endl;
+        
+        // Initialize progeny
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static) num_threads(nThreads)
+#endif 
+        for(arma::uword i=0; i<nPop; i++){
+                Progeny.col(i) = sampleInt(nCross, potCross);
+                arma::vec x = calcContr(Progeny.unsafe_col(i), nInd, 
+                                        rowPos, colPos);
+                simProgeny(i) = as_scalar(x.t()*G*x);
+                uProgeny(i) = mean(u(Progeny.col(i)));
+                calcVec(angleProgeny(i), lenProgeny(i), uProgeny(i),
+                        simProgeny(i), uMax, simMax, uMin, simMin);
+                valProgeny(i) = lenProgeny(i)-anglePenalty*fabs(angleProgeny(i)-targetAngle);
+        }
+        
+        // Select parents and best solution
+        rankProgeny = sort_index(valProgeny,"descend");
+        for(arma::uword i=0; i<nSel; i++){
+                Parents.col(i) = Progeny.col(rankProgeny(i));
+                uParents(i) = uProgeny(rankProgeny(i));
+                simParents(i) = simProgeny(rankProgeny(i));
+                angleParents(i) = angleProgeny(rankProgeny(i));
+                lenParents(i) = lenProgeny(rankProgeny(i));
+                valParents(i) = valProgeny(rankProgeny(i));
+        }
+        valBest = valParents(0);
+        angleBest = angleParents(0);
+        lenBest = lenParents(0);
+        uBest = uParents(0);
+        simBest = simParents(0);
+        Best = Parents.col(0);
+        
+        // Run GA
+        Rcpp::Rcout<<"Gen  Usefulness  Simularity  Angle  Length  Value"<<std::endl;
+        currentRun = 0;
+        for(arma::uword gen=0; gen<maxGen; gen++){
+                // Mate parents
+                crossPlan = sampHalfDialComb(nSel, nPop);
+#ifdef _OPENMP
+#pragma omp parallel for schedule(static) num_threads(nThreads)
+#endif
+                for(arma::uword i=0; i<nPop; i++){
+                        Progeny.col(i) = 
+                                mate(Parents.unsafe_col(crossPlan(0,i)),
+                                     Parents.unsafe_col(crossPlan(1,i)));
+                        arma::vec p(1,arma::fill::randu);
+                        if(as_scalar(p)<probMut){
+                                mutate(Progeny.unsafe_col(i),nMutate,potCross);
+                        }
+                        arma::vec x = calcContr(Progeny.unsafe_col(i), nInd, 
+                                                rowPos, colPos);
+                        simProgeny(i) = as_scalar(x.t()*G*x);
+                        uProgeny(i) = mean(u(Progeny.col(i)));
+                        calcVec(angleProgeny(i), lenProgeny(i), uProgeny(i),
+                                simProgeny(i), uMax, simMax, uMin, simMin);
+                        valProgeny(i) = lenProgeny(i)-anglePenalty*fabs(angleProgeny(i)-targetAngle);
+                }
+                
+                // Select parents and best solution
+                rankProgeny = sort_index(valProgeny,"descend");
+                for(arma::uword i=0; i<nSel; i++){
+                        Parents.col(i) = Progeny.col(rankProgeny(i));
+                        uParents(i) = uProgeny(rankProgeny(i));
+                        simParents(i) = simProgeny(rankProgeny(i));
+                        angleParents(i) = angleProgeny(rankProgeny(i));
+                        lenParents(i) = lenProgeny(rankProgeny(i));
+                        valParents(i) = valProgeny(rankProgeny(i));
+                }
+                if(valParents(0)>valBest){
+                        valBest = valParents(0);
+                        angleBest = angleParents(0);
+                        lenBest = lenParents(0);
+                        uBest = uParents(0);
+                        simBest = simParents(0);
+                        Best = Parents.col(0);
+                        currentRun = 0;
+                }else{
+                        ++currentRun;
+                }
+                
+                //Report status
+                if(gen%10 == 0){
+                        Rcpp::Rcout<<gen<<"  "<<uBest<<"  "<<simBest<<"  "<<angleBest<<"  "<<lenBest<<"  "<<valBest<<std::endl;
+                }
+                
+                //Terminate if solution isn't changing
+                if(currentRun>=maxRun){
+                        break;
+                }
+        }
+        
+        //Convert solution to an ordered crossing plan
+        Best = sort(Best);
+        for(arma::uword i=0; i<nCross; ++i){
+                outCrossPlan(i,0) = rowPos(Best(i));
+                outCrossPlan(i,1) = colPos(Best(i));
+        }
+        
+        return Rcpp::List::create(Rcpp::Named("crossPlan")=outCrossPlan+1, //C++ to R
+                                  Rcpp::Named("uMax")=uMax,
+                                  Rcpp::Named("uMin")=uMin,
+                                  Rcpp::Named("simMax")=simMax,
+                                  Rcpp::Named("simMin")=simMin,
+                                  Rcpp::Named("uBest")=uBest,
+                                  Rcpp::Named("simBest")=simBest,
+                                  Rcpp::Named("angleBest")=angleBest,
+                                  Rcpp::Named("lenBest")=lenBest);
+}
+
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index fb3cf993..e4da679a 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -11,6 +11,108 @@ Rcpp::Rostream<true>&  Rcpp::Rcout = Rcpp::Rcpp_cout_get();
 Rcpp::Rostream<false>& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get();
 #endif
 
+// solveRRBLUP
+Rcpp::List solveRRBLUP(const arma::mat& y, const arma::mat& X, const arma::mat& M);
+RcppExport SEXP _AlphaSimR_solveRRBLUP(SEXP ySEXP, SEXP XSEXP, SEXP MSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< const arma::mat& >::type y(ySEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type X(XSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type M(MSEXP);
+    rcpp_result_gen = Rcpp::wrap(solveRRBLUP(y, X, M));
+    return rcpp_result_gen;
+END_RCPP
+}
+// solveRRBLUPMV
+Rcpp::List solveRRBLUPMV(const arma::mat& Y, const arma::mat& X, const arma::mat& M, int maxIter, double tol);
+RcppExport SEXP _AlphaSimR_solveRRBLUPMV(SEXP YSEXP, SEXP XSEXP, SEXP MSEXP, SEXP maxIterSEXP, SEXP tolSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< const arma::mat& >::type Y(YSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type X(XSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type M(MSEXP);
+    Rcpp::traits::input_parameter< int >::type maxIter(maxIterSEXP);
+    Rcpp::traits::input_parameter< double >::type tol(tolSEXP);
+    rcpp_result_gen = Rcpp::wrap(solveRRBLUPMV(Y, X, M, maxIter, tol));
+    return rcpp_result_gen;
+END_RCPP
+}
+// solveRRBLUPMK
+Rcpp::List solveRRBLUPMK(arma::mat& y, arma::mat& X, arma::field<arma::mat>& Mlist, int maxIter);
+RcppExport SEXP _AlphaSimR_solveRRBLUPMK(SEXP ySEXP, SEXP XSEXP, SEXP MlistSEXP, SEXP maxIterSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< arma::mat& >::type y(ySEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type X(XSEXP);
+    Rcpp::traits::input_parameter< arma::field<arma::mat>& >::type Mlist(MlistSEXP);
+    Rcpp::traits::input_parameter< int >::type maxIter(maxIterSEXP);
+    rcpp_result_gen = Rcpp::wrap(solveRRBLUPMK(y, X, Mlist, maxIter));
+    return rcpp_result_gen;
+END_RCPP
+}
+// solveRRBLUP_EM
+Rcpp::List solveRRBLUP_EM(arma::mat& Y, arma::mat& X, arma::mat& M, double Vu, double Ve, double tol, int maxIter, bool useEM);
+RcppExport SEXP _AlphaSimR_solveRRBLUP_EM(SEXP YSEXP, SEXP XSEXP, SEXP MSEXP, SEXP VuSEXP, SEXP VeSEXP, SEXP tolSEXP, SEXP maxIterSEXP, SEXP useEMSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< arma::mat& >::type Y(YSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type X(XSEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type M(MSEXP);
+    Rcpp::traits::input_parameter< double >::type Vu(VuSEXP);
+    Rcpp::traits::input_parameter< double >::type Ve(VeSEXP);
+    Rcpp::traits::input_parameter< double >::type tol(tolSEXP);
+    Rcpp::traits::input_parameter< int >::type maxIter(maxIterSEXP);
+    Rcpp::traits::input_parameter< bool >::type useEM(useEMSEXP);
+    rcpp_result_gen = Rcpp::wrap(solveRRBLUP_EM(Y, X, M, Vu, Ve, tol, maxIter, useEM));
+    return rcpp_result_gen;
+END_RCPP
+}
+// solveRRBLUP_EM2
+Rcpp::List solveRRBLUP_EM2(const arma::mat& Y, const arma::mat& X, const arma::mat& M1, const arma::mat& M2, double Vu1, double Vu2, double Ve, double tol, int maxIter, bool useEM);
+RcppExport SEXP _AlphaSimR_solveRRBLUP_EM2(SEXP YSEXP, SEXP XSEXP, SEXP M1SEXP, SEXP M2SEXP, SEXP Vu1SEXP, SEXP Vu2SEXP, SEXP VeSEXP, SEXP tolSEXP, SEXP maxIterSEXP, SEXP useEMSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< const arma::mat& >::type Y(YSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type X(XSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type M1(M1SEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type M2(M2SEXP);
+    Rcpp::traits::input_parameter< double >::type Vu1(Vu1SEXP);
+    Rcpp::traits::input_parameter< double >::type Vu2(Vu2SEXP);
+    Rcpp::traits::input_parameter< double >::type Ve(VeSEXP);
+    Rcpp::traits::input_parameter< double >::type tol(tolSEXP);
+    Rcpp::traits::input_parameter< int >::type maxIter(maxIterSEXP);
+    Rcpp::traits::input_parameter< bool >::type useEM(useEMSEXP);
+    rcpp_result_gen = Rcpp::wrap(solveRRBLUP_EM2(Y, X, M1, M2, Vu1, Vu2, Ve, tol, maxIter, useEM));
+    return rcpp_result_gen;
+END_RCPP
+}
+// solveRRBLUP_EM3
+Rcpp::List solveRRBLUP_EM3(const arma::mat& Y, const arma::mat& X, const arma::mat& M1, const arma::mat& M2, const arma::mat& M3, double Vu1, double Vu2, double Vu3, double Ve, double tol, int maxIter, bool useEM);
+RcppExport SEXP _AlphaSimR_solveRRBLUP_EM3(SEXP YSEXP, SEXP XSEXP, SEXP M1SEXP, SEXP M2SEXP, SEXP M3SEXP, SEXP Vu1SEXP, SEXP Vu2SEXP, SEXP Vu3SEXP, SEXP VeSEXP, SEXP tolSEXP, SEXP maxIterSEXP, SEXP useEMSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< const arma::mat& >::type Y(YSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type X(XSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type M1(M1SEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type M2(M2SEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type M3(M3SEXP);
+    Rcpp::traits::input_parameter< double >::type Vu1(Vu1SEXP);
+    Rcpp::traits::input_parameter< double >::type Vu2(Vu2SEXP);
+    Rcpp::traits::input_parameter< double >::type Vu3(Vu3SEXP);
+    Rcpp::traits::input_parameter< double >::type Ve(VeSEXP);
+    Rcpp::traits::input_parameter< double >::type tol(tolSEXP);
+    Rcpp::traits::input_parameter< int >::type maxIter(maxIterSEXP);
+    Rcpp::traits::input_parameter< bool >::type useEM(useEMSEXP);
+    rcpp_result_gen = Rcpp::wrap(solveRRBLUP_EM3(Y, X, M1, M2, M3, Vu1, Vu2, Vu3, Ve, tol, maxIter, useEM));
+    return rcpp_result_gen;
+END_RCPP
+}
 // callFastRRBLUP
 Rcpp::List callFastRRBLUP(arma::vec y, arma::field<arma::Cube<unsigned char> >& geno, arma::Col<int>& lociPerChr, arma::uvec lociLoc, double Vu, double Ve, arma::uword maxIter, int nThreads);
 RcppExport SEXP _AlphaSimR_callFastRRBLUP(SEXP ySEXP, SEXP genoSEXP, SEXP lociPerChrSEXP, SEXP lociLocSEXP, SEXP VuSEXP, SEXP VeSEXP, SEXP maxIterSEXP, SEXP nThreadsSEXP) {
@@ -219,6 +321,52 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// solveUVM
+Rcpp::List solveUVM(const arma::mat& y, const arma::mat& X, const arma::mat& Z, const arma::mat& K);
+RcppExport SEXP _AlphaSimR_solveUVM(SEXP ySEXP, SEXP XSEXP, SEXP ZSEXP, SEXP KSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< const arma::mat& >::type y(ySEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type X(XSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type Z(ZSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type K(KSEXP);
+    rcpp_result_gen = Rcpp::wrap(solveUVM(y, X, Z, K));
+    return rcpp_result_gen;
+END_RCPP
+}
+// solveMVM
+Rcpp::List solveMVM(const arma::mat& Y, const arma::mat& X, const arma::mat& Z, const arma::mat& K, double tol, int maxIter);
+RcppExport SEXP _AlphaSimR_solveMVM(SEXP YSEXP, SEXP XSEXP, SEXP ZSEXP, SEXP KSEXP, SEXP tolSEXP, SEXP maxIterSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< const arma::mat& >::type Y(YSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type X(XSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type Z(ZSEXP);
+    Rcpp::traits::input_parameter< const arma::mat& >::type K(KSEXP);
+    Rcpp::traits::input_parameter< double >::type tol(tolSEXP);
+    Rcpp::traits::input_parameter< int >::type maxIter(maxIterSEXP);
+    rcpp_result_gen = Rcpp::wrap(solveMVM(Y, X, Z, K, tol, maxIter));
+    return rcpp_result_gen;
+END_RCPP
+}
+// solveMKM
+Rcpp::List solveMKM(arma::mat& y, arma::mat& X, arma::field<arma::mat>& Zlist, arma::field<arma::mat>& Klist, int maxIter, double tol);
+RcppExport SEXP _AlphaSimR_solveMKM(SEXP ySEXP, SEXP XSEXP, SEXP ZlistSEXP, SEXP KlistSEXP, SEXP maxIterSEXP, SEXP tolSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< arma::mat& >::type y(ySEXP);
+    Rcpp::traits::input_parameter< arma::mat& >::type X(XSEXP);
+    Rcpp::traits::input_parameter< arma::field<arma::mat>& >::type Zlist(ZlistSEXP);
+    Rcpp::traits::input_parameter< arma::field<arma::mat>& >::type Klist(KlistSEXP);
+    Rcpp::traits::input_parameter< int >::type maxIter(maxIterSEXP);
+    Rcpp::traits::input_parameter< double >::type tol(tolSEXP);
+    rcpp_result_gen = Rcpp::wrap(solveMKM(y, X, Zlist, Klist, maxIter, tol));
+    return rcpp_result_gen;
+END_RCPP
+}
 // writeASGenotypes
 void writeASGenotypes(const arma::Cube<unsigned char>& g, const arma::field<arma::uvec>& locations, const arma::uvec& allLocations, const arma::vec& snpchips, const std::vector<std::string>& names, const char missing, const std::string fname);
 RcppExport SEXP _AlphaSimR_writeASGenotypes(SEXP gSEXP, SEXP locationsSEXP, SEXP allLocationsSEXP, SEXP snpchipsSEXP, SEXP namesSEXP, SEXP missingSEXP, SEXP fnameSEXP) {
@@ -675,6 +823,12 @@ END_RCPP
 }
 
 static const R_CallMethodDef CallEntries[] = {
+    {"_AlphaSimR_solveRRBLUP", (DL_FUNC) &_AlphaSimR_solveRRBLUP, 3},
+    {"_AlphaSimR_solveRRBLUPMV", (DL_FUNC) &_AlphaSimR_solveRRBLUPMV, 5},
+    {"_AlphaSimR_solveRRBLUPMK", (DL_FUNC) &_AlphaSimR_solveRRBLUPMK, 4},
+    {"_AlphaSimR_solveRRBLUP_EM", (DL_FUNC) &_AlphaSimR_solveRRBLUP_EM, 8},
+    {"_AlphaSimR_solveRRBLUP_EM2", (DL_FUNC) &_AlphaSimR_solveRRBLUP_EM2, 10},
+    {"_AlphaSimR_solveRRBLUP_EM3", (DL_FUNC) &_AlphaSimR_solveRRBLUP_EM3, 12},
     {"_AlphaSimR_callFastRRBLUP", (DL_FUNC) &_AlphaSimR_callFastRRBLUP, 8},
     {"_AlphaSimR_callRRBLUP", (DL_FUNC) &_AlphaSimR_callRRBLUP, 8},
     {"_AlphaSimR_callRRBLUP2", (DL_FUNC) &_AlphaSimR_callRRBLUP2, 13},
@@ -685,6 +839,9 @@ static const R_CallMethodDef CallEntries[] = {
     {"_AlphaSimR_callRRBLUP_GCA2", (DL_FUNC) &_AlphaSimR_callRRBLUP_GCA2, 14},
     {"_AlphaSimR_callRRBLUP_SCA", (DL_FUNC) &_AlphaSimR_callRRBLUP_SCA, 9},
     {"_AlphaSimR_callRRBLUP_SCA2", (DL_FUNC) &_AlphaSimR_callRRBLUP_SCA2, 15},
+    {"_AlphaSimR_solveUVM", (DL_FUNC) &_AlphaSimR_solveUVM, 4},
+    {"_AlphaSimR_solveMVM", (DL_FUNC) &_AlphaSimR_solveMVM, 6},
+    {"_AlphaSimR_solveMKM", (DL_FUNC) &_AlphaSimR_solveMKM, 6},
     {"_AlphaSimR_writeASGenotypes", (DL_FUNC) &_AlphaSimR_writeASGenotypes, 7},
     {"_AlphaSimR_writeASHaplotypes", (DL_FUNC) &_AlphaSimR_writeASHaplotypes, 7},
     {"_AlphaSimR_calcGenParam", (DL_FUNC) &_AlphaSimR_calcGenParam, 3},
diff --git a/src/calcGenParam.cpp b/src/calcGenParam.cpp
index 1cff1779..49e20434 100644
--- a/src/calcGenParam.cpp
+++ b/src/calcGenParam.cpp
@@ -40,7 +40,7 @@ Rcpp::List calcGenParamE(const Rcpp::S4& trait,
     gv_d.set_size(nInd,nThreads);
     gv_d.zeros();
   }
-  arma::vec x(ploidy+1); // Genotype dossage
+  arma::vec x(ploidy+1); // Genotype dosage
   for(arma::uword i=0; i<x.n_elem; ++i)
     x(i) = double(i);
   arma::vec xa = (x-dP/2.0)*(2.0/dP);
diff --git a/src/ibd.cpp b/src/ibd.cpp
index 245f1b78..a70e1e09 100644
--- a/src/ibd.cpp
+++ b/src/ibd.cpp
@@ -212,18 +212,18 @@ arma::Mat<int> createIbdMat(arma::field<arma::field<arma::field<arma::Mat<int> >
     for(arma::uword j=0; j<ploidy; ++j){
       arma::uword stop,start=0;
       for(arma::uword k=0; k<nChr; ++k){
-        arma::uword nSeg = ibd(i)(k)(j).n_rows;
+        arma::uword nSeg = ibd(i)(chr(k))(j).n_rows;
         if(nSeg>1){
           // First segments
           for(arma::uword l=0; l<(nSeg-1); ++l){
-            stop = start + ibd(i)(k)(j)(l+1,1) - ibd(i)(k)(j)(l,1) - 1;
-            output.col(i*ploidy+j).rows(start,stop).fill(ibd(i)(k)(j)(l,0));
+            stop = start + ibd(i)(chr(k))(j)(l+1,1) - ibd(i)(chr(k))(j)(l,1) - 1;
+            output.col(i*ploidy+j).rows(start,stop).fill(ibd(i)(chr(k))(j)(l,0));
             start = stop + 1;
           }
         }
         // Last segment
-        stop = start + nLoci(chr(k)) - ibd(i)(k)(j)(nSeg-1,1);
-        output.col(i*ploidy+j).rows(start,stop).fill(ibd(i)(k)(j)(nSeg-1,0));
+        stop = start + nLoci(chr(k)) - ibd(i)(chr(k))(j)(nSeg-1,1);
+        output.col(i*ploidy+j).rows(start,stop).fill(ibd(i)(chr(k))(j)(nSeg-1,0));
         start = stop+1;
       }
     }
diff --git a/src/misc.cpp b/src/misc.cpp
index c584948b..891f9100 100644
--- a/src/misc.cpp
+++ b/src/misc.cpp
@@ -105,19 +105,19 @@ arma::uword mapIndex(arma::uword i, arma::uword j,
 // Find row given mapping index
 // k = mapping index
 // n = dimension of matrix
-arma::uword mapRow(arma::uword k, arma::uword n){
+arma::uword mapRow(const arma::uword& k, const arma::uword& n){
   return n-2-static_cast<arma::uword>(sqrt(-8*double(k) + 4*double(n)*(double(n)-1)-7)/2-0.5);
 }
 
 // Find column given mapping index
+// row = previously determined row
 // k = mapping index
 // n = dimension of matrix
-arma::uword mapCol(arma::uword k, arma::uword n){
-  arma::uword i;
-  i = mapRow(k,n);
-  return k+i+1 - n*(n-1)/2 + (n-i)*((n-i)-1)/2;
+arma::uword mapCol(const arma::uword& row, const arma::uword& k, const arma::uword& n){
+  return k+row+1 - n*(n-1)/2 + (n-row)*((n-row)-1)/2;
 }
 
+
 // Randomly samples integers without replacement
 // n number of integers to return
 // N number of integers to sample from
@@ -261,7 +261,7 @@ arma::umat sampHalfDialComb(arma::uword nLevel, arma::uword n){
   arma::umat output(n,2);
   for(arma::uword i=0; i<n; ++i){
     output(i,0) = mapRow(samples(i),nLevel);
-    output(i,1) = mapCol(samples(i),nLevel);
+    output(i,1) = mapCol(output(i,0),samples(i),nLevel);
   }
   if(fullComb>0){
     arma::umat tmp(N*fullComb,2);
@@ -269,7 +269,7 @@ arma::umat sampHalfDialComb(arma::uword nLevel, arma::uword n){
     for(arma::uword j=0; j<(N*fullComb); ++j){
       i = j%N;
       tmp(j,0) = mapRow(i,nLevel);
-      tmp(j,1) = mapCol(i,nLevel);
+      tmp(j,1) = mapCol(tmp(j,0),i,nLevel);
     }
     output = arma::join_cols(output,tmp);
   }
@@ -302,6 +302,19 @@ double choose(double n, double k){
   return (n*choose(n-1,k-1))/k;
 }
 
+// Gets the number of available threads
+//' @title Number of available threads
+//'
+//' @description
+//' Gets the number of available threads by calling the OpenMP function
+//' \code{omp_get_max_threads()}
+//'
+//' @return integer
+//'
+//' @examples
+//' getNumThreads()
+//'
+//' @export
 // [[Rcpp::export]]
 int getNumThreads(){
 #ifdef _OPENMP
diff --git a/src/misc.h b/src/misc.h
index 9be09c20..979edd60 100644
--- a/src/misc.h
+++ b/src/misc.h
@@ -1,8 +1,11 @@
 #ifndef MISC_H
 #define MISC_H
 
+arma::uword mapRow(const arma::uword& k, const arma::uword& n);
+arma::uword mapCol(const arma::uword& row, const arma::uword& k, const arma::uword& n);
 arma::uvec sampleInt(arma::uword n, arma::uword N);
 arma::uword samplePoisson(double lambda);
+arma::umat sampHalfDialComb(arma::uword nLevel, arma::uword n);
 double choose(double n, double k);
 std::bitset<8> toBits(unsigned char byte);
 unsigned char toByte(std::bitset<8> bits);