BooleanNetwork.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import cm\n",
    "%matplotlib inline\n",
    "\n",
    "\n",
    "import torch\n",
    "from torch import nn\n",
    "from torch import distributions\n",
    "from torch.nn.parameter import Parameter\n",
    "\n",
    "import models\n",
    "from models import Renorm_Dynamic\n",
    "from EI_calculation import approx_ei\n",
    "use_cuda = torch.cuda.is_available()\n",
    "device = torch.device('cuda:0') if use_cuda else torch.device('cpu')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Boolean Network"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class BnSim(object):\n",
    "    def __init__(self, number_of_vertice):\n",
    "        # 初始化节点数量\n",
    "        self.number_of_vertice = number_of_vertice\n",
    "        self.edges = self.gen_edge(False, [], 0, '4_vertice_bn')\n",
    "\n",
    "    def gen_edge(self, use_group_list, group_list, random_del_edges_ratio, adj_template_name):\n",
    "        \n",
    "        edges = np.zeros((self.number_of_vertice, self.number_of_vertice))\n",
    "        # 初始化邻接矩阵 所有连边均为1\n",
    "        \n",
    "        if use_group_list == True:\n",
    "            #使用分组列表    \n",
    "            for i in range(len(group_list)-1):\n",
    "                for j in range(self.number_of_vertice):\n",
    "                    for k in range(self.number_of_vertice):\n",
    "                        if i < len(group_list)-2:\n",
    "                            if j>=group_list[i] and j<group_list[i+1] and k>=group_list[i+1] and k<group_list[i+2]:\n",
    "                                # 依照分组列表删除不需要的连边\n",
    "                                if np.random.rand() > random_del_edges_ratio:\n",
    "                                    edges[j][k] = 1.\n",
    "                        else:\n",
    "                            if j >=group_list[i] and j<group_list[i+1] and k<group_list[1]:\n",
    "                                if np.random.rand() > random_del_edges_ratio:\n",
    "                                    edges[j][k] = 1.\n",
    "            print(edges)\n",
    "        elif use_group_list == False:\n",
    "            if adj_template_name == \"4_vertice_bn\":\n",
    "                self.number_of_vertice = 4\n",
    "                edges = np.zeros((self.number_of_vertice, self.number_of_vertice))\n",
    "                edges[0][2] = 1\n",
    "                edges[0][3] = 1\n",
    "                edges[1][2] = 1\n",
    "                edges[1][3] = 1\n",
    "                edges[2][0] = 1\n",
    "                edges[3][0] = 1\n",
    "                edges[2][1] = 1\n",
    "                edges[3][1] = 1\n",
    "                \n",
    "        \n",
    "        return edges\n",
    "\n",
    "    def sample_one_step(self, state_now, edges):\n",
    "        # 返回复杂系统下一时刻的状态\n",
    "        state_next = torch.zeros(self.number_of_vertice)\n",
    "        # 初始化复杂系统下一时刻的状态\n",
    "        for i in range(self.number_of_vertice):\n",
    "            # 对每个节点的状态进行遍历\n",
    "            n_prob = 0.\n",
    "            p_prob = 0.\n",
    "            degree = 0\n",
    "            # 初始化节点的邻居中状态为1的数量、状态为0的数量以及节点的度\n",
    "            for j in range(self.number_of_vertice):\n",
    "                # 遍历节点的所有一阶邻居，计算上述三个值：n_prob、p_prob、degree\n",
    "                if edges[i][j] == 1.:\n",
    "                    if state_now[j] == 0.:\n",
    "                        n_prob += 1\n",
    "                    elif state_now[j] == 1.:\n",
    "                        p_prob += 1\n",
    "                    degree += 1\n",
    "            if p_prob < degree:\n",
    "                #如果不是所有的邻居状态都为一，则该节点下一时刻有0.7的概率为0、0.3的概率为1\n",
    "                state_next[i] = torch.multinomial(torch.Tensor([0.7,0.3]), 1)\n",
    "            else:\n",
    "                # 如果所有的邻居状态都为1，则该节点下意识的状态一定为1\n",
    "                state_next[i] = 1.\n",
    "        return state_next\n",
    "        #返回下一时刻的状态\n",
    "\n",
    "    def sample1(self, number_of_data):\n",
    "        # 生成动力学时间序列\n",
    "        # number_of_data为生成的时间序列个数\n",
    "        # T为时间序列的长度，一般设置为2\n",
    "        # group_list为分组列表，用于指导生成邻接矩阵\n",
    "        \n",
    "        # 调用gen_edge函数生成邻接矩阵\n",
    "        n = self.number_of_vertice\n",
    "        # n为节点数量，若adj_template_name与初始化sim时的节点数目不符合，则会在self.gen_edge中更新self.number_of_vertice\n",
    "        state = torch.zeros((number_of_data, 2**n))\n",
    "        state_next = torch.zeros([number_of_data, 2**n])\n",
    "        # 初始化时间序列信息，此时各个时间序列中各个节点在各个时刻的状态均为0\n",
    "        for num in range(number_of_data):\n",
    "            #使用循环对各个时间序列进行操作\n",
    "            rands = torch.multinomial(torch.Tensor([0.5, 0.5]), n, replacement=True)\n",
    "            state[num, :]=self.encoding(rands)*1.0\n",
    "            # 随机生成初始状态，状态为0和1的概率均为50%\n",
    "            for i in range(1):\n",
    "            # 从初始状态开始向后演化(T-1)个时刻，该过程会调用sample_one_step函数\n",
    "                sn = self.sample_one_step(rands, edges=self.edges)\n",
    "                state_next[num, :] = self.encoding(sn)*1.0\n",
    "        return state, state_next, self.edges\n",
    "    def sample(self, number_of_data):\n",
    "        # 生成动力学时间序列\n",
    "        # number_of_data为生成的时间序列个数\n",
    "        # T为时间序列的长度，一般设置为2\n",
    "        # group_list为分组列表，用于指导生成邻接矩阵\n",
    "        \n",
    "        # 调用gen_edge函数生成邻接矩阵\n",
    "        n = self.number_of_vertice\n",
    "        # n为节点数量，若adj_template_name与初始化sim时的节点数目不符合，则会在self.gen_edge中更新self.number_of_vertice\n",
    "        state = torch.zeros((number_of_data, n))\n",
    "        state_next = torch.zeros([number_of_data, n])\n",
    "        # 初始化时间序列信息，此时各个时间序列中各个节点在各个时刻的状态均为0\n",
    "        for num in range(number_of_data):\n",
    "            #使用循环对各个时间序列进行操作\n",
    "            rands = torch.multinomial(torch.Tensor([0.5, 0.5]), n, replacement=True)\n",
    "            state[num, :]=rands\n",
    "            # 随机生成初始状态，状态为0和1的概率均为50%\n",
    "            for i in range(1):\n",
    "            # 从初始状态开始向后演化(T-1)个时刻，该过程会调用sample_one_step函数\n",
    "                sn = self.sample_one_step(rands, edges=self.edges)\n",
    "                state_next[num, :] = sn\n",
    "        return state, state_next, self.edges\n",
    "    def encoding(self, a):\n",
    "        if a.size()[0]==0:return -1\n",
    "        \n",
    "        out = torch.zeros(2**self.number_of_vertice)\n",
    "        v = 0\n",
    "        for i in range(a.size()[0]):\n",
    "            v = v * 2 + int(a[i].item())\n",
    "        out[v] = 1\n",
    "        return out\n",
    "\n",
    "def decimalToBinary(num, lst=[]):\n",
    "    \"\"\"This function converts decimal number\n",
    "    to binary and prints it\"\"\"\n",
    "    if num > 1:\n",
    "        lst = decimalToBinary(num // 2, lst)\n",
    "    lst.append(num % 2)\n",
    "    return lst\n",
    "def discrete_ei(input_sz, n_of_nds, netfunc):\n",
    "    input_x = torch.zeros([input_sz**n_of_nds,n_of_nds])\n",
    "    for i in range(input_sz**n_of_nds):\n",
    "        binary = decimalToBinary(i, [])\n",
    "        ll = len(binary)\n",
    "        for j in range(ll, n_of_nds):\n",
    "            binary = [0] + binary\n",
    "        xx = torch.Tensor(np.array(binary)*1.0).unsqueeze(0)\n",
    "        #print(xx)\n",
    "        input_x[i,:]=xx\n",
    "    pyx_,_,_ = netfunc(input_x)\n",
    "    pyx_ = torch.cat((torch.exp(pyx_.unsqueeze(2)),1-torch.exp(pyx_.unsqueeze(2))), 2)\n",
    "    #print(pyx_.size(),pyx_)\n",
    "    pyx_ = torch.softmax(pyx_, dim=2)\n",
    "    #print(pyx_.size(),pyx_)\n",
    "    pyx = torch.ones([input_sz**n_of_nds,input_sz**n_of_nds])\n",
    "    for i in range(input_x.size()[0]):\n",
    "        for j in range(input_x.size()[0]):\n",
    "            binary = decimalToBinary(j, [])\n",
    "            ll = len(binary)\n",
    "            probability = 1.0\n",
    "            for k in range(n_of_nds):\n",
    "                if k < n_of_nds - ll:\n",
    "                    probability *= pyx_[i,k,0]\n",
    "                else:\n",
    "                    probability *= pyx_[i,k,binary[k-n_of_nds+ll]]\n",
    "            pyx[i, j] = probability\n",
    "\n",
    "    #print(pyx.size(), pyx.sum(1), pyx)\n",
    "    logpyx=torch.log(pyx)\n",
    "    logpyx = torch.where(torch.isinf(logpyx), torch.zeros(logpyx.size()), logpyx)\n",
    "    entropy = pyx * logpyx\n",
    "    sumz = torch.sum(pyx, 0).unsqueeze(0)\n",
    "    logsumz = torch.log(sumz)\n",
    "    logsumz = torch.where(torch.isinf(logsumz), torch.zeros(logsumz.size()), logsumz)\n",
    "    logsumz = logsumz.repeat(input_sz**n_of_nds, 1)\n",
    "    #print(logsumz)\n",
    "\n",
    "    #print(logsumz)\n",
    "    first_term = torch.sum(entropy)\n",
    "    second_term = torch.sum(pyx * logsumz)\n",
    "\n",
    "    final = (first_term - second_term)/(input_sz**n_of_nds) + np.log(input_sz**n_of_nds)\n",
    "    final = final / np.log(2.0)\n",
    "    return final, final*np.log(2.0)/np.log(input_sz**n_of_nds), pyx, first_term, second_term\n",
    "def test_model(batch_size,n_of_nds, net,scale,nll,sim):\n",
    "    batch_size1 = batch_size*10\n",
    "    state,state_next,_ = sim.sample(batch_size)\n",
    "    if use_cuda:\n",
    "        state=state.cuda()\n",
    "        state_next = state_next.cuda()\n",
    "    predict, latent, latent_p = net(state) \n",
    "    predict = nn.functional.log_softmax(predict, dim=1)\n",
    "    loss = nll(predict, state_next)\n",
    "    ssp = net.encoding(state_next)\n",
    "    sigmas = torch.sqrt(torch.mean((ssp-latent_p)**2, 0))\n",
    "    sigmas_matrix = torch.diag(sigmas)\n",
    "    ei = approx_ei(scale, scale, sigmas_matrix.data, lambda x:(net.dynamics(x.unsqueeze(0))+x.unsqueeze(0)), \n",
    "                   num_samples = 1000, L=100, easy=True, device=device)\n",
    "    #one_ei = discrete_ei(input_sz, n_of_nds, lambda x: net(x))\n",
    "    #whole_sigmas = torch.sqrt(torch.mean((state_next-predict)**2, 0))\n",
    "    #inv_sigma = torch.inverse(torch.diag(whole_sigmas))\n",
    "    #whole_ei = approx_ei(scale, scale, inv_sigma.data, lambda x:net(x)[0], 1000, 10)\n",
    "    return ei, sigmas,loss.item()#, one_ei#, whole_ei\n",
    "def calc_ei_loss(latent_p, state_next, net, scale): \n",
    "    ssp = net.encoding(state_next)\n",
    "    prediction = ssp#.detach()\n",
    "    real = latent_p#.detach()\n",
    "    # detach the variables of ssp and latent_p which only optimize the dynamics NN.\n",
    "    sigmas = torch.sqrt(torch.mean((prediction - real)**2, 0))\n",
    "    sigmas_matrix = torch.diag(sigmas)\n",
    "    ei = approx_ei(scale, scale, sigmas_matrix.data, lambda x:(net.dynamics(x.unsqueeze(0))+x.unsqueeze(0)), \n",
    "                   num_samples = 1000, L=100, easy=True, device=device)\n",
    "    return ei, sigmas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def RGB_to_Hex(rgb):\n",
    "  \n",
    "    RGB = list(rgb)\n",
    "    color = '#'\n",
    "    for i in RGB:\n",
    "        num = int(i)\n",
    "        color += str(hex(num))[-2:].replace('x', '0').upper()\n",
    "    return color\n",
    "    \n",
    "def generate_colors(N=12,colormap='hsv'):\n",
    "    \n",
    "    step = max(int(255/N),1)\n",
    "    cmap = plt.get_cmap(colormap)\n",
    "    rgb_list = []\n",
    "    hex_list = []\n",
    "    for i in range(N):\n",
    "        id = step*i # cmap(int)->(r,g,b,a) in 0~1\n",
    "        id = 255 if id>255 else id\n",
    "        rgba_color = cmap(id)\n",
    "        rgb = [int(d*255) for d in rgba_color[:3]]\n",
    "        rgb_list.append(tuple(rgb))\n",
    "        hex_list.append(RGB_to_Hex(rgb))\n",
    "    return rgb_list,hex_list\n",
    "    \n",
    "rgb_list,hex_list = generate_colors(6,'cool')\n",
    "print(rgb_list)\n",
    "print(hex_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Learned Mapping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hidden_units = 64\n",
    "\n",
    "\n",
    "torch.manual_seed(1024)\n",
    "scale = 1\n",
    "batch_size =100\n",
    "nll = nn.functional.binary_cross_entropy_with_logits\n",
    "n_of_nods = 4\n",
    "sim = BnSim(number_of_vertice = n_of_nods)\n",
    "net = Renorm_Dynamic(sym_size = n_of_nods, latent_size = scale, effect_size = n_of_nods, \n",
    "                     hidden_units = hidden_units, normalized_state = False, device=device)\n",
    "net_m = Renorm_Dynamic(sym_size = n_of_nods, latent_size = n_of_nods, effect_size = n_of_nods, \n",
    "                     hidden_units = hidden_units, normalized_state = False, device=device)\n",
    "net = net.cuda() if use_cuda else net\n",
    "net_m = net_m.cuda() if use_cuda else net_m\n",
    "optimizer = torch.optim.Adam([p for p in net.parameters() if p.requires_grad==True], lr=1e-4)\n",
    "optimizer_m = torch.optim.Adam([p for p in net_m.parameters() if p.requires_grad==True], lr=1e-4)\n",
    "\n",
    "random_samples = []\n",
    "\n",
    "for t in range(500001):    \n",
    "    state,state_next,_ = sim.sample(batch_size)\n",
    "    if use_cuda:\n",
    "        state=state.cuda()\n",
    "        state_next = state_next.cuda()\n",
    "    predict, latent, latent_p = net(state)\n",
    "    predict_m, latent_m, latent_pm = net_m(state)\n",
    "    #predict = nn.functional.log_softmax(predict, dim=1)\n",
    "    loss = nll(predict, state_next)\n",
    "    loss_m = nll(predict_m, state_next)\n",
    "    optimizer.zero_grad()\n",
    "    optimizer_m.zero_grad()\n",
    "    loss.backward(retain_graph=True)\n",
    "    loss_m.backward(retain_graph=True)\n",
    "    optimizer.step()\n",
    "    optimizer_m.step()\n",
    "           \n",
    "    if t % 500 ==0:\n",
    "        ei, sigmas,_ = test_model(batch_size, n_of_nods, net, scale, nll,sim)\n",
    "        ei_m, sigmas_m, _ = test_model(batch_size, n_of_nods, net_m, n_of_nods, nll, sim)\n",
    "        print('iter %s:' % t)\n",
    "        print('Macro: loss = %.3f' % loss.item(), ', dEI= %.3f' % ei[0],\n",
    "             ', eff= %.3f' % ei[1], ', std = %.3f' % sigmas.mean().item())\n",
    "        print('Micro: loss = %.3f' % loss_m.item(), ', dEI= %.3f' % ei_m[0],\n",
    "             ', eff= %.3f' % ei_m[1], ', std = %.3f' % sigmas_m.mean().item())\n",
    "        print('Causal Emergence = %.3f' % (ei[0]-ei_m[0]))\n",
    "        \n",
    "        if t % 5000 == 0:\n",
    "            \n",
    "            xx = torch.zeros(state.size()[0],1)\n",
    "            for i in range(state.size()[0]):\n",
    "                v = 0\n",
    "                for j in range(state.size()[1]):\n",
    "                    v = v * 2 + int(state[i,j].item())\n",
    "                xx[i] = v\n",
    "            if latent.size()[1]==1:\n",
    "                yy = latent\n",
    "                if use_cuda:\n",
    "                    xx=xx.cpu()\n",
    "                    yy=yy.cpu()\n",
    "                for i in range(yy.size()[1]):\n",
    "                    plt.plot(xx.data, yy[:,i].data,'.')\n",
    "                plt.xlabel('Learned Latent State')\n",
    "                plt.ylabel('Real Latent State')\n",
    "            else:\n",
    "                lowrank=latent\n",
    "                if latent.size()[1]>2:\n",
    "                    lowrank = torch.pca_lowrank(latent, q=2)[0]\n",
    "                plotx = xx.cpu() if use_cuda else xx\n",
    "                ploty = lowrank.cpu() if use_cuda else lowrank\n",
    "                plt.scatter(ploty.data[:, 0], ploty.data[:, 1], c = plotx.view(-1).data,cmap=plt.cm.Spectral)\n",
    "            plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xx = torch.zeros(state.size()[0],1)\n",
    "for i in range(state.size()[0]):\n",
    "    v = 0\n",
    "    for j in range(state.size()[1]):\n",
    "        v = v * 2 + int(state[i,j].item())\n",
    "    xx[i] = v\n",
    "    \n",
    "yy = latent\n",
    "if use_cuda:\n",
    "    xx=xx.cpu()\n",
    "    yy=yy.cpu()\n",
    "for i in range(yy.size()[1]):\n",
    "    bools = (yy[:,i]<0)\n",
    "    plt.plot(xx[bools].data, yy[bools,i].data,'o',markersize=3, label='Class 1',color='orange')\n",
    "    bools = (yy[:,i]<4)&(yy[:,i]>2)\n",
    "    plt.plot(xx[bools].data, yy[bools,i].data,'s',markersize=4, label='Class 2',color='g')\n",
    "    bools = (yy[:,i]<2)&(yy[:,i]>0)\n",
    "    plt.plot(xx[bools].data, yy[bools,i].data,'<',markersize=4, label='Class 3',color='b')\n",
    "    bools = (yy[:,i]>10)\n",
    "    plt.plot(xx[bools].data, yy[bools,i].data,'*',markersize=5, label='Class 4',color='firebrick')\n",
    "plt.xlabel('Encoded Micro States')\n",
    "plt.ylabel('Latent Macro State')\n",
    "plt.legend()\n",
    "plt.savefig('Boolean_class.svg', dpi=600, format='svg')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Scale Search"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L = 10\n",
    "hidden_units = 64\n",
    "sigma=1\n",
    "#torch.manual_seed(2050)\n",
    "experiments = 4\n",
    "batch_size =100\n",
    "epochs=100001\n",
    "nll = nn.functional.binary_cross_entropy_with_logits\n",
    "\n",
    "err_scale = []\n",
    "ei_scale=[]\n",
    "for scale in [4,3,2,1]:\n",
    "    print('*********************',scale,'**************************')\n",
    "    err_experiment=[]\n",
    "    multi_err_experiment=[]\n",
    "    ei_experiment=[]\n",
    "    for experiment in range(experiments):\n",
    "        n_of_nods = 4\n",
    "        sim = BnSim(number_of_vertice = n_of_nods)\n",
    "        net = Renorm_Dynamic(sym_size = n_of_nods, latent_size = scale, effect_size = n_of_nods, \n",
    "                             hidden_units = hidden_units, normalized_state = False, device=device)\n",
    "        net_m = Renorm_Dynamic(sym_size = n_of_nods, latent_size = n_of_nods, effect_size = n_of_nods, \n",
    "                             hidden_units = hidden_units, normalized_state = False, device=device)\n",
    "        net = net.cuda() if use_cuda else net\n",
    "        net_m = net_m.cuda() if use_cuda else net_m\n",
    "        optimizer = torch.optim.Adam([p for p in net.parameters() if p.requires_grad==True], lr=1e-4)\n",
    "        optimizer_m = torch.optim.Adam([p for p in net_m.parameters() if p.requires_grad==True], lr=1e-4)\n",
    "\n",
    "        random_samples = []\n",
    "\n",
    "        for t in range(epochs):    \n",
    "            state,state_next,_ = sim.sample(batch_size)\n",
    "            if use_cuda:\n",
    "                state=state.cuda()\n",
    "                state_next = state_next.cuda()\n",
    "            predict, latent, latent_p = net(state)\n",
    "            predict_m, latent_m, latent_pm = net_m(state)\n",
    "            #predict = nn.functional.log_softmax(predict, dim=1)\n",
    "            loss = nll(predict, state_next)\n",
    "            loss_m = nll(predict_m, state_next)\n",
    "            optimizer.zero_grad()\n",
    "            optimizer_m.zero_grad()\n",
    "            loss.backward(retain_graph=True)\n",
    "            loss_m.backward(retain_graph=True)\n",
    "            optimizer.step()\n",
    "            optimizer_m.step()\n",
    "            if t % 500 ==0:\n",
    "                ei, sigmas,_ = test_model(batch_size, n_of_nods, net, scale, nll,sim)\n",
    "                ei_m, sigmas_m, _ = test_model(batch_size, n_of_nods, net_m, n_of_nods, nll, sim)\n",
    "                print('iter %s:' % t)\n",
    "                print('Macro: loss = %.3f' % loss.item(), ', dEI= %.3f' % ei[0],\n",
    "                     ', eff= %.3f' % ei[1], ', std = %.3f' % sigmas.mean().item())\n",
    "                print('Micro: loss = %.3f' % loss_m.item(), ', dEI= %.3f' % ei_m[0],\n",
    "                     ', eff= %.3f' % ei_m[1], ', std = %.3f' % sigmas_m.mean().item())\n",
    "                print('Causal Emergence = %.3f' % (ei[0]-ei_m[0]))\n",
    "        #test\n",
    "        ei, sigmas,_ = test_model(batch_size*10, n_of_nods, net, scale, nll, sim) \n",
    "        ei_m, sigmas_m, _ = test_model(batch_size*10, n_of_nods, net_m, n_of_nods, nll, sim)\n",
    "        err_experiment.append(loss)\n",
    "        \n",
    "        #mutual information\n",
    "        ei_experiment.append(ei[0]-ei_m[0])\n",
    "        \n",
    "    err_scale.append(err_experiment)\n",
    "    ei_scale.append(ei_experiment)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgb_list,hex_list = generate_colors(2,'cool')\n",
    "print(rgb_list)\n",
    "print(hex_list)\n",
    "#Bollean data\n",
    "ces=ei_scale\n",
    "\n",
    "scales=torch.Tensor([4,3,2,1])\n",
    "means=[]\n",
    "stds=[]\n",
    "\n",
    "for ce in ces:\n",
    "    m=np.mean(ce)\n",
    "    std=np.std(ce)\n",
    "    means.append(m)\n",
    "    stds.append(std)\n",
    "\n",
    "plt.plot(scales, means,'o')\n",
    "plt.errorbar(scales, means, stds,color=hex_list[0])\n",
    "plt.bar(scales, means, width=0.3, facecolor=hex_list[1], edgecolor='white')\n",
    "plt.title('Relationship between dCE and Scale(q)',fontsize=14) \n",
    "plt.xlabel('Scale ($q$)',fontsize=14)\n",
    "plt.ylabel('Casual Emergence (dCE)',fontsize=14)\n",
    "plt.axhline(0, color='black', linestyle='--')\n",
    "plt.savefig('Boolean_casual_emergence.svg', dpi=600, format='svg')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mutual Information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import entropy_estimators as ee"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def moving_average(interval, windowsize):\n",
    "    window = np.ones(int(windowsize)) / float(windowsize)\n",
    "    re = np.convolve(interval, window, 'same')\n",
    "    return re"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Verification for Theorem 2 when Scale(q)=2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hidden_units =32\n",
    "\n",
    "\n",
    "torch.manual_seed(1024)\n",
    "scale = 2\n",
    "batch_size =100\n",
    "nll = nn.functional.binary_cross_entropy_with_logits\n",
    "n_of_nods = 4\n",
    "sim = BnSim(number_of_vertice = n_of_nods)\n",
    "net = Renorm_Dynamic(sym_size = n_of_nods, latent_size = scale, effect_size = n_of_nods, \n",
    "                     hidden_units = hidden_units, normalized_state = False, device=device)\n",
    "net_m = Renorm_Dynamic(sym_size = n_of_nods, latent_size = n_of_nods, effect_size = n_of_nods, \n",
    "                     hidden_units = hidden_units, normalized_state = False, device=device)\n",
    "net = net.cuda() if use_cuda else net\n",
    "net_m = net_m.cuda() if use_cuda else net_m\n",
    "optimizer = torch.optim.Adam([p for p in net.parameters() if p.requires_grad==True], lr=1e-4)\n",
    "optimizer_m = torch.optim.Adam([p for p in net_m.parameters() if p.requires_grad==True], lr=1e-4)\n",
    "\n",
    "random_samples = []\n",
    "\n",
    "Isst=[]\n",
    "Illt=[]\n",
    "Isshat=[]\n",
    "Irrt=[]\n",
    "I_yt_yt1_ms=[]\n",
    "I_xt_xt1hat_ms=[]\n",
    "errs=[]\n",
    "Ts=[]\n",
    "\n",
    "for t in range(50001):    \n",
    "    state,state_next,_ = sim.sample(batch_size)\n",
    "    if use_cuda:\n",
    "        state=state.cuda()\n",
    "        state_next = state_next.cuda()\n",
    "    predict, latent, latent_p = net(state)\n",
    "    predict_m, latent_m, latent_pm = net_m(state)\n",
    "    #predict = nn.functional.log_softmax(predict, dim=1)\n",
    "    loss = nll(predict, state_next)\n",
    "    loss_m = nll(predict_m, state_next)\n",
    "    optimizer.zero_grad()\n",
    "    optimizer_m.zero_grad()\n",
    "    loss.backward(retain_graph=True)\n",
    "    loss_m.backward(retain_graph=True)\n",
    "    optimizer.step()\n",
    "    optimizer_m.step()\n",
    "           \n",
    "    if t % 500 ==0:\n",
    "       \n",
    "        I_xt_xt1=ee.mi(state.cpu(),state_next.cpu(),k=3)\n",
    "        I_yt_yt1=ee.mi(latent.cpu().detach(),latent_p.cpu().detach(),k=17)\n",
    "        I_xt_xt1hat=ee.mi(state.cpu(),predict.cpu().detach(),k=7)\n",
    "\n",
    "        #I_yt_yt1_m=np.mean([ee.mi(torch.cat([latent_m.cpu().detach()[j] for i in range(6)],0),torch.cat([latent_pm.cpu().detach()[j] for i in range(6)],0)) for j in range(batch_size)])\n",
    "        #I_xt_xt1hat_m=np.mean([ee.mi(torch.cat([state.cpu()[j] for i in range(6)],0),torch.cat([predict_m.cpu().detach()[j] for i in range(6)],0)) for j in range(batch_size)])\n",
    "\n",
    "        Isst.append(I_xt_xt1)\n",
    "        Isshat.append(I_xt_xt1hat)\n",
    "        Irrt.append(I_yt_yt1)\n",
    "\n",
    "        #I_yt_yt1_ms.append(I_yt_yt1_m)\n",
    "        #I_xt_xt1hat_ms.append(I_xt_xt1hat_m)\n",
    "\n",
    "        errs.append(np.std([I_yt_yt1,I_xt_xt1hat,I_xt_xt1]))\n",
    "        Ts.append(t)\n",
    "        \n",
    "        print(I_xt_xt1,I_yt_yt1,I_xt_xt1hat)\n",
    "\n",
    "plt.plot(Ts,Isst,label='I(xt,xt+1)')\n",
    "plt.plot(Ts,Irrt,label='I(yt,yt+1)')\n",
    "plt.plot(Ts,Isshat,label='I(xt,xthat)')\n",
    "#plt.plot(Ts,I_yt_yt1_ms,label='I(xt,xthat)m')\n",
    "#plt.plot(Ts,I_xt_xt1hat_ms,label='I(xt,xthat)m')\n",
    "plt.plot(Ts,errs,label='error')\n",
    "plt.ylim(0,3.5)\n",
    "plt.legend(loc='best')\n",
    "plt.xlabel('Iter')\n",
    "plt.ylabel('Value')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Isst_rol_mean = moving_average(Isst,10)\n",
    "Irrt_rol_mean = moving_average(Irrt,10)\n",
    "Isshat_rol_mean = moving_average(Isshat,10)\n",
    "\n",
    "#plt.plot(Ts,Isst,label='$I(x_t,x_{t+1})$',color='gold',linestyle='-.')\n",
    "#plt.plot(Ts,Irrt,label='$I(y_t,y(t+1))$',color='skyblue',linestyle='-.')\n",
    "#plt.plot(Ts,Isshat,label='$I(x_t,\\hat{x}_{t+1})$',color='thistle',linestyle='-.')\n",
    "#plt.plot(Ts,I_yt_yt1_ms,label='I(xt,xthat)m')\n",
    "#plt.plot(Ts,I_xt_xt1hat_ms,label='I(xt,xthat)m')\n",
    "#plt.plot(Ts,errs,label='error',color=hex_list[3])\n",
    "plt.plot(Ts[4:96],Isst_rol_mean[4:96],label='$I(x_t,x_{t+1})$',color='blueviolet')\n",
    "plt.plot(Ts[4:96],Irrt_rol_mean[4:96],label='$I(y_t,y(t+1))$',color='teal')\n",
    "plt.plot(Ts[4:96],Isshat_rol_mean[4:96],label='$I(x_t,\\hat{x}_{t+1})$',color='b')\n",
    "plt.ylim(0,3)\n",
    "plt.legend(loc='best')\n",
    "plt.title('Verification for Theorem 2 when Scale(q)=2',fontsize=14) \n",
    "plt.xlabel('Iter',fontsize=13)\n",
    "plt.ylabel('Mutual Information Scale(q)=2',fontsize=13)\n",
    "plt.savefig('Boolean_Mutual_Information_Scale2.svg', dpi=600, format='svg')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Verification for Theorem 6 with Scale and Mutual Information "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hidden_units = 32\n",
    "\n",
    "\n",
    "torch.manual_seed(1024)\n",
    "scale = [1,2,3,4]\n",
    "batch_size =100\n",
    "nll = nn.functional.binary_cross_entropy_with_logits\n",
    "n_of_nods = 4\n",
    "sim = BnSim(number_of_vertice = n_of_nods)\n",
    "net_1 = Renorm_Dynamic(sym_size = n_of_nods, latent_size = scale[0], effect_size = n_of_nods, \n",
    "                     hidden_units = hidden_units, normalized_state = False, device=device)\n",
    "net_2 = Renorm_Dynamic(sym_size = n_of_nods, latent_size = scale[1], effect_size = n_of_nods, \n",
    "                     hidden_units = hidden_units, normalized_state = False, device=device)\n",
    "net_3 = Renorm_Dynamic(sym_size = n_of_nods, latent_size = scale[2], effect_size = n_of_nods, \n",
    "                     hidden_units = hidden_units, normalized_state = False, device=device)\n",
    "net_4 = Renorm_Dynamic(sym_size = n_of_nods, latent_size = scale[3], effect_size = n_of_nods, \n",
    "                     hidden_units = hidden_units, normalized_state = False, device=device)\n",
    "net_1 = net_1.cuda() if use_cuda else net_1\n",
    "net_2 = net_2.cuda() if use_cuda else net_2\n",
    "net_3 = net_3.cuda() if use_cuda else net_3\n",
    "net_4 = net_4.cuda() if use_cuda else net_4\n",
    "optimizer_1 = torch.optim.Adam([p for p in net_1.parameters() if p.requires_grad==True], lr=1e-4)\n",
    "optimizer_2 = torch.optim.Adam([p for p in net_2.parameters() if p.requires_grad==True], lr=1e-4)\n",
    "optimizer_3 = torch.optim.Adam([p for p in net_3.parameters() if p.requires_grad==True], lr=1e-4)\n",
    "optimizer_4 = torch.optim.Adam([p for p in net_4.parameters() if p.requires_grad==True], lr=1e-4)\n",
    "\n",
    "random_samples = []\n",
    "\n",
    "Isshat_1=[]\n",
    "Isshat_2=[]\n",
    "Isshat_3=[]\n",
    "Isshat_4=[]\n",
    "#Isshat_m=[]\n",
    "I_xt_yt_1s=[]\n",
    "I_xt_yt_2s=[]\n",
    "I_xt_yt_3s=[]\n",
    "I_xt_yt_4s=[]\n",
    "#errs=[]\n",
    "Ts=[]\n",
    "\n",
    "for t in range(50001):    \n",
    "    state,state_next,_ = sim.sample(batch_size)\n",
    "    if use_cuda:\n",
    "        state=state.cuda()\n",
    "        state_next = state_next.cuda()\n",
    "    predict_1, latent_1, latent_p1 = net_1(state)\n",
    "    predict_2, latent_2, latent_p2 = net_2(state)\n",
    "    predict_3, latent_3, latent_p3 = net_3(state)\n",
    "    predict_4, latent_4, latent_p4 = net_4(state)\n",
    "    #predict = nn.functional.log_softmax(predict, dim=1)\n",
    "    loss_1 = nll(predict_1, state_next)\n",
    "    loss_2 = nll(predict_2, state_next)\n",
    "    loss_3 = nll(predict_3, state_next)\n",
    "    loss_4 = nll(predict_4, state_next)\n",
    "    optimizer_1.zero_grad()\n",
    "    optimizer_2.zero_grad()\n",
    "    optimizer_3.zero_grad()\n",
    "    optimizer_4.zero_grad()\n",
    "    loss_1.backward(retain_graph=True)\n",
    "    loss_2.backward(retain_graph=True)\n",
    "    loss_3.backward(retain_graph=True)\n",
    "    loss_4.backward(retain_graph=True)\n",
    "    optimizer_1.step()\n",
    "    optimizer_2.step()\n",
    "    optimizer_3.step()\n",
    "    optimizer_4.step()\n",
    "           \n",
    "    if t % 500 ==0:\n",
    "       \n",
    "        I_xt_xt1hat_1=ee.mi(state.cpu(),predict_1.cpu().detach(),k=9)\n",
    "        I_xt_yt_1=ee.mi(state.cpu(),latent_1.cpu().detach(),k=9)\n",
    "\n",
    "        I_xt_yt_2=ee.mi(state.cpu(),latent_2.cpu().detach(),k=9)\n",
    "        I_xt_xt1hat_2=ee.mi(state.cpu(),predict_2.cpu().detach(),k=9)\n",
    "\n",
    "        I_xt_xt1hat_3=ee.mi(state.cpu(),predict_3.cpu().detach(),k=9)\n",
    "        I_xt_yt_3=ee.mi(state.cpu(),latent_3.cpu().detach(),k=9)\n",
    "\n",
    "        I_xt_xt1hat_4=ee.mi(state.cpu(),predict_4.cpu().detach(),k=9)\n",
    "        I_xt_yt_4=ee.mi(state.cpu(),latent_4.cpu().detach(),k=9)\n",
    "\n",
    "        Isshat_1.append(I_xt_xt1hat_1)\n",
    "        Isshat_2.append(I_xt_xt1hat_2)\n",
    "        Isshat_3.append(I_xt_xt1hat_3)\n",
    "        Isshat_4.append(I_xt_xt1hat_4)\n",
    "        I_xt_yt_1s.append(I_xt_yt_1)\n",
    "        I_xt_yt_2s.append(I_xt_yt_2)\n",
    "        I_xt_yt_3s.append(I_xt_yt_3)\n",
    "        I_xt_yt_4s.append(I_xt_yt_4)\n",
    "\n",
    "        #errs.append(np.std([np.mean(I_yt_yt1),np.mean(I_xt_xt1hat),np.mean(I_xt_xt1)]))\n",
    "        Ts.append(t)\n",
    "        \n",
    "        print(I_xt_xt1hat_4,I_xt_yt_1,I_xt_yt_2,I_xt_yt_3,I_xt_yt_4)\n",
    "\n",
    "plt.plot(Ts,Isshat_4,label='I(xt,xt_hat)')\n",
    "#plt.plot(Ts,Isshat_m,label='I(xt,xt_hat_m)')\n",
    "plt.plot(Ts,I_xt_yt_1s,label='I(xt,yt_1)')\n",
    "plt.plot(Ts,I_xt_yt_2s,label='I(xt,yt_2)')\n",
    "plt.plot(Ts,I_xt_yt_3s,label='I(xt,yt_3)')\n",
    "plt.plot(Ts,I_xt_yt_4s,label='I(xt,yt_4)')\n",
    "#plt.plot(Ts,I_xt_xt1hat_ms,label='I(xt,xthat)m')\n",
    "#plt.plot(Ts,errs,label='error')\n",
    "plt.legend(loc='best')\n",
    "plt.xlabel('Iter')\n",
    "plt.ylabel('Value')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Isshat_4_rol_mean = moving_average(Isshat_4,10)\n",
    "Isshat_3_rol_mean = moving_average(Isshat_3,10)\n",
    "Isshat_2_rol_mean = moving_average(Isshat_2,10)\n",
    "Isshat_1_rol_mean = moving_average(Isshat_1,10)\n",
    "I_xt_yt_1s_rol_mean = moving_average(I_xt_yt_1s,10)\n",
    "I_xt_yt_2s_rol_mean = moving_average(I_xt_yt_2s,10)\n",
    "I_xt_yt_3s_rol_mean = moving_average(I_xt_yt_3s,10)\n",
    "I_xt_yt_4s_rol_mean = moving_average(I_xt_yt_4s,10)\n",
    "\n",
    "plt.plot(Ts,Isshat_3,label='I(xt,xt_hat)',color='deepskyblue')\n",
    "#plt.plot(Ts,Isshat_m,label='I(xt,xt_hat_m)')\n",
    "plt.plot(Ts,I_xt_yt_1s,label='I(xt,yt_1)',color='gold')\n",
    "plt.plot(Ts,I_xt_yt_2s,label='I(xt,yt_2)',color='lime')\n",
    "plt.plot(Ts,I_xt_yt_3s,label='I(xt,yt_3)',color='lightcoral')\n",
    "plt.plot(Ts,I_xt_yt_4s,label='I(xt,yt_4)',color='violet')\n",
    "plt.plot(Ts[6:95],Isshat_3_rol_mean[6:95],label='I(xt,xt_hat) Rolling Mean',color='b')\n",
    "plt.plot(Ts[6:95],I_xt_yt_1s_rol_mean[6:95],label='I(xt,yt_1) Rolling Mean',color='darkorange')\n",
    "plt.plot(Ts[6:95],I_xt_yt_2s_rol_mean[6:95],label='I(xt,yt_2) Rolling Mean',color='g')\n",
    "plt.plot(Ts[6:95],I_xt_yt_3s_rol_mean[6:95],label='I(xt,yt_3) Rolling Mean',color='r')\n",
    "plt.plot(Ts[6:95],I_xt_yt_4s_rol_mean[6:95],label='I(xt,yt_4) Rolling Mean',color='purple')\n",
    "#plt.plot(Ts,I_xt_xt1hat_ms,label='I(xt,xthat)m')\n",
    "#plt.plot(Ts,errs,label='error')\n",
    "plt.legend(loc='best',fontsize=7)\n",
    "plt.ylim(0,3.5)\n",
    "plt.xlabel('Iter',fontsize=14)\n",
    "plt.ylabel('Mutual Information',fontsize=14)\n",
    "#plt.savefig('Boolean Mutual Information rm.svg', dpi=600, format='svg')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgb_list,hex_list = generate_colors(12,'cubehelix')\n",
    "print(rgb_list)\n",
    "print(hex_list)\n",
    "#plt.plot(Ts,Isshat_3,label='$I(x_t,\\hat{x}_t)$',color='skyblue',linestyle='-.')\n",
    "#plt.plot(Ts,Isshat_m,label='I(xt,xt_hat_m)')\n",
    "#plt.plot(Ts,I_xt_yt_1s,label='$I(x_t,y_{t}^{q=1})$',color='orange',linestyle=':')\n",
    "#plt.plot(Ts,I_xt_yt_2s,label='$I(x_t,y_{t}^{q=2})$',color='limegreen',linestyle=':')\n",
    "#plt.plot(Ts,I_xt_yt_3s,label='$I(x_t,y_{t}^{q=3})$',color='violet',linestyle=':')\n",
    "#plt.plot(Ts,I_xt_yt_4s,label='$I(x_t,y_{t}^{q=4})$',color='deepskyblue',linestyle=':')\n",
    "plt.plot(Ts[6:95],Isshat_3_rol_mean[6:95],label='$I(x_t,\\hat{x}_{t+1})$',color='b')\n",
    "plt.plot(Ts[6:95],I_xt_yt_1s_rol_mean[6:95],label='$I(x_t,y_{t}^{q=1})$',color='darkorange')\n",
    "plt.plot(Ts[6:95],I_xt_yt_2s_rol_mean[6:95],label='$I(x_t,y_{t}^{q=2})$',color='g')\n",
    "plt.plot(Ts[6:95],I_xt_yt_3s_rol_mean[6:95],label='$I(x_t,y_{t}^{q=3})$',color='r')\n",
    "plt.plot(Ts[6:95],I_xt_yt_4s_rol_mean[6:95],label='$I(x_t,y_{t}^{q=4})$',color='purple')\n",
    "#plt.plot(Ts,I_xt_xt1hat_ms,label='I(xt,xthat)m')\n",
    "#plt.plot(Ts,errs,label='error')\n",
    "plt.legend(loc='best')\n",
    "plt.title('Verification for Theorem 6',fontsize=14) \n",
    "plt.ylim(0,3)\n",
    "plt.xlabel('Iter',fontsize=14)\n",
    "plt.ylabel('Mutual Information',fontsize=14)\n",
    "plt.savefig('Boolean_Mutual_Information_rm.svg', dpi=600, format='svg')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.7.0 ('base')",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "88279d2366fe020547cde40dd65aa0e3aa662a6ec1f3ca12d88834876c85e1a6"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}