R/cost_f.R

In jwoollia/DenHaag: Simulation of Populations with Managed Inbreeding

#' @title Recommends contributions based on relative costs of gain and inbreeding.
#' @description Carries out the standard penalised contributions algorithm, including constraints on totals contributions
#'    of males and females, and non-negative contributions. The approximate solution is obtained by the algorithm of
#'    Meuwissen. It prints out the recommendations after scaling the contributions to give the required total offspring.
#'    The function checks whether the number of parents to be considered exceeds the population, whether there are both
#'    males and female parents available. If these checks fail, error messages are printed.
#' @param ddf The data frame of 9 columns holding the population data.
#' @param tp The number of parents to be considered for contributions, assumed to be the most recent.
#' @param tc The number of offspring to be created.
#' @param cval The relative cost of inbreeding which needs to be supplied by the user.
#' @return A logical indicating the success of the algorithm or otherwise.
#' @export
#
cost_f <- function(ddf,tp,tc,cval)
{
tn <- nrow(ddf)
if(tn<tp)
  {print("Number of parents requested exceeds population!")
  return(FALSE)
  }
if( !(sum(ddf$sex==1)>0))
  {print("No males among the parents!")
  return(FALSE)
  }
if( !(sum(ddf$sex==2)>0))
  {print("No females among the parents!")
  return(FALSE)
  }
print(paste0("Recommendations for ",tc, " offspring from the most recent ",tp," parents"))
# cand is list of candidates
# q is matrix of sexes of candidates
# ebv is ebv of candidates
# nrm contains their numerator relationships
cand <- ddf[ddf$id>tn-tp,"id"]
q <- matrix(nrow=tp,ncol=2)
q[,1] <- as.double((ddf$sex[cand] < 1.5))
q[,2] <- as.double((ddf$sex[cand] > 1.5))
ebv <- ddf[cand,"ebv"]
nrm <- a_mat(ddf[c(1,6,7)])
nrm <- nrm[cand,cand]
ind <- c(1:tp)
#
nn <- 0
tt <- tp
while (!(nn==tt))
  {nn <- tt
  ai <- solve(nrm[ind,ind])
  qq <- q[ind,]
  ebvv <- ebv[ind]
  qaqi <- solve(t(qq)%*%ai%*%qq)
  lgm0 <- cval
  lgm1 <- qaqi%*% (t(qq)%*%ai%*%ebvv-c(1,1)*lgm0)
  cc <- ai%*% (ebvv-qq%*%lgm1)
  cc <- cc /matrix(2.0*lgm0,nrow=nrow(cc),ncol=ncol(cc))
  # check for negative c elements
  tt <- 0
  for (ii in 1:nn)
    {if (cc[ii]>=0)
      {tt <- tt+1
      ind[tt] <- ind[ii]
      }
    }
  ind <- ind[1:tt]
}
gc <- 0.5*(t(cc)%*%nrm[ind,ind]%*%cc)
ccdf <- data.frame("id"=cand[ind],"sex"=ddf$sex[cand[ind]],"ebv"=ebv[ind],"c"=cc,"noff"= 2.0*tc*cc)
print(paste0("Recommendations have group coancestry ",gc))
print(ccdf,row.names=FALSE)
return(TRUE)
}