Unbalanced Optimal Transport

Documented in plotGeodesic

#' Plot previously computed Wasserstein geodesics
#' @description This function generates either a sequence of images or a gif displaying the input optimal transport geodesic. 
#' @param Geodesic A list of measures in R^2 or R^3 corresponding to the discretized time steps of the geodesic. 
#' Each entry in the list must be one of the following: a \link[transport]{pp-object}; a \link[transport]{wpp-object}; or a list containing an entry named 'positions' that is an Mxd matrix.
#' @param method A string specifying the method used to plot each measure. The input "default" generates a scatterplot where the size of each point corresponds to its mass.
#' The input "bin" maps the data to a grid of prespecified size to plot it as an image. The input "binSplit" also maps the data to a grid, but afterwards the mass of each
#' pixel is split between all pixels in a specified range. This generates smoother output images.
#' @param images A boolean specifying whether image files should be generated. 
#' @param gif A boolean specifying whether a gif should be generated. To use this option the ImageMagick software and the \link[magick]{magick} R package need to be installed. 
#' (see https://cran.r-project.org/web/packages/magick/vignettes/intro.html for details and instructions.) 
#' @param File A string specifying the prefix of all generated files. 
#' @param resolution A vector with two elements specifying the resolution of the output images.
#' @param dotsize A positive number working as a multiplier on the size of the dots in the output generated by the default method.
#' @param gridsize A vector with two elements specifying the size of the grid used for "bin" and "binSplit". 
#' @param out.col A colour vector, specifying the colour scheme used in the call of the image function when plotting the output of the "bin" and "binSplit" methods.
#' See the documentation of the image function for more details.
#' @param splitrange A vector of two positive integers specifying the number of pixels the mass of a point is shared with in the "binSplit" method. 
#' The first entry controls the horizontal direction and the second controls the vertical direction.
#' @param turn A boolean specifying whether the image should be rotated to account for the output of the image function.
#' @param phi A number specifying the viewing direction in the three dimensional method. It gives the colatitude of the plot.
#' @param theta A number specifying the viewing direction in the three dimensional setting. It gives the azimuthal direction of the plot.
#' @param fps A positive number specifying the number of frames per second in the output gif.  The default adjusts the fps to output a gif file of length 10 seconds.
#' @return This function does not provide any return value, but instead generates output files in the current working directory.
#' @examples 
#' #2D-Example:
#' library(transport)
#' set.seed(420)
#' N<-2
#' supp.size<-10^2
#' L<-sqrt(supp.size)
#' d<-2
#' data.list<-vector("list",N)
#' image.list<-vector("list",N)
#' for (i in 1:N){
#'   t.vec<-seq(0,2*pi,length.out=supp.size)
#'   pos<-cbind(cos(t.vec)*runif(1,0.2,1),sin(t.vec)*runif(1,0.2,1))
#'   theta<-runif(1,0,2*pi)
#'   rotation<-matrix(c(cos(theta),sin(theta),-1*sin(theta),cos(theta)),2,2)
#'   pos<-pos%*%rotation
#'   pos<-pos+1
#'   pos<-pos/2
#'   W<-rep(1/supp.size,supp.size)
#'   data.list[[i]]<-transport::wpp(pos,W)
#' }
#' Geo<-geodesic_pos(data.list[[1]],data.list[[2]],2,seq(0,1,0.1))
#' ## Set the image and/or gif flags to TRUE to run the example. 
#' ## CRAN policy prevents examples from generating files in the working directory,
#' ## so this had to be disabled.
#' plotGeodesic(Geo,File="TestGeodesicDefault",images=FALSE,gif=FALSE)
#' plotGeodesic(Geo,method="bin",File="TestGeodesicDefaultBin",images=FALSE,gif=FALSE)
#' plotGeodesic(Geo,method="binSplit",File="TestGeodesicDefaultBinSPlit",
#' images=FALSE,gif=FALSE)
#' \donttest{
#' #3D-Example:
#' #Functions to build the example measures
#' gen_torus<-function(M,R,r){
#'  theta<-seq(0,2*pi,length.out=M)
#'  phi<-seq(0,2*pi,length.out=M)
#'  G<-expand.grid(theta,phi)
#'  x<-(R+r*cos(G[,1]))*cos(G[,2])
#'  y<-(R+r*cos(G[,1]))*sin(G[,2])
#'  z<-r*sin(G[,1])
#'  return(cbind(x,y,z))
#' }
#' sq_norm<-function(v){
#'  return(sqrt(sum(v^2)))
#' }
#' normalize<-function(v){
#'  return(v/(sq_norm(v)))
#' }
#' rotate3D<-function(pos,axis,angle){
#'  R<-matrix(0,3,3)
#'  R[1]<-cos(angle)+axis[1]^2*(1-cos(angle))
#'  R[2]<-axis[1]*axis[2]*(1-cos(angle))+axis[3]*sin(angle)
#'  R[3]<-axis[3]*axis[1]*(1-cos(angle))-axis[2]*sin(angle)
#'  R[4]<-axis[1]*axis[2]*(1-cos(angle))-axis[3]*sin(angle)
#'  R[5]<-cos(angle)+axis[2]^2*(1-cos(angle))
#'  R[6]<-axis[2]*axis[3]*(1-cos(angle))+axis[1]*sin(angle)
#'  R[7]<-axis[1]*axis[3]*(1-cos(angle))+axis[2]*sin(angle)
#'  R[8]<-axis[2]*axis[2]*(1-cos(angle))-axis[1]*sin(angle)
#'  R[9]<-cos(angle)+axis[3]^2*(1-cos(angle))
#'  return(t(diag(c(2,3,1))%*%(R%*%t(pos))))
#' }
#' ## Example
#' set.seed(123)
#' M<-40
#' U<-runif(1,0.5,1)
#' Torus<-gen_torus(M,U,min(U/2,runif(1)))
#' v<-normalize(runif(3))
#' Torus<-rotate3D(Torus,v,runif(1,0,2*pi))
#' Torus1<-Torus%*%diag(runif(3,1,3))
#' U<-runif(1,0.5,1)
#' Torus<-gen_torus(M,U,min(U/2,runif(1)))
#' v<-normalize(runif(3))
#' Torus<-rotate3D(Torus,v,runif(1,0,2*pi))
#' Torus2<-Torus%*%diag(runif(3,1,3))
#' L<-length(Torus)/3
#' Torus1<-transport::wpp(Torus1,rep(1/L,L))
#' Torus2<-transport::wpp(Torus2,rep(1/L,L))
#' geo<-geodesic_pos(Torus1,Torus2,p=2,seq(0,1,0.1))
#' ## Set the image and/or gif flags to TRUE to run the example. 
#' ## CRAN policy prevents examples from generating files in the working directory,
#' ## so this had to be disabled.
#' plotGeodesic(geo,File="3dGeodesic",images=FALSE,gif=FALSE)
#' }
#' @export
plotGeodesic<-function(Geodesic,method="default",images=FALSE,gif=FALSE,File="Geodesic",resolution=c(400,400),dotsize=2,gridsize=c(100,100),out.col= grey(0:1000/1000),splitrange=c(1,1),turn=FALSE,phi=40,theta=40,fps=NULL){
  L<-length(Geodesic)
  numbers<-numbering_gen(L)
  types<-lapply(Geodesic,type_check)
  Geodesic<-mapply(process_data,Geodesic,types,SIMPLIFY = FALSE)
  d<-dim(Geodesic[[1]]$positions)[2]
  if (gif){
    if (!requireNamespace("magick", quietly = TRUE)) {
      stop("The magick package is required to generate gifs. Install it or disable the gif option to proceed.")
    }
  }
  if (d>3){
    stop("This function does not support dimensions larger than 3.")
  }
  if (gif & (is.null(fps))){
    fps<-round(L/10)
  }
  if (!is.null(fps)){
    fps<-round(100/fps)
  }
  if (d==2){
    if (images | gif){
      FileOld<-File
      if (!images){
        dirname <- "temp"
        dirno <- 1
        while (!dir.create(dirname, showWarnings = FALSE)) {
          dirname <- paste("temp", dirno, sep = "")
          dirno <- dirno + 1
        }
        File<-paste(dirname,"/",File,sep="")
      }
      if (method=="default"){
        for (k in 1:L){
          name<-paste(File,numbers[k],".png",sep="")
          grDevices::png(name,resolution[1],resolution[2])
          graphics::plot(transport::wpp(Geodesic[[k]]$position,Geodesic[[k]]$weights),xlim=c(0,1),ylim=c(0,1),axes = FALSE,cex=dotsize)
          grDevices::dev.off()
        }
      }
      if (method=="bin"){
        Geo.array<-array(0,c(gridsize[1],gridsize[2],L))
        for (k in 1:L){
          Geo.array[,,k]<-bin2d(Geodesic[[k]]$positions,Geodesic[[k]]$weights,gridsize,turn=turn)
        }
        Geo.min<-min(Geo.array)
        Geo.max<-max(Geo.array)
        for (k in 1:L){
          name<-paste(File,numbers[k],".png",sep="")
          grDevices::png(name,resolution[1],resolution[2])
          graphics::image(Geo.array[, , k], col = out.col, axes = FALSE, 
                zlim = c(Geo.min,Geo.max))
          grDevices::dev.off()
        }
      }
      if (method=="binSplit"){
        Geo.array<-array(0,c(gridsize[1],gridsize[2],L))
        for (k in 1:L){
          I<-bin2d(Geodesic[[k]]$positions,Geodesic[[k]]$weights,gridsize,turn=turn)
          I<-smear(I,splitrange[1],splitrange[2])
          Geo.array[,,k]<-I
        }
        Geo.min<-min(Geo.array)
        Geo.max<-max(Geo.array)
        for (k in 1:L){
          name<-paste(File,numbers[k],".png",sep="")
          grDevices::png(name,resolution[1],resolution[2])
          graphics::image(Geo.array[, , k], col = out.col, axes = FALSE, 
                zlim = c(Geo.min,Geo.max))
          grDevices::dev.off()
        }
      }
      if (!images){
        FileIM<-File
        File<-FileOld
      }
    }
  }
  if (d==3){
    if (images | gif){
      FileOld<-File
      if (!images){
        dirname <- "temp"
        dirno <- 1
        while (!dir.create(dirname, showWarnings = FALSE)) {
          dirname <- paste("temp", dirno, sep = "")
          dirno <- dirno + 1
        }
        File<-paste(dirname,"/",File,sep="")
      }
      tmp<-Geodesic[[1]]$weights
      tmp<-tmp/sum(tmp)
      if (!all.equal(tmp,rep(1/length(tmp),length(tmp)))){
        warning("Weighted Point Pattern are unsupported in three dimensions. The measures are assumed to be uniformly weighted.")
      }
      for (k in 1:L){
        name<-paste(File,numbers[k],".png",sep="")
        grDevices::png(name,resolution[1],resolution[2])
        plot3D::scatter3D(Geodesic[[k]]$positions[,1],Geodesic[[k]]$positions[,2],Geodesic[[k]]$positions[,3],box=FALSE,colkey=FALSE)
        grDevices::dev.off()
      }
      if (!images){
        FileIM<-File
        File<-FileOld
      }
    }
  } 
  if (gif==TRUE){
    if (images==FALSE){
      files<-list.files(path=paste(getwd(),"/temp",sep=""),pattern=File)
      for (k in 1:length(files)){
        files[k]<-paste("temp/",files[k],sep="")
      }
    }
    else{
      files<-list.files(pattern=File)
    }
    IM<-magick::image_read(files)
    G<-magick::image_animate(IM,delay=fps)
    G.name<-paste(File,".gif",sep="")
    magick::image_write(G, G.name)
    if (!images){
      unlink(dirname, recursive = TRUE)
    }
  }
}