gravity.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tensor Field Networks\n",
    "\n",
    "Implementation of Newtonian gravity demonstration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import random\n",
    "from math import sqrt\n",
    "import tensorflow as tf\n",
    "import matplotlib\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from tensorfieldnetworks import utils \n",
    "from tensorfieldnetworks import layers\n",
    "from tensorfieldnetworks.utils import EPSILON, FLOAT_TYPE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def R(inputs, nonlin=tf.nn.relu, hidden_dim=None, output_dim=1, weights_initializer=None, biases_initializer=None):\n",
    "    with tf.variable_scope(None, \"radial_function\", values=[inputs]):\n",
    "        if weights_initializer is None:\n",
    "            weights_initializer = tf.contrib.layers.xavier_initializer()\n",
    "        if biases_initializer is None:\n",
    "            biases_initializer = tf.constant_initializer(0.)\n",
    "        input_dim = inputs.get_shape()[-1]\n",
    "        if hidden_dim is None:\n",
    "            hidden_dim = input_dim\n",
    "\n",
    "        w1 = tf.get_variable('weights1', [hidden_dim, input_dim], dtype=FLOAT_TYPE,\n",
    "                             initializer=weights_initializer)\n",
    "        b1 = tf.get_variable('biases1', [hidden_dim], dtype=FLOAT_TYPE, initializer=biases_initializer)\n",
    "\n",
    "        w2 = tf.get_variable('weights2', [output_dim, hidden_dim], dtype=FLOAT_TYPE,\n",
    "                             initializer=weights_initializer)\n",
    "        b2 = tf.get_variable('biases2', [output_dim], dtype=FLOAT_TYPE, initializer=biases_initializer)\n",
    "\n",
    "        hidden_layer = nonlin(b1 + tf.tensordot(inputs, w1, [[2], [1]]))\n",
    "        radial = b2 + tf.tensordot(hidden_layer, w2, [[2], [1]])\n",
    "\n",
    "        # [N, N, output_dim]\n",
    "        return radial\n",
    "\n",
    "\n",
    "def unit_vectors(v, axis=-1):\n",
    "    return v / utils.norm_with_epsilon(v, axis=axis, keep_dims=True)\n",
    "\n",
    "\n",
    "def Y_2(rij):\n",
    "    # rij : [N, N, 3]\n",
    "    # x, y, z : [N, N]\n",
    "    x = rij[:, :, 0]\n",
    "    y = rij[:, :, 1]\n",
    "    z = rij[:, :, 2]\n",
    "    r2 = tf.maximum(tf.reduce_sum(tf.square(rij), axis=-1), EPSILON)\n",
    "    # return : [N, N, 5]\n",
    "    output = tf.stack([x * y / r2,\n",
    "                       y * z / r2,\n",
    "                       (-tf.square(x) - tf.square(y) + 2. * tf.square(z)) / (2 * sqrt(3) * r2),\n",
    "                       z * x / r2,\n",
    "                       (tf.square(x) - tf.square(y)) / (2. * r2)],\n",
    "                      axis=-1)\n",
    "    return output\n",
    "\n",
    "\n",
    "def F_0(inputs, nonlin=tf.nn.relu, hidden_dim=None, output_dim=1,\n",
    "        weights_initializer=None, biases_initializer=None):\n",
    "    # [N, N, output_dim, 1]\n",
    "    with tf.variable_scope(None, \"F_0\", values=[inputs]):\n",
    "        return tf.expand_dims(\n",
    "            R(inputs, nonlin=nonlin, hidden_dim=hidden_dim, output_dim=output_dim,\n",
    "              weights_initializer=weights_initializer, biases_initializer=biases_initializer),\n",
    "            axis=-1)\n",
    "\n",
    "\n",
    "def F_1(inputs, rij, nonlin=tf.nn.relu, hidden_dim=None, output_dim=1,\n",
    "        weights_initializer=None, biases_initializer=None):\n",
    "    with tf.variable_scope(None, \"F_1\", values=[inputs]):\n",
    "        # [N, N, output_dim]\n",
    "        radial = R(inputs, nonlin=nonlin, hidden_dim=hidden_dim, output_dim=output_dim,\n",
    "                   weights_initializer=weights_initializer, biases_initializer=biases_initializer)\n",
    "        # Mask out for dij = 0\n",
    "        dij = tf.norm(rij, axis=-1)\n",
    "        condition = tf.tile(tf.expand_dims(dij < EPSILON, axis=-1), [1, 1, output_dim])\n",
    "        masked_radial = tf.where(condition, tf.zeros_like(radial), radial)\n",
    "        # [N, N, output_dim, 3]\n",
    "        return tf.expand_dims(unit_vectors(rij), axis=-2) * tf.expand_dims(masked_radial, axis=-1)\n",
    "\n",
    "\n",
    "def F_2(inputs, rij, nonlin=tf.nn.relu, hidden_dim=None, output_dim=1,\n",
    "        weights_initializer=None, biases_initializer=None):\n",
    "    with tf.variable_scope(None, \"F_2\", values=[inputs]):\n",
    "        # [N, N, output_dim]\n",
    "        radial = R(inputs, nonlin=nonlin, hidden_dim=hidden_dim, output_dim=output_dim,\n",
    "                   weights_initializer=weights_initializer, biases_initializer=biases_initializer)\n",
    "        # Mask out for dij = 0\n",
    "        dij = tf.norm(rij, axis=-1)\n",
    "        condition = tf.tile(tf.expand_dims(dij < EPSILON, axis=-1), [1, 1, output_dim])\n",
    "        masked_radial = tf.where(condition, tf.zeros_like(radial), radial)\n",
    "        # [N, N, output_dim, 5]\n",
    "        return tf.expand_dims(Y_2(rij), axis=-2) * tf.expand_dims(masked_radial, axis=-1)\n",
    "\n",
    "\n",
    "def filter_0(layer_input,\n",
    "             rbf_inputs,\n",
    "             nonlin=tf.nn.relu,\n",
    "             hidden_dim=None,\n",
    "             output_dim=1,\n",
    "             weights_initializer=None,\n",
    "             biases_initializer=None):\n",
    "    with tf.variable_scope(None, \"F0_to_L\", values=[layer_input]):\n",
    "        # [N, N, output_dim, 1]\n",
    "        F_0_out = F_0(rbf_inputs, nonlin=nonlin, hidden_dim=hidden_dim, output_dim=output_dim,\n",
    "                      weights_initializer=weights_initializer, biases_initializer=biases_initializer)\n",
    "        # [N, output_dim]\n",
    "        input_dim = layer_input.get_shape().as_list()[-1]\n",
    "        # Expand filter axis \"j\"\n",
    "        cg = tf.expand_dims(tf.eye(input_dim), axis=-2)\n",
    "        # L x 0 -> L\n",
    "        return tf.einsum('ijk,abfj,bfk->afi', cg, F_0_out, layer_input)\n",
    "\n",
    "\n",
    "def filter_1_output_0(layer_input,\n",
    "                      rbf_inputs,\n",
    "                      rij,\n",
    "                      nonlin=tf.nn.relu,\n",
    "                      hidden_dim=None,\n",
    "                      output_dim=1,\n",
    "                      weights_initializer=None,\n",
    "                      biases_initializer=None):\n",
    "    with tf.variable_scope(None, \"F1_to_0\", values=[layer_input]):\n",
    "        # [N, N, output_dim, 3]\n",
    "        F_1_out = F_1(rbf_inputs,\n",
    "                      rij,\n",
    "                      nonlin=nonlin,\n",
    "                      hidden_dim=hidden_dim,\n",
    "                      output_dim=output_dim,\n",
    "                      weights_initializer=weights_initializer,\n",
    "                      biases_initializer=biases_initializer)\n",
    "        # [N, output_dim, 3]\n",
    "        if layer_input.get_shape().as_list()[-1] == 1:\n",
    "            raise ValueError(\"0 x 1 cannot yield 0\")\n",
    "        elif layer_input.get_shape().as_list()[-1] == 3:\n",
    "            # 1 x 1 -> 0\n",
    "            cg = tf.expand_dims(tf.eye(3), axis=0)\n",
    "            return tf.einsum('ijk,abfj,bfk->afi', cg, F_1_out, layer_input)\n",
    "        else:\n",
    "            raise NotImplementedError(\"Other Ls not implemented\")\n",
    "\n",
    "\n",
    "def filter_1_output_1(layer_input,\n",
    "                      rbf_inputs,\n",
    "                      rij,\n",
    "                      nonlin=tf.nn.relu,\n",
    "                      hidden_dim=None,\n",
    "                      output_dim=1,\n",
    "                      weights_initializer=None,\n",
    "                      biases_initializer=None):\n",
    "    with tf.variable_scope(None, \"F1_to_1\", values=[layer_input]):\n",
    "        # [N, N, output_dim, 3]\n",
    "        F_1_out = F_1(rbf_inputs,\n",
    "                      rij,\n",
    "                      nonlin=nonlin,\n",
    "                      hidden_dim=hidden_dim,\n",
    "                      output_dim=output_dim,\n",
    "                      weights_initializer=weights_initializer,\n",
    "                      biases_initializer=biases_initializer)\n",
    "        # [N, output_dim, 3]\n",
    "        if layer_input.get_shape().as_list()[-1] == 1:\n",
    "            # 0 x 1 -> 1\n",
    "            cg = tf.expand_dims(tf.eye(3), axis=-1)\n",
    "            return tf.einsum('ijk,abfj,bfk->afi', cg, F_1_out, layer_input)\n",
    "        elif layer_input.get_shape().as_list()[-1] == 3:\n",
    "            # 1 x 1 -> 1\n",
    "            return tf.einsum('ijk,abfj,bfk->afi', utils.get_eijk(), F_1_out, layer_input)\n",
    "        else:\n",
    "            raise NotImplementedError(\"Other Ls not implemented\")\n",
    "\n",
    "\n",
    "def filter_2_output_2(layer_input,\n",
    "                      rbf_inputs,\n",
    "                      rij,\n",
    "                      nonlin=tf.nn.relu,\n",
    "                      hidden_dim=None,\n",
    "                      output_dim=1,\n",
    "                      weights_initializer=None,\n",
    "                      biases_initializer=None):\n",
    "    with tf.variable_scope(None, \"F2_to_2\", values=[layer_input]):\n",
    "        # [N, N, output_dim, 3]\n",
    "        F_2_out = F_2(rbf_inputs,\n",
    "                      rij,\n",
    "                      nonlin=nonlin,\n",
    "                      hidden_dim=hidden_dim,\n",
    "                      output_dim=output_dim,\n",
    "                      weights_initializer=weights_initializer,\n",
    "                      biases_initializer=biases_initializer)\n",
    "        # [N, output_dim, 5]\n",
    "        if layer_input.get_shape().as_list()[-1] == 1:\n",
    "            # 0 x 2 -> 2\n",
    "            cg = tf.expand_dims(tf.eye(5), axis=-1)\n",
    "            return tf.einsum('ijk,abfj,bfk->afi', cg, F_2_out, layer_input)\n",
    "        else:\n",
    "            raise NotImplementedError(\"Other Ls not implemented\")\n",
    "\n",
    "\n",
    "def self_interaction_layer_without_biases(inputs, output_dim, weights_initializer=None, biases_initializer=None):\n",
    "    # input has shape [N, C, 2L+1]\n",
    "    # input_dim is number of channels\n",
    "    if weights_initializer is None:\n",
    "        weights_initializer = tf.orthogonal_initializer()\n",
    "    if biases_initializer is None:\n",
    "        biases_initializer = tf.constant_initializer(0.)\n",
    "\n",
    "    with tf.variable_scope(None, \"self_interaction_layer\", values=[inputs]):\n",
    "        input_dim = inputs.get_shape().as_list()[-2]\n",
    "        w_si = tf.get_variable('weights', [output_dim, input_dim], dtype=FLOAT_TYPE,\n",
    "                               initializer=weights_initializer)\n",
    "        # [N, output_dim, 2l+1]\n",
    "        return tf.transpose(tf.einsum('afi,gf->aig', inputs, w_si), perm=[0, 2, 1])\n",
    "\n",
    "\n",
    "def self_interaction_layer_with_biases(inputs, output_dim, weights_initializer=None, biases_initializer=None):\n",
    "    # input has shape [N, C, 2L+1]\n",
    "    # input_dim is number of channels\n",
    "    if weights_initializer is None:\n",
    "        weights_initializer = tf.orthogonal_initializer()\n",
    "    if biases_initializer is None:\n",
    "        biases_initializer = tf.constant_initializer(0.)\n",
    "\n",
    "    with tf.variable_scope(None, \"self_interaction_layer\", values=[inputs]):\n",
    "        input_dim = inputs.get_shape().as_list()[-2]\n",
    "        w_si = tf.get_variable('weights', [output_dim, input_dim], dtype=FLOAT_TYPE,\n",
    "                               initializer=weights_initializer)\n",
    "        b_si = tf.get_variable('biases', [output_dim], initializer=biases_initializer,\n",
    "                               dtype=FLOAT_TYPE)\n",
    "\n",
    "        # [N, output_dim, 2l+1]\n",
    "        return tf.transpose(tf.einsum('afi,gf->aig', inputs, w_si) + b_si, perm=[0, 2, 1])\n",
    "\n",
    "\n",
    "def convolution(input_tensor_list, rbf, unit_vectors, weights_initializer=None, biases_initializer=None):\n",
    "    with tf.variable_scope(None, \"convolution\", values=[input_tensor_list]):\n",
    "        output_tensor_list = {0: [], 1: []}\n",
    "        for key in input_tensor_list:\n",
    "            with tf.variable_scope(None, \"L\" + str(key), values=input_tensor_list[key]):\n",
    "                for i, tensor in enumerate(input_tensor_list[key]):\n",
    "                    output_dim = tensor.get_shape().as_list()[-2]\n",
    "                    with tf.variable_scope(None, 'tensor_' + str(i), values=[tensor]):\n",
    "                        tensor = tf.identity(tensor, name=\"in_tensor\")\n",
    "                        if True:\n",
    "                            # L x 0 -> L\n",
    "                            tensor_out = filter_0(tensor,\n",
    "                                                  rbf,\n",
    "                                                  output_dim=output_dim,\n",
    "                                                  weights_initializer=weights_initializer,\n",
    "                                                  biases_initializer=biases_initializer)\n",
    "                            m = 0 if tensor_out.get_shape().as_list()[-1] == 1 else 1\n",
    "                            tensor_out = tf.identity(tensor_out, name=\"F0_to_L_out_tensor\")\n",
    "                            output_tensor_list[m].append(tensor_out)\n",
    "                        if key is 1:\n",
    "                            # L x 1 -> 0\n",
    "                            tensor_out = filter_1_output_0(tensor,\n",
    "                                                           rbf,\n",
    "                                                           unit_vectors,\n",
    "                                                           output_dim=output_dim,\n",
    "                                                           weights_initializer=weights_initializer,\n",
    "                                                           biases_initializer=biases_initializer)\n",
    "                            m = 0 if tensor_out.get_shape().as_list()[-1] == 1 else 1\n",
    "                            tensor_out = tf.identity(tensor_out, name=\"F1_to_0_out_tensor\")\n",
    "                            output_tensor_list[m].append(tensor_out)\n",
    "                        if key is 0 or key is 1:\n",
    "                            # L x 1 -> 1\n",
    "                            tensor_out = filter_1_output_1(tensor,\n",
    "                                                           rbf,\n",
    "                                                           unit_vectors,\n",
    "                                                           output_dim=output_dim,\n",
    "                                                           weights_initializer=weights_initializer,\n",
    "                                                           biases_initializer=biases_initializer)\n",
    "                            m = 0 if tensor_out.get_shape().as_list()[-1] == 1 else 1\n",
    "                            tensor_out = tf.identity(tensor_out, name=\"F1_to_1_out_tensor\")\n",
    "                            output_tensor_list[m].append(tensor_out)\n",
    "        return output_tensor_list\n",
    "\n",
    "\n",
    "def self_interaction(input_tensor_list, output_dim, weights_initializer=None, biases_initializer=None):\n",
    "    with tf.variable_scope(None, \"self_interaction\", values=[input_tensor_list]):\n",
    "        output_tensor_list = {0: [], 1: []}\n",
    "        for key in input_tensor_list:\n",
    "            with tf.variable_scope(None, \"L\" + str(key), values=input_tensor_list[key]):\n",
    "                for i, tensor in enumerate(input_tensor_list[key]):\n",
    "                    with tf.variable_scope(None, 'tensor_' + str(i), values=[tensor]):\n",
    "                        if key == 0:\n",
    "                            tensor_out = self_interaction_layer_with_biases(tensor,\n",
    "                                                                            output_dim,\n",
    "                                                                            weights_initializer=weights_initializer,\n",
    "                                                                            biases_initializer=biases_initializer)\n",
    "                        else:\n",
    "                            tensor_out = self_interaction_layer_without_biases(tensor,\n",
    "                                                                               output_dim,\n",
    "                                                                               weights_initializer=weights_initializer,\n",
    "                                                                               biases_initializer=biases_initializer)\n",
    "                        m = 0 if tensor_out.get_shape().as_list()[-1] == 1 else 1\n",
    "                        output_tensor_list[m].append(tensor_out)\n",
    "        return output_tensor_list\n",
    "\n",
    "\n",
    "def nonlinearity(input_tensor_list, nonlin=tf.nn.elu, biases_initializer=None):\n",
    "    with tf.variable_scope(None, \"nonlinearity\", values=[input_tensor_list]):\n",
    "        output_tensor_list = {0: [], 1: []}\n",
    "        for key in input_tensor_list:\n",
    "            with tf.variable_scope(None, \"L\" + str(key), values=input_tensor_list[key]):\n",
    "                for i, tensor in enumerate(input_tensor_list[key]):\n",
    "                    with tf.variable_scope(None, 'tensor_' + str(i), values=[tensor]):\n",
    "                        if key == 0:\n",
    "                            tensor_out = utils.rotation_equivariant_nonlinearity(tensor,\n",
    "                                                                                 nonlin=nonlin,\n",
    "                                                                                 biases_initializer=biases_initializer)\n",
    "                        else:\n",
    "                            tensor_out = utils.rotation_equivariant_nonlinearity(tensor,\n",
    "                                                                                 nonlin=nonlin,\n",
    "                                                                                 biases_initializer=biases_initializer)\n",
    "                        m = 0 if tensor_out.get_shape().as_list()[-1] == 1 else 1\n",
    "                        output_tensor_list[m].append(tensor_out)\n",
    "        return output_tensor_list\n",
    "\n",
    "\n",
    "def concatenation(input_tensor_list):\n",
    "    output_tensor_list = {0: [], 1: []}\n",
    "    for key in input_tensor_list:\n",
    "        with tf.variable_scope(None, \"L\" + str(key), values=input_tensor_list[key]):\n",
    "            # Concatenate along channel axis\n",
    "            # [N, channels, M]\n",
    "            output_tensor_list[key].append(tf.concat(input_tensor_list[key], axis=-2))\n",
    "    return output_tensor_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "FLOAT_TYPE = tf.float32\n",
    "EPSILON = 1e-8\n",
    "\n",
    "\n",
    "def get_eijk():\n",
    "    \"\"\"\n",
    "    Constant Levi-Civita tensor\n",
    "\n",
    "    Returns:\n",
    "        tf.Tensor of shape [3, 3, 3]\n",
    "    \"\"\"\n",
    "    eijk_ = np.zeros((3, 3, 3))\n",
    "    eijk_[0, 1, 2] = eijk_[1, 2, 0] = eijk_[2, 0, 1] = 1.\n",
    "    eijk_[0, 2, 1] = eijk_[2, 1, 0] = eijk_[1, 0, 2] = -1.\n",
    "    return tf.constant(eijk_, dtype=FLOAT_TYPE)\n",
    "\n",
    "\n",
    "def norm_with_epsilon(input_tensor, axis=None, keep_dims=False):\n",
    "    \"\"\"\n",
    "    Regularized norm\n",
    "\n",
    "    Args:\n",
    "        input_tensor: tf.Tensor\n",
    "\n",
    "    Returns:\n",
    "        tf.Tensor normed over axis\n",
    "    \"\"\"\n",
    "    return tf.sqrt(tf.maximum(tf.reduce_sum(tf.square(input_tensor), axis=axis, keep_dims=keep_dims), EPSILON))\n",
    "\n",
    "\n",
    "def ssp(x):\n",
    "    \"\"\"\n",
    "    Shifted soft plus nonlinearity.\n",
    "\n",
    "    Args:\n",
    "        x: tf.Tensor\n",
    "\n",
    "    Returns:\n",
    "        tf.Tensor of same shape as x \n",
    "   \"\"\"\n",
    "    return tf.log(0.5 * tf.exp(x) + 0.5)\n",
    "\n",
    "\n",
    "def rotation_equivariant_nonlinearity(x, nonlin=ssp, biases_initializer=None):\n",
    "    \"\"\"\n",
    "    Rotation equivariant nonlinearity.\n",
    "\n",
    "    The -1 axis is assumed to be M index (of which there are 2 L + 1 for given L).\n",
    "\n",
    "    Args:\n",
    "        x: tf.Tensor with channels as -2 axis and M as -1 axis.\n",
    "\n",
    "    Returns:\n",
    "        tf.Tensor of same shape as x with 3d rotation-equivariant nonlinearity applied.\n",
    "    \"\"\"\n",
    "    if biases_initializer is None:\n",
    "        biases_initializer = tf.constant_initializer(0.)\n",
    "    shape = x.get_shape().as_list()\n",
    "    channels = shape[-2]\n",
    "    representation_index = shape[-1]\n",
    "\n",
    "    biases = tf.get_variable('biases',\n",
    "                             [channels],\n",
    "                             dtype=FLOAT_TYPE,\n",
    "                             initializer=biases_initializer)\n",
    "\n",
    "    if representation_index == 1:\n",
    "        return nonlin(x)\n",
    "    else:\n",
    "        norm = norm_with_epsilon(x, axis=-1)\n",
    "        nonlin_out = nonlin(tf.nn.bias_add(norm, biases))\n",
    "        factor = tf.divide(nonlin_out, norm)\n",
    "        # Expand dims for representation index.\n",
    "        return tf.multiply(x, tf.expand_dims(factor, axis=-1))\n",
    "    \n",
    "\n",
    "\n",
    "def difference_matrix(geometry):\n",
    "    \"\"\"\n",
    "    Get relative vector matrix for array of shape [N, 3].\n",
    "\n",
    "    Args:\n",
    "        geometry: tf.Tensor with Cartesian coordinates and shape [N, 3]\n",
    "\n",
    "    Returns:\n",
    "        Relative vector matrix with shape [N, N, 3]\n",
    "    \"\"\"\n",
    "    # [N, 1, 3]\n",
    "    ri = tf.expand_dims(geometry, axis=1)\n",
    "    # [1, N, 3]\n",
    "    rj = tf.expand_dims(geometry, axis=0)\n",
    "    # [N, N, 3]\n",
    "    rij = ri - rj\n",
    "    return rij\n",
    "\n",
    "\n",
    "def distance_matrix(geometry):\n",
    "    \"\"\"\n",
    "    Get relative distance matrix for array of shape [N, 3].\n",
    "\n",
    "    Args:\n",
    "        geometry: tf.Tensor with Cartesian coordinates and shape [N, 3]\n",
    "\n",
    "    Returns:\n",
    "        Relative distance matrix with shape [N, N]\n",
    "    \"\"\"\n",
    "    # [N, N, 3]\n",
    "    rij = difference_matrix(geometry)\n",
    "    # [N, N]\n",
    "    dij = norm_with_epsilon(rij, axis=-1)\n",
    "    return dij\n",
    "\n",
    "\n",
    "def random_rotation_matrix(numpy_random_state):\n",
    "    \"\"\"\n",
    "    Generates a random 3D rotation matrix from axis and angle.\n",
    "\n",
    "    Args:\n",
    "        numpy_random_state: numpy random state object\n",
    "\n",
    "    Returns:\n",
    "        Random rotation matrix.\n",
    "    \"\"\"\n",
    "    rng = numpy_random_state\n",
    "    axis = rng.randn(3)\n",
    "    axis /= np.linalg.norm(axis) + EPSILON\n",
    "    theta = 2 * np.pi * rng.uniform(0.0, 1.0)\n",
    "    return rotation_matrix(axis, theta)\n",
    "\n",
    "\n",
    "def rotation_matrix(axis, theta):\n",
    "    return scipy.linalg.expm(np.cross(np.eye(3), axis * theta))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# radial basis functions\n",
    "rbf_low = 0.\n",
    "rbf_high = 2.\n",
    "rbf_count = 30\n",
    "rbf_spacing = (rbf_high - rbf_low) / rbf_count\n",
    "centers = tf.lin_space(rbf_low, rbf_high, rbf_count)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING:tensorflow:From C:\\Users\\jacob\\Documents\\ece594n_hw\\tensorfieldnetworks\\tensorfieldnetworks\\utils.py:32: calling reduce_sum (from tensorflow.python.ops.math_ops) with keep_dims is deprecated and will be removed in a future version.\n",
      "Instructions for updating:\n",
      "keep_dims is deprecated, use keepdims instead\n"
     ]
    }
   ],
   "source": [
    "# [N, 3]\n",
    "r = tf.placeholder(FLOAT_TYPE, shape=[None, 3])\n",
    "\n",
    "# [N, 1, 1]\n",
    "masses = tf.placeholder(FLOAT_TYPE, shape=[None, 1, 1])\n",
    "\n",
    "# [N, N, 3]\n",
    "rij = utils.difference_matrix(r)\n",
    "\n",
    "# [N, N]\n",
    "dij = utils.distance_matrix(r)\n",
    "\n",
    "# rbf : [N, N, rbf_count]\n",
    "gamma = 1. / rbf_spacing\n",
    "rbf = tf.exp(-gamma * tf.square(tf.expand_dims(dij, axis=-1) - centers))\n",
    "\n",
    "# shifted softplus\n",
    "nonlin = utils.ssp\n",
    "\n",
    "hidden_dim = rbf_count\n",
    "output_dim = 1\n",
    "\n",
    "# single layer, 0 -> 1\n",
    "with tf.variable_scope(None, \"layer1\", values=[masses]):\n",
    "    # [N, 1, 3]\n",
    "    output = layers.filter_1_output_1(masses, rbf, rij, nonlin=nonlin, \n",
    "                                      hidden_dim=hidden_dim, output_dim=output_dim, \n",
    "                                      weights_initializer=tf.glorot_uniform_initializer(),\n",
    "                                      biases_initializer=tf.glorot_uniform_initializer())\n",
    "    \n",
    "# accel : [N, 1, 3]\n",
    "accel = tf.placeholder(FLOAT_TYPE, shape=[None, 1, 3])\n",
    "\n",
    "# loss : []\n",
    "loss = tf.nn.l2_loss(accel - output)\n",
    "\n",
    "optim = tf.train.AdamOptimizer(learning_rate=1.e-3)\n",
    "\n",
    "train_op = optim.minimize(loss)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for plotting\n",
    "\n",
    "# x : []\n",
    "x = tf.placeholder(FLOAT_TYPE, shape=(None))\n",
    "\n",
    "# rbf_plot : [rbf_count]\n",
    "rbf_plot = tf.exp(-gamma * tf.square(x - centers))\n",
    "\n",
    "for v in tf.global_variables():\n",
    "    if 'layer1/F1_to_1/F_1/radial_function/weights1' in v.name and \"Adam\" not in v.name:\n",
    "        w1 = v\n",
    "    elif 'layer1/F1_to_1/F_1/radial_function/biases1' in v.name and \"Adam\" not in v.name:\n",
    "        b1 = v\n",
    "    elif 'layer1/F1_to_1/F_1/radial_function/weights2' in v.name and \"Adam\" not in v.name:\n",
    "        w2 = v\n",
    "    elif 'layer1/F1_to_1/F_1/radial_function/biases2' in v.name and \"Adam\" not in v.name:\n",
    "        b2 = v\n",
    "        \n",
    "hidden_layer_plot = nonlin(b1 + tf.tensordot(rbf_plot, w1, [[0], [1]]))\n",
    "R_plot = b2 + tf.tensordot(hidden_layer_plot, w2, [[0], [1]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def accelerations(points, masses=None):\n",
    "    \"\"\"\n",
    "    inputs:\n",
    "    -points: a list of 3-tuples of point coordinates\n",
    "    -masses: a list (of equal length N) of masses\n",
    "    \n",
    "    returns: \n",
    "    -shape [N, 3] numpy array of accelerations under Newtonian gravity\n",
    "    \"\"\"\n",
    "    accels = []\n",
    "    if masses is None:\n",
    "        masses = [1.0 for _ in range(len(points))]\n",
    "    for i, ri_ in enumerate(points):\n",
    "        accel_vec = np.array((0., 0., 0.))\n",
    "        for j, rj_ in enumerate(points):\n",
    "            rij_ = ri_ - rj_\n",
    "            dij_ = np.linalg.norm(rij_)\n",
    "            if (ri_ != rj_).any():\n",
    "                accel_update = -rij_ / (np.power(dij_, 3) + EPSILON) * masses[j]\n",
    "                accel_vec += accel_update\n",
    "        accels.append(accel_vec)\n",
    "    assert len(accels) == len(points)\n",
    "    return np.array(accels)\n",
    "\n",
    "\n",
    "def random_points_and_masses(max_points=10, min_mass=0.5, max_mass=2.0, \n",
    "                             max_coord=rbf_high, min_separation=0.5):\n",
    "    \"\"\"\n",
    "    returns:\n",
    "    -shape [N, 3] numpy array of points, where N is between 2 and max_points\n",
    "    -shape [N] numpy array of masses\n",
    "    \"\"\"\n",
    "    num_points = random.randint(2, max_points)\n",
    "    candidate_points = []\n",
    "    for point in range(num_points):\n",
    "        candidate_points.append(\n",
    "            np.array([random.uniform(-max_coord, max_coord) for _ in range(3)]))\n",
    "    \n",
    "    # remove points that are closer than min_separation\n",
    "    output_points = []\n",
    "    for point in candidate_points:\n",
    "        include_point = True\n",
    "        for previous_point in output_points:\n",
    "            if np.linalg.norm(point - previous_point) < min_separation:\n",
    "                include_point = False\n",
    "        if include_point:\n",
    "            output_points.append(point)\n",
    "    \n",
    "    points_ = np.array(output_points)\n",
    "    masses_ = np.random.rand(len(output_points)) * (max_mass - min_mass) + min_mass\n",
    "    return points_, masses_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'xrange' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-9-2516faad2555>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m     13\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mstep\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mmax_steps\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     14\u001b[0m     rand_points, rand_masses = random_points_and_masses(max_points=10, \n\u001b[1;32m---> 15\u001b[1;33m                                                         min_separation=min_separation)\n\u001b[0m\u001b[0;32m     16\u001b[0m     \u001b[0mrand_masses_m\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mreshape\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mrand_masses\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m[\u001b[0m\u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mrand_masses\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     17\u001b[0m     sess.run(train_op, feed_dict={r: rand_points, \n",
      "\u001b[1;32m<ipython-input-7-925b69458547>\u001b[0m in \u001b[0;36mrandom_points_and_masses\u001b[1;34m(max_points, min_mass, max_mass, max_coord, min_separation)\u001b[0m\n\u001b[0;32m     33\u001b[0m     \u001b[0mnum_points\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mrandom\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mrandint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mmax_points\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     34\u001b[0m     \u001b[0mcandidate_points\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 35\u001b[1;33m     \u001b[1;32mfor\u001b[0m \u001b[0mpoint\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mxrange\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnum_points\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     36\u001b[0m         candidate_points.append(\n\u001b[0;32m     37\u001b[0m             np.array([random.uniform(-max_coord, max_coord) for _ in range(3)]))\n",
      "\u001b[1;31mNameError\u001b[0m: name 'xrange' is not defined"
     ]
    }
   ],
   "source": [
    "max_steps = 1001\n",
    "validation_size = 1000\n",
    "print_freq = 1000\n",
    "\n",
    "F1_y_vals = []\n",
    "x_vals = [index * rbf_high / rbf_count for index in range(rbf_count + 1)]\n",
    "\n",
    "sess = tf.Session()\n",
    "sess.run(tf.global_variables_initializer())\n",
    "\n",
    "# training\n",
    "min_separation = 0.5\n",
    "for step in range(max_steps):\n",
    "    rand_points, rand_masses = random_points_and_masses(max_points=10, \n",
    "                                                        min_separation=min_separation)\n",
    "    rand_masses_m = np.reshape(rand_masses, [len(rand_masses), 1, 1])\n",
    "    sess.run(train_op, feed_dict={r: rand_points, \n",
    "                                  masses: rand_masses_m, \n",
    "                                  accel: np.expand_dims(accelerations(rand_points,\n",
    "                                                                      rand_masses), axis=-2)})\n",
    "    \n",
    "    # print performance on different random point set every print_freq steps\n",
    "    if step % print_freq == 0:\n",
    "        loss_sum = 0.\n",
    "        for _ in range(validation_size):\n",
    "            validation_points, validation_masses = random_points_and_masses(max_points=50,\n",
    "                                                                            min_separation=min_separation)\n",
    "            validation_masses_m = np.reshape(validation_masses, [len(validation_masses), 1, 1])\n",
    "            validation_loss = sess.run(loss, feed_dict={r: validation_points, \n",
    "                                                        masses: validation_masses_m, \n",
    "                                                        accel: np.expand_dims(\n",
    "                                                            accelerations(validation_points, \n",
    "                                                                          validation_masses), axis=-2)})\n",
    "            loss_sum += validation_loss\n",
    "        print(\"Step %d: validation loss = %.3f\" % (step, loss_sum / validation_size))\n",
    "\n",
    "        # for plotting radial function\n",
    "        y_vals = []\n",
    "        for x_val in x_vals:\n",
    "            y_vals.append(sess.run(R_plot, feed_dict={x: x_val}))\n",
    "        F1_y_vals.append((step, y_vals))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVoAAADFCAYAAAAVFjpUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJztnXd4VVXW/z8rBAglBCEgQoBQBAUCiBQRFVAEhUFwxEFER8aCdXjxZ+/ldSwjoqNYXhwUxBlEsaAoijTHQjEgiIAICAxNadKbJOv3xz43uWk3J8k9aazP8+wn55y9z97fe3Kzss8ua4mqYhiGYQRHTEkLMAzDKO+YoTUMwwgYM7SGYRgBY4bWMAwjYMzQGoZhBIwZWsMwjIAxQ2sYhhEwZmgNwzACplCGVkT+EG0hhmEY5ZXC9mg7RVWFYRhGOUZsC65hGEawxPopJCJnAsnh5VX1jYA0GYZhlCvyNbQiMhFoBiwB0rzLCpihNQzD8EG+QwcishJopTbGYBiGUSj8TIb9ANQLWohhGEZ5xc8YbSKwQkQWAkdCF1X1osBUGYZhlCP8GNqHgxZhGIZRnvG1vEtETiRz7exCVd0WqCrDMIxyRL5jtCLyJ2AhcCnwJ2CBiAwKWphhGEZ5wc+qg6XA+aFerIjUAWaqarti0GcYhlHm8bPqICbbUMFOn/cZhmEY+JsM+1REPgMmeeeDgU+Ck2QYhlG+8DsZdgnQzTv9UlXfD1RVAUlMTNTk5OSSllGqWLVqFQAtW7YsYSWGUXZZtGjRDlWtU9R6yoVTmY4dO2pqampJyyhV9OjRA4C5c+eWqA7DKMuIyCJV7VjUevIcOhCRr1T1LBHZh/NtkJEFqKrWKGrjhmEYxwN5GlpVPcv7GV98cgzDMMofftbRTvRzzTAMw8gdP6sOWoefiEgscHowcoxocdttt5W0BMMwPCKN0d4D3AtUEZG9ocvAUWBsMWgzikD//v1LWoJhGB55Dh2o6hPe+OzTqlrDS/GqWltV7ylGjUYhWLVqVcYSL8MwShY/O7wWikhC6EREaorIwAA1hdq5QERWicgaEbk76PbKG9dffz3XX399ScswDAN/hvYhVd0TOlHV3cBDwUkCEakAvAhcCLQChohIqyDbNAzDCApfvg5yueYrqGMR6AysUdWfVfUo8BYwIOA2DcMwAsGPoU0VkdEi0sxLo4FFAetqAGwMO9/kXTMMwyhz+DG0f8WtNJjspSPAzUGK8oOIDBeRVBFJ3b59e0nLMQzDyJN8hwBU9QBQ3JNRm4GGYedJ3rUMVHUs3jKzjh07ln2HDVHm/vvvL2kJhlF2+P13OHwYjh6FuDioVi2q1edraEWkBXA7kBxeXlXPjaqSrHwLnCwiTXAG9jLg8gDbK3f06tWrpCUYhn/274fNm52xq14dmjXLmr9yJcyd6/JbtYI+fbLmT58Ob74JR45A375w9dVZ88eMgX/8wxnSESMgtKFH1RnZ22+HF15w1+67z91/9GjUPp6fSa13gFeAfwJpUWs5Aqp6TERuAT4DKgCvqery4mi7vLBkyRIA2rdvX8JKjHLBwYPO2O3fD7Gx0LmzM1DHjrm0Zg2MGwcHDkD9+nDddZl5aWmwYAE8+qgzlKedBg89BOnpLi89HWbNgiefdG2dfTbcdZczgqE0Ywa8+KLL79nT6VB196rCJ5/Av//t8vfudRrT0jLTjBlOI7iyaWlO/++/u2vffZf5WefNc/dHET+1HVPVl6Paqg9U9RPMwXihGTlyJGBuEo9L0tKgQoXM8/R0+O03mDwZduxwxm7oUNdjO3LE/dy0Ce6/3xnU+Hh44AF3/dgxZ4w2bIBnnnH11a0LN96Ytc2NG+G119xxgwZQtWrW/J9/hi1b3PHatc6wAohATIy7P8SOHbBihcsLpd27M/P37nV1heeHDCbAoUPuvEIFqFzZ/awR5mwwPh5at4aKFV2qVAl++QVSU91xx47wl7+440ceKfjzzwU/hvYjEbkJeB83EQaAqu6KigLDMPJHFd57D3buhG3bXI/x8GHXgzxwwBnSK66AffuckRkzxhnRQ4dcucOH4amnXF0h4xLOwYOwbp07PnDA1VGpknuNj41194SIiYHzz8+8HhvrjGfI0FarBtde667HxjpD9+23MNHzRVW3rjPqMTEuAXz5pev1VqkC3bvDHXdk1Zea6rTExTlDeMUVWfN79XJDBpUrQ/PmrtccTp8+8OCDLr9mTUhIyJrfrRu88or/30cB8ROccV0ul1VVmwYjqeCY4++cmOPvUsr+/a731LCh+6MHZ9QOHnTGc+1aZ0xfeMEZupAhPXAARo50xhPgzjudUQqhCn/7m+vNgnuNr1nTlalSxdXVrVtm+a1bndGqXNnlHTrkenrgyh88mFX3jh3OuFavDklJMGlS1vydO+H1111Ptl49+OMfs+YfOuR6oVWquDpqlA131oE7/g6hqk2K2ohhHBds3Qo//eRekzt2dEbn0CFntA4ehOHDYfVqV/b++6FWLXc99No7e7YzWOAmd5KSXO+wWjXXA6tZE3791eWfcw6ceqrLq1rVpVdfdRpE4MIL4aSTsuq79lpnWBMSXNvhvdpq1VyvMT4+0+CGk5iYdRwzO7VruwmlvKhSJecE13GEn1UHf87tuqq+EX05hlGKCPU0Qz3KtWvhxx/dK3bLls5gHTqUmV55JfP1+/LL4eSTs9aXnp617iZNnIGsUsUZumnTMg3txRe7HmQ4333nepa1arlX4/r1s+YvXJjZW4zJZYn8q6/m/VlF4HTzfhoUfsZoO4UdxwHnAYsBM7SlmMcff7ykJZRu0tPda/zu3S7t2QPbt8OqVW52unp1OOGErPe89x4sW+aOL7vMzX5XqZL5in7yyZmG9pRT4Kqr3PWQMV2zxk061asHZ50Ff/hD1vpfeMFNQCUmQtNcRub+/vfInykpqXDPwggcP0MHfw0/F5GaON8DRinmzDPPLGkJpYODB50BDU+//OJ6juFjnOBe3b/80h0PHQqDBmW+uler5nq1IUPbrJkbCghnxQrX42zYEDp0gEaNsuaPG+cmhERy19q1a9E/r1EqKcxisQOAjduWcr755hvgODK4hw87AxpuULdtc8YxxI8/wuefw65d0K8fPPxw5gx0zZpuWVLI0MbEuPWc4XTp4q41beoMaXbuuCPnbHk4ockv47jDzxjtR2RGwY3BuS18O0hRRtG59957gXK66uDIETfps3Wrm8nessX1UA8fhvXrnZE94QS48kqoU8ctJ6pTxxnRyZNdHfv2uQmrcFq1gsaN3SRTp045mmXQIJcMo4BECmVTWVWPAKPCLh8DNqjqpsCVGQa4BfNbt7rtmeFGNbQsMSHBza63a+euh9aKNm8OU6ZkrSs02RMbm3l/ON26OUNtGFEmUo92HtABuFZVrywmPcbxzr59bpfQxo1u4mjLlsy1oVWrunHQdeuc8V282E1ahThyxBnR0JbQPXuyLkw/6SQ3xtqyZdYF+IYRMJEMbSURuRw4U0T+mD1TVd8LTpZxXJCe7taFhgzrxo2ZWy1jY93ypTPOcJNLSUnO0CYmul1Q4Hq5LVtm1le5Mgwb5pY3deiQ+371Nm0C/1iGkZ1IhvYGYChQE8geUlUBM7RGwUhLcz3RDRvcK/p//5u50yk+3s3S168PS5fCzJnwz3+6ZVLhnH02fPihO/7qq6yGFiKvFTWMEiJPQ6uqXwFfiUiqqo4rRk1GFHjuuedKWoJ7hd+82RnWDRtcjzXkei4x0fUuGzd2qUYNt+zpkkvcelVwBjW7ob3uOujd2+2Hb2Vh5IyyQb6+DsoC5uuglJCe7gzrunUubdzojC3AiSc6g5qc7Hqumze73VTZ3ThOmOBe/8FNTn31VXF+AsPIQrH5OjDKJjNnzgQCdgCu6sZYQ4Z1w4bMoYB69dzyqeRkZ2BDmwOWLnVelr780m0xnTEja539+rntpxdd5I4NoxxQPnq0DRtq6rhx0Lat6znltfPmOCIQ712qbglVyLCuW+d6peCGApo0cSk5Oac/0hD//a/LD33vfv7Z3WMYpZDAe7QiksvWl0xUdXFRG8+j3adxk29HgbXAX1R1d8SbKlaE+fPhm2/c4vS2bSElJafPSaPg7NnjjGHIsO7b564nJLiJqJBxze727uBB50rv8suzbnVt1AguuMD1ZPv3d5sMDKOcE2no4JkIeQoEFTPsc+AeL5zNU8A9wF0R76hVy7loW74cvv/ezVjPnOl6Tm3bukmTuLiA5JYz9u93KwJChnWX59+9WrVMo9qkidt5ldebw9ixzsnyr786V3xXZluGPXp05kYDwzgOiLTqoGdxCglrN3zQbj7gb89j1apu22SnTs44LFvmjO6HH8LHH0OLFm63UHKyM8w2vOBe33ftcq/zoRRy01e5sntWnTu7vf116vh/ZmvXZvpNffXVnIY2+0oCwyjn+JoME5E2OB8HGd3CYvJHezUwucB31arllv+cc45bt/n99663u3Kly69RwxmR5GTXO6tZ8/gwvOnpzvFKuGHdv9/lVaniXus7dHCTV/Xr5+7T1A933ul8s9aoAQMHOoN+PDxfw8gDP6FsHgJ64AztJ8CFwFeqWmjvGiIyE6iXS9Z9qjrVK3Mf0BH4o+YiUkSGA8MBGjVqdPqGDRsiNxqayAm9Fq9fn+nZKSEhcxKnfn1nqKMcBbO4WbVsGezcSctq1Vzv8pdf3M+QN/+aNZ1hbdzY/UxMLJgxVIUPPoCXX4apU3O6HFywwPkfsCEbowwTrckwP4Z2GdAO+E5V24nIicCbqnp+xBuLIkpkGHA9cJ6qHsyneOHW0ao6L08hw7thQ2acpJgYNwaZmJg11alT+gzH0aNuwmrXrkyD+ssvbptq6HcbF+eWW514otvO2qhR0WM2DRzoDCzAqFFw221Fq88wSiHFuY72kKqmi8gxEakBbAMaFrXhvBCRC4A7ge5+jGwRGnIrFOrWdeOQqs5/6bZtznnz9u3u55o1mU5NINPzfrhD6FCqWjXrcWFevVUzQzwfPZoZTmXPHpf27s083rMn56x9rVpw4ol8dOAA1KxJ/yFDXI892q/uvXtnGtoxY1zgwPAQ14ZhZODH0KZ6URVeBRYB+3GevYJiDFAZ+FyccZivqjcE2J5DxPX4Tjwx6/X0dNc73LEjM+3e7XqQGzc6I5jXW0FMjDM+4Sn7tbS0TIMaSpHeMqpUyXRUHdq6GjqvWzfDufQz3jra/jfeGIWHkwvXXusMbN++cNddZmQNIwJ+Qtnc5B2+IiKfAjVU9fugBKlq86DqLhQxMS7CZ+3aOR2YgDOKhw5lDQsdinqalpZ/io1164BDqVKlnMch41qjRtbIpcXBsWPw4otwzTVZXRJWquQmGcv4WLZhFAeRNiycoqo/5rZxQUQ6BLVhocwhkhnuuU6dklYTXdauhSuucJtBli9362PDMSNrGL6I9JdyG3AduW9cCHLDglFamDfPGVlw62GHDHGRXw3DKBCRNixc5/20v6zjlaFD4aOPnNvCRx5x65INwygweS7vyi2qQjilKcKCuUnMycaNGwFo2LCIC0R27XLL30LxtgzjOKI4lneFoirUBc4EZnvnPYFvsAgLpZoCG9ivv4a5c+G++7Jer1XLJcMwCk2koYO/AIjIDKCVqm71zk8CxheLOqPQTPbCag8ePDj/wv/5j1umdeCAW872wAMBqzOM4ws/K+obhoysx69Ao4D0GFHi5Zdf5uWXX/ZXeMyYzO3IY8a4tcKGYUQNP+tzZonIZ8Ak73wwMDM4SUax88YbbpfZsmUwe7bbbmwYRtTws2HhFm9i7Gzv0lhVfT9YWUaxEhfnHMT88otFOzCMAPC14txbYWCTX+WFHTty9lqrVDEjaxgBke8YrYicISLfish+ETkqImkisrc4xBkB8P77zqBOn17SSgzjuMFPj3YMcBnwDs4/7J+BFkGKMorOlClTcl786CO49FLnY+Hii52xtZ1ehhE4vvz4qeoaoIKqpqnq68AFwcoyikpiYiKJ2YcH2rVz/mjBef5qYf8vDaM48GNoD4pIJWCJiPxdRG71eZ9RgowfP57x48dnvdioEcyZA336uM0JDRqUhDTDOO7wE2GhMW7tbCXgViABeMnr5ZYKbAtuTnp4/mjnzp1bojoMoyxTLBEWRKQC8LiqDgUOA48UtUGjmLHAiIZR4kQcAlDVNKCxN3RQrIjIbSKiImKr5wvL3r3QtSv8+GNJKzGM4xo/qw5+Br4WkQ+BA6GLqjo6KFEi0hDoDfw3qDbKPWlpLrz64cNw2mnw73+7lQaGYRQ7fia11gLTvLLxYSlInsUFaIw8gGzkzd69cOSIO65Y0RlbwzBKBD9bcIt1XFZEBgCbVXWp2Nhioflk3jz44QcXRPH22yE5uaQlGcZxS4kEfRKRmUC9XLLuA+7FDRvkV8dwYDhAo0bmTCw7VatWdWHUU1Ndj9YwjBIj3+VdxYmIpACzgIPepSRgC9BZVX/J6z5b3pWTl156CYCbbropn5KGYeRFtJZ3laqNB6q6TFXrqmqyqiYDm4AOkYysEcbcubBvHwBvv/02b7/9dsnqMQwDiBxu/AUiTEap6ohAFBmFY/Vq6NcPTjwRJk4saTWGYYQRqUebCiwC4oAOwGovtcftEgscr2dr7v7zIy0NrrwSDh50gRRvvLGkFRmGEUaehlZVJ6jqBKAt0ENVX1DVF4DzcMbWKC1UqAAjRkBCgpv4mjChpBUZhhGGn1UHJwA1gF3eeXXvmlGauPxyOOssmDfP1swaRinDj6F9EvhOROYAApwDPBykKKOQNGrkEuZMxjBKE342LLwuItOBLt6lu2wVgGEYhn/8hLIRoBfQTlWnApVEpHPgyoz8mTfPbbXNhVGjRjFq1KhiFmQYRm74WUf7EtAVGOKd7wNeDEyR4Y9Dh2DAABcx4Y47MtbPhpg2bRrTpk0rIXGGYYTjx9B2UdWbcf5oUdXfKKblXUYEJk6E7dtdj3byZBcy3DCMUokfQ/u75wBcAUSkDpAeqCojfxIToWVLd3zrrebPwDBKMX4M7fPA+0BdEfkb8BXweKCqjPz54x9hxQr48EPnocswjFKLn1UH/xKRRbiNCgIMVNWVgSsz8icmBvr3zzWrSpUqxSzGMIy88OsmcTWwN1ReRBqpqkU/KMVMnz69pCUYhuGRr6EVkb8CD+Ei4abherWK25prFDe7d0PNmiWtIiK///47mzZt4vDhwyUtxTB8ERcXR1JSEhUDmuvw06P9H6Clqu4MRIHhn7VroU0b+NOf4LbboG3e/+v+93//F4AHHniguNRlsGnTJuLj40lOTsaiZBilHVVl586dbNq0iSZNmgTShp/JsI3AnkBaNwrGs8+6YItvvAF33hmx6KxZs5g1a1YxCcvK4cOHqV27thlZo0wgItSuXTvQNzC/UXDnisjHwJHQxSCj4Bq5oApr1mSe33FHyWnxgRlZoywR9PfVj6H9r5cqUUwbFbxx4ZtxY8Ifq2rk7tvxgAh8+iksWADvvAPnnlvSigzD8ElEQ+ttVIhX1duLSQ8i0hMYgPOtcERE6hZX22WCLl1cMgyjzBDR0Kpqmoh0Ky4xHjcCT6rqEU/DtmJuv1xQu3btkpZQLvnggw/4+OOP2bt3L9dccw29e+cbsNkwfE2GLRGRD0XkShH5YygFqKkFcLaILBCRL0SkU4BtlVveffdd3n333ZKWUaq5+uqrqVu3Lm3atMmRd8MNN/D111/nuD5w4EBeffVVXnnlFSZPnlwcMo1ygB9DGwfsBM4F+nvpD0VpVERmisgPuaQBuF52LeAM4A7gbcllpFpEhotIqoikbt++vShySjfTpsGFF8Ls2W5CzIgaw4YN49NPP801b/78+Zxxxhk5rqelpQHw2GOPcfPNNweqzyg/+NmC+5doN6qqvfLKE5EbgfdUVYGFIpIOJAJZrKmqjgXGAnTs2LH8WqCnn4b//MdNhI0a5dbP+uCee+4B4IknnghSXZnmnHPOYf369Tmur1y5khYtWlChQgUALr30UmrVqsXSpUvp168f+/fv58ILL6RDhw7FrNgoq/hx/B0nIjeLyEsi8looBajpA6Cn13YL3EqH4zMS7po18OWX7jg2FgYP9n3rvHnzmDdvXkDCyjfTp0/nggsuyDhftmwZJ554IvPnzychIYGZM2cyZcoUXnnllRJUaZQl/AwdTATqAX2AL4AknPPvoHgNaCoiPwBvAVd5vdvjj+bNYdUquOkmuPpqSEoqaUWF4+GH3fI0EXecndtuy8x/5pmc+cOHZ+aPHeu72V69etGmTZscaerUqRHv++yzzzIM7eHDh9m1axcPPvggACNGjGDRokW88sor3HDDDb61GMc3ftbRNlfVS0VkgKpOEJF/A18GJUhVjwJXBFV/mePkk+FFC2hRGGbOnFngew4ePMju3bupX78+AMuXL6dLly7Exvr1v2QYOfHl+Nv7uVtE2gAJgK1tNcolc+bMoWfPnhnny5Yto20EnxKG4Qc//6bHisgJwAPAh0B14MFAVRlFJqk0DTM8/HDuQwYhnnkm9yGDEGPHFmjIwC9Dhgxh7ty57Nixg6SkJB555BG+++47Bg0alFFm2bJldO5ssUiNoiHlYfizY8eOmpqaWtIyosuGDS7wYoyfl47SxcqVKzn11FNLWkah6NChAwsWLAjMXZ5Resnteysii1S1Y1Hr9rPq4EQRGSci073zViJyTVEbNiKQlgYdOkC9ejB0aJ4hxY3os3jxYjOyRtTx010aD3wG1PfOfwJGBiXIABYtgl27XJTb2bMhPr7AVYwcOZKRI+3XZBilAT+GNlFV38aLfKuqx3BetYyg2LwZ6tRxx717u2VNBWTJkiUsWbIkysIMwygMfibDDohIbTLDjZ+BOQIPlosvhgEDYOlSCyNuGOUAP4b2/+FWGzQTka+BOsCgyLcYRSYmBk47raRVGIYRBfz4OlgsIt2BlrjAjKtU9fd8bjMMwzA88jS0EVwhthARVPW9gDQZUaBFixYlLcEwDI9IPdr+EfIUMEMbBM895yIodOrkHMkUkrEBLPA3DKNw5LnqQFX/EiFdXZwijxvWr4dbb4Uzz4TGjd16WqNQ/O1vf6N169a0bduW9u3bs2DBAgCee+45Dh48GJU2fvzxR7p27UrlypUZNWpUlrxPP/2Uli1b0rx5c5588smM6+vWraNLly40b96cwYMHc/ToUQCOHDnC4MGDad68OV26dMnVfaNRdil7247KMzNmZB63aweeP9TCMHz4cIYPHx4FUWWPefPmMW3aNBYvXsz333/PzJkzadiwIRBdQ1urVi2ef/55br89a0i9tLQ0br75ZqZPn86KFSuYNGkSK1asAOCuu+7i1ltvZc2aNZxwwgmMGzcOgHHjxnHCCSewZs0abr31Vu66666oaDRKB2ZoSxOtW2e6Qzz//CJV9dNPP/HTTz9FSVjZYuvWrSQmJlK5cmUAEhMTqV+/Ps8//zxbtmyhZ8+eGY5jZsyYQdeuXenQoQOXXnop+/fvByA5OZk777yTlJQUOnfuzJrwUO8edevWpVOnTjl2ki1cuJDmzZvTtGlTKlWqxGWXXcbUqVNRVWbPnp3hS+Gqq67igw8+AGDq1KlcddVVAAwaNIhZs2ZRHrbHG45Ik2GXquo7ItJEVdcVp6jjlm7dXFKFY8dKWk10+PRT+OWX6NZZrx6EOebOTu/evXn00Udp0aIFvXr1YvDgwXTv3p0RI0YwevRo5syZQ2JiIjt27OCxxx5j5syZVKtWjaeeeorRo0dn+J5NSEhg2bJlvPHGG4wcOZJp06b5krd58+aMHjQ4Bz8LFixg586d1KxZM8PlYlJSEps3b85xT2xsLAkJCezcuZPExMRCPSKjdBGpR3uP99Mi/BU3IrZRoQhUr16dRYsWMXbsWOrUqcPgwYMZP358jnLz589nxYoVdOvWjfbt2zNhwgQ2bNiQkT9kyJCMnxatwigKkaa1d4rIDKCJiHyYPVNVLwpCkIi0B17BBYU8BtykqguDaMsoBiL0PIOkQoUK9OjRgx49epCSksKECRMYNmxYljKqyvnnn8+kSZNyrSM8Jmgu8UHzpEGDBmzcuDHjfNOmTTRo0IDatWuze/dujh07RmxsbMb18HuSkpI4duwYe/bssZDx5YhIPdp+OL+zO4BncklB8XfgEVVt77X/9wDbKre0b9+e9u3bl7SMEmHVqlWsXr0643zJkiU0btwYgPj4ePbtc5GYzjjjDL7++uuM8dcDBw5kGdcOhROfPHkyXbt29d1+p06dWL16NevWrePo0aO89dZbXHTRRYgIPXv2ZMqUKQBMmDCBAQMGAHDRRRcxYcIEAKZMmcK5555bIONulHJUNWIC6ng/qwPV8ytf1ITzFDbYOx4C/Du/e04//XQt0+zerdqpk+rdd6vOmVPSaorMihUrSrT91NRU7dq1q5566qmakpKiF198sW7fvl1VVZ9//nlt0aKF9ujRQ1VVZ82apR07dtSUlBRNSUnRqVOnqqpq48aN9c4779SUlBTt2LGjrl69Okc7W7du1QYNGmh8fLwmJCRogwYNdM+ePaqq+vHHH+vJJ5+sTZs21cceeyzjnrVr12qnTp20WbNmOmjQID18+LCqqh46dEgHDRqkzZo1006dOunatWsDfUZGTnL73gKpGgW7lq/jby98zUSgFm4L7nZcwMQfgjD8InKqZ2wF1+M+U1U3RLqnzDv+/uAD50gG4PTToSx/Fsq24+8QycnJpKam2mTUcUSQjr99hbIB/p+qzvEa7uFdO7OwjYrITFxk3ezcB5wH3Kqq74rIn4BxQK9c6hgODAdo1KhRYaWUDsLXz/buHZUqr7jCxbd88803o1KfYRiFx4+hrRYysgCqOldEqhWlUVXNYThDiMgbwP94p+8A/8yjjrE4g0/Hjh3L9oLDJ5+EPn2cwR04MCpVbtq0KSr1HK/YziwjmvgxtD+LyAO44QNwocB/Dk4SW4DuwFzgXGB1xNLlgRo1nP9Zb2LEMIzyhR9DezXwCM6JjAJfeteC4jrgHyISCxzGGx4wDMMoq/jxR/sbMKIYtITa+wo4vbjaMwzDCJrC++FCALK0AAATZUlEQVQzis6xY257alJS1KsuyLpPwzCCxZzKlCQLFkDDhs6ZzBNPRLXqJ554gieiXGdZQkQyVl4AHDt2jDp16vCHP/wBgA8//DCL+8Lc2LJlS4YDmNJA9erVgfx17d69m5deeqm4ZPnmnXfe4dRTT81w6FMU+vbty+7duyOWGT9+PFu2bClyW9EgX0MrIt38XDMKQWhZ14oV8HOQ84vHH9WqVeOHH37g0KFDAHz++ecZ213B7cS6++67I9ZRv379jF1cpYn8dJVWQztu3DheffVV5syZk3/hfPjkk0+oWbNmxDJlytACL/i8ZhSUPXvAc+UXrfWzIS655BIuueSSqNZZ1ujbty8ff/wxAJMmTcpwEgPuj/CWW24BYNiwYYwYMYIzzzyTpk2bZhix9evX06ZNm4zyAwcO5Pzzzyc5OZkxY8YwevRoTjvtNM444wx27doFQI8ePQhtntmxYwfJyckFuj+cdevW0bVrV1JSUrj//vszrofrWr58OZ07d6Z9+/a0bduW1atXc/fdd7N27Vrat2/PHXfcwf79+znvvPPo0KEDKSkpTJ06NaOeU089leuuu47WrVvTu3fvjH9Ma9asoVevXrRr144OHTqwdu1aAJ5++mk6depE27Zteeihh3J97pMmTSIlJYU2bdpk+NV99NFH+eqrr7jmmmu44447spSfO3cu55xzDv369aNly5bccMMNpKen51kXuA0lO3bsyPMzTJkyhdTUVIYOHUr79u05dOgQd999N61ataJt27Y5fAgHTl5bxoCuwG3ARlwk3FB6GFgajW1p0UplegvugQOqn36q+ttvUa22e/fu2r1796jW6ZfsWxlDWsLTiy++qKqqBw4cyDX/9ddfV1XV7du358jzQ7Vq1XTp0qV6ySWX6KFDh7Rdu3Y6Z84c7devn6qqvv7663rzzTerqupVV12lgwYN0rS0NF2+fLk2a9ZMVVXXrVunrVu3zijfrFkz3bt3r27btk1r1KihL7/8sqqqjhw5Up999tmMz/rtt99maG/cuHGB7g+nf//+OmHCBFVVHTNmjFarVi2HrltuuUXffPNNVVU9cuSIHjx4MEu+qurvv/+esTV4+/bt2qxZM01PT9d169ZphQoV9LvvvlNV1UsvvVQnTpyoqqqdO3fW9957T1Xd9uADBw7oZ599ptddd52mp6drWlqa9uvXT7/44ossmjdv3qwNGzbUbdu26e+//649e/bU999/P8ezCWfOnDlauXJlXbt2rR47dkx79eql77zzTsS6GjdurNu3b4/4GcLb27Fjh7Zo0ULT09NVVfW3XP7egtyCG6lHWwnn3yAWiA9Le7Fw49GjalW3WSGf1yCj4LRt25b169czadIk+vbtG7HswIEDiYmJoVWrVvz666+5lunZsyfx8fHUqVOHhIQE+vd3YfVSUlJ8bXAo6P1ff/11Ri/8yiuvzLXOrl278vjjj/PUU0+xYcMGqlSpkqOMqnLvvffStm1bevXqxebNmzM+Y5MmTTKcD51++umsX7+effv2sXnzZi72toXHxcVRtWpVZsyYwYwZMzjttNPo0KEDP/74YxbnPQDffvstPXr0oE6dOsTGxjJ06FD+85//5PtsOnfuTNOmTalQoQJDhgzhq6++8l1Xbp8hOwkJCcTFxXHNNdfw3nvvUbVq1Xw1RZM8Vx2o6hfAFyIyXvPxNWAYkZg7d26eeVWrVo2Yn5iYGDE/Py666CJuv/125s6dy86dO/MsF4rGAOQZ2SC8TExMTMZ5TEwMxzxH7bGxsRmvvYcPHy7w/dnJz4PX5ZdfTpcuXfj444/p27cv//d//0fTpk2zlPnXv/7F9u3bWbRoERUrViQ5OTlDW7imChUqZAwd5Iaqcs8993D99ddH1FQYsn/Ogngu8/MZYmNjWbhwIbNmzWLKlCmMGTOG2bNnF15wAfEzRjteRGZnT4ErM4wocPXVV/PQQw+RkpJSLO0lJyezaNEigCJPpHXr1o233noLcMYyN37++WeaNm3KiBEjGDBgAN9//30WV5AAe/bsoW7dulSsWJE5c+ZkcW6eG/Hx8SQlJWWE2Tly5AgHDx6kT58+vPbaaxnhfjZv3sy2bduy3Nu5c2e++OILduzYQVpaGpMmTaJ79+75ftaFCxeybt060tPTmTx5MmeddVah6wr/HKHnsH//fvbs2UPfvn159tlnWbp0qe96ooGfdbTho8ZxwCU4h9xGYUlLg08+gRYtoGnTQKIpnHfeeVGvsyySlJTEiBHFtt+G22+/nT/96U+MHTuWfv36Famuf/zjH1x++eU89dRTGX5rs/P2228zceJEKlasSL169bj33nupVasW3bp1o02bNlx44YXcdddd9O/fn5SUFDp27Mgpp5ySb9sTJ07k+uuv58EHH6RixYq888479O7dm5UrV2as0a5evTpvvvkmdevWzbjvpJNO4sknn6Rnz56oKv369ctTezidOnXilltuYc2aNfTs2ZOLL76YmJiYQtUVYtiwYdxwww1UqVKF6dOnM2DAAA4fPoyqMnr0aN/1RIN83STmepPIQlXtHICeQlHm3CSuXw9Nmrjj+vXBixtVXigPbhKN4mPu3LmMGjXKd0y2oChRN4kiUivsNAa3PTahqA0f14RHp802nmYYRvnDz9DBIpwzGcENGawDrglSVLmnUiU47zxncFu2DKSJCy+8EIDp06cHUr9hRItQbLfyjB+nMk2KQ8hxRY8eLgF4M9TRJtLssWEYxYufoYM44CbgLDLdJL6iqocj3mj4I6Z8uptQVQsuaJQZCjNXVRD8/JW/AbTGbbsd4x1PjHiHcVwTFxfHzp07A//yGkY0UFV27txJXFxcYG34GaNto6qtws7niMiKoAQZZZ+kpCQ2bdrE9u3bS1qKYfgiLi6OpADclYbwY2gXi8gZqjofQES6AEVaSyUil+J8JpwKdFbV1LC8e3CTbWnACFX9rChtlTq2bIGpU90a2lNOgTCPUtEk5A6wJKhYsSJNmtjQvmGE8GNoTwe+EZH/eueNgFUisgxQVW1biHZ/AP4I/F/4RRFpBVyGG56oD8wUkRaqmlaINkonCxfCTTe54z594NNPA2mm2L0TGYaRJ34M7QXRblRVV0Ku+5kHAG+p6hFgnYisAToD86KtocQIX0PbokXJ6TAMo9jwY2gfU9UsroNEZGL2a1GiATA/7HyTd6380KYNXHWVM7innRZYM6F1iUVxyGIYRnTwY2hbh5940WnzDZ4oIjOBerlk3aeqU/3Ji1j/cDIj5B4RkR+KWmeUSAR2+Co5bx5cHWRAYRJFxJ+W4PH/XILHtOSktOiA0qUlKjuK8jS03qTUvUAVEdmL2xkGcBQYm1/FqtqrEHo2Aw3DzpO8a7nVPzakQ0RSo7EfORqYltwxLblTWrSUFh1Q+rREo54819Gq6hOqGg88rao1VDXeS7VV9Z5oNJ4LHwKXiUhlEWkCnAwsDKgtwzCMYsHP0MF0ETkn+0VVzd9teh6IyMW4DRB1gI9FZImq9lHV5SLyNrAC51fh5nK14sAwjOMSP4Y2PJJaHG4VwCLg3MI2qqrvA+/nkfc34G8FrDLfoYxixLTkjmnJndKipbTogHKopcD+aEWkIfCcqh7fIVYNwzB8UhiPJptwO7oMwzAMH+RraEXkBRF53ktjcN67FgcvLaP9C0RklYisEZG7c8mvLCKTvfwFIpIclnePd32ViPQpBi3/T0RWiMj3IjJLRBqH5aWJyBIvfVgMWoaJyPawNq8Ny7tKRFZ76aqAdTwbpuEnEdkdlhftZ/KaiGzLa6mfOJ73tH4vIh3C8qL2THxqGeppWCYi34hIu7C89d71JUWd9faho4eI7An7PTwYlhfxdxuAljvCdPzgfT9qeXlReyZefQ1FZI7397pcRP4nlzLR+77kF48cuCosDQW6RSPOuZ8EVADWAk1x4c+XAq2ylbkJ57YR3Pbdyd5xK698ZaCJV0+FgLX0BKp6xzeGtHjn+4v5uQwDxuRyby3gZ+/nCd7xCUHpyFb+r8BrQTwTr75zgA7AD3nk9wWm45YqngEsiPYzKYCWM0NtABeGtHjn64HEYnomPYBpRf3dRkNLtrL9gdlBPBOvvpOADt5xPPBTLn9DUfu++Bk6mIyb/FoEvKuqX/u4J1p0Btao6s+qehR4C7dNN5wBwATveApwnogIYdt5VXUdENrOG5gWVZ2jqge90/m4dcBB4Oe55EUf4HNV3aWqvwGfU/ht1gXVMQSYVMi28kXdSphdEYoMAN5Qx3ygpoicRHSfiS8tqvqN1xYE+F3x8UzyoijfsWhoCfq7slVVF3vH+4CV5NyFGrXvS56GVkRiReTvuDHZCTi/tBtF5O8iEv2wrbnTANgYdp7bltyMMqp6DNgD1PZ5b7S1hHMN7r9hiDgRSRWR+SIysAg6CqLlEu+VZ4q4ScyC3BtNHXjDKE2A8FD10XwmfshLb7S/KwUl+3dFgRkiskjcDsig6SoiS0VkuoiEdoKW2DMRkao4w/Vu2OXAnom44cbTgAXZsqL2fYm0vOtpXJe6iWfxEZEawCgv5RjTMBwicgXQEQgPQt9YVTeLSFNgtogsU9W1Acr4CJikqkdE5HrcP8tCL8mLApcBUzTruujifialDhHpiTO0Z4VdPst7LnWBz0XkRy3CuvV8WIz7PewXkb7AB7iNQiVJf+BrVQ3v/QbyTESkOs6gj1TVvUWtLy8iDR38AbguZGQBPCE34sYuigM/W3Izyojzw5AA7PR5b7S1ICK9gPuAi9R5IQNAVTd7P38G5uL+gwamRVV3hrX/TzL9U0TzuRSkrsvI9ioY5Wfih7z0Rvu74gsRaYv73QxQ1Z2h62HPZRtuvXlRhrwioqp7VXW/d/wJUFFEEimhZ+IR6bsStWfivZm/C/xLVd/LpUj0vi8RBot/KkxeNBOux/0z7pUzNCDfOluZm8k6Gfa2d9yarJNhP1O0yTA/Wk7DTSCcnO36CUBl7zgRWE0RJhZ8ajkp7PhiYL5mDuSv8zSd4B3XCkqHV+4U3GSGBPVMwupNJu+Jn35kndxYGO1nUgAtjXDzBmdmu14NiA87/ga4IEAd9UK/F5zx+q/3fHz9bqOpxctPwI3jVgv4mQhuOPS5CGWi9n2JJOQD4M+5XL8C+LCoD7wAD6QvbkZwLc7zF8CjuB4juN1q73hf2oVA07B77/PuWwVcWAxaZgK/Aku89KF3/UxgmfdlXQZcUwxangCWe23OAU4Ju/dq73mtAf4SpA7v/GHgyWz3BfFMJgFbgd9x42bXADcAN3j5ArzoaV0GdAzimfjU8k/gt7DvSqp3van3TJZ6v7/7AtZxS9j3ZD5hhj+3322QWrwyw3CT2OH3RfWZeHWGgs1+H/Y76BvU9yXPnWEi0gB4DziEW3EAbtyxCnCxel15wzAMIzL5bsEVkXPJ9Em7QlVnBa7KMAyjHFFgXweGYRhGwSiMrwPDMAyjAJihNQzDCBgztIZhGAFjhtYIjDDvXD+IyEciUrOA9z8sIrd7x496m0EilR8vIoPyKZMc8h4lIh1F5Pl8yl5eEM2GkRtmaI0gOaSq7VW1DW4R+s2FrUhVH1TVmdGTBqqaqqojIhRJBszQGkXGDK1RXMzDc7whItXF+etd7PkYzfAKJSL3ifNb+xVhoZ7De6si8qCIfOv1lMd63tryRERO95ymLCXM2Hu+WKd5x93DfKF+JyLxwJPA2d61W70e7pee7sUicmZYPXM95z0/isi/QppEpJM4f7NLRWShiMSLSAURedr7DN97viiMcowZWiNwRKQCcB4uyjHAYdymlw44H77PeE6WT8dto26P26XTKY8qx6hqJ6+nXAXnlyMSrwN/VdV2EcrcjgsG2h44G7dR527gS69X/iywDTjf0z0YCB92OA0YifOD3BToJiKVcG5G/8dru5dX7zXAHlXt5H3G68RFfTbKKWZojSCpIiJLgF+AE3F+O8FtbXxcRL7HbVtu4OWfDbyvqgfVOTDKK+pCT3HRNJbhPJK1zqMc3rhwTc309DQxj6JfA6NFZIRX/lguZSoCr3rtvoMzqiEWquomVU3HbedMxvXIt6rqt5DhwOUY0Bv4s/dsFuDcepa0xywjQMzQGkFyyOshNsYZ19Br+1BcqPnTvfxfcT4r8kVE4oCXgEGqmgK86vfeSKjqk8C1uB7y1yJySi7FbvW0tsNtR68Ulnck7DiNyC5IBdfDbu+lJqo6o0gfwCjVmKE1Akdd1IkRwG1hriy3qervnj/WUGy1/wADRaSKN0baP5fqQkZ1h+dLNOIqA1XdDewWkZC/16G5lRORZqq6TFWfAr7FeRzbh/PJHCIB10NNB67EhXuJxCrgJBHp5LUR733+z4AbPTd9iEgLEamWT11GGSbSf13DiBqq+p03VDAE+BfwkfcKngr86JVZLCKTcV6atuEMXvZ6dovIq8APuCGJHGVy4S/AayKiQF49x5Ge0U/HeYia7h2neZNo43E96XdF5M/Ap8CBfD7zUREZDLwgIlVw47O9cJ67koHF3qTZdqA4IkwYJYT5OjAMwwgYGzowDMMIGDO0hmEYAWOG1jAMI2DM0BqGYQSMGVrDMIyAMUNrGIYRMGZoDcMwAsYMrWEYRsD8f4Szt2sml7udAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f8eecf2c910>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot radial function and compare with -1/r^2\n",
    "\n",
    "radial_fig = plt.figure(figsize=(5,2.5))\n",
    "ax = radial_fig.add_subplot(1,1,1)\n",
    "\n",
    "min_index_cutoff = int((1 / sqrt(40.) - rbf_low) / rbf_spacing)\n",
    "ax.plot(x_vals[min_index_cutoff:], [-1 / r_**2 for r_ in x_vals[min_index_cutoff:]], \"r:\", lw=3, label=\"$-1/r^2$\")\n",
    "\n",
    "for step, y_vals in F1_y_vals[1:]:\n",
    "    line, = ax.plot(x_vals, y_vals, 'r', alpha=0.5, label=\"Step {}\".format(step))\n",
    "    line.set_ydata(y_vals)\n",
    "\n",
    "ax.plot([min_separation, min_separation], [10,-50], 'k--', label=\"Minimum distance of points\")\n",
    "ax.set_ylabel(\"Output of learned radial function.\")\n",
    "ax.set_xlabel(\"Radial distance\")\n",
    "ax.set_xlim(0., 2.0)\n",
    "ax.set_ylim(-10, 1.)\n",
    "handles, labels = ax.get_legend_handles_labels()\n",
    "ax.legend(handles, labels)\n",
    "radial_fig.canvas.draw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ece594n_hw4",
   "language": "python",
   "name": "ece594n_hw4"
  },
  "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}