From aafc0df1c1d202ed30b8f70d8dc72a39c8eb53b4 Mon Sep 17 00:00:00 2001 From: Kotaro Tanahashi Date: Mon, 3 Oct 2022 16:01:25 +0900 Subject: [PATCH] Release/1.3 (#189) - optimize polynomial data structure in the compilation - optimize the construction of sum expressions - optimize the to_qubo method to run faster - fix hash collision to avoid a crash in robin_hood --- .circleci/config.yml | 649 ------------------------ .github/workflows/build_and_upolad.yaml | 3 +- .gitignore | 1 + CMakeLists.txt | 6 +- README.rst | 2 +- benchmark/benchmark.py | 127 +++++ benchmark/benchmark_tsp.py | 75 --- docs/reference/model.rst | 6 +- external/robin_hood.cmake | 2 +- pyqubo/array.py | 2 +- pyqubo/integer/one_hot_enc_integer.py | 2 + pyqubo/package_info.py | 4 +- setup.py | 93 ++-- src/abstract_syntax_tree.hpp | 137 +++-- src/compiler.hpp | 219 +------- src/expand.hpp | 214 ++++++++ src/linkedlist.hpp | 18 + src/main.cpp | 46 +- src/model.hpp | 327 ++++++------ src/poly.hpp | 186 +++++++ src/product.hpp | 67 +++ src/test.cpp | 194 +++++++ src/variables.hpp | 134 +++++ tests/test_array.py | 107 ---- tests/test_constraint.py | 9 - tests/test_express.py | 61 +-- tests/test_integer.py | 24 +- tests/test_logical_constraint.py | 9 - tests/test_model.py | 22 + 29 files changed, 1389 insertions(+), 1357 deletions(-) delete mode 100644 .circleci/config.yml create mode 100644 benchmark/benchmark.py delete mode 100644 benchmark/benchmark_tsp.py create mode 100644 src/expand.hpp create mode 100644 src/linkedlist.hpp create mode 100644 src/poly.hpp create mode 100644 src/product.hpp create mode 100644 src/test.cpp create mode 100644 src/variables.hpp diff --git a/.circleci/config.yml b/.circleci/config.yml deleted file mode 100644 index 57dfbeea..00000000 --- a/.circleci/config.yml +++ /dev/null @@ -1,649 +0,0 @@ -version: 2.1 - -orbs: - win: circleci/windows@2.2.0 - -jobs: - -################################################################################################## -# Linux -################################################################################################## - - test-linux-36: &linux-test-template - docker: - - image: quay.io/pypa/manylinux1_x86_64:latest - - environment: - PYTHON_PATH: /opt/python/cp36-cp36m - - working_directory: ~/repo - - steps: - - - checkout - - - run: &set-python-path - name: set python path - command: | - echo 'export PATH=$PYTHON_PATH/bin/:$PATH' >> $BASH_ENV - source $BASH_ENV - python -V - - - run: &create-virtualenv - name: create virtualenv - command: | - python --version - python -m pip install virtualenv - python -m virtualenv env - - - restore_cache: - keys: - - v2-dependencies-{{ checksum "requirements.txt" }}-{{ .Environment.CIRCLE_JOB }} - - - run: &install-dependencies-template - name: install dependencies - command: | - # Install dependencies such that pre-release versions are not installed when executing `setup.py install` - . env/bin/activate - python --version - pip install -r requirements.txt - - - save_cache: - paths: - - ./env - key: v2-dependencies-{{ checksum "requirements.txt" }}-{{ .Environment.CIRCLE_JOB }} - - - run: &install-package-template - name: install package - command: | - . env/bin/activate - python --version - python setup.py install - - - run: &run-tests-template - name: run unittests - command: | - . env/bin/activate - python --version - coverage run -m unittest discover - - - test-linux-doctest: - <<: *linux-test-template - - steps: - - - checkout - - - run: *set-python-path - - - run: *create-virtualenv - - - restore_cache: - keys: - - v2-dependencies-{{ checksum "requirements.txt" }}-{{ checksum "requirements_doctest.txt" }}-{{ .Environment.CIRCLE_JOB }} - - - run: - name: install dependencies - command: | - . env/bin/activate - python --version - pip install -r requirements.txt - pip install -r requirements_doctest.txt - - - save_cache: - paths: - - ./env - key: v2-dependencies-{{ checksum "requirements.txt" }}-{{ checksum "requirements_doctest.txt" }}-{{ .Environment.CIRCLE_JOB }} - - - run: *install-package-template - - - run: *run-tests-template - - - run: - name: test doc build - command: | - . env/bin/activate - sphinx-build -W -b html docs docs/_build/html - - - run: - name: run doctest - command: | - . env/bin/activate - make doctest - - test-linux-35: - <<: *linux-test-template - environment: - PYTHON_PATH: /opt/python/cp35-cp35m - - test-linux-37: - <<: *linux-test-template - environment: - PYTHON_PATH: /opt/python/cp37-cp37m - - test-linux-38: - <<: *linux-test-template - environment: - PYTHON_PATH: /opt/python/cp38-cp38 - - test-linux-39: - <<: *linux-test-template - environment: - PYTHON_PATH: /opt/python/cp39-cp39 - - test-linux32-35: - <<: *linux-test-template - docker: - - image: quay.io/pypa/manylinux1_i686:latest - environment: - PYTHON_PATH: /opt/python/cp35-cp35m - - test-linux32-36: - <<: *linux-test-template - docker: - - image: quay.io/pypa/manylinux1_i686:latest - environment: - PYTHON_PATH: /opt/python/cp36-cp36m - - test-linux32-37: - <<: *linux-test-template - docker: - - image: quay.io/pypa/manylinux1_i686:latest - environment: - PYTHON_PATH: /opt/python/cp37-cp37m - - test-linux32-38: - <<: *linux-test-template - docker: - - image: quay.io/pypa/manylinux1_i686:latest - environment: - PYTHON_PATH: /opt/python/cp38-cp38 - - test-linux32-39: - <<: *linux-test-template - docker: - - image: quay.io/pypa/manylinux1_i686:latest - environment: - PYTHON_PATH: /opt/python/cp39-cp39 - - deploy-linux: &deploy-linux-template - docker: - - image: quay.io/pypa/manylinux1_x86_64:latest - - environment: - GIT_URL: << pipeline.project.git_url >> - GIT_REVISION: << pipeline.git.revision >> - GIT_TAG: << pipeline.git.tag >> - - working_directory: ~/repo - - steps: - - run: - name: checkout directory - command: | - # For the bug of CircleCI that checkout doesn't work when the tag is added - echo $GIT_URL - echo $GIT_REVISION - echo $GIT_TAG - cd ../ - rm -rf repo - git clone $GIT_URL repo - cd repo - git checkout -b workbranch $GIT_REVISION - - - run: - name: build wheels - command: | - for PYBIN in /opt/python/*/bin; do - echo $PYBIN - if "${PYBIN}/python" -c "import sys; sys.exit(sys.version_info>=(3, 5) and sys.version_info<(3, 10))"; then continue; fi; - "${PYBIN}/pip" install -r requirements.txt - #"${PYBIN}/pip" wheel . -w ./wheelhouse - PATH=$PYBIN:$PATH python setup.py bdist_wheel -d ./wheelhouse - "${PYBIN}/python" setup.py sdist -d ./dist - done - - - run: - name: bundle shared libraries into wheels - command: | - ls ./wheelhouse/ - for whl in ./wheelhouse/pyqubo*.whl; do - auditwheel repair "$whl" -w ./dist - done - - - store_artifacts: - path: ./dist - - - run: &init-pypirc - name: init .pypirc - command: | - echo -e "[pypi]" >> ~/.pypirc - echo -e "username = $PYPI_USERNAME" >> ~/.pypirc - - # for testpypi - #echo -e "repository = https://test.pypi.org/legacy/" >> ~/.pypirc - #echo -e "password = $TEST_PYPI_PASSWORD" >> ~/.pypirc - - # for actual pypi - echo -e "password = $PYPI_PASSWORD" >> ~/.pypirc - - - run: - name: python setup - command: | - echo 'export PATH=/opt/python/cp36-cp36m/bin/:$PATH' >> $BASH_ENV - source $BASH_ENV - python -V - - - run: *create-virtualenv - - - run: &upload-wheel - name: install twine and deploy - command: | - . env/bin/activate - pip install twine==1.15.0 - twine upload --skip-existing ./dist/* - - deploy-linux32: - <<: *deploy-linux-template - docker: - - image: quay.io/pypa/manylinux1_i686:latest - -################################################################################################## -# Mac OSX -################################################################################################## - - test-osx-36: &test-osx-template - macos: - xcode: "11.2.0" - environment: - PYTHON: 3.6.5 - HOMEBREW_NO_AUTO_UPDATE: 1 - - working_directory: ~/repo - - steps: - - checkout - - - restore_cache: - keys: - - v6-brew_cache-{{ .Environment.CIRCLE_JOB }}-xcode11.2.1 - - - run: &brew-update - name: brew update - command: | - # Only when the brew is not updated and Python=3.9, we update the brew - # since the older pyenv doesn't contain Python3.9. - - brew --version - brew_version=`echo "$(brew -v)" | head -n 1 | awk '{print $2}'` - if [ $PYTHON = 3.9.1 -a $brew_version = 2.1.16 ]; then - brew update - fi - brew --version - - - save_cache: - key: v6-brew_cache-{{ .Environment.CIRCLE_JOB }}-xcode11.2.1 - paths: - - /Users/distiller/Library/Caches/Homebrew - - /usr/local/Homebrew - - - run: &install-cmake-pyenv - name: install cmake - command: | - brew --version - brew install cmake - brew install pyenv - - - restore_cache: - keys: - - v4-pyenv-{{ .Environment.CIRCLE_JOB }}-xcode11.2.1 - - - run: &install-python - name: install python - command: | - pyenv install $PYTHON -s - - - save_cache: - paths: - - ~/.pyenv - key: v4-pyenv-{{ .Environment.CIRCLE_JOB }}-xcode11.2.1 - - - run: &create-virtualenv-using-pyenv - name: create virtualenv - command: | - export PYENV_ROOT="$HOME/.pyenv" - export PATH="$PYENV_ROOT/bin:$PATH" - eval "$(pyenv init --path)" - if command -v pyenv 1>/dev/null 2>&1; then - eval "$(pyenv init -)" - fi - pyenv local $PYTHON - which python - python -m pip install virtualenv - python -m virtualenv env - - - restore_cache: - keys: - - v2-dependencies-{{ checksum "requirements.txt" }}-{{ .Environment.CIRCLE_JOB }} - - - run: *install-dependencies-template - - - save_cache: - paths: - - ./env - key: v2-dependencies-{{ checksum "requirements.txt" }}-{{ .Environment.CIRCLE_JOB }} - - - run: *install-package-template - - - run: *run-tests-template - - - run: - name: codecov - command: | - . env/bin/activate - # skip here when local run - if ! ([ -z "$CIRCLE_PROJECT_REPONAME" ] || [ -z "$CIRCLE_PROJECT_USERNAME" ]); then - codecov - fi - - test-osx-35: - <<: *test-osx-template - environment: - PYTHON: 3.5.8 - HOMEBREW_NO_AUTO_UPDATE: 1 - - test-osx-37: - <<: *test-osx-template - environment: - PYTHON: 3.7.5 - HOMEBREW_NO_AUTO_UPDATE: 1 - - test-osx-38: - <<: *test-osx-template - environment: - PYTHON: 3.8.0 - HOMEBREW_NO_AUTO_UPDATE: 1 - - test-osx-39: - <<: *test-osx-template - environment: - PYTHON: 3.9.1 - HOMEBREW_NO_AUTO_UPDATE: 0 - - deploy-osx-36: &deploy-osx-template - macos: - xcode: "11.2.0" - - environment: - PYTHON: 3.6.5 - HOMEBREW_NO_AUTO_UPDATE: 1 - MACOSX_DEPLOYMENT_TARGET: 10.9 - - working_directory: ~/repo - - steps: - - checkout - - - restore_cache: - keys: - - v6-brew_cache-{{ .Environment.CIRCLE_JOB }}-xcode11.2.1 - - - run: *brew-update - - - save_cache: - key: v6-brew_cache-{{ .Environment.CIRCLE_JOB }}-xcode11.2.1 - paths: - - /Users/distiller/Library/Caches/Homebrew - - /usr/local/Homebrew - - - run: *install-cmake-pyenv - - - restore_cache: - keys: - - v4-pyenv-{{ .Environment.CIRCLE_JOB }}-xcode11.2.0 - - - run: *install-python - - - save_cache: - paths: - - ~/.pyenv - key: v4-pyenv-{{ .Environment.CIRCLE_JOB }}-xcode11.2.0 - - - run: *create-virtualenv-using-pyenv - - - restore_cache: - keys: - - v2-dependencies-{{ checksum "requirements.txt" }}-{{ .Environment.CIRCLE_JOB }} - - - run: *install-dependencies-template - - - save_cache: - paths: - - ./env - key: v2-dependencies-{{ checksum "requirements.txt" }}-{{ .Environment.CIRCLE_JOB }} - - - run: *install-package-template - - - run: *run-tests-template - - - run: - name: create bdist_wheel - command: | - . env/bin/activate - python setup.py bdist_wheel - - - store_artifacts: - path: ./dist - - - run: *init-pypirc - - - run: *upload-wheel - - deploy-osx-35: - <<: *deploy-osx-template - environment: - PYTHON: 3.5.8 - HOMEBREW_NO_AUTO_UPDATE: 1 - MACOSX_DEPLOYMENT_TARGET: 10.9 - - deploy-osx-37: - <<: *deploy-osx-template - environment: - PYTHON: 3.7.5 - HOMEBREW_NO_AUTO_UPDATE: 1 - MACOSX_DEPLOYMENT_TARGET: 10.9 - - deploy-osx-38: - <<: *deploy-osx-template - environment: - PYTHON: 3.8.0 - HOMEBREW_NO_AUTO_UPDATE: 1 - MACOSX_DEPLOYMENT_TARGET: 10.9 - - deploy-osx-39: - <<: *deploy-osx-template - environment: - PYTHON: 3.9.1 - HOMEBREW_NO_AUTO_UPDATE: 0 - MACOSX_DEPLOYMENT_TARGET: 10.9 - - -################################################################################################## -# Windows -################################################################################################## - - test-win-36: &test-win-template - executor: - name: win/default - shell: bash.exe - - environment: - PYTHON: 3.6.5 - PYTHON_PATH: C:/Python36/ - - steps: - - checkout - - - run: &win-install-camke-python - name: install cmake and python - command: | - choco install -y cmake - - # Python 3.7 is already installed - if [ $PYTHON != '3.7.3' ]; then - choco install -y python --version "${PYTHON}" - fi - - - run: - name: install pyqubo - command: | - export PATH=$PATH:"C:\Program Files\CMake\bin" - $PYTHON_PATH/python -m pip install -r requirements.txt - $PYTHON_PATH/python setup.py install - - - run: - name: run unittests - command: | - $PYTHON_PATH/python -m unittest - - test-win-37: - <<: *test-win-template - environment: - PYTHON: 3.7.3 - PYTHON_PATH: C:/tools/miniconda3 - - test-win-38: - <<: *test-win-template - environment: - PYTHON: 3.8.5 - PYTHON_PATH: C:/Python38/ - - test-win-39: - <<: *test-win-template - environment: - PYTHON: 3.9.1 - PYTHON_PATH: C:/Python39/ - - test-win-35: - <<: *test-win-template - environment: - PYTHON: 3.5.4 - PYTHON_PATH: C:/Python35/ - - deploy-win-36: &deploy-win-template - executor: - name: win/default - shell: bash.exe - - environment: - PYTHON: 3.6.5 - PYTHON_PATH: C:/Python36/ - - steps: - - checkout - - - run: *win-install-camke-python - - - run: - name: create wheel and test install - command: | - export PATH=$PATH:"C:\Program Files\CMake\bin" - $PYTHON_PATH/python -m pip install -r requirements.txt - $PYTHON_PATH/python -m pip install wheel twine - $PYTHON_PATH/python setup.py bdist_wheel - $PYTHON_PATH/python -m pip install pyqubo --no-index -f dist - - - run: *init-pypirc - - - run: - name: upload wheel - command: | - $PYTHON_PATH/python -m twine upload --skip-existing ./dist/* - shell: bash.exe - - deploy-win-37: - <<: *deploy-win-template - environment: - PYTHON: 3.7.3 - PYTHON_PATH: C:/tools/miniconda3 - - deploy-win-38: - <<: *deploy-win-template - environment: - PYTHON: 3.8.5 - PYTHON_PATH: C:/Python38/ - - deploy-win-39: - <<: *deploy-win-template - environment: - PYTHON: 3.9.1 - PYTHON_PATH: C:/Python39/ - - deploy-win-35: - <<: *deploy-win-template - environment: - PYTHON: 3.5.4 - PYTHON_PATH: C:/Python35/ - - -workflows: - tests: - jobs: - - test-linux-doctest - - - test-osx-35 - - test-osx-36 - - test-osx-37 - - test-osx-38 - - test-osx-39 - - - test-win-35 - - test-win-36 - - test-win-37 - - test-win-38 - - test-win-39 - - - test-linux-35 - - test-linux-36 - - test-linux-37 - - test-linux-38 - - test-linux-39 - - - test-linux32-35 - - test-linux32-36 - - test-linux32-37 - - test-linux32-38 - - test-linux32-39 - - deploy: - jobs: - - - deploy-osx-35: &deploy-template - filters: - tags: - only: /^[0-9]+(\.[0-9]+){2}?/ - branches: - ignore: /.*/ - - deploy-osx-36: - <<: *deploy-template - - deploy-osx-37: - <<: *deploy-template - - deploy-osx-38: - <<: *deploy-template - - deploy-osx-39: - <<: *deploy-template - - deploy-win-35: - <<: *deploy-template - - deploy-win-36: - <<: *deploy-template - - deploy-win-37: - <<: *deploy-template - - deploy-win-38: - <<: *deploy-template - - deploy-win-39: - <<: *deploy-template - - deploy-linux: - <<: *deploy-template - - deploy-linux32: - <<: *deploy-template diff --git a/.github/workflows/build_and_upolad.yaml b/.github/workflows/build_and_upolad.yaml index 6abf4d3a..28e9fb6c 100644 --- a/.github/workflows/build_and_upolad.yaml +++ b/.github/workflows/build_and_upolad.yaml @@ -172,7 +172,6 @@ jobs: path: dist - uses: pypa/gh-action-pypi-publish@release/v1 with: - user: __token__ password: ${{ secrets.PYPI_PASSWORD }} upload_test_pypi: @@ -185,7 +184,7 @@ jobs: path: dist - uses: pypa/gh-action-pypi-publish@release/v1 with: - user: __token__ password: ${{ secrets.PYPI_TEST_PASSWORD }} repository_url: https://test.pypi.org/legacy/ skip_existing: true + verbose: true diff --git a/.gitignore b/.gitignore index 3cf89d0c..96493832 100644 --- a/.gitignore +++ b/.gitignore @@ -41,6 +41,7 @@ __pycache__/ # Distribution / packaging .Python build/ +build_test/ develop-eggs/ dist/ downloads/ diff --git a/CMakeLists.txt b/CMakeLists.txt index 397c8bbf..9db55834 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -18,7 +18,11 @@ include(external/eigen.cmake) include(external/pybind11.cmake) include(external/robin_hood.cmake) -pybind11_add_module(cpp_pyqubo src/main.cpp) +if(DEFINED ENV{build_test}) + add_executable(cpp_pyqubo src/test.cpp) +else() + pybind11_add_module(cpp_pyqubo src/main.cpp) +endif() target_compile_definitions(cpp_pyqubo PRIVATE VERSION_INFO=${PYQUBO_VERSION_INFO}) target_compile_features(cpp_pyqubo PRIVATE cxx_std_17) diff --git a/README.rst b/README.rst index f0fc9984..7e0c9c5e 100644 --- a/README.rst +++ b/README.rst @@ -45,7 +45,7 @@ By calling ``model.to_qubo()``, we get the resulting QUBO. >>> H = (4*s1 + 2*s2 + 7*s3 + s4)**2 >>> model = H.compile() >>> qubo, offset = model.to_qubo() ->>> pprint(qubo) +>>> pprint(qubo) # doctest: +SKIP {('s1', 's1'): -160.0, ('s1', 's2'): 64.0, ('s1', 's3'): 224.0, diff --git a/benchmark/benchmark.py b/benchmark/benchmark.py new file mode 100644 index 00000000..56462eef --- /dev/null +++ b/benchmark/benchmark.py @@ -0,0 +1,127 @@ +from pyqubo import Array, Constraint, Placeholder +import logging +import time +import argparse +import numpy as np +from timeout_decorator import timeout, TimeoutError +from memory_profiler import memory_usage + + +parser = argparse.ArgumentParser() + +logging.basicConfig(level=logging.INFO) +logger = logging.getLogger("benchmark_tsp") + + +def number_partition_with_timeout(n, timeout_sec): + + @timeout(timeout_sec) + def number_partition(n): + t0 = time.time() + s = Array.create('s', n, 'SPIN') + numbers = np.random.randint(0, 10, n) + H = sum([si * ni for si, ni in zip(s, numbers)])**2 + + # Compile model + t1 = time.time() + model = H.compile() + t2 = time.time() + qubo, offset = model.to_qubo(index_label=False) + t3 = time.time() + print("len(qubo)", len(qubo)) + + return t1-t0, t2-t1, t3-t2 + + return number_partition(n) + +def tsp_with_timeout(n_city, timeout_sec): + + @timeout(timeout_sec) + def tsp(n_city): + t0 = time.time() + x = Array.create('c', (n_city, n_city), 'BINARY') + use_for_loop=False + + # Constraint not to visit more than two cities at the same time. + time_const = 0.0 + for i in range(n_city): + # If you wrap the hamiltonian by Const(...), this part is recognized as constraint + time_const += Constraint((sum(x[i, j] for j in range(n_city)) - 1)**2, label="time{}".format(i)) + + # Constraint not to visit the same city more than twice. + city_const = 0.0 + for j in range(n_city): + city_const += Constraint((sum(x[i, j] for i in range(n_city)) - 1)**2, label="city{}".format(j)) + + # distance of route + feed_dict = {} + + if use_for_loop: + distance = 0.0 + for i in range(n_city): + for j in range(n_city): + for k in range(n_city): + # we set the constant distance + d_ij = 10 + distance += d_ij * x[k, i] * x[(k + 1) % n_city, j] + + else: + distance = [] + for i in range(n_city): + for j in range(n_city): + for k in range(n_city): + # we set the constant distance + d_ij = 10 + distance.append(d_ij * x[k, i] * x[(k + 1) % n_city, j]) + distance = sum(distance) + + print("express done") + + # Construct hamiltonian + A = Placeholder("A") + H = distance + + feed_dict["A"] = 1.0 + + # Compile model + t1 = time.time() + model = H.compile() + t2 = time.time() + qubo, offset = model.to_qubo(index_label=False, feed_dict=feed_dict) + t3 = time.time() + + print("len(qubo)", len(qubo)) + + return t1-t0, t2-t1, t3-t2 + + return tsp(n_city) + + + +def measure(problem, step, init_size, max_size, timeout): + if problem == "tsp": + f = tsp_with_timeout + elif problem == "number_partition": + f = number_partition_with_timeout + + for n in range(init_size, max_size+step, step): + try: + max_memory, (express_time, compile_time, to_qubo_time) = memory_usage((f, (n, timeout)), max_usage=True, retval=True) + logger.info("Memory usage is {} MB for n={}".format(max_memory, n)) + logger.info("Elapsed time is {} sec (expression: {} sec, compile: {} sec, to_qubo {} sec), for n={}".format( + express_time+compile_time+to_qubo_time, express_time, compile_time, to_qubo_time, n)) + + except TimeoutError as e: + logger.error("TimeoutError: Elapsed time exceeded {} sec for n_city={}".format(timeout, n)) + break + + +if __name__=="__main__": + parser.add_argument('-p', '--problem', type=str) + parser.add_argument('-m', '--max_size', type=int) + parser.add_argument('-i', '--init_size', type=int) + parser.add_argument('-s', '--step', type=int) + parser.add_argument('-t', '--timeout', type=int) + args = parser.parse_args() + #number_partition_with_timeout(2000, timeout_sec=10) + measure(args.problem, args.step, args.init_size, args.max_size, args.timeout) diff --git a/benchmark/benchmark_tsp.py b/benchmark/benchmark_tsp.py deleted file mode 100644 index 9d540142..00000000 --- a/benchmark/benchmark_tsp.py +++ /dev/null @@ -1,75 +0,0 @@ -from pyqubo import Array, Constraint, Placeholder -import logging -import time -import argparse -from timeout_decorator import timeout, TimeoutError -from memory_profiler import memory_usage - - -parser = argparse.ArgumentParser() - -logging.basicConfig(level=logging.INFO) -logger = logging.getLogger("benchmark_tsp") - -def tsp_with_timeout(n_city, timeout_sec): - - @timeout(timeout_sec) - def tsp(n_city): - t0 = time.time() - x = Array.create('c', (n_city, n_city), 'BINARY') - - # Constraint not to visit more than two cities at the same time. - time_const = 0.0 - for i in range(n_city): - # If you wrap the hamiltonian by Const(...), this part is recognized as constraint - time_const += Constraint((sum(x[i, j] for j in range(n_city)) - 1)**2, label="time{}".format(i)) - - # Constraint not to visit the same city more than twice. - city_const = 0.0 - for j in range(n_city): - city_const += Constraint((sum(x[i, j] for i in range(n_city)) - 1)**2, label="city{}".format(j)) - - # distance of route - distance = 0.0 - for i in range(n_city): - for j in range(n_city): - for k in range(n_city): - # we set the constant distance - d_ij = 10 - distance += d_ij * x[k, i] * x[(k+1)%n_city, j] - - # Construct hamiltonian - A = Placeholder("A") - H = distance + A * (time_const + city_const) - - # Compile model - t1 = time.time() - model = H.compile() - qubo, offset = model.to_qubo(index_label=True, feed_dict={"A": 2.0}) - t2 = time.time() - - return t1-t0, t2-t1 - - return tsp(n_city) - - -def measure(step, init_size, max_size, timeout): - for n_city in range(init_size, max_size+step, step): - try: - max_memory, (express_time, compile_time) = memory_usage((tsp_with_timeout, (n_city, timeout)), max_usage=True, retval=True) - logger.info("Memory usage is {} MB for n_city={}".format(max_memory, n_city)) - logger.info("Elapsed time is {} sec (expression: {} sec, compile: {} sec), for n_city={}".format( - express_time+compile_time, express_time, compile_time, n_city)) - - except TimeoutError as e: - logger.error("TimeoutError: Elapsed time exceeded {} sec for n_city={}".format(timeout, n_city)) - break - - -if __name__=="__main__": - parser.add_argument('-m', '--max_size', type=int) - parser.add_argument('-i', '--init_size', type=int) - parser.add_argument('-s', '--step', type=int) - parser.add_argument('-t', '--timeout', type=int) - args = parser.parse_args() - measure(args.step, args.init_size, args.max_size, args.timeout) diff --git a/docs/reference/model.rst b/docs/reference/model.rst index 91dba975..f431fd40 100644 --- a/docs/reference/model.rst +++ b/docs/reference/model.rst @@ -38,7 +38,7 @@ Model >>> model.to_qubo(index_label=True) # doctest: +SKIP ({(2, 2): 3.0, (1, 1): 2.0, (0, 0): 1.0}, 0.0) >>> model.variables - ['a', 'b', 'c'] + ['c', 'a', 'b'] This indicaretes the mapping of indices and labels as 'c'->0, 'a'->1, 'b'->2 @@ -96,7 +96,7 @@ Model >>> pprint(model.to_qubo(index_label=True)) # doctest: +SKIP ({(0, 0): 3.0, (0, 2): 1.0, (1, 1): 0.0, (1, 2): 1.0, (2, 2): 0.0}, 0.0) >>> model.variables - ['x', 'y', 'z'] + ['z', 'x', 'y'] .. py:method:: to_ising(index_label=False, feed_dict=None) @@ -129,7 +129,7 @@ Model >>> pprint(model.to_ising(index_label=True)) # doctest: +SKIP ({0: 1.75, 1: 0.25, 2: 0.5}, {(0, 2): 0.25, (1, 2): 0.25}, 2.0) >>> model.variables - ['x', 'y', 'z'] + ['z', 'x', 'y'] .. py:method:: to_bqm(index_label=False, feed_dict=None) diff --git a/external/robin_hood.cmake b/external/robin_hood.cmake index b0e9b129..d2060682 100644 --- a/external/robin_hood.cmake +++ b/external/robin_hood.cmake @@ -3,7 +3,7 @@ include(FetchContent) FetchContent_Declare( robin_hood GIT_REPOSITORY https://github.com/martinus/robin-hood-hashing - GIT_TAG 3.11.2 + GIT_TAG 3.11.3 ) FetchContent_GetProperties(robin_hood) diff --git a/pyqubo/array.py b/pyqubo/array.py index 695f0118..5b7cfb1e 100644 --- a/pyqubo/array.py +++ b/pyqubo/array.py @@ -64,7 +64,7 @@ class Array: >>> array[:, 1] # = array[(slice(None), 1)] Array([Binary('x1'), Binary('x3')]) - >>> sum(array[:, 1]) + >>> sum(array[:, 1]) # doctest: +SKIP (Binary('x1') + Binary('x3')) Use list or tuple to select a subset of the array. diff --git a/pyqubo/integer/one_hot_enc_integer.py b/pyqubo/integer/one_hot_enc_integer.py index 9302dda8..3df01ec6 100644 --- a/pyqubo/integer/one_hot_enc_integer.py +++ b/pyqubo/integer/one_hot_enc_integer.py @@ -65,6 +65,8 @@ def __init__(self, label, value_range, strength): express = SubH(lower + sum(i*x for i, x in enumerate(self.array)), label=label) penalty = self.constraint * strength + self._express = express + super().__init__( label=label, value_range=value_range, diff --git a/pyqubo/package_info.py b/pyqubo/package_info.py index 964eb50b..838a27f2 100644 --- a/pyqubo/package_info.py +++ b/pyqubo/package_info.py @@ -1,11 +1,11 @@ # (major, minor, patch, prerelease) -VERSION = (1, 2, 0, "") +VERSION = (1, 3, 0, "") __shortversion__ = '.'.join(map(str, VERSION[:3])) __version__ = '.'.join(map(str, VERSION[:3])) + "".join(VERSION[3:]) __package_name__ = 'pyqubo' -__contact_names__ = 'Recruit Communications Co., Ltd.' +__contact_names__ = 'Recruit Co., Ltd.' __contact_emails__ = 'rco_pyqubo@ml.cocorou.jp' __homepage__ = 'https://pyqubo.readthedocs.io/en/latest/' __repository_url__ = 'https://github.com/recruit-communications/pyqubo' diff --git a/setup.py b/setup.py index 7254c227..74606da9 100644 --- a/setup.py +++ b/setup.py @@ -1,13 +1,11 @@ import os +import re import subprocess import sys -import platform -import re -from setuptools import setup, Extension +from setuptools import Extension, setup from setuptools.command.build_ext import build_ext - # Convert distutils Windows platform specifiers to CMake -A arguments PLAT_TO_CMAKE = { "win32": "Win32", @@ -17,21 +15,11 @@ } -class PackageInfo(object): - def __init__(self, info_file): - with open(info_file) as f: - exec(f.read(), self.__dict__) - self.__dict__.pop('__builtins__', None) - - def __getattribute__(self, name): - return super(PackageInfo, self).__getattribute__(name) - - -package_info = PackageInfo(os.path.join('pyqubo', 'package_info.py')) - - +# A CMakeExtension needs a sourcedir instead of a file list. +# The name must be the _single_ output extension from the CMake build. +# If you need multiple extensions, see scikit-build. class CMakeExtension(Extension): - def __init__(self, name, sourcedir=''): + def __init__(self, name, sourcedir=""): Extension.__init__(self, name, sources=[]) self.sourcedir = os.path.abspath(sourcedir) @@ -40,11 +28,12 @@ class CMakeBuild(build_ext): def build_extension(self, ext): extdir = os.path.abspath(os.path.dirname(self.get_ext_fullpath(ext.name))) - # required for auto-detection of auxiliary "native" libs + # required for auto-detection & inclusion of auxiliary "native" libs if not extdir.endswith(os.path.sep): extdir += os.path.sep - cfg = "Debug" if self.debug else "Release" + debug = int(os.environ.get("DEBUG", 0)) if self.debug is None else self.debug + cfg = "Debug" if debug else "Release" # CMake lets you override the generator - we need to check this. # Can be set with Conda-Build, for example. @@ -54,12 +43,18 @@ def build_extension(self, ext): # EXAMPLE_VERSION_INFO shows you how to pass a value into the C++ code # from Python. cmake_args = [ - "-DCMAKE_LIBRARY_OUTPUT_DIRECTORY={}".format(extdir), - "-DPYTHON_EXECUTABLE={}".format(sys.executable), - "-DPYQUBOC_VERSION_INFO={}".format(self.distribution.get_version()), - "-DCMAKE_BUILD_TYPE={}".format(cfg), # not used on MSVC, but no harm + f"-DCMAKE_LIBRARY_OUTPUT_DIRECTORY={extdir}", + f"-DPYTHON_EXECUTABLE={sys.executable}", + f"-DCMAKE_BUILD_TYPE={cfg}", # not used on MSVC, but no harm ] build_args = [] + # Adding CMake arguments set as environment variable + # (needed e.g. to build for ARM OSx on conda-forge) + if "CMAKE_ARGS" in os.environ: + cmake_args += [item for item in os.environ["CMAKE_ARGS"].split(" ") if item] + + # In this example, we pass in the version to C++. You might not need to. + cmake_args += [f"-DEXAMPLE_VERSION_INFO={self.distribution.get_version()}"] if self.compiler.compiler_type != "msvc": # Using Ninja-build since it a) is available as a wheel and b) @@ -67,8 +62,17 @@ def build_extension(self, ext): # exported for Ninja to pick it up, which is a little tricky to do. # Users can override the generator with CMAKE_GENERATOR in CMake # 3.15+. - if not cmake_generator: - cmake_args += ["-GNinja"] + if not cmake_generator or cmake_generator == "Ninja": + try: + import ninja # noqa: F401 + + ninja_executable_path = os.path.join(ninja.BIN_DIR, "ninja") + cmake_args += [ + "-GNinja", + f"-DCMAKE_MAKE_PROGRAM:FILEPATH={ninja_executable_path}", + ] + except ImportError: + pass else: @@ -87,18 +91,16 @@ def build_extension(self, ext): # Multi-config generators have a different way to specify configs if not single_config: cmake_args += [ - "-DCMAKE_LIBRARY_OUTPUT_DIRECTORY_{}={}".format(cfg.upper(), extdir) + f"-DCMAKE_LIBRARY_OUTPUT_DIRECTORY_{cfg.upper()}={extdir}" ] build_args += ["--config", cfg] - - if platform.system() == 'Darwin': + if sys.platform.startswith("darwin"): # Cross-compile support for macOS - respect ARCHFLAGS if set archs = re.findall(r"-arch (\S+)", os.environ.get("ARCHFLAGS", "")) if archs: cmake_args += ["-DCMAKE_OSX_ARCHITECTURES={}".format(";".join(archs))] - # Set CMAKE_BUILD_PARALLEL_LEVEL to control the parallel build level # across all generators. if "CMAKE_BUILD_PARALLEL_LEVEL" not in os.environ: @@ -106,24 +108,35 @@ def build_extension(self, ext): # using -j in the build_ext call, not supported by pip or PyPA-build. if hasattr(self, "parallel") and self.parallel: # CMake 3.12+ only. - build_args += ["-j{}".format(self.parallel)] + build_args += [f"-j{self.parallel}"] + + build_temp = os.path.join(self.build_temp, ext.name) + if not os.path.exists(build_temp): + os.makedirs(build_temp) + + subprocess.check_call(["cmake", ext.sourcedir] + cmake_args, cwd=build_temp) + subprocess.check_call(["cmake", "--build", "."] + build_args, cwd=build_temp) + + +class PackageInfo(object): + def __init__(self, info_file): + with open(info_file) as f: + exec(f.read(), self.__dict__) + self.__dict__.pop('__builtins__', None) + + def __getattribute__(self, name): + return super(PackageInfo, self).__getattribute__(name) - if not os.path.exists(self.build_temp): - os.makedirs(self.build_temp) - subprocess.check_call( - ["cmake", ext.sourcedir] + cmake_args, cwd=self.build_temp - ) - subprocess.check_call( - ["cmake", "--build", "."] + build_args, cwd=self.build_temp - ) +package_info = PackageInfo(os.path.join('pyqubo', 'package_info.py')) + packages = ['pyqubo', 'pyqubo.integer', 'pyqubo.utils'] install_requires = [ 'numpy>=1.17.3', - 'dimod>=0.9.14, <0.11', + 'dimod>=0.9.14, <0.12', 'dwave-neal>=0.5.7', 'Deprecated>=1.2.12', #'jij-cimod>=1.0.1' diff --git a/src/abstract_syntax_tree.hpp b/src/abstract_syntax_tree.hpp index d4df0de2..2955c5a1 100644 --- a/src/abstract_syntax_tree.hpp +++ b/src/abstract_syntax_tree.hpp @@ -6,8 +6,11 @@ #include #include #include +#include #include +#include "linkedlist.hpp" + namespace pyqubo { enum class expression_type { @@ -53,64 +56,74 @@ namespace std { } namespace pyqubo { - class add_operator final : public expression { - std::vector> _children; + + using add_list = LinkedList>; + + class add_operator final : public expression { public: - add_operator(const std::shared_ptr& lhs, const std::shared_ptr& rhs) noexcept : _children{lhs, rhs} { - ; - } + pyqubo::add_list* node; - const auto& children() const noexcept { - return _children; + add_operator(const std::shared_ptr& lhs, const std::shared_ptr& rhs) noexcept { + this->node = new pyqubo::add_list(lhs); + this->node->next = new pyqubo::add_list(rhs); } - auto add_child(const std::shared_ptr& expression) noexcept { - _children.emplace_back(expression); + auto create_node(const std::shared_ptr add, const std::shared_ptr new_child){ + auto new_node = new pyqubo::add_list(new_child, add->node); + return new_node; } + add_operator(const std::shared_ptr add, const std::shared_ptr child): + node(create_node(add, child)){} + pyqubo::expression_type expression_type() const noexcept override { return expression_type::add_operator; } std::string to_string() const noexcept override { - return "(" + - std::accumulate(std::begin(_children), std::end(_children), std::string(), [](const auto& acc, const auto& child) { - return acc + (std::size(acc) > 0 ? " + " : "") + child->to_string(); - }) + - ")"; + std::string s; + s += "("; + pyqubo::add_list* next_node = this->node; + int i = 0; + std::string sep = " + "; + while(next_node != nullptr){ + s += next_node->value->to_string(); + s += sep; + next_node = next_node->next; + i++; + } + for(int j=0;j(0); - boost::hash_combine(result, "+"); - - for (const auto& child : _children) { - boost::hash_combine(result, std::hash()(*child)); - } - return result; } bool equals(const std::shared_ptr& other) const noexcept override { - if (!expression::equals(other)) { - return false; - } - - const auto& other_add_operator = std::static_pointer_cast(other); - - if (std::size(_children) != std::size(other_add_operator->_children)) { - return false; - } - - for (auto i = 0; i < static_cast(std::size(_children)); ++i) { - if (!_children[i]->equals(other_add_operator->_children[i])) { + if(expression::equals(other)){ + pyqubo::add_list* other_next_node = std::static_pointer_cast(other)->node; + pyqubo::add_list* next_node = this->node; + while(next_node != nullptr && other_next_node != nullptr){ + if(!next_node->value->equals(other_next_node->value)){ + return false; + } + next_node = next_node->next; + other_next_node = other_next_node->next; + } + if(next_node == nullptr && other_next_node == nullptr){ + return true; + }else{ return false; } + }else{ + return false; } - - return true; } }; @@ -246,7 +259,7 @@ namespace pyqubo { }; class sub_hamiltonian : public variable { - std::shared_ptr _expression; + public: sub_hamiltonian(const std::shared_ptr& expression, const std::string& name) noexcept : variable(name), _expression(expression) { @@ -267,16 +280,17 @@ namespace pyqubo { std::size_t hash() const noexcept override { auto result = variable::hash(); - boost::hash_combine(result, "sub_hamiltonian"); boost::hash_combine(result, std::hash()(*_expression)); - return result; } bool equals(const std::shared_ptr& other) const noexcept override { return variable::equals(other) && _expression->equals(std::static_pointer_cast(other)->_expression); } + + protected: + std::shared_ptr _expression; }; class constraint final : public sub_hamiltonian { @@ -297,7 +311,7 @@ namespace pyqubo { } std::string to_string() const noexcept override { - return "Constraint(" + expression()->to_string() + ", '" + name() + "')"; // conditionは文字列化できない……。 + return "Constraint(" + expression()->to_string() + ", '" + name() + "')"; } std::size_t hash() const noexcept override { @@ -403,38 +417,59 @@ namespace pyqubo { inline std::shared_ptr operator+(const std::shared_ptr& lhs, const std::shared_ptr& rhs) noexcept { if (lhs->expression_type() == expression_type::numeric_literal && rhs->expression_type() == expression_type::numeric_literal) { - return std::make_shared(std::static_pointer_cast(lhs)->value() + std::static_pointer_cast(rhs)->value()); + double left_value = std::static_pointer_cast(lhs)->value(); + double right_value = std::static_pointer_cast(rhs)->value(); + return std::make_shared(left_value + right_value); } - if (lhs->expression_type() == expression_type::numeric_literal && std::static_pointer_cast(lhs)->value() == 0) { + + /*if (lhs->expression_type() == expression_type::numeric_literal && std::static_pointer_cast(lhs)->value() == 0) { return rhs; } if (rhs->expression_type() == expression_type::numeric_literal && std::static_pointer_cast(rhs)->value() == 0) { return lhs; - } + }*/ + return std::make_shared(lhs, rhs); } inline std::shared_ptr operator*(const std::shared_ptr& lhs, const std::shared_ptr& rhs) noexcept { + /*if(lhs->expression_type() != expression_type::numeric_literal && rhs->expression_type() != expression_type::numeric_literal){ + return std::make_shared(lhs, rhs); + }else if (lhs->expression_type() == expression_type::numeric_literal){ + if (rhs->expression_type() == expression_type::numeric_literal) { + return std::make_shared(std::static_pointer_cast(lhs)->value() * std::static_pointer_cast(rhs)->value()); + }else if (std::static_pointer_cast(lhs)->value() == 1) { + return rhs; + }else{ + return std::make_shared(lhs, rhs); + } + } else if (rhs->expression_type() == expression_type::numeric_literal){ + if (std::static_pointer_cast(rhs)->value() == 1) { + return lhs; + }else{ + return std::make_shared(lhs, rhs); + } + }*/ + // constant times constant if (lhs->expression_type() == expression_type::numeric_literal && rhs->expression_type() == expression_type::numeric_literal) { - return std::make_shared(std::static_pointer_cast(lhs)->value() * std::static_pointer_cast(rhs)->value()); - } - - if (lhs->expression_type() == expression_type::numeric_literal && std::static_pointer_cast(lhs)->value() == 1) { - return rhs; - } - - if (rhs->expression_type() == expression_type::numeric_literal && std::static_pointer_cast(rhs)->value() == 1) { - return lhs; + return std::make_shared( + std::static_pointer_cast(lhs)->value() * std::static_pointer_cast(rhs)->value()); + }else{ + return std::make_shared(lhs, rhs); } + + //return std::make_shared(lhs, rhs); + } + inline std::shared_ptr multiply_express(const std::shared_ptr& lhs, const std::shared_ptr& rhs) noexcept { return std::make_shared(lhs, rhs); } template - Result visit(Functor& functor, const std::shared_ptr& expression) noexcept { + Result visit(Functor& functor, const std::shared_ptr& expression) { switch (expression->expression_type()) { case expression_type::add_operator: return functor(std::static_pointer_cast(expression)); diff --git a/src/compiler.hpp b/src/compiler.hpp index 9144f1b6..0a6897e1 100644 --- a/src/compiler.hpp +++ b/src/compiler.hpp @@ -14,210 +14,29 @@ #include "abstract_syntax_tree.hpp" #include "model.hpp" +#include "expand.hpp" +#include "product.hpp" +#include "variables.hpp" namespace pyqubo { - // Expand to polynomial. - - // TODO: ペナルティを最後ではなく途中で足し合わせられるか検討する。もし途中で足し合わせられるなら、戻り値が一つになって嬉しい。 - - class expand final { - robin_hood::unordered_map _sub_hamiltonians; - robin_hood::unordered_map>> _constraints; - variables* _variables; - - public: - auto operator()(const std::shared_ptr& expression, variables* variables) noexcept { - _sub_hamiltonians = {}; - _constraints = {}; - _variables = variables; - - auto [polynomial, penalty] = visit>(*this, expression); - - return std::tuple{polynomial + penalty, _sub_hamiltonians, _constraints}; - } - - auto operator()(const std::shared_ptr& add_operator) noexcept { - // 汚いコードでごめんなさい。パフォーマンスのためなので、ご容赦を。 - - const auto& merge = [](auto& polyominal, const auto& other) { - for (const auto& [product, coefficient] : other) { - const auto [it, emplaced] = polyominal.emplace(product, coefficient); - - if (!emplaced) { - it->second = it->second + coefficient; - } - } - }; - - auto polynomial = pyqubo::polynomial{}; - auto penalty = pyqubo::polynomial{}; - - for (const auto& child: add_operator->children()) { - const auto [child_polynomial, child_penalty] = visit>(*this, child); - - merge(polynomial, child_polynomial); - merge(penalty, child_penalty); - } - - return std::tuple{polynomial, penalty}; - } - - auto operator()(const std::shared_ptr& mul_operator) noexcept { - const auto [l_polynomial, l_penalty] = visit>(*this, mul_operator->lhs()); - const auto [r_polynomial, r_penalty] = visit>(*this, mul_operator->rhs()); - - return std::tuple{l_polynomial * r_polynomial, l_penalty + r_penalty}; - } - - auto operator()(const std::shared_ptr& binary_variable) noexcept { - return std::tuple{polynomial{{{_variables->index(binary_variable->name())}, std::make_shared(1)}}, polynomial{}}; - } - - auto operator()(const std::shared_ptr& spin_variable) noexcept { - return std::tuple{polynomial{{{_variables->index(spin_variable->name())}, std::make_shared(2)}, {{}, std::make_shared(-1)}}, polynomial{}}; - } - - auto operator()(const std::shared_ptr& place_holder_variable) noexcept { - return std::tuple{polynomial{{{}, place_holder_variable}}, polynomial{}}; - } - - auto operator()(const std::shared_ptr& sub_hamiltonian) noexcept { - const auto [polynomial, penalty] = visit>(*this, sub_hamiltonian->expression()); - - _sub_hamiltonians.emplace(sub_hamiltonian->name(), polynomial); - - return std::tuple{polynomial, penalty}; - } - - auto operator()(const std::shared_ptr& constraint) noexcept { - const auto [polynomial, penalty] = visit>(*this, constraint->expression()); - - _constraints.emplace(constraint->name(), std::pair{polynomial, constraint->condition()}); - - return std::tuple{polynomial, penalty}; - } - - auto operator()(const std::shared_ptr& with_penalty) noexcept { - const auto [e_polynomial, e_penalty] = visit>(*this, with_penalty->expression()); - const auto [p_polynomial, p_penalty] = visit>(*this, with_penalty->penalty()); - - return std::tuple{e_polynomial, e_penalty + p_penalty + p_polynomial}; - } - - auto operator()(const std::shared_ptr& user_defined_expression) noexcept { - return visit>(*this, user_defined_expression->expression()); - } - - auto operator()(const std::shared_ptr& numeric_literal) noexcept { - return std::tuple{polynomial{{{}, numeric_literal}}, polynomial{}}; - } - }; - - // Convert to quadratic polynomial. - - inline std::optional> find_replacing_pair(const pyqubo::polynomial& polynomial) noexcept { - auto counts = [&] { - auto result = std::map, int>{}; - - for (const auto& [product, _] : polynomial) { - if (std::size(product.indexes()) <= 2) { - continue; - } - - for (auto it_1 = std::begin(product.indexes()); it_1 != std::prev(std::end(product.indexes())); ++it_1) { - for (auto it_2 = std::next(it_1); it_2 != std::end(product.indexes()); ++it_2) { - const auto [it, emplaced] = result.emplace(std::pair{*it_1, *it_2}, 1); - - if (!emplaced) { - it->second++; - } - } - } - } - - return result; - }(); - - if (std::size(counts) == 0) { - return std::nullopt; - } - - const auto it = std::max_element(std::begin(counts), std::end(counts), [](const auto& count_1, const auto& count_2) { - return count_1.second < count_2.second; - }); - - return it->first; - } - - inline auto convert_to_quadratic(const pyqubo::polynomial& polynomial, double strength, variables* variables) noexcept { - auto result = polynomial; - - for (;;) { - const auto replacing_pair = find_replacing_pair(result); - - if (!replacing_pair) { - break; - } - - const auto replacing_pair_index = variables->index(variables->name(replacing_pair->first) + " * " + variables->name(replacing_pair->second)); - - // replace. - - for (;;) { - auto it = std::find_if(std::begin(result), std::end(result), [&](const auto& term) { - return std::binary_search(std::begin(term.first.indexes()), std::end(term.first.indexes()), replacing_pair->first) && std::binary_search(std::begin(term.first.indexes()), std::end(term.first.indexes()), replacing_pair->second); - }); - - if (it == std::end(result)) { - break; - } - - const auto indexes = [&] { - auto result = pyqubo::indexes{}; - - std::copy_if(std::begin(it->first.indexes()), std::end(it->first.indexes()), std::back_inserter(result), [&](const auto& index) { - return index != replacing_pair->first && index != replacing_pair->second; - }); - - result.emplace_back(replacing_pair_index); - - return result; - }(); - const auto expression = it->second; - - result.erase(it); - result.emplace(product(indexes), expression); - } - - // insert. - - const auto emplace_term = [](pyqubo::polynomial& polynomial, const pyqubo::product& product, const std::shared_ptr& coefficient) { - const auto [it, emplaced] = polynomial.emplace(product, coefficient); - - if (!emplaced) { - it->second = it->second + coefficient; - } - }; - - // clang-format off - emplace_term(result, product{replacing_pair_index }, std::make_shared(strength * 3)); - emplace_term(result, product{replacing_pair->first, replacing_pair_index }, std::make_shared(strength * -2)); - emplace_term(result, product{replacing_pair->second, replacing_pair_index }, std::make_shared(strength * -2)); - emplace_term(result, product{replacing_pair->first, replacing_pair->second}, std::make_shared(strength * 1)); - // clang-format on - } - - return result; - } - - // Compile. - - inline auto compile(const std::shared_ptr& expression, double strength) noexcept { + // Compile. + inline auto compile(const std::shared_ptr& express, const std::shared_ptr& strength) noexcept { auto variables = pyqubo::variables(); - const auto [polynomial, sub_hamiltonians, constraints] = expand()(expression, &variables); - const auto quadratic_polynomial = convert_to_quadratic(polynomial, strength, &variables); + const auto [polynomial, sub_hamiltonians, constraints] = expand()(express, &variables); + + //std::cout << "compile" << polynomial.to_string() << std::endl; + const auto quadratic_polynomial = convert_to_quadratic(*(polynomial.get_terms()), strength, &variables); + + + /*for(auto [key, val]: quadratic_polynomial){ + std::cout << "quadratic_polynomial " << key.to_string() << ", " << val->to_string() << std::endl; + } + for(auto [key, val]: sub_hamiltonians){ + std::cout << "sub_hamiltonians " << key << ", " << val.to_string() << std::endl; + }*/ + return model(quadratic_polynomial, sub_hamiltonians, constraints, variables); } -} +} \ No newline at end of file diff --git a/src/expand.hpp b/src/expand.hpp new file mode 100644 index 00000000..3e1a4879 --- /dev/null +++ b/src/expand.hpp @@ -0,0 +1,214 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "abstract_syntax_tree.hpp" +#include "product.hpp" +#include "variables.hpp" +#include "poly.hpp" + +namespace pyqubo { + // Expand to polynomial. + + class expand final { + robin_hood::unordered_map _sub_hamiltonians; + robin_hood::unordered_map>> _constraints; + variables* _variables; + + public: + auto operator()(const std::shared_ptr& expression, variables* variables) noexcept { + _sub_hamiltonians = {}; + _constraints = {}; + _variables = variables; + + auto [polynomial, penalty] = visit>(*this, expression); + polynomial = polynomial + penalty; + return std::tuple{polynomial, _sub_hamiltonians, _constraints}; + } + + auto operator()(const std::shared_ptr& add_operator) noexcept { + + auto polynomial = pyqubo::poly(); + auto penalty = pyqubo::poly(); + + add_list* next_node = add_operator->node; + auto [child_polynomial, child_penalty] = visit>(*this, next_node->value); + polynomial = polynomial + child_polynomial; + penalty = penalty + child_penalty; + + next_node = next_node->next; + while(next_node != nullptr){ + auto [child_polynomial_tmp, child_penalty_tmp] = visit>(*this, next_node->value); + polynomial = polynomial + child_polynomial_tmp; + penalty = penalty + child_penalty_tmp; + next_node = next_node->next; + } + + return std::tuple{polynomial, penalty}; + } + + auto operator()(const std::shared_ptr& mul_operator) noexcept { + auto [l_polynomial, l_penalty] = visit>(*this, mul_operator->lhs()); + auto [r_polynomial, r_penalty] = visit>(*this, mul_operator->rhs()); + return std::tuple{l_polynomial * r_polynomial, l_penalty + r_penalty}; + } + + auto operator()(const std::shared_ptr& binary_variable) noexcept { + auto p1 = poly(std::make_shared(1), new product({_variables->index(binary_variable->name())})); + return std::tuple{ + p1, + poly() + }; + } + + auto operator()(const std::shared_ptr& spin_variable) noexcept { + polynomial* p = new polynomial{{ + {{_variables->index(spin_variable->name())}, std::make_shared(2)}, + {{}, std::make_shared(-1)} + }}; + return std::tuple{poly(p), poly()}; + } + + auto operator()(const std::shared_ptr& place_holder_variable) noexcept { + return std::tuple{ + poly(place_holder_variable, new product({})), + poly() + }; + } + + auto operator()(const std::shared_ptr& sub_hamiltonian) noexcept { + const auto [polynomial, penalty] = visit>(*this, sub_hamiltonian->expression()); + _sub_hamiltonians.emplace(sub_hamiltonian->name(), polynomial.copy()); + return std::tuple{polynomial, penalty}; + } + + auto operator()(const std::shared_ptr& constraint) noexcept { + const auto [polynomial, penalty] = visit>(*this, constraint->expression()); + _constraints.emplace(constraint->name(), std::pair{polynomial.copy(), constraint->condition()}); + return std::tuple{polynomial, penalty}; + } + + auto operator()(const std::shared_ptr& with_penalty) noexcept { + auto [e_polynomial, e_penalty] = visit>(*this, with_penalty->expression()); + auto [p_polynomial, p_penalty] = visit>(*this, with_penalty->penalty()); + return std::tuple{e_polynomial, p_polynomial}; + } + + auto operator()(const std::shared_ptr& user_defined_expression) noexcept { + return visit>(*this, user_defined_expression->expression()); + } + + auto operator()(const std::shared_ptr& numeric_literal) noexcept { + return std::tuple{ + poly(numeric_literal, new product({})), + poly() + }; + } + }; + + // Convert to quadratic polynomial. + + inline std::optional> find_replacing_pair(const pyqubo::polynomial& polynomial) noexcept { + auto counts = [&] { + auto result = std::map, int>{}; + + for (const auto& [product, _] : polynomial) { + if (std::size(product.indexes()) <= 2) { + continue; + } + + for (auto it_1 = std::begin(product.indexes()); it_1 != std::prev(std::end(product.indexes())); ++it_1) { + for (auto it_2 = std::next(it_1); it_2 != std::end(product.indexes()); ++it_2) { + const auto [it, emplaced] = result.emplace(std::pair{*it_1, *it_2}, 1); + + if (!emplaced) { + it->second++; + } + } + } + } + + return result; + }(); + + if (std::size(counts) == 0) { + return std::nullopt; + } + + const auto it = std::max_element(std::begin(counts), std::end(counts), [](const auto& count_1, const auto& count_2) { + return count_1.second < count_2.second; + }); + + return it->first; + } + + inline auto convert_to_quadratic(const pyqubo::polynomial& polynomial, std::shared_ptr strength, variables* variables) noexcept { + auto result = polynomial; + + for (;;) { + const auto replacing_pair = find_replacing_pair(result); + + if (!replacing_pair) { + break; + } + + const auto replacing_pair_index = variables->index(variables->name(replacing_pair->first) + " * " + variables->name(replacing_pair->second)); + + // replace. + + for (;;) { + auto it = std::find_if(std::begin(result), std::end(result), [&](const auto& term) { + return std::binary_search(std::begin(term.first.indexes()), std::end(term.first.indexes()), replacing_pair->first) && std::binary_search(std::begin(term.first.indexes()), std::end(term.first.indexes()), replacing_pair->second); + }); + + if (it == std::end(result)) { + break; + } + + const auto indexes = [&] { + auto result = pyqubo::indexes{}; + + std::copy_if(std::begin(it->first.indexes()), std::end(it->first.indexes()), std::back_inserter(result), [&](const auto& index) { + return index != replacing_pair->first && index != replacing_pair->second; + }); + + result.emplace_back(replacing_pair_index); + + return result; + }(); + const auto expression = it->second; + result.erase(it); + result.emplace(product(indexes), expression); + } + + // insert. + + const auto emplace_term = [](pyqubo::polynomial& polynomial, const pyqubo::product& product, const std::shared_ptr& coefficient) { + const auto [it, emplaced] = polynomial.emplace(product, coefficient); + + if (!emplaced) { + it->second = it->second + coefficient; + } + }; + + // clang-format off + emplace_term(result, product{replacing_pair_index }, std::make_shared(3) * strength); + emplace_term(result, product{replacing_pair->first, replacing_pair_index }, std::make_shared(-2) * strength); + emplace_term(result, product{replacing_pair->second, replacing_pair_index }, std::make_shared(-2) * strength); + emplace_term(result, product{replacing_pair->first, replacing_pair->second}, strength); + // clang-format on + } + + return result; + } +} diff --git a/src/linkedlist.hpp b/src/linkedlist.hpp new file mode 100644 index 00000000..93f1ca1d --- /dev/null +++ b/src/linkedlist.hpp @@ -0,0 +1,18 @@ +#pragma once +#include + + +template +class LinkedList { +public: + T value; + LinkedList* next = nullptr; + LinkedList(T _value): value(_value){} + LinkedList(T _value, LinkedList* _next): value(_value), next(_next){} + + ~LinkedList(){ + if(next != nullptr){ + delete next; + } + } +}; diff --git a/src/main.cpp b/src/main.cpp index acff1bf4..cd2c28d9 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -11,6 +11,7 @@ #include #include "abstract_syntax_tree.hpp" +#include "expand.hpp" #include "compiler.hpp" namespace py = pybind11; @@ -18,6 +19,7 @@ using namespace py::literals; PYBIND11_MODULE(cpp_pyqubo, m) { m.doc() = "pyqubo C++ binding"; + py::class_>(m, "Base") .def("__add__", [](const std::shared_ptr& expression, const std::shared_ptr& other) { @@ -72,9 +74,14 @@ PYBIND11_MODULE(cpp_pyqubo, m) { }) .def( "compile", [](const std::shared_ptr& expression, double strength) { - return pyqubo::compile(expression, strength); + return pyqubo::compile(expression, std::make_shared(strength)); }, py::arg("strength") = 5) + .def( + "compile", [](const std::shared_ptr& expression, const std::shared_ptr& placeholder_strength) { + return pyqubo::compile(expression, placeholder_strength); + }, + py::arg("strength")) .def("__hash__", [](const pyqubo::expression& expression) { // 必要？ return std::hash()(expression); }) @@ -83,17 +90,23 @@ PYBIND11_MODULE(cpp_pyqubo, m) { .def("__repr__", &pyqubo::expression::to_string); py::class_, pyqubo::expression>(m, "Add") - .def("__iadd__", [](std::shared_ptr& add_operator, const std::shared_ptr& other) { - add_operator->add_child(other); - - return add_operator; + .def("__add__", [](const std::shared_ptr& add_operator, const std::shared_ptr& other) { + return std::make_shared(add_operator, other); }) - .def("__iadd__", [](std::shared_ptr& add_operator, double other) { - add_operator->add_child(std::make_shared(other)); - - return add_operator; + .def("__add__", [](const std::shared_ptr& add_operator, double other) { + return std::make_shared(add_operator, std::make_shared(other)); + }) + .def("__radd__", [](const std::shared_ptr& add_operator, double other) { + return std::make_shared(add_operator, std::make_shared(other)); + }) + .def("__sub__", [](const std::shared_ptr& add_operator, const std::shared_ptr& other) { + return std::make_shared(add_operator, std::make_shared(-1) * other); + }) + .def("__sub__", [](const std::shared_ptr& add_operator, double other) { + return std::make_shared(add_operator, std::make_shared(-other)); }); + py::class_, pyqubo::expression>(m, "Binary") .def(py::init()); @@ -110,7 +123,10 @@ PYBIND11_MODULE(cpp_pyqubo, m) { .def(py::init&, const std::string&, const std::function&>(), py::arg("hamiltonian"), py::arg("label"), py::arg("condition") = py::cpp_function([](double x) { return x == 0; })); py::class_, pyqubo::expression>(m, "WithPenalty") - .def(py::init&, const std::shared_ptr&, const std::string&>()); + .def(py::init&, const std::shared_ptr&, const std::string&>()) + .def_property_readonly("express", &pyqubo::with_penalty::expression) + .def_property_readonly("penalty", &pyqubo::with_penalty::penalty); + py::class_, pyqubo::expression>(m, "UserDefinedExpress") .def(py::init&>()); @@ -170,7 +186,9 @@ PYBIND11_MODULE(cpp_pyqubo, m) { }(); return solution.sample().at(name_and_indexes); - }); + }) + .def("value", &pyqubo::solution::evaluate) + .def("__repr__", &pyqubo::solution::to_string); py::class_(m, "Model") .def_property_readonly("variables", &pyqubo::model::variable_names) .def( @@ -190,9 +208,9 @@ PYBIND11_MODULE(cpp_pyqubo, m) { .def( "to_qubo", [](const pyqubo::model& model, bool index_label, const std::unordered_map& feed_dict) { if (!index_label) { - return py::cast(model.to_bqm(feed_dict, cimod::Vartype::BINARY).to_qubo()); - } else { - return py::cast(model.to_bqm(feed_dict, cimod::Vartype::BINARY).to_qubo()); + return py::cast(model.to_qubo_string(feed_dict)); + }else{ + return py::cast(model.to_qubo_int(feed_dict)); } }, py::arg("index_label") = false, py::arg("feed_dict") = std::unordered_map{}) diff --git a/src/model.hpp b/src/model.hpp index 28c9f5b5..d3cd5131 100644 --- a/src/model.hpp +++ b/src/model.hpp @@ -17,172 +17,108 @@ #include #include "abstract_syntax_tree.hpp" +#include "expand.hpp" +#include "product.hpp" +#include "variables.hpp" -namespace pyqubo { - class variables final { - robin_hood::unordered_map _indexes; - robin_hood::unordered_map _names; - - public: - variables() noexcept : _indexes{}, _names{} { - ; - } - - auto index(const std::string& variable_name) noexcept { - const auto [it, emplaced] = _indexes.emplace(variable_name, std::size(_indexes)); - - if (emplaced) { - _names.emplace(it->second, variable_name); - } - - return it->second; - } - - const auto& name(int index) const noexcept { - return _names.find(index)->second; - } - - auto names() const noexcept { - auto result = std::vector(std::size(_names)); - - for (const auto& [index, name] : _names) { - result[index] = name; - } - - return result; - } - }; - - using indexes = boost::container::small_vector; - - class product final { - pyqubo::indexes _indexes; - std::size_t _hash; - - public: - product(const pyqubo::indexes& indexes) noexcept : _indexes(indexes), _hash(boost::hash_range(std::begin(indexes), std::end(indexes))) { - ; - } - - product(std::initializer_list init) noexcept : product(pyqubo::indexes(init)) { - ; - } - - const auto& indexes() const noexcept { - return _indexes; - } - - friend std::hash; - }; - - inline auto operator*(const product& product_1, const product& product_2) noexcept { - return product([&] { - auto result = indexes{}; - - std::set_union(std::begin(product_1.indexes()), std::end(product_1.indexes()), - std::begin(product_2.indexes()), std::end(product_2.indexes()), - std::back_inserter(result)); - - return result; - }()); - } - - inline bool operator==(const product& product_1, const product& product_2) noexcept { - return product_1.indexes() == product_2.indexes(); - } -} - -namespace std { - template <> - struct hash { - auto operator()(const pyqubo::product& product) const noexcept { - return product._hash; - } - }; -} namespace pyqubo { - // std::variantを使用してzeroやmonomialな場合の処理削減をやってみたのですが、パフォーマンスは向上しませんでした。なので、unordered_map一本でやります。 - // よく考えれば、要素数が0の場合の処理とかはunordered_mapの中でやっていそうですし。。。 - - using polynomial = robin_hood::unordered_map>; - - inline auto operator+(const polynomial& polynomial_1, const polynomial& polynomial_2) noexcept { - auto result = polynomial_1; - - for (const auto& [product, coefficient] : polynomial_2) { - const auto [it, emplaced] = result.emplace(product, coefficient); - - if (!emplaced) { - it->second = it->second + coefficient; - } - } - return result; - } - - inline auto operator*(const polynomial& polynomial_1, const polynomial& polynomial_2) noexcept { - auto result = polynomial{}; - - for (const auto& [product_1, coefficient_1] : polynomial_1) { - for (const auto& [product_2, coefficient_2] : polynomial_2) { - const auto [it, emplaced] = result.emplace(product_1 * product_2, coefficient_1 * coefficient_2); - - if (!emplaced) { - it->second = it->second + coefficient_1 * coefficient_2; - } - } - } - - return result; - } class evaluate final { std::unordered_map _feed_dict; public: - evaluate(const std::unordered_map& feed_dict) noexcept : _feed_dict(feed_dict) { + evaluate(const std::unordered_map& feed_dict) : _feed_dict(feed_dict) { ; } - auto operator()(const std::shared_ptr& expression) const noexcept { + auto operator()(const std::shared_ptr& expression) const { return visit(*this, expression); } - auto operator()(const std::shared_ptr& add_operator) const noexcept { - return std::accumulate(std::begin(add_operator->children()), std::end(add_operator->children()), 0.0, [&](const auto& acc, const auto& child) { + auto operator()(const std::shared_ptr& add_operator) const { + double acc = 0.0; + add_list* next_node = add_operator->node; + while(next_node != nullptr){ + acc = acc + visit(*this, next_node->value); + next_node = next_node->next; + } + return acc; + /*return std::accumulate(std::begin(add_operator->children()), std::end(add_operator->children()), 0.0, [&](const auto& acc, const auto& child) { return acc + visit(*this, child); - }); + });*/ } - auto operator()(const std::shared_ptr& mul_operator) const noexcept { + auto operator()(const std::shared_ptr& mul_operator) const { return visit(*this, mul_operator->lhs()) * visit(*this, mul_operator->rhs()); } - auto operator()(const std::shared_ptr& place_holder_variable) const noexcept { - return _feed_dict.at(place_holder_variable->name()); + auto operator()(const std::shared_ptr& place_holder_variable) const { + + auto found = _feed_dict.find(place_holder_variable->name()); + if(found != _feed_dict.end()){ + return found->second; + }else{ + std::string err_msg = "the value of " + place_holder_variable->name() + " is not provided in feed_dict."; + throw std::invalid_argument(err_msg); + return 0.0; + } } - auto operator()(const std::shared_ptr& user_defined_expression) const noexcept { + auto operator()(const std::shared_ptr& user_defined_expression) const { return visit(*this, user_defined_expression->expression()); } - auto operator()(const std::shared_ptr& numeric_literal) const noexcept { + auto operator()(const std::shared_ptr& numeric_literal) const { return numeric_literal->value(); } }; + class solution final { std::unordered_map _sample; double _energy; std::unordered_map _sub_hamiltonians; std::unordered_map> _constraints; + const std::unordered_map _feed_dict; + const std::string _vartype; + variables _variables; public: - solution(const std::unordered_map sample, double energy, const std::unordered_map& sub_hamiltonians, const std::unordered_map>& constraints) noexcept : _sample(sample), _energy(energy), _sub_hamiltonians(sub_hamiltonians), _constraints(constraints) { + solution( + const std::unordered_map sample, + double energy, + const std::unordered_map& sub_hamiltonians, + const std::unordered_map>& constraints, + const std::unordered_map feed_dict, + const std::string vartype, + variables variables + ) noexcept : + _sample(sample), + _energy(energy), + _sub_hamiltonians(sub_hamiltonians), + _constraints(constraints), + _feed_dict(feed_dict), + _vartype(vartype), + _variables(variables) { ; } + std::string to_string(){ + std::string s = "DecodedSolution({"; + int counter = 0; + for(auto [k, v]: _sample){ + s += k + ":" + std::to_string(v); + if(counter != _sample.size() - 1){ + s += ", "; + } + counter ++; + } + s += "}, energy=" + std::to_string(_energy) + ")"; + return s; + } + const auto& sample() const noexcept { return _sample; } @@ -198,12 +134,43 @@ namespace pyqubo { const auto& constraints() const noexcept { return _constraints; } + + auto evaluate(const std::shared_ptr& expression) { + + const auto [polynomial, sub_hamiltonians, constraints] = pyqubo::expand()(expression, &_variables); + //std::cout << _variables.to_string(); + //std::cout << "compile" << polynomial.to_string() << std::endl; + auto& poly_terms = *polynomial.terms; + + const auto evaluate = pyqubo::evaluate(_feed_dict); + const auto evaluate_polynomial = [&](const auto& poly_terms) { + return std::accumulate(std::begin(poly_terms), std::end(poly_terms), 0.0, [&](const auto acc, const auto& term) { + return acc + + std::accumulate(std::begin(term.first.indexes()), std::end(term.first.indexes()), 1, [&](const auto acc, const auto& index) { + const auto value = _sample.at(_variables.name(index)); + return acc * (_vartype == "BINARY" ? value : (value + 1) / 2); + }) * evaluate(term.second); + }); + }; + const auto energy = evaluate_polynomial(poly_terms); + + // check constraints + for (const auto& [name, pair] : constraints) { + const auto& [polynomial, condition] = pair; + auto& const_poly_terms = *polynomial.terms; + const auto const_energy = evaluate_polynomial(const_poly_terms); + if(const_energy > 0){ + throw std::runtime_error("constraint: " + name + " is broken."); + } + } + return energy; + } }; class model final { - polynomial _quadratic_polynomial; - robin_hood::unordered_map _sub_hamiltonians; - robin_hood::unordered_map>> _constraints; + const polynomial _quadratic_polynomial; + robin_hood::unordered_map _sub_hamiltonians; // コンパイル中にpolyのコピーをしたかチェック + robin_hood::unordered_map>> _constraints; variables _variables; static auto to_cimod_vartype(const std::string vartype) noexcept { @@ -211,7 +178,7 @@ namespace pyqubo { } public: - model(const polynomial& quadratic_polynomial, const robin_hood::unordered_map& sub_hamiltonians, const robin_hood::unordered_map>>& constraints, const variables& variables) noexcept : _quadratic_polynomial(quadratic_polynomial), _sub_hamiltonians(sub_hamiltonians), _constraints(constraints), _variables(variables) { + model(const polynomial &quadratic_polynomial, const robin_hood::unordered_map& sub_hamiltonians, const robin_hood::unordered_map>>& constraints, const variables& variables) noexcept : _quadratic_polynomial(quadratic_polynomial), _sub_hamiltonians(sub_hamiltonians), _constraints(constraints), _variables(variables) { ; } @@ -220,7 +187,8 @@ namespace pyqubo { } template - auto to_bqm_parameters(const std::unordered_map& feed_dict) const noexcept { // 不格好でごめんなさい。PythonのBinaryQuadraticModelを作成可能にするために、このメンバ関数でBinaryQuadraticModelの引数を生成します。 + auto to_bqm_parameters(const std::unordered_map& feed_dict) const { // 不格好でごめんなさい。PythonのBinaryQuadraticModelを作成可能にするために、このメンバ関数でBinaryQuadraticModelの引数を生成します。 + //throw std::runtime_error("test to_qubo."); const auto evaluate = pyqubo::evaluate(feed_dict); auto linear = cimod::Linear{}; @@ -252,12 +220,81 @@ namespace pyqubo { } template - auto to_bqm(const std::unordered_map& feed_dict, cimod::Vartype vartype) const noexcept { + auto to_bqm(const std::unordered_map& feed_dict, cimod::Vartype vartype) const { const auto [linear, quadratic, offset] = to_bqm_parameters(feed_dict); return cimod::BinaryQuadraticModel(linear, quadratic, offset, vartype); } + auto to_qubo_int(const std::unordered_map& feed_dict) const { + + const auto evaluate = pyqubo::evaluate(feed_dict); + auto quadratic = cimod::Quadratic{}; + auto offset = 0.0; + + for (const auto& [product, coefficient] : _quadratic_polynomial) { + const auto coefficient_value = evaluate(coefficient); + + switch (std::size(product.indexes())) { + case 0: { + offset = coefficient_value; // 次数が0の項は1つにまとめられるので、この処理は最大で1回しか実行されません。なので、+=ではなくて=を使用しています。 + break; + } + case 1: { + if(coefficient_value != 0.0){ + quadratic.emplace(std::pair{product.indexes()[0], product.indexes()[0]}, coefficient_value); + } + break; + } + case 2: { + if(coefficient_value != 0.0){ + quadratic.emplace(std::pair{product.indexes()[0], product.indexes()[1]}, coefficient_value); + } + break; + } + default: + throw std::runtime_error("invalid term."); + } + } + + return std::make_tuple(quadratic, offset); + } + + auto to_qubo_string(const std::unordered_map& feed_dict) const { + + const auto evaluate = pyqubo::evaluate(feed_dict); + + auto quadratic = cimod::Quadratic{}; + auto offset = 0.0; + + for (const auto& [product, coefficient] : _quadratic_polynomial) { + const auto coefficient_value = evaluate(coefficient); + + switch (std::size(product.indexes())) { + case 0: { + offset = coefficient_value; // 次数が0の項は1つにまとめられるので、この処理は最大で1回しか実行されません。なので、+=ではなくて=を使用しています。 + break; + } + case 1: { + if(coefficient_value != 0.0){ + quadratic.emplace(std::pair{_variables.name(product.indexes()[0]), _variables.name(product.indexes()[0])}, coefficient_value); + } + break; + } + case 2: { + if(coefficient_value != 0.0){ + quadratic.emplace(std::pair{_variables.name(product.indexes()[0]), _variables.name(product.indexes()[1])}, coefficient_value); + } + break; + } + default: + throw std::runtime_error("invalid term."); // ここには絶対にこないはず。 + } + } + + return std::make_tuple(quadratic, offset); + } + template auto energy(const std::unordered_map& sample, const std::string& vartype, const std::unordered_map& feed_dict) const noexcept { return to_bqm(feed_dict, to_cimod_vartype(vartype)).energy([&] { @@ -277,7 +314,7 @@ namespace pyqubo { } template - auto decode_sample(const std::unordered_map& sample, const std::string& vartype, const std::unordered_map& feed_dict) const noexcept { + auto decode_sample(const std::unordered_map& sample, const std::string& vartype, const std::unordered_map& feed_dict) const { const auto evaluate = pyqubo::evaluate(feed_dict); const auto evaluate_polynomial = [&](const auto& polynomial, const auto& sample) { return std::accumulate(std::begin(polynomial), std::end(polynomial), 0.0, [&](const auto acc, const auto& term) { @@ -297,7 +334,7 @@ namespace pyqubo { auto result = std::unordered_map{}; for (const auto& [name, polynomial] : _sub_hamiltonians) { - result.emplace(name, evaluate_polynomial(polynomial, sample)); + result.emplace(name, evaluate_polynomial(*polynomial.get_terms(), sample)); } return result; @@ -307,17 +344,20 @@ namespace pyqubo { for (const auto& [name, pair] : _constraints) { const auto& [polynomial, condition] = pair; - const auto energy = evaluate_polynomial(polynomial, sample); + const auto energy = evaluate_polynomial(*polynomial.get_terms(), sample); result.emplace(name, std::pair{condition(energy), energy}); } return result; - }()); + }(), + feed_dict, + vartype, + _variables); } template - auto decode_samples(const std::vector>& samples, const std::string& vartype, const std::unordered_map& feed_dict) const noexcept { + auto decode_samples(const std::vector>& samples, const std::string& vartype, const std::unordered_map& feed_dict) const { auto result = std::vector{}; std::transform(std::begin(samples), std::end(samples), std::back_inserter(result), [&](const auto& sample) { @@ -329,7 +369,7 @@ namespace pyqubo { }; template <> - inline auto model::to_bqm_parameters(const std::unordered_map& feed_dict) const noexcept { // メンバ関数を特殊化するときは、クラスの外に書かなければなりません。。。 + inline auto model::to_bqm_parameters(const std::unordered_map& feed_dict) const { // メンバ関数を特殊化するときは、クラスの外に書かなければなりません。。。 const auto evaluate = pyqubo::evaluate(feed_dict); auto linear = cimod::Linear{}; @@ -361,7 +401,7 @@ namespace pyqubo { } template <> - inline auto model::decode_sample(const std::unordered_map& sample, const std::string& vartype, const std::unordered_map& feed_dict) const noexcept { + inline auto model::decode_sample(const std::unordered_map& sample, const std::string& vartype, const std::unordered_map& feed_dict) const { const auto evaluate = pyqubo::evaluate(feed_dict); const auto evaluate_polynomial = [&](const auto& polynomial, const auto& sample) { return std::accumulate(std::begin(polynomial), std::end(polynomial), 0.0, [&](const auto acc, const auto& term) { @@ -389,7 +429,7 @@ namespace pyqubo { auto result = std::unordered_map{}; for (const auto& [name, polynomial] : _sub_hamiltonians) { - result.emplace(name, evaluate_polynomial(polynomial, sample)); + result.emplace(name, evaluate_polynomial(*polynomial.get_terms(), sample)); } return result; @@ -399,12 +439,15 @@ namespace pyqubo { for (const auto& [name, pair] : _constraints) { const auto& [polynomial, condition] = pair; - const auto energy = evaluate_polynomial(polynomial, sample); + const auto energy = evaluate_polynomial(*polynomial.get_terms(), sample); result.emplace(name, std::pair{condition(energy), energy}); } return result; - }()); + }(), + feed_dict, + vartype, + _variables); } } diff --git a/src/poly.hpp b/src/poly.hpp new file mode 100644 index 00000000..e7eb66b7 --- /dev/null +++ b/src/poly.hpp @@ -0,0 +1,186 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include "product.hpp" + + +namespace pyqubo { + + enum class poly_type { + multi_poly, + single_poly + }; + + using Coeff = std::shared_ptr; + + class poly { + public: + poly_type _poly_type = poly_type::single_poly; + pyqubo::polynomial* terms = nullptr; + product* prd = nullptr; + Coeff coeff; + + poly(pyqubo::polynomial* terms): terms(terms){ + _poly_type = poly_type::multi_poly; + } + + poly(Coeff coeff): coeff(coeff){ + this->prd = new product({}); + } + + poly(){ + this->coeff = std::make_shared(0); + this->prd = new product({}); + } + + poly(Coeff coeff, product* prd): coeff(coeff), prd(prd){ + } + + poly copy() const { + if(_poly_type == poly_type::single_poly){ + return poly(this->coeff, this->prd); + }else{ + return poly(new pyqubo::polynomial(*terms)); + } + } + + // polynomial is composed of only numeric constant or not + bool is_numeric() const { + return _poly_type == poly_type::single_poly && prd->indexes().size() == 0; + } + + poly to_multi() const { + if(_poly_type == poly_type::single_poly){ + return poly(new pyqubo::polynomial({{*prd, coeff}})); + }else{ + return *this; + } + } + + pyqubo::polynomial* get_terms() const { + if(_poly_type == poly_type::single_poly){ + return new pyqubo::polynomial({{*prd, coeff}}); + }else{ + return terms; + } + } + + size_t size() const { + if(_poly_type == poly_type::single_poly){ + return 1; + }else{ + return terms->size(); + } + } + + std::string to_string() const { + if(_poly_type == poly_type::single_poly){ + return "single_poly(" + coeff->to_string() + "," + this->prd->to_string() + ")"; + }else{ + std::stringstream ss; + ss << "multi_poly("; + for(const auto& [name, value]: *terms){ + ss << "["; + ss << name.to_string(); + ss << ","; + ss << value->to_string(); + ss << "],"; + } + ss << ")"; + std::string s = ss.str(); + return s; + } + } + }; + + auto multiply_multi_multi(const poly& poly_1, const poly& poly_2){ + auto result = new polynomial{}; + for (const auto& [product_1, coefficient_1] : *poly_1.terms) { + for (const auto& [product_2, coefficient_2] : *poly_2.terms) { + const auto [it, emplaced] = result->emplace(product_1 * product_2, coefficient_1 * coefficient_2); + if (!emplaced) { + it->second = it->second + coefficient_1 * coefficient_2; + } + } + } + return poly(result); + } + + auto& add_multi_single(poly& multi_poly, poly& single_poly){ + //std::cout << "add_multi_single: " << multi_poly.to_string() << " <= " << single_poly.to_string() << "\n"; + const auto [it, emplaced] = multi_poly.terms->emplace(*(single_poly.prd), single_poly.coeff); + if(!emplaced){ + it->second = it->second + single_poly.coeff; + } + return multi_poly; + } + + auto& add_multi_multi(const poly& multi_poly_1, const poly& multi_poly_2) { + //printf("add_multi_multi"); + for (const auto& [product, coefficient] : *multi_poly_2.terms) { + const auto [it, emplaced] = multi_poly_1.terms->emplace(product, coefficient); + if (!emplaced) { + it->second = it->second + coefficient; + } + } + return multi_poly_1; + } + + auto operator*(const poly& poly_1, const poly& poly_2) noexcept { + if(poly_1._poly_type == poly_type::single_poly && poly_2._poly_type == poly_type::single_poly){ + /*if(poly_1.is_numeric() && poly_2.is_numeric()){ + return poly(poly_1.coeff * poly_2.coeff, new product({})); + }else{ + return poly(poly_1.coeff * poly_2.coeff, multiply(poly_1.prd, poly_2.prd)); + }*/ + return poly(poly_1.coeff * poly_2.coeff, multiply(poly_1.prd, poly_2.prd)); + }else if(poly_1._poly_type == poly_type::multi_poly && poly_2._poly_type == poly_type::single_poly){ + return multiply_multi_multi(poly_1, poly_2.to_multi()); + }else if(poly_1._poly_type == poly_type::single_poly && poly_2._poly_type == poly_type::multi_poly){ + return multiply_multi_multi(poly_1.to_multi(), poly_2); + }else{ + return multiply_multi_multi(poly_1, poly_2); + } + } + + auto operator+(poly& poly_1, poly& poly_2) noexcept { + + if(poly_1._poly_type == poly_type::single_poly && poly_2._poly_type == poly_type::single_poly){ + //printf("single poly %s + single poly %s\n", poly_1.coeff->to_string().c_str(), poly_2.coeff->to_string().c_str()); + if(poly_1.prd->equals(*poly_2.prd)){ + return poly(poly_1.coeff + poly_2.coeff, poly_1.prd); + }else{ + auto terms = new pyqubo::polynomial({{*poly_1.prd, poly_1.coeff}, {*poly_2.prd, poly_2.coeff}}); + return poly(terms); + } + }else if(poly_1._poly_type == poly_type::multi_poly && poly_2._poly_type == poly_type::single_poly){ + //std::cout << "operator+ add_multi_single\n"; + return add_multi_single(poly_1, poly_2); + + }else if(poly_1._poly_type == poly_type::single_poly && poly_2._poly_type == poly_type::multi_poly){ + //std::cout << "operator+ add_multi_single\n"; + return add_multi_single(poly_2, poly_1); + }else{ + //std::cout << "operator+ add_multi_multi\n"; + if(poly_1.size() > poly_2.size()){ + return add_multi_multi(poly_1, poly_2); + }else{ + return add_multi_multi(poly_2, poly_1); + } + } + } +} \ No newline at end of file diff --git a/src/product.hpp b/src/product.hpp new file mode 100644 index 00000000..9ccc1af8 --- /dev/null +++ b/src/product.hpp @@ -0,0 +1,67 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +namespace pyqubo { + using indexes = boost::container::small_vector; + class product final { + pyqubo::indexes _indexes; + std::size_t _hash; + + public: + product(const pyqubo::indexes& indexes) noexcept : _indexes(indexes), _hash(create_hash(indexes)) { + ; + } + + std::size_t create_hash(const pyqubo::indexes& indexes){ + std::size_t seed = robin_hood::hash_int(4711); + for(int v: indexes){ + // Combine algorithm is same as boost::hash_conbine except hash function. + // We use robin_hood::hash_int since boost::hash doesn't work well. + seed ^= robin_hood::hash_int(v) + 0x9e3779b9 + (seed << 6) + (seed >> 2); + } + return seed; + } + + product(std::initializer_list init) noexcept : product(pyqubo::indexes(init)) { + ; + } + + std::string to_string() const { + std::stringstream ss; + ss << "["; + for(size_t i = 0; i < _indexes.size(); ++i){ + if(i != 0) ss << ","; + ss << _indexes[i]; + } + ss << "]"; + std::string s = ss.str(); + return s; + } + + const auto& indexes() const noexcept { + return _indexes; + } + + bool equals(pyqubo::product other){ + return this->indexes() == other.indexes(); + } + + friend std::hash; + }; + + +} \ No newline at end of file diff --git a/src/test.cpp b/src/test.cpp new file mode 100644 index 00000000..373e2764 --- /dev/null +++ b/src/test.cpp @@ -0,0 +1,194 @@ +#include +#include +#include +#include +#include +#include +#include + +#include "abstract_syntax_tree.hpp" +#include "expand.hpp" +#include "compiler.hpp" +#include "poly.hpp" +#include "product.hpp" +#include "variables.hpp" + + +std::shared_ptr create_binary(std::string s){ + return std::make_shared(s); +} + +std::shared_ptr create_numeric(double a){ + return std::make_shared(a); +} + +namespace pyqubo { + +void test_poly1(){ + auto variables = pyqubo::variables(); + int index1 = variables.index("x_1"); + int index2 = variables.index("x_1"); + auto p1 = new pyqubo::product({index1}); + auto p2 = new pyqubo::product({index2}); + auto sp1 = pyqubo::poly(create_numeric(2), p1); + auto sp2 = pyqubo::poly(create_numeric(3), p2); + auto sp3 = sp1 + sp2; + + const product cp1 = *p1; + const product cp2 = *p2; + //auto h = std::hash()(cp1); + std::cout << "p1: " << std::to_string(std::hash()(cp1)) << ", p2: " << std::to_string(std::hash()(cp2)) << std::endl; + if(cp1 == cp2){ + std::cout << "p1==p2" << std::endl; + }else{ + std::cout << "p1!=p2" << std::endl; + } + std::cout << sp3.to_string() << std::endl; + +} + +void test_poly2(){ + auto sp1 = pyqubo::poly(create_numeric(2), new pyqubo::product({1})); + auto sp2 = pyqubo::poly(create_numeric(3), new pyqubo::product({2})); + auto sp3 = pyqubo::poly(create_numeric(4), new pyqubo::product({3})); + auto sp0 = pyqubo::poly(); + + auto sp4 = sp1 * sp2; + auto sp5 = sp1 + sp2; + auto sp6 = sp1 + sp3; + auto sp7 = sp5 * sp6; + + std::cout << sp3.to_string() << std::endl; + std::cout << sp4.to_string() << std::endl; + std::cout << sp7.to_string() << std::endl; + + auto sp8 = sp7 + sp5; + std::cout << sp8.to_string() << std::endl; + + auto sp9 = sp0 + sp8; + std::cout << sp9.to_string() << std::endl; +} + +void test_poly3(int n){ + auto variables = pyqubo::variables(); + clock_t t0 = clock(); + auto sp0 = pyqubo::poly(); + for(int i=0; i < n; i++){ + for(int j=0; j < n; j++){ + for(int k=0; k < n; k++){ + int index1 = variables.index("xxx_" + std::to_string(i) + "_" + std::to_string(j)); + int index2 = variables.index("xxx_" + std::to_string(j) + "_" + std::to_string(k)); + auto sp1 = pyqubo::poly(create_numeric(2), new pyqubo::product({index1})); + auto sp2 = pyqubo::poly(create_numeric(2), new pyqubo::product({index2})); + auto sp3 = sp1 * sp2; + sp0 = sp0 + sp3; + } + } + } + std::cout << sp0.to_string() << std::endl; + printf("size: %zu\n", sp0.size()); + clock_t t1 = clock(); + const double expression_time = static_cast(t1 - t0) / CLOCKS_PER_SEC * 1000.0; + printf("time %f\n", expression_time); +} + + +void test_express(int n){ + clock_t t0 = clock(); + auto a = create_binary("a"); + auto b = create_binary("b"); + std::shared_ptr c = std::make_shared(a, b); + + for(int i=0; i < n; i++){ + for(int j=0; j < n; j++){ + for(int k=0; k < n; k++){ + auto bi = create_binary("xxx_" + std::to_string(i) + "_" + std::to_string(j)); + auto bj = create_binary("xxx_" + std::to_string(j) + "_" + std::to_string(k)); + auto coeff = std::make_shared(2); + auto p = coeff*bi*bj; + c = std::make_shared(c, p); + //c->add_child(coeff*bi); + } + } + } + //std::cout << c->to_string() << std::endl; + std::cout << "express done\n"; + clock_t t1 = clock(); + + auto model = pyqubo::compile(c, create_numeric(5)); + std::cout << "compile done\n"; + + clock_t t2 = clock(); + const double expression_time = static_cast(t1 - t0) / CLOCKS_PER_SEC * 1000.0; + const double compile_time = static_cast(t2 - t1) / CLOCKS_PER_SEC * 1000.0; + printf("express time: %f, compile time %f\n", expression_time, compile_time); + const auto [qubo, offset] = model.to_qubo_string({}); + std::cout << "size:" << qubo.size() << "\n"; + for(const auto& [key, value]: qubo){ + //std::cout << std::get<0>(key) << "," << std::get<1>(key) << "," << value << std::endl; + } +} + +void test_express2(int n){ + clock_t t0 = clock(); + + for(int i=0; i < n; i++){ + auto a = create_binary("a"); + pyqubo::compile(a, create_numeric(5)); + } + //std::cout << a->to_string() << std::endl; + + clock_t t1 = clock(); + + clock_t t2 = clock(); + const double expression_time = static_cast(t1 - t0) / CLOCKS_PER_SEC * 1000.0; + const double compile_time = static_cast(t2 - t1) / CLOCKS_PER_SEC * 1000.0; + printf("express time: %f, compile time %f\n", expression_time, compile_time); + //auto qubo = model.to_qubo(); +} + +void test_express3(){ + + auto a = create_binary("a"); + auto b = create_binary("b"); + auto c = create_binary("c"); + + auto H = a * b * c; + auto m = pyqubo::compile(H, create_numeric(5)); + //std::cout << ->to_string() << std::endl; + auto qubo = m.to_qubo_string({}); +} + +void test_express4(){ + + auto a = create_binary("a"); + auto b = create_binary("b"); + auto c = create_binary("c"); + auto n2 = create_numeric(2); + auto n1 = create_numeric(-1); + auto cc = (a + b + n1); + auto subh_var = std::make_shared(a + n2 * b, "my_subh"); + auto with_p = std::make_shared(subh_var, cc * cc, "my_with_penalty"); + + auto H = subh_var + create_numeric(-3); + auto m = pyqubo::compile(H, create_numeric(5)); + //std::cout << ->to_string() << std::endl; + auto qubo = m.to_qubo_string({}); +} + +} + +int main(int argc, char* argv[]){ + if(argc > 1){ + int n = atoi(argv[1]); + printf("hello world\n"); + pyqubo::test_express(n); + } + + //pyqubo::test_poly1(); + //pyqubo::test_poly2(); + //pyqubo::test_express3(); + pyqubo::test_express4(); + + return 0; +} diff --git a/src/variables.hpp b/src/variables.hpp new file mode 100644 index 00000000..8c8e2d96 --- /dev/null +++ b/src/variables.hpp @@ -0,0 +1,134 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +namespace std { + template <> + struct hash { + auto operator()(const pyqubo::product& product) const noexcept { + return product._hash; + } + }; +} + +namespace pyqubo { + using polynomial = robin_hood::unordered_map>; + //using polynomial = std::unordered_map>; + + class variables final { + robin_hood::unordered_map _indexes; + robin_hood::unordered_map _names; + + public: + variables() noexcept : _indexes{}, _names{} { + ; + } + + std::string to_string() const { + std::string s = "variables("; + for(auto [name, index]: _indexes){ + s += name + "->" + std::to_string(index) + "\n"; + } + s += ")"; + return s; + } + + auto index(const std::string& variable_name) noexcept { + const auto [it, emplaced] = _indexes.emplace(variable_name, std::size(_indexes)); + + if (emplaced) { + _names.emplace(it->second, variable_name); + } + + return it->second; + } + + const auto& name(int index) const noexcept { + return _names.find(index)->second; + } + + auto names() const noexcept { + auto result = std::vector(std::size(_names)); + + for (const auto& [index, name] : _names) { + result[index] = name; + } + + return result; + } + }; + + inline auto multiply(const product* product_1, const product* product_2) noexcept { + return new product([&] { + auto result = indexes{}; + + std::set_union(std::begin(product_1->indexes()), std::end(product_1->indexes()), + std::begin(product_2->indexes()), std::end(product_2->indexes()), + std::back_inserter(result)); + + return result; + }()); + } + + inline auto operator*(const product& product_1, const product& product_2) noexcept { + return product([&] { + auto result = indexes{}; + + std::set_union(std::begin(product_1.indexes()), std::end(product_1.indexes()), + std::begin(product_2.indexes()), std::end(product_2.indexes()), + std::back_inserter(result)); + + return result; + }()); + } + + inline bool operator==(const product& product_1, const product& product_2) noexcept { + return product_1.indexes() == product_2.indexes(); + } + + inline auto operator+(const polynomial& polynomial_1, const polynomial& polynomial_2) noexcept { + auto result = polynomial_1; + + for (const auto& [product, coefficient] : polynomial_2) { + const auto [it, emplaced] = result.emplace(product, coefficient); + + if (!emplaced) { + it->second = it->second + coefficient; + } + } + + return result; + } + + inline auto operator*(const polynomial& polynomial_1, const polynomial& polynomial_2) noexcept { + auto result = polynomial{}; + + for (const auto& [product_1, coefficient_1] : polynomial_1) { + for (const auto& [product_2, coefficient_2] : polynomial_2) { + const auto [it, emplaced] = result.emplace(product_1 * product_2, coefficient_1 * coefficient_2); + + if (!emplaced) { + it->second = it->second + coefficient_1 * coefficient_2; + } + } + } + + return result; + } + +} + diff --git a/tests/test_array.py b/tests/test_array.py index f8483d6e..7e8a415b 100644 --- a/tests/test_array.py +++ b/tests/test_array.py @@ -73,40 +73,6 @@ def test_array_transpose(self): [Binary('x[0][2]'), Binary('x[1][2]')]]) self.assertTrue(array.T == expected) - def test_array_add(self): - array1 = Array.create('x', (2, 2), 'BINARY') - array2 = Array.create('y', (2, 2), 'BINARY') - expected = Array([[Binary('x[0][0]') + Binary('y[0][0]'), Binary('x[0][1]') + Binary('y[0][1]')], - [Binary('x[1][0]') + Binary('y[1][0]'), Binary('x[1][1]') + Binary('y[1][1]')]]) - self.assertTrue(array1 + array2 == expected) - - expected2 = Array([[Binary('x[0][0]') + 1, Binary('x[0][1]') + 1], - [Binary('x[1][0]') + 1, Binary('x[1][1]') + 1]]) - self.assertTrue(array1 + 1 == expected2) - self.assertTrue(1 + array1 == expected2) - self.assertTrue(array1 + np.ones((2, 2)) == expected2) - - self.assertRaises(TypeError, lambda: array1 + "str") - array3 = Array.create('z', (2, 3), 'BINARY') - self.assertRaises(ValueError, lambda: array1 + array3) - - def test_array_subtract(self): - array1 = Array.create('x', (2, 2), 'BINARY') - array2 = Array.create('y', (2, 2), 'BINARY') - expected = Array([[Binary('x[0][0]') - Binary('y[0][0]'), Binary('x[0][1]') - Binary('y[0][1]')], - [Binary('x[1][0]') - Binary('y[1][0]'), Binary('x[1][1]') - Binary('y[1][1]')]]) - self.assertTrue(array1 - array2 == expected) - - expected2 = Array([[Binary('x[0][0]') - 1, Binary('x[0][1]') - 1], - [Binary('x[1][0]') - 1, Binary('x[1][1]') - 1]]) - self.assertTrue(array1 - 1 == expected2) - # expected3 = Array([[1 - Binary('x[0][0]'), 1 - Binary('x[0][1]')], - # [1 - Binary('x[1][0]'), 1 - Binary('x[1][1]')]]) - expected3 = Array([[Binary('x[0][0]') * -1 + 1, Binary('x[0][1]') * -1 + 1], - [Binary('x[1][0]') * -1 + 1, Binary('x[1][1]') * -1 + 1]]) - self.assertTrue(1 - array1 == expected3) - self.assertTrue(array1 - np.ones((2, 2)) == expected2) - def test_array_multiply(self): array1 = Array.create('x', (2, 2), 'BINARY') array2 = Array.create('y', (2, 2), 'BINARY') @@ -129,79 +95,6 @@ def test_array_division(self): self.assertRaises(ValueError, lambda: 1 / array1) self.assertRaises(ValueError, lambda: array1 / array1) - def test_array_dot(self): - - # inner product of 1-D arrays - array_a = Array([Binary('a'), Binary('b')]) - array_b = Array([Binary('c'), Binary('d')]) - expected = ((Binary('a') * Binary('c')) + (Binary('b') * Binary('d'))) - self.assertTrue(expected == array_a.dot(array_b)) - - # dot product of 1-D array and N-D array - array_a = Array([[Binary('a'), Binary('b')], [Binary('c'), Binary('d')]]) - array_b = Array([Binary('e'), Binary('f')]) - expected = Array([((Binary('a') * Binary('e')) + (Binary('b') * Binary('f'))), - ((Binary('c') * Binary('e')) + (Binary('d') * Binary('f')))]) - self.assertTrue(expected == array_a.dot(array_b)) - - # both are 2-D arrays - array_a = Array([[Binary('a'), Binary('b')], [Binary('c'), Binary('d')]]) - array_b = Array([[Binary('e'), Binary('f')], [Binary('g'), Binary('h')]]) - expected = Array([[((Binary('a') * Binary('e')) + (Binary('b') * Binary('g'))), - ((Binary('a') * Binary('f')) + (Binary('b') * Binary('h')))], - [((Binary('c') * Binary('e')) + (Binary('d') * Binary('g'))), - ((Binary('c') * Binary('f')) + (Binary('d') * Binary('h')))]]) - self.assertTrue(expected == array_a.dot(array_b)) - - # array_a is a N-D array and array_b is an M-D array (where N, M>=2) - array_a = Array.create('a', shape=(3, 2, 4), vartype='BINARY') - array_b = Array.create('b', shape=(5, 4, 3), vartype='BINARY') - i, j, k, m = (1, 1, 3, 2) - self.assertTrue(array_a.dot(array_b)[i, j, k, m] - == sum(array_a[i, j, :] * array_b[k, :, m])) - - # array_a is 1-D array and array_b is a list - array_a = Array([Binary('a'), Binary('b')]) - array_b = [3, 4] - self.assertTrue((Binary('a') * 3 + (Binary('b') * 4)) == array_a.dot(array_b)) - - # test validation - self.assertRaises(TypeError, lambda: array_a.dot(1)) - - def test_array_matmul(self): - - # the either of the arguments is 1-D array, - array_a = Array([[Binary('a'), Binary('b')], [Binary('c'), Binary('d')]]) - array_b = Array([Binary('e'), Binary('f')]) - expected = Array([((Binary('a') * Binary('e')) + (Binary('b') * Binary('f'))), - ((Binary('c') * Binary('e')) + (Binary('d') * Binary('f')))]) - self.assertTrue(array_a.matmul(array_b) == expected) - - # both arguments are 2-D array - array_a = Array([[Binary('a'), Binary('b')], [Binary('c'), Binary('d')]]) - array_b = Array([[Binary('e'), Binary('f')], [Binary('g'), Binary('h')]]) - expected = Array([[((Binary('a') * Binary('e')) + (Binary('b') * Binary('g'))), - ((Binary('a') * Binary('f')) + (Binary('b') * Binary('h')))], - [((Binary('c') * Binary('e')) + (Binary('d') * Binary('g'))), - ((Binary('c') * Binary('f')) + (Binary('d') * Binary('h')))]]) - self.assertTrue(array_a.matmul(array_b) == expected) - - # either argument is N-D (where N > 2) - array_a = Array.create('a', shape=(2, 2, 3), vartype='BINARY') - array_b = Array.create('b', shape=(3, 2), vartype='BINARY') - self.assertTrue((array_a.matmul(array_b))[0] == array_a[0].matmul(array_b)) - - # array_a is 2-D array and array_b is a list - array_a = Array([[Binary('a'), Binary('b')], [Binary('c'), Binary('d')]]) - array_b = [3, 4] - expected = Array([((Binary('a') * Num(3)) + (Binary('b') * Num(4))), ((Binary('c') * Num(3)) + (Binary('d') * Num(4)))]) - self.assertTrue(expected == array_a.matmul(array_b)) - - # test validation - array_a = Array.create('a', shape=(2, 2, 3), vartype='BINARY') - array_b = Array.create('b', shape=(3, 3, 2), vartype='BINARY') - self.assertRaises(AssertionError, lambda: array_a.matmul(array_b)) - def test_array_reshape(self): array = Array.create('a', shape=(2, 3), vartype='BINARY') reshaped = array.reshape((3, 2)) diff --git a/tests/test_constraint.py b/tests/test_constraint.py index 83f904e4..ba342d25 100644 --- a/tests/test_constraint.py +++ b/tests/test_constraint.py @@ -65,15 +65,6 @@ def test_xor(self): self.assertTrue(model.energy({"a": 0, "b": 0, "c": 1, "aux_xor": 1}, vartype="BINARY") > 0) self.assertTrue(model.energy({"a": 1, "b": 1, "c": 1, "aux_xor": 1}, vartype="BINARY") > 0) - def test_equality(self): - xor1 = XorConst(Binary("a"), Binary("b"), Binary("c"), label="xor") - xor2 = XorConst(Binary("a"), Binary("b"), Binary("c"), label="xor") - xor3 = XorConst(Binary("b"), Binary("c"), Binary("a"), label="xor") - or1 = OrConst(Binary("a"), Binary("b"), Binary("c"), label="xor") - self.assertTrue(xor1 + 1 == xor2 + 1) - self.assertTrue(xor1 == xor2) - self.assertFalse(xor1 == or1) - self.assertFalse(xor1 == xor3) if __name__ == '__main__': unittest.main() diff --git a/tests/test_express.py b/tests/test_express.py index 6c5a0e78..d27d6003 100644 --- a/tests/test_express.py +++ b/tests/test_express.py @@ -4,46 +4,20 @@ class TestExpress(unittest.TestCase): - def test_equality(self): - self.assertEqual(Binary("a"), Binary("a")) - self.assertNotEqual(Binary("a"), Binary("b")) - - exp1 = Binary("a") + 2 - exp2 = Binary("a") + 2 - exp3 = Binary("b") + 2 - self.assertEqual(exp1, exp2) - self.assertNotEqual(exp1, exp3) - - mul1 = Placeholder("p") * Binary("b") - mul2 = Placeholder("p") * Binary("b") - mul3 = Placeholder("p") * Binary("a") - self.assertEqual(mul1, mul2) - self.assertNotEqual(mul1, mul3) - - class CustomPenalty(WithPenalty): - def __init__(self, hamiltonian, penalty, label="label"): - super().__init__(hamiltonian, penalty, label) - - a, b = Binary("a"), Binary("b") - p1 = 1 + CustomPenalty(a+b, a*b) - p2 = 1 + CustomPenalty(a+b, a*b) - p3 = 1 + CustomPenalty(a*b, a+b) - self.assertEqual(p1, p2) - self.assertNotEqual(p1, p3) - - subh1 = SubH(a+b, label="c1") - subh2 = SubH(a+b, label="c1") - subh3 = SubH(a+b, label="c2") - self.assertEqual(subh1, subh2) - self.assertNotEqual(subh1, subh3) - - def compile_check(self, exp, expected_qubo, expected_offset, feed_dict={}): model = exp.compile(strength=5) qubo, offset = model.to_qubo(feed_dict=feed_dict) assert_qubo_equal(qubo, expected_qubo) self.assertEqual(offset, expected_offset) + def test_compile_binary2(self): + a, b = Binary("a"), Binary("b") + exp = a + 2*a + b -1 + # expected_qubo = {('a', 'a'): 1.0, ('a', 'b'): 1.0, ('b', 'b'): 0.0} + expected_qubo = {('a', 'a'): 3.0, ('b', 'b'): 1.0} + expected_offset = -1 + self.compile_check(exp, expected_qubo, expected_offset) + def test_compile_binary(self): a, b = Binary("a"), Binary("b") exp = 1 + a*b + a - 2 @@ -57,6 +31,7 @@ def test_compile_spin(self): expected_qubo = {('a', 'a'): 4.0, ('b', 'b'): -2.0} expected_offset = -2.0 self.compile_check(exp, expected_qubo, expected_offset) + def test_compile_expand_add(self): a, b = Binary("a"), Binary("b") @@ -65,7 +40,8 @@ def test_compile_expand_add(self): expected_qubo = {('a', 'a'): 1.0, ('b', 'b'): -1.0} expected_offset = 0.0 self.compile_check(exp, expected_qubo, expected_offset) - + + def test_compile_div(self): a, b = Binary("a"), Binary("b") exp = a*b / 2 + 1 @@ -74,7 +50,7 @@ def test_compile_div(self): expected_offset = 1.0 q, offset = exp.compile().to_qubo() self.compile_check(exp, expected_qubo, expected_offset) - + def test_compile_power(self): a, b = Binary("a"), Binary("b") exp = (a+b)**3 @@ -82,13 +58,13 @@ def test_compile_power(self): expected_offset = 0.0 q, offset = exp.compile().to_qubo() self.compile_check(exp, expected_qubo, expected_offset) - + def test_compile_neg(self): exp = -Binary("a") expected_qubo = {('a', 'a'): -1.0} expected_offset = 0.0 self.compile_check(exp, expected_qubo, expected_offset) - + def test_compile_placeholder(self): a, b = Binary("a"), Binary("b") p, q, r = Placeholder("p"), Placeholder("q"), Placeholder("r") @@ -98,7 +74,7 @@ def test_compile_placeholder(self): feed_dict={"p": 3, "q": 2, "r": 2} q, offset = exp.compile().to_qubo(feed_dict=feed_dict) self.compile_check(exp, expected_qubo, expected_offset, feed_dict) - + def test_compile_subh(self): a, b = Binary("a"), Binary("b") p = Placeholder("p") @@ -107,16 +83,18 @@ def test_compile_subh(self): expected_offset = 3 feed_dict={"p": 3} self.compile_check(exp, expected_qubo, expected_offset, feed_dict) - + def test_compile_constraint(self): a, b, c = Binary("a"), Binary("b"), Binary("c") exp = Constraint(a*b*c, label="constraint") + #exp=a*b*c # expected_qubo = {('a', '0*1'): -10.0, ('b', '0*1'): -10.0, ('0*1', '0*1'): 15.0, ('a', 'a'): 0.0, ('a', 'b'): 5.0, ('c', '0*1'): 1.0, ('b', 'b'): 0.0, ('c', 'c'): 0.0} + model = exp.compile() expected_qubo = {('a', 'a * b'): -10.0, ('b', 'a * b'): -10.0, ('a * b', 'a * b'): 15.0, ('a', 'b'): 5.0, ('c', 'a * b'): 1.0} expected_offset = 0 self.compile_check(exp, expected_qubo, expected_offset, feed_dict={}) - + def test_compile_with_penalty(self): class CustomPenalty(WithPenalty): def __init__(self, hamiltonian, penalty, label): @@ -129,6 +107,7 @@ def __init__(self, hamiltonian, penalty, label): expected_offset = 0.0 feed_dict={"p": 2} self.compile_check(custom_penalty, expected_qubo, expected_offset, feed_dict) + if __name__ == '__main__': unittest.main() diff --git a/tests/test_integer.py b/tests/test_integer.py index 53edc16f..6af36d47 100644 --- a/tests/test_integer.py +++ b/tests/test_integer.py @@ -29,7 +29,12 @@ def test_one_hot_enc_integer(self): decoded = model.decode_sampleset( sampleset, feed_dict={"s": 10.0}) best = min(decoded, key=lambda x: x.energy) - self.assertTrue(best.subh["a"]==3) + self.assertTrue(best.value(a) == 3) + + # will raise runtime error if the constraint is broken, when evaluating the value of the one-hot-integer object. + worst = max(decoded, key=lambda x: x.energy) + self.assertRaises(RuntimeError, lambda: worst.value(a)) + self.assertTrue(a.value_range == (0, 4)) # expected_q = {('a[0]', 'a[1]'): 20.0, @@ -50,20 +55,20 @@ def test_one_hot_enc_integer(self): # expected_offset = 19 # assert_qubo_equal(q, expected_q) # self.assertTrue(offset == expected_offset) - + def test_one_hot_enc_integer_equal(self): a = OneHotEncInteger("a", (0, 4), strength=Placeholder("s")) b = OneHotEncInteger("b", (0, 4), strength=Placeholder("s")) M = 2.0 - H = (a + b - 5) ** 2 + M*(a.equal_to(3)-1)**2 + H = (a + b - 5) ** 2 + M * (a.equal_to(3) - 1)**2 model = H.compile() q, offset = model.to_qubo(feed_dict={"s": 10.0}) sampleset = dimod.ExactSolver().sample_qubo(q) decoded = model.decode_sampleset( sampleset, feed_dict={"s": 10.0}) best = min(decoded, key=lambda x: x.energy) - self.assertTrue(best.subh["a"]==3) - self.assertTrue(best.subh["b"]==2) + self.assertTrue(best.value(a) == 3) + self.assertTrue(best.value(b) == 2) self.assertTrue(best.subh["a_const"]==0) self.assertTrue(best.subh["b_const"]==0) self.assertEqual(len(best.constraints(only_broken=True)), 0) @@ -113,8 +118,8 @@ def test_log_enc_integer(self): sampleset = dimod.ExactSolver().sample_qubo(q) decoded = model.decode_sampleset(sampleset) best = min(decoded, key=lambda x: x.energy) - self.assertTrue(best.subh["a"] == 2) - self.assertTrue(best.subh["b"] == 3) + self.assertTrue(best.value(a) == 2) + self.assertTrue(best.value(b) == 3) self.assertTrue(a.value_range == (0, 4)) self.assertTrue(b.value_range == (0, 4)) @@ -129,8 +134,8 @@ def test_unary_enc_integer(self): sampleset = dimod.ExactSolver().sample_qubo(q) decoded = model.decode_sampleset(sampleset) best = min(decoded, key=lambda x: x.energy) - self.assertTrue(best.subh["a"] == 1) - self.assertTrue(best.subh["b"] == 2) + self.assertTrue(best.value(a) == 1) + self.assertTrue(best.value(b) == 2) self.assertTrue(a.value_range == (0, 3)) self.assertTrue(b.value_range == (0, 3)) @@ -156,6 +161,7 @@ def test_unary_enc_integer(self): # ('b[1]', 'b[1]'): -7.0, # ('b[2]', 'b[2]'): -7.0} # assert_qubo_equal(q, expected_q) + if __name__ == '__main__': unittest.main() diff --git a/tests/test_logical_constraint.py b/tests/test_logical_constraint.py index 6f1f9af7..af0a51f2 100644 --- a/tests/test_logical_constraint.py +++ b/tests/test_logical_constraint.py @@ -68,15 +68,6 @@ def test_xor(self): self.assertTrue(model.energy({"a": 0, "b": 0, "c": 1, "aux_xor": 1}, vartype="BINARY") > 0) self.assertTrue(model.energy({"a": 1, "b": 1, "c": 1, "aux_xor": 1}, vartype="BINARY") > 0) - def test_equality(self): - xor1 = XorConst(Binary("a"), Binary("b"), Binary("c"), label="xor") - xor2 = XorConst(Binary("a"), Binary("b"), Binary("c"), label="xor") - xor3 = XorConst(Binary("b"), Binary("c"), Binary("a"), label="xor") - or1 = OrConst(Binary("a"), Binary("b"), Binary("c"), label="xor") - self.assertTrue(xor1 + 1 == xor2 + 1) - self.assertTrue(xor1 == xor2) - self.assertFalse(xor1 == or1) - self.assertFalse(xor1 == xor3) if __name__ == '__main__': unittest.main() diff --git a/tests/test_model.py b/tests/test_model.py index 1a89de8d..5a930ed8 100644 --- a/tests/test_model.py +++ b/tests/test_model.py @@ -244,6 +244,25 @@ def test_placeholders(self): self.assertTrue(best_sample.array("S", 2) == 1) self.assertTrue(np.isclose(best_sample.energy, -1.8)) + def test_compiler_strength(self): + a, b, c, d = Binary("a"), Binary("b"), Binary("c"), Binary("d") + H = a * b * c - 2 * c * d - 3 * a * b + model = H.compile(strength=Placeholder("s")) + feed_dict = {"s": 10.0} + qubo, offset = model.to_qubo(feed_dict=feed_dict) + sampler = dimod.ExactSolver() + sampleset = sampler.sample_qubo(qubo) + decoded_samples = model.decode_sampleset(sampleset, feed_dict=feed_dict) + best_sample = min(decoded_samples, key=lambda s: s.energy) + + self.assertTrue(qubo[('a', 'a * b')] == -20) + self.assertTrue(best_sample.value(b) == 1) + self.assertTrue(best_sample.value(a) == 1) + self.assertTrue(best_sample.value(b) == 1) + self.assertTrue(best_sample.value(c) == 1) + self.assertTrue(best_sample.value(d) == 1) + self.assertTrue(np.isclose(best_sample.energy, -4)) + def test_constraint(self): sampler = dimod.ExactSolver() x = Array.create('x',shape=(3),vartype="BINARY") @@ -257,6 +276,7 @@ def test_constraint(self): if sol.energy==1.0: self.assertEqual(sol.subh['C1'], 1.0) + def test_higher_order(self): x = Array.create('x', 5, 'BINARY') exp = x[0]*x[1]*x[2]*x[3] @@ -277,5 +297,7 @@ def test_higher_order(self): e = model.energy(sample, vartype='BINARY') self.assertEqual(e, 10.0) + + if __name__ == '__main__': unittest.main()