Z3による汎用大域的最適化 ref: http://qiita.com/kotauchisunsun/items/90bce940a1c08ca1b7ae

file0.txt

Z3 is a theorem prover from Microsoft Research. It is licensed under the MIT license.

If you are not familiar with Z3, you can start here.

Z3 can be built using Visual Studio, a Makefile or using CMake. It provides bindings for several programming languages.

See the release notes for notes on various stable releases of Z3.

file1.py

from z3 import *

x = Real("x")
y = Real("y")
c = Real("c")

s = Solver()

s.add(
    x > 0,
    y > 0,
    c > 0,
    c * c == 2.0,
    x*x + y*y == 1,
    y == -x + c
)

print s.check()
print s.model()

file2.txt

$ python z3_test.py 
> sat
> [x = 0.7071067811?, c = 1.4142135623?, y = 0.7071067811?]

file3.txt

c > 0
c * c == 2.0

file4.py

import sys
from random import gauss,seed

seed(0)
n = int(sys.argv[1])

for i in range(n):
    print i,gauss(2.0 * i + 3.0, 0.1)

file5.txt

$python data.py 50 > data.csv

file6.py

from z3 import *

if __name__ == "__main__":
    data = map(lambda x:map(float,x.split()),open("data.csv")) # 1.

    s = Solver()
    a = Real("a")
    b = Real("b")
    delta = Real("delta")

    deltas = []
    for i,(x,y) in enumerate(data):
        d = Real("d%d"%i)
        deltas.append(d)
        s.add(
            d == y - ( a * x + b )  # 2.
        )
        s.check()

    s.add(delta == sum(d * d for d in deltas)) # 3.

    print s.check() # 4.
    print s.model()

    max_delta = float(s.model()[delta].as_decimal(15)[:-1]) # 5.
    min_delta = 0                                           # 6.

    while 1:
        tmp_delta = (max_delta + min_delta) / 2.0           #7.
        s.push()

        s.add( delta < tmp_delta )   # 8.

        if s.check() == sat: # 9.
            print "sat:",tmp_delta,min_delta,max_delta
            s.push()
            max_delta = tmp_delta
        else: # 10.
            print "unsat:",tmp_delta,min_delta,max_delta
            s.pop()
            min_delta = tmp_delta

        if max_delta - min_delta < 1.0e-6: # 11.
            break

    print s.check()
    model = s.model()

    print delta,model[delta].as_decimal(15)
    print a,model[a].as_decimal(15)
    print b,model[b].as_decimal(15)

file7.txt

delta 0.604337347396090?
a 2.001667022705078?
b 2.963314133483435?

file8.py

from z3 import *

if __name__ == "__main__":
    data = map(lambda x:map(float,x.split()),open("data.csv")) # 1.

    s = Solver()
    a = Real("a")
    b = Real("b")
    epsilon = Real("epsilon")

    s.add(epsilon >= 0.0) # 2.

    for i,(x,y) in enumerate(data):
        s.add(
            y - epsilon <=  ( a * x + b ) , (a*x+b) <= y + epsilon # 3.
        )
        s.check() 

    print s.check() # 4.
    print s.model()

    max_epsilon = float(s.model()[epsilon].as_decimal(15)[:-1]) # 5.
    min_epsilon = 0 # 6.

    while 1:
        tmp_epsilon = (max_epsilon + min_epsilon) / 2.0 # 7.
        s.push()

        s.add( epsilon < tmp_epsilon ) # 8.

        if s.check() == sat: # 9.
            print "sat:",tmp_epsilon,min_epsilon,max_epsilon
            s.push()
            max_epsilon = tmp_epsilon
        else: # 10.
            print "unsat:",tmp_epsilon,min_epsilon,max_epsilon
            s.pop()
            min_epsilon = tmp_epsilon

        if max_epsilon - min_epsilon < 1.0e-6: # 11.
            break

    print s.check()
    model = s.model()

    print epsilon,model[epsilon].as_decimal(15)
    print a,model[a].as_decimal(15)
    print b,model[b].as_decimal(15)

file9.txt

epsilon 0.223683885406060?
a 2.000628752115151?
b 3.006668013951515?