TREX_wrapper.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## T-REX wrapper\n",
    "\n",
    "First we import functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from tqdm import tqdm\n",
    "from newick import read\n",
    "from datasets.loading import load_data\n",
    "import os\n",
    "import ipdb\n",
    "\n",
    "os.environ['DATAPATH']=\"./data\"\n",
    "os.environ['SAVEPATH']=\"./embeddings\"\n",
    "\n",
    "def sim2metric(S):\n",
    "    \"\"\"\n",
    "    Input S is a n-by-n np array, representing similarities.\n",
    "    Output D turn S into a metric.\n",
    "    It should satisfy:\n",
    "    1) Sij = 1 iff Dij = 0\n",
    "    2) Sij = -1 iff Dij = max(D)\n",
    "    3) basic metric properties (triangle ineq...etc)\n",
    "    \n",
    "    Now we choose D = 1-S\n",
    "    \"\"\"\n",
    "    return 1-S\n",
    "\n",
    "def save_dist2phylip(D,path):\n",
    "    n = D.shape[0]\n",
    "    f = open(path, \"w\")\n",
    "    f.write(str(n)+'\\n')\n",
    "    for i in tqdm(range(n)):\n",
    "        f.write(str(i)+' ')\n",
    "        np.savetxt(f,D[i],fmt='%.4f',newline=\" \")\n",
    "        f.write('\\n')\n",
    "    f.close()\n",
    "    print(\"Matrix saved\")\n",
    "\n",
    "def newick2dist(trees,n):\n",
    "    # This convert a tree from newick format to distance matrix (n-by-n)\n",
    "    # Note that all internal nodes are named \"None\"\n",
    "    # All leaves node have name in np.arange(n)\n",
    "    \n",
    "    D = np.zeros((n,n))\n",
    "    D = Fill_D_from_tree(trees[0],D)\n",
    "    return D\n",
    "\n",
    "def Fill_D_from_tree(subtree,D):\n",
    "    All_leaves = np.arange(D.shape[0])    \n",
    "\n",
    "    num_children = len(subtree.descendants)\n",
    "    for i in range(num_children):\n",
    "        sub_leaves = np.array(subtree.descendants[i].get_leaf_names()).astype(int)\n",
    "        # We simply set the negative branch weights to zero following the common practice\n",
    "        if subtree.descendants[i].length>0:\n",
    "            D[np.ix_(sub_leaves,np.setdiff1d(All_leaves,sub_leaves))] += subtree.descendants[i].length\n",
    "            D[np.ix_(np.setdiff1d(All_leaves,sub_leaves),sub_leaves)] += subtree.descendants[i].length\n",
    "    \n",
    "    \n",
    "        if subtree.descendants[i].name == None:\n",
    "            D = Fill_D_from_tree(subtree.descendants[i],D)\n",
    "    \n",
    "    return D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "from scipy.sparse import csr_matrix, lil_matrix\n",
    "from scipy.sparse.csgraph import shortest_path\n",
    "import newick\n",
    "\n",
    "def newick2W(trees,n):\n",
    "    # This convert a tree from newick format to distance csr matrix, including interior nodes ((2n-1)-by-(2n-1))\n",
    "    # Note that all internal nodes are named \"None\"\n",
    "    # All leaves node have name in np.arange(n)\n",
    "    \n",
    "    W = csr_matrix(((2*n-1),(2*n-1)))\n",
    "    assert len(trees[0].descendants) == 2 # An internal/root node should always have 2 children!\n",
    "    \n",
    "    W,_,_ = Fill_W_from_tree(trees[0],W,n)\n",
    "    return W\n",
    "\n",
    "def Fill_W_from_tree(subtree,W,cur_idx):   \n",
    "    # Check if we visit current node\n",
    "    if subtree.name == None:\n",
    "        # Not visit yet, give id\n",
    "        subtree.name = str(cur_idx)\n",
    "        cur_idx += 1\n",
    "    # Now add all adjacent relation to W\n",
    "    if not subtree.is_leaf:\n",
    "        Flag = False\n",
    "        for i in range(len(subtree.descendants)):\n",
    "            if subtree.descendants[i].name == None:\n",
    "                subtree.descendants[i].name = str(cur_idx)\n",
    "                cur_idx += 1\n",
    "                Flag = True\n",
    "\n",
    "            W[int(subtree.descendants[i].name),int(subtree.name)] = subtree.descendants[i].length\n",
    "            W[int(subtree.name),int(subtree.descendants[i].name)] = subtree.descendants[i].length\n",
    "\n",
    "            if Flag:\n",
    "                W,subtree.descendants[i],cur_idx = Fill_W_from_tree(subtree.descendants[i],W,cur_idx)\n",
    "    \n",
    "    return W,subtree,cur_idx\n",
    "\n",
    "def generate_rand_TreeMetric_V2(N,option='unweight',noise_ratio=0.0):\n",
    "    # Generate (weighted) random binary tree in newick format.\n",
    "    T = generate_tree_newick(np.arange(N),option=option)\n",
    "    # Make it into a tree data structure in python\n",
    "    T = newick.loads(T)\n",
    "    # Get the distance matrix for leaves in T\n",
    "    W = newick2W(T,N)\n",
    "    \n",
    "    # Add noise. We randomly add noise_ratio*(2n-2) fake connections between nodes (including internal nodes).\n",
    "    m_Noise = int(noise_ratio*(2*N-2))\n",
    "    for _ in range(m_Noise):\n",
    "        idx0,idx1 = np.random.choice(N,size=(2),replace=False)\n",
    "        if option == 'unweight':\n",
    "            W[idx0,idx1] = 1.0\n",
    "            W[idx1,idx0] = 1.0\n",
    "        else:\n",
    "            W[idx0,idx1] = np.random.random()\n",
    "            W[idx1,idx0] = np.random.random()\n",
    "    \n",
    "    # Now construct the shortest path based on this noisy observation\n",
    "    Noisy_D = get_leaves_Dmatrix(W,N)\n",
    "    \n",
    "    return Noisy_D\n",
    "\n",
    "def generate_rand_TreeMetric(N,option='unweight',noise_scale=1.0,p=2):\n",
    "    # Generate (weighted) random binary tree in newick format.\n",
    "    T = generate_tree_newick(np.arange(N),option='unweight')\n",
    "    # Make it into a tree data structure in python\n",
    "    T = newick.loads(T)\n",
    "    # Get the distance matrix for leaves in T\n",
    "    D = newick2dist(T,N)\n",
    "    \n",
    "    # Add noise (noise_scale*uniform([0,1])) to each leaf distance. \n",
    "    Noise = np.around(np.triu(noise_scale*np.random.random(size=D.shape),1),5)\n",
    "    Noise = Noise + Noise.transpose()\n",
    "    Noisy_D = D + Noise\n",
    "    Noise_norm = (np.abs(Noise)**p).sum().sum()**(1/p)\n",
    "    print('Noise_norm:{0:4f}'.format(Noise_norm))\n",
    "    \n",
    "    return Noisy_D\n",
    "\n",
    "def generate_tree_newick(L,option='unweight'):\n",
    "    # base case\n",
    "    if len(L) == 1: \n",
    "        return str(L[0])\n",
    "    else:\n",
    "        perm_L = np.random.permutation(L)\n",
    "        split = np.random.randint(1,len(L))\n",
    "        left = perm_L[:split]\n",
    "        right = perm_L[split:]\n",
    "    # recursion\n",
    "    if option == 'unweight':\n",
    "        return str((generate_tree_newick(left)+':1.0', generate_tree_newick(right)+':1.0')).replace(' ','').replace(\"'\",\"\")\n",
    "    else:\n",
    "        left_weight = ':{0:.5f}'.format(np.random.rand())\n",
    "        right_weight = ':{0:.5f}'.format(np.random.rand())\n",
    "        return str((generate_tree_newick(left,option=option)+left_weight, generate_tree_newick(right,option=option)+right_weight)).replace(' ','').replace(\"'\",\"\")\n",
    "\n",
    "def newick2dist(trees,n):\n",
    "    # This convert a tree from newick format to distance matrix (n-by-n)\n",
    "    # Note that all internal nodes are named \"None\"\n",
    "    # All leaves node have name in np.arange(n)\n",
    "    \n",
    "    D = np.zeros((n,n))\n",
    "    assert len(trees[0].descendants) == 2 # An internal/root node should always have 2 children!\n",
    "    \n",
    "    \n",
    "    D = Fill_D_from_tree(trees[0],D)\n",
    "    return D\n",
    "\n",
    "def Fill_D_from_tree(subtree,D):\n",
    "    All_leaves = np.arange(D.shape[0])    \n",
    "\n",
    "    num_children = len(subtree.descendants)\n",
    "    for i in range(num_children):\n",
    "        sub_leaves = np.array(subtree.descendants[i].get_leaf_names()).astype(int)\n",
    "        # We simply set the negative branch weights to zero following the common practice\n",
    "        if subtree.descendants[i].length>0:\n",
    "            D[np.ix_(sub_leaves,np.setdiff1d(All_leaves,sub_leaves))] += subtree.descendants[i].length\n",
    "            D[np.ix_(np.setdiff1d(All_leaves,sub_leaves),sub_leaves)] += subtree.descendants[i].length\n",
    "    \n",
    "    \n",
    "        if subtree.descendants[i].name == None:\n",
    "            D = Fill_D_from_tree(subtree.descendants[i],D)\n",
    "    \n",
    "    return D\n",
    "\n",
    "def generate_rand_Emetric(num_nodes,E_rank=2):\n",
    "    X = np.random.normal(size=(num_nodes,E_rank))\n",
    "    X = X/np.linalg.norm(X,axis=1,keepdims=True)*np.random.rand(num_nodes,1)\n",
    "    D_metric = np.zeros((num_nodes,num_nodes))\n",
    "    all_pairs = generate_all_pairs(num_nodes)\n",
    "    D_metric[all_pairs[:,0],all_pairs[:,1]] = np.linalg.norm(X[all_pairs[:,0],:]-X[all_pairs[:,1],:],axis=1)\n",
    "    D_metric[all_pairs[:,1],all_pairs[:,0]] = D_metric[all_pairs[:,0],all_pairs[:,1]]\n",
    "    return D_metric\n",
    "\n",
    "def generate_Emetric(x):\n",
    "    D_metric = np.zeros((x.shape[0],x.shape[0]))\n",
    "    all_pairs = generate_all_pairs(x.shape[0])\n",
    "    D_metric[all_pairs[:,0],all_pairs[:,1]] = np.linalg.norm(x[all_pairs[:,0],:]-x[all_pairs[:,1],:],axis=1)\n",
    "    D_metric[all_pairs[:,1],all_pairs[:,0]] = D_metric[all_pairs[:,0],all_pairs[:,1]]\n",
    "    # Round first to prevent numerical issue...\n",
    "    D_metric = np.round(D_metric,5)\n",
    "    return D_metric\n",
    "\n",
    "def sim2metric(S):\n",
    "    \"\"\"\n",
    "    Input S is a n-by-n np array, representing similarities.\n",
    "    Output D turn S into a metric.\n",
    "    It should satisfy:\n",
    "    1) Sij = 1 iff Dij = 0\n",
    "    2) Sij = -1 iff Dij = max(D)\n",
    "    3) basic metric properties (triangle ineq...etc)\n",
    "    \n",
    "    Now we choose D = 1-S\n",
    "    \"\"\"\n",
    "    return 1-S\n",
    "\n",
    "def generate_tree(L):\n",
    "    # base case\n",
    "    if len(L) == 1: \n",
    "        return L[0]\n",
    "    else:\n",
    "        perm_L = np.random.permutation(L)\n",
    "        split = np.random.randint(1,len(L))\n",
    "        left = perm_L[:split]\n",
    "        right = perm_L[split:]\n",
    "    # recursion\n",
    "    return (generate_tree(left), generate_tree(right))\n",
    "\n",
    "def plot_rand_binary_tree(brt,n,plot=True,normalization=False,d_mtx=None):\n",
    "    w_star = None\n",
    "    A = np.zeros((2*n-1,2*n-1))\n",
    "    # Note that we keep the leave node id the same!!!\n",
    "    # Thus the internal node id start at n.\n",
    "    # For clarity, we set node id n to be the root node\n",
    "    A,next_idx = Fill_adj(brt,n,A)\n",
    "    \n",
    "    assert next_idx == 2*n-1\n",
    "    \n",
    "    if normalization:\n",
    "        # Apply Puoya's normalized weight\n",
    "        A_T=np.zeros((n*(n-1)/2,2*n-2))\n",
    "        # First construct an edge name dictionary\n",
    "        w_dict = -np.ones((2*n-1,2*n-1))\n",
    "        w_dict,next_idx,cur_w_idx = Fill_w_dict(brt,n,w_dict)\n",
    "        assert next_idx == 2*n-1\n",
    "        assert cur_w_idx == 2*n-2\n",
    "        \n",
    "        # Now go through all shortest path for leaves and get d_vec\n",
    "        G = nx.from_numpy_matrix(np.matrix(A))\n",
    "        All_path = nx.shortest_path(G)\n",
    "        count = 0\n",
    "        d_vec = np.zeros(n*(n-1)/2)\n",
    "        for i in range(n):\n",
    "            for j in range(i+1,n):\n",
    "                d_vec[count] = d_mtx[i,j]\n",
    "                path = All_path[i][j]\n",
    "                for k in range(len(path)-1):\n",
    "                    A_T[count,int(w_dict[path[k],path[k+1]])] = 1\n",
    "                count += 1\n",
    "        \n",
    "        # optimal weight\n",
    "        w_star = np.matmul(np.matmul(np.linalg.pinv(np.matmul(A_T.transpose(),A_T)),A_T.transpose()),d_vec)\n",
    "        # modify adj matrix A\n",
    "        for i in range(2*n-1):\n",
    "            for j in range(i+1,2*n-1):\n",
    "                if A[i,j] == 1:\n",
    "                    A[i,j] = w_star[int(w_dict[i,j])]\n",
    "                    A[j,i] = w_star[int(w_dict[i,j])]\n",
    "    \n",
    "    # Compute d, the shortest path for leave nodes\n",
    "    G = nx.from_numpy_matrix(np.matrix(A))\n",
    "    d = nx.algorithms.shortest_paths.dense.floyd_warshall_numpy(G)\n",
    "    d = d[:n]\n",
    "    d = d[:,:n]\n",
    "    \n",
    "    if plot:\n",
    "        plot_Tree_metric(A,n,n)\n",
    "    \n",
    "    return d,w_star\n",
    "\n",
    "def Fill_adj(brt,cur_idx,A):\n",
    "    left = brt[0]\n",
    "    if isinstance(left,np.int64):\n",
    "        A[cur_idx,left] = 1\n",
    "        A[left,cur_idx] = 1\n",
    "        next_idx = cur_idx+1\n",
    "    else:\n",
    "        cur_idx_left = cur_idx+1\n",
    "        A[cur_idx,cur_idx_left] = 1\n",
    "        A[cur_idx_left,cur_idx] = 1\n",
    "        A,next_idx = Fill_adj(left,cur_idx_left,A)\n",
    "    \n",
    "    right = brt[1]\n",
    "    if isinstance(right,np.int64):\n",
    "        A[cur_idx,right] = 1\n",
    "        A[right,cur_idx] = 1\n",
    "    else:\n",
    "        cur_idx_right = next_idx\n",
    "        A[cur_idx,cur_idx_right] = 1\n",
    "        A[cur_idx_right,cur_idx] = 1\n",
    "        A,next_idx = Fill_adj(right,cur_idx_right,A)\n",
    "\n",
    "    return A,next_idx\n",
    "    \n",
    "def Fill_w_dict(brt,cur_idx,w_dict,cur_w_idx=0):\n",
    "    left = brt[0]\n",
    "    if isinstance(left,np.int64):\n",
    "        w_dict[cur_idx,left] = cur_w_idx\n",
    "        w_dict[left,cur_idx] = cur_w_idx\n",
    "        next_idx = cur_idx+1\n",
    "        cur_w_idx += 1\n",
    "    else:\n",
    "        cur_idx_left = cur_idx+1\n",
    "        w_dict[cur_idx,cur_idx_left] = cur_w_idx\n",
    "        w_dict[cur_idx_left,cur_idx] = cur_w_idx\n",
    "        cur_w_idx += 1\n",
    "        w_dict,next_idx,cur_w_idx = Fill_w_dict(left,cur_idx_left,w_dict,cur_w_idx)\n",
    "    \n",
    "    right = brt[1]\n",
    "    if isinstance(right,np.int64):\n",
    "        w_dict[cur_idx,right] = cur_w_idx\n",
    "        w_dict[right,cur_idx] = cur_w_idx\n",
    "        cur_w_idx += 1\n",
    "    else:\n",
    "        cur_idx_right = next_idx\n",
    "        w_dict[cur_idx,cur_idx_right] = cur_w_idx\n",
    "        w_dict[cur_idx_right,cur_idx] = cur_w_idx\n",
    "        cur_w_idx += 1\n",
    "        w_dict,next_idx,cur_w_idx = Fill_w_dict(right,cur_idx_right,w_dict,cur_w_idx)\n",
    "\n",
    "    return w_dict,next_idx,cur_w_idx\n",
    "\n",
    "def plot_Tree_metric(W,N_leaves,root_id=-1,Return_metric=False,plot=True,option='scipy'):\n",
    "    \n",
    "    print('Start computing shortest path')\n",
    "    \n",
    "    # First to remove potential zero rows\n",
    "    W = W[~np.all(W == 0, axis=1)]\n",
    "    W = W[:,~np.all(W == 0, axis=0)]\n",
    "    \n",
    "#     # Replace all negative entries to 0 (if any)\n",
    "#     W[W<0] = 0.0\n",
    "    \n",
    "#     # Rounding for better vitualization\n",
    "#     W = np.around(W,3)\n",
    "    \n",
    "    G = nx.from_numpy_matrix(np.matrix(W))\n",
    "    \n",
    "    if plot:\n",
    "        # Use spring_layout to handle positioning of graph\n",
    "        layout = nx.kamada_kawai_layout(G)\n",
    "\n",
    "        # # Use a list for node_sizes\n",
    "        # sizes = [1000,400,200]\n",
    "\n",
    "        # # Use a list for node colours\n",
    "        # # Here we use different color to distiguish leave nodes and steiner nodes\n",
    "        # # By default, we use 'tab:blue' for leave nodes and 'tab:red' for steiner nodes\n",
    "        N_all = W.shape[0]\n",
    "        color_map = []\n",
    "        for _ in range(N_leaves):\n",
    "            color_map.append('tab:blue')\n",
    "        for _ in range(N_all-N_leaves):\n",
    "            color_map.append('tab:red')\n",
    "\n",
    "        if root_id>-1:\n",
    "            color_map[root_id] = 'tab:green'\n",
    "\n",
    "        # Draw the graph using the layout - with_labels=True if you want node labels.\n",
    "        nx.draw(G, layout, with_labels=True, node_color= color_map)\n",
    "\n",
    "        # Get weights of each edge and assign to labels\n",
    "        labels = nx.get_edge_attributes(G, \"weight\")\n",
    "\n",
    "        # Draw edge labels using layout and list of labels\n",
    "        nx.draw_networkx_edge_labels(G, pos=layout, edge_labels=labels)\n",
    "    \n",
    "    # Also compute the shortest distance matrix for leaves\n",
    "    if Return_metric:\n",
    "        if option == 'networkx':\n",
    "            n = N_leaves\n",
    "            d = nx.algorithms.shortest_paths.dense.floyd_warshall_numpy(G)\n",
    "            d = d[:n]\n",
    "            d = d[:,:n]\n",
    "            return np.array(d)\n",
    "        elif option == 'scipy':\n",
    "            n = N_leaves\n",
    "            graph = csr_matrix(W) \n",
    "            d = shortest_path(csgraph=graph,directed=False)\n",
    "            d = d[:n]\n",
    "            d = d[:,:n]\n",
    "            return d\n",
    "    else:\n",
    "        return None\n",
    "    \n",
    "def is_ultrametric(M):\n",
    "    # First check if all entries are non-negative!\n",
    "    if M[M<0].any():\n",
    "        print(\"Input metric has negative distance!\")\n",
    "        return False\n",
    "    n = M.shape[0]\n",
    "    \n",
    "    for i in range(n):\n",
    "        for j in range(i):\n",
    "            for k in range(j):\n",
    "                if M[i,j] > max(M[i,k],M[j,k]):\n",
    "                    return False\n",
    "                if M[i,k] > max(M[i,j],M[j,k]):\n",
    "                    return False \n",
    "                if M[j,k] > max(M[i,j],M[i,k]):\n",
    "                    return False\n",
    "    \n",
    "    return True\n",
    "\n",
    "def is_treemetric(M):\n",
    "    # First check if all entries are non-negative!\n",
    "    if M[M<0].any():\n",
    "        print(\"Input metric has negative distance!\")\n",
    "        return False\n",
    "    \n",
    "    # Also compute Gromov hyperbolicity (delta)\n",
    "    n = M.shape[0]\n",
    "    Delta = 0\n",
    "    Flag = True\n",
    "    for i in tqdm(range(n)):\n",
    "        for j in range(i):\n",
    "            for k in range(j):\n",
    "                for m in range(k):\n",
    "                    ARRAY = np.array([i,j,k,m])\n",
    "                    for elmt in list(itertools.permutations(ARRAY)):\n",
    "                        if M[elmt[0],elmt[1]]+M[elmt[2],elmt[3]] > max(M[elmt[0],elmt[2]]+M[elmt[1],elmt[3]],M[elmt[0],elmt[3]]+M[elmt[1],elmt[2]]):\n",
    "                            diff = M[elmt[0],elmt[1]]+M[elmt[2],elmt[3]] - max(M[elmt[0],elmt[2]]+M[elmt[1],elmt[3]],M[elmt[0],elmt[3]]+M[elmt[1],elmt[2]])\n",
    "                            if diff > Delta:\n",
    "                                Delta = diff\n",
    "                            Flag = False\n",
    "                            \n",
    "                    \n",
    "    return Flag, Delta\n",
    "\n",
    "def get_leaves_Dmatrix(W,N_leaves):\n",
    "    # Input W is a csr_matrix.\n",
    "    n = N_leaves \n",
    "    d = shortest_path(csgraph=W,directed=False)\n",
    "    d = d[:n]\n",
    "    d = d[:,:n]\n",
    "    return d\n",
    "\n",
    "def Linkage2dist(Z):\n",
    "    # First, build a dictionary for mapping all nodes to actual cluters\n",
    "    n = Z.shape[0]+1\n",
    "    # Include all leaves as singleton first\n",
    "    W_Size = int(max(Z[:,0].max(),Z[:,1].max()))+2\n",
    "    W = csr_matrix((W_Size,W_Size))    \n",
    "\n",
    "    for i in range(Z.shape[0]):\n",
    "        idx0 = Z[i,0]\n",
    "        idx1 = Z[i,1]\n",
    "        length = Z[i,2]\n",
    "        pid = i+n\n",
    "        W[idx0,pid] = length\n",
    "        W[idx1,pid] = length\n",
    "        W[pid,idx0] = length\n",
    "        W[pid,idx1] = length\n",
    "\n",
    "    d = shortest_path(csgraph=W,directed=False)\n",
    "    d = d[:n]\n",
    "    d = d[:,:n]\n",
    "    return d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Save the dissimilarity matrix for real world datasets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DATASET = 'zoo'\n",
    "CURV = '100' # set it to be the curvature you use for using HyperAid.\n",
    "SEED = '0' # Change if necessary.\n",
    "x, y_true, similarities = load_data(DATASET)\n",
    "D_metric = sim2metric(similarities)\n",
    "np.save('./embeddings/{}/seed_{}_curv_{}/D_target_{}_{}.npy'.format(DATASET,SEED,CURV,DATASET,SEED),D_metric) #"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Save the raw input metric and trained hyperbolic metric in phylip format. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "Note that D_metric_$dataset_$seed.npy is D_hyp\n",
    "D_hyp_$dataset_$seed.phylip is its phylip format\n",
    "D_metric_$dataset_$seed.phylip is the taget metric in phylip format\n",
    "\"\"\"\n",
    "\n",
    "x, y_true, similarities = load_data(DATASET)\n",
    "D_metric = sim2metric(similarities)\n",
    "\n",
    "D_hyp_path = './embeddings/{}/seed_{}_curv_{}/D_metric_segmentation_{}.npy'.format(DATASET,SEED,CURV,DATASET,SEED)\n",
    "D_hyp = np.load(D_hyp_path)\n",
    "\n",
    "D_metric_phylo_path = './embeddings/{}/seed_{}_curv_{}/D_metric.phylip'.format(DATASET,SEED,CURV)\n",
    "D_hyp_phylo_path = './embeddings/{}/seed_{}_curv_{}/D_hyp.phylip'.format(DATASET,SEED,CURV)\n",
    "save_dist2phylip(D_metric,D_metric_phylo_path)\n",
    "save_dist2phylip(D_hyp,D_hyp_phylo_path)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Applying T-REX [website](http://www.trex.uqam.ca/index.php?action=home)\n",
    "\n",
    "1) First click \"NJ\" on the left panel.\n",
    "2) Choose \"Distance matrix\".\n",
    "3) Click \"File\" then upload the .phylip file.\n",
    "4) Choose the desired tree reconstruction method, NJ or UNJ.\n",
    "5) After redirecting to the result page, click \"Fitting statistics\".\n",
    "6) Copy the numbers right after \"TREE METRIC (ADDITIVE DISTANCE) MATRIX (AD)\", save it to the same folder with name `D_hyp_{DEC_method}.phylip` if the input file is `D_hyp.phylip` and `D_metric_{}.phylip` for the other case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DEC = 'NJ' # choose 'NJ' or 'UNJ'\n",
    "\n",
    "Hyp_dec_path = './embeddings/{}/seed_{}_curv_{}/D_hyp_{}.phylip'.format(DATASET,SEED,CURV,DEC)\n",
    "Metric_dec_path = './embeddings/{}/seed_{}_curv_{}/D_metric_{}.phylip'.format(DATASET,SEED,CURV,DEC)\n",
    "\n",
    "def read_phylip2npy(Hyp_NJ_path):\n",
    "    f = open(Hyp_NJ_path,'r')\n",
    "    Array = np.loadtxt(f,)\n",
    "    Array = Array[:,1:]\n",
    "#     print(Array)\n",
    "    f.close()\n",
    "    return Array\n",
    "\n",
    "D_metric_dec_output = read_phylip2npy(Metric_dec_path)\n",
    "loss_metric_dec = ((D_metric_dec_output-D_metric)**2).sum().sum()**(1/2)\n",
    "D_hyp_dec_output = read_phylip2npy(Hyp_dec_path)\n",
    "loss_hyp_dec = ((D_hyp_dec_output-D_metric)**2).sum().sum()**(1/2)\n",
    "print('l2 loss for Direct: ',loss_metric_dec)\n",
    "print('l2 loss for HyperAid: ',loss_hyp_dec)\n",
    "print('Gain (%): ',(loss_metric_dec/loss_hyp_dec-1)*100)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "hyphc",
   "language": "python",
   "name": "hyphc"
  },
  "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}