259 lines
5.5 KiB
Julia
259 lines
5.5 KiB
Julia
module OptimMixture
|
|
|
|
import DelimitedFiles as DF
|
|
import BenchmarkTools as BT
|
|
import UnicodePlots as UP
|
|
import LinearAlgebra: ⋅
|
|
|
|
const ENERGIES = [2803, 2811, 2819, 2826, 2834, 2842, 2850, 2858, 2866, 2874,
|
|
2882, 2890, 2897, 2905, 2913, 2921, 2929, 2937, 2945, 2953,
|
|
2961, 2969, 2977, 2985, 2993, 3001, 3009, 3018, 3026, 3034,
|
|
3042, 3050]
|
|
|
|
@kwdef struct Slab
|
|
energy::Int
|
|
data::Matrix{Int}
|
|
end
|
|
|
|
function load_slabs(width::Int=512, height::Int=512)
|
|
slabs = Slab[]()
|
|
i_min, i_max = typemax(Int), 0
|
|
for e in ENERGIES
|
|
slab = Slab(energy=e, data=DF.readdlm("test_sample/HeLa_F-SRS_512x512_2803cm-1.txt", ',', Int))
|
|
push!(slabs, slab)
|
|
|
|
i_min, i_max = min(i_min, minimum(slab.data)), max(i_max, maximum(slab.data))
|
|
|
|
return slabs, i_min, i_max
|
|
end
|
|
end
|
|
|
|
function Q_matrix_transpose(X::Matrix{Float64}, Y::Matrix{Float64}, N::Int)
|
|
return sum(abs2, Y - view(X, 1:N, :) * view(X, (N+1):size(X, 1), :)')
|
|
end
|
|
|
|
function Q_matrix_reshape(X::Matrix{Float64}, Y::Matrix{Float64}, N::Int)
|
|
N_plus_M, m = size(X)
|
|
return sum(abs2, Y - view(X, 1:N, :) * reshape(view(X, (N+1):N_plus_M, :), m, N_plus_M - N))
|
|
end
|
|
|
|
function E_loop(X::Matrix{Float64}, Y::Matrix{Float64})
|
|
s = 0
|
|
N, M = size(Y)
|
|
m = size(X, 2)
|
|
for i in 1:N
|
|
for j in N+1:N+M
|
|
x = 0
|
|
for k in 1:m
|
|
@inbounds x += X[i, k] * X[j, k]
|
|
end
|
|
@inbounds s += (x - Y[i, j-N])^2
|
|
end
|
|
end
|
|
return s
|
|
end
|
|
|
|
function Q_loop!(dst::Matrix{Float64}, X::Matrix{Float64}, N::Int)
|
|
N_plus_M, m = size(X)
|
|
for i in 1:N
|
|
for j in N+1:N_plus_M
|
|
x = 0
|
|
for k in 1:m
|
|
@inbounds x += X[i, k] * X[j, k]
|
|
end
|
|
@inbounds dst[i, j-N] = x
|
|
end
|
|
end
|
|
end
|
|
|
|
function DE_loop!(dst::Matrix{Float64}, X::Matrix{Float64}, Y::Matrix{Float64})
|
|
dst_W = zeros(size(Y))
|
|
DE_loop!(dst, dst_W, X, Y)
|
|
end
|
|
|
|
function dot_prod_dual_loop(W::Matrix{Float64}, X::Matrix{Float64}, V::Matrix{Float64})
|
|
_, m = size(X)
|
|
N, M = size(W)
|
|
|
|
s = 0
|
|
|
|
for k in 1:m
|
|
for i in 1:N
|
|
x = 0
|
|
for j in 1:M
|
|
@inbounds x += W[i, j] * X[N+j, k]
|
|
end
|
|
s += x * V[i, k]
|
|
end
|
|
for j in 1:M
|
|
x = 0
|
|
for i in 1:N
|
|
@inbounds x += W[i, j] * X[i, k]
|
|
end
|
|
s += x * V[N+j, k]
|
|
end
|
|
end
|
|
|
|
return s
|
|
end
|
|
|
|
function dot_prod_primal_loop(W::Matrix{Float64}, X::Matrix{Float64}, V::Matrix{Float64})
|
|
_, m = size(X)
|
|
N, M = size(W)
|
|
|
|
s = 0
|
|
|
|
for i in 1:N
|
|
for j in 1:M
|
|
x = 0
|
|
for k in 1:m
|
|
@inbounds x += X[i, k] * V[N+j, k] + X[N+j, k] * V[i, k]
|
|
|
|
end
|
|
s += x * W[i, j]
|
|
end
|
|
end
|
|
|
|
return s
|
|
end
|
|
|
|
function DE_loop!(dst::Matrix{Float64}, dst_W::Matrix{Float64}, X::Matrix{Float64}, Y::Matrix{Float64})
|
|
_, m = size(X)
|
|
N, M = size(Y)
|
|
# dst is the same size as X
|
|
# we need storage of size Y
|
|
Q_loop!(dst_W, X, N)
|
|
# compute W = Q(X) - Y
|
|
dst_W .-= Y
|
|
for k in 1:m
|
|
for i in 1:N
|
|
x = 0
|
|
for j in 1:M
|
|
@inbounds x += dst_W[i, j] * X[N+j, k]
|
|
end
|
|
@inbounds dst[i, k] = x
|
|
end
|
|
for j in 1:M
|
|
x = 0
|
|
for i in 1:N
|
|
@inbounds x += dst_W[i, j] * X[i, k]
|
|
end
|
|
@inbounds dst[N+j, k] = x
|
|
end
|
|
end
|
|
dst .*= 2
|
|
end
|
|
|
|
function DQ_V_loop!(dst::Matrix{Float64}, X::Matrix{Float64}, V::Matrix{Float64}, N::Int)
|
|
N_plus_M, m = size(X)
|
|
for i in 1:N
|
|
for j in N+1:N_plus_M
|
|
x = 0
|
|
for k in 1:m
|
|
@inbounds x += (V[i, k] * X[j, k] + X[i, k] * V[j, k])
|
|
end
|
|
@inbounds dst[i, j-N] = x
|
|
end
|
|
end
|
|
end
|
|
|
|
function Q_loop_transposed(X::Matrix{Float64}, Y::Matrix{Float64}, N::Int)
|
|
s = 0
|
|
m, N_plus_M = size(X)
|
|
for i in 1:N
|
|
for j in N+1:N_plus_M
|
|
x = 0
|
|
for k in 1:m
|
|
@inbounds x += X[k, i] * X[k, j]
|
|
end
|
|
@inbounds s += (Y[j-N, i] - x)^2
|
|
end
|
|
end
|
|
return s
|
|
end
|
|
|
|
function test_dot_prod(N::Int, M::Int, m::Int)
|
|
X = rand(N + M, m)
|
|
W = rand(N, M)
|
|
V = rand(N + M, m) .- 0.5
|
|
|
|
@show dot_prod_primal_loop(W, X, V)
|
|
@show dot_prod_dual_loop(W, X, V)
|
|
end
|
|
|
|
function test_DE(N::Int, M::Int, m::Int)
|
|
X = rand(N + M, m)
|
|
Y = rand(N, M)
|
|
V = rand(N + M, m) .- 0.5
|
|
|
|
EX = E_loop(X, Y)
|
|
grad_EX = zeros(size(X))
|
|
DE_loop!(grad_EX, X, Y)
|
|
|
|
eps = [10.0^k for k in 1:-1:-12]
|
|
errors = zeros(size(eps))
|
|
|
|
for i in eachindex(eps)
|
|
EX_V_true = E_loop(X + eps[i] * V, Y)
|
|
EX_V_approx = EX + eps[i] * sum(grad_EX .* V)
|
|
@show EX_V_true
|
|
@show EX_V_approx
|
|
errors[i] = abs(EX_V_true - EX_V_approx)
|
|
end
|
|
|
|
@show eps
|
|
@show errors
|
|
|
|
foo = UP.lineplot(eps, errors, xscale=:log10, yscale=:log10)
|
|
UP.lineplot!(foo, eps, eps)
|
|
println(foo)
|
|
end
|
|
|
|
function test_DQ(N::Int, M::Int, m::Int)
|
|
X = rand(N + M, m) * 255
|
|
V = rand(N + M, m)
|
|
|
|
QX = zeros(N, M)
|
|
Q_loop!(QX, X, N)
|
|
|
|
dst_true = zeros(N, M)
|
|
dst_approx = zeros(N, M)
|
|
|
|
eps = [10.0^k for k in 1:-1:-8]
|
|
errors = zeros(size(eps))
|
|
|
|
for i in eachindex(eps)
|
|
Q_loop!(dst_true, X + eps[i] * V, N)
|
|
DQ_V_loop!(dst_approx, X, eps[i] * V, N)
|
|
dst_approx .+= QX
|
|
errors[i] = sum(abs, dst_approx - dst_true)
|
|
end
|
|
|
|
@show eps
|
|
@show errors
|
|
|
|
println(UP.lineplot(eps, errors, xscale=:log10, yscale=:log10))
|
|
end
|
|
|
|
function test(N::Int, M::Int, m::Int)
|
|
Y = rand(N, M) * 255
|
|
X = rand(N + M, m) * 255
|
|
|
|
X_pre = copy(X)
|
|
X_pre[N+1:N+M] .= vec(X_pre[N+1:N+M]')
|
|
|
|
r_matrix = Q_matrix_transpose(X, Y, N)
|
|
r_loop = Q_loop(X, Y, N)
|
|
|
|
@show r_matrix, r_loop, abs(r_matrix - r_loop) / r_loop
|
|
end
|
|
|
|
end # module
|
|
|
|
# OptimMixture.test(512 * 512, 20, 2)
|
|
# OptimMixture.test_DQ(512 * 512, 20, 2)
|
|
OptimMixture.test_DE(512 * 512, 20, 2)
|
|
# OptimMixture.test_DE(20, 2, 2)
|
|
# OptimMixture.test_dot_prod(512^2, 20, 2)
|
|
|
|
# vim: ts=2:sw=2:sts=2
|