Итак, наконец-то я решил это (не лучшее решение, но оно все равно может кому-то помочь). Я решил начать с числового дифференцирования (центральных различий) и не использовать локальную сплайн / полиномиальную подгонку, как предлагалось @ 42- .
# our dataset consists of points (x, y, z), where z = f(x, y)
z_dx <- function(x, y) { # first partial derivative dz/dx
pos_x <- which(df_x == x)
if (pos_x == 1 | pos_x == length(df_x)) return(NA)
xf <- df_x[pos_x + 1]
xb <- df_x[pos_x - 1]
return((df[which(df$x == xf & df$y == y), "z"] - df[which(df$x == xb & df$y == y), "z"]) / (xf - xb))
}
z_dy <- function(x, y) { # first partial derivative dz/dy
pos_y <- which(df_y == y)
if (pos_y == 1 | pos_y == length(df_y)) return(NA)
yf <- df_y[pos_y + 1]
yb <- df_y[pos_y - 1]
return((df[which(df$y == yf & df$x == x), "z"] - df[which(df$y == yb & df$x == x), "z"]) / (yf - yb))
}
z_dx2 <- function(x, y) { # second partial derivative d2z/dx2
pos_x <- which(df_x == x)
if (pos_x <= 2 | pos_x >= (length(df_x) - 1)) return(NA)
xf <- df_x[pos_x + 1]
xb <- df_x[pos_x - 1]
return((df[which(df$x == xf & df$y == y), "dx"] - df[which(df$x == xb & df$y == y), "dx"]) / (xf - xb))
}
z_dy2 <- function(x, y) { # second partial derivative d2z/dy2
pos_y <- which(df_y == y)
if (pos_y <= 2 | pos_y >= (length(df_y) - 1)) return(NA)
yf <- df_y[pos_y + 1]
yb <- df_y[pos_y - 1]
return((df[which(df$y == yf & df$x == x), "dy"] - df[which(df$y == yb & df$x == x), "dy"]) / (yf - yb))
}
z_dxdy <- function(x, y) { # second partial derivative d2z/dxdy
pos_y <- which(df_y == y)
if (pos_y <= 2 | pos_y >= (length(df_y) - 1)) return(NA)
yf <- df_y[pos_y + 1]
yb <- df_y[pos_y - 1]
return((df[which(df$y == yf & df$x == x), "dx"] - df[which(df$y == yb & df$x == x), "dx"]) / (yf - yb))
}
# read dataframe and set proper column names
df <- read.table("map.dat", header = T)
names(df) <- c("x", "y", "z")
# extract unique abscissae
df_x <- unique(df$x)
df_y <- unique(df$y)
# calculate first derivatives numerically
df$dx <- apply(df, 1, function(x) z_dx(x["x"], x["y"]))
df$dy <- apply(df, 1, function(x) z_dy(x["x"], x["y"]))
# set tolerance limits, so we can filter out every non-stationary point later
tolerance_x <- 30.0
tolerance_y <- 10.0
# select only stationary points
df$stat <- apply(df, 1, function(x) ifelse(is.na(x["dx"]) | is.na(x["dy"]), FALSE, ifelse(abs(x["dx"]) <= tolerance_x & abs(x["dy"]) <= tolerance_y, TRUE, FALSE)))
stationaries <- df[which(df$stat == TRUE), ]
# compute second derivatives for selected points
stationaries$dx2 <- apply(stationaries, 1, function(x) z_dx2(x["x"], x["y"]))
stationaries$dy2 <- apply(stationaries, 1, function(x) z_dy2(x["x"], x["y"]))
stationaries$dxdy <- apply(stationaries, 1, function(x) z_dxdy(x["x"], x["y"]))
# let's make classifying function...
stat_type <- function(dx2, dy2, dxdy) {
if (is.na(dx2) | is.na(dy2) | is.na(dxdy)) return("unk")
if (dx2 * dy2 - dxdy^2 < 0) return("saddle")
if (dx2 * dy2 - dxdy^2 > 0) {
if (dx2 < 0 & dy2 < 0) return("max")
else if (dx2 > 0 & dy2 > 0) return("min")
else return("unk")
}
else return("unk")
}
# ...and apply it to our stationary points
stationaries$type <- apply(stationaries, 1, function(x) stat_type(x["dx2"], x["dy2"], x["dxdy"]))
stationaries <- stationaries[which(stationaries$type != "unk"), ] # leave out all non-stationarities
# set upper cutoff (in the units of z = f(x, y)) - points with the greatest value of z can form very large areas where derivatives are so close to zero that we can't discriminate them (this is often the case when our map is built by numerical integration with predefined grid so final result is discrete function of two variables)
cutwidth <- 2.5
stationaries <- stationaries[which(stationaries$z <= (max(stationaries$z) - cutwidth)), ] # select only points which are lying below maximum by selected cutoff
# create dataframes for minima, maxima and saddles
minima <- stationaries[which(stationaries$type == "min"), ]
maxima <- stationaries[which(stationaries$type == "max"), ]
saddles <- stationaries[which(stationaries$type == "saddle"), ]
Все становится проще, если ваша карта имеет постоянный размер сетки - тогда нет необходимости вычислять размер шага для дифференциации для каждой точки