TUKR実装-福永

Lecture_TUKR/TUKR.py

@@ -0,0 +1,181 @@
+import numpy as np
+# from tqdm import tqdm #プログレスバーを表示させてくれる
+import jax,jaxlib
+import jax.numpy as jnp
+
+
+class UKR:
+    def __init__(self, X, latent_dim, sigma1, sigma2, prior='random', Uinit=None, Vinit=None):
+        #--------初期値を設定する．---------
+        self.X = X
+        #ここから下は書き換えてね
+        self.xsamples = X.shape[0]
+        self.ysamples = X.shape[1]
+        self.ob_dim = X.shape[2]
+        self.sigma1 = sigma1
+        self.sigma2 = sigma2
+        self.latent_dim = latent_dim
+
+
+        if Uinit is None:
+            if prior == 'random': #一様事前分布のとき
+                self.U = np.random.uniform(low=0.0, high=1.0, size=(self.xsamples, self.latent_dim))
+
+
+            else: #ガウス事前分布のとき
+                self.U = np.random.normal(0, self.sigma1, (self.xsamples, self.latent_dim))
+
+        else: #Zの初期値が与えられた時
+            self.U = Uinit
+        if Vinit is None:
+            if prior == 'random': #一様事前分布のとき
+                self.V = np.random.normal(0, self.sigma2, (self.ysamples, self.latent_dim))
+
+            else: #ガウス事前分布のとき
+                self.V = np.random.normal(self.ysamples*self.latent_dim).reshape(self.ysamples, self.latent_dim)
+        else: #Zの初期値が与えられた時
+            self.V = Vinit
+
+        self.history = {}
+    def f(self, U, V):
+        Ud = np.sum((U[:, None, :] - self.U[None, :, :])**2, axis=2)
+        UH = -1*(Ud/(2*self.sigma1**2))
+        Uh = jnp.exp(UH)
+
+        Vd = np.sum((V[:, None, :] - self.V[None, :, :]) ** 2, axis=2)
+        VH = -1 * (Vd / (2 * self.sigma2 ** 2))
+        Vh = jnp.exp(VH)
+
+        bunshi = jnp.einsum('KI,MJ,IJD -> KMD', Uh, Vh, self.X)
+
+
+        bunbo = jnp.einsum('KI,MJ -> KM', Uh, Vh)
+        newbunbo = bunbo[:, :, None]
+
+
+
+        f = bunshi/newbunbo
+        # print(U.shape)
+        #写像の計算
+        return f
+
+
+    def E(self, U, V,  X, alpha, norm): #目的関数の計算
+        E = jnp.sum((X - self.f(U, V))**2)
+        UR = alpha * jnp.sum(jnp.abs(U ** norm))
+        VR = alpha * jnp.sum(jnp.abs(V ** norm))
+        E = E/self.xsamples/self.ysamples + UR/self.xsamples + VR/self.ysamples
+
+        return E
+
+
+
+    def fit(self, nb_epoch: int, eta: float, alpha: float, norm: float):
+        # 学習過程記録用
+        self.history['u'] = np.zeros((nb_epoch, self.xsamples, self.latent_dim))
+        self.history['v'] = np.zeros((nb_epoch, self.ysamples, self.latent_dim))
+        self.history['f'] = np.zeros((nb_epoch, self.xsamples, self.ysamples, self.ob_dim))
+        self.history['error'] = np.zeros(nb_epoch)
+
+        for epoch in range(nb_epoch):
+            dEdu = jax.grad(self.E, argnums=0)(self.U, self.V, self.X, alpha, norm)
+            self.U = self.U - eta * dEdu
+            dEdv = jax.grad(self.E, argnums=1)(self.U, self.V, self.X, alpha, norm)
+            self.V = self.V - eta * dEdv
+
+           # U,Vの更新
+
+
+
+
+            # 学習過程記録用
+            self.history['u'][epoch] = self.U
+            self.history['v'][epoch] = self.V
+            self.history['f'][epoch] = self.f(self.U, self.V)
+            self.history['error'][epoch] = self.E(self.U, self.V, self.X, alpha, norm)
+
+    #--------------以下描画用(上の部分が実装できたら実装してね)---------------------
+    def calc_approximate_fu(self, resolution): #fのメッシュ描画用，resolution:一辺の代表点の数
+        nb_epoch = self.history['u'].shape[0]
+        self.history['y'] = np.zeros((nb_epoch, resolution ** self.latent_dim, self.ob_dim))
+        for epoch in range(nb_epoch):
+            uzeta = self.create_uzeta(self.U)
+            vzeta = self.create_vzeta(self.V)
+            Y = self.f(uzeta, vzeta)
+            # print(self.history['y'][epoch])
+            self.history['y'][epoch] = Y
+        return self.history['y']
+
+    # def calc_approximate_fv(self, resolution): #fのメッシュ描画用，resolution:一辺の代表点の数
+    #     nb_epoch = self.history['v'].shape[0]
+    #     self.history['y'] = np.zeros((nb_epoch, resolution ** self.latent_dim, self.ob_dim))
+    #     for epoch in range(nb_epoch):
+    #         vzeta = create_vzeta(self.V, resolution)
+    #         Y = self.f(vzeta, self.history['v'][epoch])
+    #         self.history['y'][epoch] = Y
+    #     return self.history['y']
+
+
+    def create_uzeta(self, Z): #fのメッシュの描画用に潜在空間に代表点zetaを作る．
+
+        u_y = np.linspace(np.amin(Z), np.amax(Z), self.xsamples).reshape(-1, 1)
+        # UYY = np.meshgrid(u_y)
+        # A = np.amin(Z)
+        # print(self.U.shape)
+
+        # zeta = np.concatenate([uxx[:, None], uyy[:, None]], axis=1)
+        return u_y
+
+    def create_vzeta(self, Z): #fのメッシュの描画用に潜在空間に代表点zetaを作る．
+
+        v_y = np.linspace(np.min(Z), np.max(Z), self.ysamples).reshape(-1, 1)
+        # UYY = np.meshgrid(u_y)
+
+
+        # zeta = np.concatenate([uxx[:, None], uyy[:, None]], axis=1)
+        return v_y
+
+# def create_vzeta(V, resolution):
+#     v_x = np.linspace(np.min(V), np.max(V), resolution)
+#     v_y = np.linspace(np.min(V), np.max(V), resolution)
+#     VXX, VYY = np.meshgrid(v_x, v_y)
+#     vxx = VXX.reshape(-1)
+#     vyy = VYY.reshape(-1)
+#
+#     vzeta = np.concatenate([vxx[:, None], vyy[:, None]], axis=1)
+#
+#     return vzeta
+
+if __name__ == '__main__':
+    from Lecture_UKR.data import create_kura
+    from Lecture_UKR.data import create_rasen
+    # from Lecture_UKR.data import create_2d_sin_curve
+    from Lecture_TUKR.visualizer import visualize_history
+    from Lecture_TUKR.data_scratch import load_kura_tsom
+
+    #各種パラメータ変えて遊んでみてね．
+    epoch = 100 #学習回数
+    sigma1 = 0.5 #カーネルの幅
+    sigma2 = 0.5
+    eta = 5 #学習率
+    latent_dim = 1 #潜在空間の次元
+    alpha = 0.001
+    norm = 2
+    seed = 4
+    np.random.seed(seed)
+
+    #入力データ（詳しくはdata.pyを除いてみると良い）
+    xsamples = 20 #データ数
+    ysamples = 10
+    X = load_kura_tsom(xsamples, ysamples) #鞍型データ　ob_dim=3, 真のL=2
+
+    #X = create_rasen(nb_samples) #らせん型データ　ob_dim=3, 真のL=1
+    # X = create_2d_sin_curve(nb_samples) #sin型データ　ob_dim=2, 真のL=1
+    ukr = UKR(X, latent_dim, sigma1, sigma2, prior='random')
+    ukr.fit(epoch, eta, alpha, norm)
+    #visualize_history(X, ukr.history['f'], ukr.history['z'], ukr.history['error'], save_gif=False, filename="tmp")
+
+    #----------描画部分が実装されたらコメントアウト外す----------
+    # ukr.calc_approximate_fu(resolution=10)
+    # ukr.calc_approximate_fv(resolution=10)
+    visualize_history(X, ukr.history['f'], ukr.history['u'], ukr.history['v'], ukr.history['error'], save_gif=True, filename="/Users/furukawashuushi/Desktop/3ヶ月コースGIF/TUKR")

Lecture_TUKR/data_scratch.py

@@ -2,19 +2,18 @@
 import matplotlib.pyplot as plt
 
 def load_kura_tsom(xsamples, ysamples, missing_rate=None,retz=False):
-    z1 = 
-    z2 = 
+    z1 = np.linspace(-1, 1, xsamples)
+    z2 = np.linspace(-1, 1, ysamples)
 
-    z1_repeated, z2_repeated = 
-    x1 = 
-    x2 = 
-    x3 = 
+    z1_repeated, z2_repeated = np.meshgrid(z1, z2, indexing = 'ij')
+    x1 = z1_repeated #+ np.random.normal(loc=0.0, scale=0.1, size=(xsamples, ysamples))
+    x2 = z2_repeated
+    x3 = x1**2 - x2**2
     #ノイズを加えたい時はここをいじる,locがガウス分布の平均、scaleが分散,size何個ノイズを作るか
     #このノイズを加えることによって三次元空間のデータ点は上下に動く
 
-    x = 
-    truez = 
-
+    x = np.concatenate((x1[:, :, np.newaxis],x2[:, :, np.newaxis], x3[:, :, np.newaxis]), axis=2)
+    truez = np.concatenate((x1[:, :, np.newaxis],x2[:, :, np.newaxis], x3[:, :, np.newaxis]), axis=2)
     if missing_rate == 0 or missing_rate == None:
         if retz:
             return x, truez
@@ -25,13 +24,13 @@ def load_kura_tsom(xsamples, ysamples, missing_rate=None,retz=False):
     import matplotlib.pyplot as plt
     from mpl_toolkits.mplot3d import Axes3D
 
-    xsamples = 
-    ysamples = 
+    xsamples = 10
+    ysamples = 5
+
+    x, truez = load_kura_tsom(xsamples, ysamples, retz=True)
 
-    x, truez = load_kura_tsom() 
-
     fig = plt.figure(figsize=[5, 5])
     ax_x = fig.add_subplot(projection='3d')
-    ax_x.scatter()
+    ax_x.scatter(x[:, :, 0].flatten(), x[:, :, 1].flatten(), x[:, :, 2].flatten())
     ax_x.set_title('Generated three-dimensional data')
     plt.show()

Lecture_TUKR/visualizer.py

@@ -0,0 +1,101 @@
+import numpy as np
+from matplotlib import pyplot as plt
+from matplotlib.animation import FuncAnimation
+
+STEP = 150
+
+
+def visualize_history(X, Y_history, U_history, V_history, error_history, save_gif=False, filename="tmp"):
+    input_dim, latent_dim1, latent_dim2 = X.shape[2], U_history[0].shape[1], V_history[0].shape[1]
+    input_projection_type = '3d' if input_dim > 2 else 'rectilinear'
+
+    fig = plt.figure(figsize=(10, 8))
+    gs = fig.add_gridspec(3, 3)
+    input_ax = fig.add_subplot(gs[0:2, 0], projection=input_projection_type)
+    latent_ax1 = fig.add_subplot(gs[0:2, 1], aspect='equal')
+    latent_ax2 = fig.add_subplot(gs[0:2, 2], aspect='equal')
+    error_ax = fig.add_subplot(gs[2, :])
+    num_epoch = len(Y_history)
+
+    # if input_dim == 3 and latent_dim1 == 2:
+    #     resolution = int(np.sqrt(Y_history.shape[1]))
+    #     if Y_history.shape[1] == resolution ** 2:
+    #         Y_history = np.array(Y_history).reshape((num_epoch, resolution, resolution, input_dim))
+    #
+    # if input_dim == 3 and latent_dim2 == 2:
+    #     resolution = int(np.sqrt(Y_history.shape[1]))
+    #     if Y_history.shape[1] == resolution ** 2:
+    #         Y_history = np.array(Y_history).reshape((num_epoch, resolution, resolution, input_dim))
+    observable_drawer = [None, None, draw_observable_2D,
+                         draw_observable_3D][input_dim]
+
+    latent_drawer1 = [None, draw_latent_1D, draw_latent_2D][latent_dim1]
+    latent_drawer2 = [None, draw_latent_1D, draw_latent_2D][latent_dim2]
+
+    ani = FuncAnimation(
+        fig,
+        update_graph,
+        frames=num_epoch,  # // STEP,
+        repeat=True,
+        interval=50,
+        fargs=(observable_drawer, latent_drawer1, latent_drawer2, X, Y_history, U_history, V_history, error_history, fig,
+               input_ax, latent_ax1, latent_ax2, error_ax, num_epoch))
+    plt.show()
+    if save_gif:
+        ani.save(f"{filename}.gif", writer='ffmpeg')
+
+
+def update_graph(epoch, observable_drawer, latent_drawer1,latent_drawer2, X, Y_history,
+                 U_history, V_history, error_history, fig, input_ax, latent_ax1, latent_ax2, error_ax, num_epoch):
+    fig.suptitle(f"epoch: {epoch}")
+    input_ax.cla()
+    #  input_ax.view_init(azim=(epoch * 400 / num_epoch), elev=30)
+    latent_ax1.cla()
+    latent_ax2.cla()
+    error_ax.cla()
+
+    Y, U, V= Y_history[epoch], U_history[epoch], V_history[epoch]
+    colormap = X[:, :, 0]
+    colormap1 = U[:, 0]
+    colormap2 = V[:, 0]
+    # print(X.shape)
+
+    observable_drawer(input_ax, X, Y, colormap)
+    latent_drawer1(latent_ax1, U,  colormap1)
+    latent_drawer2(latent_ax2, V, colormap2)
+    draw_error(error_ax, error_history, epoch)
+
+
+def draw_observable_3D(ax, X, Y, colormap):
+    # print(X.shape,type(colormap))
+    # print(colormap.shape)
+    ax.scatter(X[:, :, 0], X[:, :, 1], X[:, :, 2], c=colormap)
+    # ax.set_zlim(-1, 1)
+    if len(Y.shape) == 3:
+        ax.plot_wireframe(Y[:, :, 0], Y[:, :, 1], Y[:, :, 2], color='black')
+        # ax.scatter(Y[:, :, 0], Y[:, :, 1], Y[:, :, 2], color='black')
+    else:
+        ax.plot(Y[:, 0], Y[:, 1], Y[:, 2], color='black')
+# ax.plot(Y[:, 0], Y[:, 1], Y[:, 2], color='black')
+# ax.plot_wireframe(Y[:, :, 0], Y[:, :, 1], Y[:, :, 2], color='black')
+
+
+def draw_observable_2D(ax, X, Y, colormap):
+    ax.scatter(X[:, 0], X[:, 1], c=colormap)
+    ax.plot(Y[:, 0], Y[:, 1], c='black')
+
+
+def draw_latent_2D(ax, Z, colormap):
+    ax.set_xlim(-1.1, 1.1)
+    ax.set_ylim(-1.1, 1.1)
+    ax.scatter(Z[:, 0], Z[:, 1], c=colormap)
+
+
+def draw_latent_1D(ax, Z, colormap):
+    ax.scatter(Z, np.zeros(Z.shape), c=colormap)
+    ax.set_ylim(-1, 1)
+
+def draw_error(ax, error_history, epoch):
+    ax.set_title("error_function", fontsize=8)
+    ax.plot(error_history, label='誤差関数')
+    ax.scatter(epoch, error_history[epoch], s=55, marker="*")

Lecture_UKR/data.py

@@ -47,4 +47,30 @@ def create_2d_sin_curve(nb_samples, noise_scale=0.01):
     ax.set_ylabel("x2")
     ax.set_xlabel("x3")
 
-    plt.show()
+    plt.show()
+
+
+    # 機械学習の課題
+    # import random
+    # A=0
+    # for a in range(1000000):
+    #     d = random.randint(1, 6)
+    #     if (d == 1) | (d == 2):
+    #         A+=1
+    # print(str(100*A/1000000)+"%")
+    #
+    # bag = 2  # 袋の総数
+    # bag_x = 1  # 袋Xの数
+    # bag_y = 1  # 袋Yの数
+    #
+    # ball_in_bag_x = [4, 6]  # 袋Xの中身(赤玉が4個、白玉が6個)
+    # ball_in_bag_y = [5, 2]  # 袋Yの中身(赤玉が5個、白玉が2個)
+    #
+    #
+    # def jouhou(bag, bag_sum, ball_in_bag):
+    #     ans = (ball_in_bag[0] / sum(ball_in_bag)) * (bag / bag_sum)
+    #     return ans
+    #
+    #
+    # j = jouhou(bag_x, bag, ball_in_bag_x)
+    # print(str(j*100) + "%")