## Time-stamp: "Last modified 2024-04-11 15:28:06 delucia"

## Compile the fortran code. On Linux, if a fortran compiler is
## installed, it should work out of the box

setwd("/home/work/simR/Rphree/SpecSim/")
system("R CMD SHLIB SS_Sim.f SS_Vario.f -o SSlib.so")
dyn.load("SSlib.so")

## Load the functions
source("SS_Fun_sim.R")
source("SS_Fun_vario.R")
setwd("/home/work/simR/Rphree/NAAICE_tug_input/build")

library(lattice)
library(PoetUtils)

PlotField2 <- function (data, grid, nx, ny, contour = TRUE, nlevels = 12, breaks, 
    palette = "heat.colors", rev.palette = TRUE, scale = TRUE, 
    plot.axes = TRUE, ...) 
{
    if (!missing(grid)) {
        xc <- unique(sort(grid$cell$XCOORD))
        yc <- unique(sort(grid$cell$YCOORD))
        nx <- length(xc)
        ny <- length(yc)
        if (!length(data) == nx * ny) 
            stop("Wrong nx, ny or grid")
    }
    else {
        xc <- seq(1, nx)
        yc <- seq(1, ny)
    }
    z <- matrix(round(data, 6), ncol = nx, nrow = ny, byrow = TRUE)
    pp <- t(z[rev(seq(1, nrow(z))), ])
    if (missing(breaks)) {
        breaks <- pretty(data, n = nlevels)
    }
    breakslen <- length(breaks)
    colors <- do.call(palette, list(n = breakslen - 1))
    if (rev.palette) 
        colors <- rev(colors)
    if (scale) {
        par(mfrow = c(1, 2))
        nf <- layout(matrix(c(1, 2), 1, 2, byrow = TRUE), widths = c(4, 
            1))
    }
    par(las = 1, mar = c(5, 5, 3, 1))
    image(xc, yc, pp, las = 1, asp = 1, breaks = breaks, col = colors,
          axes = FALSE, ann = plot.axes, ...)
    if (plot.axes) {
        axis(1)
        axis(2)
    }
    if (contour) 
        contour(unique(sort(xc)), unique(sort(yc)), pp, breaks = breaks, 
            add = TRUE)
    if (scale) {
        par(las = 1, mar = c(5, 1, 5, 5))
        PlotImageScale(data, breaks = breaks, add.axis = FALSE, 
            axis.pos = 4, col = colors)
        axis(4, at = breaks)
    }
    invisible(pp)
}


## generate a grid. At the moment it must be a cartesian square grid
## with number of cells a power of two, e.g. 128, 256, ...
n <- 256            ## we want a 256x256 grid

## Simple structure describing the simulation grid
gr <- SS.Grid(nodes=n, side=100) 

## fix the seed of RNG for reproducibility
set.seed(64276)

## Extract the random numbers (uniform distribution between 0 and 1)
## once and for all
randoms <- runif(n**2)

## Anisotropic SpectralSimulations with 2 different isotropic correlation lengths
## using the same random numbers. Here "range" is in meters!
a1 <- SS.Sim(vec.rn = randoms, nugget=0., sill=1, type="spherical", range=20, ratio=8, angle = -30, grid=gr, normalize=TRUE)

## Compute the experimental variograms ONLY in X and Y direction.
## "lag" refers to the nodes on one side of the grid, not meters!
variog1 <- SS.Vario(a1, grid=gr, lag=2, nstep=72)

## Visualize

# cairo_pdf("Spherical_40_80.pdf", height=9, width=9)
par(mfrow=c(1,2))
SS.Plot(a1, grid=gr, main="Range 20 m, anisotropy ratio 4", axes=FALSE)
axis(1, at=c(-100, 0, 100))
axis(2, at=c(-100, 0, 100))
SS.VarioPlot(vario=variog1, main="Experimental Variograms", lwd=2) 
lines(variog1$htrue, SS.Spherical(variog1$htrue, range=20, sill=1), col="blue") 
lines(variog1$htrue, SS.Spherical(variog1$htrue, range=20/4, sill=1), col="red")


sim1 <- matrix(a1, n, n)
levelplot(sim1)

cut1 <- sim1[1:200, 1:200]

levelplot(cut1)
sd(cut1)
mean(cut1)

ax <- 1E-6*(cut1 - min(cut1))/max(cut1)+1E-7
ay <- 1E-7*(cut1 - min(cut1))/max(cut1)+1E-7

sd(ax)
mean(ax)

sd(ay)
mean(ay)

range(ax)
range(ay)



PlotField2(log10(ax), nx=200, ny=200, contour = FALSE, xaxs="i", yaxs="i", nlevels=6, plot.axes = FALSE)

title(expression(log[10](alpha[x])))

levelplot(log10(t(ax)))

levelplot(log10(ay))

data.table::fwrite(ax, "alpha_x.csv", col.names = FALSE)
data.table::fwrite(ay, "alpha_y.csv", col.names = FALSE)

init <- c(110.01240000000783,55.50867875567298,-1.2162968960604385e-9,4.455329999159673e-7,2.0000000000005577e-12,0.0006151616266928334,0.0006147160946983908)

names <- c("H","O","Charge","Ba","Cl","S_6_","Sr")

inmat <- matrix(rep(init, 200*200), nrow = 200*200, byrow = TRUE)

colnames(inmat) <- names

data.table::fwrite(inmat, "barite_200_init.csv", col.names = TRUE)

library(PoetUtils)

out <- fread("../../build/barite_200_output.csv")

cairo_pdf("barite_200_field_alphax.pdf", width = 9, height = 6, family="serif")
PlotField2(log10(ax), nx=200, ny=200, contour = FALSE, nlevels=10,
           palette = "heat.colors", rev.palette = FALSE,
           plot.axes = FALSE, main=expression(log[10](Sr)))
dev.off()

cairo_pdf("barite_200_field_Ba.pdf", width = 9, height = 6, family="serif")
PlotField2(log10(out$Ba), nx=200, ny=200, contour = FALSE, nlevels=10,
           palette = "terrain.colors", plot.axes = FALSE, main=expression(log[10](Sr)))
dev.off()

PlotField2(log10(out$H), nx=200, ny=200, contour = FALSE, palette = "terrain.colors", nlevels=10, plot.axes = FALSE, main=expression(log[10](Sr)))

PlotField2(out$S_6_, nx=200, ny=200, contour = FALSE, nlevels=10, plot.axes = FALSE, main=expression(log[10](Ba)))

###################### barite large

alphamat <- matrix(1e-6, nrow = 1000, ncol = 1000)
data.table::fwrite(alphamat, "barite_large_alpha.csv", col.names = FALSE)


init <- c(110.0124, 55.508678, -1.216296896e-9, 4.45533e-7, 0, 0.0006151616, 0.0006147161)
nams <- c("H","O","Charge","Ba","Cl","S_6_","Sr")
names(init) <- nams

n <- 1000
inmat <- matrix(rep(init, 1000*1000-n), nrow = 1000*1000 - n, byrow = TRUE)
colnames(inmat) <- nams

bounds <- c(111.0124, 55.50622, -3.0E-07, 1, 2, 0.01, 0.001)
names(bounds) <- nams
bmat <- matrix(rep(bounds, n), nrow = n, byrow = TRUE)
colnames(bmat) <- nams

set.seed(9342)
inds <- sample(seq_len(1E6), size=1E6)

dim(inmat)
dim(bmat)

totmat <- rbind(bmat, inmat)
dim(totmat)

fmat <- totmat[inds, ]


xc <- rep(seq(1, 1000), times = 1000)
yc <- rep(seq(1, 1000), each = 1000)
iplot <- which(fmat[,"Cl"]==2)

plot(xc[iplot], 1000-yc[iplot], pch=4, cex=0.5, xaxs="i",yaxs="i")


cairo_pdf("barite_large_init_locs.pdf", width = 6, height = 6, family="serif")
par(mar=c(1,1,0.5,0.5))
plot(xc[iplot], yc[iplot], pch=4, cex=0.5, las=1, xaxs="i", yaxs="i",
     xlab="Easting", ylab="Northig", asp=1, axes = FALSE)
box()
dev.off()

data.table::fwrite(fmat, "barite_large_init.csv", col.names = TRUE)


out <- fread("../../build/barite_large_output.csv")

x11()

cairo_pdf("barite_large_field_Ba.pdf", width = 8, height = 6, family="serif")
##par(mar=c(4,4,0.5,0.5))
PlotField2(log10(out$Ba), nx=1000, ny=1000, contour = FALSE, nlevels=10,
           palette = "terrain.colors", plot.axes = FALSE)
box()
dev.off()


init

round(rbind(init, bounds), 3)



##### surfex, grid 500x250

nx <- 200
ny <- 100

alphamat <- matrix(1.1e-12, nrow = ny, ncol = nx)
data.table::fwrite(alphamat, "surfex_alpha.csv", col.names = FALSE)

a <- read.table(
    textConnection(
"H             1.11e+02 120.0,    
O             5.55e+01  55.1,       
Charge       -2.0e-13   8.0e-17,     
C(-4)         2.0e-16   2.0e-15,     
C(4)          2.0e-03   0.2,       
Ca            2.0e-01   0.03,        
Cl            3.0e-01   0.5,         
Fe(2)         1.4e-04   0.0002,      
Fe(3)         1.3e-09   2.0e-08,     
H(0)          6.0e-12   2.0e-11,     
K             2.0e-03   1.0e-05,     
Mg            1.0e-02   0.2,         
Na            2.0e-01   0.3,         
S(-2)         5.9e-10   0,           
S(2)          8.3e-15   8.3e-12,     
S(4)          2.1e-14   5.1e-14,     
S(6)          1.6e-02   0.026,       
Sr            4.5e-04   0.045,       
U(4)          2.5e-09   2.5e-08,     
U(5)          1.6e-10   1.6e-10,     
U(6)          2.3e-07   1.0e-05" 
))
vals <- t(a[, 2:3])

nams <- c("H", "O", "Charge", "CH4", "C", "Ca", "Cl", "Fe2", "Fe3",
          "H0", "K", "Mg", "Na", "HS2", "S2", "S4", "S6","Sr", "U4",
          "U5", "U6")
colnames(vals) <- nams
rownames(vals) <- c("IC", "BC")
vals

length(nams)

dim(vals)
t(vals)

ic <- matrix(rep(vals[1,], nx*ny), nrow = nx*ny, byrow = TRUE)

colnames(ic) <- nams

data.table::fwrite(ic, "surfex_init.csv", col.names = TRUE)


dput(vals[2,])

out <- fread("../../build/surfex_output.csv")

x11()

cairo_pdf("surfex_field_U6.pdf", width = 10, height = 6, family="serif")
out <- fread("../../build/surfex_output.csv")
PlotField2(log10(out$U6), nx=nx, ny=ny, contour = FALSE, nlevels=10,
           palette = "cm.colors", plot.axes = FALSE, rev.palette = FALSE)
dev.off()

cairo_pdf("surfex_field_Na.pdf", width = 10, height = 6, family="serif")
out <- fread("../../build/surfex_output.csv")
PlotField2(out$Na, nx=nx, ny=ny, contour = FALSE, nlevels=10,
           palette = "topo.colors", plot.axes = FALSE, rev.palette = TRUE)
dev.off()



PlotField2(log10(out$U4), nx=nx, ny=ny, contour = FALSE, nlevels=10,
           palette = "terrain.colors", plot.axes = FALSE)
out$U4[seq(50*200, 51*200)]




PlotField2(out$Na, nx=nx, ny=ny, contour = FALSE, nlevels=10,
           palette = "terrain.colors", plot.axes = FALSE)
out$Na[seq(50*200, 51*200)]
out$Cl[seq(50*200, 51*200)]


PlotField2(log10(out$Fe2), nx=nx, ny=ny, contour = FALSE, nlevels=12,
           palette = "terrain.colors", plot.axes = FALSE)

PlotField2((out$Cl), nx=nx, ny=ny, contour = FALSE, nlevels=10,
           palette = "terrain.colors", plot.axes = FALSE)


dev.off()

##### debug, grid 200x100

nx <- 200
ny <- 100

alphamat <- matrix(1.1e-11, nrow = ny, ncol = nx)
data.table::fwrite(alphamat, "debug_alpha.csv", col.names = FALSE)

a <- read.table(
    textConnection(
"H             1.11e+02 120.0,    
O             5.55e+01  55.1,       
Charge       -2.0e-13   8.0e-17,     
C(-4)         2.0e-16   2.0e-15,     
C(4)          2.0e-03   0.2,       
Ca            2.0e-01   0.03,        
Cl            3.0e-01   0.5,         
Fe(2)         1.4e-04   0.0002,      
Fe(3)         1.3e-09   2.0e-08,     
H(0)          6.0e-12   2.0e-11,     
K             2.0e-03   1.0e-05,     
Mg            1.0e-02   0.2,         
Na            2.0e-01   0.3,         
S(-2)         5.9e-10   0,           
S(2)          8.3e-15   8.3e-12,     
S(4)          2.1e-14   5.1e-14,     
S(6)          1.6e-02   0.026,       
Sr            4.5e-04   0.045,       
U(4)          2.5e-09   2.5e-08,     
U(5)          1.6e-10   1.6e-10,     
U(6)          2.3e-07   1.0e-05" 
))
vals <- t(a[, 2:3])

nams <- c("H", "O", "Charge", "CH4", "C", "Ca", "Cl", "Fe2", "Fe3",
          "H0", "K", "Mg", "Na", "HS2", "S2", "S4", "S6","Sr", "U4",
          "U5", "U6")
colnames(vals) <- nams
rownames(vals) <- c("IC", "BC")
vals

length(nams)

dim(vals)
t(vals)

ic <- matrix(rep(vals[1,], nx*ny), nrow = nx*ny, byrow = TRUE)

colnames(ic) <- nams

data.table::fwrite(ic, "debug_init.csv", col.names = TRUE)


dput(vals[2,])

out <- fread("debug_output.csv")

x11()

## cairo_pdf("debug_field_U6.pdf", width = 10, height = 6, family="serif")
# out <- fread("./debug_output.csv")
PlotField2(log10(out$U6), nx=nx, ny=ny, contour = FALSE, nlevels=10,
           palette = "cm.colors", plot.axes = FALSE, rev.palette = FALSE)
##dev.off()

## cairo_pdf("debug_field_Na.pdf", width = 10, height = 6, family="serif")
out2 <- fread("./debug_output.csv")
PlotField2(out2$Na, nx=nx, ny=ny, contour = FALSE, nlevels=10,
           palette = "topo.colors", plot.axes = FALSE, rev.palette = TRUE)
## dev.off()
out2$Na[seq(1, 200)]

range(out2$Na)


PlotField2(log10(out$U4), nx=nx, ny=ny, contour = FALSE, nlevels=10,
           palette = "terrain.colors", plot.axes = FALSE)
out$U4[seq(50*200, 51*200)]




PlotField2(out2$Fe3, nx=nx, ny=ny, contour = FALSE, nlevels=10,
           palette = "terrain.colors", plot.axes = FALSE)

out2$Na[seq(20*200, 21*200)]

out$Cl[seq(50*200, 51*200)]


PlotField2(log10(out$Fe2), nx=nx, ny=ny, contour = FALSE, nlevels=12,
           palette = "terrain.colors", plot.axes = FALSE)

PlotField2((out$Cl), nx=nx, ny=ny, contour = FALSE, nlevels=10,
           palette = "terrain.colors", plot.axes = FALSE)


dev.off()

library(Rcpp)

cppFunction("static const std::vector<int> gen_seq(int from, int to) {
  int vsize = to - from + 1;
  std::vector<int> vec(vsize);
  for (int i =  0; i < vsize; i++) {
    vec[i] = i+from;
  }
  return vec;
}")

gen_seq(24, 49)

cppFunction("static const std::vector<int> gen_vec(int elems) {
  std::vector<int> vec(elems);
  for (int i = 0; i < elems; i++) {
    vec[i] = i;
  }
  return vec;
}")


gen_vec(20)