New: GraphProblem

Project.toml

@@ -6,18 +6,22 @@ version = "0.1.1"
 [deps]
 BitBasis = "50ba71b6-fa0f-514d-ae9a-0916efc90dcf"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 MLStyle = "d8e11817-5142-5d16-987a-aa16d5891078"
 PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 
 [compat]
 BitBasis = "0.9"
 DocStringExtensions = "0.9"
+Documenter = "1.8.0"
 Graphs = "1"
 InteractiveUtils = "1"
 MLStyle = "0.4"
 PrettyTables = "2"
+Printf = "1.11.0"
 julia = "1.10"
 
 [extras]

src/ProblemReductions.jl

@@ -7,12 +7,14 @@ using BitBasis
 using MLStyle
 using InteractiveUtils: subtypes
 
+using Printf
+
 export @bit_str
 export TruthTable
 export HyperGraph, UnitDiskGraph, GridGraph, PlanarGraph, SimpleGraph
-export @bv_str, StaticElementVector, StaticBitVector, statictrues, staticfalses, onehotv
+export @bv_str, StaticElementVector, StaticBitVector, statictrues, staticfalses, onehotv, hamming_distance
 export num_variables, num_flavors, variables, flavors, weights, set_weights, is_weighted, energy, weight_type, problem_size, configuration_space_size, constraints
-export UnitWeight
+export UnitWeight,ZeroWeight
 
 # models
 export BooleanExpr, Circuit, Assignment, simple_form, CircuitSAT, @circuit, booleans, ¬, ∨, ∧, ⊻, is_literal, is_cnf, is_dnf
@@ -29,10 +31,11 @@ export QUBO
 export Factoring
 export Matching, is_matching
 export MaximalIS
-export PaintShop
+export PaintShop,num_paint_shop_color_switch, paint_shop_coloring_from_config, paint_shop_from_pairs
+export OpenPitMining, is_valid_mining, print_mining
 
 # rules
-export target_problem, AbstractProblem, ConstraintSatisfactionProblem, AbstractReductionResult, reduceto, extract_solution, extract_multiple_solutions, reduce_size
+export target_problem, AbstractProblem, ConstraintSatisfactionProblem, AbstractSatisfiabilityProblem, AbstractReductionResult, reduceto, extract_solution, extract_multiple_solutions, reduce_size
 export LogicGadget, truth_table
 export ReductionSATTo3SAT
 export ReductionCircuitToSpinGlass, ReductionMaxCutToSpinGlass, ReductionSpinGlassToMaxCut, ReductionVertexCoveringToSetCovering, ReductionSatToColoring,

src/bitvector.jl

@@ -169,6 +169,13 @@ function Base.count_ones(x::StaticBitVector)
     sum(v->count_ones(v), x.data)
 end
 
+"""
+    hamming_distance(x::StaticBitVector, y::StaticBitVector)
+
+Calculate the Hamming distance between two static bit vectors.
+"""
+hamming_distance(x::StaticBitVector, y::StaticBitVector) = count_ones(x ⊻ y)
+
 """
 Constructing a static bit vector.
 """

src/models/OpenPitMining.jl

@@ -0,0 +1,189 @@
+"""
+$TYPEDEF
+
+The open pit mining problem.
+This problem can be solved in polynomial time with the pseudoflow algorithm.
+
+Positional arguments
+-------------------------------
+* `rewards` is a matrix of rewards.
+* `blocks` are the locations of the blocks.
+
+Example
+-----------------------------------
+```jldoctest; setup=:(using GenericTensorNetworks)
+julia> rewards =  [-4  -7  -7  -17  -7  -26;
+         0  39  -7   -7  -4    0;
+         0   0   1    8   0    0;
+         0   0   0    0   0    0;
+         0   0   0    0   0    0;
+         0   0   0    0   0    0];
+
+julia> gp = GenericTensorNetwork(OpenPitMining(rewards));
+
+julia> res = solve(gp, SingleConfigMax())[]
+(21.0, ConfigSampler{12, 1, 1}(111000100000))ₜ
+
+julia> is_valid_mining(rewards, res.c.data)
+true
+
+julia> print_mining(rewards, res.c.data)
+     -4      -7      -7     -17      -7     -26 
+      ◼      39      -7      -7      -4       ◼ 
+      ◼       ◼       1       8       ◼       ◼ 
+      ◼       ◼       ◼       ◼       ◼       ◼ 
+      ◼       ◼       ◼       ◼       ◼       ◼ 
+      ◼       ◼       ◼       ◼       ◼       ◼
+```
+
+You will the the mining is printed as green in an colored REPL.
+"""
+struct OpenPitMining{ET} <: ConstraintSatisfactionProblem{ET}
+    rewards::Matrix{ET}
+    blocks::Vector{Tuple{Int,Int}}  # non-zero locations
+    function OpenPitMining(rewards::Matrix{ET}, blocks::Vector{Tuple{Int,Int}}) where ET
+        for (i, j) in blocks
+            checkbounds(rewards, i, j)
+        end
+        new{ET}(rewards, blocks)
+    end
+end
+function OpenPitMining(rewards::Matrix{ET}) where ET
+    # compute block locations
+    blocks = Tuple{Int,Int}[]
+    for i=1:size(rewards, 1), j=i:size(rewards,2)-i+1
+        push!(blocks, (i,j))
+    end
+    OpenPitMining(rewards, blocks)
+end
+
+flavors(::Type{<:OpenPitMining}) = [0, 1]
+function set_weights(c::OpenPitMining, get_weights)
+    rewards = copy(c.rewards)
+    for (w, b) in zip(get_weights, c.blocks)
+        rewards[b...] = w
+    end
+    OpenPitMining(rewards, c.blocks)
+end
+
+"""
+    is_valid_mining(rewards::AbstractMatrix, config)
+
+Return true if `config` (a boolean mask for the feasible region) is a valid mining of `rewards`.
+"""
+function is_valid_mining(rewards::AbstractMatrix, config)
+    blocks = get_blocks(rewards)
+    assign = Dict(zip(blocks, config))
+    for block in blocks
+        if block[1] != 1 && !iszero(assign[block])
+            if iszero(assign[(block[1]-1, block[2]-1)]) ||
+                iszero(assign[(block[1]-1, block[2])]) ||
+                iszero(assign[(block[1]-1, block[2]+1)])
+                return false
+            end
+        end
+    end
+    return true
+end
+function get_blocks(rewards)
+    blocks = Tuple{Int,Int}[]
+    for i=1:size(rewards, 1), j=i:size(rewards,2)-i+1
+        push!(blocks, (i,j))
+    end
+    return blocks
+end
+
+"""
+    print_mining(rewards::AbstractMatrix, config)
+
+Printing the mining solution in a colored REPL.
+"""
+function print_mining(rewards::AbstractMatrix{T}, config) where T
+    k = 0
+    for i=1:size(rewards, 1)
+        for j=1:size(rewards, 2)
+            if j >= i && j <= size(rewards,2)-i+1
+                k += 1
+                if T <: Integer
+                    str = Printf.@sprintf " %6i " rewards[i,j]
+                else
+                    str = Printf.@sprintf " %6.2F " rewards[i,j]
+                end
+                if iszero(config[k])
+                    printstyled(str; color = :red)
+                else
+                    printstyled(str; color = :green)
+                end
+            else
+                str = Printf.@sprintf " %6s " "◼"
+                printstyled(str; color = :black)
+            end
+        end
+        println()
+    end
+end
+
+function _open_pit_mining_branching!(rewards::AbstractMatrix{T}, mask::AbstractMatrix{Bool}, setmask::AbstractMatrix{Bool}, idx::Int) where T
+    # find next
+    idx < 1 && return zero(T)
+    i, j = divrem(idx-1, size(mask, 2)) .+ 1 # row-wise!
+    while i > j || size(mask, 1)-i+1 < j || setmask[i, j]  # skip not allowed or already decided
+        idx -= 1
+        idx < 1 && return zero(T)
+        i, j = divrem(idx-1, size(mask, 2)) .+ 1 # row-wise!
+    end
+    if rewards[i, j] < 0 # do not mine!
+        setmask[i, j] = true
+        return _open_pit_mining_branching!(rewards, mask, setmask, idx-1)
+    else
+        _mask = copy(mask)
+        _setmask = copy(setmask)
+        # CASE 1: try mine current block
+        # set mask
+        reward0 = set_recur!(mask, setmask, rewards, i, j)
+        reward1 = _open_pit_mining_branching!(rewards, mask, setmask, idx-1) + reward0
+
+        # CASE 1: try do not mine current block
+        # unset mask
+        _setmask[i, j] = true
+        reward2 = _open_pit_mining_branching!(rewards, _mask, _setmask, idx-1)
+
+        # choose the right branch
+        if reward2 > reward1
+            copyto!(mask, _mask)
+            copyto!(setmask, _setmask)
+            return reward2
+        else
+            return reward1
+        end
+    end
+end
+
+function set_recur!(mask, setmask, rewards::AbstractMatrix{T}, i, j) where T
+    reward = zero(T)
+    for k=1:i
+        start = max(1, j-(i-k))
+        stop = min(size(mask, 2), j+(i-k))
+        @inbounds for l=start:stop
+            if !setmask[k,l]
+                mask[k,l] = true
+                setmask[k,l] = true
+                reward += rewards[k,l]
+            end
+        end
+    end
+    return reward
+end
+
+"""
+    open_pit_mining_branching(rewards::AbstractMatrix)
+
+Solve the open pit mining problem with the naive branching algorithm.
+NOTE: open pit mining can be solved in polynomial time, but this one runs in exponential time.
+"""
+function open_pit_mining_branching(rewards::AbstractMatrix{T}) where T
+    idx = length(rewards)
+    mask = falses(size(rewards))
+    rewards = _open_pit_mining_branching!(rewards, mask, falses(size(rewards)), idx)
+    return rewards, mask
+end

src/models/Paintshop.jl

@@ -50,12 +50,24 @@ struct PaintShop{LT} <: ConstraintSatisfactionProblem{Int}
         new{eltype(sequence)}(sequence, isfirst)
     end
 end
+function paint_shop_from_pairs(pairs::AbstractVector{Tuple{Int,Int}})
+    n = length(pairs)
+    @assert sort!(vcat(collect.(pairs)...)) == collect(1:2n)
+    sequence = zeros(Int, 2*n)
+    @inbounds for i=1:n
+        sequence[pairs[i]] .= i
+    end
+    return PaintShop(sequence)
+end
 
 variables(gp::PaintShop) = unique(gp.sequence)
 flavors(::Type{<:PaintShop}) = [0, 1]
 problem_size(c::PaintShop) = (; sequence_length=length(c.sequence))
 Base.:(==)(a::PaintShop, b::PaintShop) = a.sequence == b.sequence && a.isfirst == b.isfirst
 
+# weights interface
+set_weights(c::PaintShop, weights) = c # constant UnitWeight
+
 # constraints interface
 function energy_terms(c::PaintShop)
     # constraints on alphabets with the same color
@@ -72,6 +84,15 @@ end
 
 @nohard_constraints PaintShop
 
+"""
+    num_paint_shop_color_switch(sequence::AbstractVector, coloring)
+
+Returns the number of color switches.
+"""
+function num_paint_shop_color_switch(sequence::AbstractVector, coloring)
+    return count(i->coloring[i] != coloring[i+1], 1:length(sequence)-1)
+end
+
 """
     paint_shop_coloring_from_config(p::PaintShop, config)
 

src/models/models.jl

@@ -200,6 +200,16 @@ end
 Base.getindex(::UnitWeight, i) = 1
 Base.size(w::UnitWeight) = (w.n,)
 
+"""
+    UnitWeight <: AbstractVector{Int}
+
+The unit weight vector of length `n`.
+"""
+struct ZeroWeight <: AbstractVector{Int}
+    n::Int
+end
+Base.getindex(::ZeroWeight, i) = 0
+Base.size(w::ZeroWeight) = (w.n,)
 """
     energy_terms(problem::AbstractProblem) -> Vector{LocalConstraint}
 
@@ -248,3 +258,4 @@ include("Factoring.jl")
 include("Matching.jl")
 include("MaximalIS.jl")
 include("Paintshop.jl")
+include("OpenPitMining.jl")