Multilevel Monte Carlo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Multilevel Monte Carlo Solver (Tensorflow v.1.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import packages, reset the computational graph, start interactive Tensorflow session and set random seed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import tensorflow as tf\n",
    "from timeit import default_timer as timer\n",
    "import numpy as np\n",
    "tf.reset_default_graph()\n",
    "sess = tf.InteractiveSession()\n",
    "tf.set_random_seed(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define Parameter of the PDE (example: constant coefficient parabolic equation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with tf.variable_scope('PDE-Parameter'):\n",
    "    d=10 #dimension\n",
    "    mu=tf.zeros((1,d),name='Mu') #explicit solution available\n",
    "    sigma=tf.multiply(1/d,tf.ones((d,d)),name='Sigma') #Sigma TRANSPOSED "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define Parameter of the Algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "T=1.0 #time, where the solution is evaluated\n",
    "a=0.0 \n",
    "b=1.0 #approx. the solution of the PDE in [a,b]^d\n",
    "l_max=1 #2^(l_max) is the number of time steps in the finest level\n",
    "l_min=0 #2^(l_min) is the number of time steps in the coarsest level\n",
    "n_hidlayer=3 #number of hidden layers for the neural nets\n",
    "n_neurons=[16,8,4] #number of neurons for the hidden layers\n",
    "K=256 #multiplicator for the batch-size\n",
    "n_valid=10000 #number of samples for validation of the expected value/loss\n",
    "n_acc=100 #number of samples for measuring the accuracy (by Monte-Carlo-Simulation)\n",
    "valid_steps=500 #every valid_steps the loss is computed \n",
    "start_lr=9e-5 #starting learning rate of the gradient descent\n",
    "decay_steps=2000 #learning rate decays exponentially \n",
    "w_initializer=tf.contrib.layers.variance_scaling_initializer() #initializer of the weights (alternativ: tf.contrib.layers.xavier_initializer() for Sigmoid, tf.constant_initializer(0), tf.contrib.layers.variance_scaling_initializer() for ReLu, tf.truncated_normal_initializer(0,1))\n",
    "b_initializer=tf.zeros_initializer() #initializer of the biases (alternativ: tf.zeros_initializer(), tf.constant_initializer(0.01) for ReLu)\n",
    "activation=tf.nn.elu #activation function (alternativ: tf.nn.sigmoid, tf.nn.elu, tf.nn.relu, tf.nn.tanh)\n",
    "adam_eps=1e-5 #epsilon of the Adam optimizer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Functions for defining a neuronal network, attaching summaries to the Tensors (for TensorBoard visualization), feeding with data, calculating the loss and optimizing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def variable_summaries(var): \n",
    "    with tf.variable_scope('Summaries'):\n",
    "        mean = tf.reduce_mean(var)\n",
    "        tf.summary.scalar('Average', mean)\n",
    "        tf.summary.histogram('Histogram', var)\n",
    "\n",
    "def phi(x,lvl): \n",
    "    with tf.variable_scope('Euler-Scheme'):\n",
    "        with tf.variable_scope('Fine_Realization'):\n",
    "            N_fine=tf.constant(2**lvl,name='Steps_fine',dtype=tf.int32) \n",
    "            h_fine=tf.divide(T,tf.cast(N_fine,tf.float32),name='Step-Size_fine')\n",
    "            dw_fine=tf.random_normal(shape=[N_fine,tf.shape(x)[0],d],mean=0.0,stddev=tf.sqrt(h_fine),name='DW')\n",
    "            count=tf.constant(0,name='Count')\n",
    "            def scheme_fine(i,realisation):\n",
    "                realisation+=mu*h_fine+tf.matmul(dw_fine[i,:,:],sigma)\n",
    "                return [i+1,realisation]\n",
    "            _, y_fine=tf.while_loop(lambda i,realisation: i<N_fine, scheme_fine, loop_vars=[count,x],name='Euler-Loop_fine')\n",
    "        if lvl==l_min:\n",
    "            phi=tf.reduce_sum(y_fine**2, axis=1, keepdims=True, name='Phi')\n",
    "            return phi\n",
    "        else:\n",
    "            with tf.variable_scope('Coarse_Realization'):\n",
    "                N_coarse=tf.constant(2**(lvl-1),name='Steps_coarse',dtype=tf.int32)\n",
    "                h_coarse=tf.divide(T,tf.cast(N_coarse,tf.float32),name='Step-Size_coarse')\n",
    "                dw_coarse=dw_fine[0::2,:,:]+dw_fine[1::2,:,:]\n",
    "                def scheme_coarse(i,realisation):\n",
    "                    realisation+=mu*h_coarse+tf.matmul(dw_coarse[i,:,:],sigma)\n",
    "                    return [i+1,realisation]                 \n",
    "                _, y_coarse=tf.while_loop(lambda i,realisation: i<N_coarse, scheme_coarse, loop_vars=[count,x],name='Euler-Loop_coarse')\n",
    "            phi=tf.subtract(tf.reduce_sum(y_fine**2, axis=1, keepdims=True),tf.reduce_sum(y_coarse**2, axis=1, keepdims=True),name='Phi')\n",
    "            return phi\n",
    "            \n",
    "def nn(input_layer, num_hidlayer, num_neurons, level, weight_initializer, bias_initializer, start_learn_rate, training):\n",
    "    name_suffix=str(level)\n",
    "    with tf.variable_scope('Network_'+name_suffix):           \n",
    "        with tf.variable_scope('Target'):\n",
    "            is_validation=tf.placeholder(tf.bool,name='Is_Validation') #\n",
    "            z=tf.cond(is_validation,lambda: tf.placeholder(tf.float32,[None,1],name='Z-Input'),lambda: phi(input_layer,level),name='Network-Target')\n",
    "        with tf.variable_scope('Normalization'):\n",
    "            prev_output=input_layer-(a+b)/2 \n",
    "        for n in range(num_hidlayer):\n",
    "            with tf.variable_scope('Hidden_Layer%d' %(n+1)):\n",
    "                prev_output=tf.contrib.layers.fully_connected(prev_output,num_neurons[n],activation_fn=activation,normalizer_fn=tf.contrib.layers.batch_norm,normalizer_params={'is_training':training,'updates_collections':'updates'},weights_initializer=weight_initializer)\n",
    "                variable_summaries(prev_output)\n",
    "        with tf.variable_scope('Output_Layer'):\n",
    "            output=tf.contrib.layers.fully_connected(prev_output,1,activation_fn=None,normalizer_fn=None,weights_initializer=weight_initializer,biases_initializer=bias_initializer)\n",
    "            variable_summaries(output)\n",
    "        with tf.variable_scope('Losses'):\n",
    "            delta=tf.clip_by_value(z-output,-100.0,100.0) \n",
    "            loss=tf.reduce_mean(delta**2, name='Loss') \n",
    "            tf.summary.scalar('Loss-Summary', loss)\n",
    "    with tf.name_scope('Train_'+name_suffix):            \n",
    "        global_step=tf.Variable(0, trainable=False, name='Global_Step') \n",
    "        learn_rate=tf.train.exponential_decay(start_learn_rate, global_step,decay_steps, 0.85, staircase=True, name='Learn_Rate')\n",
    "        optimizer=tf.train.AdamOptimizer(learn_rate,epsilon=adam_eps, name='Adam') \n",
    "        var_list=tf.get_collection(tf.GraphKeys.TRAINABLE_VARIABLES, scope='Network_'+name_suffix)\n",
    "        update_ops = tf.get_collection('updates', scope='Network_'+name_suffix)\n",
    "        with tf.control_dependencies(update_ops):\n",
    "            training=optimizer.minimize(loss,global_step=global_step,var_list=var_list,name='Minimizer')\n",
    "        summaries=tf.get_collection(tf.GraphKeys.SUMMARIES, scope='Network_'+name_suffix)\n",
    "        merged=tf.summary.merge(summaries, name='Merged') \n",
    "    return output, z\n",
    "\n",
    "def uniform_input():\n",
    "    batch_size=tf.placeholder(tf.int32,shape=[],name='Batchsize')\n",
    "    return tf.random_uniform([batch_size,d],minval=a,maxval=b,name='Xi-Input')\n",
    "\n",
    "def build_model():\n",
    "    with tf.variable_scope('Input'):\n",
    "        is_training=tf.placeholder(tf.bool,name='Is_Training') \n",
    "        nn_input=tf.cond(is_training,lambda: uniform_input(),lambda: tf.placeholder(tf.float32,[None,d],name='X-Input'),name='Network-Input')\n",
    "    out=[nn(nn_input,n_hidlayer,n_neurons,l, w_initializer, b_initializer, start_lr, is_training) for l in range(l_min,l_max+1)]\n",
    "    nn_outputs=[lst[0] for lst in out]\n",
    "    phi_outputs=[lst[1] for lst in out]\n",
    "    with tf.variable_scope('Accuracy'):\n",
    "        u_real=tf.placeholder(tf.float32,[None,1],name='U_Real') \n",
    "        u_approx=tf.add_n(nn_outputs,name='U_Approximation') \n",
    "        abs_diff=tf.abs(u_approx-u_real,name='Absolute-Error')\n",
    "        max_error=tf.reduce_max(abs_diff,name='Max-Error')\n",
    "        tf.summary.scalar('Max_Error-Summary', max_error)\n",
    "        l2_error=tf.sqrt(tf.reduce_mean(abs_diff**2)*(b-a)**d,name='L2-Error')\n",
    "        tf.summary.scalar('L2_Error-Summary', l2_error)\n",
    "        summaries_acc=tf.get_collection(tf.GraphKeys.SUMMARIES, scope='Accuracy')\n",
    "        merged_acc=tf.summary.merge(summaries_acc, name='Merged_Acc') \n",
    "\n",
    "def valid_accuracy_data(num_validation,num_accuracy):\n",
    "    dictionary={'Input/Is_Training:0': True, 'Input/Network-Input/Batchsize:0': num_accuracy}\n",
    "    accuracy_data=[sess.run('Input/Network-Input/Merge:0',feed_dict=dictionary)]\n",
    "    dictionary['Input/Network-Input/Batchsize:0']=num_validation\n",
    "    dictionary.update({'Network_'+str(l)+'/Target/Is_Validation:0': False for l in range(l_min,l_max+1)})\n",
    "    fetch=['Input/Network-Input/Merge:0']\n",
    "    fetch.append(['Network_'+str(l)+'/Target/Network-Target/Merge:0' for l in range(l_min,l_max+1)])\n",
    "    validation_data=sess.run(fetch,feed_dict=dictionary)\n",
    "    dictionary['Input/Is_Training:0']=False    \n",
    "    del dictionary['Input/Network-Input/Batchsize:0']\n",
    "    MC_mean=np.empty((num_accuracy,1))\n",
    "    for sample in range(num_accuracy):\n",
    "        dictionary.update({'Input/Network-Input/X-Input:0': np.tile(accuracy_data[0][sample],(10000,1))})\n",
    "        MC_mean[sample]=sum([np.mean(multilevel_real) for multilevel_real in sess.run(fetch[1],feed_dict=dictionary)])\n",
    "    accuracy_data.append(MC_mean)\n",
    "    return validation_data, accuracy_data\n",
    "    \n",
    "def trainNN(level,n_iterations,validation_data,accuracy_data,batch_sizes):\n",
    "    name_suffix=str(level)\n",
    "    scope='Network_'+name_suffix+'/'\n",
    "    scope_train='Train_'+name_suffix+'/'\n",
    "    scope_target=scope+'Target/'\n",
    "    scope_z_target=scope_target+'Network-Target/'\n",
    "    valid_dictionary={'Input/Network-Input/X-Input:0': validation_data[0], scope_z_target+'Z-Input:0': validation_data[1][level-l_min], 'Input/Is_Training:0': False, scope_target+'Is_Validation:0': True}\n",
    "    accuracy_dictionary={'Input/Network-Input/X-Input:0': accuracy_data[0], 'Input/Is_Training:0': False, 'Accuracy/U_Real:0': accuracy_data[1]}\n",
    "    glob_iterations=sess.run(scope_train+'Global_Step:0')\n",
    "    for iteration in range(glob_iterations,glob_iterations+n_iterations): \n",
    "        if ((iteration)%valid_steps)==0: \n",
    "            summary, rate, valid_loss=sess.run([scope_train+'Merged/Merged:0',scope_train+'Learn_Rate:0', scope+'Losses/Loss:0'], feed_dict=valid_dictionary)\n",
    "            summary_acc, l2_err, max_err=sess.run(['Accuracy/Merged_Acc/Merged_Acc:0','Accuracy/L2-Error:0', 'Accuracy/Max-Error:0'], feed_dict=accuracy_dictionary)\n",
    "            writer[level-l_min].add_summary(summary, iteration) \n",
    "            writer[level-l_min].add_summary(summary_acc, iteration)  \n",
    "            print('Network: %d, Iteration: %d, Loss: %.8f, Max. Error: %.4f, L2-Error: %.4f, Learning Rate: %.1E' %(level,iteration,valid_loss,max_err,l2_err,rate))\n",
    "        sess.run(scope_train+'Minimizer:0', feed_dict={'Input/Is_Training:0': True, scope_target+'Is_Validation:0': False, 'Input/Network-Input/Batchsize:0': batch_sizes})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Build the neuronal networks, initialize variables and prepare summaries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "build_model()\n",
    "print('Model build!')\n",
    "valid_data, acc_data=valid_accuracy_data(n_valid,n_acc)\n",
    "print('Data for validation and accuracy measurements computed!')\n",
    "tf.global_variables_initializer().run() \n",
    "print('Variables initialized!')\n",
    "count=0 #find a new directory for the summary logs\n",
    "while True:\n",
    "    count+=1\n",
    "    if not tf.gfile.Exists('logs/log'+str(count)):\n",
    "        dir='logs/log'+str(count)\n",
    "        break\n",
    "writer = [tf.summary.FileWriter(dir+'/network_'+str(l)) for l in range(l_min,l_max+1)] \n",
    "writer[0].add_graph(sess.graph) \n",
    "print('Tensorboard summaries prepared!')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Train the networks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "n_iter=30000 #number of gradient descents\n",
    "start = timer()\n",
    "for l in range(l_min,l_max+1):    \n",
    "    trainNN(l,n_iter,valid_data,acc_data,K*2**(l_max+l_min-l))\n",
    "end = timer()\n",
    "print('Elapsed minutes to train the networks: ',(end - start)/60)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Optional: Train further selected networks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "l=0\n",
    "n_iter=20000\n",
    "batch_sz=K*2**(l_max+l_min-l)\n",
    "trainNN(l,n_iter,valid_data,acc_data,batch_sz)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical experiments"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test accuracy against analytic solution (only for mu=0!!)  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def real_value(X):\n",
    "    sigma_matrix=np.transpose(sess.run('PDE-Parameter/Sigma:0'))\n",
    "    return (np.sum(X**2,axis=1,keepdims=True)+T*np.trace(sigma_matrix@np.transpose(sigma_matrix)))\n",
    "\n",
    "def max_l2_analytic_error(n_tests): \n",
    "    dictionary={'Input/Is_Training:0': True, 'Input/Network-Input/Batchsize:0': n_tests}\n",
    "    xi=sess.run('Input/Network-Input/Merge:0',feed_dict=dictionary)\n",
    "    del dictionary['Input/Network-Input/Batchsize:0']\n",
    "    real_u=real_value(xi)\n",
    "    dictionary.update({'Input/Network-Input/X-Input:0': xi, 'Input/Is_Training:0': False, 'Accuracy/U_Real:0': real_u})\n",
    "    max_error, l2_error=sess.run(['Accuracy/Max-Error:0','Accuracy/L2-Error:0'],feed_dict=dictionary)\n",
    "    return max_error, l2_error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "max_analytic, l2_analytic=max_l2_analytic_error(500000)\n",
    "print('Max. Error: %.4f L2_Error: %.4f' %(max_analytic,l2_analytic))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test accuracy at given points "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def u_eval(X):\n",
    "    return sess.run('Accuracy/U_Approximation:0', feed_dict={'Input/Network-Input/X-Input:0': X, 'Input/Is_Training:0': False})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_0=np.zeros((1,d))\n",
    "x_0[0,0]=0.2\n",
    "tests=[np.random.uniform(a,b,(1,d)),np.random.uniform(a,b,(1,d)),np.ones((1,d)),x_0,np.zeros((1,d)),0.1*np.ones((1,d)),0.2*np.ones((1,d)),0.3*np.ones((1,d)),0.4*np.ones((1,d)),0.5*np.ones((1,d)),0.6*np.ones((1,d)),np.array([np.arange(1,d+1,1)**(-1.0)])]\n",
    "for count, test in enumerate(tests, start=1):\n",
    "    print('\\n test: ',count,', x=',test,'\\n')\n",
    "    print('Real value (only for mu=0): ',real_value(test))\n",
    "    print('Neuronal Network approx.:   ',u_eval(test))\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Close the tensorflow session"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sess.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Optional: Show the Neuronal networks and Summaries on tensorboard"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!tensorboard --logdir=\"logs\" --port=6006"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "now open: http://localhost:6006/ or http://PC-NAME:6006/ (z.B.: http://Julius-PC:6006/) restart the kernel afterwards"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.5.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}