Pages

Tuesday, 22 November 2016

Constraint Programming with python-constraint

I have recently been looking at constraint programming and I thought I would post some of my code as I go. This is using the python-constraint package.

Doing things with python-constraint is pretty easy. To constrain a variable to have a certain value, we can use InSetConstraint. To make sure all the variables in a set are different, we use the AllDifferentConstraint. These are enough to solve sudoku.

Sudoku

I don't think Sudoku needs any introduction. The constraints are fairly easy to specify: each row must be all different, each col must be all different and each 3x3 block must be all different. Also we have to ensure that the given starting state is used, for this we use equality constraints. We use 1 variable for each number in the sudoku square, so 81 variables all up.

from constraint import *
problem = Problem()
problem.addVariables(range(9*9), range(1,10))

# specify the given starting sudoku problem here, dots are blank squares
sudoku = '\
.42986.1.\
....2...5\
......82.\
4.1..793.\
.683.254.\
.758..1.2\
.39......\
7...4....\
.1.67325.'

# Make sure our solution matches the specified starting thing
for i,char in enumerate(sudoku):
    if char != '.':  problem.addConstraint(InSetConstraint([int(char)]),[i])
    
# add row and col constraints
for i in range(9): 
    problem.addConstraint(AllDifferentConstraint(),range(9*i,9*i+9))
for i in range(9): 
    problem.addConstraint(AllDifferentConstraint(),range(i,81,9))
# add block constraints
for i in [0,3,6,27,30,33,54,57,60]: 
    ind = []
    for j in range(3):
        for k in range(3):
            ind.append(i + j*9 + k)
    problem.addConstraint(AllDifferentConstraint(),ind)

# pretty print the output
a = problem.getSolution()
for i in range(9):
    for j in range(9): print a[i*9+j],
    print ''

The Golomb Ruler

A bit of info on Golomb Rulers: here. If you know how to improve the code below, leave a comment. Golomb rulers are described by their order, which is the number of marks on their edge. The length of a Golomb ruler is the distance between the outer two marks and represents the longest distance it can measure.

There is no requirement that a Golomb ruler measures all distances up to their length. However, if a ruler does measure all distances, it is classified as a perfect Golomb ruler. There are no perfect Golomb rulers above order 4.

Finally, a Golomb ruler is described as optimal if no shorter ruler of the same order exists. The code below first looks for perfect rulers, then increments length each time until an optimal ruler is found. It starts to slow down by order = 7, it takes a few minutes.

This problem is harder than the sudoku problem, we have to use a set of auxiliary variables to represent the pairwise differences, then apply the AllDifferentConstraint to the auxiliary variables.

from constraint import *
from itertools import combinations

order = 7 # the number of marks on the Golomb ruler
diffs = [] # the set of pairwise differences between marks
for i in combinations(range(order),2): diffs.append(i)
length = len(diffs)

def lessthan(a,b): return a < b
def diffeq(a,b,res): return (b-a)==res     

while True:
    print 'trying length: '+str(length)
    problem = Problem()
    problem.addVariables(range(order), range(0,length+1))
    problem.addVariables(diffs,    range(1,length+1))

    # make sure the first mark is 0 and last is length
    problem.addConstraint(InSetConstraint([0]),[0])
    problem.addConstraint(InSetConstraint([length]),[order-1])

    # ensure that the marks are in increasing order
    # make sure all possible pairwise relations are constrained
    for i in diffs: problem.addConstraint(lessthan,[i[0],i[1]])

    # ensure that the differences reflect the marks
    for i in diffs: problem.addConstraint(diffeq,[i[0],i[1],i])
    
    # all the differences must be different    
    problem.addConstraint(AllDifferentConstraint(),diffs)
   
    # pretty print the output
    solutions = problem.getSolutions()
    for ruler in solutions:
        for mark in range(order): print ruler[mark],
        print ''
    if len(solutions) > 0: break # if we find some solutions, quit
    length = length + 1

example output for order = 7, note the shortest ruler is length 25, though a perfect ruler would be length 21:

trying length: 21
trying length: 22
trying length: 23
trying length: 24
trying length: 25
0 4 9 15 22 23 25
0 3 4 12 18 23 25
0 2 3 10 16 21 25
0 2 5 14 18 24 25
0 2 6 9 14 24 25
0 2 7 13 21 22 25
0 2 7 15 21 24 25
0 1 4 10 18 23 25
0 1 7 11 20 23 25
0 1 11 16 19 23 25

Magic Squares

Magic Squares seem to be very common examples when doing constraint programming, so I may as well include the code for them. Info on Magic Squares can be found here. This code is quite slow, I thought it would be much faster. This may be because python-constraint is just kind of slow, or maybe the problem is harder than I thought.

from constraint import *

problem = Problem()
size = 4
problem.addVariables(range(size*size), range(1,size*size+1))
total = size*(size*size + 1)/2

# make sure rows add to total
for i in range(size): 
    problem.addConstraint(ExactSumConstraint(total),range(i*size,i*size+size))
for i in range(size): 
    problem.addConstraint(ExactSumConstraint(total),range(i,size*size,size))
# add the diagonal constraints
diag1,diag2 = [],[]
for i in range(size):
    diag1.append(i*size + i)
    diag2.append(i*size + (size-i-1))

problem.addConstraint(ExactSumConstraint(total),diag1)
problem.addConstraint(ExactSumConstraint(total),diag2)
problem.addConstraint(AllDifferentConstraint())
   
# pretty print the output
solution = problem.getSolution()
for i in range(size):
    for j in range(size):
        print solution[i*size + j],
    print ''

Monday, 14 November 2016

Terraforming Venus

I thought I would look at some of the figures involved in terraforming Venus, discussions of which can be found all over the place e.g. here. The numbers from this page come from here. Venus is the planet next to earth as you go towards the Sun. It has a year of 224.7 Earth days, and a day length of 243 Earth days (longer than the year!). It also rotates in the opposite direction to most other planets, likely due to some large collision early in its life.

I'll quote a little bit from wikipedia: Venus has an extremely dense atmosphere composed of 96.5% carbon dioxide, 3.5% nitrogen, and traces of other gases, most notably sulfur dioxide. The mass of its atmosphere is 93 times that of Earth's, whereas the pressure at its surface is about 92 times that at Earth's—a pressure equivalent to that at a depth of nearly 1 kilometre under Earth's oceans. The density at the surface is 65 kg/m3, 6.5% that of water or 50 times as dense as Earth's atmosphere at 20 °C at sea level. The CO2-rich atmosphere generates the strongest greenhouse effect in the Solar System, creating surface temperatures of at least 735 K (462 °C).

Terraforming Venus would involve removing most of the sulfuric acid and carbon dioxide, needless to say this is a monumental task that we don't currently have the technology for. Some of the truly staggering figures:

Total Mass of Venus Atmosphere = 4.80E+20 kg
Percent Atmosphere CO2 = 96.50%
Total Mass of CO2 = 4.63E+20 kg
For comparison:
Total Mass of Earth's Atmosphere = 5.10E+18 kg
Mass Ratio Venus to Earth = 94.118

Utilizing the Bosch reaction (CO2 + 2H2 -> C + 2H2O) , combining hydrogen with carbon dioxide to make carbon graphite and water could be used to remove Venus' carbon dioxide. This reaction requires the introduction of iron, cobalt or nickel as a catalyst and requires a temperature level of 530-730 degrees Celsius. Since the reaction releases some heat, it could be self sustaining if you could supply enough hydrogen. A problem with this is that that the production of elemental carbon tends to foul the catalyst's surface, which is detrimental to the reaction's efficiency [ref].

Molecular Weight of CO2 = 44
Molecular Weight of 2*H2 = 4
Total Initial Molecular Weight = 48
Molecular Weight of C (graphite) = 12
Molecular Weight of 2*H2O = 36
Total Final Molecular Weight = 48

You could use water from e.g. comets to get hydrogen for the reaction. You would have to use Electrolysis to split the water and initiate the reaction. Alternatively you could reuse the water output from the Bosch reaction , recover the Hydrogen and continue. But this way you wouldn't have oceans on Venus at the end.

If you were to provide hydrogen for the reaction (4E19 kg), this massive reaction would create an ocean nearly as large as 1/4 of the Earth’s ocean:

Ratio of 2*H2O to CO2 (36 / 44) = 0.818
Resultant Mass of H2O = 3.79E+20 kg
Density of H2O = 1.000 g / cm^3
= 1,000,000.000 g / M^3
= 1,000.000 kg / M^3
Volume of Resultant H2O = 3.79E+17 M^3
= 3.79E+08 kM^3
Area of Venus Surface = 4.602E+08 kM^2
Average Depth of H2O = 0.824 kM
Total Volume of Earth's Oceans = 1.300E+09 kM^3
Average Depth of Earth's Oceans = 3.682 kM

It would also result in the deposition of a layer of graphite with an average thickness over the entire surface of Venus roughly equal to a 40 story building:

Ratio of graphite to CO2 (12 / 44) = 0.273
Resultant Mass of graphite = 1.26E+20 kg
Density of C (graphite) = 2.230 g / cm^3
= 2,230,000.000 g / M^3
= 2,230.000 kg / M^3
Volume of Resultant C (graphite) = 5.66E+16 M^3
= 5.66E+07 kM^3
Area of Venus Surface = 4.602E+08 kM^2
Average Depth of C (graphite) = 0.123 kM

So after converting all the carbon dioxide to graphite, you would be left with ~120 meter thick blanket of carbon over the planet. Though this may not be terrible, carbon makes a pretty good base for soils. There would still be a great deal of atmospheric pressure due to all the oxygen though, and a spark would result in a mighty conflagration that would undo all your work. A better idea may be to trap the carbon dioxide as carbonate rocks, but for this you would need large quantities of calcium and/or magnesium (About 8×10^20 kg of calcium or 5×10^20 kg of magnesium would be required). This would have the advantage of removing a lot of the oxygen, which would lower the atmospheric pressure. Unfortunately Magnesium carbonate begins to decompose at 350 degrees C, so Calcium might be a better bet. You would also need to remove all the sulphuric acid since it would dissolve all your carbonate rocks, using up the sulphuric acid but producing more carbon dioxide.

Of course you have to power all this as well, electrolysis doesn't come cheap. Solar panels can work well on Venus because the sunlight is 4 times stronger. Perhaps some genetically engineered plant...