function parallelTemper(M = 2000, T = [0.1, 0.2, 0.3, 0.4]; d = 100, alpha0 = T[1])
x = 2*rand(Bernoulli(1/2), d, I) .- 1
if rand() < alpha0 # parallel step
# println(sum(abs.(y2.-y1)))
# idx = sample(1:I, 2, replace = false) # not neigbor
# rho = pdfIsing(x[:,idx[2]], T[idx[1]]) * pdfIsing(x[:,idx[1]], T[idx[2]]) / (pdfIsing(x[:,idx[1]], T[idx[1]]) * pdfIsing(x[:,idx[2]], T[idx[2]]))
rho = logpdfIsing(x[:,idx[2]], T[idx[1]]) + logpdfIsing(x[:,idx[1]], T[idx[2]]) - (logpdfIsing(x[:,idx[1]], T[idx[1]]) + logpdfIsing(x[:,idx[2]], T[idx[2]]))
x[:,idx[1]] .= x[:, idx[2]]
res[m, :] = sum(x, dims=1)