emd_sample.r

library(lpSolve)

# 欧几里得距离
euclid_dist <- function(f1, f2) {
    return(sqrt(sum((f1 - f2)^2)))
}

# 计算EMD
emd <- function(dist, w1, w2) {
    # 准备lp.transport()参数
    costs <- dist
    row.signs <- rep("<", length(w1))
    row.rhs <- w1
    col.signs <- rep(">", length(w2))
    col.rhs <- w2
    # 解决运输问题
    t <- lp.transport(costs, "min", row.signs, row.rhs, col.signs, col.rhs)
    # 获得最优解
    flow <- t$solution
    # 计算最少工作量
    work <- sum(flow * dist)
    # 对EMD结果作归一化
    e <- work / sum(flow)
    return(e)
}

# 特征量
f1 = matrix(c(100, 40, 22, 211, 20, 2, 32, 190, 150, 2, 100, 100), 4, 3, byrow=T)
f2 = matrix(c(0, 0, 0, 50, 100, 80, 255, 255, 255), 3, 3, byrow=T)

# 权重（要用整数）
w1 = c(4, 3, 2, 1)
w2 = c(5, 3, 2)

n1 = length(f1[,1])
n2 = length(f2[,1])

# 创建一个距离矩阵
dist = matrix(0, n1, n2)
for (i in 1:n1) {
    for (j in 1:n2) {
        dist[i, j] = euclid_dist(f1[i,], f2[j,])
    }
}

# 从权重和距离矩阵得到EMD
e = emd(dist, w1, w2)
cat(sprintf("emd = %f\n", e))