Skip to content

Commit

Permalink
Fix: temperature was not correctly applied.
Browse files Browse the repository at this point in the history
  • Loading branch information
@shinmorino
shinmorino committed Nov 7, 2018
1 parent 70b48d0 commit 25de2d8
Show file tree
Hide file tree
Showing 8 changed files with 50 additions and 47 deletions.
24 changes: 12 additions & 12 deletions sqaodc/cpu/CPUBipartiteGraphAnnealer.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -406,10 +406,10 @@ void CPUBipartiteGraphAnnealer<real>::annealOneStepColoringParallel(real G, real


template<class real>
void CPUBipartiteGraphAnnealer<real>::annealOneStepSANaive(real kT, real beta) {
void CPUBipartiteGraphAnnealer<real>::annealOneStepSANaive(real kT, real _) {
throwErrorIfQNotSet();

real Tnorm = kT * beta;
real invKT = real(1.) / kT;

sq::Random &random = random_[0];
int N = N0_ + N1_;
Expand All @@ -420,7 +420,7 @@ void CPUBipartiteGraphAnnealer<real>::annealOneStepSANaive(real kT, real beta) {
real qyx = matQ0_(y, x);
real sum = J_.transpose().row(x).dot(matQ1_.row(y));
real dE = real(2.) * qyx * (h0_(x) + sum);
real threshold = (dE < real(0.)) ? real(1.) : std::exp(-dE * Tnorm);
real threshold = (dE < real(0.)) ? real(1.) : std::exp(-dE * invKT);
if (threshold > random.random<real>())
matQ0_(y, x) = - qyx;
}
Expand All @@ -429,7 +429,7 @@ void CPUBipartiteGraphAnnealer<real>::annealOneStepSANaive(real kT, real beta) {
real qyx = matQ1_(y, x);
real sum = J_.row(x).dot(matQ0_.row(y));
real dE = real(2.) * qyx * (h1_(x) + sum);
real threshold = (dE < real(0.)) ? real(1.) : std::exp(-dE * Tnorm);
real threshold = (dE < real(0.)) ? real(1.) : std::exp(-dE * invKT);
if (threshold > random.random<real>())
matQ1_(y, x) = - qyx;
}
Expand All @@ -444,11 +444,11 @@ void CPUBipartiteGraphAnnealer<real>::annealOneStepSANaive(real kT, real beta) {
template<class real, class T> static inline
void tryFlipSA(sq::EigenMatrixType<real> &qAnneal, int im,
const sq::EigenMatrixType<real> &dEmat, const sq::EigenRowVectorType<real> &h,
const T &J, sq::SizeType N, real Tnorm, sq::Random &random) {
const T &J, sq::SizeType N, real invKT, sq::Random &random) {
for (int iq = 0; iq < N; ++iq) {
real q = qAnneal(im, iq);
real dE = real(2.) * q * (h[iq] + dEmat(im, iq));
real thresh = dE < real(0.) ? real(1.) : std::exp(- dE * Tnorm);
real thresh = dE < real(0.) ? real(1.) : std::exp(- dE * invKT);
if (thresh > random.random<real>())
qAnneal(im, iq) = -q;
}
Expand All @@ -458,7 +458,7 @@ template<class real>
void CPUBipartiteGraphAnnealer<real>::
annealHalfStepSAColoring(int N, EigenMatrix &qAnneal,
const EigenRowVector &h, const EigenMatrix &J,
const EigenMatrix &qFixed, real Tnorm) {
const EigenMatrix &qFixed, real invKT) {
#ifdef _OPENMP
EigenMatrix dEmat(qFixed.rows(), J.rows());
// dEmat = qFixed * J.transpose(); // For debug
Expand All @@ -475,7 +475,7 @@ annealHalfStepSAColoring(int N, EigenMatrix &qAnneal,
sq::Random &random = random_[threadNum];
# pragma omp for
for (int im = 0; im < m_; ++im)
tryFlipSA(qAnneal, im, dEmat, h, J, N, Tnorm, random);
tryFlipSA(qAnneal, im, dEmat, h, J, N, invKT, random);
}
#else
abort_("Must not reach here.");
Expand All @@ -484,13 +484,13 @@ annealHalfStepSAColoring(int N, EigenMatrix &qAnneal,


template<class real>
void CPUBipartiteGraphAnnealer<real>::annealOneStepSAColoring(real kT, real beta) {
void CPUBipartiteGraphAnnealer<real>::annealOneStepSAColoring(real kT, real _) {
throwErrorIfQNotSet();
clearState(solSolutionAvailable);

real Tnorm = kT * beta;
annealHalfStepSAColoring(N1_, matQ1_, h1_, J_, matQ0_, Tnorm);
annealHalfStepSAColoring(N0_, matQ0_, h0_, J_.transpose(), matQ1_, Tnorm);
real invKT = real(1.) / kT;
annealHalfStepSAColoring(N1_, matQ1_, h1_, J_, matQ0_, invKT);
annealHalfStepSAColoring(N0_, matQ0_, h0_, J_.transpose(), matQ1_, invKT);
}


Expand Down
6 changes: 3 additions & 3 deletions sqaodc/cpu/CPUBipartiteGraphAnnealer.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -83,9 +83,9 @@ class CPUBipartiteGraphAnnealer : public sq::BipartiteGraphAnnealer<real> {
/* simulated annealing */
void annealHalfStepSAColoring(int N, EigenMatrix &qAnneal,
const EigenRowVector &h, const EigenMatrix &J,
const EigenMatrix &qFixed, real Tnorm);
void annealOneStepSANaive(real G, real beta);
void annealOneStepSAColoring(real G, real beta);
const EigenMatrix &qFixed, real invKT);
void annealOneStepSANaive(real kT, real _);
void annealOneStepSAColoring(real kT, real _);


void syncBits();
Expand Down
12 changes: 6 additions & 6 deletions sqaodc/cpu/CPUDenseGraphAnnealer.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -326,7 +326,7 @@ void CPUDenseGraphAnnealer<real>::annealOneStepColoringParallel(real G, real bet

template<class real> inline static
void tryFlipSA(sq::MatrixType<real> &matQ, int y, const sq::VectorType<real> &h, const sq::MatrixType<real> &J,
sq::Random &random, real Tnorm) {
sq::Random &random, real invKT) {
int N = J.rows;
int x = random.randInt(N);
real qyx = matQ(y, x);
Expand All @@ -337,17 +337,17 @@ void tryFlipSA(sq::MatrixType<real> &matQ, int y, const sq::VectorType<real> &h,
#else
real sum = dot_naive(J.rowPtr(x), matQ.rowPtr(y), N);
#endif
real dE = real(2.) * qyx * (h(x) + sum);
real threshold = (dE < real(0.)) ? real(1.) : std::exp(-dE * Tnorm);
real dE = real(2.) * qyx * (h(x) + 2. * sum);
real threshold = (dE < real(0.)) ? real(1.) : std::exp(-dE * invKT);
if (threshold > random.random<real>())
matQ(y, x) = - qyx;
}

template<class real>
void CPUDenseGraphAnnealer<real>::annealOneStepSANaive(real kT, real beta) {
void CPUDenseGraphAnnealer<real>::annealOneStepSANaive(real kT, real _) {
throwErrorIfQNotSet();

real Tnorm = kT * beta;
real invKT = real(1.) / kT;

#ifndef _OPENMP
sq::Random &random = random_[0];
Expand All @@ -357,7 +357,7 @@ void CPUDenseGraphAnnealer<real>::annealOneStepSANaive(real kT, real beta) {
#endif
for (int y = 0; y < m_; ++y) {
for (int loop = 0; loop < sq::IdxType(N_); ++loop)
tryFlipSA(matQ_, y, h_, J_, random, Tnorm);
tryFlipSA(matQ_, y, h_, J_, random, invKT);
}
clearState(solSolutionAvailable);
}
Expand Down
2 changes: 1 addition & 1 deletion sqaodc/cpu/CPUDenseGraphAnnealer.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -66,7 +66,7 @@ class CPUDenseGraphAnnealer : public sq::DenseGraphAnnealer<real> {
void annealColoredPlane(real G, real beta);
void annealColoredPlaneParallel(real G, real beta);
/* simulated annealing */
void annealOneStepSANaive(real kT, real beta);
void annealOneStepSANaive(real kT, real _);

void syncBits();

Expand Down
14 changes: 7 additions & 7 deletions sqaodc/cuda/CUDABipartiteGraphAnnealer.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -513,22 +513,22 @@ __device__ __forceinline__ static void
deviceTryFlipSA(int gidx, int gidy, real *d_qAnneal, sq::SizeType qAnnealStride,
int N, int m, const real *d_Emat, sq::SizeType EmatStride,
const real *d_h, const real *d_realRand,
real Tnorm) {
real invKT) {
int iq = gidx;

if ((iq < N) && (gidy < m)) {
int im = gidy;
real q = d_qAnneal[im * qAnnealStride + iq];
real dE = real(2.) * q * (d_h[iq] + d_Emat[im * qAnnealStride + iq]);
real thresh = dE < real(0.) ? real(1.) : exp(-dE * Tnorm);
real thresh = dE < real(0.) ? real(1.) : exp(-dE * invKT);
if (thresh > d_realRand[N * gidy + iq])
d_qAnneal[im * qAnnealStride + iq] = - q;
}
}

template<class real> void CUDABipartiteGraphAnnealer<real>::
tryFlipSA(DeviceMatrix *d_qAnneal, const DeviceMatrix &d_Jq, int N, int m,
const DeviceVector &d_h, const real *d_realRand, real Tnorm) {
const DeviceVector &d_h, const real *d_realRand, real invKT) {

real *d_qAnneal_data = d_qAnneal->d_data;
sq::SizeType qAnnealStride = d_qAnneal->stride;
Expand All @@ -541,7 +541,7 @@ tryFlipSA(DeviceMatrix *d_qAnneal, const DeviceMatrix &d_Jq, int N, int m,
auto flipOp0 = [=]__device__(int gidx, int gidy) {
deviceTryFlipSA(gidx, gidy,
d_qAnneal_data, qAnnealStride, N, m, d_Emat, EmatStride,
d_h_data, d_realRand, Tnorm);
d_h_data, d_realRand, invKT);
};
transform2d(flipOp0, N, m, dim3(64, 2), stream);
}
Expand All @@ -557,18 +557,18 @@ void CUDABipartiteGraphAnnealer<real>::annealOneStepSA(real kT, real beta) {
d_randReal_.generate<real>(d_random_, nRequiredRandNum);
const real *d_randNum;

real Tnorm = kT * beta;
real invKT = real(1.) / kT;

/* 1st */
calculate_Jq(&d_Jq1_, d_J_, opNone, d_matq0_);
d_randNum = d_randReal_.acquire<real>(N1_ * m_);
tryFlipSA(&d_matq1_, d_Jq1_, N1_, m_, d_h1_, d_randNum, Tnorm);
tryFlipSA(&d_matq1_, d_Jq1_, N1_, m_, d_h1_, d_randNum, invKT);
DEBUG_SYNC;

/* 2nd */
calculate_Jq(&d_Jq0_, d_J_, opTranspose, d_matq1_);
d_randNum = d_randReal_.acquire<real>(N0_ * m_);
tryFlipSA(&d_matq0_, d_Jq0_, N0_, m_, d_h0_, d_randNum, Tnorm);
tryFlipSA(&d_matq0_, d_Jq0_, N0_, m_, d_h0_, d_randNum, invKT);
DEBUG_SYNC;
}

Expand Down
16 changes: 8 additions & 8 deletions sqaodc/cuda/CUDADenseGraphAnnealer.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -406,7 +406,7 @@ tryFlipSQAKernel(char *d_q, sq::SizeType qStride, const real *d_Jq,

int neibour0 = (y == 0) ? m - 1 : y - 1;
int neibour1 = (y == m - 1) ? 0 : y + 1;
real dE = twoDivM * (real)qyx * (d_Jq[y] + d_h[x]);
real dE = twoDivM * (real)qyx * (real(2.) * d_Jq[y] + d_h[x]);
dE -= (real)qyx * (d_q[qStride * neibour0 + x] + d_q[qStride * neibour1 + x]) * coef;
real threshold = (dE < real(0.)) ? real(1.) : exp(-dE * beta);
if (threshold > d_random[y])
Expand Down Expand Up @@ -508,7 +508,7 @@ __global__ static void
tryFlipSAKernel(char *d_q, sq::SizeType qStride, const real *d_Jq,
const real *d_h,
const int *d_x, const real *d_random, sq::SizeType m,
const real Tnorm) {
const real invKT) {

int gid = blockDim.x * blockIdx.x + threadIdx.x;
if (gid < m) {
Expand All @@ -517,7 +517,7 @@ tryFlipSAKernel(char *d_q, sq::SizeType qStride, const real *d_Jq,
char qyx = d_q[qStride * y + x];

real dE = real(2.) * (real)qyx * (d_Jq[y] + d_h[x]);
real threshold = (dE < real(0.)) ? real(1.) : exp(-dE * Tnorm);
real threshold = (dE < real(0.)) ? real(1.) : exp(-dE * invKT);
if (threshold > d_random[y])
d_q[qStride * y + x] = - qyx;
}
Expand All @@ -526,19 +526,19 @@ tryFlipSAKernel(char *d_q, sq::SizeType qStride, const real *d_Jq,
template<class real> void CUDADenseGraphAnnealer<real>::
annealOneStepSA(DeviceBitMatrix *d_matq, const DeviceVector &d_Jq, const int *d_x,
const real *d_random,
const DeviceVector &d_h, const DeviceMatrix &d_J, real Tnorm) {
const DeviceVector &d_h, const DeviceMatrix &d_J, real invKT) {
dim3 blockDim(128);
dim3 gridDim(std::max(divru(m_, blockDim.x), 1));
cudaStream_t stream = devStream_->getCudaStream();
#if 0
tryFlipSAKernel<real><<<gridDim, blockDim, 0, stream>>>(d_matq->d_data, d_matq->stride,
d_Jq.d_data, d_h.d_data,
d_x, d_random, m_,
Tnorm);
invKT);
#else
void *args[] = {(void*)&d_matq->d_data, (void*)&d_matq->stride,
(void*)&d_Jq.d_data, (void*)&d_h.d_data, (void*)&d_x, (void*)&d_random,
(void*)&m_, (void*)&Tnorm, NULL};
(void*)&m_, (void*)&invKT, NULL};
cudaLaunchKernel((void*)tryFlipSAKernel<real>, gridDim, blockDim, args, 0, stream);
#endif
DEBUG_SYNC;
Expand All @@ -555,12 +555,12 @@ void CUDADenseGraphAnnealer<real>::annealOneStepSA(real kT, real beta) {
if (!realNumBuffer_.available(m_ * N_))
realNumBuffer_.generate<real>(d_random_, N_ * m_ * nRunsPerRandGen_);

real Tnorm = kT * beta;
real invKT = real(1.) / kT;
for (int idx = 0; idx < N_; ++idx) {
const int *d_flipPos = flipPosBuffer_.acquire<int>(m_);
const real *d_random = realNumBuffer_.acquire<real>(m_);
calculate_Jq(&d_Jq_, d_J_, d_matq_, d_flipPos);
annealOneStepSA(&d_matq_, d_Jq_, d_flipPos, d_random, d_h_, d_J_, Tnorm);
annealOneStepSA(&d_matq_, d_Jq_, d_flipPos, d_random, d_h_, d_J_, invKT);
}
}

Expand Down
18 changes: 10 additions & 8 deletions sqaodpy/sqaod/py/bipartite_graph_annealer.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -336,12 +336,13 @@ def anneal_one_step_coloring(self, G, beta) :
self._anneal_half_step_coloring(N0, q0, h0, J.T, q1, G, beta, m)

# simulated annealing
def anneal_one_step_sa_naive(self, kT, beta) :
def anneal_one_step_sa_naive(self, kT, _) :
""" (sqaod.py only) sqaod.algorithm.sa_naive version of SA """
h0, h1, J, c, q0, q1 = self._vars()
N0, N1 = self.get_problem_size()
m = self._m
N = N0 + N1
invKT = 1. / kT

for im in range(m) :
q0m, q1m = q0[im],q1[im]
Expand All @@ -351,36 +352,37 @@ def anneal_one_step_sa_naive(self, kT, beta) :
if (iq < N0) :
q = q0m[iq]
dE = 2. * q * (h0[iq] + np.dot(J.T[iq], q1m))
thresh = 1 if dE < 0 else np.exp(-dE * kT * beta)
thresh = 1 if dE < 0 else np.exp(-dE * invKT)
if thresh > np.random.rand():
q0m[iq] = -q
else :
iq -= N0
q = q1m[iq]
dE = 2. * q * (h1[iq] + np.dot(J[iq], q0m))
thresh = 1 if dE < 0 else np.exp(-dE * kT * beta)
thresh = 1 if dE < 0 else np.exp(-dE * invKT)
if thresh > np.random.rand():
q1m[iq] = -q

def _anneal_half_step_sa_coloring(self, N, qAnneal, h, J, qFixed, kT, beta, m) :
def _anneal_half_step_sa_coloring(self, N, qAnneal, h, J, qFixed, invKT, m) :
""" (sqaod.py only) try to flip spins in a colored plane for SA. """
dEmat = np.matmul(J, qFixed.T)
for im in range(m) :
qAnnealm = qAnneal[im]
for iq in range(0, N) :
q = qAnnealm[iq]
dE = 2. * q * (h[iq] + dEmat[iq, im])
thresh = 1 if dE < 0 else np.exp(-dE * kT * beta)
thresh = 1 if dE < 0 else np.exp(-dE * invKT)
if thresh > np.random.rand():
qAnnealm[iq] = -q

def anneal_one_step_sa_coloring(self, kT, beta) :
def anneal_one_step_sa_coloring(self, kT, _) :
""" (sqaod.py only) sqaod.algorithm.sa_coloring version of SA """
h0, h1, J, c, q0, q1 = self._vars()
N0, N1 = self.get_problem_size()
m = self._m
self._anneal_half_step_sa_coloring(N1, q1, h1, J, q0, kT, beta, m)
self._anneal_half_step_sa_coloring(N0, q0, h0, J.T, q1, kT, beta, m)
invKT = 1. / kT
self._anneal_half_step_sa_coloring(N1, q1, h1, J, q0, invKT, m)
self._anneal_half_step_sa_coloring(N0, q0, h0, J.T, q1, invKT, m)



Expand Down
5 changes: 3 additions & 2 deletions sqaodpy/sqaod/py/dense_graph_annealer.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -319,10 +319,11 @@ def anneal_one_step_coloring(self, G, beta) :
self.anneal_colored_plane(G, beta, 0)
self.anneal_colored_plane(G, beta, 1)

def anneal_one_step_sa_naive(self, kT, beta) :
def anneal_one_step_sa_naive(self, kT, _) :
""" (sqaod.py only) sqaod.algorithm.sa_naive version of SA """
h, J, c, q = self._vars()
N = self._N
invKT = 1. / kT

for iq in range(self._m) :
qm = q[iq]
Expand All @@ -331,7 +332,7 @@ def anneal_one_step_sa_naive(self, kT, beta) :
qx = qm[x]
sum = np.dot(J[x], qm); # diagnoal elements in J are zero.
dE = 2. * qx * (h[x] + sum)
threshold = 1. if (dE <= 0.) else np.exp(-dE * kT * beta)
threshold = 1. if (dE <= 0.) else np.exp(- dE * invKT)
if threshold > np.random.rand():
qm[x] = - qx

Expand Down

0 comments on commit 25de2d8

Please sign in to comment.