diff --git a/dimod/discrete/cppdqm.pxd b/dimod/discrete/cppdqm.pxd new file mode 100644 index 000000000..4a21f715a --- /dev/null +++ b/dimod/discrete/cppdqm.pxd @@ -0,0 +1,100 @@ +# distutils: language = c++ +# cython: language_level=3 +# +# Copyright 2019 D-Wave Systems Inc. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# +# ============================================================================= + +from dimod.bqm.cppbqm cimport AdjVectorBQM as cppAdjVectorBQM +from dimod.bqm.common cimport VarIndex, Bias +from libcpp cimport bool +from libcpp.pair cimport pair +from libcpp.vector cimport vector + +cimport numpy as np + + +cdef extern from "dimod/adjvectordqm.h" namespace "dimod" nogil: + + cdef cppclass AdjVectorDQM[V, B]: + ctypedef V variable_type + ctypedef B bias_type + ctypedef size_t size_type + + cppAdjVectorBQM[variable_type, bias_type] bqm_ + vector[variable_type] case_starts_ + vector[vector[variable_type]] adj_ + + # constructors + # cython cannot handle templated constructors, so we call out the types + # explicitly + + AdjVectorDQM() except + + AdjVectorDQM(const AdjVectorDQM & ) except + + + # the actual signature is more general, but we already have a large + # number of these so we'll add them as needed + # template + # AdjVectorDQM(io_variable_type *case_starts, size_type num_variables, + # io_bias_type *linear_biases, size_type num_cases, + # io_variable_type *irow, io_variable_type *icol, + # io_bias_type *quadratic_biases, size_type num_interactions) + AdjVectorDQM(np.uint32_t*, size_type, np.uint32_t*, size_type, np.uint32_t*, np.uint32_t*, np.uint32_t*, size_type) + AdjVectorDQM(np.uint32_t*, size_type, np.uint64_t*, size_type, np.uint32_t*, np.uint32_t*, np.uint64_t*, size_type) + AdjVectorDQM(np.uint32_t*, size_type, np.int32_t*, size_type, np.uint32_t*, np.uint32_t*, np.int32_t*, size_type) + AdjVectorDQM(np.uint32_t*, size_type, np.int64_t*, size_type, np.uint32_t*, np.uint32_t*, np.int64_t*, size_type) + AdjVectorDQM(np.uint32_t*, size_type, np.float32_t*, size_type, np.uint32_t*, np.uint32_t*, np.float32_t*, size_type) + AdjVectorDQM(np.uint32_t*, size_type, np.float64_t*, size_type, np.uint32_t*, np.uint32_t*, np.float64_t*, size_type) + AdjVectorDQM(np.uint64_t*, size_type, np.uint32_t*, size_type, np.uint64_t*, np.uint64_t*, np.uint32_t*, size_type) + AdjVectorDQM(np.uint64_t*, size_type, np.uint64_t*, size_type, np.uint64_t*, np.uint64_t*, np.uint64_t*, size_type) + AdjVectorDQM(np.uint64_t*, size_type, np.int32_t*, size_type, np.uint64_t*, np.uint64_t*, np.int32_t*, size_type) + AdjVectorDQM(np.uint64_t*, size_type, np.int64_t*, size_type, np.uint64_t*, np.uint64_t*, np.int64_t*, size_type) + AdjVectorDQM(np.uint64_t*, size_type, np.float32_t*, size_type, np.uint64_t*, np.uint64_t*, np.float32_t*, size_type) + AdjVectorDQM(np.uint64_t*, size_type, np.float64_t*, size_type, np.uint64_t*, np.uint64_t*, np.float64_t*, size_type) + AdjVectorDQM(np.int32_t*, size_type, np.uint32_t*, size_type, np.int32_t*, np.int32_t*, np.uint32_t*, size_type) + AdjVectorDQM(np.int32_t*, size_type, np.uint64_t*, size_type, np.int32_t*, np.int32_t*, np.uint64_t*, size_type) + AdjVectorDQM(np.int32_t*, size_type, np.int32_t*, size_type, np.int32_t*, np.int32_t*, np.int32_t*, size_type) + AdjVectorDQM(np.int32_t*, size_type, np.int64_t*, size_type, np.int32_t*, np.int32_t*, np.int64_t*, size_type) + AdjVectorDQM(np.int32_t*, size_type, np.float32_t*, size_type, np.int32_t*, np.int32_t*, np.float32_t*, size_type) + AdjVectorDQM(np.int32_t*, size_type, np.float64_t*, size_type, np.int32_t*, np.int32_t*, np.float64_t*, size_type) + AdjVectorDQM(np.int64_t*, size_type, np.uint32_t*, size_type, np.int64_t*, np.int64_t*, np.uint32_t*, size_type) + AdjVectorDQM(np.int64_t*, size_type, np.uint64_t*, size_type, np.int64_t*, np.int64_t*, np.uint64_t*, size_type) + AdjVectorDQM(np.int64_t*, size_type, np.int32_t*, size_type, np.int64_t*, np.int64_t*, np.int32_t*, size_type) + AdjVectorDQM(np.int64_t*, size_type, np.int64_t*, size_type, np.int64_t*, np.int64_t*, np.int64_t*, size_type) + AdjVectorDQM(np.int64_t*, size_type, np.float32_t*, size_type, np.int64_t*, np.int64_t*, np.float32_t*, size_type) + AdjVectorDQM(np.int64_t*, size_type, np.float64_t*, size_type, np.int64_t*, np.int64_t*, np.float64_t*, size_type) + + # methods + + bool self_loop_present() except + + bool connection_present(variable_type, variable_type) except + + void connect_variables(variable_type, variable_type) except + + size_type num_variables() except + + size_type num_variable_interactions() except + + size_type num_cases(variable_type) except + + size_type num_case_interactions() except + + bias_type & linear_case(variable_type, variable_type) except + + void get_linear[io_bias_type](variable_type, io_bias_type*) except + + void set_linear[io_bias_type](variable_type, io_bias_type*) except + + pair[bias_type, bool] get_quadratic_case(variable_type, variable_type, variable_type, variable_type) except + + bool set_quadratic_case(variable_type, variable_type, variable_type, variable_type, bias_type) except + + bool get_quadratic[io_bias_type](variable_type, variable_type, io_bias_type*) except + + bool set_quadratic[io_bias_type](variable_type, variable_type, io_bias_type*) except + + double get_energy[io_variable_type](io_variable_type*) except + + void to_coo[io_variable_type, io_bias_type](io_variable_type*, io_bias_type*, io_variable_type*, io_variable_type*, io_bias_type*) except + + + # shapeable methods + + variable_type add_variable(variable_type) except + diff --git a/dimod/discrete/cydiscrete_quadratic_model.pxd b/dimod/discrete/cydiscrete_quadratic_model.pxd index 76f842fdc..2a4d8e127 100644 --- a/dimod/discrete/cydiscrete_quadratic_model.pxd +++ b/dimod/discrete/cydiscrete_quadratic_model.pxd @@ -20,7 +20,7 @@ from libcpp.vector cimport vector cimport numpy as np -from dimod.bqm.cppbqm cimport AdjVectorBQM as cppAdjVectorBQM +from dimod.discrete.cppdqm cimport AdjVectorDQM as cppAdjVectorDQM from dimod.bqm.common cimport Integral32plus, Numeric, Numeric32plus @@ -36,9 +36,7 @@ ctypedef fused Unsigned: cdef class cyDiscreteQuadraticModel: - cdef cppAdjVectorBQM[CaseIndex, Bias] bqm_ - cdef vector[CaseIndex] case_starts_ # len(adj_) + 1 - cdef vector[vector[VarIndex]] adj_ + cdef cppAdjVectorDQM[VarIndex, Bias] dqm_ cdef readonly object dtype cdef readonly object case_dtype diff --git a/dimod/discrete/cydiscrete_quadratic_model.pyx b/dimod/discrete/cydiscrete_quadratic_model.pyx index 41687879b..07840b63a 100644 --- a/dimod/discrete/cydiscrete_quadratic_model.pyx +++ b/dimod/discrete/cydiscrete_quadratic_model.pyx @@ -28,14 +28,12 @@ from dimod.utilities import asintegerarrays, asnumericarrays cdef class cyDiscreteQuadraticModel: def __init__(self): - self.case_starts_.push_back(0) - self.dtype = np.float64 self.case_dtype = np.int64 @property def adj(self): - return self.adj_ + return self.dqm_.adj_ cpdef Py_ssize_t add_variable(self, Py_ssize_t num_cases) except -1: """Add a discrete variable. @@ -51,69 +49,41 @@ cdef class cyDiscreteQuadraticModel: if num_cases <= 0: raise ValueError("num_cases must be a positive integer") - cdef VarIndex v = self.adj_.size() # index of new variable - - self.adj_.resize(v+1) - - cdef Py_ssize_t i - for i in range(num_cases): - self.bqm_.add_variable() - - self.case_starts_.push_back(self.bqm_.num_variables()) - - return v + return self.dqm_.add_variable(num_cases) def copy(self): cdef cyDiscreteQuadraticModel dqm = type(self)() - dqm.bqm_ = self.bqm_ - dqm.case_starts_ = self.case_starts_ - dqm.adj_ = self.adj_ + dqm.dqm_ = self.dqm_ dqm.dtype = self.dtype dqm.case_dtype = self.dtype + return dqm @cython.boundscheck(False) @cython.wraparound(False) cpdef Bias[:] energies(self, CaseIndex[:, :] samples): - + if samples.shape[1] != self.num_variables(): raise ValueError("Given sample(s) have incorrect number of variables") cdef Py_ssize_t num_samples = samples.shape[0] cdef VarIndex num_variables = samples.shape[1] - cdef Bias[:] energies = np.zeros(num_samples, dtype=self.dtype) + cdef Bias[:] energies = np.empty(num_samples, dtype=self.dtype) cdef Py_ssize_t si, vi cdef CaseIndex cu, case_u, cv, case_v cdef VarIndex u, v - for si in range(num_samples): # this could be parallelized + for si in range(num_samples): for u in range(num_variables): case_u = samples[si, u] if case_u >= self.num_cases(u): raise ValueError("invalid case") - cu = self.case_starts_[u] + case_u - - energies[si] += self.bqm_.get_linear(cu) - - for vi in range(self.adj_[u].size()): - v = self.adj_[u][vi] - - # we only care about the lower triangle - if v > u: - break - - case_v = samples[si, v] - - cv = self.case_starts_[v] + case_v - - out = self.bqm_.get_quadratic(cu, cv) - - if out.second: - energies[si] += out.first + for si in range(num_samples): + energies[si] = self.dqm_.get_energy(& samples[si, 0]) return energies @@ -160,68 +130,15 @@ cdef class cyDiscreteQuadraticModel: if irow[qi] == icol[qi]: raise ValueError("quadratic data contains a self-loop") - cdef cyDiscreteQuadraticModel dqm = cls() - - # set the BQM - if num_interactions: - dqm.bqm_ = cppAdjVectorBQM[CaseIndex, Bias]( - &irow[0], &icol[0], &quadratic_biases[0], - num_interactions, True) - - # add the linear biases - while dqm.bqm_.num_variables() < num_cases: - dqm.bqm_.add_variable() - for ci in range(num_cases): - dqm.bqm_.set_linear(ci, linear_biases[ci]) - - # set the case starts - dqm.case_starts_.resize(case_starts.shape[0] + 1) - for v in range(case_starts.shape[0]): - dqm.case_starts_[v] = case_starts[v] - dqm.case_starts_[case_starts.shape[0]] = dqm.bqm_.num_variables() - - # and finally the adj. This is not really the memory bottleneck so - # we can build an intermediate (unordered) set version - cdef vector[unordered_set[VarIndex]] adjset - adjset.resize(num_variables) - u = 0 - for ci in range(dqm.bqm_.num_variables()): - - # we've been careful so don't need ui < case_starts.size() - 1 - while ci >= dqm.case_starts_[u+1]: - u += 1 - - span = dqm.bqm_.neighborhood(ci) - - v = 0 - while span.first != span.second: - cj = deref(span.first).first - - # see above note - while cj >= dqm.case_starts_[v+1]: - v += 1 - - adjset[u].insert(v) - - inc(span.first) - - # now put adjset into adj - dqm.adj_.resize(num_variables) - for v in range(num_variables): - dqm.adj_[v].insert(dqm.adj_[v].begin(), - adjset[v].begin(), adjset[v].end()) - sort(dqm.adj_[v].begin(), dqm.adj_[v].end()) - - # do one last final check for self-loops within a variable - for v in range(num_variables): - for ci in range(dqm.case_starts_[v], dqm.case_starts_[v+1]): - span2 = dqm.bqm_.neighborhood(ci, dqm.case_starts_[v]) - if span2.first == span2.second: - continue - if deref(span2.first).first < dqm.case_starts_[v+1]: - raise ValueError("A variable has a self-loop") + cdef cyDiscreteQuadraticModel cyDQM = cls() - return dqm + cyDQM.dqm_ = cppAdjVectorDQM[VarIndex, Bias](&case_starts[0], num_variables, &linear_biases[0], + num_cases, &irow[0], &icol[0], &quadratic_biases[0], num_interactions) + + if cyDQM.dqm_.self_loop_present(): + raise ValueError("A variable has a self-loop") + + return cyDQM @classmethod def from_numpy_vectors(cls, case_starts, linear_biases, quadratic): @@ -252,13 +169,9 @@ cdef class cyDiscreteQuadraticModel: cdef Py_ssize_t num_cases = self.num_cases(v) - biases = np.empty(num_cases, dtype=np.float64) + biases = np.empty(num_cases, dtype=self.dtype) cdef Bias[:] biases_view = biases - - cdef Py_ssize_t c - for c in range(num_cases): - biases_view[c] = self.bqm_.get_linear(self.case_starts_[v] + c) - + self.dqm_.get_linear(v, & biases_view[0]) return biases @cython.boundscheck(False) @@ -271,18 +184,17 @@ cdef class cyDiscreteQuadraticModel: raise ValueError("case {} is invalid, variable only supports {} " "cases".format(case, self.num_cases(v))) - return self.bqm_.get_linear(self.case_starts_[v] + case) + return self.dqm_.linear_case(v, case) def get_quadratic(self, VarIndex u, VarIndex v, bint array=False): # check that the interaction does in fact exist - if u < 0 or u >= self.adj_.size(): + if u < 0 or u >= self.num_variables(): raise ValueError("unknown variable") - if v < 0 or v >= self.adj_.size(): + if v < 0 or v >= self.num_variables(): raise ValueError("unknown variable") - it = lower_bound(self.adj_[u].begin(), self.adj_[u].end(), v) - if it == self.adj_[u].end() or deref(it) != v: + if not self.dqm_.connection_present(u,v): raise ValueError("there is no interaction between given variables") cdef CaseIndex ci @@ -291,32 +203,23 @@ cdef class cyDiscreteQuadraticModel: if array: # build a numpy array - quadratic = np.zeros((self.num_cases(u), self.num_cases(v)), + quadratic = np.empty((self.num_cases(u), self.num_cases(v)), dtype=self.dtype) quadratic_view = quadratic - - for ci in range(self.case_starts_[u], self.case_starts_[u+1]): - - span = self.bqm_.neighborhood(ci, self.case_starts_[v]) - - while (span.first != span.second and deref(span.first).first < self.case_starts_[v+1]): - case_u = ci - self.case_starts_[u] - case_v = deref(span.first).first - self.case_starts_[v] - quadratic_view[case_u, case_v] = deref(span.first).second - - inc(span.first) + self.dqm_.get_quadratic(u, v, & quadratic_view[0, 0]) else: # store in a dict quadratic = {} - for ci in range(self.case_starts_[u], self.case_starts_[u+1]): + for ci in range(self.dqm_.case_starts_[u], self.dqm_.case_starts_[u+1]): - span = self.bqm_.neighborhood(ci, self.case_starts_[v]) + span = self.dqm_.bqm_.neighborhood( + ci, self.dqm_.case_starts_[v]) - while (span.first != span.second and deref(span.first).first < self.case_starts_[v+1]): - case_u = ci - self.case_starts_[u] - case_v = deref(span.first).first - self.case_starts_[v] + while (span.first != span.second and deref(span.first).first < self.dqm_.case_starts_[v+1]): + case_u = ci - self.dqm_.case_starts_[u] + case_v = deref(span.first).first - self.dqm_.case_starts_[v] quadratic[case_u, case_v] = deref(span.first).second inc(span.first) @@ -329,7 +232,7 @@ cdef class cyDiscreteQuadraticModel: @cython.wraparound(False) cpdef Bias get_quadratic_case(self, VarIndex u, CaseIndex case_u, - VarIndex v, CaseIndex case_v) except? -45.3: + VarIndex v, CaseIndex case_v) except? -45.3: if case_u < 0 or case_u >= self.num_cases(u): raise ValueError("case {} is invalid, variable only supports {} " @@ -339,11 +242,7 @@ cdef class cyDiscreteQuadraticModel: raise ValueError("case {} is invalid, variable only supports {} " "cases".format(case_v, self.num_cases(v))) - - cdef CaseIndex cu = self.case_starts_[u] + case_u - cdef CaseIndex cv = self.case_starts_[v] + case_v - - return self.bqm_.get_quadratic(cu, cv).first + return self.dqm_.get_quadratic_case(u, case_u, v, case_v).first @cython.boundscheck(False) @cython.wraparound(False) @@ -352,30 +251,26 @@ cdef class cyDiscreteQuadraticModel: the total number of cases in the DQM. """ if v < 0: - return self.bqm_.num_variables() + return self.dqm_.bqm_.num_variables() if v >= self.num_variables(): raise ValueError("unknown variable {}".format(v)) - return self.case_starts_[v+1] - self.case_starts_[v] + return self.dqm_.num_cases(v) cpdef Py_ssize_t num_case_interactions(self): """The total number of case interactions.""" - return self.bqm_.num_interactions() + return self.dqm_.num_case_interactions() @cython.boundscheck(False) @cython.wraparound(False) cpdef Py_ssize_t num_variable_interactions(self) except -1: """The total number of case interactions.""" - cdef Py_ssize_t num = 0 - cdef Py_ssize_t v - for v in range(self.num_variables()): - num += self.adj_[v].size() - return num // 2 + return self.dqm_.num_variable_interactions() cpdef Py_ssize_t num_variables(self): """The number of discrete variables in the DQM.""" - return self.adj_.size() + return self.dqm_.num_variables() @cython.boundscheck(False) @cython.wraparound(False) @@ -387,14 +282,12 @@ cdef class cyDiscreteQuadraticModel: raise ValueError('Recieved {} bias(es) for a variable of degree {}' ''.format(biases.shape[0], self.num_cases(v))) - cdef Py_ssize_t c - for c in range(biases.shape[0]): - self.bqm_.set_linear(self.case_starts_[v] + c, biases[c]) + self.dqm_.set_linear(v, & biases[0]) @cython.boundscheck(False) @cython.wraparound(False) cpdef Py_ssize_t set_linear_case(self, VarIndex v, CaseIndex case, Bias b) except -1: - + # self.num_cases checks that the variable is valid if case < 0: @@ -403,14 +296,15 @@ cdef class cyDiscreteQuadraticModel: raise ValueError("case {} is invalid, variable only supports {} " "cases".format(case, self.num_cases(v))) - self.bqm_.set_linear(self.case_starts_[v] + case, b) + cdef Bias *bias_ptr = &(self.dqm_.linear_case(v, case)) + bias_ptr[0] = b def set_quadratic(self, VarIndex u, VarIndex v, biases): # check that the interaction does in fact exist - if u < 0 or u >= self.adj_.size(): + if u < 0 or u >= self.num_variables(): raise ValueError("unknown variable") - if v < 0 or v >= self.adj_.size(): + if v < 0 or v >= self.num_variables(): raise ValueError("unknown variable") if u == v: raise ValueError("there cannot be a quadratic interaction between " @@ -435,32 +329,14 @@ cdef class cyDiscreteQuadraticModel: raise ValueError("case {} is invalid, variable only supports {} " "cases".format(case_v, self.num_cases(v))) - cu = self.case_starts_[u] + case_u - cv = self.case_starts_[v] + case_v + self.dqm_.set_quadratic_case(u, case_u, v, case_v, bias) - self.bqm_.set_quadratic(cu, cv, bias) + self.dqm_.connect_variables(u,v) else: - - biases_view = np.asarray(biases, dtype=self.dtype).reshape(num_cases_u, num_cases_v) - - for case_u in range(biases_view.shape[0]): - cu = self.case_starts_[u] + case_u - for case_v in range(biases_view.shape[1]): - cv = self.case_starts_[v] + case_v - - bias = biases_view[case_u, case_v] - - if bias: - self.bqm_.set_quadratic(cu, cv, bias) - # track in adjacency - low = lower_bound(self.adj_[u].begin(), self.adj_[u].end(), v) - if low == self.adj_[u].end() or deref(low) != v: - # need to add - self.adj_[u].insert(low, v) - self.adj_[v].insert( - lower_bound(self.adj_[v].begin(), self.adj_[v].end(), u), - u) + biases_view = np.asarray(biases, dtype=self.dtype).reshape( + num_cases_u, num_cases_v) + self.dqm_.set_quadratic(u, v, & biases_view[0, 0]) cpdef Py_ssize_t set_quadratic_case(self, VarIndex u, CaseIndex case_u, @@ -481,20 +357,7 @@ cdef class cyDiscreteQuadraticModel: raise ValueError("there cannot be a quadratic interaction between " "two cases in the same variable") - - cdef CaseIndex cu = self.case_starts_[u] + case_u - cdef CaseIndex cv = self.case_starts_[v] + case_v - - self.bqm_.set_quadratic(cu, cv, bias) - - # track in adjacency - low = lower_bound(self.adj_[u].begin(), self.adj_[u].end(), v) - if low == self.adj_[u].end() or deref(low) != v: - # need to add - self.adj_[u].insert(low, v) - self.adj_[v].insert( - lower_bound(self.adj_[v].begin(), self.adj_[v].end(), u), - u) + self.dqm_.set_quadratic_case(u, case_u, v, case_v, bias) @cython.boundscheck(False) @cython.wraparound(False) @@ -502,32 +365,13 @@ cdef class cyDiscreteQuadraticModel: Unsigned[:] irow, Unsigned[:] icol, Bias[:] qdata): # we don't do array length checking so be careful! This can segfault # if the given arrays are incorrectly sized - - cdef Py_ssize_t vi - for vi in range(self.num_variables()): - starts[vi] = self.case_starts_[vi] - - cdef Py_ssize_t ci = 0 - cdef Py_ssize_t qi = 0 - for ci in range(self.bqm_.num_variables()): - ldata[ci] = self.bqm_.linear(ci) - - span = self.bqm_.neighborhood(ci) - while span.first != span.second and deref(span.first).first < ci: - - irow[qi] = ci - icol[qi] = deref(span.first).first - qdata[qi] = deref(span.first).second - - inc(span.first) - qi += 1 - + self.dqm_.to_coo( & starts[0], & ldata[0], & irow[0], & icol[0], & qdata[0]) def to_numpy_vectors(self): - + cdef Py_ssize_t num_variables = self.num_variables() cdef Py_ssize_t num_cases = self.num_cases() - cdef Py_ssize_t num_interactions = self.bqm_.num_interactions() + cdef Py_ssize_t num_interactions = self.num_case_interactions() # use the minimum sizes of the various index types. We combine for # variables and cases and exclude int8 to keep the total number of diff --git a/dimod/include/dimod/adjvectordqm.h b/dimod/include/dimod/adjvectordqm.h new file mode 100644 index 000000000..8fe197337 --- /dev/null +++ b/dimod/include/dimod/adjvectordqm.h @@ -0,0 +1,344 @@ +// Copyright 2020 D-Wave Systems Inc. +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#ifndef DIMOD_ADJVECTORDQM_H_ +#define DIMOD_ADJVECTORDQM_H_ + +#include +#include +#include +#include + +#include "dimod/adjvectorbqm.h" +#include "dimod/utils.h" + +namespace dimod { + +template +class AdjVectorDQM { + public: + using bias_type = B; + using variable_type = V; + using size_type = std::size_t; + + AdjVectorBQM bqm_; + std::vector case_starts_; // len(adj_) + 1 + std::vector> adj_; + + AdjVectorDQM() { case_starts_.push_back(0); } + + template + explicit AdjVectorDQM(const DQM &dqm) { + bqm_ = dqm.bqm_; + case_starts_.insert(case_starts_.begin(), dqm.case_starts_.begin(), + dqm.case_starts_.end()); + adj_.resize(dqm.adj_.size()); + variable_type num_variables = dqm.num_variables(); + for (variable_type v = 0; v < num_variables; v++) { + adj_[v].insert(adj_[v].begin(), dqm.adj_[v].begin(), + dqm.adj_[v].end()); + } + } + + template + AdjVectorDQM(io_variable_type *case_starts, size_type num_variables, + io_bias_type *linear_biases, size_type num_cases, + io_variable_type *irow, io_variable_type *icol, + io_bias_type *quadratic_biases, size_type num_interactions) { + // Set the BQM, linear biases will be added separately. + if (num_interactions) { + bqm_ = AdjVectorBQM( + irow, icol, quadratic_biases, num_interactions, true); + } + + // Accounting for the cases/variables at the end without interaction. + while (bqm_.num_variables() < num_cases) { + bqm_.add_variable(); + } + assert(bqm_.num_variables() == num_cases); + + // Add the linear biases. + for (variable_type ci = 0; ci < num_cases; ci++) { + bqm_.set_linear(ci, linear_biases[ci]); + } + + // Set the case starts. + case_starts_.resize(num_variables + 1); + for (variable_type v = 0; v < num_variables; v++) { + case_starts_[v] = case_starts[v]; + } + case_starts_[num_variables] = num_cases; + + // Fill the adjacency list for variables. + std::vector> adjset; + adjset.resize(num_variables); + variable_type u = 0; + variable_type num_total_cases = bqm_.num_variables(); + for (variable_type ci = 0; ci < num_total_cases; ci++) { + while (ci >= case_starts_[u + 1]) { + u++; + } + auto span = bqm_.neighborhood(ci); + variable_type v = 0; + while (span.first != span.second) { + variable_type cj = (span.first)->first; + while (cj >= case_starts_[v + 1]) { + v++; + } + adjset[u].insert(v); + span.first++; + } + } + + adj_.resize(num_variables); +#pragma omp parallel for + for (variable_type v = 0; v < num_variables; v++) { + adj_[v].insert(adj_[v].begin(), adjset[v].begin(), adjset[v].end()); + std::sort(adj_[v].begin(), adj_[v].end()); + } + } + + template + void to_coo(io_variable_type *case_starts, io_bias_type *linear_biases, + io_variable_type *irow, io_variable_type *icol, + io_bias_type *quadratic_biases) { + variable_type num_variables = this->num_variables(); + variable_type num_total_cases = bqm_.num_variables(); + + for (variable_type v = 0; v < num_variables; v++) { + case_starts[v] = case_starts_[v]; + } + + size_type qi = 0; + for (variable_type ci = 0; ci < num_total_cases; ci++) { + linear_biases[ci] = bqm_.get_linear(ci); + auto span = bqm_.neighborhood(ci); + while ((span.first != span.second) && ((span.first)->first < ci)) { + irow[qi] = ci; + icol[qi] = (span.first)->first; + quadratic_biases[qi] = (span.first)->second; + span.first++; + qi++; + } + } + } + + bool self_loop_present() { + variable_type num_variables = this->num_variables(); + for (variable_type v = 0; v < num_variables; v++) { + for (variable_type ci = case_starts_[v], + ci_end = case_starts_[v + 1]; + ci < ci_end; ci++) { + auto span = bqm_.neighborhood(ci, case_starts_[v]); + if ((span.first != span.second) && + ((span.first)->first < case_starts_[v + 1])) { + return true; + } + } + } + return false; + } + + bool connection_present(variable_type u, variable_type v) { + bool connected = true; + auto it = std::lower_bound(adj_[u].begin(), adj_[u].end(), v); + if (it == adj_[u].end() || *it != v) { + connected = false; + } + return connected; + } + + void connect_variables(variable_type u, variable_type v) { + auto low = std::lower_bound(adj_[u].begin(), adj_[u].end(), v); + if (low == adj_[u].end() || *low != v) { + adj_[u].insert(low, v); + adj_[v].insert(std::lower_bound(adj_[v].begin(), adj_[v].end(), u), + u); + } + } + + size_type num_variables() { return adj_.size(); } + + size_type num_variable_interactions() { + size_type num = 0; + variable_type num_variables = this->num_variables(); + for (variable_type v = 0; v < num_variables; v++) { + num += adj_[v].size(); + } + return (num / 2); + } + + size_type num_cases(variable_type v) { + assert(v >= 0 && v < this->num_variables()); + return (case_starts_[v + 1] - case_starts_[v]); + } + + size_type num_case_interactions() { return bqm_.num_interactions(); } + + variable_type add_variable(variable_type num_cases) { + assert(num_cases > 0); + variable_type v = adj_.size(); + adj_.resize(v + 1); + for (variable_type n = 0; n < num_cases; n++) { + bqm_.add_variable(); + } + case_starts_.push_back(bqm_.num_variables()); + return v; + } + + bias_type &linear_case(variable_type v, variable_type case_v) { + assert(v >= 0 && v < this->num_variables()); + assert(case_v >= 0 && case_v < this->num_cases(v)); + return bqm_.linear(case_starts_[v] + case_v); + } + + template + void get_linear(variable_type v, io_bias_type *biases) { + assert(v >= 0 && v < this->num_variables()); + variable_type num_cases_v = this->num_cases(v); + for (variable_type case_v = 0; case_v < num_cases_v; case_v++) { + biases[case_v] = bqm_.linear(case_starts_[v] + case_v); + } + } + + template + void set_linear(variable_type v, io_bias_type *biases) { + assert(v >= 0 && v < this->num_variables()); + variable_type num_cases_v = this->num_cases(v); + for (variable_type case_v = 0; case_v < num_cases_v; case_v++) { + bqm_.linear(case_starts_[v] + case_v) = biases[case_v]; + } + } + + std::pair get_quadratic_case(variable_type u, + variable_type case_u, + variable_type v, + variable_type case_v) { + assert(u >= 0 && u < this->num_variables()); + assert(case_u >= 0 && case_u < this->num_cases(u)); + assert(v >= 0 && v < this->num_variables()); + assert(case_v >= 0 && case_v < this->num_cases(v)); + variable_type cu = case_starts_[u] + case_u; + variable_type cv = case_starts_[v] + case_v; + return bqm_.get_quadratic(cu, cv); + } + + // Check if boolean type is still okay + bool set_quadratic_case(variable_type u, variable_type case_u, + variable_type v, variable_type case_v, + bias_type bias) { + assert(u >= 0 && u < this->num_variables()); + assert(case_u >= 0 && case_u < this->num_cases(u)); + assert(v >= 0 && v < this->num_variables()); + assert(case_v >= 0 && case_v < this->num_cases(v)); + if (u == v) { + return false; + } + variable_type cu = case_starts_[u] + case_u; + variable_type cv = case_starts_[v] + case_v; + bqm_.set_quadratic(cu, cv, bias); + connect_variables(u, v); + return true; + } + + // Returns false if there is no interaction among the variables. + template + bool get_quadratic(variable_type u, variable_type v, + io_bias_type *quadratic_biases) { + assert(u >= 0 && u < this->num_variables()); + assert(v >= 0 && v < this->num_variables()); + if (!connection_present(u, v)) { + return false; + } + variable_type num_cases_u = num_cases(u); + variable_type num_cases_v = num_cases(v); + memset(quadratic_biases, 0, + (size_t)num_cases_u * (size_t)num_cases_v * + sizeof(io_bias_type)); +#pragma omp parallel for + for (variable_type case_u = 0; case_u < num_cases_u; case_u++) { + auto span = bqm_.neighborhood(case_starts_[u] + case_u, + case_starts_[v]); + while (span.first != span.second && + (span.first)->first < case_starts_[v + 1]) { + variable_type case_v = (span.first)->first - case_starts_[v]; + quadratic_biases[case_u * num_cases_v + case_v] = + (span.first)->second; + span.first++; + } + } + return true; + } + + template + bool set_quadratic(variable_type u, variable_type v, io_bias_type *biases) { + assert(u >= 0 && u < this->num_variables()); + assert(v >= 0 && v < this->num_variables()); + if (u == v) { + return false; + } + variable_type num_cases_u = num_cases(u); + variable_type num_cases_v = num_cases(v); + // This cannot be parallelized since the vectors cannot be reshaped in + // parallel. + bool inserted = false; + for (variable_type case_u = 0; case_u < num_cases_u; case_u++) { + variable_type cu = case_starts_[u] + case_u; + for (variable_type case_v = 0; case_v < num_cases_v; case_v++) { + variable_type cv = case_starts_[v] + case_v; + bias_type bias = biases[case_u * num_cases_v + case_v]; + if (bias != (bias_type)0) { + bqm_.set_quadratic(cu, cv, bias); + inserted = true; + } + } + } + + if (inserted) { + connect_variables(u, v); + return true; + } else { + return false; + } + } + + template + double get_energy(io_variable_type *sample) { + variable_type num_variables = this->num_variables(); + double sample_energy = 0; + for (variable_type u = 0; u < num_variables; u++) { + variable_type case_u = sample[u]; + assert(case_u < num_cases(u)); + variable_type cu = case_starts_[u] + case_u; + sample_energy += bqm_.get_linear(cu); + for (variable_type vi = 0; vi < adj_[u].size(); vi++) { + variable_type v = adj_[u][vi]; + // We only care about lower triangle. + if (v > u) { + break; + } + variable_type case_v = sample[v]; + variable_type cv = case_starts_[v] + case_v; + auto out = bqm_.get_quadratic(cu, cv); + if (out.second) { + sample_energy += out.first; + } + } + } + return sample_energy; + } +}; +} // namespace dimod + +#endif // DIMOD_ADJVECTORDQM_H_