Penalty model lp

.circleci/config.yml

@@ -52,6 +52,12 @@ jobs:
             . env/bin/activate
             pip install penaltymodel_mip/
 
+      - run: &install-lp-template
+          name: install penaltymodel-lp
+          command: |
+            . env/bin/activate
+            pip install penaltymodel_lp/
+
       - run: &core-tests-template
           name: core tests
           command: |
@@ -70,6 +76,12 @@ jobs:
             . env/bin/activate
             coverage run -a -m unittest discover -s penaltymodel_maxgap/
 
+      - run: &lp-tests-template
+          name: lp tests
+          command: |
+            . env/bin/activate
+            coverage run -a -m unittest discover -s penaltymodel_lp/
+
       - run: &mip-tests-template
           name: mip tests
           command: |
@@ -172,6 +184,8 @@ jobs:
 
       - run: *install-mip-template
 
+      - run: *install-lp-template
+
       - run: *core-tests-template
 
       - run: *cache-tests-template
@@ -180,6 +194,8 @@ jobs:
 
       - run: *mip-tests-template
 
+      - run: *lp-tests-template
+
       - run: *integration-tests-template
 
   test-osx-3.6:
@@ -350,6 +366,40 @@ jobs:
 
       - run: *twine-template
 
+  deploy-lp:
+    docker:
+      - image: circleci/python:3.6-jessie
+
+    working_directory: ~/repo
+
+    steps:
+      - checkout
+
+      - run: *create-virtualenv-template
+
+      - run: *install-core-template
+
+      - run: *install-lp-template
+
+      - run:
+          name: verify version matches tag
+          command: |
+            . env/bin/activate
+            [[ "$(pip show penaltymodel-lp 2>/dev/null | grep Version)" == "Version: $(echo "$CIRCLE_TAG" | sed -E 's/lp-([0-9]+\.[0-9]+\.[0-9]+(\.dev[0-9]*)?)$/\1/')" ]]
+
+      - run: *pypirc-file-template
+
+      - run:
+          name: build bdist and sdist
+          command: |
+            . env/bin/activate
+            cd penaltymodel_lp
+            python setup.py sdist
+            python setup.py bdist_wheel
+            mv dist/* ../dist/
+
+      - run: *twine-template
+
 workflows:
   version: 2
   test:
@@ -388,3 +438,9 @@ workflows:
               only: /^mip-[0-9]+(\.[0-9]+)*(\.dev([0-9]+)?)?$/
             branches:
               ignore: /.*/
+      - deploy-lp:
+          filters:
+            tags:
+              only: /^lp-[0-9]+(\.[0-9]+)*(\.dev([0-9]+)?)?$/
+            branches:
+              ignore: /.*/

appveyor.yml

@@ -19,6 +19,7 @@ install:
   - "%PYTHON%\\python.exe -m pip install penaltymodel_cache\\"
   - "%PYTHON%\\python.exe -m pip install penaltymodel_maxgap\\"
   - "%PYTHON%\\python.exe -m pip install penaltymodel_mip\\"
+  - "%PYTHON%\\python.exe -m pip install penaltymodel_lp\\"
 
 build: off
 
@@ -27,4 +28,5 @@ test_script:
   - "%PYTHON%\\python.exe -m unittest discover -s penaltymodel_cache"
   - "%PYTHON%\\python.exe -m unittest discover -s penaltymodel_maxgap"
   - "%PYTHON%\\python.exe -m unittest discover -s penaltymodel_mip"
+  - "%PYTHON%\\python.exe -m unittest discover -s penaltymodel_lp"
   - "%PYTHON%\\python.exe -m unittest discover tests/"

penaltymodel_lp/penaltymodel/__init__.py

@@ -0,0 +1,2 @@
+import pkgutil
+__path__ = pkgutil.extend_path(__path__, __name__)

penaltymodel_lp/penaltymodel/lp/__init__.py

@@ -0,0 +1,5 @@
+from penaltymodel.lp.generation import *
+import penaltymodel.lp.generation
+
+from penaltymodel.lp.interface import *
+import penaltymodel.lp.interface

penaltymodel_lp/penaltymodel/lp/generation.py

@@ -0,0 +1,166 @@
+import dimod
+from itertools import product
+import numpy as np
+from scipy.optimize import linprog
+
+#TODO: put these values in a common penaltymodel folder
+MIN_LINEAR_BIAS = -2
+MAX_LINEAR_BIAS = 2
+MIN_QUADRATIC_BIAS = -1
+MAX_QUADRATIC_BIAS = 1
+DEFAULT_GAP = 2
+
+
+def get_item(dictionary, tuple_key, default_value):
+    """Grab values from a dictionary using an unordered tuple as a key.
+
+    Dictionary should not contain None, 0, or False as dictionary values.
+
+    Args:
+        dictionary: Dictionary that uses two-element tuple as keys
+        tuple_key: Unordered tuple of two elements
+        default_value: Value that is returned when the tuple_key is not found in the dictionary
+    """
+    u, v = tuple_key
+
+    # Grab tuple-values from dictionary
+    tuple1 = dictionary.get((u, v), None)
+    tuple2 = dictionary.get((v, u), None)
+
+    # Return the first value that is not {None, 0, False}
+    return tuple1 or tuple2 or default_value
+
+
+def _get_lp_matrix(spin_states, nodes, edges, offset_weight, gap_weight):
+    """Creates an linear programming matrix based on the spin states, graph, and scalars provided.
+    LP matrix:
+        [spin_states, corresponding states of edges, offset_weight, gap_weight]
+
+    Args:
+        spin_states: Numpy array of spin states
+        nodes: Iterable
+        edges: Iterable of tuples
+        offset_weight: Numpy 1-D array or number
+        gap_weight: Numpy 1-D array or a number
+    """
+    if len(spin_states) == 0:
+        return None
+
+    # Set up an empty matrix
+    n_states = len(spin_states)
+    m_linear = len(nodes)
+    m_quadratic = len(edges)
+    matrix = np.empty((n_states, m_linear + m_quadratic + 2))   # +2 columns for offset and gap
+
+    # Populate linear terms (i.e. spin states)
+    if spin_states.ndim == 1:
+        spin_states = np.expand_dims(spin_states, 1)
+    matrix[:, :m_linear] = spin_states
+
+    # Populate quadratic terms
+    node_indices = dict(zip(nodes, range(m_linear)))
+    for j, (u, v) in enumerate(edges):
+        u_ind = node_indices[u]
+        v_ind = node_indices[v]
+        matrix[:, j + m_linear] = np.multiply(matrix[:, u_ind], matrix[:, v_ind])
+
+    # Populate offset and gap columns, respectively
+    matrix[:, -2] = offset_weight
+    matrix[:, -1] = gap_weight
+    return matrix
+
+
+#TODO: check table is not empty (perhaps this check should be in bqm.stitch or as a common
+# penaltymodel check)
+#TODO: Check if len(table.keys()[0]) == len(decision_variables)
+
+#TODO: Adding min_gap into code
+#TODO: Checking min_gap < max_gap (occurs at bqm.stitch or penaltymodel?)
+def generate_bqm(graph, table, decision_variables,
+                 linear_energy_ranges=None, quadratic_energy_ranges=None):
+
+    # Check for auxiliary variables in the graph
+    if len(graph) != len(decision_variables):
+        raise ValueError('Penaltymodel-lp does not handle problems with auxiliary variables')
+
+    if not linear_energy_ranges:
+        linear_energy_ranges = {}
+
+    if not quadratic_energy_ranges:
+        quadratic_energy_ranges = {}
+
+    # Simplify graph naming
+    # Note: nodes' and edges' order determine the column order of the LP
+    nodes = decision_variables
+    edges = graph.edges
+
+    # Set variable names for lengths
+    m_linear = len(nodes)                   # Number of linear biases
+    m_quadratic = len(edges)                # Number of quadratic biases
+    n_noted = len(table)                    # Number of spin combinations specified in the table
+    n_unnoted = 2**m_linear - n_noted       # Number of spin combinations that were not specified
+
+    # Linear programming matrix for spin states specified by 'table'
+    noted_states = table.keys() if isinstance(table, dict) else table
+    noted_states = list(noted_states)
+    noted_matrix = _get_lp_matrix(np.asarray(noted_states), nodes, edges, 1, 0)
+
+    # Linear programming matrix for spins states that were not specified by 'table'
+    # Note: When spin
+    spin_states = product([-1, 1], repeat=m_linear) if m_linear > 1 else [-1, 1]
+    unnoted_states = [state for state in spin_states if state not in noted_states]
+    unnoted_matrix = _get_lp_matrix(np.asarray(unnoted_states), nodes, edges, 1, -1)
+    if unnoted_matrix is not None:
+        unnoted_matrix *= -1   # Taking negative in order to flip the inequality
+
+    # Constraints
+    if isinstance(table, dict):
+        noted_bound = np.asarray([table[state] for state in noted_states])
+        unnoted_bound = np.full((n_unnoted, 1), -1 * min(table.values()))  # -1 for flipped inequality
+    else:
+        noted_bound = np.zeros((n_noted, 1))
+        unnoted_bound = np.zeros((n_unnoted, 1))
+
+    # Bounds
+    # Note: max_gap's max(..or..) is to support python2.7. TODO: ideally, use max([..], default)
+    linear_range = (MIN_LINEAR_BIAS, MAX_LINEAR_BIAS)
+    quadratic_range = (MIN_QUADRATIC_BIAS, MAX_QUADRATIC_BIAS)
+
+    bounds = [linear_energy_ranges.get(node, linear_range) for node in nodes]
+    bounds += [get_item(quadratic_energy_ranges, edge, quadratic_range) for edge in edges]
+
+    # Note: Since ising has {-1, 1}, the largest possible gap is [-largest_bias, largest_bias],
+    #   hence that 2 * sum(largest_biases)
+    max_gap = 2 * sum(max(abs(lbound), abs(ubound)) for lbound, ubound in bounds)
+    bounds.append((None, None))     # Bound for offset
+    bounds.append((0, max_gap))     # Bound for gap. TODO: bound with min_gap as well
+
+    # Cost function
+    cost_weights = np.zeros((1, m_linear + m_quadratic + 2))
+    cost_weights[0, -1] = -1     # Only interested in maximizing the gap
+
+    # Returns a Scipy OptimizeResult
+    result = linprog(cost_weights.flatten(), A_eq=noted_matrix, b_eq=noted_bound,
+                     A_ub=unnoted_matrix, b_ub=unnoted_bound, bounds=bounds)
+
+    #TODO: propagate scipy.optimize.linprog's error message?
+    if not result.success:
+        raise ValueError('Penaltymodel-lp is unable to find a solution.')
+
+    # Split result
+    x = result.x
+    h = x[:m_linear]
+    j = x[m_linear:-2]
+    offset = x[-2]
+    gap = x[-1]
+
+    if gap <= 0:
+        raise ValueError('Penaltymodel-lp is unable to find a solution.')
+
+    # Create BQM
+    bqm = dimod.BinaryQuadraticModel.empty(dimod.SPIN)
+    bqm.add_variables_from((v, bias) for v, bias in zip(nodes, h))
+    bqm.add_interactions_from((u, v, bias) for (u, v), bias in zip(edges, j))
+    bqm.add_offset(offset)
+
+    return bqm, gap

penaltymodel_lp/penaltymodel/lp/interface.py

@@ -0,0 +1,48 @@
+import dimod
+
+from six import iteritems
+
+import penaltymodel.core as pm
+
+from penaltymodel.lp.generation import generate_bqm
+
+__all__ = ['get_penalty_model']
+
+
+@pm.penaltymodel_factory(100)
+def get_penalty_model(specification):
+    """Factory function for penaltymodel-lp.
+
+    Args:
+        specification (penaltymodel.Specification): The specification
+            for the desired penalty model.
+
+    Returns:
+        :class:`penaltymodel.PenaltyModel`: Penalty model with the given specification.
+
+    Raises:
+        :class:`penaltymodel.ImpossiblePenaltyModel`: If the penalty cannot be built.
+
+    Parameters:
+        priority (int): -100
+
+    """
+    # check that the feasible_configurations are spin
+    feasible_configurations = specification.feasible_configurations
+    if specification.vartype is dimod.BINARY:
+        feasible_configurations = {tuple(2 * v - 1 for v in config): en
+                                   for config, en in iteritems(feasible_configurations)}
+
+    # convert ising_quadratic_ranges to the form we expect
+    ising_quadratic_ranges = specification.ising_quadratic_ranges
+    quadratic_ranges = {(u, v): ising_quadratic_ranges[u][v] for u, v in specification.graph.edges}
+
+    try:
+        bqm, gap = generate_bqm(specification.graph, feasible_configurations,
+                                specification.decision_variables,
+                                linear_energy_ranges=specification.ising_linear_ranges,
+                                quadratic_energy_ranges=quadratic_ranges)
+    except ValueError:
+        raise pm.exceptions.FactoryException("Specification is for too large of a model")
+
+    return pm.PenaltyModel.from_specification(specification, bqm, gap, 0.0)

penaltymodel_lp/requirements.txt

@@ -0,0 +1,4 @@
+dimod==0.7.2
+dwavebinarycsp==0.0.6
+numpy==1.15.3
+scipy==1.1.0

penaltymodel_lp/setup.py

@@ -0,0 +1,21 @@
+from setuptools import setup
+
+FACTORY_ENTRYPOINT = 'penaltymodel_factory'
+
+install_requires = ['dimod>=0.6.0,<0.8.0',
+                    'penaltymodel>=0.15.0,<0.16.0',
+                    'scipy>=0.15.0,<1.1.0',
+                    'numpy>=0.0.0,<1.16.0'
+                    ]
+
+packages = ['penaltymodel',
+            'penaltymodel.lp',
+            ]
+
+setup(
+    name="penaltymodel-lp",
+    install_requires=install_requires,
+    entry_points={
+        FACTORY_ENTRYPOINT: ['lp = penaltymodel.lp:get_penalty_model']
+    }
+)