OpenGL绘制3D梯度下降动画（小批量梯度下降算法）

OpenGL绘制3D梯度下降动画（小批量梯度下降算法）

import copy
import random
import time
import math
from OpenGL.GL import *
from OpenGL.GLU import *
from OpenGL.GLUT import *
from OpenGL.arrays import vbo

import AutodiffEngine.autodiff as ad
import numpy as np
import pandas as pd
import threading
import three_d.freeTypeFont as glFreeType
from AutodiffEngine.gradient import gradient_color

IS_PERSPECTIVE = True  # 透视投影
VIEW = np.array([-0.8, 0.8, -0.8, 0.8, 0.5, 50.0])  # 视景体的left/right/bottom/top/near/far六个面
SCALE_K = np.array([1.0, 1.0, 1.0])  # 模型缩放比例
EYE = np.array([0.0, 1.0, 1.0])  # 眼睛的位置（默认z轴的正方向）
LOOK_AT = np.array([0.0, 0.0, 0.0])  # 瞄准方向的参考点（默认在坐标原点）
EYE_UP = np.array([0.0, 1.0, 0.0])  # 定义对观察者而言的上方（默认y轴的正方向）
WIN_W, WIN_H = 640, 480  # 保存窗口宽度和高度的变量
LEFT_IS_DOWNED = False  # 鼠标左键被按下
MOUSE_X, MOUSE_Y = 0, 0  # 考察鼠标位移量时保存的起始位置

input_colors = ["#b93f43", "#e4b1b6", "#037ee4"]
colors = gradient_color(input_colors, color_sum=300)
colors = np.array(colors)[:255, :] / 255
R = colors[:, 0]
G = colors[:, 1]
B = colors[:, 2]

# R = np.hstack((np.linspace(200, 255, 127), np.linspace(255, 47, 128)))[::-1] / 255
# G = np.hstack((np.linspace(3, 255, 127), np.linspace(255, 47, 128)))[::-1] / 255
# B = np.hstack((np.linspace(33, 200, 127), np.hstack((np.linspace(220, 255, 25), np.linspace(255, 150, 103)))))[
#     ::-1] / 255
A = np.ones(255, dtype=np.float) * 0.5
TAG_XYZ = [0.0, -0.2, 0.0]
loss_mean, loss_std = None, None
tracks = []  # 轨迹


def getposture():
    global EYE, LOOK_AT

    dist = np.sqrt(np.power((EYE - LOOK_AT), 2).sum())
    if dist > 0:
        phi = np.arcsin((EYE[1] - LOOK_AT[1]) / dist)
        theta = np.arcsin((EYE[0] - LOOK_AT[0]) / (dist * np.cos(phi)))
    else:
        phi = 0.0
        theta = 0.0

    return dist, phi, theta


DIST, PHI, THETA = getposture()  # 眼睛与观察目标之间的距离、仰角、方位角


def init():
    global our_font
    our_font = glFreeType.font_data("msyh.ttc", 20)
    glClearColor(0.0, 0.0, 0.0, 1.0)  # 设置画布背景色。注意：这里必须是4个参数
    glEnable(GL_DEPTH_TEST)  # 开启深度测试，实现遮挡关系
    glDepthFunc(GL_LEQUAL)  # 设置深度测试函数（GL_LEQUAL只是选项之一）


def getColor(scale):
    global R, G, B, A
    scale = min(scale, 254)
    return R[scale], B[scale], B[scale], A[scale]


def start_draw():
    global w0, w1, loss, TAG_XYZ, vbo_vertices, vbo_indices, tracks, scope
    XYZ = copy.deepcopy(TAG_XYZ)
    tracks.append(XYZ)

    # 绘制方形轨迹
    size = 10.0
    glPointSize(size)
    glBegin(GL_POINTS)
    glColor4f(0.0, 1.0, 0.0, 1.0)
    # glVertex3f(XYZ[0], XYZ[1] + size / 5, XYZ[2])
    glVertex3f(XYZ[0], XYZ[1], XYZ[2])
    glEnd()

    # 绘制球形轨迹
    # glPushMatrix()
    # glTranslatef(XYZ[0], XYZ[1], XYZ[2])
    # glColor4f(0, 1, 1, 1)
    # glutWireSphere(scope * 0.08, 10, 10)
    # # glutSolidSphere(scope*0.08, 10, 10)
    # glPopMatrix()

    # 绘制曲面
    vbo_vertices.bind()
    # glInterleavedArrays(GL_V3F, 0, None)
    glInterleavedArrays(GL_C3F_V3F, 0, None)
    vbo_vertices.unbind()
    vbo_indices.bind()
    # glDrawElements(GL_QUADS, int(vbo_indices.size / 4), GL_UNSIGNED_INT, None)
    glDrawElements(GL_TRIANGLES, int(vbo_indices.size / 4), GL_UNSIGNED_INT, None)
    vbo_indices.unbind()

    # 绘制轨迹
    glEnable(GL_BLEND)
    # glEnable(GL_LINE_SMOOTH)
    glLineWidth(2.0)
    glBegin(GL_LINES)  # 开始绘制线段（世界坐标系）
    glColor4f(0.9, 0.9, 0.9, 1.0)  # 设置当前颜色为红色不透明
    for i in range(len(tracks) - 1):
        glVertex3f(tracks[i][0], float(tracks[i][1]), tracks[i][2])
        glVertex3f(tracks[i + 1][0], float(tracks[i + 1][1]), tracks[i + 1][2])
    glEnd()


def draw():
    global IS_PERSPECTIVE, VIEW
    global EYE, LOOK_AT, EYE_UP
    global SCALE_K
    global WIN_W, WIN_H
    global scope

    # 清除屏幕及深度缓存
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)

    # 设置投影（透视投影）
    glMatrixMode(GL_PROJECTION)
    glLoadIdentity()

    if WIN_W > WIN_H:
        if IS_PERSPECTIVE:
            glFrustum(VIEW[0] * WIN_W / WIN_H, VIEW[1] * WIN_W / WIN_H, VIEW[2], VIEW[3], VIEW[4],
                      VIEW[5])  # glFrustum() 用来设置平行投影
        else:
            glOrtho(VIEW[0] * WIN_W / WIN_H, VIEW[1] * WIN_W / WIN_H, VIEW[2], VIEW[3], VIEW[4],
                    VIEW[5])  # glOrtho() 用来设置平行投影
    else:
        if IS_PERSPECTIVE:
            glFrustum(VIEW[0], VIEW[1], VIEW[2] * WIN_H / WIN_W, VIEW[3] * WIN_H / WIN_W, VIEW[4], VIEW[5])
        else:
            glOrtho(VIEW[0], VIEW[1], VIEW[2] * WIN_H / WIN_W, VIEW[3] * WIN_H / WIN_W, VIEW[4], VIEW[5])

    # 设置模型视图
    glMatrixMode(GL_MODELVIEW)
    glLoadIdentity()

    # 几何变换
    glScale(SCALE_K[0], SCALE_K[1], SCALE_K[2])

    # 设置视点
    gluLookAt(
        EYE[0], EYE[1], EYE[2],
        LOOK_AT[0], LOOK_AT[1], LOOK_AT[2],
        EYE_UP[0], EYE_UP[1], EYE_UP[2]
    )

    # 设置视口
    glViewport(0, 0, WIN_W, WIN_H)

    # ---------------------------------------------------------------
    glLineWidth(1.0)
    glBegin(GL_LINES)  # 开始绘制线段（世界坐标系）

    # 以红色绘制x轴
    glColor4f(1.0, 0.0, 0.0, 1.0)  # 设置当前颜色为红色不透明
    glVertex3f(-scope, 0.0, 0.0)  # 设置x轴顶点（x轴负方向）
    glVertex3f(scope, 0.0, 0.0)  # 设置x轴顶点（x轴正方向）
    glVertex3f(scope, 0.0, 0.0)  # x轴正方向箭头
    glVertex3f(scope - 0.1 * scope, -0.05 * scope, 0.0)  # x轴正方向箭头

    # 以绿色绘制y轴
    glColor4f(0.0, 1.0, 0.0, 1.0)  # 设置当前颜色为绿色不透明
    glVertex3f(0.0, -scope, 0.0)  # 设置y轴顶点（y轴负方向）
    glVertex3f(0.0, scope, 0.0)  # 设置y轴顶点（y轴正方向）
    glVertex3f(0.0, scope, 0.0)  # y轴正方向箭头
    glVertex3f(-0.05 * scope, scope - 0.1 * scope, 0.0)  # y轴正方向箭头

    # 以蓝色绘制z轴
    glColor4f(0.0, 0.0, 1.0, 1.0)  # 设置当前颜色为蓝色不透明
    glVertex3f(0.0, 0.0, -scope)  # 设置z轴顶点（z轴负方向）
    glVertex3f(0.0, 0.0, scope)  # 设置z轴顶点（z轴正方向）
    glVertex3f(0.0, 0.0, scope)  # z轴正方向箭头
    glVertex3f(0.0, -0.05 * scope, scope - 0.1 * scope)  # z轴正方向箭头

    glEnd()  # 结束绘制线段

    start_draw()
    # ---------------------------------------------------------------
    glutSwapBuffers()  # 切换缓冲区，以显示绘制内容


def reshape(width, height):
    global WIN_W, WIN_H

    WIN_W, WIN_H = width, height
    glutPostRedisplay()


def mouseclick(button, state, x, y):
    global SCALE_K
    global LEFT_IS_DOWNED
    global MOUSE_X, MOUSE_Y

    MOUSE_X, MOUSE_Y = x, y
    if button == GLUT_LEFT_BUTTON:
        LEFT_IS_DOWNED = state == GLUT_DOWN
    elif button == 3:
        SCALE_K *= 1.05
        glutPostRedisplay()
    elif button == 4:
        SCALE_K *= 0.95
        glutPostRedisplay()


def mousemotion(x, y):
    global LEFT_IS_DOWNED
    global EYE, EYE_UP
    global MOUSE_X, MOUSE_Y
    global DIST, PHI, THETA
    global WIN_W, WIN_H

    if LEFT_IS_DOWNED:
        dx = MOUSE_X - x
        dy = y - MOUSE_Y
        MOUSE_X, MOUSE_Y = x, y

        PHI += 2 * np.pi * dy / WIN_H
        PHI %= 2 * np.pi
        THETA += 2 * np.pi * dx / WIN_W
        THETA %= 2 * np.pi
        r = DIST * np.cos(PHI)

        EYE[1] = DIST * np.sin(PHI) * scope
        EYE[0] = r * np.sin(THETA) * scope
        EYE[2] = r * np.cos(THETA) * scope

        if 0.5 * np.pi < PHI < 1.5 * np.pi:
            EYE_UP[1] = -1.0
        else:
            EYE_UP[1] = 1.0
        glutPostRedisplay()


def keydown(key, x, y):
    global DIST, PHI, THETA
    global EYE, LOOK_AT, EYE_UP
    global IS_PERSPECTIVE, VIEW

    if key in [b'x', b'X', b'y', b'Y', b'z', b'Z']:
        if key == b'x':  # 瞄准参考点 x 减小
            LOOK_AT[0] -= 0.1
        elif key == b'X':  # 瞄准参考 x 增大
            LOOK_AT[0] += 0.1
        elif key == b'y':  # 瞄准参考点 y 减小
            LOOK_AT[1] -= 0.1
        elif key == b'Y':  # 瞄准参考点 y 增大
            LOOK_AT[1] += 0.1
        elif key == b'z':  # 瞄准参考点 z 减小
            LOOK_AT[2] -= 0.1
        elif key == b'Z':  # 瞄准参考点 z 增大
            LOOK_AT[2] += 0.1

        DIST, PHI, THETA = getposture()
        glutPostRedisplay()
    elif key == b'\r':  # 回车键，视点前进
        EYE = LOOK_AT + (EYE - LOOK_AT) * 0.9
        DIST, PHI, THETA = getposture()
        glutPostRedisplay()
    elif key == b'\x08':  # 退格键，视点后退
        EYE = LOOK_AT + (EYE - LOOK_AT) * 1.1
        DIST, PHI, THETA = getposture()
        glutPostRedisplay()
    elif key == b' ':  # 空格键，切换投影模式
        IS_PERSPECTIVE = not IS_PERSPECTIVE
        glutPostRedisplay()


def logistic_prob(_w):
    def wrapper(_x):
        return 1 / (1 + np.exp(-np.sum(_x * _w)))

    return wrapper


def test_accuracy(_w, _X, _Y):
    prob = logistic_prob(_w)
    correct = 0
    total = len(_Y)
    for i in range(len(_Y)):
        x = _X[i]
        y = _Y[i]
        p = prob(x)
        if p >= 0.5 and y == 1.0:
            correct += 1
        elif p < 0.5 and y == 0.0:
            correct += 1
    print("总数：%d, 预测正确：%d" % (total, correct))


def gen_2d_data(n):
    np.random.seed(3)
    x_data = np.random.random([n, 1])
    y_data = np.ones(n)
    for i in range(n):
        d = x_data[i]
        if d[0] < 0.5:
            y_data[i] = 0
    x_data_with_bias = np.ones([n, 2])
    x_data_with_bias[:, 1:] = x_data
    return x_data_with_bias, y_data


def gen_2d_data2(n):
    np.random.seed()
    x_data = np.random.random([n, 2]) * 4 - 2

    w = np.random.random_integers(-5, 5, [2, 1])/1.0
    b = np.random.random_integers(-100, 100, [1])/100
    print(f"W = {w} b = {b}")
    hidden = np.matmul(x_data, w)[:, 0] + b
    y_data = 1 / (1 + np.exp(-hidden))
    #
    # y_data = np.where(np.sum(x_data, axis=1) > 0, 0, 1)

    # y_data = 1 - np.round(np.mean(np.square(x_data), axis=1), decimals=0)
    return x_data, y_data

def gen_2d_data3(n):
    np.random.seed()
    x_data1 = np.random.rand(int(n/2), 2)*0.2+0.2
    x_data2 = np.random.rand(int(n/2), 2)*0.2+0.7
    x_data = np.concatenate((x_data1, x_data2))

    y_data1 = np.zeros_like(x_data1[:, 0])
    y_data2 = np.ones_like(x_data1[:, 0])
    y_data = np.concatenate((y_data1, y_data2))

    c = list(zip(x_data, y_data))
    np.random.shuffle(c)
    x_data, y_data = zip(*c)
    x_data, y_data = np.array(x_data), np.array(y_data)
    return x_data, y_data


def cost(theta0, theta1):
    global data_x, data_y
    J = 0
    m = len(data_x)
    hidden = np.dot(data_x, np.array([theta0, theta1]))
    sigmoid = 1 / (1 + np.exp(-hidden))
    J = np.sum(np.square(sigmoid - data_y))
    # J = np.sum(-data_y * ad.log(sigmoid+0.01) - (1 - data_y) * ad.log(1 - sigmoid+0.01))
    return J


def move(args):
    glutPostRedisplay()
    glutTimerFunc(33, move, 1)


def train():
    global X_val, Y_val, w_val, TAG_XYZ, N, loss, loss_mean, loss_std
    x = ad.Variable(name='x')
    w = ad.Variable(name='w')
    y = ad.Variable(name='y')

    # 注意，以下实现某些情况会有很大的数值误差，
    # 所以一般真实系统实现会提供高阶算子，从而减少数值误差

    hidden = ad.matmul(x, w)
    sigmoid = 1 / (1 + ad.exp(-hidden))
    # L = -y * ad.log(sigmoid+0.01) - (1 - y) * ad.log(1 - sigmoid+0.01)
    L = ad.reduce_sum(ad.square(sigmoid - y))

    w_grad, = ad.gradients(L, [w])
    executor = ad.Executor([L, sigmoid, w_grad])

    test_accuracy(w_val, X_val, Y_val)
    batch = 10
    learning_rate = 0.1
    max_iters = 5000
    for iteration in range(max_iters):
        for i in range(int(N / batch)):
            x_val = X_val[i * batch:i * batch + batch, :]
            y_val = Y_val[i * batch:i * batch + batch]
            L_val, predict, w_grad_val = executor.run(feed_dict={w: w_val, x: x_val, y: y_val})
            w_val -= learning_rate * w_grad_val
            TAG_XYZ[0] = float(w_val[0])
            TAG_XYZ[2] = float(w_val[1])
            l = cost(float(w_val[0]), float(w_val[1]))
            TAG_XYZ[1] = (l - loss_mean) / loss_std * 5
            print(f"iter = {iteration}, w = {w_val}, loss={l}")
            time.sleep(0.01)
    test_accuracy(w_val, X_val, Y_val)


def deal_vbo(w0, w1, loss):
    s = time.time()
    vertices = []
    indices = []
    loss_flatten = loss.flatten()
    loss_flatten.sort()
    edge_max = loss_flatten[-20]
    edge_min = loss_flatten[20]
    for i in range(len(w0)):
        for j in range(len(w1)):
            point = [float(w0[i]), float(loss[i][j]), float(w1[j])]
            rgba = getColor(int((float(loss[i][j]) - loss.min()) / (loss.max() - loss.min()) * 255))
            if float(loss[i][j]) < edge_min:
                rgba = (1.0, 0.0, 0.0, 0.5)
            if float(loss[i][j]) > edge_max:
                rgba = (0.0, 1.0, 0.0, 0.5)
            vertices += rgba[:-1]
            vertices += point
            if i < len(w0) - 1 and j < len(w1) - 1:
                # indices += [i * len(w0) + j, (i + 1) * len(w0) + j, (i + 1) * len(w0) + j + 1, i * len(w0) + j + 1]
                indices += [i * len(w0) + j,
                            (i + 1) * len(w0) + j,
                            (i + 1) * len(w0) + j + 1
                            ]
                indices += [i * len(w0) + j,
                            (i + 1) * len(w0) + j + 1,
                            i * len(w0) + j + 1
                            ]
    print("ddd", time.time() - s)
    return np.array(vertices, np.float32), np.array(indices, np.int)


if __name__ == "__main__":
    N = 100
    # X_val, Y_val = gen_2d_data(N)
    X_val, Y_val = gen_2d_data2(N)
    data_x = X_val
    data_y = Y_val
    # w_val = np.zeros(2)
    # w_val = np.array([-0.0, -1.0])
    w_val = np.random.random_integers(-10, 10, [2])/1.0  # 初始化权重
    print(w_val)

    scope = 15
    EYE *= scope
    SCALE_K *= 1 / scope
    num = 200
    w0 = np.linspace(-scope, scope, num)
    w1 = np.linspace(-scope, scope, num)
    loss = np.empty(shape=(num, num))
    # w0, w1 = np.meshgrid(w0, w1)
    for i in range(num):
        for j in range(num):
            loss[i, j] = cost(w0[i], w1[j])
    loss_mean = loss.mean()
    loss_std = loss.std()
    loss = (loss - loss_mean) / loss_std * 5

    vertices, indices = deal_vbo(w0, w1, loss)

    vbo_vertices = vbo.VBO(vertices)
    vbo_indices = vbo.VBO(indices, target=GL_ELEMENT_ARRAY_BUFFER)

    glutInit()
    glutInitDisplayMode(GLUT_DOUBLE | GLUT_ALPHA | GLUT_DEPTH | GLUT_RGBA)

    glutInitWindowSize(WIN_W, WIN_H)
    glutInitWindowPosition(300, 200)
    glutCreateWindow('梯度下降动画'.encode("gbk"))

    init()  # 初始化画布
    glutDisplayFunc(draw)  # 注册回调函数draw()
    glutTimerFunc(33, move, 1)
    glutReshapeFunc(reshape)  # 注册响应窗口改变的函数reshape()
    glutMouseFunc(mouseclick)  # 注册响应鼠标点击的函数mouseclick()
    glutMotionFunc(mousemotion)  # 注册响应鼠标拖拽的函数mousemotion()
    glutKeyboardFunc(keydown)  # 注册键盘输入的函数keydown()

    t = threading.Thread(target=train, args=())
    t.start()

    glutMainLoop()  # 进入glut主循环

源码: Github.