Root-finding animation script

zetanim.py

"""
Animates output of

   build/examples/real_roots 0 0 50 -verbose > outlog.txt
   build/examples/real_roots 0 47 50 -verbose > outlog2.txt

assuming that the following hack of a patch has been applied to arb:

diff --git a/arb_calc/isolate_roots.c b/arb_calc/isolate_roots.c
index ac425fb..4981601 100644
--- a/arb_calc/isolate_roots.c
+++ b/arb_calc/isolate_roots.c
@@ -54,8 +54,21 @@ check_block(arb_calc_func_t func, void * param, const arf_interval_t block,
     arb_init(x);
 
     arf_interval_get_arb(x, block, prec);
+
+    if (arb_calc_verbose)
+    {
+        printf("block ");
+        arf_interval_printd(block, 15);
+    }
+
     func(t, x, param, 1, prec);
 
+    if (arb_calc_verbose)
+    {
+        printf("   value   ");
+        arb_printd(t, 30); printf("    ");
+    }
+
     result = BLOCK_UNKNOWN;
 
     if (arb_is_positive(t) || arb_is_negative(t))
@@ -75,6 +88,11 @@ check_block(arb_calc_func_t func, void * param, const arf_interval_t block,
         }
     }
 
+    if (arb_calc_verbose)
+    {
+        printf("status %d\n", result);
+    }
+
     arb_clear(t + 0);
     arb_clear(t + 1);
     arb_clear(x);
@@ -115,6 +133,7 @@ isolate_roots_recursive(arf_interval_ptr * blocks, int ** flags,
     else
     {
         *eval_count -= 1;
+
         status = check_block(func, param, block, asign, bsign, prec);
 
         if (status != BLOCK_NO_ZERO)

"""
import numpy
import scipy
import scipy.misc
import os

xmin = 0.0
xmax = 50.0
ymin = -8.0
ymax = 8.0

nx = 1280
ny = 720

#os.system("rm *.png")

frame = 1

os.system("convert -background white -fill black -size 1280x720 -pointsize 72 "
"-gravity center label:'"
"Using interval arithmetic to isolate\\n\\n"
"roots of the Z-function on [0,50]."
"' card0.png")

os.system("convert -background white -fill black -size 1280x720 -pointsize 64 "
"-gravity center label:'"
"Roots of Z(t) precisely match nontrivial\\n\\n"
"zeros of the Riemann zeta function:\\n\\n\\n"
"Z(t) = C(t) zeta(1/2+it) where C(t) != 0"
"' card0a.png")

os.system(u"convert -background white -fill black -size 1280x720 -pointsize 44 "
"-gravity center label:'"
"Root criteria for a generic function f(t) on \\[a,b\\]:\\n\\n"
"Exactly 0 roots (RED/BLUE): the enclosure of f(\\[a,b\\])\\n"
"is strictly positive (blue) or negative (red).\\n\\n"
"Exactly 1 root (GREEN): sgn(f(a)) sgn(f(b)) = -1 and the\\n"
"enclosure of the derivative of f on \\[a,b\\] does not contain 0.\\n\\n"
"Unknown (MAGENTA): the enclosure of f(\\[a,b\\]) contains 0,\\n"
"but the GREEN test failed. The algorithm splits\\n"
"\\[a,b\\] in half and recurses."
"' card0b.png")

os.system("convert -background white -fill black -size 1280x720 -pointsize 24 "
"-gravity center -pointsize 72 label:'Depth-first refinement\\n\\n(the actual order used by\\nthe root isolation algorithm).' "
"card1.png")

os.system("convert -background white -fill black -size 1280x720 -pointsize 24 "
"-gravity center -pointsize 72 label:'The same subdivisions, but\\nre-animated in breadth-first order.' "
"card2.png")

os.system("convert -background white -fill black -size 1280x720 -pointsize 24 "
"-gravity center -pointsize 72 label:'The same subdivisions, but\\nre-animated in breadth-first order.' "
"card2.png")

os.system("convert -background white -fill black -size 1280x720 -pointsize 24 "
"-gravity center -pointsize 72 label:'"
"Looking only at the interval \\[47,50\\],\\n"
"which contains two roots\\n"
"at 48.005... and 49.773...\\n\\n"
"Depth-first."
"' "
"card3.png")

os.system("convert -background white -fill black -size 1280x720 -pointsize 24 "
"-gravity center -pointsize 72 label:'The same subdivisions, but\\nre-animated in breadth-first order.' "
"card4.png")

os.system("convert card0.png -define png:color-type=2 card0.png")
os.system("convert card0a.png -define png:color-type=2 card0a.png")
os.system("convert card0b.png -define png:color-type=2 card0b.png")
os.system("convert card1.png -define png:color-type=2 card1.png")
os.system("convert card2.png -define png:color-type=2 card2.png")
os.system("convert card3.png -define png:color-type=2 card3.png")
os.system("convert card4.png -define png:color-type=2 card4.png")

for i in range(160):
    os.system("cp card0.png frame%04i.png" % frame)
    frame += 1

for i in range(240):
    os.system("cp card0a.png frame%04i.png" % frame)
    frame += 1

for i in range(720):
    os.system("cp card0b.png frame%04i.png" % frame)
    frame += 1

for i in range(160):
    os.system("cp card1.png frame%04i.png" % frame)
    frame += 1

datas = []
for line in open("outlog.txt", "r").readlines():
    if line.startswith("block"):
        line = line.split()
        ba = float(line[1].lstrip("[").rstrip(","))
        bb = float(line[2].rstrip("]"))
        vm = -float(line[4])   # y axis is flipped
        vr = float(line[6])
        st = int(line[8])
        datas.append((ba, bb, vm, vr, st))

def showthings():
    buf = scipy.zeros((nx, ny, 3), dtype=numpy.uint8)
    flagged = scipy.zeros(nx, dtype=numpy.bool)
    global frame, datas
    for (ba, bb, vm, vr, st) in datas:
        print "frame", frame
        xi1 = int((ba - xmin) * nx / (xmax - xmin))
        xi2 = int((bb - xmin) * nx / (xmax - xmin))
        try:
            yi1 = int((vm - vr - ymin) * ny / (ymax - ymin))
            yi2 = int((vm + vr - ymin) * ny / (ymax - ymin))
        except:
            yi1 = 0
            yi2 = ny
        if st == 1:
            xi1, xi2 = min(xi1-1, xi2), max(xi2, xi1+1)
            flagged[xi1:xi2] = 1
        xi1 = max(xi1, 0)
        xi2 = min(xi2, nx)
        yi1 = max(yi1, 0)
        yi2 = min(yi2, ny)
        # print ba, bb, vm, vr, st, "    ", xi1, xi2
        for i in xrange(xi1,xi2):
            if flagged[i]:
                buf[i,:,:] = (0, 224, 0)
                buf[i,yi1:yi2,:] = (0, 255, 0)
            else:
                buf[i,:yi1,:] = 255
                buf[i,yi2:,:] = 255

                if st == 2:
                    buf[i,yi1:yi2,:] = (255, 0, 255)
                elif st == 0:
                    if -vm > 0.0:
                        buf[i,yi1:yi2,:] = (0, 0, 255)
                    else:
                        buf[i,yi1:yi2,:] = (255, 0, 0)

        buf[:,ny/2,:] = 0
        scipy.misc.imsave("frame%04i.png" % frame, numpy.swapaxes(buf, 0, 1))
        frame += 1
    for i in range(120):
        os.system("cp frame%04i.png frame%04i.png" % (frame-1, frame))
        frame += 1

showthings()

for i in range(160):
    os.system("cp card2.png frame%04i.png" % frame)
    frame += 1

datas.sort(key=lambda t: t[1]-t[0], reverse=True)

showthings()

for i in range(240):
    os.system("cp card3.png frame%04i.png" % frame)
    frame += 1

xmin = 47.0
xmax = 50.0
ymin = -2.0
ymax = 2.0
datas = []
for line in open("outlog2.txt", "r").readlines():
    if line.startswith("block"):
        line = line.split()
        ba = float(line[1].lstrip("[").rstrip(","))
        bb = float(line[2].rstrip("]"))
        vm = -float(line[4])   # y axis is flipped
        vr = float(line[6])
        st = int(line[8])
        datas.append((ba, bb, vm, vr, st))

showthings()

for i in range(160):
    os.system("cp card4.png frame%04i.png" % frame)
    frame += 1

datas.sort(key=lambda t: t[1]-t[0], reverse=True)

showthings()

os.system("avconv -r 48 -i frame%04d.png -r 48 -c:v h264 -crf 1 test2.flv")