content: {- TooManyCells.Spatial.Relationships
  Gregory W. Schwartz
  
  Collects the functions pertaining to quantifying relationships between points.
  -}
  
  {-# LANGUAGE QuasiQuotes #-}
  {-# LANGUAGE BangPatterns #-}
  {-# LANGUAGE OverloadedStrings #-}
  
  module TooManyCells.Spatial.Relationships
      ( spatialRelationshipsR
      ) where
  
  -- Remote
  import BirchBeer.Types (LabelMap (..))
  import Data.Bool (bool)
  import Data.Maybe (isJust)
  import Language.R as R
  import Language.R.QQ (r)
  
  -- Local
  import TooManyCells.File.Types (OutputDirectory (..))
  import TooManyCells.Matrix.Types (SingleCells (..), ProjectionMap (..), X (..), Y (..))
  import TooManyCells.Matrix.Utility (sparseMatToSparseRMat)
  import TooManyCells.Spatial.Utility (scToRPat)
  import qualified TooManyCells.Spatial.Types as Too
  
  spatialRelationshipsR :: OutputDirectory
                        -> Too.PCFCrossFlag
                        -> ProjectionMap
                        -> Maybe LabelMap
                        -> SingleCells
                        -> Maybe Too.TopDistances
                        -> [Too.Mark]
                        -> IO ()
  spatialRelationshipsR (OutputDirectory outDir) pcfCrossFlag' pm lm sc td marks = R.runRegion $ do
    let isLabelMarks = bool 0 1 . isJust $ lm :: Double
        pcfCrossFlagDouble = bool 0 1 $ Too.unPCFCrossFlag pcfCrossFlag' :: Double
        (tVal, tQFlag) = case td of
                          Nothing -> (10, 1 :: Double)
                          Just (Too.TopQuantile x) -> (x :: Double, 1)
                          Just (Too.TopDistance x) -> (x :: Double, 0)
  
    pat <- scToRPat lm pm sc marks
  
    [r|
      suppressMessages(library(spatstat))
      suppressMessages(library(reshape2))
      suppressMessages(library(plyr))
  
      # Find index of first switch from above to below or equal value.
      findNegSwap = function(x, xs) {
        if(all(xs$obs > x)) {
          return(xs$r[length(xs$r)])
        } else if(all(xs$obs <= x)) {
          return(xs$r[1])
        } else {
          for (i in c(1:length(xs$r))) {
            if (i == length(xs$r)) {
              return(xs$r[i])
            } else if ((xs$obs[i] > x) && (xs$obs[i + 1] <= x)) {
              return(xs$r[i + 1])
            }
          }
        }
      }
  
      # Find index of first switch from below or equal value to above.
      findPosSwap = function(x, xs) {
        if (all(xs$obs <= x)) {
          return(xs$r[length(xs$r)])
        } else if (all(xs$obs > x)) {
          return(xs$r[1])
        } else {
          for (i in c(1:length(xs$r))) {
            if (i == length(xs$r)) {
              return(xs$r[i])
            } else if (xs$obs[i] <= x && xs$obs[i + 1] > x) {
              return(xs$r[i + 1])
            }
          }
        }
      }
  
      # Find the longest length of a continuously positive (or negative / equal to) stretch.
      findLongestLength = function(posFlag, x, xs) {
        if(posFlag) {
          lengthInfo = rle(xs$obs > x)
        } else {
          lengthInfo = rle(xs$obs <= x)
        }
        stretches = lengthInfo$lengths[which(lengthInfo$values)]
        if(length(stretches) == 0) {
          return(0)
        } else {
          return(xs$r[max(stretches)] - xs$r[1])
        }
      }
  
      # Find the first r value of the maximum value.
      findMaxPos = function(xs) {
        return(xs$r[min(which(xs$obs == max(xs$obs)))])
      }
  
      # Find the first r value of the maximum value with maximum value.
      findMaxPosWithVal = function(as) {
        # Do not use 0 values
        xs = as[-1,]
        r = findMaxPos(xs)
        val = max(xs$obs)
  
        return(val / r)
      }
  
      # Find the first r value of the minimum value.
      findMinPos = function(xs) {
        return(xs$r[min(which(xs$obs == min(xs$obs)))])
      }
  
      topDistances = function(xs) {
        if(tQFlag_hs == 1) {
          return(xs$obs[1:round(length(xs$obs) / tVal_hs)])
        } else {
          return(xs$obs[1:max(which(xs$r <= tVal_hs))])
        }
      }
  
      # Find statistics for a single curve
      getCurveStats = function (outFolder, df) {
        statsDf = data.frame(meanCorr = mean(df$obs))
        statsDf$maxCorr = max(df$obs)
        statsDf$minCorr = min(df$obs)
        statsDf$meanCorr = mean(df$obs)
        statsDf$topMaxCorr = max(topDistances(df))
        statsDf$topMeanCorr = mean(topDistances(df))
        statsDf$topMinCorr = min(topDistances(df))
        statsDf$negSwap = findNegSwap(1, df)
        statsDf$posSwap = findPosSwap(1, df)
        statsDf$longestPosLength = findLongestLength(TRUE, 1, df)
        statsDf$longestNegLength = findLongestLength(FALSE, 1, df)
        statsDf$maxPosWithVal = findMaxPosWithVal(df)
        statsDf$logMaxPosWithVal = log(findMaxPosWithVal(df))
        statsDf$maxPos = findMaxPos(df)
        statsDf$minPos = findMinPos(df)
        statsDf$label = basename(outFolder)
        statsDf$topVal = tVal_hs
        statsDf$topQuantile = tQFlag_hs
  
        return(statsDf)
      }
  
      crossCorrCurveStats = function (outFolder) {
  
        # Cross correlation function (several columns of marks).
        if(pcfCrossFlagDouble_hs == 1) {
          crossFn = envelope(pat_hs, pcfcross, funargs=marks)
          crossFn$obs = crossFn$iso
        } else {
          crossFn = markcrosscorr(pat_hs)
          for(i in c(1:length(crossFn$fns))) {
            crossFn$fns[[i]]$obs = crossFn$fns[[i]]$iso
          }
        }
        plot(crossFn)
  
        # Object
        outObject = crossFn
        outObject$label = basename(outFolder)
        saveRDS(outObject, file = file.path(outFolder,"crosscorr.rds"))
  
        # Statistics
        if(pcfCrossFlagDouble_hs == 1) {
          statsDf = getCurveStats(outFolder, crossFn)
        } else {
          varDf = melt(crossFn$which)
          varDf = varDf[order(varDf$value),]
          statsDf = cbind(varDf, do.call(rbind, lapply(crossFn$fns, function (x) getCurveStats(outFolder, x))))
        }
  
        # To get the sample size, different for quantitative vs. nominal marks
        markVars = marks(pat_hs)
        if (is.data.frame(markVars)) {
          statsDf$n = nrow(markVars)
        } else {
          statsDf$n = length(markVars[markVars != "Other"])
        }
  
        # Output statistics
        write.csv(statsDf, file.path(outFolder, "stats.csv"), quote = FALSE, row.names = FALSE)
  
        # Curve
        if(pcfCrossFlagDouble_hs == 1) {
          curveDf = data.frame(x = crossFn$r, y = crossFn$obs, label = basename(outFolder))
        } else {
          curveDf = adply(melt(crossFn$which), 1, function(xs) {
            data.frame(x = crossFn$fns[[xs$value]]$r, y = crossFn$fns[[xs$value]]$obs, label = basename(outFolder))
          })
        }
        write.csv(curveDf, file.path(outFolder, "curve.csv"), quote = FALSE, row.names = FALSE)
  
        return ()
      }
  
      mainRelationships = function(isLabelMarks, outFolder) {
        message(paste("Processing", outFolder))
        dir.create(outFolder, showWarnings = FALSE, recursive = TRUE)
  
        message("Plotting point process")
        pdf(file = file.path(outFolder, "basic_plot.pdf"))
        tryCatch(plot(pat_hs), error = function (e) message(paste("Could not plot basic_plot.pdf:", e)))
        dev.off()
  
        message("Plotting mark correlation function")
        # Mark correlation function.
        pdf(file = file.path(outFolder, "mark_correlation_function.pdf"))
        tryCatch(plot(markcorr(pat_hs)), error = function (e) message(paste("Could not plot mark_correlation_function.pdf:", e)))
        dev.off()
  
        # Mark variogram (lower value, more similar mark values at a distance).
        if (!isLabelMarks) {
          message("Plotting variogram")
          pdf(file = file.path(outFolder, "mark_variogram.pdf"))
          tryCatch(plot(markvario(pat_hs)), error = function (e) message(paste("Could not plot mark_variogram:", e)))
          dev.off()
        }
  
        # Envelope.
        message("Plotting envelope")
        pdf(file = file.path(outFolder, "envelope.pdf"))
        tryCatch(plot(envelope(pat_hs)), error = function (e) message(paste("Could not plot envelope", e)))
        dev.off()
  
        # Mark cross correlation
        message("Plotting mark cross correlation")
        pdf(file = file.path(outFolder, "cross_correlation_function.pdf"))
        tryCatch(crossCorrCurveStats(outFolder), error = function (e) message(paste("Could not compute mark cross correlation statistics:", e)))
        dev.off()
  
        return()
  
      }
  
      mainRelationships(isLabelMarks_hs, outDir_hs)
  
    |]
  
    pure ()
  