이번엔 소나 데이터로 저주파 통과 필터 만듬

내용은 단순하다

x_k = alpha * x_k-1 + (1 - alpha) * x_k

평균 이동 필터가 새로 들어온 값을 1/n 만큼 반영하다보니

필터 결과가 뒤쳐짐

 

저주파 통고가 필터에서는 alpha를 설정해줘서 1/n만큼이 아닌

위의 경우 (1 - alpha) 만큼 반영시켜서 좀 반응이 빠름

 

 

이전에 구현한 소나 데이터 로드함수 그대로 쓰고

import scipy.io

sonar = scipy.io.loadmat("SonarAlt.mat")

def get_sonar():
    sonar_arr = sonar["sonarAlt"][0]
    sonar_len = len(sonar_arr)
    for i in range(sonar_len):
        yield sonar_arr[i]

for val in get_sonar():
    print(val)

 

매트랩 코드 보니 조금 다른데?

그냥 이동 평균 구현한거 고쳐서 만듬

avg 계산은 없고 그대로 쓴다.

이건 아닌거가아서 이동 평균 구현한걸 lpf로 만들어버림

 

 

 

구현은 했는데 위에 이동 평균이랑 별 차이를 모르겠다 이동평균이랑 lpf랑 같이 섞어봄

import matplotlib.pyplot as plt

def lpf_filter(buf_list, alpha):
    prev_avg = sum(buf_list[:-1]) / len(buf_list[:-1])
    avg = prev_avg * alpha +  (1 - alpha) * buf_list[-1]
    return avg


idx_list, buf_list, val_list, avg_list = [], [], [], []
buf_len = 10
alpha = 0.7


for idx, val in enumerate(get_sonar()):
    if idx > 500:
        break
    buf_list.append(val)
    if len(buf_list) > buf_len:
        buf_list = buf_list[1:]
        avg = lpf_filter(buf_list, alpha)
        idx_list.append(idx)
        val_list.append(val)
        avg_list.append(avg)

plt.plot(idx_list, val_list, label="raw val")
plt.plot(idx_list, avg_list, label="filtered val")

 

 

alpha 값을 0.9로 했을때는 

이동 평균이랑 별 차이가 없었음 

alpha를 0.3으로 지정해서 새 값이 0.7만큼 가중을 시켰을때

실제 데이터에 좀 더 가깝고 반응이 빨랏음

 

500까지 하면 너무 새밀해서

200으로 샘플을 줄여 플로팅

import matplotlib.pyplot as plt

def lpf_filter(buf_list, alpha):
    prev_avg = sum(buf_list[:-1]) / len(buf_list[:-1])
    avg = prev_avg * alpha +  (1 - alpha) * buf_list[-1]
    return avg

def mov_avg_filter(buf_list):
    avg = sum(buf_list) / len(buf_list)
    return avg


idx_list, buf_list, val_list, lpf_list, mov_avg_list = [], [], [], [], []
buf_len = 10
alpha = 0.3


for idx, val in enumerate(get_sonar()):
    if idx > 200:
        break
    buf_list.append(val)
    if len(buf_list) > buf_len:
        buf_list = buf_list[1:]
        lpf = lpf_filter(buf_list, alpha)
        mov_avg = mov_avg_filter(buf_list)
        idx_list.append(idx)
        val_list.append(val)
        lpf_list.append(lpf)
        mov_avg_list.append(mov_avg)

plt.plot(idx_list, val_list, label="raw val")
plt.plot(idx_list, lpf_list, label="lpf val")
plt.plot(idx_list, mov_avg_list, label="mov avg val")
plt.legend()

 

진짜 쩌는 책이다

개발자를 위핸 수학 지대넓얕 이라 해도 될듯

 

 

많은 곳에서 사용되는 2, 3차원에서 이동, 회전이라던가

스칼라, 벡터장

심볼릭 방정식 푸는 내용

오일러 법을 이용한 물체 이동

최적의 궤적을 구하기 위한 경사 상승법

푸리에 변환

머신러닝까지

 

딴대 보면 보통 c++이나 matlab으로 구현되어있을법한 것들을

가장 가볍게 하기 좋은 언어인 python으로 구현되어있는 것도 모잘라

 

plt로 거의 대부분 다 시각화를 해줘서 이해하기 직관적이다.

수학적인 내용도 있지만 다른 수학책에 비하면 정말 난잡한 내용은 발라내서 깔끔하다

 

차근 차근따라가며 수학자의 사고방식이 이런걸까 싶은 느낌을 받을수 있다.

 

 

 

https://www.yes24.com/Product/Goods/106045879

 

프로그래머를 위한 수학 : 파이썬으로 하는 3D 그래픽스, 시뮬레이션, 머신러닝 - 예스24

개발자에게 필요 없는 수학은 없다롱런하고 싶은 프로그래머를 위한 핵심 비법!`3D 그래픽스, 시뮬레이션, 머신러닝`세 분야의 공통점은 수학적 논리와 사고방식이 필요하다는 것!프로그래밍에

www.yes24.com

 

 

 

 

 

이번에는 평균 이동 필터를 파이썬으로 구현하려고한다.

근데 여기서 이용하는 소나 데이터가 메트랩 확장자로 되어있다.

 

찾아보니까 사이파이쓰면 된다칸다

 

1. mat file 소나 데이터 읽기

sicpy io loadmat으로 mat 읽고

yield로 데이터 뽑아내도록 구현함

 

import scipy.io

sonar = scipy.io.loadmat("SonarAlt.mat")

def get_sonar():
    sonar_arr = sonar["sonarAlt"][0]
    sonar_len = len(sonar_arr)
    for i in range(sonar_len):
        yield sonar_arr[i]

for val in get_sonar():
    print(val)

 

2, 평균 필터 구현하기

대충 버퍼 사이즈 5라 치자

 

avg x_k = (x_k-4 + x_k-3 + x_k-2 + x_k-1 + x_k) / 5

 

avg x_k-1 같은 경우

= (x_k-4-1 + x_k-3-1 + x_k-2-1 + x_k-1-1 + x_k-1) / 5

= (x_k-5 + x_k-4 + x_k-3 + x_k-2 + x_k-1) / 5

이다 

 

avg x_k - avg x_k-1 로 우항 정리하면

= (x_k - x_k-5) / 5임

 

avg x_k = avg x_k-1 + (x_k - x_k-5) / 5

가 된다.

 

버퍼 사이즈를 5대신 n으로 놓자

avg x_k = avg x_k-1 + (x_k - x_k-n) / n다

 

다시 말을 정리하면

지금 x 평균 = 한번 전 x 평균 + (지금 x 값 + n번 전 x 평균) / n

 

근대 생각해보니까 파이썬에선 리스트 조절만하면 되지 이럴 필요가 없잔아--

위 내용은 무시해도 되고 

 

버퍼 사이즈 10 채운 상태서 평균 계산

다음 시간에선 오래된 값비우고, 새값 넣기 -> 평균계산

으로 루프 돌리면 된다.

 

 

대충 평균 이동필터 구현하긴 했는데

좀 지저분하게 만들어서 찜찜하다

idx 는 x축

x_list는 평균 필터 버퍼 리스트

xsaved는 소나 값 전체 모음

avg 는 평균 값 전체 모음 리스트

동작은 맞으니 pass

import matplotlib.pyplot as plt
def mov_avg_filter(x_list):
    avg = sum(x_list) / len(x_list)
    return avg
idx_list, x_list, xsaved_list, avg_list = [], [], [], []
buf_len = 10
for idx, val in enumerate(get_sonar()):
    if idx > 500:
        break
    idx_list.append(idx)
    x_list.append(val)
    xsaved_list.append(val)
    if len(x_list) < buf_len:
        avg_list.append(val)
    else:
        if len(x_list) > buf_len:
            x_list = x_list[1:]
        avg = mov_avg_filter(x_list)
        avg_list.append(avg)
plt.plot(idx_list, xsaved_list, label="sonar val")
plt.plot(idx_list, avg_list, label="filtered val")
plt.legend()

 

칼만필터는 어렵지않아

내가 이 책을 봤던게 19년도 쯤이었던거 같은데

그때 정말 설명을 너무 친절하게 해주셔서 감탄하면서 봤었다

당시엔 matlab있으니까 그냥했는데

지금은 메트랩이 없엉

 

또 2판이 새로나오면서 파티클필터도 추가되기도 했고

이참에 파이썬으로 작성해봐야되겟다 싶엇는데 구현하게 됨

 

 

1. 볼트 값 가져오기

평균 필터 파트의 get volt는

직류분 14.4에다가

정규분포 노이즈로 4 *  randn 사용

파이썬으로 구현하면 이럼

import numpy as np

def get_volt():
    w = 0 + 4 * np.random.normal()
    z = 14.4 + w
    return z

print(get_volt())

 

2. 평균 필터와 테스트 로직 구현하기

메트랩 코드를 파이썬으로 만들다보니

persistent 때매 좀 지저분 해지긴했지만 동작은 잘된다

 

import matplotlib.pyplot as plt
def avg_filter(prev_avg, k, x):

    if prev_avg == 0:
        return x, 1
    alpha = (k - 1) / k
    avg = alpha * prev_avg + (1 - alpha) * x
    k = k + 1
    return avg, k
t = np.arange(0, 10, 0.2)
sample_num = len(t)
avg_list, xm_list = [], []
prev_avg, k = 0, 0
for _ in range(sample_num):
    xm = get_volt()
    xm_list.append(xm)

    avg, k = avg_filter(prev_avg, k, xm)
    prev_avg = avg
    avg_list.append(avg)
plt.plot(t, xm_list, color="red", label="voltage")
plt.plot(t, avg_list, color="blue", label="filtered voltage")
plt.legend()

 

 

 

2차원 공간 발사체의 최적 발사 각도 구하기

(1) 발사체 궤적 플로팅

발사체의 최대거리 함수 r(theta)가 있을시, 거리를 최대로 하는 theta를 찾는 문제를 다룬다고 하자

일단 발사체 궤적과 플로팅 함수는 다음과 같다

 

vx, vz는 속도를 xy축으로 분해한 성분으로

궤적 루프는 z가 0이 될때 즉 발사체가 땅으로 떨어지면 종료한다.

import matplotlib.pyplot as plt
from matplotlib import cm
import numpy as np
from math import sin, cos, pi

def trajectory(theta, speed=20, height=0, dt=0.01, g=-9.81):
    vx = speed * cos(pi * theta / 180)
    vz = speed * sin(pi * theta / 180)
    t, x, z = 0, 0, height
    ts, xs, zs = [t], [x], [z]
    while z >= 0:
        t += dt
        vz += g * dt
        x += vx * dt
        z += vz * dt
        ts.append(t)
        xs.append(x)
        zs.append(z)
    return ts, xs, zs

def plot_trajectories(*trajs, show_seconds=False):
    for traj in trajs:
        xs, zs = traj[1], traj[2]
        plt.plot(xs, zs)
        if show_seconds:
            second_indices = []
            second = 0
            for i, t in enumerate(traj[0]):
                if t >= second:
                    second_indices.append(i)
                    second +=1
            plt.scatter([xs[i] for i in second_indices], [zs[i] for i in second_indices])
    xl = plt.xlim()
    plt.plot(plt.xlim(), [0, 0], c='k')
    plt.xlim(*xl)

    width = 7
    coords_height = (plt.ylim()[1] - plt.ylim()[0])
    coords_width = (plt.xlim()[1] - plt.xlim()[0])
    plt.gcf().set_size_inches(width, width * coords_height / coords_width)

 

 

일단 발사각 45, 60인 경우 발사 궤적은 이런식인데

45도일때 가장 멀리 나간다.

plot_trajectories(
    trajectory(45),
    trajectory(60)
)

 

 

(2) 발사체 궤적 성질보기

최대 발사각 구하려면 지금 처럼 모든 각도에 대한 궤적을 플로팅해서 가장 멀리가는걸 찾거나

r(theta)를 최대로 하는 theta를 최적화 기법을 찾아내야한다.

 

일단 궤적 성질이라하면 도착거리, 최대 높이, 체공 시간 정도가 있는데

궤적 함수로부터 다음과같이 뽑아낼 수 있음

 

def landong_position(traj):
    return traj[1][-1]

def hang_time(traj):
    return traj[0][-1]

def max_height(traj):
    return max(traj[2])

traj_45 = trajectory(theta=45)

print(landong_position(traj_45))
print(hang_time(traj_45))
print(max_height(traj_45))

 

(3) 여러 각도에서 도착 거리 보기

위에서 도착거리, 체공시간, 최대 높이에 대한 함수를 구현했다.

r(theta)는 발사체의 도착 거리이므로

여러 각도 값을 주어 산점도로 플로팅해보자

 

아까 구현한 궤적 함수의 경우 시간에 따른 도착 거리지만

지금 플로팅한것은 발사각 세타에 대한 높이를 나타낸다.

대충 봐도 45도쯤에서 발사 거리가 최대로 나온다.

angles = range(0, 90, 5)
landong_positions = [landong_position(trajectory(theta)) for theta in angles]
plt.scatter(angles, landong_positions)

 

(4) 여러 발사 각도에서 채공시간 보기

 

1초간 발사체 위치는 아까 플로팅 궤적 함수에서 옵션이 있어 true를 주면 이런식으로 나온다.

 

 

아깐 theta별 거리에 대한 함수를 만들어서 플로팅했다면

이번에는 theta별 채공 시간에 대해 만들어서 플로팅해보자

 

45도 즘에서 최대 거리가 됫던것과는 다르게

90도쯤에서 가장 채공시간이 높게 나온다.

test_angles = range(0, 181, 5)
hang_times = [hang_time(trajectory(theta)) for theta in test_angles]
plt.scatter(test_angles, hang_times)

 

 (5) 지표와 각도가 주어질때 산점도 플로팅 함수

아까 발사체 지표가 도착거리, 채공시간, 최대 높이 3가지가 있었는데

방금 구현한것처럼 지표를 주면 바로 산점도 띄우는 함수를 만들자

 

이런식으로 채공 시간을 플로팅했고,

최대 높이나 거리도 플로팅 가능

def plot_trajectories(metric, thetas, **settings):
    plt.scatter(thetas, [metric(trajectory(theta, **settings)) for theta in thetas])

plot_trajectories(hang_time, range(0, 181, 5))

 

최대 거리의 경우 발사 시작점 높이가 10인 가정시 다음과 같음

plot_trajectories(landong_position, range(0,90, 5), height=10)

 

(6) 함수로 최적 거리구하기

위에서 한내용은 대충 오일러 방법으로 구하고, 그래프 플로팅해서 눈대중으로 찾아보았다.

이번에는 함수로 구해봄

일단 45도때 궤적을 플롯해보자

지금까지 계속 보던 포물선 그래프가 나온다.

발사 거리는 z(t) = 0인 시점이 되는 t로 정해짐

z 방향 속도 z'(t) = 초기 속도 + 중력 속도(가속도 아님)로 정해지고 이를 적분하면

z 방향 거리 z(t)가 아래와 같이 나온다.

z(t)가 t에 대한 이차식이 되었으므로 근의 공식으로 해를 구할수 있음

t = -2vz/g를 

발사거리 r = vx * dt = |v| cos(theta) dt에 대입하여

각도별 거리

r(theta) = - 2|v|**2/g * sin(theta) * cos(theta)로 정리가 가능

 

다시 정리하자

1. z축 거리가 0이 되는 시간 t를 구함, z(t) = 0이되는 t를 구함

2. x축 거리 함수에 대입하여 theta별 거리 r(t)를 구함

 

이를 플로팅하면 아까 궤적 함수를 플로팅했을때랑 동일한 형태 결과가 나온다.

최대 거리를 구하려면? r(theta) = 발사 거리 함수가 주어졌을때

r'(theta) = 0, 극점이 되는 theta를 구하면 된다.

 

trj = trajectory(45)
ts, zs = trj[0], trj[2]
plt.plot(ts, zs)


def z(t):
    return 20 * sin(45*pi/180) * t + (-9.81/2)*t**2

def plot_function(f, xmin, xmax, **kwargs):
    ts = np.linspace(xmin, xmax, 1000)
    plt.plot(ts, [f(t) for t in ts], **kwargs)

plt.plot(ts,zs)
plot_function(z,0,2.9)



def r(theta):
    return (-2 * 20 * 20/-9.81) * sin(theta * pi / 180) * cos(theta * pi / 180)
plot_function(r, 0, 90)

 

3치원 궤적 다루기

(1) 3차원 궤적 플로팅

위에서 한 경우는 발사각 theta가 주어졌을때 xz축에서의 거리와 높이를 다루었음

이번엔 y축 phi를 추가하여 xyz 3d 공간에서의 궤적을 다룸

 

r(t) = 거리 했던거랑은 달리

이런 식으로 표현하면 되나 싶음

r(t, theta, phi) = (x, y, z)

def trajectory3d(theta, phi, speed= 20, height=0, dt=0.01, g=-9.81):
    vx = speed * cos(pi * theta / 180) * cos(pi * phi / 180)
    vy = speed * cos(pi * theta / 180) * sin(pi * phi / 180)
    vz = speed * sin(pi * theta / 180)
    t, x, y, z = 0, 0, 0, height
    ts, xs, ys, zs = [t], [x], [y], [z]
    while z >= 0:
        t += dt
        vz += g * dt
        x += vx * dt
        y += vy * dt
        z += vz * dt
        ts.append(t)
        xs.append(x)
        ys.append(y)
        zs.append(z)
    return ts, xs, ys, zs

def plot_trajectory3d(traj):
    fig = plt.gcf()
    fig.set_size_inches(7, 7)
    ax = fig.gca(projection='3d')
    ax.plot(traj[1], traj[2], traj[3])

plot_trajectory3d(trajectory3d(45, 45))

 

 

(2) 지형을 추가한 3차원 궤적 플로팅

위 3차원 궤적 플로팅에서는 지형이 존재하지 않음

3d 궤적함수와 플롯팅 함수에 지형 요소를 추가

 

지형 함수로는 릿지 함수 사용

(x^2 - 5 y ^2)/2500

 

궤적 계산시 z >= 지형(x, y)인 경우에만 루프하여 궤적 추적

3차원 플로팅시 지형 서피스도 띄움

 

traj(20, 0)

traj(20, 270)인 두 경우에 대한 3차원 궤적

 

def flat_ground(x, y):
    return 0

def ridge(x, y):
    return (x ** 2 - 5 * y **2) / 2500

def trajectory3d(theta, phi, speed=20, height=0, dt=0.01, g=-9.81, elevation=flat_ground):
    vx = speed * cos(pi * theta / 180) * cos(pi * phi / 180)
    vy = speed * cos(pi * theta / 180) * sin(pi * phi / 180)
    vz = speed * sin(pi * theta / 180)
    t, x, y, z = 0, 0, 0, height
    ts, xs, ys, zs = [t], [x], [y], [z]
    while z >= elevation(x, y):
        t += dt
        vz += g * dt
        x += vx * dt
        y += vy * dt
        z += vz * dt
        ts.append(t)
        xs.append(x)
        ys.append(y)
        zs.append(z)
    return ts, xs, ys, zs

def plot_trajectories_3d(*trajs, elevation=flat_ground, bounds=None, zbounds=None, shadows=False):
    fig, ax = plt.gcf(), plt.gca(projection='3d')
    fig.set_size_inches(7, 7)
    if not bounds:
        xmin = min([x for traj in trajs for x in traj[1]])
        xmax = max([x for traj in trajs for x in traj[1]])
        ymin = min([x for traj in trajs for x in traj[2]])
        ymax = max([x for traj in trajs for x in traj[2]])

        padding_x = 0.1 * (xmax-xmin)
        padding_y = 0.1 * (ymax-ymin)
        xmin -= padding_x
        xmax += padding_x
        ymin -= padding_y
        ymax += padding_y
    else:
        xmin, xmax, ymin, ymax = bounds
    plt.plot([xmin, xmax], [0, 0], [0, 0], c='k')
    plt.plot([0, 0], [ymin, ymax], [0, 0], c='k')
    g = np.vectorize(elevation)
    ground_x = np.linspace(xmin, xmax, 20)
    ground_y = np.linspace(ymin, ymax, 20)
    ground_x, ground_y = np.meshgrid(ground_x, ground_y)
    ground_z = g(ground_x, ground_y)
    ax.plot_surface(ground_x, ground_y, ground_z, cmap=cm.coolwarm, 
                    alpha=0.5,linewidth=0, antialiased=True)
    
    for traj in trajs:
        ax.plot(traj[1], traj[2], traj[3])

    if zbounds:
        ax.set_zlim(*zbounds)

plot_trajectories_3d(
    trajectory3d(20, 0, elevation=ridge),
    trajectory3d(20, 270, elevation=ridge),
    bounds=[0,40,-40,0],
    elevation=ridge
)

 

(3) 3차원 공간에서 발사체 궤적 거리 구하기

2d 차원에서 theta를 x축에 놓고 거리를 y축에 놓아 그래프로 표현했었음

이번에는 3d에서 다룸

고도 함수 h(x, y) = (x^2 - 5y^2) / 2500이므로

h(x, y) = Bx^2 - Cy^2로 놓으면

B = 1/2500

C= 5/2500이 된다.

 

발사체의 고도는 시간 t에 따라 변하므로 h(t)라 일단 하고

x, y는 vx * t, vy * t로 구할수 있다.

h(x, y)에 이를 대입하면

 

h(vx * t, vy * t) = B vx^2  t^2 - C vy^2 t^2 가 된다.

vx, vy는 나두고 t에 대한 이차방정식으로 정리후 근의 공식으로 t에 대한 해를 구하면

 

t = -vz / (g/s - B vx^2 + cvy^2)가 온다.

다시 정리해보자

 

지금까지 릿지 경사에서 채공시간을 구하는 과정을 다루었다.

1. 고도 함수 h(x, y)가 있다

2. 발사체 x, y는 vx * t, vy * t로 구한다.

3. 발사체 xy 좌표를 고도함수에 대입한다.

4. x, y에 대한 고도함수를 t에 대한 이차식으로 정리한다.

5. h(vx * t, vy * t) = 0가 되는 t를 근의 공식으로 구한다. <- 채공시간

 

채공 시간을 구하였으면 xy 평면에서의 이동 거리를 구하면 된다.

vxy * 채공시간을 하면 되는데

vxy = 피타고라스정리 vx, vy로 구하면 됨

 

 

2차원에서는 theta 한 변수만 하면됬지만

지금은 theta와 phi 두 변수를 다루므로 따로 플로팅 x

from math import sqrt
B = 0.001
C = 0.005
v = 20
g = -9.81

def velocity_component(v, theta, phi):
    vx = v * cos(theta * pi / 180) * cos(phi * pi / 180)
    vy = v * cos(theta * pi / 180) * sin(phi * pi / 180)
    vz = v * sin(theta * pi / 180)
    return vx, vy, vz

def landing_distance(theta, phi):
    vx, vy, vz = velocity_component(v, theta, phi)
    v_xy = sqrt(vx**2 + vy**2)
    a = (g/2) - B * vx**2 + C * vy **2
    b = vz
    landing_time = -b/a
    landing_distance = v_xy * landing_time
    return landing_distance

 

 

(4) 항력 요소를 발사체 궤적에 가미하기

실제 발사체는 공기저항을 받으므로 위 식그대로 날라가지는 않음

vx, vy, vz요소에 항력 영향을 가중치 0.1로하여 가미하면 다음과 같음

 

def trajectory3d(theta, phi, speed=20, height=0, dt=0.01, g=-9.81, elevation=flat_ground,drag=0):
    vx = speed * cos(pi * theta / 180) * cos(pi * phi / 180)
    vy = speed * cos(pi * theta / 180) * sin(pi * phi / 180)
    vz = speed * sin(pi * theta / 180)
    t, x, y, z = 0, 0, 0, height
    ts, xs, ys, zs = [t], [x], [y], [z]
    while z >= elevation(x, y):
        t += dt
        vx -= (drag * vx) * dt
        vy -= (drag * vy) * dt
        vz += (g - (drag * vz)) * dt
        x += vx * dt
        y += vy * dt
        z += vz * dt
        ts.append(t)
        xs.append(x)
        ys.append(y)
        zs.append(z)
    return ts, xs, ys, zs

plot_trajectories_3d(
    trajectory3d(20, -20, elevation=ridge),
    trajectory3d(20, -20, elevation=ridge, drag=0.1),
    bounds=[0,40,-40,0],
    elevation=ridge
)

 

(5) xy 방향각에 따른 발사체 거리 스칼라 필드를 히트맵으로 시각화

아까 landing_distance 함수 구하긴 했는데 theta, phi에 따라 거리를 어떻게 시각화하지는 못했는데

거리 스칼라장을 히트맵으로 보여줄수 있음

히트맵으로 플로팅한 결과는 다음과 같다

 

발사각(y or x 축 회전) theta는 0 ~ 90도 사이

발사 방향(z축 회전) phi는 0 ~ 360도

 

릿지 지형에서 theta는 40, phi는 100, 270즘에서 가장 멀리 날라간다.

def scalar_field_heatmap(f, xmin, xmax, ymin, ymax, xsteps=100, ysteps=100):
    fig = plt.figure()
    fig.set_size_inches(7, 7)
    fv = np.vectorize(f)
    X = np.linspace(xmin, xmax, xsteps)
    Y = np.linspace(ymin, ymax, ysteps)
    X, Y = np.meshgrid(X, Y)
    z = fv(X, Y)
    
    fig, ax = plt.subplots()
    c = ax.pcolormesh(X, Y, z, cmap='plasma')
    ax.axis([X.min(), X.max(), Y.min(), Y.max()])
    fig.colorbar(c, ax=ax)

scalar_field_heatmap(landing_distance, 0, 90, 0, 360)

 

(6) 스칼라장을 히트맵 대신 3d에서 플로팅하기

스칼라장을 히트맵대신 서피스로 플로팅한 결과

봉우리 2개가 나오는데

이 두 지점에선 기울기가 0이 됨

 

def plot_scalar_field(f, xmin, xmax, ymin, ymax, xsteps=100, ysteps=100, c=None,cmap=cm.coolwarm, alpha=1,antialiased=False):
    fig = plt.gcf()
    fig.set_size_inches(7, 7)
    ax = fig.gca(projection='3d')
    fv = np.vectorize(f)
    X = np.linspace(xmin, xmax, xsteps)
    Y = np.linspace(ymin, ymax, ysteps)
    X, Y = np.meshgrid(X, Y)
    Z = fv(X, Y)
    surf = ax.plot_surface(X, Y, Z, cmap=cmap, color=c, alpha=alpha,
                          linewidth =0, antialiased=antialiased)

plot_scalar_field(landing_distance, 0, 90, 0, 360)

 

 

경사 상승을 이용한 최적화

(1) 거리 히트맵과 경사 벡터장

경사 상승이라니 딥러닝 공부할때 생각나네

딥러닝에서는 비용 함수가 최소가 되는 방향으로 학습을 시켰는데

여기서는 비행체 거리가 최대가 되는 방향을 찾아서 경사 상승이라 하나보다

 

일단 경사와 관련된 함수부터 구현하자

도함수 공식으로 시컨트 슬로프 함수 만들고

시컨트 슬로프로 근사 미분 함수를

근사 미분 함수로 근사 그라디언트 xy 방향 편미분 반환하는 함수 만든다..

그리고 이 함수와 도착거리함수로 도착 거리 그라디언트 함수 구현

def secant_slope(f, xmin, xmax):
    return (f(xmax) - f(xmin)) / (xmax - xmin)

def approx_derivative(f, x, dx=1e-6):
    return secant_slope(f, x-dx, x+dx)

def approx_gradient(f, x0, y0, dx=1e-6):
    partial_x = approx_derivative(lambda x: f(x, y0), x0, dx=dx)
    partial_y = approx_derivative(lambda y: f(x0, y), y0, dx=dx)
    return (partial_x, partial_y)

def landing_distance_gradient(theta, phi):
    return approx_gradient(landing_distance, theta, phi)

 

벡터장 플로팅 함수로 히트맵에다가 화살표들을 올림 

 

from vectors import to_polar, to_cartesian

def plot_vector_field(f, xmin, xmax, ymin, ymax, xsteps=10, ysteps=10, color='k'):
    X, Y = np.meshgrid(np.linspace(xmin,xmax,xsteps), np.linspace(ymin, ymax, ysteps))
    U = np.vectorize(lambda x, y : f(x, y)[0])(X, Y)
    V = np.vectorize(lambda x, y : f(x, y)[1])(X, Y)
    plt.quiver(X, Y, U, V, color=color)
    fig = plt.gcf()

scalar_field_heatmap(landing_distance, 0, 90, 0, 360)
plot_vector_field(landing_distance_gradient, 0, 90, 0, 360, xsteps=10, ysteps=10, color='k')
plt.xlabel('theta')
plt.ylabel('phi')
plt.gcf().set_size_inches(9, 7)

 

(2) 경사 상승법 구현하기

위에서 근사치 그라디언트 함수를 작성하여 벡터장 플로팅까지 했다

이번에는 이 근사 그라디언트로 경사 방향을 계속 올라 최대 지점까지 가는 함수 구현

 

36, 83을 줄시

끝점 37.58, 89.9에 도달함

 

from vectors import length

def gradient_ascent(f, xstart, ystart, tolerance=1e-6):
    x = xstart
    y = ystart
    grad = approx_gradient(f, x, y)
    while length(grad) > tolerance:
        x += grad[0]
        y += grad[1]
        grad = approx_gradient(f, x, y)
    return x, y

gradient_ascent(landing_distance, 36, 83)

 

 

(3) 경사 상승 궤적 그리기

그라디언트 어센트 함수로 시작점을 주면 최적의 지점까지 찾는건 잘 됬는데

궤적을 드로잉해보면 다음과 같음

 

def gradient_ascent_points(f, xstart, ystart, tolerance=1e-6):
    x = xstart
    y = ystart
    xs, ys = [x], [y]
    grad = approx_gradient(f, x, y)
    while length(grad) > tolerance:
        x += grad[0]
        y += grad[1]
        grad = approx_gradient(f, x, y)
        xs.append(x)
        ys.append(y)
    return xs, ys

scalar_field_heatmap(landing_distance, 35, 40, 80, 100)

plt.plot(*gradient_ascent_points(landing_distance, 36, 83), c='k')
plt.scatter([36,37.58114751557887],[83,89.99992616039857],c='k',s=75)
plt.xlabel('theta')
plt.ylabel('phi')
plt.gcf().set_size_inches(9, 7)

 

 

 

 

 

드디어 장을 다뤄보게 됬다

 

벡터장 플로팅 하기

좌표 x,y를 입력받아 반환하는 다음 함수를 주면 벡터를 반환

f(x, y) = (-2y, x)

 

import numpy as np
import matplotlib.pyplot as plt

def f(x, y):
    return (-2*y, x)

def plot_vector_field(f, xmin, xmax, ymin, ymax, xstep=1, ystep=1):
    X, Y = np.meshgrid(np.arange(xmin, xmax, xstep), np.arange(ymin, ymax, ystep))
    U = np.vectorize(lambda x, y : f(x, y)[0])(X, Y)
    V = np.vectorize(lambda x, y : f(x, y)[1])(X, Y)
    plt.quiver(X, Y, U, V, color='red')
    fig = plt.gcf()
    fig.set_size_inches(7,7)


plot_vector_field(f, -5, 5, -5, 5)

 

 

x, y = 1, 2 인경우

f(1, 2) = (-5, 1)인 벡터를 반환

meshgrid xy 둘다 -5 ~ 5 범위로 step은 1로 해서

이런 식으로 반환

 

 

-5, -5가 f에 들어갔다고 하자

f(-5, -5) = (10, -5)가 된다.

U, V 행렬 첫 값을 보면 생각한대로 나온다.

그리고 plt.quiver()는 2d 화살표 장을 플로팅 해주는 함수라고 한다.

 

 

위 사진에서 잘렸는데

x,y = -5, -5에서 벡터는 (10, -5)인데 사진 밖에 나가서 보이질 않는다

대신 x, y = -3, -3이라고 하면 벡터는 (6, -3)

-3, -3위치의 화살표를 보면 대충 느낌은 맞는거같다

 

 

이번엔 함수를 f(x,y) = (-x, -y)로 바꿔보면 이런식으로 중앙을 향하는 벡터장이 나온다.

 

 

 

벡터장의 중심을 -2, 4로 이동시키기

함수 -2-x, 4-y로 하면 벡터장의 중심이 -2, 4로 이동한 결과가 나온다.

 

 

스칼라장 플로팅하기

이번에는 스칼라장을 띄워봄

포텐셜 함수 u(x, y) = 0.5 * (x^2 + y^2)를 플로팅하면

 

 

from matplotlib.pyplot import cm
def plot_scalar_field(
            f, xmin, xmax, ymin, ymax, xstep=0.25, ystep=0.25, 
            c=None, cmap=cm.coolwarm, alpha=1, antialiased=False
        ):
    fig = plt.figure()
    fig.set_size_inches(7, 7)
    ax = fig.gca(projection='3d')
    fv = np.vectorize(f)
    X = np.arange(xmin, xmax, xstep)
    Y = np.arange(ymin, ymax, ystep)
    X, Y = np.meshgrid(X, Y)
    Z = fv(X, Y)
    surf = ax.plot_surface(X, Y, Z, cmap=cmap, color=c, alpha=alpha,
                             linewidth=0, antialiased=antialiased)
def u(x, y):
    return 0.5 * (x**2 + y**2)
plot_scalar_field(u, -5, 5, -5, 5)

 

 

히트맵으로 플로팅하면 이런식

def scalar_field_heatmap(f, xmin, xmax, ymin, ymax, xstep=0.1, ystep=0.1):
    fig = plt.figure()
    fig.set_size_inches(7,7)
    fv = np.vectorize(f)
    X = np.arange(xmin, xmax, xstep)
    Y = np.arange(ymin, ymax, ystep)
    X, Y = np.meshgrid(X, Y)
    z = fv(X, Y)
    fig, ax = plt.subplots()
    c = ax.pcolormesh(X, Y, z, cmap='plasma')
    ax.axis([X.min(), X.max(), Y.min(), Y.max()])
    fig.colorbar(c, ax=ax)
    plt.xlabel('x')
    plt.ylabel('y')
scalar_field_heatmap(u, -5, 5, -5, 5)

 

이번엔 스칼라장을 등고선으로 플로팅

def scalar_field_contour(f, xmin, xmax, ymin, ymax, levels=None):
    fv = np.vectorize(f)
    X = np.arange(xmin, xmax, 0.1)
    Y = np.arange(ymin, ymax, 0.1)
    X, Y = np.meshgrid(X, Y)
    Z = fv(X, Y)
    fig, ax = plt.subplots()
    CS = ax.contour(X, Y, Z, levels=levels)
    ax.clabel(CS, inline=1, fontsize=10, fmt='%1.1f')
    plt.xlabel('x')
    plt.ylabel('y')
    fig.set_size_inches(7, 7)

scalar_field_contour(u, -10, 10, -10, 10, levels=[10, 20, 30, 40, 50, 60])

 

 

지금까지 벡터장과 스칼라장을 플로팅 시켰고

이번에는 우주선 게임에 중력장을 추가

 

 

우주선 게임에 중력장 추가하기

진행할때마다 내용이 달라져서 처음부터 구현함

 

 

 

1. 폴리곤 모델 구현하기

이번 구현 게임은 화면 벽에 닿으면 반대로 튀어나갈수 있게 bounce 여부를 추가

맨앞에 bounce 있으며 bounce 켜짐 여부에 따라 반대편에 튀어나오게 할지, 벽에 부딪히면 튕겨지도록 move에 구현됨

 

그외 차이는 변수에 gravity와 mass가 추가되고

move 함수 내에 가속도를 나타내는 추진력 벡터와 중력장 파트가 추가됨

최종 가속도 = 추력 가속도 + 중력 가속도

 

import pygame
import vectors
from math import pi, sqrt, cos, sin, atan2
from random import randint, uniform
from linear_solver import do_segments_intersect
import numpy as np

boundce = True

class PolygonModel():
    def __init__(self, points):
        self.points = points
        self.rotation_angle = 0
        self.gravity = 0
        self.mass = 1
        self.x = 0
        self.y = 0
        self.vx = 0
        self.vy = 0
        self.angular_velocity = 0
        self.draw_center = False
    
    def transformed(self):
        rotated = [vectors.rotate2d(self.rotation_angle, v) for v in self.points]
        return [vectors.add((self.x, self.y), v) for v in rotated]

    def move(self, milliseconds, thrust_vector=(0, 0), gravity_sources=[]):
        ax, ay = thrust_vector
        gx, gy = gravitational_field(gravity_source, self.x, self.y)
        ax += gx
        ay += gy
        self.vx += ax * milliseconds / 1000
        self.vy += ay * milliseconds / 1000
        
        dx, dy = self.vx * milliseconds / 1000, self.vy * milliseconds / 1000
        self.x, self.y = vectors.add((self.x, self.y), (dx, dy))
        
        if boundce:
            if self.x < -10 or self.x > 10:
                self.vx = - self.vx
            if self.y < -10 or self.y > 10:
                self.vy = - self.vy
        else:
            if self.x < -10:
                self.x += 20
            if self.y < -10:
                self.y += 20
            if self.x > 10:
                self.x -= 20
            if self.y > 10:
                self.y -= 20
        self.rotation_angle += self.angular_velocity * milliseconds / 1000
    
    def segments(self):
        point_count = len(self.points)
        points = self.transformed()
        return [(points[i], points[(i+1)%point_count]) for i in range(0, point_count)]
    
    def does_collide(self, other_poly):
        for other_segment in other_poly.segments():
            if self.does_intersect(other_segment):
                return True
        return False

    def does_intersect(self, other_segment):
        for segment in self.segments():
            if do_segments_intersect(other_segment, segment):
                return True
        return False

 

2. 중력장 구현하기

폴리곤 모델로 블랙홀 클래스 구현

블랙홀은 폴리곤 모델 move 함수의 gravity sources로 들어감

블랙홀 1개만이 아닌 여러개 영향을 줄수 있음

 

gravitational_field에서는 (x, y) 벡터 반환하나

블랙홀들의 x축 중력, y축 중력들을 합하여 반환해줌

 

class BlackHole(PolygonModel):
    def __init__(self, gravity):
        vs = [vectors.to_cartesian((0.5, 2 * pi * i / 20)) for i in range(0, 20)]
        super().__init__(vs)
        self.gravity = gravity
    

def gravitational_field(sources, x, y):
    fields = [vectors.scale(-source.gravity, (x-source.x, y-source.y))
              for source in sources]
    return vectors.add(*fields)

black_hole = BlackHole(0.1)
black_hole.x, black_hole.y = 0, 0
black_holes = [
    black_hole
]

 

 

 

 

3. 메인함수 제외한 나머지 구현

 

먼저 우주선, 소행성 구현

변경사항은 우주선 시작점 빼곤 없음

class Ship(PolygonModel):
    def __init__(self):
        super().__init__([(0.5, 0,), (-0.25, 0.25), (-0.25, -0.25)])
    def laser_segment(self):
        dist = 20 * sqrt(2)
        x, y = self.transformed()[0]
        return (x, y), (x + dist * cos(self.rotation_angle), y + dist*sin(self.rotation_angle))

class Asteroid(PolygonModel):
    def __init__(self):
        sides = randint(5, 9)
        vs = [vectors.to_cartesian((uniform(0.5,1.0), 2 * pi * i/sides))
              for i in range(0, sides)]
        super().__init__(vs)
        self.vx = uniform(-1, 1)
        self.vy = uniform(-1, 1)
        self.angular_velocity = uniform(-pi/2, pi/2)


ship = Ship()
ship.x = 7
ship.y = 3
asteroid_count = 10
default_asteroids = [Asteroid() for _ in range(0, asteroid_count)]
for ast in default_asteroids:
    ast.x = randint(-9, 9)
    ast.y = randint(-9, 9)

 

나머지는 컬러, 화면크기, 드로잉 관련 함수

 

BLACK = (0, 0, 0)
WHITE = (255, 255, 255)
BLUE =  (0, 0, 255)
GREEN = (0, 255, 0)
RED =   (255, 0, 0)
LIGHT_GRAY =  (240, 240, 240)
DARK_GRAY = (128, 128, 128)
width, height = 400, 400

def to_pixels(x, y):
    return (width/2 + width * x / 20, height/2 - height * y / 20)

def draw_poly(screen, polygon_model, color=BLACK, fill=False):
    pixel_points = [to_pixels(x, y) for x, y in polygon_model.transformed()]
    if fill:
        pygame.draw.polygon(screen, color, pixel_points, 0)
    else:
        pygame.draw.lines(screen, color, True, pixel_points, 2)


def draw_segment(screen, v1, v2, color=RED):
    pygame.draw.line(screen, color, to_pixels(*v1), to_pixels(*v2), 2)

def draw_grid(screen):
    for x in range(-9, 10):
        draw_segment(screen, (x, -10), (x, 10), color=LIGHT_GRAY)
    for y in range(-9, 10):
        draw_segment(screen, (-10, y), (10, y), color=LIGHT_GRAY)
    draw_segment(screen, (-10, 0), (10, 0), color=DARK_GRAY)
    draw_segment(screen, (0, -10), (0, 10), color=DARK_GRAY)

thrust = 3

 

 

 

 

메인함수도 크게 바뀐건 없지만

우주선, 소행성 무브에 추력 벡터와 중력 소스가 추가되었으며

추력 벡터를 계산하는 과정

블랙홀을 드로잉하는 과정등이 추가되었음

 

전후진 입력없으면 추력 벡터는 (0, 0)이지만 

입력있을시 우주선 방향에 크기 3인 벡터를 x, y축 방향으로 분해해서 추력 벡터로 무브에 전달

 

 

바운스도 잘되고 블랙홀을 중심으로 움직인다.

def main(asteroids=default_asteroids):
    pygame.init()
    screen = pygame.display.set_mode([width, height])
    pygame.display.set_caption('asteroids')
    done = False
    clock = pygame.time.Clock()

    while not done:
        clock.tick()
        for event in pygame.event.get():
            if event.type == pygame.QUIT:
                done = True
        milliseconds = clock.get_time()
        keys = pygame.key.get_pressed()

        if keys[pygame.K_LEFT]:
            ship.rotation_angle += milliseconds * (2 * pi / 1000)
        if keys[pygame.K_RIGHT]:
            ship.rotation_angle -= milliseconds * (2 * pi / 1000)
        for ast in asteroids:
            ast.move(milliseconds, (0, 0), gravity_sources=black_holes)
        
        ship_thrust_vector = (0, 0)
        if keys[pygame.K_UP]:
            ship_thrust_vector = vectors.to_cartesian((thrust, ship.rotation_angle))
        elif keys[pygame.K_DOWN]:
            ship_thrust_vector = vectors.to_cartesian((-thrust, ship.rotation_angle))
        ship.move(milliseconds, ship_thrust_vector, black_holes)
        laser = ship.laser_segment()


        screen.fill(WHITE)
        draw_grid(screen)
        if keys[pygame.K_SPACE]:
            draw_segment(screen, *laser)
        draw_poly(screen, ship)
        for bh in black_holes:
            draw_poly(screen, bh, fill=True)
        for ast in asteroids:
            if keys[pygame.K_SPACE] and ast.does_intersect(laser):
                asteroids.remove(ast)
            else:
                draw_poly(screen, ast, color=GREEN)
        pygame.display.flip()
    pygame.quit()

if __name__ == "__main__":
    main()

 

 

 

 

이번엔 블랙홀을 2개로 늘린경우

소행성들이 두 블랙홀의 영향을 받아 움직인다.

 

black_hole = BlackHole(0.1)
black_hole.x, black_hole.y = 0, 0


black_hole2 = BlackHole(0.2)
black_hole.x, black_hole.y = -4, -4

black_holes = [
    black_hole,
    black_hole2
]

 

 

 

 

 

 

 

 

 

 

 

전체 코드

 

import pygame
import vectors
from math import pi, sqrt, cos, sin, atan2
from random import randint, uniform
from linear_solver import do_segments_intersect
import numpy as np

boundce = True

class PolygonModel():
    def __init__(self, points):
        self.points = points
        self.rotation_angle = 0
        self.gravity = 0
        self.mass = 1
        self.x = 0
        self.y = 0
        self.vx = 0
        self.vy = 0
        self.angular_velocity = 0
        self.draw_center = False
    
    def transformed(self):
        rotated = [vectors.rotate2d(self.rotation_angle, v) for v in self.points]
        return [vectors.add((self.x, self.y), v) for v in rotated]

    def move(self, milliseconds, thrust_vector=(0, 0), gravity_sources=[]):
        ax, ay = thrust_vector
        gx, gy = gravitational_field(gravity_sources, self.x, self.y)
        ax += gx
        ay += gy
        self.vx += ax * milliseconds / 1000
        self.vy += ay * milliseconds / 1000
        
        dx, dy = self.vx * milliseconds / 1000, self.vy * milliseconds / 1000
        self.x, self.y = vectors.add((self.x, self.y), (dx, dy))
        
        if boundce:
            if self.x < -10 or self.x > 10:
                self.vx = - self.vx
            if self.y < -10 or self.y > 10:
                self.vy = - self.vy
        else:
            if self.x < -10:
                self.x += 20
            if self.y < -10:
                self.y += 20
            if self.x > 10:
                self.x -= 20
            if self.y > 10:
                self.y -= 20
        self.rotation_angle += self.angular_velocity * milliseconds / 1000
    
    def segments(self):
        point_count = len(self.points)
        points = self.transformed()
        return [(points[i], points[(i+1)%point_count]) for i in range(0, point_count)]
    
    def does_collide(self, other_poly):
        for other_segment in other_poly.segments():
            if self.does_intersect(other_segment):
                return True
        return False

    def does_intersect(self, other_segment):
        for segment in self.segments():
            if do_segments_intersect(other_segment, segment):
                return True
        return False





class BlackHole(PolygonModel):
    def __init__(self, gravity):
        vs = [vectors.to_cartesian((0.5, 2 * pi * i / 20)) for i in range(0, 20)]
        super().__init__(vs)
        self.gravity = gravity
    

def gravitational_field(sources, x, y):
    fields = [vectors.scale(-source.gravity, (x-source.x, y-source.y))
              for source in sources]
    return vectors.add(*fields)

black_hole = BlackHole(0.1)
black_hole.x, black_hole.y = 0, 0


black_hole2 = BlackHole(0.2)
black_hole.x, black_hole.y = -4, -4

black_holes = [
    black_hole,
    black_hole2
]




class Ship(PolygonModel):
    def __init__(self):
        super().__init__([(0.5, 0,), (-0.25, 0.25), (-0.25, -0.25)])
    def laser_segment(self):
        dist = 20 * sqrt(2)
        x, y = self.transformed()[0]
        return (x, y), (x + dist * cos(self.rotation_angle), y + dist*sin(self.rotation_angle))

class Asteroid(PolygonModel):
    def __init__(self):
        sides = randint(5, 9)
        vs = [vectors.to_cartesian((uniform(0.5,1.0), 2 * pi * i/sides))
              for i in range(0, sides)]
        super().__init__(vs)
        self.vx = uniform(-1, 1)
        self.vy = uniform(-1, 1)
        self.angular_velocity = uniform(-pi/2, pi/2)


ship = Ship()
ship.x = 7
ship.y = 3
asteroid_count = 10
default_asteroids = [Asteroid() for _ in range(0, asteroid_count)]
for ast in default_asteroids:
    ast.x = randint(-9, 9)
    ast.y = randint(-9, 9)


BLACK = (0, 0, 0)
WHITE = (255, 255, 255)
BLUE =  (0, 0, 255)
GREEN = (0, 255, 0)
RED =   (255, 0, 0)
LIGHT_GRAY =  (240, 240, 240)
DARK_GRAY = (128, 128, 128)
width, height = 400, 400

def to_pixels(x, y):
    return (width/2 + width * x / 20, height/2 - height * y / 20)

def draw_poly(screen, polygon_model, color=BLACK, fill=False):
    pixel_points = [to_pixels(x, y) for x, y in polygon_model.transformed()]
    if fill:
        pygame.draw.polygon(screen, color, pixel_points, 0)
    else:
        pygame.draw.lines(screen, color, True, pixel_points, 2)

def draw_segment(screen, v1, v2, color=RED):
    pygame.draw.line(screen, color, to_pixels(*v1), to_pixels(*v2), 2)

def draw_grid(screen):
    for x in range(-9, 10):
        draw_segment(screen, (x, -10), (x, 10), color=LIGHT_GRAY)
    for y in range(-9, 10):
        draw_segment(screen, (-10, y), (10, y), color=LIGHT_GRAY)
    draw_segment(screen, (-10, 0), (10, 0), color=DARK_GRAY)
    draw_segment(screen, (0, -10), (0, 10), color=DARK_GRAY)

thrust = 3


def main(asteroids=default_asteroids):
    pygame.init()
    screen = pygame.display.set_mode([width, height])
    pygame.display.set_caption('asteroids')
    done = False
    clock = pygame.time.Clock()

    while not done:
        clock.tick()
        for event in pygame.event.get():
            if event.type == pygame.QUIT:
                done = True
        milliseconds = clock.get_time()
        keys = pygame.key.get_pressed()

        if keys[pygame.K_LEFT]:
            ship.rotation_angle += milliseconds * (2 * pi / 1000)
        if keys[pygame.K_RIGHT]:
            ship.rotation_angle -= milliseconds * (2 * pi / 1000)
        for ast in asteroids:
            ast.move(milliseconds, (0, 0), gravity_sources=black_holes)
        
        ship_thrust_vector = (0, 0)
        if keys[pygame.K_UP]:
            ship_thrust_vector = vectors.to_cartesian((thrust, ship.rotation_angle))
        elif keys[pygame.K_DOWN]:
            ship_thrust_vector = vectors.to_cartesian((-thrust, ship.rotation_angle))
        ship.move(milliseconds, ship_thrust_vector, black_holes)
        laser = ship.laser_segment()


        screen.fill(WHITE)
        draw_grid(screen)
        if keys[pygame.K_SPACE]:
            draw_segment(screen, *laser)
        draw_poly(screen, ship)
        for bh in black_holes:
            draw_poly(screen, bh, fill=True)
        for ast in asteroids:
            if keys[pygame.K_SPACE] and ast.does_intersect(laser):
                asteroids.remove(ast)
            else:
                draw_poly(screen, ast, color=GREEN)
        pygame.display.flip()
    pygame.quit()

if __name__ == "__main__":
    main()

 

 

 

대수식을 심볼릭 프로그래밍으로 풀기 위해

자료 구조로 만들기 진행

 

제곱, 수, 심볼 변수 정의하기

설명할 내용없이 단순함

class Power():
    def __init__(self, base, exponent):
        self.base = base
        self.exponent = exponent
    
class Number():
    def __init__(self, number):
        self.number = number

class Variable():
    def __init__(self, symbol):
        self.symbol = symbol

Power(Variable("x"), Number(2))

 

 

곱 클래스 만들기

실제 연산하는 내용은 없지만 단순하게 이것도 구현

이건 3 x^2

class Product():
    def __init__(self, exp1, exp2):
        self.exp1 = exp1
        self.exp2 = exp2

Product(Number(3), Power(Variable("x"), Number(2)))

 

 

약간 복잡한 함수 정의하기

Function은 함수

Apply는 함수와 해당 함수 인자

다음 함수식은

(3x^2 + x) * sin(x)

class Sum():
    def __init__(self, *exps):
        self.exps = exps
    
class Function():
    def __init__(self, name):
        self.name = name

class Apply():
    def __init__(self, function, argument):
        self.function = function
        self.argument = argument

f_expression = Product(
    Sum(
        Product(
            Number(3),
            Power(
                Variable("x"),
                Number(2)
            )
        ),
        Variable("x")
    ),
    Apply(
        Function("sin"),
        Variable('x')
    )
)

 

cos(x^3 - 5)는 다음과 같이 표현

 

Apply(
    Function("cos"),
    Sum(
        Power(
            Variable("x"),
            Number("3")
        ),
        Number(-5)
    )
)

 

 

 

로그 연산 표현해보기

ln(y^z)는 원래 쉽게 표현할수 있으나

지금 구현한 함수로는 다음과 같이 표현가능

from math import log
def f(y, z):
    return log(y ** z)

Apply(
    Function("ln"),
    Power(
        Variable("y"),
        Variable("z")
    )
)

 

 

 

 

 

 

나눗셈 구현하기

(a + b) / 2같은 경우는 다음과 같이 표현

class Quotient():
    def __init__(self, numerator, denominator):
        self.numerator = numerator
        self.denominator = denominator

Quotient(
    Sum(
        Variable("a"),
        Variable("b")
    ),
    Number(2)
)

 

 

차 연산 구현해서 이차 방정식 판별식 구현하기

이차 방정식 y = ax^2 + bx + c의 판별식 D = b  ^ 2 - 4 ac를 다음과 같이 표현 가능

 

class Difference():
    def __init__(self, exp1, exp2):
        self.exp1 = exp1
        self.exp2 = exp2

Difference(
    Power(
        Variable("b"),
        Number(2)
    ),
    Product(
        Number(4),
        Product(
            Variable("a"),
            Variable("c")
        )
    )
)

 

 

차연산 표현하기

차연산도 간단하게 이런식으로

class Negative():
    def __init__(self, exp):
        self.exp = exp

Negative(
    Sum(
        Power(
            Variable("x"),
            Number(2)
        ),
        Variable("y")
    )
)

 

 

 

 

2차 방정식 근의 공식 표현하기

근의 공식은 길어지지만 이런식으로

A = Variable('a')
B = Variable('b')
C = Variable('c')
Sqrt = Function('sqrt')

Quotient(
    Sum(
        Negative(B),
        Apply(
            Sqrt,
            Difference(
                Power(
                    B, Number(2)
                ),
                Product(
                    Number(4),
                    Product(A, C)
                )
            )
        )
    ),
    Product(
        Number(2),
        A
    )
)

 

식에 있는 모든 변수 찾아내기

def f(x):

   return x^3

기존의 파이썬 함수는 결과만 반환하지 안에 뭐가 들어있는지 모르므로

앞서 동작은 안하지만 대수 방정식 구조를 만듬

일단 기존에 만든 것에 어떤 변수가 있는지 반환하는 함수 만듬

 

일단 변수만 찾아내는 함수니까 넣은 변수는 잘 찾는거같다.

def distinct_variables(exp):
    if isinstance(exp, Variable):
        return (set(exp.symbol))
    elif isinstance(exp, Number):
        return set()
    elif isinstance(exp, Sum):
        return set().union(*[distinct_variables(exp) for exp in exp.exps])
    elif isinstance(exp, Product):
        return distinct_variables(exp.exp1).union(distinct_variables(exp.exp2))
    elif isinstance(exp, Power):
        return distinct_variables(exp.base).union(distinct_variables(exp.exponent))
    elif isinstance(exp, Apply):
        return distinct_variables(exp.argument)
    else:
        raise TypeError("Not a valid expression")

print(distinct_variables(Variable("z")))
print(distinct_variables(Number(3)))
print(distinct_variables(Sum(Number(3), Number(4))))
print(distinct_variables(Sum(Number(3), Variable("x"))))
print(distinct_variables(Sum(Variable("x"), Variable("y"))))
print(distinct_variables(Power(Variable("x"), Number(3))))
print(distinct_variables(Product(Power(Variable("x"), Number(3)), Variable("y"))))

 

 

방정식 동작하게 하기

아까 구현한 변수, 곱같은 클래스에 evaluate 함수 추가

evaluate로 변수를 바인딩해서 계산 수행

from abc import ABC, abstractmethod

class Expression(ABC):
    @abstractmethod
    def evaluate(self, **bindings):
        pass

class Number(Expression):
    def __init__(self, number):
        self.number = number
    def evaluate(self, **bindings):
        return self.number

class Variable(Expression):
    def __init__(self, symbol):
        self.symbol = symbol
    def evaluate(self, **bindings):
        try:
            return bindings[self.symbol]
        except:
            raise KeyError("Variable '{}' is not bound.".format(self.symbol))

class Product(Expression):
    def __init__(self, exp1, exp2):
        self.exp1 = exp1
        self.exp2 = exp2
    def evaluate(self, **bindings):
        return self.exp1.evaluate(**bindings) * self.exp2.evaluate(**bindings)

Product(Variable('x'), Variable('y')).evaluate(x=2, y=5)

 

 

타 연산 구현하기

방금 변수와 곱을 구현했고

나머지 연산들을 구현해보자

sin, cos, ln 같은 초월함수는 따로 Apply에서 바인딩 되도록 하고

합, 제곱, 차, 나눗셈 연산은 기본 연산대로 evaluate 추가

 

import math
from math import sin, cos, log

_function_bindings = {
    "sin":math.sin,
    "cos":math.cos,
    "ln":math.log
}

class Apply(Expression):
    def __init__(self, function, argument):
        self.function = function
        self.argument = argument
    def evaluate(self, **bindings):
        return _function_bindings[self.function.name](self.argument.evaluate(**bindings))


class Sum(Expression):
    def __init__(self, *exps):
        self.exps = exps
    def evaluate(self, **bindings):
        return sum([exp.evaluate(**bindings) for exp in self.exps])

class Power(Expression):
    def __init__(self, base, exponent):
        self.base = base
        self.exponent = exponent
    def evaluate(self, **bindings):
        return self.base.evaluate(**bindings) ** self.exponent.evaluate(**bindings)

class Difference(Expression):
    def __init__(self, exp1, exp2):
        self.exp1 = exp1
        self.exp2 = exp2
    def evaluate(self, **bindings):
        return self.exp1.evaluate(**bindings) - self.exp2.evaluate(**bindings)

class Quotient(Expression):
    def __init__(self, numerator, denominator):
        self.numerator = numerator
        self.denominator = denominator
    def evaluate(self, **bindings):
        return self.numerate.evaluate(**bindings) / self.denominator.evaluate(**bindings)

 

 

 

(3 x ^2 + x) * sin(x)

이 식을 심볼릭 방법으로 한것과

그냥 파이썬 함수에 대입한  결과를 보면 동일하게 잘 연산되었음

f_expression = Product(
                Sum(
                    Product(
                        Number(3),
                        Power(
                            Variable("x"),
                            Number(2)
                        )
                    ),
                    Variable("x")
                ),
                Apply(
                    Function("sin"),
                    Variable("x")
                )
            )

print(f_expression.evaluate(x=5))

def f(x):
    return (3 * x ** 2 + x) * sin(x)
print(f(5))

 

 

표현식 전개하기

지금까지 구한 표현식을 그대로 출력하면 클래스랑 주소를 알려주는데

내용을 전개해주면 좋으니 우선 추상 클래스 expression을 다음과 같이 구현

class Expression(ABC):
    @abstractmethod
    def evaluate(self, **bindings):
        pass
    @abstractmethod
    def expand(self):
        pass
    @abstractmethod
    def display(self):
        pass
    def __repr__(self):
        return self.display()

 

 

 

지금까지 구현한

숫자, 변수, 연산들을 다시 구현해보자

먼저 가장 간단한 숫자와 변수

아까와는 달리 주소가 아니라 클래스명과 값이 잘 나온다

class Number(Expression):
    def __init__(self, number):
        self. number = number
    def evaluate(self, **bindings):
        return self.number
    def expand(self):
        return self
    def display(self):
        return "Number({})".format(self.number)

class Variable(Expression):
    def __init__(self, symbol):
        self.symbol = symbol
    def evaluate(self, **bindings):
        return bindings[self.symbol]
    def expand(self):
        return self
    def display(self):
        return "Variable(\"{}\")".format(self.symbol)
print(Number(4))
print(Variable("x"))

 

합, 차연산 구현하고 동작시켜보면 이런식

 

 

class Sum(Expression):
    def __init__(self, *exps):
        self.exps = exps
    def evaluate(self, **bindings):
        return sum([exp.evaluate(**bindings) for exp in self.exps])
    def expand(self):
        return (Sum(*[exp.expand() for exp in self.exps]))
    def display(self):
        return "Sum({})".format(",".join([e.display() for e in self.exps]))

class Difference(Expression):
    def __init__(self, exp1, exp2):
        self.exp1 = exp1
        self.exp2 = exp2
    def evaluate(self, **bindings):
        return self.exp1.evaluate(**bindings) - self.exp2.evaluate(**bindings)
    def expand(self):
        return self
    def display(self):
        return "Difference({},{})".format(self.exp1.display(), self.exp2.display())

print(Sum(Number(3), Number(5)))
print(Sum(Variable("x"), Number(3)))
print(Difference(Number(3), Number(5)))
print(Difference(Variable("x"), Number(3)))
print(Difference(Variable("x"), Number(3)).evaluate(x=4))
print(Sum(Variable("x"), Difference(Number(3), Number(5))))

tmp_exp = Sum(Variable("x"), Number(3))
print(tmp_exp)
tmp_exp.expand()

 

 

 

 

나머지 연산도 구현

class Product(Expression):
    def __init__(self, exp1, exp2):
        self.exp1 = exp1
        self.exp2 = exp2
    def evaluate(self, **bindings):
        return self.exp1.evaluate(**bindings) * self.exp2.evaluate(**bindings)
    def expand(self):
        expanded1 = self.exp1.expand()
        expanded2 = self.exp2.expand()
        if isinstance(expanded1, Sum):
            return Sum(*[Product(e, expanded2).expand() for e in expanded1.exps])
        elif isinstance(expanded2, Sum):
            return Sum(*[Product(expanded1, e) for e in expanded2.exps])
        else:
            return Product(expanded1, expanded2)
    def display(self):
        return "Product({},{})".format(self.exp1.display(), self.exp2.display())

class Quotient(Expression):
    def __init__(self, numerator, denominator):
        self.numerator = numerator
        self.denominator = denominator
    def evaluate(self, **bindings):
        return self.numerator.evaluate(**bindings) / self.denominator.evaluate(**bindings)
    def expand(self):
        return self
    def display(self):
        return "Quotient({},{})".format(self.numerator.display(),self.denominator.display())

class Negative(Expression):
    def __init__(self, exp):
        self.exp = exp
    def evaluate(self, **bindings):
        return - self.exp.evaluate(**bindings)
    def expand(self):
        return self
    def display(self):
        return "Negative({})".format(self.exp.display())
    
class Power(Expression):
    def __init__(self, base, exponent):
        self.base = base
        self.exponent = exponent
    def evaluate(self, **bindings):
        return self.base.evaluate(**bindings) ** self.exponent.evaluate(**bindings)
    def expand(self):
        return self
    def display(self):
        return "Power({},{})".format(self.base.display(),self.exponent.display())

class Function():
    def __init__(self, name, make_latex=None):
        self.name = name
        self.make_latex = make_latex
    def latex(self, arg_latex):
        if self.make_latex:
            return self.make_latex(arg_latex)
        else:
            return " \\operatorname{{ {} }} \\left( {} \\right)".format(self.name, arg_latex)

class Apply(Expression):
    def __init__(self, function, argument):
        self.function = function
        self.argument = argument
    def evaluate(self, **bindings):
        return _function_bindings[self.function.name](self.argument.evaluate(**bindings))
    def expand(self):
        return Apply(self.function, self.argument.expand())
    def display(self):
        return "Apply(Function(\"{}\"),{})".format(self.function.name, self.argument.display())

 

식 간단하게 만들고 그냥 출력했을때랑

전개했을때 차이는 이런식

(a + b) * (y + z)

= ay + az + by + bz

출력 결과는 길어서 좀 그렇긴한데 결과도 잘 나오긴 함

Y = Variable('y')
Z = Variable('z')
A = Variable('a')
B = Variable('b')
Product(Sum(A, B), Sum(Y, Z))


Product(Sum(A, B), Sum(Y, Z)).expand()




f_expression = Product(
                Sum(
                    Product(
                        Number(3),
                        Power(
                            Variable("x"),
                            Number(2)
                        )
                    ),
                    Variable("x")
                ),
                Apply(
                    Function("sin"),
                    Variable("x")
                )
            )
f_expression.expand()

 

미분, 적분 내용까지 가긴 너무 힘드니 대충 sympy로 마무리

 

Sympy

기호로 수학 푸는 파이썬 lib

지금까지 구현한거랑 거의 비슷

숫자 그대로 써도되고, 심볼 쓸수있고

아래 내용은 방정식에 대입, 미분, 적분 결과

from sympy import *
from sympy.core.core import *

Mul(Symbol('y'), Add(3, Symbol('x')))



y = Symbol('y')
x = Symbol('x')
print(y * (3+x))
print(y * (3+x).subs(x, 1)) # substitude
print((x**2).diff(x)) # derivative
print((3*x**2).integrate(x)) # integrate

 

 

이번에는 오일러 방법으로 지난번에 만든 게임에 이동을 추가함

사실 지난번에 만든 게임에는 우주선 이동을 안넣엇지 운석 이동을 시키면서 이미 오일러 방법을 사용했었다.

 

 

변화 거리 = 속도 * 소요시간(밀리세컨드) / 1000

기존 위치 = 이전 위치 + 변화 거리

 

 

 

그래도 게임을 개선한건 맞으니 그냥 쭉 정리하면

오일러 방법으로 무브는 하는데

벽을 넘어가면 반대편에 나오도록 수정

 

class PolygonModel():
    def __init__(self, points):
        self.points = points
        self.rotation_angle = 0
        self.x = 0
        self.y = 0
        self.vx = 0
        self.vy = 0
        self.angular_velocity = 0
        self.draw_center = True
    
    def transformed(self):
        rotated = [vectors.rotate2d(self.rotation_angle, v) for v in self.points]
        return [vectors.add((self.x, self.y), v) for v in rotated]
    
    def move(self, milliseconds):
        dx, dy = self.vx * milliseconds / 1000.0, self.vy * milliseconds / 1000
        self.x, self.y = vectors.add((self.x, self.y), (dx, dy))
        if self.x < -10:
            self.x += 20
        if self.y < -10:
            self.y += 20
        if self.x > 10:
            self.x -= 20
        if self.y > 10:
            self.y -= 20
        self.rotation_angle += self.angular_velocity * milliseconds / 1000
    
    def segments(self):
        point_count = len(self.points)
        points = self.transformed()
        return [(points[i], points[(i+1)%point_count]) for i in range(0, point_count)]

    def does_collide(self, other_poly):
        for other_segment in other_poly.segments():
            if self.does_intersect(other_segment):
                return True
        return False

    def does_intersect(self, other_segment):
        for segment in self.segments():
            if do_segments_intersect(other_segment, segment):
                return True
        return False

 

 

우주선, 운석도 변경사항은 없음

 

class Ship(PolygonModel):
    def __init__(self):
        super().__init__([(0.5, 0), (-0.25, 0.25), (-0.25, -0.25)])
    
    def laser_segment(self):
        dist = 20 * sqrt(2)
        x, y = self.transformed()[0]
        return (x, y), (x+dist*cos(self.rotation_angle), y + dist * sin(self.rotation_angle))

class Asteroid(PolygonModel):
    def __init__(self):
        sides = randint(5, 9)
        vs = [vectors.to_cartesian((uniform(0.5, 1.0), 2 * pi * i / sides))
            for i in range(0, sides)]
        super().__init__(vs)
        self.vx = uniform(-1, 1)
        self.vy = uniform(-1, 1)
        self.angular_velocity = uniform(-pi/2, pi/2)

 

 

폴리 그리기 함수에서 중심점 파트랑

그리드 그리기 함수가 추가

 

BLACK = (0, 0, 0)
WHITE = (255, 255, 255)
BLUE =  (0, 0, 255)
GREEN = (0, 255, 0)
RED =   (255, 0, 0)
LIGHT_GRAY =  (240, 240, 240)
DARK_GRAY = (128, 128, 128)
width, height = 400, 400

def to_pixel(x, y):
    return (width / 2 + width * x / 20, height/2 - height * y / 20)

def draw_poly(screen, polygon_model, color=BLACK):
    pixel_points = [to_pixel(x, y) for x, y in polygon_model.transformed()]
    pygame.draw.lines(screen, color, True, pixel_points, 2)
    if polygon_model.draw_center:
        cx, cy = to_pixel(polygon_model.x, polygon_model.y)
        pygame.draw.circle(screen, BLACK, (int(cx), int(cy)), 4, 4)
    
def draw_segment(screen, v1, v2, color=RED):
    pygame.draw.line(screen, color, to_pixel(*v1), to_pixel(*v2), 2)

def draw_grid(screen):
    for x in range(-9, 10):
        draw_segment(screen, (x, -10), (x, 10), color=LIGHT_GRAY)
    for y in range(-9, 10):
        draw_segment(screen, (-10, y), (10, y), color=LIGHT_GRAY)
    
    draw_segment(screen, (-10, 0), (10, 0), color=DARK_GRAY)
    draw_segment(screen, (0, -10), (0, 10), color=DARK_GRAY)

 

 

 

 

메인 함수 진입 전에

우주선이랑, 운석, 가속도 변수 초기화

 

 

메인함수가 지난번이랑 다른건

우주선 좌우전후, 그리드 그리기가 추가된게 다다

 

ship = Ship()
asteroid_count = 10
default_asteroids = [Asteroid() for _ in range(0, asteroid_count)]
for ast in default_asteroids:
    ast.x = randint(-9, 9)
    ast.y = randint(-9, 9)


acceleration = 3
def main(asteroids=default_asteroids):
    pygame.init()
    screen = pygame.display.set_mode([width, height])
    done = False
    clock = pygame.time.Clock()

    while not done:
        clock.tick()

        for event in pygame.event.get():
            if event.type == pygame.QUIT:
                done = True
        milliseconds = clock.get_time()
        keys = pygame.key.get_pressed()

        for ast in asteroids:
            ast.move(milliseconds)
        
        if keys[pygame.K_LEFT]:
            ship.rotation_angle += milliseconds * (2 * pi / 1000)
        if keys[pygame.K_RIGHT]:
            ship.rotation_angle -= milliseconds * (2 * pi / 1000)
        if keys[pygame.K_UP]:
            ax = acceleration * cos(ship.rotation_angle)
            ay = acceleration * sin(ship.rotation_angle)
            ship.vx += ax * milliseconds/1000
            ship.vy += ay * milliseconds/1000
        if keys[pygame.K_DOWN]:
            ax = acceleration * cos(ship.rotation_angle)
            ay = acceleration * sin(ship.rotation_angle)
            ship.vx -= ax * milliseconds/1000
            ship.vy -= ay * milliseconds/1000
        ship.move(milliseconds)

        laser = ship.laser_segment()
        screen.fill(WHITE)
        draw_grid(screen)
        if keys[pygame.K_SPACE]:
            draw_segment(screen, *laser)
        draw_poly(screen, ship)
        for ast in asteroids:
            if keys[pygame.K_SPACE] and ast.does_intersect(laser):
                asteroids.remove(ast)
            else:
                draw_poly(screen, ast, color=GREEN)
        pygame.display.flip()
    pygame.quit()

 

 

 

 

동작 결과는 이럼

우주선도 움직이니까 이젠 게임같긴하다

 

 

 

 

 

import pygame
import vectors
from math import pi, sqrt, cos, sin, atan2
from random import randint, uniform
from linear_solver import do_segments_intersect
import numpy as np
import sys

class PolygonModel():
    def __init__(self, points):
        self.points = points
        self.rotation_angle = 0
        self.x = 0
        self.y = 0
        self.vx = 0
        self.vy = 0
        self.angular_velocity = 0
        self.draw_center = True
    
    def transformed(self):
        rotated = [vectors.rotate2d(self.rotation_angle, v) for v in self.points]
        return [vectors.add((self.x, self.y), v) for v in rotated]
    
    def move(self, milliseconds):
        dx, dy = self.vx * milliseconds / 1000.0, self.vy * milliseconds / 1000
        self.x, self.y = vectors.add((self.x, self.y), (dx, dy))
        if self.x < -10:
            self.x += 20
        if self.y < -10:
            self.y += 20
        if self.x > 10:
            self.x -= 20
        if self.y > 10:
            self.y -= 20
        self.rotation_angle += self.angular_velocity * milliseconds / 1000
    
    def segments(self):
        point_count = len(self.points)
        points = self.transformed()
        return [(points[i], points[(i+1)%point_count]) for i in range(0, point_count)]

    def does_collide(self, other_poly):
        for other_segment in other_poly.segments():
            if self.does_intersect(other_segment):
                return True
        return False

    def does_intersect(self, other_segment):
        for segment in self.segments():
            if do_segments_intersect(other_segment, segment):
                return True
        return False

class Ship(PolygonModel):
    def __init__(self):
        super().__init__([(0.5, 0), (-0.25, 0.25), (-0.25, -0.25)])
    
    def laser_segment(self):
        dist = 20 * sqrt(2)
        x, y = self.transformed()[0]
        return (x, y), (x+dist*cos(self.rotation_angle), y + dist * sin(self.rotation_angle))

class Asteroid(PolygonModel):
    def __init__(self):
        sides = randint(5, 9)
        vs = [vectors.to_cartesian((uniform(0.5, 1.0), 2 * pi * i / sides))
            for i in range(0, sides)]
        super().__init__(vs)
        self.vx = uniform(-1, 1)
        self.vy = uniform(-1, 1)
        self.angular_velocity = uniform(-pi/2, pi/2)




BLACK = (0, 0, 0)
WHITE = (255, 255, 255)
BLUE =  (0, 0, 255)
GREEN = (0, 255, 0)
RED =   (255, 0, 0)
LIGHT_GRAY =  (240, 240, 240)
DARK_GRAY = (128, 128, 128)
width, height = 400, 400

def to_pixel(x, y):
    return (width / 2 + width * x / 20, height/2 - height * y / 20)

def draw_poly(screen, polygon_model, color=BLACK):
    pixel_points = [to_pixel(x, y) for x, y in polygon_model.transformed()]
    pygame.draw.lines(screen, color, True, pixel_points, 2)
    if polygon_model.draw_center:
        cx, cy = to_pixel(polygon_model.x, polygon_model.y)
        pygame.draw.circle(screen, BLACK, (int(cx), int(cy)), 4, 4)
    
def draw_segment(screen, v1, v2, color=RED):
    pygame.draw.line(screen, color, to_pixel(*v1), to_pixel(*v2), 2)

def draw_grid(screen):
    for x in range(-9, 10):
        draw_segment(screen, (x, -10), (x, 10), color=LIGHT_GRAY)
    for y in range(-9, 10):
        draw_segment(screen, (-10, y), (10, y), color=LIGHT_GRAY)
    
    draw_segment(screen, (-10, 0), (10, 0), color=DARK_GRAY)
    draw_segment(screen, (0, -10), (0, 10), color=DARK_GRAY)


ship = Ship()
asteroid_count = 10
default_asteroids = [Asteroid() for _ in range(0, asteroid_count)]
for ast in default_asteroids:
    ast.x = randint(-9, 9)
    ast.y = randint(-9, 9)


acceleration = 3
def main(asteroids=default_asteroids):
    pygame.init()
    screen = pygame.display.set_mode([width, height])
    done = False
    clock = pygame.time.Clock()

    while not done:
        clock.tick()

        for event in pygame.event.get():
            if event.type == pygame.QUIT:
                done = True
        milliseconds = clock.get_time()
        keys = pygame.key.get_pressed()

        for ast in asteroids:
            ast.move(milliseconds)
        
        if keys[pygame.K_LEFT]:
            ship.rotation_angle += milliseconds * (2 * pi / 1000)
        if keys[pygame.K_RIGHT]:
            ship.rotation_angle -= milliseconds * (2 * pi / 1000)
        if keys[pygame.K_UP]:
            ax = acceleration * cos(ship.rotation_angle)
            ay = acceleration * sin(ship.rotation_angle)
            ship.vx += ax * milliseconds/1000
            ship.vy += ay * milliseconds/1000
        if keys[pygame.K_DOWN]:
            ax = acceleration * cos(ship.rotation_angle)
            ay = acceleration * sin(ship.rotation_angle)
            ship.vx -= ax * milliseconds/1000
            ship.vy -= ay * milliseconds/1000
        ship.move(milliseconds)

        laser = ship.laser_segment()
        screen.fill(WHITE)
        draw_grid(screen)
        if keys[pygame.K_SPACE]:
            draw_segment(screen, *laser)
        draw_poly(screen, ship)
        for ast in asteroids:
            if keys[pygame.K_SPACE] and ast.does_intersect(laser):
                asteroids.remove(ast)
            else:
                draw_poly(screen, ast, color=GREEN)
        pygame.display.flip()
    pygame.quit()

if __name__ == "__main__":
    main()

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

ㅡㅇㅂ재ㅔㄹㅈ배렂벫재

잠깐 딴짓했다가 다시 할 기분이 잘 안들었는데 겨우 시작

 

 

유량 다루기 시작 - 먼저 그래프 드로잉 함수부터

plot function은 이름 그대로 함수와 시간 구간 주면 그대로 드로잉 해주는 함수

여기서 사용한 함수는

유체량과 유체흐름률 두개가 있음

유체량은 (t-4)^3 / 64 + 3.3

유체흐름률은 유체량 함수의 도함수 3 * (t - 4)^2 / 64

0 ~ 10초 사이 유체량 함수를 드로잉한 결과는 다음과 같음

import matplotlib.pyplot as plt
import numpy as np

def plot_function(f, tmin, tmax, tlabel=None, xlabel=None, axes=False, **kwargs):
    ts = np.linspace(tmin, tmax, 1000)
    if tlabel:
        plt.xlabel(tlabel, fontsize=18)
    if xlabel:
        plt.ylabel(xlabel, fontsize=18)
    plt.plot(ts, [f(t) for t in ts], **kwargs)
    if axes:
        total_t = tmax-tmin
        plt.plot([tmin-total_t,tmax+total_t/10], [0, 0], c='k', linewidth=1)
        plt.xlim(tmin-total_t/10,tmax+total_t/10)
        xmin, xmax = plt.ylim()
        plt.plot([0, 0], [xmin, xmax], c='k', linewidth=1)
        plt.ylim(xmin, xmax)
def plot_volume(f,tmin,tmax,axes=False,**kwargs):
    plot_function(f, tmin, tmax, tlabel='time (hr)', xlabel='volume (bbl)', axes=axes, **kwargs)
def plot_flow_rate(f, tmin, tmax, axes=False, **kwargs):
    plot_function(f, tmin, tmax, tlabel='time (hr)', xlabel='flow rate(bbl/hr)', axes=axes, **kwargs)

def volume(t):
    return (t-4)**3 / 64 + 3.3
def flow_rate(t):
    return 3 * (t - 4) ** 2 / 64

plot_volume(volume, 0, 10)

 

 

특정 시간에서 유체의 량과 유체 량의 평균 변화를 구해보자

유체량에 대한 함수 volume에다가 시간 4, 9일때 량을 구하면 3.3, 5.2 가 나옴

평균 유체 흐름량 함수를 구현해보면 0.39가 나온다.

 

 

 

def average_flow_rate(v, t1, t2):
    return (v(t2) - v(t1))/(t2 - t1)

print(volume(4))
print(volume(9))
print(average_flow_rate(volume, 4, 9))

 

 

음의 변화율 다루기

아까 구현한 볼륨 함수는 점점 증가하는 형태였으나 이번엔

시간에 따라 내려가도록 구현

0 ~ 4 초사이 평균 흐름율은 -0.8

플로팅 결과는 저렇다

 

def decreasing_volume(t):
    if t < 5:
        return 10 - (t ** 2) / 5
    else:
        return 0.2 * (10 - t) ** 2

print(average_flow_rate(decreasing_volume, 0, 4))
plot_volume(decreasing_volume, 0, 10)

 

 

시컨트 라인 그리기

12시에 오도메트리가 77600 마일인데, 4시 30분에는 77900 마일이 됬다. 평균 속도 변화율은 어떨까?

으악 삼각함수도 잘 모르는데 시컨트라니

잠깐 찾아보니 삼각함수의 역수란다 이게 왜쓰는지 모르겠는데

코사인은 밑변/빗변이니 시컨트는 빗변/밑변 

흐름은 잘 모르겠지만 대충 시컨트 라인은 x1, x2사이 주어진 함수에 대한 직선 방정식을 구하는 함순거같음

 

def secant_line(f, x1, x2):
    def line(x):
        return f(x1) + (x - x1) * (f(x2) - f(x1)) / (x2 - x1)
    return line

def plot_secant(f, x1, x2, color='k'):
    line = secant_line(f, x1, x2)
    plot_function(line,x1,x2,c=color)
    plt.scatter([x1, x2], [f(x1), f(x2)], c=color)

plot_volume(volume, 0, 10)
plot_secant(volume, 4, 8)

 

 

 

 

 

 

구간별 평균 변화율

이번에는 1초 간격 평균 변화율을 구해보자

코드는 여전히 단순

def interval_flow_rate(v, t1, t2, dt):
    return [(t, average_flow_rate(v, t, t+dt)) for t in np.arange(t1, t2, dt)]

interval_flow_rate(volume, 0, 10, 1)

 

 

 

 

구간별 평균 유체 변화율 산점도 그리기

방금 작성한 구간별 평균 변화율로 산점도 드로잉하는 함수

 

def plot_interval_flow_rates(volume, t1, t2, dt):
    series = interval_flow_rate(volume, t1, t2, dt)
    times = [t for (t, _) in series]
    rates = [q for (_, q) in series]
    plt.scatter(times, rates)

plot_interval_flow_rates(volume, 0, 10, 1)

 

 

 

유체 흐름률과 구간별 평균 흐름률 플로팅

유체 흐름률은 유체량 함수의 도함수로 부터 구함

구간별 흐름률은 유체량에서 구간별 평균 변화율로 구함

모양은 비슷하나 후자가 살짝 뒤처짐

 

plot_flow_rate(flow_rate, 0, 10)
plot_interval_flow_rates(volume, 0, 10, 1)

 

 

하지만 구간을 1 -> 1/3으로 줄이면 원 함수에 더 가까워짐

 

plot_flow_rate(flow_rate, 0, 10)
plot_interval_flow_rates(volume, 0, 10, 1/3)

 

 

감소 하는 함수의 구간별 평균 흐름율 플로팅하면 이럼

plot_interval_flow_rates(decreasing_volume, 0, 10, 0.5)

 

직선의 경우 구간별 평균 흐름률

def linear_volume_function(t):
    return 5 * t + 3

plot_interval_flow_rates(linear_volume_function, 0, 10, 0.25)

 

 

 

유체량의 평균 변화율은 1에 가까워질수록 어떻게 될까?

아까 구현한 평균 유체 변화율로 1에 가까워지게 해보면 특정 값에 점점 가까워진다.

print(average_flow_rate(volume, 0.5, 1.5))
print(average_flow_rate(volume, 0.9, 1.1))
print(average_flow_rate(volume, 0.99, 1.01))
print(average_flow_rate(volume, 0.999, 1.001))
print(average_flow_rate(volume, 0.9999, 1.0001))

 

 

 

 

 

 

 

시컨트 라인과 탄젠트

계속 시컨트 시컨트하길래 잠깐 딴데 찾아보니 이렇게 설명한다

시컨트라인은 곡선 두 점 지나는 선

접선은 곡선의 한점을 지나는 선

https://ko.strephonsays.com/chord-and-vs-secant-and-vs-tangent-1750

 

 

 

순간 변화율 함수 구현하기

아깐 평균 변화율 구간을 1을 중심으로 좁힘으로서 특정 수에 수렴하는 것을 봤는데

이번에는 아예 순간변화율을 구하는 함수 구현됨

먼저 톨레랑스 용인도? 얼마까지 에러 봐줄래 하는 변수를 잡아줌

디폴트로 디지트 6했으니 10 ^ -6까지 오차는 용인한다.

h = 1로 놓고 위아래 1차이와 평균 변화율 구함

2 * digit 그러니까 12회 까지 반복수행하는대 오차가 10 ^-6 보다 크지않으면 반올림하고 반환

루프마다 h는 1/10배로 줄어듬

def instantaneous_flow_rate(v, t, digits=6):
    tolerance = 10 ** (-digits)
    h = 1
    approx = average_flow_rate(v, t-h, t+h)
    for i in range(0, 2 * digits):
        h = h / 10
        next_approx = average_flow_rate(v, t-h, t+h)
        if abs(next_approx - approx) < tolerance:
            return round(next_approx, digits)
        else:
            approx = next_approx
    raise Exception("derivative did not converge")

instantaneous_flow_rate(volume, 1)

 

 

 

 

구현한 순간 변화율 함수가 옳은지 플로팅 해보기

먼저 기존에 볼륨 도함수인 변화율 함수 플로팅

이번엔 순간 변화율 구하는 함수에 볼륨 함수를 집어넣은 뒤 플로팅

둘다 동일함

 

 

def get_flow_rate_function(v):
    def flow_rate_function(t):
        return instantaneous_flow_rate(volume, t)
    return flow_rate_function

plot_function(flow_rate, 0, 10)

plot_function(get_flow_rate_function(volume),0,10)

 

 

짧은 시간 동안 유체량 변화 정도 

1초 동안 유체량 변화 정도 계산시  0.1875나 실제로는 1.09가 나온다.

대신 구간을 1초 대신 0.01로 줄여보면 0.001875 실제로는 0.00186으로

구간을 좁히면 실제 값에 근사하게 나옴

 

def small_volume_change(q, t, dt):
    return q(t) * dt

print(small_volume_change(flow_rate, 2, 1))
print(volume(3) - volume(2))
print("#######")
print(small_volume_change(flow_rate, 2, 0.01))
print(volume(2.01) - volume(2))

 

 

 

볼륨 변화량 구하기

계속 유체라고 했는데 변수 명 그대로 볼륨이라고 해야겠다.

0 ~ 10 초간 0.1초 간격으로 볼륨 변화량을 합치면 4.3289가 나온다.

하지만 실제 볼륨10, 볼륨0으로 구해보면 4.375가 뜬다.

이번엔 구간을 0.0001로 줄이면 4.3749로 실제 값에 가까워짐

def volume_change(q, t1, t2, dt):
    return (sum(small_volume_change(q, t, dt) for t in np.arange(t1, t2, dt)))

print(volume_change(flow_rate, 0, 10, 0.1))
print(volume(10) - volume(0))
print(volume_change(flow_rate, 0, 10, 0.0001))

 

 

볼륨 플로팅하기

이번엔 실제 볼륨 함수와 근사화한 볼륨 함수를 플로팅함

볼륨 함수 근사화는 아까 구현한 볼륨 변화 사용

근사시 구간을 0.5로 해서 좀 안맞아보이긴 한데 거의 비슷하게 나오긴 함

 

근사 구간을 0.1로하면 아까 보다 원함수에 가까워짐

 

구간을 0.01로 더 줄이면 더가까워짐

def approximate_volume(q, v0, dt, T):
    return v0 + volume_change(q, 0, T, dt)

def approximate_volume_function(q, v0, dt):
    def volume_function(T):
        return approximate_volume(q, v0, dt, T)
    return volume_function

plot_volume(approximate_volume_function(flow_rate, 2.3, 0.5), 0, 10)
plot_volume(volume, 0, 10)



plot_volume(approximate_volume_function(flow_rate, 2.3, 0.1), 0, 10)
plot_volume(volume, 0, 10)




plot_volume(approximate_volume_function(flow_rate, 2.3, 0.01), 0, 10)
plot_volume(volume, 0, 10)

 

 

 

 

 

 

원래 같으면 마음에 드는 영상 보고난지 얼마 안되서 적긴 해야되는데

잠깐 다른거 하느라 흐름을 놓침

 

이제와서 쓰자니 잘 기억도 안나고 그렇지만

다시 흐름 잡기 준비하고자 작성

 

내용은 이전에 해석학 공부하면서 봤던 사실들을 반복했지만

설명은 훨씬 간단하고 듣기 담백했었다.

 

특히 초월함수는 대수 방정식으로 표현못하는 경우 초월함수랬는데

테일러 급수는 초월함수를 방정식으로 표현 가능했다고한다.

그게 어떻게 되나 싶었는데 진짜 급수 형태로 표현해냈더라

 

연속 미분 가능하다는 의미도 한번 더보니까

여전히 잘 모르긴 하지만 뭐 안본것보단 나았고

4분짜리라 다시 보고 적어도 되긴한데 봣던거 또보고 하긴 귀찬

 

 

 

 

+ Recent posts