讨论/综合讨论/新人报到/
新人报到

能不能看一下还有什么可以优化的

如题。

代码片段


import os
import os.path as osp
from numpy import abs, sqrt, log
from numpy.random import randint, uniform, choice
from tqdm import trange
from PIL import Image, ImageDraw

if not osp.exists("./output"):
    os.mkdir("output") 

K = 1024
MAXITERATIONS = K * K
MAXCUTOUT = MAXITERATIONS
NEXAMPLES = 1000000


def createattractor():
    global xe, ye, dx, dy
    for n in range(NEXAMPLES):
        lyapunov = 0
        a, x, y = [], [], []
        
        # Initialize coefficients for this attractor
        r = 4
        k = choice(12, r, replace=False)
        for i in range(12):
            a.append(0)
        for i in k:    
            while a[i] == 0:
                a[i] = float(randint(-20, 21) / 10)

        # Calculate the attractor
        drawit = True
        x0 = 0
        y0 = 0
        x.append(x0)
        y.append(y0)
        
        d0 = -1
        while d0 <= 0:
            dx = float(uniform(-1, 1, 1) / 1e5)
            dy = float(uniform(-1, 1, 1) / 1e5)
            d0 = sqrt(dx * dx + dy * dy)
        xe = x[0] + dx
        ye = y[0] + dy

        for i in trange(MAXITERATIONS + MAXCUTOUT, ascii=True):
            # Calculate next term
            def f1(x, y):
                temp = a[0] + a[1]*x + a[2]*y + a[3]*x*x +\
                       a[4]*x*y + a[5]*y*y
                return float(temp)

            def f2(x, y):
                temp = a[6] + a[7]*x + a[8]*y + a[9]*x*x +\
                       a[10]*x*y + a[11]*y*y
                return float(temp)
            
            x.append(f1(x[i-1], y[i-1]))
            y.append(f2(x[i-1], y[i-1]))
            xenew = f1(xe, ye)
            yenew = f2(xe, ye)

            # Check Avalid
            b = abs(a)
            q1 = b[3]+b[4]+b[5]+b[9]+b[10]+b[11]
            q2 = b[2]+b[4]+b[5]
            q3 = b[7]+b[9]+b[10]
            q4 = b[0]+b[6]
            K0 = (q1 * q2 * q3 * q4 == 0)
            if K0:
                drawit = False
                print("unwanted attractor")
                break

            # Does the series tend to infinity
            if x[i] < -1e10 or y[i] < -1e10 or\
               x[i] > 1e10 or y[i] > 1e10:
                drawit = False
                print("infinite attractor")
                break

            # Does the series tend to a point
            dx = x[i] - x[i-1]
            dy = y[i] - y[i-1]
            if abs(dx) < 1e-10 and abs(dy) < 1e-10:
                drawit = False
                print("point attractor")
                break

            # Calculate the lyapunov exponents
            if i >= MAXCUTOUT:
                dx = x[i] - xenew
                dy = y[i] - yenew
                dd = sqrt(dx * dx + dy * dy)
                lyapunov += log(abs(dd / d0))
                xe = x[i] + d0 * dx / dd
                ye = y[i] + d0 * dy / dd
            
        # Classify the series according to lyapunov
        if drawit:
            if abs(lyapunov) < K * K:
                print("neutrally stable")
                drawit = False
            elif lyapunov < 0:
                print("periodic {}".format(lyapunov))
                drawit = False 
            else:
                print("chaotic {}".format(lyapunov))
            
        # Save the image
        if drawit:
            xmin = +1e10
            xmax = -1e10
            ymin = +1e10
            ymax = -1e10
            for i in trange(MAXCUTOUT, MAXITERATIONS + MAXCUTOUT,
                            ascii=True):
                # Update the bounds
                xmin = min(xmin, x[i])
                xmax = max(xmax, x[i])
                ymin = min(ymin, y[i])
                ymax = max(ymax, y[i])
            if xmax - xmin < 1e-10 or ymax - ymin < 1e-10:
                continue             
            saveattractor(n, r, a, xmin, xmax, ymin, ymax,
                          x, y, lyapunov)


def saveattractor(n, r, a, xmin, xmax, ymin, ymax, x, y,
                  lyapunov):
    D = K

    # Save the parameters
    if not osp.exists("./output/{}".format(r)):
        os.mkdir("output/{}".format(r))
    with open("./output/{}/{}.txt".format(r, n), "w") as f:
        f.write("{}<=x<={} {}<=y<={}\n".format(xmin, xmax,
                                                         ymin, ymax))
        f.write("x0={} y0={}\n".format(x[0], y[0]))
        for i in range(len(a)):
            f.write("a[{}]={}\n".format(i, a[i]))
        f.write("negelegant={}\n".format(r))
        f.write("lyapunov={}\n".format(lyapunov / MAXITERATIONS))
        f.close()
    
    # Save the image
    image = Image.new("RGB", (D, D), "white")

    draw = ImageDraw.Draw(image)
    for i in trange(MAXCUTOUT, MAXITERATIONS + MAXCUTOUT, ascii=True):
        ix = D * (x[i] - xmin + 1e-10) / (xmax - xmin + 2e-10)
        iy = D * (y[i] - ymin + 1e-10) / (ymax - ymin + 2e-10)
        draw.point([ix, iy], fill="black")
    
    image.save("./output/{}/{}.png".format(r, n), "PNG")
    print("saved attractor {}-{}\n".format(r, n))


createattractor()
os.system("pause")

展开讨论
Grotex发起于 2020-05-18
共 0 个讨论
无讨论