Temporary example for testing optim replacement (KMbranch)

.Rbuildignore

@@ -8,3 +8,5 @@
 ^CRAN-SUBMISSION$
 ^doc$
 ^Meta$
+^examples/temp/
+

DESCRIPTION

@@ -27,7 +27,7 @@ Description: Provides users with an EZ-to-use platform for representing data
 License: MIT + file LICENSE
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.2
+RoxygenNote: 7.3.3
 VignetteBuilder: knitr
 Config/testthat/edition: 3
 Suggests: 
@@ -47,7 +47,10 @@ Imports:
     graphics,
     grDevices,
     plotrix,
+    Rcpp,
     splines,
     stats,
     withr
 NeedsCompilation: yes
+LinkingTo: 
+    Rcpp

NAMESPACE

@@ -43,6 +43,7 @@ export(rotate)
 export(samples)
 export(sqrtManhattan)
 export(translate_axes)
+importFrom(Rcpp,sourceCpp)
 importFrom(grDevices,chull)
 importFrom(grDevices,colorRampPalette)
 importFrom(grDevices,grey)

R/RcppExports.R

@@ -0,0 +1,11 @@
+# Generated by using Rcpp::compileAttributes() -> do not edit by hand
+# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+
+alfunc <- function(BVEC, X, Y, M, MU, LAMBDA, CONST1, CONST2, U, V) {
+    .Call(`_biplotEZ_alfunc`, BVEC, X, Y, M, MU, LAMBDA, CONST1, CONST2, U, V)
+}
+
+optimize_spline <- function(BVEC, X, Y, M, MU, LAMBDA, CONST1, CONST2, U, V, TAU, FTOL, TINY, ITMAX) {
+    .Call(`_biplotEZ_optimize_spline`, BVEC, X, Y, M, MU, LAMBDA, CONST1, CONST2, U, V, TAU, FTOL, TINY, ITMAX)
+}
+

R/biplotEZ-package.R

@@ -54,5 +54,7 @@
 #' @importFrom graphics points
 #' @importFrom graphics text
 #' @importFrom grDevices colorRampPalette
+#' @importFrom Rcpp sourceCpp
+#' @useDynLib biplotEZ, .registration = TRUE
 ## usethis namespace: end
 NULL

R/plot2D.R

@@ -1087,13 +1087,189 @@
 #' @useDynLib biplotEZ, .registration = TRUE
 #'
 #' @noRd
+
+
+
+# biplot.spline.axis <- function(j, X, Y, means, sd, 
+#                                n.int, spline.control, dmeth=0, ... )
+# {
+#   n <- nrow(X)
+#   p <- ncol(X)
+#   if (n > 103)
+#     {  
+#       my.sample <- sample (1:n, size=103, replace=F)
+#       X <- X[my.sample,]
+#       Y <- Y[my.sample,]
+#       n <- nrow(X)
+#     }
+
+#   tau <- spline.control$tau
+#   nmu <- spline.control$nmu
+#   u <- spline.control$u
+#   v <- spline.control$v
+#   lambda <- spline.control$lambda
+#   smallsigma <- spline.control$smallsigma
+#   bigsigma <- spline.control$bigsigma
+#   gamma <- spline.control$gamma
+#   bigsigmaactivate <- spline.control$bigsigmaactivate
+#   eps <- spline.control$eps
+#   tiny <- spline.control$tiny
+#   itmax <- spline.control$itmax
+#   ftol <- spline.control$ftol
+
+#   cat ("Calculating spline axis for variable", j, "\n")
+#   if(dmeth==1) stop("dmeth should be equal to zero or integer greater than 1 \n")  
+#   Ytilde <- scale(scale(Y, center=FALSE, scale=1/sd), center=-1*means, scale=FALSE)
+
+#   ytilde <- Ytilde[,j]
+#   mutilde <- seq(from=min(ytilde),to=max(ytilde),length.out=nmu)
+#   y <- Y[,j]
+#   rangey <- max(y)-min(y)
+#   mu <- seq(from=min(y)-.3*rangey,to=max(y)+.3*rangey,length.out=nmu)
+#   markers <- (pretty(ytilde)-means[j])/sd[j]
+#   mu <- sort(c(mu,markers))
+#   mu <- unique(mu)
+#   nmu <- length(mu)
+
+#   if (v>0)
+#   {
+#     knots <- seq.int(from=0,to=1,length.out=v+2)[-c(1,v+2)]
+#     knots <- stats::quantile(y,knots)
+#     M <- splines::bs(mu,knots=knots,degree=u,intercept=FALSE)
+#   } else M <- splines::bs(mu,df=u+v,degree=u,intercept=FALSE)
+#   M <- scale(M,scale=FALSE,center=M[which.min(abs(mu)),]) # To ensure that the spline passes through the origin at the calibration which represents the mean of the variable
+#   Breg <- t(solve(t(X)%*%X)%*%t(X)%*%y)
+#   Zreg <- mu%*%Breg/sum(Breg^2)
+#   Bvec <- as.vector(solve(t(M)%*%M)%*%t(M)%*%Zreg)  # Closest to regression biplot
+
+#   const1 <- sum(y^2)
+#   const2 <- sum(X^2)/(n*p)
+#   TotalNumberOfLossFunctionCalls <- 0
+
+#   optimtouse <- function(Bvec)
+#   {
+#     timetemp <- proc.time()[3]
+#     LOSS <- 1.0
+#     LOSS1 <- 1.0
+#     Ind <- rep(1,n)
+#     pred <- rep(0,nmu)
+#     deltmp <- 0
+#     tau <- tau
+#     #.5 # the choice of tau seems to affect perfomance quite substantially.
+#     # tau is used to specify the points on the inital simplex.
+#     Ay <- rep(0,(u+v)*p+1)
+#     TEMPVK <- rep(0,(u+v)*p)
+#     iter1 <- 0
+#     iter <- 0
+#     ERRO <- 0
+
+#     # Prepare for Fortran subroutine
+#     storage.mode(X) <- "double"
+#     storage.mode(Ind) <- "integer"
+#     storage.mode(mu) <- "double"
+#     storage.mode(pred) <- "double"
+#     storage.mode(y) <- "double"
+#     storage.mode(M) <- "double"
+#     storage.mode(Bvec) <- "double"
+#     storage.mode(Ay) <- "double"
+#     storage.mode(TEMPVK) <- "double"
+
+
+#     returned_data <-.Fortran('L',LOSS=as.double(LOSS),X=X,n=as.integer(n),p=as.integer(p),nmu=as.integer(nmu),Ind=Ind,
+#                             mu=mu,pred=pred,lambda=as.double(lambda),y=y,const1=as.double(const1),const2=as.double(const2),u=as.integer(u),
+#                             v=as.integer(v),M=M,Bvec=Bvec,tau=as.double(tau),Ay=Ay,TEMPVEK=TEMPVK,iter=as.integer(iter),
+#                             ftol=as.double(ftol),LOSS1=as.double(LOSS1),iter1=as.integer(iter1),fout = as.integer(ERRO),
+#                             const3=as.double(tiny), itmax=as.integer(itmax))
+#     if(returned_data$fout > 0)
+#     {
+#       cat("Fout is: ", returned_data$fout, "\n")
+#       warning("Increase itmax for Fortran \n")
+#     }
+
+#     B <- matrix(returned_data$Bvec,ncol=p)
+#     Z <- M%*%B 
+
+#     aa <- list(BestValue=returned_data$LOSS,BestSolution=returned_data$Bvec,ConvergenceCode=returned_data$fout, iter1=returned_data$iter1,
+#              iter=returned_data$iter,TimeTaken=proc.time()[3]-timetemp)
+#     aa
+#   }
+#   EuclidDist2 <- function (X, Y) 
+#   {
+#     n <- nrow(X)
+#     m <- nrow(Y)
+#     bx <- rowSums(X^2)
+#     by <- rowSums(Y^2)
+#     outer(bx, by, FUN = "+") - 2 * X %*% t(Y)
+#   }
+
+#   ### Variable initialisation
+#   outBestValues <- rep(NA,gamma+1)
+#   outBestSolutions <- matrix(nrow=2*(u+v),ncol=gamma+1)
+#   outTimeTaken <- rep(NA,gamma+1) # Is made one element longer at each iteration.
+#   BestSolutionsFrequency <- rep(NA,gamma+1)
+#   BestSolutionsIndices <- rep(NA,gamma+1) # Is made one element longer at each iteration.
+#   SquaredDistancesBetweenBestSolutions <- matrix(nrow=gamma+1,ncol=gamma+1)
+
+#   ### Initial coefficients closest to regression biplot
+#   temp <- optimtouse(Bvec)
+#   outBestValues[1] <- temp$BestValue
+#   outBestSolutions[,1] <- temp$BestSolution
+#   outTimeTaken[1] <- temp$TimeTaken
+#   BestSolutionsFrequency[1] <- 1
+#   BestSolutionsIndices[1] <- 1
+#   DistinctSolutions <- 1
+#   PreviousBestSolution <- NA
+#   nSameSolutionConsecutively <- 0
+#   BigSigmaActivations <- NULL
+
+#   test.iter <- temp$iter
+#   test.iter1 <- temp$iter1
+
+#   ### Last best coefficients perturbed
+#   for (gammacounter in 2:(gamma+1))
+#   {
+#     if (nSameSolutionConsecutively>=bigsigmaactivate)
+#     {
+#       temp <- optimtouse(outBestSolutions[,which.min(outBestValues)]+stats::rnorm((u+v)*2,mean=0,sd=bigsigma))
+#       BigSigmaActivations <- c(BigSigmaActivations,gammacounter)
+#     }
+#     else temp <- optimtouse(outBestSolutions[,which.min(outBestValues)]+stats::rnorm((u+v)*2,mean=0,sd=smallsigma))
+#     outTimeTaken[gammacounter] <- temp$TimeTaken
+#     tempSquaredDistances <- EuclidDist2(matrix(temp$BestSolution,nrow=1),t(outBestSolutions[,1:DistinctSolutions]))
+#     if (any(tempSquaredDistances<eps))
+#     {
+#       BestSolutionsIndices[gammacounter] <- tempAA<-which.min(tempSquaredDistances)
+#       BestSolutionsFrequency[tempAA] <- BestSolutionsFrequency[tempAA]+1
+#       if (!is.na(PreviousBestSolution) && tempAA==PreviousBestSolution) nSameSolutionConsecutively<-nSameSolutionConsecutively+1
+#       else
+#       {
+#         PreviousBestSolution <- tempAA
+#         nSameSolutionConsecutively <- 0
+#       }
+#     }
+#     else
+#     {
+#       DistinctSolutions <- DistinctSolutions+1
+#       outBestValues[DistinctSolutions] <- temp$BestValue
+#       outBestSolutions[,DistinctSolutions] <- temp$BestSolution
+#       BestSolutionsFrequency[DistinctSolutions] <- 1
+#       BestSolutionsIndices[gammacounter] <- DistinctSolutions
+#       SquaredDistancesBetweenBestSolutions[1:(DistinctSolutions-1),DistinctSolutions]<-tempSquaredDistances
+#       nSameSolutionConsecutively <- 0
+#     }
+#   }
+#   axis.points <- cbind(M%*%matrix(outBestSolutions[,which.min(outBestValues)],ncol=2), mu, 0)
+
+#   for (i in 1:nrow(axis.points)) if (any(zapsmall(axis.points[i,3] - markers) == 0)) axis.points[i, 4] <- 1
+#   axis.points[,3] <- axis.points[,3]*sd[j] + means[j]
+#   axis.points
+# }
+
 biplot.spline.axis <- function(j, X, Y, means, sd, 
-                               n.int, spline.control, dmeth=0, ... )
-{
+                               n.int, spline.control, dmeth=0, ... ){
   n <- nrow(X)
   p <- ncol(X)
-  if (n > 103)
-    {  
+  if (n > 103){  
       my.sample <- sample (1:n, size=103, replace=F)
       X <- X[my.sample,]
       Y <- Y[my.sample,]
@@ -1128,8 +1304,7 @@ biplot.spline.axis <- function(j, X, Y, means, sd,
   mu <- unique(mu)
   nmu <- length(mu)
 
-  if (v>0)
-  {
+  if (v>0){
     knots <- seq.int(from=0,to=1,length.out=v+2)[-c(1,v+2)]
     knots <- stats::quantile(y,knots)
     M <- splines::bs(mu,knots=knots,degree=u,intercept=FALSE)
@@ -1142,56 +1317,44 @@ biplot.spline.axis <- function(j, X, Y, means, sd,
   const1 <- sum(y^2)
   const2 <- sum(X^2)/(n*p)
   TotalNumberOfLossFunctionCalls <- 0
-
-  optimtouse <- function(Bvec)
-  {
+  optimtouse <- function(Bvec) {
     timetemp <- proc.time()[3]
-    LOSS <- 1.0
-    LOSS1 <- 1.0
-    Ind <- rep(1,n)
-    pred <- rep(0,nmu)
-    deltmp <- 0
-    tau <- tau
-    #.5 # the choice of tau seems to affect perfomance quite substantially.
-    # tau is used to specify the points on the inital simplex.
-    Ay <- rep(0,(u+v)*p+1)
-    TEMPVK <- rep(0,(u+v)*p)
-    iter1 <- 0
-    iter <- 0
-    ERRO <- 0
-
-    # Prepare for Fortran subroutine
-    storage.mode(X) <- "double"
-    storage.mode(Ind) <- "integer"
-    storage.mode(mu) <- "double"
-    storage.mode(pred) <- "double"
-    storage.mode(y) <- "double"
-    storage.mode(M) <- "double"
-    storage.mode(Bvec) <- "double"
-    storage.mode(Ay) <- "double"
-    storage.mode(TEMPVK) <- "double"
 
+      # Use Rcpp optimize_spline function
+    returned_data <- optimize_spline(
+      BVEC = Bvec, 
+      X = X, 
+      Y = y, 
+      M = M, 
+      MU = mu,
+      LAMBDA = lambda, 
+      CONST1 = const1, 
+      CONST2 = const2,
+      U = u, 
+      V = v, 
+      TAU = tau, 
+      FTOL = ftol, 
+      TINY = tiny, 
+      ITMAX = itmax
+    )
 
-    returned_data <-.Fortran('L',LOSS=as.double(LOSS),X=X,n=as.integer(n),p=as.integer(p),nmu=as.integer(nmu),Ind=Ind,
-                            mu=mu,pred=pred,lambda=as.double(lambda),y=y,const1=as.double(const1),const2=as.double(const2),u=as.integer(u),
-                            v=as.integer(v),M=M,Bvec=Bvec,tau=as.double(tau),Ay=Ay,TEMPVEK=TEMPVK,iter=as.integer(iter),
-                            ftol=as.double(ftol),LOSS1=as.double(LOSS1),iter1=as.integer(iter1),fout = as.integer(ERRO),
-                            const3=as.double(tiny), itmax=as.integer(itmax))
-    if(returned_data$fout > 0)
-    {
-      cat("Fout is: ", returned_data$fout, "\n")
-      warning("Increase itmax for Fortran \n")
+    if(returned_data$ERRO > 0) {
+      cat("  Warning: Error code =", returned_data$ERRO, "\n")
     }
 
-    B <- matrix(returned_data$Bvec,ncol=p)
-    Z <- M%*%B 
+    aa <- list(
+      BestValue = returned_data$LOSS,
+      BestSolution = returned_data$BVEC,
+      ConvergenceCode = returned_data$ERRO,
+      iter1 = returned_data$ITER1,
+      iter = returned_data$ITER,
+      TimeTaken = proc.time()[3] - timetemp
+    )
 
-    aa <- list(BestValue=returned_data$LOSS,BestSolution=returned_data$Bvec,ConvergenceCode=returned_data$fout, iter1=returned_data$iter1,
-             iter=returned_data$iter,TimeTaken=proc.time()[3]-timetemp)
     aa
   }
-  EuclidDist2 <- function (X, Y) 
-  {
+
+  EuclidDist2 <- function (X, Y) {
     n <- nrow(X)
     m <- nrow(Y)
     bx <- rowSums(X^2)
@@ -1223,29 +1386,24 @@ biplot.spline.axis <- function(j, X, Y, means, sd,
   test.iter1 <- temp$iter1
 
   ### Last best coefficients perturbed
-  for (gammacounter in 2:(gamma+1))
-  {
-    if (nSameSolutionConsecutively>=bigsigmaactivate)
-    {
+  for (gammacounter in 2:(gamma+1)){
+    if (nSameSolutionConsecutively>=bigsigmaactivate){
       temp <- optimtouse(outBestSolutions[,which.min(outBestValues)]+stats::rnorm((u+v)*2,mean=0,sd=bigsigma))
       BigSigmaActivations <- c(BigSigmaActivations,gammacounter)
     }
     else temp <- optimtouse(outBestSolutions[,which.min(outBestValues)]+stats::rnorm((u+v)*2,mean=0,sd=smallsigma))
     outTimeTaken[gammacounter] <- temp$TimeTaken
     tempSquaredDistances <- EuclidDist2(matrix(temp$BestSolution,nrow=1),t(outBestSolutions[,1:DistinctSolutions]))
-    if (any(tempSquaredDistances<eps))
-    {
+    if (any(tempSquaredDistances<eps)){
       BestSolutionsIndices[gammacounter] <- tempAA<-which.min(tempSquaredDistances)
       BestSolutionsFrequency[tempAA] <- BestSolutionsFrequency[tempAA]+1
       if (!is.na(PreviousBestSolution) && tempAA==PreviousBestSolution) nSameSolutionConsecutively<-nSameSolutionConsecutively+1
-      else
-      {
+      else{
         PreviousBestSolution <- tempAA
         nSameSolutionConsecutively <- 0
       }
     }
-    else
-    {
+    else{
       DistinctSolutions <- DistinctSolutions+1
       outBestValues[DistinctSolutions] <- temp$BestValue
       outBestSolutions[,DistinctSolutions] <- temp$BestSolution
@@ -1263,6 +1421,9 @@ biplot.spline.axis <- function(j, X, Y, means, sd,
 }
 
 
+
+
+
 #' Plot nonlinear axes on biplots
 #'
 #' @param z.axes list containing all the info to draw axis.