From c2f4e26c51145c7f6fb23688bc2497042152a31d Mon Sep 17 00:00:00 2001 From: Alexander Condello Date: Wed, 26 May 2021 15:55:43 -0700 Subject: [PATCH 1/2] Comment out roof duality tests This commit should be reverted when the dwave-preprocessing optional dependency is added. --- tests/test_fixedvariablecomposite.py | 120 +++++++++++++-------------- tests/test_roof_duality.py | 18 ++-- 2 files changed, 66 insertions(+), 72 deletions(-) diff --git a/tests/test_fixedvariablecomposite.py b/tests/test_fixedvariablecomposite.py index 25a38c2b6..e28fe6f3f 100644 --- a/tests/test_fixedvariablecomposite.py +++ b/tests/test_fixedvariablecomposite.py @@ -26,85 +26,85 @@ from dimod import SampleSet -@dimod.testing.load_sampler_bqm_tests(FixedVariableComposite(ExactSolver())) -@dimod.testing.load_sampler_bqm_tests(FixedVariableComposite(dimod.NullSampler())) -class TestFixedVariableComposite(unittest.TestCase): +# @dimod.testing.load_sampler_bqm_tests(FixedVariableComposite(ExactSolver())) +# @dimod.testing.load_sampler_bqm_tests(FixedVariableComposite(dimod.NullSampler())) +# class TestFixedVariableComposite(unittest.TestCase): - def test_instantiation_smoketest(self): - sampler = FixedVariableComposite(ExactSolver()) +# def test_instantiation_smoketest(self): +# sampler = FixedVariableComposite(ExactSolver()) - dtest.assert_sampler_api(sampler) +# dtest.assert_sampler_api(sampler) - def test_sample(self): - bqm = BinaryQuadraticModel({1: -1.3, 4: -0.5}, - {(1, 4): -0.6}, - 0, - vartype=Vartype.SPIN) +# def test_sample(self): +# bqm = BinaryQuadraticModel({1: -1.3, 4: -0.5}, +# {(1, 4): -0.6}, +# 0, +# vartype=Vartype.SPIN) - fixed_variables = {1: -1} - sampler = FixedVariableComposite(ExactSolver()) - response = sampler.sample(bqm, fixed_variables=fixed_variables) +# fixed_variables = {1: -1} +# sampler = FixedVariableComposite(ExactSolver()) +# response = sampler.sample(bqm, fixed_variables=fixed_variables) - self.assertEqual(response.first.sample, {4: -1, 1: -1}) - self.assertAlmostEqual(response.first.energy, 1.2) +# self.assertEqual(response.first.sample, {4: -1, 1: -1}) +# self.assertAlmostEqual(response.first.energy, 1.2) - def test_empty_bqm(self): - bqm = BinaryQuadraticModel({1: -1.3, 4: -0.5}, - {(1, 4): -0.6}, - 0, - vartype=Vartype.SPIN) +# def test_empty_bqm(self): +# bqm = BinaryQuadraticModel({1: -1.3, 4: -0.5}, +# {(1, 4): -0.6}, +# 0, +# vartype=Vartype.SPIN) - fixed_variables = {1: -1, 4: -1} - sampler = FixedVariableComposite(ExactSolver()) - response = sampler.sample(bqm, fixed_variables=fixed_variables) - self.assertIsInstance(response, SampleSet) +# fixed_variables = {1: -1, 4: -1} +# sampler = FixedVariableComposite(ExactSolver()) +# response = sampler.sample(bqm, fixed_variables=fixed_variables) +# self.assertIsInstance(response, SampleSet) - def test_empty_fix(self): - linear = {1: -1.3, 4: -0.5} - quadratic = {(1, 4): -0.6} +# def test_empty_fix(self): +# linear = {1: -1.3, 4: -0.5} +# quadratic = {(1, 4): -0.6} - sampler = FixedVariableComposite(ExactSolver()) - response = sampler.sample_ising(linear, quadratic) - self.assertIsInstance(response, SampleSet) +# sampler = FixedVariableComposite(ExactSolver()) +# response = sampler.sample_ising(linear, quadratic) +# self.assertIsInstance(response, SampleSet) - self.assertEqual(response.first.sample, {4: 1, 1: 1}) - self.assertAlmostEqual(response.first.energy, -2.4) +# self.assertEqual(response.first.sample, {4: 1, 1: 1}) +# self.assertAlmostEqual(response.first.energy, -2.4) -@dimod.testing.load_sampler_bqm_tests(RoofDualityComposite(ExactSolver())) -@dimod.testing.load_sampler_bqm_tests(RoofDualityComposite(dimod.NullSampler())) -class TestRoofDualityComposite(unittest.TestCase): - def test_construction(self): - sampler = RoofDualityComposite(dimod.ExactSolver()) - dtest.assert_sampler_api(sampler) +# @dimod.testing.load_sampler_bqm_tests(RoofDualityComposite(ExactSolver())) +# @dimod.testing.load_sampler_bqm_tests(RoofDualityComposite(dimod.NullSampler())) +# class TestRoofDualityComposite(unittest.TestCase): +# def test_construction(self): +# sampler = RoofDualityComposite(dimod.ExactSolver()) +# dtest.assert_sampler_api(sampler) - def test_3path(self): - sampler = RoofDualityComposite(dimod.ExactSolver()) - sampleset = sampler.sample_ising({'a': 10}, {'ab': -1, 'bc': 1}) +# def test_3path(self): +# sampler = RoofDualityComposite(dimod.ExactSolver()) +# sampleset = sampler.sample_ising({'a': 10}, {'ab': -1, 'bc': 1}) - # all should be fixed, so should just see one - self.assertEqual(len(sampleset), 1) - self.assertEqual(set(sampleset.variables), set('abc')) +# # all should be fixed, so should just see one +# self.assertEqual(len(sampleset), 1) +# self.assertEqual(set(sampleset.variables), set('abc')) - def test_triangle(self): - sampler = RoofDualityComposite(dimod.ExactSolver()) +# def test_triangle(self): +# sampler = RoofDualityComposite(dimod.ExactSolver()) - bqm = dimod.BinaryQuadraticModel.from_ising({}, {'ab': -1, 'bc': -1, 'ac': -1}) +# bqm = dimod.BinaryQuadraticModel.from_ising({}, {'ab': -1, 'bc': -1, 'ac': -1}) - # two equally good solutions - sampleset = sampler.sample(bqm) +# # two equally good solutions +# sampleset = sampler.sample(bqm) - self.assertEqual(set(sampleset.variables), set('abc')) - dimod.testing.assert_response_energies(sampleset, bqm) +# self.assertEqual(set(sampleset.variables), set('abc')) +# dimod.testing.assert_response_energies(sampleset, bqm) - def test_triangle_sampling_mode_off(self): - sampler = RoofDualityComposite(dimod.ExactSolver()) +# def test_triangle_sampling_mode_off(self): +# sampler = RoofDualityComposite(dimod.ExactSolver()) - bqm = dimod.BinaryQuadraticModel.from_ising({}, {'ab': -1, 'bc': -1, 'ac': -1}) +# bqm = dimod.BinaryQuadraticModel.from_ising({}, {'ab': -1, 'bc': -1, 'ac': -1}) - # two equally good solutions, but with sampling mode off it will pick one - sampleset = sampler.sample(bqm, sampling_mode=False) +# # two equally good solutions, but with sampling mode off it will pick one +# sampleset = sampler.sample(bqm, sampling_mode=False) - self.assertEqual(set(sampleset.variables), set('abc')) - self.assertEqual(len(sampleset), 1) # all should be fixed - dimod.testing.assert_response_energies(sampleset, bqm) +# self.assertEqual(set(sampleset.variables), set('abc')) +# self.assertEqual(len(sampleset), 1) # all should be fixed +# dimod.testing.assert_response_energies(sampleset, bqm) diff --git a/tests/test_roof_duality.py b/tests/test_roof_duality.py index a5478bbec..6c3b93c71 100644 --- a/tests/test_roof_duality.py +++ b/tests/test_roof_duality.py @@ -17,17 +17,11 @@ import dimod -try: - from dimod import fix_variables -except ImportError: - cpp = False -else: - cpp = True +from dimod import fix_variables -@unittest.skipUnless(cpp, "no cpp extensions built") -class TestFixVariables(unittest.TestCase): - def test_3path(self): - bqm = dimod.BinaryQuadraticModel.from_ising({'a': 10}, {'ab': -1, 'bc': 1}) - fixed = dimod.fix_variables(bqm) - self.assertEqual(fixed, {'a': -1, 'b': -1, 'c': 1}) +# class TestFixVariables(unittest.TestCase): +# def test_3path(self): +# bqm = dimod.BinaryQuadraticModel.from_ising({'a': 10}, {'ab': -1, 'bc': 1}) +# fixed = dimod.fix_variables(bqm) +# self.assertEqual(fixed, {'a': -1, 'b': -1, 'c': 1}) From b4fa72127b24d51e1c5d57bdc25a30a74961eb04 Mon Sep 17 00:00:00 2001 From: Alexander Condello Date: Wed, 26 May 2021 15:56:31 -0700 Subject: [PATCH 2/2] Remove boost and roof duality code Closes https://github.com/dwavesystems/dimod/issues/618 Closes https://github.com/dwavesystems/dimod/issues/519 --- .circleci/config.yml | 30 - MANIFEST.in | 1 - README.rst | 5 - dimod/__init__.py | 9 +- dimod/roof_duality/_fix_variables.pyx | 57 - dimod/roof_duality/fix_variables.py | 56 +- dimod/roof_duality/src/compressed_matrix.hpp | 715 ---------- dimod/roof_duality/src/fix_variables.cpp | 1297 ------------------ dimod/roof_duality/src/fix_variables.hpp | 45 - setup.py | 5 - testscpp/Makefile | 4 +- testscpp/tests/test_roofduality.cpp | 44 - 12 files changed, 5 insertions(+), 2263 deletions(-) delete mode 100644 dimod/roof_duality/_fix_variables.pyx delete mode 100644 dimod/roof_duality/src/compressed_matrix.hpp delete mode 100644 dimod/roof_duality/src/fix_variables.cpp delete mode 100644 dimod/roof_duality/src/fix_variables.hpp delete mode 100644 testscpp/tests/test_roofduality.cpp diff --git a/.circleci/config.yml b/.circleci/config.yml index c254569da..08f3fdb7f 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -24,11 +24,6 @@ jobs: name: fix git checkout command: git reset --hard - - run: - name: install boost - command: | - yum install -y boost-devel - - run: name: build wheels command: | @@ -62,12 +57,6 @@ jobs: steps: - checkout - # in the future we'd like to not need boost here - - run: &linux-install-boost-template - name: install boost - command: | - sudo apt-get install libboost-dev - - run: name: build sdist command: | @@ -174,11 +163,6 @@ jobs: command: | brew install pyenv - - run: &install-boost-osx-template - name: install boost - command: | - brew install boost - - restore_cache: keys: - pyenv-{{ .Environment.CIRCLE_JOB }}-xcode12.2.0 @@ -247,8 +231,6 @@ jobs: command: | brew install pyenv - - run: *install-boost-osx-template - - restore_cache: keys: - pyenv-{{ .Environment.CIRCLE_JOB }}-xcode12.2.0 @@ -310,16 +292,10 @@ jobs: pip install -r requirements.txt pip install -r tests\requirements.txt - - run: &win-install-boost-template - name: install boost - command: | - nuget install boost -ExcludeVersion -OutputDirectory . - - run: name: build extension command: | env\Scripts\activate.ps1 - $env:CL+=' -Iboost\lib\native\include\' python setup.py build_ext --inplace - run: @@ -354,13 +330,10 @@ jobs: - run: *win-install-dependencies-template - - run: *win-install-boost-template - - run: name: create wheel command: | env\Scripts\activate.ps1 - $env:CL+=' -Iboost\lib\native\include\' python setup.py bdist_wheel - store_artifacts: @@ -383,8 +356,6 @@ jobs: steps: - checkout - - run: *linux-install-boost-template - - run: name: install doxygen command: sudo apt-get install doxygen @@ -429,7 +400,6 @@ jobs: steps: - checkout - - run: *linux-install-boost-template - run: name: run cpp tests command: | diff --git a/MANIFEST.in b/MANIFEST.in index 1241775d9..01cbc831b 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -1,4 +1,3 @@ include pyproject.toml -recursive-include dimod/roof_duality *.cpp *.hpp recursive-include dimod/include *h recursive-include dimod *.pyx *.pxd *.pyx.src diff --git a/README.rst b/README.rst index dee67c452..323808eeb 100644 --- a/README.rst +++ b/README.rst @@ -87,11 +87,6 @@ in place: pip install -r requirements.txt python setup.py build_ext --inplace -Note that installation from source requires that your system have the Boost_ -C++ libraries installed. - -.. _Boost: https://www.boost.org/ - .. installation-end-marker License diff --git a/dimod/__init__.py b/dimod/__init__.py index 404a96303..49d92d10d 100644 --- a/dimod/__init__.py +++ b/dimod/__init__.py @@ -23,14 +23,7 @@ from dimod.reference import * import dimod.reference -try: - import dimod.roof_duality._fix_variables as _ -except ImportError: - pass -else: - # we only import fix_variables function into the top level namespace if the c++ extension - # is built - from dimod.roof_duality import fix_variables +from dimod.roof_duality import fix_variables from dimod.binary_quadratic_model import BinaryQuadraticModel, BQM import dimod.binary_quadratic_model diff --git a/dimod/roof_duality/_fix_variables.pyx b/dimod/roof_duality/_fix_variables.pyx deleted file mode 100644 index 55afa93d7..000000000 --- a/dimod/roof_duality/_fix_variables.pyx +++ /dev/null @@ -1,57 +0,0 @@ -# distutils: language = c++ -# cython: language_level=3 -# -# Copyright 2018 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 libcpp.map cimport map -from libcpp.utility cimport pair -from libcpp.vector cimport vector - -from dimod.vartypes import Vartype - -cdef extern from "fix_variables.hpp" namespace "fix_variables_": - vector[pair[int, int]] fixQuboVariablesMap(map[pair[int, int], double] QMap, - int QSize, int method) except + - - -def fix_variables_wrapper(bqm, method): - """bqm should be binary and linear indexed - - method 1 -> roof-duality & strongly connected components - method 2 -> roof-duality only - - """ - - if bqm.vartype is not Vartype.BINARY: - raise ValueError("bqm must be BINARY") - if not all(v in bqm.linear for v in range(len(bqm))): - raise ValueError("bqm must be integer-labeled") - if not isinstance(method, int): - raise TypeError("method should be an int") - if method < 1 or method > 2: - raise ValueError("method should 1 or 2") - - cdef map[pair[int, int], double] QMap - for (u, v), bias in bqm.quadratic.items(): - QMap[pair[int, int](u, v)] = bias - for v, bias in bqm.linear.items(): - QMap[pair[int, int](v, v)] = bias - - fixed = fixQuboVariablesMap(QMap, len(bqm), int(method)) - - return {int(v - 1): int(val) for v, val in fixed} diff --git a/dimod/roof_duality/fix_variables.py b/dimod/roof_duality/fix_variables.py index 46b200ccd..50f4ee0c4 100644 --- a/dimod/roof_duality/fix_variables.py +++ b/dimod/roof_duality/fix_variables.py @@ -39,31 +39,6 @@ def fix_variables(bqm, sampling_mode=True): Returns: dict: Variable assignments for some variables of the specified binary quadratic model. - Examples: - This example creates a binary quadratic model with a single ground state - and fixes the model's single variable to the minimizing assignment. - - >>> bqm = dimod.BinaryQuadraticModel.from_ising({'a': 1.0}, {}) - >>> dimod.fix_variables(bqm) - {'a': -1} - - This example has two ground states, :math:`a=b=-1` and :math:`a=b=1`, with - no variable having a single value for all ground states, so neither variable - is fixed. - - >>> bqm = dimod.BinaryQuadraticModel.empty(dimod.SPIN) - >>> bqm.add_interaction('a', 'b', -1.0) - >>> dimod.fix_variables(bqm) # doctest: +SKIP - {} - - This example turns sampling model off, so variables are fixed to an assignment - that attains the ground state. - - >>> bqm = dimod.BinaryQuadraticModel.empty(dimod.SPIN) - >>> bqm.add_interaction('a', 'b', -1.0) - >>> dimod.fix_variables(bqm, sampling_mode=False) # doctest: +SKIP - {'a': 1, 'b': 1} - .. [BHT] Boros, E., P.L. Hammer, G. Tavares. Preprocessing of Unconstraint Quadratic Binary Optimization. Rutcor Research Report 10-2006, April, 2006. @@ -71,32 +46,5 @@ def fix_variables(bqm, sampling_mode=True): (2002), pp. 155-225 """ - try: - from dimod.roof_duality._fix_variables import fix_variables_wrapper - except ImportError: - raise ImportError("c++ extension roof_duality is not built") - - if sampling_mode: - method = 2 # roof-duality only - else: - method = 1 # roof-duality and strongly connected components - - linear = bqm.linear - if all(v in linear for v in range(len(bqm))): - # we can work with the binary form of the bqm directly - fixed = fix_variables_wrapper(bqm.binary, method) - else: - try: - inverse_mapping = dict(enumerate(sorted(linear))) - except TypeError: - # in python3 unlike types cannot be sorted - inverse_mapping = dict(enumerate(linear)) - mapping = {v: i for i, v in inverse_mapping.items()} - - fixed = fix_variables_wrapper(bqm.relabel_variables(mapping, inplace=False).binary, method) - fixed = {inverse_mapping[v]: val for v, val in fixed.items()} - - if bqm.vartype is Vartype.SPIN: - return {v: 2*val - 1 for v, val in fixed.items()} - else: - return fixed + # todo: add temporary dependency on dwave-preprocessing + raise NotImplementedError diff --git a/dimod/roof_duality/src/compressed_matrix.hpp b/dimod/roof_duality/src/compressed_matrix.hpp deleted file mode 100644 index dbaa773b7..000000000 --- a/dimod/roof_duality/src/compressed_matrix.hpp +++ /dev/null @@ -1,715 +0,0 @@ -/** -# Copyright 2018 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 COMPRESSED_MATRIX_HPP_INCLUDED -#define COMPRESSED_MATRIX_HPP_INCLUDED - -/** -* \file compressed_matrix.hpp -* -* Compressed matrix C++ header file. -*/ - -#include -#include -#include - -/** -* \namespace sapi -* \brief Namespace for compressed matrix. -*/ -namespace compressed_matrix -{ -// declaration - -/** -* \class CompressedMatrix -* \brief Sparse matrix. -* -* This class stores a (sparse) matrix in compressed-row format. -* -* The matrix is represented by three vectors: \p rowOffsets, \p colIndices, and \p values. -* \p rowOffsets has (fixed) length equal to the number of rows in the matrix, plus one. -* For any row index \p r and \p k such that rowOffsets[r] <= k < rowOffsets[r+1], -* the entry at postition (r, colIndices[k]) has value values[k]. Column indices are listed in -* increasing order within each row (ie. values in each row are listed from left to right). Missing column indices -* indicate a value of zero. -* -* Note that nothing prevents you from explicitly storing a zero value but this is wasteful. -*/ -// forward declarations -template -class CompressedMatrixIterator; - -template -class CompressedMatrixConstIterator; - -template -class CompressedMatrix -{ - friend class CompressedMatrixIterator; - friend class CompressedMatrixConstIterator; -public: - /** - * \brief Construct an empty CompressedMatrix of given size. - * - * \param numRows number of rows - * \param numCols number of columns - */ - CompressedMatrix(int numRows = 0, int numCols = 0); - - /** - * \brief Construct a CompressedMatrix from a map. - * - * For all valid row indices \p r and column indices \p c, the following is true of the new CompressedMatrix: - * \code get(r,c) == m[std::make_pair(r,c)] \endcode - * - * \param numRows number of rows - * \param numCols number of columns - * \param m sparse matrix represented as a map - * \throw ValueException if any key of \p m has an invalid row or column index - */ - CompressedMatrix(int numRows, int numCols, - const std::map, T>& m); - - /** - * \brief Construct a CompressedMatrix from raw compressed-row sparse matrix data. - * - * \param numRows number of rows in the matrix - * \param numCols number of columns in the matrix - * \param rowOffsets row offset data - * \param colIndices column index data - * \param values matrix values - * - * Requirements: - * - rowOffsets.size() == numRows + 1 - * - rowOffsets[0] == 0 - * - rowOffsets[i] <= rowOffsets[i + 1] for all valid \p i - * - colIndices[i] < numCols for all valid \p k - * - colIndices[k] < colIndices[k + 1] whenever rowIndices[i] <= k < rowIndices[i + 1] - * for all valid \p i - * - colIndices.size() == values.size() - */ - CompressedMatrix(int numRows, int numCols, - std::vector rowOffsets, std::vector colIndices, std::vector values); - - /** - * \return number of rows - */ - int numRows() const; - - /** - * \return number of columns - */ - int numCols() const; - - /** - * \brief Get the number of nonzero entries in the matrix. - * - * Note that if you set an entry to zero, it still counts as a "nonzero entry." Really, this function counts - * the number of not-zero-by-default entries. - * - * \return number of nonzero entries - */ - int nnz() const; - - /** - * \brief Look up a matrix entry. - * - * \param row row index of the entry - * \param col column index of the entry - * \return matrix entry (defaults to zero if not found) - * \throw ValueException if \p row is not less than numRows() or \p col is not less than numCols() - */ - T get(int row, int col) const; - - /** - * \brief Get a reference to matrix entry. - * - * Note that this function \em creates an entry (with value zero) if none exists! Consequently, this function can: - * - be slow, and - * - create unwanted zero entries. - * Use the get() for pure lookups! - * - * To minimize shuffling of data, fill the matrix from left to right, then top to bottom. - * - * \param row row index of the matrix entry - * \param col column index of the matrix entry - * \return reference to the matrix entry - * \throw ValueException if \p row is not less than numRows() or \p col is not less than numCols() - */ - T& operator()(int row, int col); - - /** - * \brief Get row offset data. - * - * This vector always has size numRows + 1. - * - * \return row offset data - */ - const std::vector& rowOffsets() const { return rowOffsets_; } - - /** - * \brief Get column index data. - * - * Values from index rowOffsets()[r] up to (but not including) index rowOffsets()[r + 1] are in strictly - * increasing order for any valid row index r. - * - * \return column index data - */ - const std::vector& colIndices() const { return colIndices_; } - - /** - * \return matrix entry values data - */ - const std::vector& values() const { return values_; } - - /** - * \brief Reserves space in colIndices and values. - * - * Calling this function before adding a large number of new entries can save on the number of memory reallocations. - * - * \param minValues new minimum size of colIndices and values. - */ - void reserve(int minValues); - - /** - * \p iterator type - */ - typedef CompressedMatrixIterator iterator; - - /** - * \p const_iterator type - */ - typedef CompressedMatrixConstIterator const_iterator; - - /** - * \Returns an iterator referring to the first element in the CompressedMatrix - */ - iterator begin(); - - /** - * \Returns a const_iterator referring to the first element in the CompressedMatrix - */ - const_iterator begin() const; - - /** - * \Returns an iterator referring to the past-the-end element in the CompressedMatrix - */ - iterator end(); - - /** - * \Returns a const_iterator referring to the past-the-end element in the CompressedMatrix - */ - const_iterator end() const; - -private: - int numRows_; - int numCols_; - - std::vector rowOffsets_; - std::vector colIndices_; - std::vector values_; -}; - - -/** -* \brief Converts a CompressedMatrix to a map. -* \param cm the CompressedMatrix to convert -* \return map m where m[make_pair(r,c)] == cm.get(r,c) for all r,c -*/ -template -std::map, T> compressedMatrixToMap(const CompressedMatrix& cm); - -/** -* \class CompressedMatrixIterator -* \brief Sparse matrix iterator. -* -* This class creates a sparse matrix iterator. -*/ -template -class CompressedMatrixIterator -{ -public: - friend class CompressedMatrixConstIterator; - /** - * \brief Construct a CompressedMatrix iterator - * - * \param cm a CompressedMatrix pointer - * \param currRow current row - * \param currIndex current index of colIndices and values - */ - CompressedMatrixIterator(CompressedMatrix* cm, int currRow = 0, int currIndex = 0); - - /** - * \brief Increment the iterator by one - */ - void operator++(); - - /** - * \brief Dereference the iterator - * - * \return a double reference - */ - T& operator*() const; - - /** - * \brief Get the row of the iterator - * - * \return row number - */ - int row() const; - - /** - * \brief Get the column of the iterator - * - * \return column number - */ - int col() const; - - /** - * \brief Compare whether two iterators are the same - * - * \return true if two iterators are the same; false otherwise - */ - bool operator!=(const CompressedMatrixIterator& iter) const; - -private: - - CompressedMatrix* cm_; - int currRow_; - int currIndex_; -}; - -/** -* \class CompressedMatrixConstIterator -* \brief Sparse matrix const_iterator. -* -* This class creates a sparse matrix const_iterator. -*/ -template -class CompressedMatrixConstIterator -{ -public: - /** - * \brief Construct a CompressedMatrix const_iterator - * - * \param cm a CompressedMatrix pointer - * \param currRow current row - * \param currIndex current index of colIndices and values - */ - CompressedMatrixConstIterator(const CompressedMatrix* cm, int currRow = 0, int currIndex = 0); - - /** - * \brief Construct a CompressedMatrixConstIterator from a CompressedMatrixIterator - * - * \param iter a CompressedMatrixIterator - */ - CompressedMatrixConstIterator(const CompressedMatrixIterator& iter); - - /** - * \brief Increment the iterator by one - */ - void operator++(); - - /** - * \brief Dereference the iterator - * - * \return a double value - */ - T operator*() const; - - /** - * \brief Get the row of the iterator - * - * \return row number - */ - int row() const; - - /** - * \brief Get the column of the iterator - * - * \return column number - */ - int col() const; - - /** - * \brief Compare whether two iterators are the same - * - * \return true if two iterators are the same; false otherwise - */ - bool operator!=(const CompressedMatrixConstIterator& iter) const; - -private: - - const CompressedMatrix* cm_; - int currRow_; - int currIndex_; -}; - - -//base exception class -class CompressedMatrixException -{ -public: - CompressedMatrixException(const std::string& m = "compressed matrix exception") : message(m) {} - const std::string& what() const {return message;} -private: - std::string message; -}; - -} //namespace compressed_matrix - - -// definition - -template -compressed_matrix::CompressedMatrix::CompressedMatrix(int numRows, int numCols) : numRows_(numRows), numCols_(numCols) -{ - rowOffsets_.resize(numRows_ + 1, 0); -} - -template -compressed_matrix::CompressedMatrix::CompressedMatrix(int numRows, int numCols, const std::map, T>& m) : numRows_(numRows), numCols_(numCols) -{ - int offset = 0; - int currRow = 0; - rowOffsets_.resize(numRows_ + 1); - for (typename std::map, T>::const_iterator it = m.begin(), end = m.end(); it != end; ++it) - { - if (it->first.first >= numRows_ || it->first.second >= numCols_) - throw CompressedMatrixException("map contains invalid index for compressed matrix constructor."); - - if (currRow <= it->first.first) - { - for (int i = currRow; i <= it->first.first; i++) - rowOffsets_[i] = offset; - currRow = it->first.first + 1; - } - - colIndices_.push_back(it->first.second); - values_.push_back(it->second); - ++offset; - } - for (int i = currRow; i < rowOffsets_.size(); i++) - rowOffsets_[i] = offset; -} - -template -compressed_matrix::CompressedMatrix::CompressedMatrix(int numRows, int numCols, std::vector rowOffsets, std::vector colIndices, std::vector values) : numRows_(numRows), numCols_(numCols) -{ - if (rowOffsets.size() != numRows_ + 1) - throw CompressedMatrixException("row offset vector's size is incorrect."); - if (colIndices.size() != values.size()) - throw CompressedMatrixException("column indices vector's size and values vector's size are different."); - - rowOffsets_ = rowOffsets; - colIndices_ = colIndices; - values_ = values; -} - -template -int compressed_matrix::CompressedMatrix::numRows() const -{ - return numRows_; -} - -template -int compressed_matrix::CompressedMatrix::numCols() const -{ - return numCols_; -} - -template -int compressed_matrix::CompressedMatrix::nnz() const -{ - return static_cast(values_.size()); -} - -template -T compressed_matrix::CompressedMatrix::get(int row, int col) const -{ - if (row >= numRows_ || col >= numCols_) - throw CompressedMatrixException("row or col contains invalid index for compressed matrix get."); - - int start = rowOffsets_[row]; - int end = rowOffsets_[row + 1]; - - if (start == end) - return 0; - - int low = start; - int high = end - 1; - int mid; - while (low <= high) - { - mid = low + (high - low) / 2; - if (colIndices_[mid] == col) - return values_[mid]; - else if (colIndices_[mid] < col) - low = mid + 1; - else - high = mid - 1; - } - - return 0; -} - - -template -T& compressed_matrix::CompressedMatrix::operator()(int row, int col) -{ - if (row >= numRows_ || col >= numCols_) - throw CompressedMatrixException("row or col contains invalid index for compressed matrix operator()."); - - int start = rowOffsets_[row]; - int end = rowOffsets_[row + 1]; - - if (start == end) //not found - { - for (int i = row + 1; i < rowOffsets_.size(); i++) - ++rowOffsets_[i]; - colIndices_.insert(colIndices_.begin() + start, col); - values_.insert(values_.begin() + start, 0); - return values_[start]; - } - - int low = start; - int high = end - 1; - int mid; - while (low <= high) - { - mid = low + (high - low) / 2; - if (colIndices_[mid] == col) - return values_[mid]; //found, return the reference - else if (colIndices_[mid] < col) - low = mid + 1; - else - high = mid - 1; - } - - //not found - for (int i = row + 1; i < rowOffsets_.size(); i++) - ++rowOffsets_[i]; - if (colIndices_[mid] < col) - ++mid; - colIndices_.insert(colIndices_.begin() + mid, col); - values_.insert(values_.begin() + mid, 0); - return values_[mid]; -} - -template -void compressed_matrix::CompressedMatrix::reserve(int minValues) -{ - colIndices_.reserve(minValues); - values_.reserve(minValues); -} - -template -typename compressed_matrix::CompressedMatrix::const_iterator compressed_matrix::CompressedMatrix::begin() const -{ - //find the last 0 in rowOffsets_[]; there will always be at least one 0 in rowOffsets_[], so doesn't need to check if rowOffsets_[low] is 0 or not - int low = 0; - int high = numRows_; - int mid; - while (low < high) - { - mid = low + (high - low + 1) / 2; - if (rowOffsets_[mid] == 0) - low = mid; - else - high = mid - 1; - } - - return CompressedMatrixConstIterator(this, low, 0); -} - -template -typename compressed_matrix::CompressedMatrix::const_iterator compressed_matrix::CompressedMatrix::end() const -{ - return CompressedMatrixConstIterator(this, numRows_); -} - -template -typename compressed_matrix::CompressedMatrix::iterator compressed_matrix::CompressedMatrix::begin() -{ - //find the last 0 in rowOffsets_[]; there will always be at least one 0 in rowOffsets_[], so doesn't need to check if rowOffsets_[low] is 0 or not - int low = 0; - int high = numRows_; - int mid; - while (low < high) - { - mid = low + (high - low + 1) / 2; - if (rowOffsets_[mid] == 0) - low = mid; - else - high = mid - 1; - } - - return CompressedMatrixIterator(this, low, 0); -} - -template -typename compressed_matrix::CompressedMatrix::iterator compressed_matrix::CompressedMatrix::end() -{ - return CompressedMatrixIterator(this, numRows_); -} - - -template -std::map, T> compressed_matrix::compressedMatrixToMap(const CompressedMatrix& cm) -{ - std::map, T> ret; - for (int i = 0; i < cm.numRows(); i++) - { - int start = cm.rowOffsets()[i]; - int end = cm.rowOffsets()[i + 1]; - for (int j = start; j < end; j++) - ret.insert(std::make_pair(std::make_pair(i, cm.colIndices()[j]), cm.values()[j])); - } - - return ret; -} - - -template -compressed_matrix::CompressedMatrixIterator::CompressedMatrixIterator(CompressedMatrix* cm, int currRow, int currIndex) : cm_(cm), currRow_(currRow), currIndex_(currIndex) -{ - -} - -template -void compressed_matrix::CompressedMatrixIterator::operator++() -{ - ++currIndex_; - if (currIndex_ == cm_->rowOffsets_[currRow_ + 1]) - { - //++currRow_; - //while (currRow_numRows_ && cm_->rowOffsets_[currRow_]==cm_->rowOffsets_[currRow_+1]) - // ++currRow_; - //find the last element which is equal to cm_->rowOffsets_[currRow_+1] - int target = cm_->rowOffsets_[currRow_ + 1]; - int low = currRow_ + 1; - int high = cm_->numRows_; - int mid; - while (low < high) - { - mid = low + (high - low + 1) / 2; - if (cm_->rowOffsets_[mid] == target) - low = mid; - else - high = mid - 1; - } - - currRow_ = low; - } -} - -template -T& compressed_matrix::CompressedMatrixIterator::operator*() const -{ - return cm_->values_[currIndex_]; -} - -template -int compressed_matrix::CompressedMatrixIterator::row() const -{ - return currRow_; -} - -template -int compressed_matrix::CompressedMatrixIterator::col() const -{ - return cm_->colIndices_[currIndex_]; -} - -template -bool compressed_matrix::CompressedMatrixIterator::operator!=(const CompressedMatrixIterator& iter) const -{ - if (cm_ != iter.cm_) - throw CompressedMatrixException("compare iterators from two different CompressedMatrix objects."); - - return (currRow_ != iter.currRow_); -} - -template -compressed_matrix::CompressedMatrixConstIterator::CompressedMatrixConstIterator(const CompressedMatrix* cm, int currRow, int currIndex) : cm_(cm), currRow_(currRow), currIndex_(currIndex) -{ - -} - -template -compressed_matrix::CompressedMatrixConstIterator::CompressedMatrixConstIterator(const compressed_matrix::CompressedMatrixIterator& iter) : cm_(iter.cm_), currRow_(iter.currRow_), currIndex_(iter.currIndex_) -{ - -} - -template -void compressed_matrix::CompressedMatrixConstIterator::operator++() -{ - ++currIndex_; - if (currIndex_ == cm_->rowOffsets_[currRow_ + 1]) - { - //++currRow_; - //while (currRow_numRows_ && cm_->rowOffsets_[currRow_]==cm_->rowOffsets_[currRow_+1]) - // ++currRow_; - //find the last element which is equal to cm_->rowOffsets_[currRow_+1] - int target = cm_->rowOffsets_[currRow_ + 1]; - int low = currRow_ + 1; - int high = cm_->numRows_; - int mid; - while (low < high) - { - mid = low + (high - low + 1) / 2; - if (cm_->rowOffsets_[mid] == target) - low = mid; - else - high = mid - 1; - } - - currRow_ = low; - } -} - -template -T compressed_matrix::CompressedMatrixConstIterator::operator*() const -{ - return cm_->values_[currIndex_]; -} - -template -int compressed_matrix::CompressedMatrixConstIterator::row() const -{ - return currRow_; -} - -template -int compressed_matrix::CompressedMatrixConstIterator::col() const -{ - return cm_->colIndices_[currIndex_]; -} - -template -bool compressed_matrix::CompressedMatrixConstIterator::operator!=(const CompressedMatrixConstIterator& iter) const -{ - if (cm_ != iter.cm_) - throw CompressedMatrixException("compare iterators from two different CompressedMatrix objects."); - - return (currRow_ != iter.currRow_); -} - -#endif // COMPRESSED_MATRIX_HPP_INCLUDED diff --git a/dimod/roof_duality/src/fix_variables.cpp b/dimod/roof_duality/src/fix_variables.cpp deleted file mode 100644 index da956cfa3..000000000 --- a/dimod/roof_duality/src/fix_variables.cpp +++ /dev/null @@ -1,1297 +0,0 @@ -/** -# Copyright 2018 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. -# -# ================================================================================================ -*/ -#include "fix_variables.hpp" -#include "compressed_matrix.hpp" - -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include - -//for debugging only -#include -#include -#include -#include - -namespace -{ - -class compClass -{ -public: - bool operator()(const std::pair& a, const std::pair& b) - { - if (a.second != b.second) - return !(a.second < b.second); - else - return a.first < b.first; - } -}; - -bool compareAbs(double a, double b) -{ - return std::abs(a) < std::abs(b); -} - - -//the index is 1-based, according to the paper -struct Posiform -{ - std::vector, long long int> > quadratic; - std::vector > linear; - long long int cst; - int numVars; -}; - -struct SC -{ - std::vector original; - std::vector positive; - std::vector negative; -}; - -struct SCRet -{ - std::vector S; - compressed_matrix::CompressedMatrix G; -}; - -Posiform BQPToPosiform(const compressed_matrix::CompressedMatrix& Q, bool makeRandom = false) -{ - int numVariables = Q.numRows(); - Posiform ret; - ret.numVars = numVariables; - std::vector q(numVariables); - - //q = diag(Q); - for (int i = 0; i < numVariables; i++) - q[i] = Q.get(i, i); - - //tmpQ = triu(Q,1); - std::vector tmpQRowOffsets(numVariables + 1); - std::vector tmpQColIndices; tmpQColIndices.reserve(Q.nnz()); - std::vector tmpQValues; tmpQValues.reserve(Q.nnz()); - int index = 0; - for (int i = 0; i < Q.numRows(); i++) - { - int start = Q.rowOffsets()[i]; - int end = Q.rowOffsets()[i + 1]; - tmpQRowOffsets[i] = index; - for (int j = start; j < end; j++) - { - int r = i; - int c = Q.colIndices()[j]; - if (r != c) - { - tmpQColIndices.push_back(c); - tmpQValues.push_back(Q.values()[j]); //Q.values()[j] is Q(r, c) - ++index; - } - } - } - tmpQRowOffsets.back() = index; //number of non-zero elements - - compressed_matrix::CompressedMatrix tmpQ(numVariables, numVariables, tmpQRowOffsets, tmpQColIndices, tmpQValues); - - - if (!makeRandom) //make a deterministic posiform from Q - { - //cPrime = q + sum(tmpQ.*(tmpQ<0),2); and quadratic term - std::vector cPrime = q; - for (int i = 0; i < tmpQ.numRows(); i++) - { - int start = tmpQ.rowOffsets()[i]; - int end = tmpQ.rowOffsets()[i + 1]; - for (int j = start; j < end; j++) - { - int r = i; - int c = tmpQ.colIndices()[j]; - if (tmpQ.values()[j] < 0) - { - cPrime[i] += tmpQ.values()[j]; - ret.quadratic.push_back(std::make_pair(std::make_pair(r + 1, -(c + 1)), -tmpQ.values()[j])); //+1 makes it 1-based - } - else if (tmpQ.values()[j] > 0) - ret.quadratic.push_back(std::make_pair(std::make_pair(r + 1, c + 1), tmpQ.values()[j])); //+1 makes it 1-based - } - } - - //constant term and linear term - ret.cst = 0; - for (int i = 0; i < numVariables; i++) - { - if (cPrime[i] < 0) - { - ret.cst += cPrime[i]; //constant term - ret.linear.push_back(std::make_pair(-(i + 1), -cPrime[i])); //+1 makes it 1-based - } - else if (cPrime[i] > 0) - ret.linear.push_back(std::make_pair(i + 1, cPrime[i])); //+1 makes it 1-based - } - } - else - { - std::vector linear = q; - std::set > negPairs; - for (int i = 0; i < tmpQ.numRows(); i++) - { - int start = tmpQ.rowOffsets()[i]; - int end = tmpQ.rowOffsets()[i + 1]; - for (int j = start; j < end; j++) - { - int r = i; - int c = tmpQ.colIndices()[j]; - if (tmpQ.values()[j] < 0) - { - if (rand() % 2 == 0) - negPairs.insert(std::make_pair(-(r + 1), c + 1)); //+1 makes it 1-based - else - negPairs.insert(std::make_pair(r + 1, -(c + 1))); //+1 makes it 1-based - } - } - } - - for (int i = 0; i < tmpQ.numRows(); i++) - { - int start = tmpQ.rowOffsets()[i]; - int end = tmpQ.rowOffsets()[i + 1]; - for (int j = start; j < end; j++) - { - int r = i; - int c = tmpQ.colIndices()[j]; - if (tmpQ.values()[j] > 0) - ret.quadratic.push_back(std::make_pair(std::make_pair(r + 1, c + 1), tmpQ.values()[j])); - else - { - if (negPairs.find(std::make_pair(r + 1, -(c + 1)))!=negPairs.end()) - ret.quadratic.push_back(std::make_pair(std::make_pair(r + 1, -(c + 1)), -tmpQ.values()[j])); - else - ret.quadratic.push_back(std::make_pair(std::make_pair(-(r + 1), c + 1), -tmpQ.values()[j])); - } - } - } - - for (std::set >::iterator it = negPairs.begin(); it != negPairs.end(); ++it) - { - if (it->first > 0) - linear[it->first - 1] += tmpQ(it->first - 1, -it->second - 1); - else - linear[it->second - 1] += tmpQ(-it->first - 1, it->second - 1); - } - - ret.cst = 0; - - for (int i = 0; i < numVariables; i++) - { - if (linear[i] < 0) - { - ret.cst += linear[i]; - ret.linear.push_back(std::make_pair(-(i + 1), -linear[i])); - } - else if (linear[i] > 0) - ret.linear.push_back(std::make_pair(i + 1, linear[i])); - } - } - - return ret; -} - -compressed_matrix::CompressedMatrix posiformToImplicationNetwork_1(const Posiform& p) -{ - int n = p.numVars; - int numVertices = 2 * n + 2; - int source = 0; - int sink = n + 1; - - //n = p.numVars - //0: source; - //1 to n: x_1 to x_n - //n+1: sink - //n+2 to 2*n+1: \overline_x_1 to \overline_x_n - std::map, long long int> m; - for (int i = 0; i < p.linear.size(); i++) - { - int v = p.linear[i].first; - - long long int capacity = p.linear[i].second; // originally p.linear[i].second/2 - if (v > 0) - { - m[std::make_pair(source, v + n + 1)] = capacity; - m[std::make_pair(v, sink)] = capacity; - } - else - { - v = std::abs(v); - m[std::make_pair(source, v)] = capacity; - m[std::make_pair(v + n + 1, sink)] = capacity; - } - } - - for (int i = 0; i < p.quadratic.size(); i++) - { - int v_1 = p.quadratic[i].first.first; - int v_2 = p.quadratic[i].first.second; - - long long int capacity = p.quadratic[i].second; // originally p.quadratic[i].second/2 - if (v_1 < 0) - { - v_1 = std::abs(v_1); - m[std::make_pair(v_1 + n + 1, v_2 + n + 1)] = capacity; - m[std::make_pair(v_2, v_1)] = capacity; - } - else if (v_2 < 0) - { - v_2 = std::abs(v_2); - m[std::make_pair(v_1, v_2)] = capacity; - m[std::make_pair(v_2 + n + 1, v_1 + n + 1)] = capacity; - } - else - { - m[std::make_pair(v_1, v_2 + n + 1)] = capacity; - m[std::make_pair(v_2, v_1 + n + 1)] = capacity; - } - } - - return compressed_matrix::CompressedMatrix(numVertices, numVertices, m); -} - -compressed_matrix::CompressedMatrix posiformToImplicationNetwork_2(const Posiform& p) -{ - int n = p.numVars; - int numVertices = 2 * n + 2; - int source = n; - int sink = 2 * n + 1; - - //n = p.numVars - //0 to n-1: x_1 to x_n - //n: source; - //n+1 to 2*n: \overline_x_1 to \overline_x_n - //2*n+1: sink - std::map, long long int> m; - for (int i = 0; i < p.linear.size(); i++) - { - int v = p.linear[i].first; - long long int capacity = p.linear[i].second; // originally p.linear[i].second/2 - if (v > 0) - { - m[std::make_pair(source, v + n)] = capacity; - m[std::make_pair(v - 1, sink)] = capacity; - } - else - { - v = std::abs(v); - m[std::make_pair(source, v - 1)] = capacity; - m[std::make_pair(v + n, sink)] = capacity; - } - } - - for (int i = 0; i < p.quadratic.size(); i++) - { - int v_1 = p.quadratic[i].first.first; - int v_2 = p.quadratic[i].first.second; - long long int capacity = p.quadratic[i].second; // originally p.quadratic[i].second/2 - if (v_1 < 0) - { - v_1 = std::abs(v_1); - m[std::make_pair(v_1 + n, v_2 + n)] = capacity; - m[std::make_pair(v_2 - 1, v_1 - 1)] = capacity; - } - else if (v_2 < 0) - { - v_2 = std::abs(v_2); - m[std::make_pair(v_1 - 1, v_2 - 1)] = capacity; - m[std::make_pair(v_2 + n, v_1 + n)] = capacity; - } - else - { - m[std::make_pair(v_1 - 1, v_2 + n)] = capacity; - m[std::make_pair(v_2 - 1, v_1 + n)] = capacity; - } - } - - return compressed_matrix::CompressedMatrix(numVertices, numVertices, m); -} - -//this version returns R instead of F -compressed_matrix::CompressedMatrix maxFlow(const compressed_matrix::CompressedMatrix& A) -{ - int numVertices = A.numRows(); - int numVariables = numVertices / 2 - 1; - - //clock_t curr_1 = clock(); - //clock_t curr_2; - - using namespace boost; - - typedef adjacency_list_traits Traits; - typedef adjacency_list, property > > > Graph; //for edge capacity is long long int - //typedef adjacency_list, property > > > Graph; //for edge capacity is double - - Graph g; - - property_map::type capacity = get(edge_capacity, g); - property_map::type reverse_edge = get(edge_reverse, g); - property_map::type residual_capacity = get(edge_residual_capacity, g); - - std::vector verts(numVertices); - for (int i = 0; i < numVertices; ++i) - verts[i] = add_vertex(g); - - Traits::vertex_descriptor s = verts[0]; - Traits::vertex_descriptor t = verts[numVariables + 1]; - - for (int i = 0; i < A.numRows(); i++) - { - int start = A.rowOffsets()[i]; - int end = A.rowOffsets()[i + 1]; - for (int j = start; j < end; j++) - { - int r = i; - int c = A.colIndices()[j]; - long long int cap = A.values()[j]; - - Traits::edge_descriptor e1, e2; - bool in1, in2; - boost::tie(e1, in1) = add_edge(verts[r], verts[c], g); - boost::tie(e2, in2) = add_edge(verts[c], verts[r], g); - - capacity[e1] = cap; - capacity[e2] = 0; - reverse_edge[e1] = e2; - reverse_edge[e2] = e1; - } - } - - - //curr_2 = clock(); - //mexPrintf("inside maxFlow int version: Time elapsed_constructing_graph_for_max_flow: %f\n", ((double)curr_2 - curr_1) / CLOCKS_PER_SEC); - //curr_1 = curr_2; - - push_relabel_max_flow(g, s, t); - - //curr_2 = clock(); - //mexPrintf("inside maxFlow int version: Time elapsed_for_boost_max_flow: %f\n", ((double)curr_2 - curr_1) / CLOCKS_PER_SEC); - //curr_1 = curr_2; - - std::vector RRowOffsets(numVertices+1); - std::vector RColIndices; - std::vector RValues; - int offset = 0; - int currRow = 0; - graph_traits::vertex_iterator u_iter, u_end; - graph_traits::out_edge_iterator ei, e_end; - for (boost::tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) - for (boost::tie(ei, e_end) = out_edges(*u_iter, g); ei != e_end; ++ei) - { - if (capacity[*ei] > 0) - { - if (currRow <= *u_iter) - { - for (int i = currRow; i <= *u_iter; i++) - RRowOffsets[i] = offset; - currRow = static_cast((*u_iter) + 1); - } - - RColIndices.push_back(static_cast(target(*ei, g))); - RValues.push_back(residual_capacity[*ei]); - ++offset; - } - } - for (int i = currRow; i < RRowOffsets.size(); i++) - RRowOffsets[i] = offset; - - compressed_matrix::CompressedMatrix R(numVertices, numVertices, RRowOffsets, RColIndices, RValues); //residual - - //curr_2 = clock(); - //mexPrintf("inside maxFlow int version: Time elapsed_get_R: %f\n", ((double)curr_2 - curr_1) / CLOCKS_PER_SEC); - //curr_1 = curr_2; - - std::map, long long int> rm = compressed_matrix::compressedMatrixToMap(R); - std::map, long long int> rmMissing; - - //curr_2 = clock(); - //mexPrintf("inside maxFlow int version: Time elapsed_copy_R_to_map: %f\n", ((double)curr_2 - curr_1) / CLOCKS_PER_SEC); - //curr_1 = curr_2; - - for (int i = 0; i < R.numRows(); i++) - { - int start = R.rowOffsets()[i]; - int end = R.rowOffsets()[i + 1]; - for (int j = start; j < end; j++) - { - int r = i; - int c = R.colIndices()[j]; - long long int Arc = A.get(r, c); - long long int Rcr = R.get(c, r); - if (Arc != 0 && R.values()[j] + Rcr - Arc != 0) //it->second is: R(r, c); here it means R(r, c)+R(c, r)!=A(r, c), so R(c, r) is missing - rmMissing.insert(std::make_pair(std::make_pair(c, r), Arc - R.values()[j])); - } - } - - for (std::map, long long int>::iterator it = rmMissing.begin(), end = rmMissing.end(); it != end; ++it) - rm.insert(*it); - - //curr_2 = clock(); - //mexPrintf("inside maxFlow int version: Time elapsed_get_missing_elements_in_R: %f\n", ((double)curr_2 - curr_1) / CLOCKS_PER_SEC); - //curr_1 = curr_2; - - return compressed_matrix::CompressedMatrix(numVertices, numVertices, rm); -} - -std::vector bfs_for_method_2(const compressed_matrix::CompressedMatrix& g, int s) -{ - std::queue q; - q.push(s); - int numVertices = g.numRows(); - std::vector visited(numVertices, -1); - visited[s] = 1; - - while (!q.empty()) - { - int curr = q.front(); q.pop(); - int start = g.rowOffsets()[curr]; - int end = g.rowOffsets()[curr + 1]; - for (int j = start; j < end; j++) - { - int c = g.colIndices()[j]; - if (g.values()[j] != 0 && visited[c] == -1) - { - visited[c] = 1; - q.push(c); - } - } - } - - return visited; -} - -compressed_matrix::CompressedMatrix makeResidualSymmetric(const compressed_matrix::CompressedMatrix& R) -{ - //clock_t curr_1 = clock(); - //clock_t curr_2; - - int numVertices = R.numRows(); - int numVariables = numVertices / 2 - 1; - std::map, long long int> MsymR; - for (int i = 0; i < R.numRows(); i++) - { - int start = R.rowOffsets()[i]; - int end = R.rowOffsets()[i + 1]; - for (int j = start; j < end; j++) - { - int r = i; - int c = R.colIndices()[j]; - if (r != c) - { - MsymR[std::make_pair(r, c)] += R.values()[j]; - int compR; - if (r <= numVariables) - compR = r + (numVariables + 1); - else - compR = r - (numVariables + 1); - int compC; - if (c <= numVariables) - compC = c + (numVariables + 1); - else - compC = c - (numVariables + 1); - MsymR[std::make_pair(compC, compR)] += R.values()[j]; - } - } - } - - //curr_2 = clock(); - //mexPrintf("inside makeResidualSym: Time elapsed_get_MsymR: %f\n", ((double)curr_2 - curr_1) / CLOCKS_PER_SEC); - //curr_1 = curr_2; - - //remove extra 0s in MsymR, required because there will possibly be 0 in MsymR and the 0 will add non-exist edges in stronglyConnectedComponents() - //cause the connected components results incorrct !!! - std::map, long long int> MsymRWithoutZero; - for (std::map, long long int>::const_iterator it = MsymR.begin(), end = MsymR.end(); it != end; ++it) - { - if (it->second != 0 && it->first.first != numVariables + 1 && it->first.second != 0) //add clearing R here !!! - MsymRWithoutZero.insert(*it); - } - - //curr_2 = clock(); - //mexPrintf("inside makeResidualSym: Time elapsed_remove_zero_in_MsymR: %f\n", ((double)curr_2 - curr_1) / CLOCKS_PER_SEC); - //curr_1 = curr_2; - - return compressed_matrix::CompressedMatrix(numVertices, numVertices, MsymRWithoutZero); -} - -SCRet stronglyConnectedComponents(const compressed_matrix::CompressedMatrix& R) -{ - //clock_t curr_1 = clock(); - //clock_t curr_2; - - int numVertices = R.numRows(); - int numVariables = numVertices / 2 - 1; - - using namespace boost; - - typedef adjacency_list_traits < vecS, vecS, directedS > Traits; - typedef adjacency_list < vecS, vecS, directedS, - property < vertex_name_t, std::string, - property < vertex_index_t, long, - property < vertex_color_t, boost::default_color_type, - property < vertex_distance_t, long, - property < vertex_predecessor_t, Traits::edge_descriptor > > > > > - > Graph; - - typedef graph_traits::vertex_descriptor Vertex; - - Graph G; - - std::vector root(numVertices); - for (int i = 0; i < numVertices; ++i) - root[i] = add_vertex(G); - - std::vector component(num_vertices(G)), discover_time(num_vertices(G)); - std::vector color(num_vertices(G)); - - for (int i = 0; i < R.numRows(); i++) - { - int start = R.rowOffsets()[i]; - int end = R.rowOffsets()[i + 1]; - for (int j = start; j < end; j++) - { - int r = i; - int c = R.colIndices()[j]; - add_edge(root[r], root[c], G); - } - } - - //curr_2 = clock(); - //mexPrintf("Time elapsed_constructing_graph_for_strongly_connected_components_graph: %f\n", ((double)curr_2 - curr_1) / CLOCKS_PER_SEC); - //curr_1 = curr_2; - - //int num = strong_components(G, &component[0], root_map(&root[0]).color_map(&color[0]).discover_time_map(&discover_time[0])); - int num = strong_components(G, boost::make_iterator_property_map(component.begin(), boost::get(boost::vertex_index, G)), - root_map(boost::make_iterator_property_map(root.begin(), boost::get(boost::vertex_index, G))).color_map(boost::make_iterator_property_map(color.begin(), boost::get(boost::vertex_index, G))).discover_time_map(boost::make_iterator_property_map(discover_time.begin(), boost::get(boost::vertex_index, G)))); - - //curr_2 = clock(); - //mexPrintf("Time elapsed_boost_strongly_connected_components_graph: %f\n", ((double)curr_2 - curr_1) / CLOCKS_PER_SEC); - //curr_1 = curr_2; - - //here original's content is from 0 to 2*n+1, 0: source, 1 to n, variables, n+1: sink, n+2 to 2*n+1: bar variables - //positive: from 1 to n (1-based variables index) - //negative: from 1 to n (1-based variables index) - std::vector vsc(num); - for (int i = 0; i < component.size(); i++) - { - int whichComponent = component[i]; - vsc[whichComponent].original.push_back(i); - if (vsc[whichComponent].original.back() <= numVariables) - vsc[whichComponent].positive.push_back(vsc[whichComponent].original.back()); - else - vsc[whichComponent].negative.push_back(vsc[whichComponent].original.back() - (numVariables + 1)); - } - - std::map, long long int> M; - - //speed up the program a lot !!! - //build the graph G for strongly connected components - for (int i = 0; i < R.numRows(); i++) - { - int start = R.rowOffsets()[i]; - int end = R.rowOffsets()[i + 1]; - for (int j = start; j < end; j++) - { - int r = i; - int c = R.colIndices()[j]; - int sc1 = component[r]; - int sc2 = component[c]; - if (sc1 != sc2 && R.values()[j] > 0 && M.find(std::make_pair(sc1, sc2)) == M.end()) - M.insert(std::make_pair(std::make_pair(sc1, sc2), 1)); - } - } - - //curr_2 = clock(); - //mexPrintf("Time elapsed_after_boost_scc: %f\n", ((double)curr_2 - curr_1) / CLOCKS_PER_SEC); - //curr_1 = curr_2; - - SCRet ret; - ret.S = vsc; - ret.G = compressed_matrix::CompressedMatrix(num, num, M); - - return ret; -} - -std::vector classifyStronglyConnectedComponents(const std::vector& S) -{ - std::vector ret(S.size()); //0: self-complement, 1: non-self-complement - - for (int i = 0; i < S.size(); i++) - { - bool flag = true; - std::set pos(S[i].positive.begin(), S[i].positive.end()); - for (int j = 0; j < S[i].negative.size(); j++) - { - if (pos.find(S[i].negative[j]) != pos.end()) - { - flag = false; - break; - } - } - - if (flag) - ret[i] = 1; - else - ret[i] = 0; - } - - return ret; -} - -void shortestPath(int src, const compressed_matrix::CompressedMatrix& g, std::vector& visited) -{ - std::queue q; - q.push(src); - visited.resize(g.numRows(), 0); - visited[src] = 1; - - while (!q.empty()) - { - int curr = q.front(); q.pop(); - - int start = g.rowOffsets()[curr]; - int end = g.rowOffsets()[curr + 1]; - for (int j = start; j < end; j++) - { - int c = g.colIndices()[j]; - if (g.values()[j] != 0 && !visited[c]) - { - visited[c] = 2; - q.push(c); - } - } - } -} - -std::vector > fixVariables(const std::vector& classifiedSC, const SCRet& scRet, int numVariables) -{ - std::vector > ret; - - int x0Location; - int x0BarLocation; - std::vector I; - - for (int i = 0; i < classifiedSC.size(); i++) - { - if (classifiedSC[i] == 1) //1: non self-complement - { - if (scRet.S[i].original[0] == 0) - x0Location = i; - else if (scRet.S[i].original[0] == numVariables + 1) - x0BarLocation = i; - else - I.push_back(i); - } - } - - //calculate if there is a path from x0 to other nodes - std::vector visited; - shortestPath(x0Location, scRet.G, visited); //speed up the program - - //find complement pairs - std::vector positive(numVariables + 1, -1); - std::vector negative(numVariables + 1, -1); - for (int i = 0; i < I.size(); i++) - { - for (int k = 0; k < scRet.S[I[i]].positive.size(); k++) - positive[scRet.S[I[i]].positive[k]] = I[i]; - for (int k = 0; k < scRet.S[I[i]].negative.size(); k++) - negative[scRet.S[I[i]].negative[k]] = I[i]; - } - - std::map complementPairs; - for (int i = 0; i < positive.size(); i++) - { - if (positive[i] != -1 && negative[i] != -1) - { - complementPairs[positive[i]] = negative[i]; - complementPairs[negative[i]] = positive[i]; - } - } - - std::queue q; - - std::map, long long int> GTransMap; - for (int i = 0; i < scRet.G.numRows(); i++) - { - int start = scRet.G.rowOffsets()[i]; - int end = scRet.G.rowOffsets()[i + 1]; - for (int j = start; j < end; j++) - GTransMap.insert(std::make_pair(std::make_pair(scRet.G.colIndices()[j], i), scRet.G.values()[j])); - } - compressed_matrix::CompressedMatrix GTrans(scRet.G.numRows(), scRet.G.numRows(), GTransMap); - - std::vector outDegrees(scRet.G.numRows(), -1); - for (int i = 0; i < classifiedSC.size(); i++) - { - if (classifiedSC[i] == 1) //non self-complement components - { - outDegrees[i] = scRet.G.rowOffsets()[i + 1] - scRet.G.rowOffsets()[i]; - if (outDegrees[i] == 0 && i != x0Location && i != x0BarLocation) //need to push the original outdegree 0 nodes into q !!! - q.push(i); - } - } - outDegrees[x0Location] = -1; - outDegrees[x0BarLocation] = -1; - - for (int i = 0; i < visited.size(); i++) - { - if (visited[i] == 2) //exclude x0 - { - for (int k = 0; k < scRet.S[i].positive.size(); k++) - ret.push_back(std::make_pair(scRet.S[i].positive[k], 1)); - - for (int k = 0; k < scRet.S[i].negative.size(); k++) - ret.push_back(std::make_pair(scRet.S[i].negative[k], 0)); - - int complement = complementPairs[i]; - - outDegrees[i] = -1; - outDegrees[complement] = -1; - } - } - - - //decrease the outdegrees of node which has outgoing edges to i and to complement - //push node which has 0 outdegrees into the queue - for (int i = 0; i < visited.size(); i++) - { - if (visited[i] == 2) //exclude x0 - { - int start = GTrans.rowOffsets()[i]; - int end = GTrans.rowOffsets()[i + 1]; - for (int k = start; k < end; k++) - { - if (outDegrees[GTrans.colIndices()[k]] > 0) - { - --outDegrees[GTrans.colIndices()[k]]; - if (outDegrees[GTrans.colIndices()[k]] == 0) - q.push(GTrans.colIndices()[k]); - } - } - - int complement = complementPairs[i]; - start = GTrans.rowOffsets()[complement]; - end = GTrans.rowOffsets()[complement + 1]; - for (int k = start; k < end; k++) - { - if (outDegrees[GTrans.colIndices()[k]] > 0) - { - --outDegrees[GTrans.colIndices()[k]]; - if (outDegrees[GTrans.colIndices()[k]] == 0) - q.push(GTrans.colIndices()[k]); - } - } - } - } - - - while (!q.empty()) - { - int curr = q.front(); q.pop(); - if (outDegrees[curr] == 0) - { - outDegrees[curr] = -1; - int complement = complementPairs[curr]; - outDegrees[complement] = -1; - - //fixed all variables in component curr - for (int k = 0; k < scRet.S[curr].positive.size(); k++) - ret.push_back(std::make_pair(scRet.S[curr].positive[k], 1)); - - for (int k = 0; k < scRet.S[curr].negative.size(); k++) - ret.push_back(std::make_pair(scRet.S[curr].negative[k], 0)); - - //decrease the outdegrees of node which has outgoing edges to i and to complement - //push node which has 0 outdegrees into the queue - int start = GTrans.rowOffsets()[curr]; - int end = GTrans.rowOffsets()[curr + 1]; - for (int k = start; k < end; k++) - { - if (outDegrees[GTrans.colIndices()[k]] > 0) - { - --outDegrees[GTrans.colIndices()[k]]; - if (outDegrees[GTrans.colIndices()[k]] == 0) - q.push(GTrans.colIndices()[k]); - } - } - - start = GTrans.rowOffsets()[complement]; - end = GTrans.rowOffsets()[complement + 1]; - for (int k = start; k < end; k++) - { - if (outDegrees[GTrans.colIndices()[k]] > 0) - { - --outDegrees[GTrans.colIndices()[k]]; - if (outDegrees[GTrans.colIndices()[k]] == 0) - q.push(GTrans.colIndices()[k]); - } - } - } - } - - std::sort(ret.begin(), ret.end(), compClass()); - - return ret; -} - -compressed_matrix::CompressedMatrix computeNewQAndOffset(const compressed_matrix::CompressedMatrix& Q, const std::vector >& fixed, double& offset) -{ - //fixed now is 1-based - if (!fixed.empty()) - { - int numVariables = Q.numRows(); - - std::map mI; - for (int i = 0; i < fixed.size(); i++) - mI[fixed[i].first - 1] = i; //-1 to make it 0-based - - std::vector J; - std::map mJ; - int cnt = 0; - for (int i = 0; i < numVariables; i++) - { - if (mI.find(i) == mI.end()) - { - J.push_back(i); - mJ[i] = cnt++; - } - } - - std::map, double> QII; - std::map, double> QJJ; - std::map, double> QIJ; - std::map, double> QJI; - - for (int i = 0; i < Q.numRows(); i++) - { - int start = Q.rowOffsets()[i]; - int end = Q.rowOffsets()[i + 1]; - for (int j = start; j < end; j++) - { - int r = i; - int c = Q.colIndices()[j]; - if (mI.find(r) != mI.end() && mI.find(c) != mI.end()) - QII[std::make_pair(mI[r], mI[c])] = Q.values()[j]; - if (mJ.find(r) != mJ.end() && mJ.find(c) != mJ.end()) - QJJ[std::make_pair(mJ[r], mJ[c])] = Q.values()[j]; - if (mI.find(r) != mI.end() && mJ.find(c) != mJ.end()) - QIJ[std::make_pair(mI[r], mJ[c])] = Q.values()[j]; - if (mJ.find(r) != mJ.end() && mI.find(c) != mI.end()) - QJI[std::make_pair(mJ[r], mI[c])] = Q.values()[j]; - } - } - - //off_set = x0'*Q_II*x0; - offset = 0; - for (std::map, double>::const_iterator it = QII.begin(), end = QII.end(); it != end; ++it) - { - int r = it->first.first; - int c = it->first.second; - offset += fixed[r].second * it->second * fixed[c].second; - } - - //x0'*Q_IJ - std::vector tmp2(numVariables - fixed.size(), 0); - for (std::map, double>::const_iterator it = QIJ.begin(), end = QIJ.end(); it != end; ++it) - { - int r = it->first.first; - int c = it->first.second; - tmp2[c] += it->second * fixed[r].second; - } - - //Q_JI*x0 - std::vector tmp3(numVariables - fixed.size(), 0); - for (std::map, double>::const_iterator it = QJI.begin(), end = QJI.end(); it != end; ++it) - { - int r = it->first.first; - int c = it->first.second; - tmp3[r] += it->second * fixed[c].second; - } - - for (int i = 0; i < (int)tmp3.size(); i++) - tmp2[i] += tmp3[i]; - - cnt = 0; - for (int i = 0; i < (int)mJ.size() * (int)mJ.size(); i += numVariables - static_cast(fixed.size()) + 1) - { - int r = i / static_cast(mJ.size()); - int c = i % mJ.size(); - QJJ[std::make_pair(r, c)] += tmp2[cnt++]; - } - - std::map, double> mfixedQ; - for (std::map, double>::const_iterator it = QJJ.begin(), end = QJJ.end(); it != end; ++it) - { - int r = it->first.first; - int c = it->first.second; - if (it->second != 0) - mfixedQ[std::make_pair(J[r], J[c])] = it->second; - } - - return compressed_matrix::CompressedMatrix(numVariables, numVariables, mfixedQ); - } - else - { - offset = 0; - return Q; - } -} - -std::vector > applyImplication(const compressed_matrix::CompressedMatrix& A) -{ - int numVertices = A.numRows(); - int numVariables = numVertices / 2 - 1; - - //debuging only - //clock_t curr_1 = clock(); - //clock_t curr_2; - - using namespace boost; - - typedef adjacency_list_traits Traits; - typedef adjacency_list, property > > > Graph; //for edge capacity is long long int - //typedef adjacency_list, property > > > Graph; //for edge capacity is double - - Graph g; - - property_map::type capacity = get(edge_capacity, g); - property_map::type reverse_edge = get(edge_reverse, g); - property_map::type residual_capacity = get(edge_residual_capacity, g); - - std::vector verts(numVertices); - for (int i = 0; i < numVertices; ++i) - verts[i] = add_vertex(g); - - Traits::vertex_descriptor s = verts[numVariables]; - Traits::vertex_descriptor t = verts[2 * numVariables + 1]; - - for (int i = 0; i < A.numRows(); i++) - { - int start = A.rowOffsets()[i]; - int end = A.rowOffsets()[i + 1]; - for (int j = start; j < end; j++) - { - int r = i; - int c = A.colIndices()[j]; - long long int cap = A.values()[j]; - - Traits::edge_descriptor e1, e2; - bool in1, in2; - boost::tie(e1, in1) = add_edge(verts[r], verts[c], g); - boost::tie(e2, in2) = add_edge(verts[c], verts[r], g); - - capacity[e1] = cap; - capacity[e2] = 0; - reverse_edge[e1] = e2; - reverse_edge[e2] = e1; - } - } - - //curr_2 = clock(); - //mexPrintf("inside applyImplication int version: Time elapsed_building_graph: %f\n", ((double)curr_2 - curr_1) / CLOCKS_PER_SEC); - //curr_1 = curr_2; - - push_relabel_max_flow(g, s, t); - - //curr_2 = clock(); - //mexPrintf("inside applyImplication int version: Time elapsed_max_flow: %f\n", ((double)curr_2 - curr_1) / CLOCKS_PER_SEC); - //curr_1 = curr_2; - - //residual - std::vector RRowOffsets(numVertices+1); - std::vector RColIndices; - std::vector RValues; - //flow - std::vector FRowOffsets(numVertices+1); - std::vector FColIndices; - std::vector FValues; - - int offset = 0; - int currRow = 0; - graph_traits::vertex_iterator u_iter, u_end; - graph_traits::out_edge_iterator ei, e_end; - for (boost::tie(u_iter, u_end) = vertices(g); u_iter != u_end; ++u_iter) - for (boost::tie(ei, e_end) = out_edges(*u_iter, g); ei != e_end; ++ei) - { - if (capacity[*ei] > 0) - { - if (currRow <= *u_iter) - { - for (int i = currRow; i <= *u_iter; i++) - { - RRowOffsets[i] = offset; - FRowOffsets[i] = offset; - } - currRow = static_cast((*u_iter) + 1); - } - RColIndices.push_back(static_cast(target(*ei, g))); - RValues.push_back(residual_capacity[*ei]); - FColIndices.push_back(static_cast(target(*ei, g))); - FValues.push_back(capacity[*ei] - residual_capacity[*ei]); - ++offset; - } - } - for (int i = currRow; i < RRowOffsets.size(); i++) - RRowOffsets[i] = offset; - for (int i = currRow; i < FRowOffsets.size(); i++) - FRowOffsets[i] = offset; - - compressed_matrix::CompressedMatrix R(numVertices, numVertices, RRowOffsets, RColIndices, RValues); //residual - compressed_matrix::CompressedMatrix F(numVertices, numVertices, FRowOffsets, FColIndices, FValues); //flow - - //curr_2 = clock(); - //mexPrintf("inside applyImplication int version: Time elapsed_get_RF: %f\n", ((double)curr_2 - curr_1) / CLOCKS_PER_SEC); - //curr_1 = curr_2; - - //resid = (A>0).*R + F'.*(A'>0); - std::map, long long int> residM; - for (int i = 0; i < A.numRows(); i++) - { - int start = A.rowOffsets()[i]; - int end = A.rowOffsets()[i + 1]; - for (int j = start; j < end; j++) - { - int r = i; - int c = A.colIndices()[j]; - long long int RValue = R.get(r, c); - if (A.values()[j] > 0 && RValue != 0 && r != 2 * numVariables + 1 && c != numVariables) - residM[std::make_pair(r, c)] += RValue; - - long long int FTValue = F.get(r, c); //not F(c, r) !!! - if (A.values()[j] > 0 && FTValue != 0 && c != 2 * numVariables + 1 && r != numVariables) - residM[std::make_pair(c, r)] += FTValue; - } - } - - compressed_matrix::CompressedMatrix resid(numVertices, numVertices, residM); - - //curr_2 = clock(); - //mexPrintf("inside applyImplication int version: Time elapsed_finally_get_resid: %f\n", ((double)curr_2 - curr_1) / CLOCKS_PER_SEC); - //curr_1 = curr_2; - - std::vector forced = bfs_for_method_2(resid, numVariables); - - //curr_2 = clock(); - //mexPrintf("inside applyImplication int version: Time elapsed_bfs_for_method_2: %f\n", ((double)curr_2 - curr_1) / CLOCKS_PER_SEC); - //curr_1 = curr_2; - - std::vector > fixed; - for (int i = 0; i < forced.size(); i++) - { - if (forced[i] > 0) - { - if (i <= numVariables - 1) - fixed.push_back(std::make_pair(i + 1, 1)); //1-based - else if (i >= numVariables + 1 && i <= 2 * numVariables) - fixed.push_back(std::make_pair(i - numVariables, 0)); //1-based - } - } - - std::sort(fixed.begin(), fixed.end(), compClass()); - - //curr_2 = clock(); - //mexPrintf("inside applyImplication int version: Time elapsed_fix_vars_and_sorting: %f\n", ((double)curr_2 - curr_1) / CLOCKS_PER_SEC); - //curr_1 = curr_2; - - return fixed; -} - -} // anonymous namespace - -namespace fix_variables_ -{ - -FixVariablesResult fixQuboVariables(const compressed_matrix::CompressedMatrix& Q, int method) -{ - //Q needs to be a square matrix - if (Q.numRows() != Q.numCols()) - throw std::invalid_argument("Q's size is not correct."); - - if (!(method == 1 || method == 2)) - throw std::invalid_argument("method must be an integer of 1 or 2."); - - FixVariablesResult ret; - - //check if Q is empty - if (Q.numRows() == 0 || Q.numCols() == 0) - { - ret.offset = 0; - return ret; - } - - int numVariables = Q.numRows(); - - //clock_t curr_1 = clock(); - //clock_t curr_2; - - //uTriQ = triu(Q) + tril(Q,-1)'; //make upper triangular - std::map, double> uTriQMap; - std::set usedVariables; - for (int i = 0; i < Q.numRows(); i++) - { - int start = Q.rowOffsets()[i]; - int end = Q.rowOffsets()[i + 1]; - for (int j = start; j < end; j++) - { - int r = i; - int c = Q.colIndices()[j]; - uTriQMap[std::make_pair(std::min(r, c), std::max(r, c))] += Q.values()[j]; - usedVariables.insert(r); - usedVariables.insert(c); - } - } - - compressed_matrix::CompressedMatrix uTriQ(numVariables, numVariables, uTriQMap); - - double maxAbsValue = 0; - - if (!uTriQ.values().empty()) - maxAbsValue = std::fabs(*std::max_element(uTriQ.values().begin(), uTriQ.values().end(), compareAbs)); - - double ratio = 1; - - if (maxAbsValue != 0) - ratio = static_cast(std::numeric_limits::max()) / maxAbsValue; - - ratio /= static_cast(1LL << 10); - - if (ratio < 1) - ratio = 1; - - std::map, long long int> uTriQMapLLI; - - for (std::map, double>::iterator it = uTriQMap.begin(), end = uTriQMap.end(); it != end; ++it) - uTriQMapLLI[it->first] = static_cast(it->second * ratio); - - compressed_matrix::CompressedMatrix uTriQLLI(numVariables, numVariables, uTriQMapLLI); - - //curr_2 = clock(); - //mexPrintf("Time elapsed_make upper triangular: %f\n", ((double)curr_2 - curr_1) / CLOCKS_PER_SEC); - //curr_1 = curr_2; - - Posiform p = BQPToPosiform(uTriQLLI); - - //curr_2 = clock(); - //mexPrintf("Time elapsed_BQPToPosiform: %f\n", ((double)curr_2 - curr_1) / CLOCKS_PER_SEC); - //curr_1 = curr_2; - - if (method == 1) - { - compressed_matrix::CompressedMatrix A = posiformToImplicationNetwork_1(p); - - //curr_2 = clock(); - //mexPrintf("Time elapsed_posiformToImplicationNetwork_1: %f\n", ((double)curr_2 - curr_1) / CLOCKS_PER_SEC); - //curr_1 = curr_2; - - compressed_matrix::CompressedMatrix R = maxFlow(A); //use this - - //curr_2 = clock(); - //mexPrintf("Time elapsed_maxFlow: %f\n", ((double)curr_2 - curr_1) / CLOCKS_PER_SEC); - //curr_1 = curr_2; - - compressed_matrix::CompressedMatrix symR = makeResidualSymmetric(R); //use this - - //curr_2 = clock(); - //mexPrintf("Time elapsed_makeRSym: %f\n", ((double)curr_2 - curr_1) / CLOCKS_PER_SEC); - //curr_1 = curr_2; - - //add clearing R in makeResidualSym, so here just use symR directly !!! - SCRet scRet = stronglyConnectedComponents(symR); - - - //curr_2 = clock(); - //mexPrintf("Time elapsed_stronglyConnectedComponents: %f\n", ((double)curr_2 - curr_1) / CLOCKS_PER_SEC); - //curr_1 = curr_2; - - std::vector classifiedSC = classifyStronglyConnectedComponents(scRet.S); - - //curr_2 = clock(); - //mexPrintf("Time elapsed_classifyStrongComponents: %f\n", ((double)curr_2 - curr_1) / CLOCKS_PER_SEC); - //curr_1 = curr_2; - - ret.fixedVars = fixVariables(classifiedSC, scRet, numVariables); - - //curr_2 = clock(); - //mexPrintf("Time elapsed_fixVarsUsingOutDegree: %f\n", ((double)curr_2 - curr_1) / CLOCKS_PER_SEC); - //curr_1 = curr_2; - } - else if (method == 2) - { - compressed_matrix::CompressedMatrix A = posiformToImplicationNetwork_2(p); - - //curr_2 = clock(); - //mexPrintf("Time elapsed_posiformToImplicationNetwork_2: %f\n", ((double)curr_2 - curr_1) / CLOCKS_PER_SEC); - //curr_1 = curr_2; - - ret.fixedVars = applyImplication(A); - - //curr_2 = clock(); - //mexPrintf("Time elapsed_applyImplication: %f\n", ((double)curr_2 - curr_1) / CLOCKS_PER_SEC); - //curr_1 = curr_2; - } - - if (ret.fixedVars.size() > numVariables) - throw std::logic_error("ret.fixedVars has wrong size."); - - // remove unused variables from ret.fixedVars - std::vector > updatedFixedVars; - for (int i = 0; i < ret.fixedVars.size(); ++i) - { - if (usedVariables.find(ret.fixedVars[i].first - 1) != usedVariables.end()) // -1 to make it 0-based since usedVariables is 0-based - updatedFixedVars.push_back(ret.fixedVars[i]); - } - - ret.fixedVars = updatedFixedVars; - - ret.newQ = computeNewQAndOffset(uTriQ, ret.fixedVars, ret.offset); - - //curr_2 = clock(); - //mexPrintf("Time elapsed_computeNewQAndOffset: %f\n", ((double)curr_2 - curr_1) / CLOCKS_PER_SEC); - //curr_1 = curr_2; - - return ret; -} - - -std::vector > fixQuboVariablesMap(std::map, double> QMap, int QSize, int mtd) -{ - compressed_matrix::CompressedMatrix QInput(QSize, QSize, QMap); - - FixVariablesResult ret = fixQuboVariables(QInput, mtd); - - return ret.fixedVars; -} - - -} // namespace fix_variables_ diff --git a/dimod/roof_duality/src/fix_variables.hpp b/dimod/roof_duality/src/fix_variables.hpp deleted file mode 100644 index da6dc166c..000000000 --- a/dimod/roof_duality/src/fix_variables.hpp +++ /dev/null @@ -1,45 +0,0 @@ -/** -# Copyright 2018 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 FIX_VARIABLES_HPP_INCLUDED -#define FIX_VARIABLES_HPP_INCLUDED - -#include -#include - -#include "compressed_matrix.hpp" - -namespace fix_variables_ -{ - -struct FixVariablesResult -{ - std::vector > fixedVars; //1-based - compressed_matrix::CompressedMatrix newQ; //0-based - double offset; -}; - -FixVariablesResult fixQuboVariables(const compressed_matrix::CompressedMatrix& Q, int method); - -std::vector > fixQuboVariablesMap(std::map, double> QMap, int QSize, int method); - -} // namespace fix_variables_ - -#endif // FIX_VARIABLES_HPP_INCLUDED - - - diff --git a/setup.py b/setup.py index 17ee114f8..b47bf5055 100644 --- a/setup.py +++ b/setup.py @@ -57,11 +57,6 @@ def build_extensions(self): 'dimod/bqm/common.pyx', 'dimod/discrete/cydiscrete_quadratic_model.pyx', 'dimod/cyvariables.pyx', - # roofduality needs to be treated differently until deprecation - Extension("dimod.roof_duality._fix_variables", - ['dimod/roof_duality/_fix_variables.pyx', - 'dimod/roof_duality/src/fix_variables.cpp'], - include_dirs=['dimod/roof_duality/src/']), ], annotate=bool(os.getenv('CYTHON_ANNOTATE', False)), nthreads=int(os.getenv('CYTHON_NTHREADS', 0)), diff --git a/testscpp/Makefile b/testscpp/Makefile index 2afe7dbe9..586be73ba 100644 --- a/testscpp/Makefile +++ b/testscpp/Makefile @@ -11,11 +11,11 @@ tests_parallel: test_main_parallel.out test_main: test_main.cpp g++ -std=c++11 -Wall -c test_main.cpp - g++ -std=c++11 -Wall test_main.o tests/*.cpp $(ROOT)/dimod/roof_duality/src/fix_variables.cpp -o test_main -I $(SRC) + g++ -std=c++11 -Wall test_main.o tests/*.cpp -o test_main -I $(SRC) test_main_parallel: test_main.cpp g++ -std=c++11 -fopenmp -Wall -c test_main.cpp -o test_main_parallel.o - g++ -std=c++11 -fopenmp -Wall test_main_parallel.o tests/*.cpp $(ROOT)/dimod/roof_duality/src/fix_variables.cpp -o test_main_parallel -I $(SRC) + g++ -std=c++11 -fopenmp -Wall test_main_parallel.o tests/*.cpp -o test_main_parallel -I $(SRC) catch2: git submodule init diff --git a/testscpp/tests/test_roofduality.cpp b/testscpp/tests/test_roofduality.cpp deleted file mode 100644 index 9c8e41bbd..000000000 --- a/testscpp/tests/test_roofduality.cpp +++ /dev/null @@ -1,44 +0,0 @@ -// 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. -// -// ============================================================================= - -#include "../../dimod/roof_duality/src/fix_variables.hpp" -#include "../Catch2/single_include/catch2/catch.hpp" - -namespace fix_variables_ { - -TEST_CASE("Test roof_duality's fixQuboVariables()", "[roofduality]") { - SECTION("Test invalid cases") { - auto nonsquare_matrix = - compressed_matrix::CompressedMatrix(2, 3); - REQUIRE_THROWS_AS(fixQuboVariables(nonsquare_matrix, 2), - std::invalid_argument); - - auto matrix = compressed_matrix::CompressedMatrix(2, 2); - REQUIRE_THROWS_AS(fixQuboVariables(matrix, 0), std::invalid_argument); - } - - SECTION("Test empty case") { - auto empty_matrix = compressed_matrix::CompressedMatrix(0, 0); - FixVariablesResult result = fixQuboVariables(empty_matrix, 2); - - REQUIRE(result.offset == 0); - REQUIRE(result.fixedVars.empty()); - REQUIRE(result.newQ.numCols() == 0); - REQUIRE(result.newQ.numRows() == 0); - } -} - -} // namespace fix_variables_