{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Example 1: Sandstone Model" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# Importing\n", "import theano.tensor as T\n", "import theano\n", "import sys, os\n", "sys.path.append(\"../GeMpy\")\n", "sys.path.append(\"../\")\n", "# Importing GeMpy modules\n", "import gempy as GeMpy\n", "\n", "# Reloading (only for development purposes)\n", "import importlib\n", "importlib.reload(GeMpy)\n", "\n", "# Usuful packages\n", "import numpy as np\n", "import pandas as pn\n", "\n", "import matplotlib.pyplot as plt\n", "\n", "# This was to choose the gpu\n", "os.environ['CUDA_LAUNCH_BLOCKING'] = '1'\n", "\n", "# Default options of printin\n", "np.set_printoptions(precision = 6, linewidth= 130, suppress = True)\n", "\n", "#%matplotlib inline\n", "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# Importing the data from csv files and settign extent and resolution\n", "geo_data = GeMpy.create_data([696000,747000,6863000,6950000,-20000, 200],[ 50, 50, 50],\n", " path_f = os.pardir+\"/input_data/a_Foliations.csv\",\n", " path_i = os.pardir+\"/input_data/a_Points.csv\")\n", "\n", "# Assigning series to formations as well as their order (timewise)\n", "GeMpy.set_data_series(geo_data, {\"EarlyGranite_Series\":geo_data.formations[-1], \n", " \"BIF_Series\":(geo_data.formations[0], geo_data.formations[1]),\n", " \"SimpleMafic_Series\":geo_data.formations[2]}, \n", " order_series = [\"EarlyGranite_Series\",\n", " \"BIF_Series\",\n", " \"SimpleMafic_Series\"], verbose=0)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": true }, "outputs": [], "source": [ "GeMpy.data_to_pickle(geo_data, 'sandstone')" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "I am in the setting\n", "float32\n", "I am here\n", "[2, 2]\n" ] } ], "source": [ "inter = GeMpy.InterpolatorInput(geo_data)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([1, 2, 3, 4])" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "inter.interpolator.tg.n_formation.get_value()" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([], shape=(100, 0), dtype=float64)" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import numpy as np\n", "np.zeros((100,0))" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "100.0" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "100000/1000" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAakAAAFgCAYAAAABy4YnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd4VFX6wPHvnZ7eExJIIAECSkITCB1NQARXkJXiqqss\nrrqCDTu4qwsqKvpbdV3LIipFVxEFREWRgEgJRES69ASSQEIa6ZOZTPn9ETMS0jHJDMn7eR4fyZ1z\n733nzsx97zn3nHMVu91uRwghhHBBKmcHIIQQQtRFkpQQQgiXJUlKCCGEy5IkJYQQwmVd1knKYrGQ\nkZGBxWJxdihCCCFawGWdpLKyskhISCArK8vZoQghhGgBl3WSEkII0bZJkhJCCOGyJEkJIYRwWZKk\nhBBCuCxJUkIIIVyWJCkhhBAuS5KUEEIIlyVJSohmkJSUxJYtWzCbzfWWu5SHDtjtdlJTU8nNzb3U\n8IS4bEmSEqIZ9OnTh8ceeww/Pz+uu+46Xn75ZXbv3o3Vaq1WbvPmzRQVFdW7LbvdzsmTJ3nvvfe4\n/fbb6dy5M0899RT+/v4t+RaEcEkaZwcgRFvg4eHB119/zbBhw1i/fj3r168HwM/Pj2uuuYaEhATi\n4+P57LPPyM/P56abbnKsW5WUNm/ezObNm/nhhx/IyMhwvD5u3DiWLFmCSiXXlKL9kSQlRDMJDAxk\n/fr1DBkyxDFV1/nz51m1ahWrVq0CwM3NjcLCQkeS2r17N5MnT+bUqVO1bnPIkCGsXLkSnU7XKu9B\nCFcjl2ZCNKMuXbrw7bff4u3tXeO1G264AaPRyPr16x3NgFdddRVr166lQ4cONcrHxMTw1Vdf4eHh\n0eJxC+GqpCbVTthsNnJyckhPT6ewsJCEhIQW25fVaKRg/wEsRUVovL3x7R2L2s0NALPJQurxXMpK\nzbh76IjsHohO37xfw/r2fymaGnOfPn1Ys2YN1113XbWOFF999RUAubm5JCcnM3ToUA4ePMgDDzxQ\nY5LkLl26sH79+jZ9H6o1vgvi8qfYL6W7kYvIyMggISGBjRs30qlTJ2eH41RlZWUcPXqU9PR00tPT\nycjIcPw7PT2dM2fOOE6Y8156k8cevBs3vQar1UpBQQF5eXlkZmVjxJ+CXAsaILKTL0atijKzFW8P\nHX26BwGw73gORaVmxzK3C04s5xI3kr3xe6ym307Oar2O4IRrOOvVnb270jGbLRjLLVhtdrQ6NYOH\ndWHQsMga7+nCZGN3N5DWQU+JYkavcsNaGEB5OTViOJe4kawNGyktK8JitaJRq/Fw96bDmARCRjc9\nMe9JTmPvrnQsFb91gNBo1fQdGE7v3kG1JsOqk++aL1Yx5x/3YbfbiYuLIzk52bGNMWPGcOWVV/Kf\n//zHUavq3Lkzp0+fJjg4mO3bt9OtW7daY2rOk3vVMTaez+eMrYjSyCC8vP2JDe6BQWu4pG02Rn3H\ntV9cRLPuq7kvWkTrkiTlRHa7HavVikbz2wnGaLLUmwTqYrVa+fjjj5k3bx4nTpyos5xao8PDOwDs\nNmyWckpLihzdortfMYQbrr2f8rJiikvOU1JaQJmxiFJrCSUVpXS/ahy+wZ3xMPwWj06rZvTACMbE\ndeZc4kYy162v9v6W//wTOSUllCgelBoCsaq1mC12FBQURYWiKKhUKjp2DuatN5/By8sLqJ7sik0l\nFJtLsKgVjkX58nOIO6U5xZSnqFAKPMFuIcjDhPXMMY7s309mcRFuWi2ju0dRbDLh5+7On/peRfSk\nmwgZneA4xrkFRvIKywn0NRDg4+Y41uUV5exPO0BK4gmyzhjQafVo3N1RVCpsdjvGcgu2khL88/bT\nw5qBxg5mrNg8DGT3GEahrRNWS+Ux3bT1cz794k3WrUniuglxtX4moaGdiOnXF5OtlN3bd5GYuIHB\ngwb/9tn+epItyy8g+aSF9CIdCgoGvZpyqwlUNiJivRkff1WNxHLh96m8JI+irKPodRpKyko4sn0z\nuQeOcrbgPHa7jb8NGYhVo+JIuCcnOvkwOHI4tw8dV+37d6nfzwvtSU7jp6RTdb4+YGiXZktUdV00\n+Y26mtJOMVKLuwy0+yT1e350f//HM/ywbQduHj4EBgTQo2snOoQEERAQgL+/PwEBAY7/3Oq4cps9\nezbvvvsuvr6+aHTuVKBDq/dA7+aJ3uCBm4cnvbp3Iq5PV26++WZ8fHzqjclisbB8+XLmz59f5834\nuqkAW42lbm5ejP3jHNSdYgHw8dTh7VH9Rv74AWF0WLes2skA4GReLnd++jGF5eV17rVzeAzXj3uA\nG24fzXXDIsnakMibL73M6YJ8ckuLyS0todhkoqjcRF5ZGSaLtc5tXcxdq+X1G8cT6u2Ft6cf5Tfd\nR+K+c+QWGCkqNWO3g6KAl7sOg16NT3g2Hc7sIuJEAan64dhQgwJqlQa7wZNCuxad2YjBXMqx1D0s\n/vp13LUafPV63N180HmFYPDywcPLF7vNRklJAaczjmK328k/f67eWAcNHkf3kT0JGh5IlFsMAzsM\npGP6QQq2/MD5/GJSrcFku3fGDpRrNZTqLKhUKjQaNSoV0LmQUSNjiI8ait1u57Nvf+Zfr71G2vF9\nlBTmUlqUj9ViqnXfCnDhiSC8UwjDn/8LvuYeXNt9OIUlZg6m5JJz3oiXuw6NWsFmhwqLlehwPwbH\nhjp+O7m5uTz44IP4+fkRGhpKWFiY4/8B/kF8u+oENqsdu82G1WjEbrWiqNWo3dxQVCo0WjW33hVX\nZ9Jo7G/24oumKmmEkEYIKh8/tN5e2Ox2TBYb/uG+dIvtcEmJV7Scdv1JbEg+TeKuNMwXNDms3nzC\nUTNoaN0iz0EcOLCYwryzdZbTarXMmjWLZ555Bl9f3xqvv/zyy+zcuZOdO3fWuY0zqX0YNGRBgwkK\nKnuTBQYGMnXqVFau/IzU1JQG1/lNzQQVFNiFmyb9HW/vYFLsduyKQlGpGU83LSqV4ih3cEMSgeUm\nFKVymd1u59T5fJJOpeLt7l1rktLr3Em45i/07zMWRVGxZVsqI2KCyNm0mb5hHVmUnER2SUmN9RQg\nIDKU3FOZ1c+stVAUWLTzJ3qFBNEvrCOl33xPrnskhSW/JVOL1UZuoRFNUAZBRw8SnlpIka4TVtQo\nAHawWC1QUohO647eXIZitxHduTdjusey/vh+SioqoKQEcs406kjXjFMh80wK7j/6oPMNZGennRi/\n30Ov01loNSrKLAq5vh2x2WykZZ1g7/EkrBoNo8b+GYvVhmIF1WkvvjuyBbPJxEf/t4aPl76D1VL/\n4OIqFx5GjUrhicFxbDNZyVEd4/2vvVBsGkeZ/KJy9Fp15Xp2yM43cuBkLoG+bo7fTkJCAnfeeWet\n+1KrNXh7+OHj4cPIftfRL3owGdmnsCt2unbvB95epB7PpUdMzc4kjf3NWo1Gsjd+X23dvNJSdheq\n0IeFAmArLsKo0lJUZsFut5ObX8aBk7ms9jU06hwgWke7TVIbkk+zLim1xnJzhdWxvLYvqclk4r8f\nruPzLxPJPHUIs6mszn388Y9/5KWXXqp2b8FisbBv3z62bt3Ktm3b2LZtG+fO1X6F7eUXwsgJ99Kt\n90hSCzWUmywYLrrCKy8v54MPPmDHjh3s2LGj3qY+AJVaS6+4ceRlpnI29UC9ZQFyck+xcdNi/jhx\nLl52KFIqT0xGkwUPN+1v2y0vpay8Ag+3yhrWG9u38sGu5Lo2S/dugxg/Zibe3oGOZdYKK/vXb8fN\nZEanUaPX1Px6+hoMPDxqKOXX9OKgpoIDa7aTuq3u91FqrmBX+hl2pZ9hCXth7dfoDF54+HbAt0NX\nOnSLwyekO1p3A3qfdHodrRxoa1HpADuVDZOVSsoKSc/aS8a5k5w+d5L07FMYTaW17lelUtOr1xBO\nnthHmbG4zviq2O120tOPkp5+lM2bP0Hv4cYeT0+u8A8jwtsfi1sABw5/xv7jyZwvziXItwP33/wM\nZSiVzbWKgq0CTGcNbNXuwr3zUK6epGXf9jXknj3Z4P4vZLHZue+ztSifKxj8PAmK7EtIp5H4hHTF\n4OkHKJSbrahUClp1ZQfhwhIT9vJ83lj0A8vezaHgXAoqlQqbreaFj9VqITSgI+OHTsHD4MlH69/h\nZMZh5k5/hYrCQgDKSmsm16b8Zgv2H3DU6ktMJpbv3sVnB/bx8K0voq96nxYbZYUl2DV6x7Z0pRWU\nuGnqPQeI1tUuk5TRZCFxV1q9ZRJ3pTGib0cMeg379+9nyZIl7Nixg59//rnBqW86RPRk+ftvMzrh\n6mrLP/74Y+6++25KaqkdXEijMzAo4Vb6j5qCRlf5AzJXWNl3PIe4mNDqZTUaHnvsMUpLfztZenl5\nERcXR3DElXz37deOk1SHiJ4EhnUjYfJsANKP72HzmrfIyzxBl859KCw8x/mCLEDBzc0L7HaM5cX4\n+nYA7FX1CgCstupVGJPWrdqyPqFhAET6B9Av8kq++2UfJcYiPNy8GTvmb1zZc4Sj1lXFrlIwni/A\nDQjx9KLootrXFcFBPDJqGP7ubhw0WfEM82PI3X9gQMLNHF31KfsP7KtxLIM8PSgylmO6YOYHc3kx\n5qxizmcdJ3XvtwC4+fhz9wM3obHYQQGNrbJZzGa38cn6RZzOPM754rzaP7CL9AiPJfN8BgcObKuj\nxMWNazWZSo0cLjVy+FxOjddUiorbxs9CrzVgsdgwaVXY7XYURSH7VA6rH30dY2HtybMp7HY7xvxi\n0vK3krZ7KwAB4TEMnjyf0oIsirJTKMlNpeDcSQrOncTciGR8Vb9BDIu+ji4durN+5yo2/vQVFmsF\n00b/Fa2m8gLHUlyE4aJhYU39zVqKirDZ7Xyy92feTd5BgdHItEFX4+kZ4Chvs9lRXXS3Q7Hb0Zqs\nVLhpqm1POE+7PPr7judUay6ozYVJITU1lVdffdXxmqKoCAyLIrTzlQSERvH9568B4OUbzLDr/0qP\nfgl4hfSosc3w8HBKSkpQFIWYmBhGjBjB8OHDKbT6ce+fxwHQ86oxDL/+Ljx9g2qsX1TL1aVGo2Hq\n1KlA5cDPIUOGcMUVV6BWq9lx4CyrP11GePf+jJhwL8Edu2Gz/fa+w7v348aZ/+H07g1kn9hLZOc+\n/PTz10y84Qk6hvXEbrdzviATtUoDKFx4xNSq6gkm27cT6uIjjr/jOnfmqxl30dHHlzKLjTV7kujf\nYxjxY2fi5laz2dKuKFTo1bj5VTaJ6jQaXrp+PIt27uDnM5nc2OsKbunfG82vsy6U69WOdf2DOzPi\nrtncsWkln+89QHZJCaFeXmxNPc2f+sZybY9ulJjM7MvK4YtyNVmpGZSez8RcXvlZ2KwW7HYrPspv\n23Q3n0Nx64lKUZOVl+FIUHqtgfDgSDqHRNE5uDNjSOdvm9ZTaDbTreMV/GHoFAJ9gnnm/QcB6BDS\nmaxzpyuPmVrDiFF/Jq7PeH7ev44t21Zg+rUmPmDgdURG9SZ5x1eUlJzHaCzAaDTWOE4AYUERdAio\nvAeruuAcawcMAZ6YSirXM3h4ExLek5DwHoSE98DslsPe79dw7pfKeALc3ThvNGKz15029d7umEtM\nKCoNNosJg2cgZYXZbP5gZq2x6d28CO7UnZHD48g9m+KYeaNfv34sWLCA3h4+rFp9Eisqenbpzfrk\n1QT5dmBQr5GObahsVvzKMoHfOk809Ter8fZGpShE+vtTYDTia3BjwoAEqlK+zWbHDtguulACUH69\n2KrrwlC0rnaZpGo72ddXbsiQIdxwww0MGTIEjU8kp4p80OkrO0KcOvIjOr07AxJuof/IyY6aT237\nGDBgAF999RVDhw7Fz8/Psfz/3lpOSERPrr7xfkK7XFlnPBd3Vqjy/vvv17o8PEDN+Fvn0Cl6oKPW\nolKpq5Xx83bDb9QNaLuP5NixHfx1+mvo9JVJRFEU/P0qa0R2oPjX37OiUOPGstrNjc7jxpC7YQMA\nBo2Wjj6VCafYWMpTN0xH320cJk3t3ZrNHlp0eg29xw4kZc9WrCYzg8IjWfxjMk/Gj2B4ZGfMFhtg\nx6JSkRHsXhmLXY3eHALYKesfydwAP06fLyDC14dpfWM4npsPVHagcB9xJUO011NQWHmys1itWCwW\nSs+fxaw6hVH3W81NgwU/Uyr5hm5cfdV4bHYr4SFRBPuEovm1RhlWdISSM2UEGAw8MWAwhj63YVdV\nHpdpNz5AUEQ33NQ6Xvy/uwgLjeaG62fj3bkruhIzgwZNILbfMJJ3rqO4OA+VWkPXbn3p2rUPiqKg\n1h0n6tAeUrJLOFdkZGhEND8qHUj8aS3D+45xxGm74BxrV2yY/EoZ88TN+OsHEOzZt1pt1Wg4g29/\nHd+/+DGac4UM7RLBTbE9+cf67/nnmFG8u/Nnfkj9NaEqCjpfT8YtvI/yw31R7FrKCrPBbsfdJwSN\n3h21Ro9PcBT+Hbrh26ErXbpeQUhoRxRFYfzQSOY9cis9evTg2Wef5aabbkKlUnFuQyIRnCOVUCLD\nonn01ufJOZ+F+oLvZQTnUIzVE0NTf7O+vWM5u2YtgyO68OakyZzKz8dPr3EkqaqEbNHU/E3ZL7gA\na+x+Rctpl0mqrpN9XeWCg4NZu3YtADsPZnJ2w1FHGb2bJ3fMXY6Hl3+t617IYDBw/fXX11g+afzV\nnLNEUGGtuwlIp1U7xik1VqewEO75y7Ra2/GrXPtrm3vihmNc2XM4UNmUd2HTnVqlUKxXY//1Stbb\nQ1et0wTA6IERdIrrjFarqtHlNywwgFnTriExx58j+zNRLmhisSsKZg8tZg8t4wdG4OHjSXDCNWSu\nW4+iqHh69BhUKosjDqvNzpEuPlg0lTUq97JIVPbKr7HvkDGc8vqB8CMqFIuNLv5+dPH3w6pRkdYz\ngLOevfEoNFColP7aq0+FSq3BKzACVGGkB+xigKYAjbWyauFvqux0MqBHHDZF/euZTQ22CsKKjtOh\n5ARFbu4sih+NWlHIKk3hrFc0dpVC99g4FEVFdvoJJl1/N716/QGzlwGzmwZ9aQUqmwa93oORoyZj\nsVpQqX5NfYqCTWXjbMdy+pzyYFBYCOZgG3Y79Pdwo/OEx1FUv97nURTMmt8mjSnwz8auttEjrhe6\ntL5YKqp/RnpTMGp3HaP/8Wd6pBTgs+cUbjodz49LQKNS8cDIwcSEhfDfpF0MHzUAn6nDseZ0RIUe\nFPDw/a0jQ/ydi9AZPFGrFDRqFYoCwYEejqTorldx9913c8stt1QbYqHx9iaCynuwaYTQKbgLHYMq\nv4NqbERwjgjOofGqPmNHU3+zajc3x/doaJdI4iI6Y6eAk3TCigoFMGvdsHNRk/OvNfqm7le0HPU/\n//nPfzo7iEtVVFTEsmXLuOOOO2qdhqYuAT5ubN17psZ9lQvptGpuHtMDjab6zFEXr+vlG+SoVTW0\nbl38fH0AhePpBXWWGTu4Cz06N332ga6dfFGrFE5nFVd7vzqtmrGDuzAmrjNdO/mCXkNqZiGK2YYj\n/9hBpVaweuux++jRalS4GTR4uWtr3Q6AZ1QUgcOHog8MxKNzOH79+xM+bQpe0dH0jg3F7KYhNaeE\nCo0Ki0FDubcOtaeuxjYUtQpjWjpeWgMKCmZbBTatihPRARyK8Eaxq/Eo64qHMcoRw59GxlES7EVy\nkJESDw3F/m5kd/bh9OAIBgwdRw/P3pzIqDzGJrMVRals9ql8ryrsagXFo5AOBSaqKiBu1gICrRlo\nbWXoK8oo9/RErTURkXcEtc2CXq1CpVR2rvCynMemUcjzDcai1eCt9yQ8JIore/YnJCqAdHMFVRvW\nmCsTjV2xVk4cq+A4XRYEnMPoVYbeFExQQeW9JZsd3CqKADtmnTd2RcGo11ChUbCrbJz3P0dJUC6K\nojC5TwLd/LvU+D4pVJ58K3T5WMJ96ODnjc95E8oFF0fhgQH4DIzmhNqOn+9QrLmdfl23MvSq46L+\ntQai1ahQUPDx1Dnu3ei0am4ZewVXXdWvxqS4+gB/8rYn4W0tIowc3DHhrRgJ4TzRpOFHCWq9jvBp\nU1BdkNwu5Td74fcImw3Vr/WnQpU3Wm9viqiZgMyeOqy/Jqmm/o5Fy2i346Tq6ilUZfzQyDp79vye\ndRuK6eLutRcOlv09ymsZW1Kjp6DJwu5fsshIzUeDQmS4D0ZNwzNONPXGcmNiAbCWl1Ow7wCW4iLs\nHm6kh2gppgKDyg1bkT9GI7WuX24xcfDcEYpMpXjrPYgJ6Ynh1x5cVce4apyUxWrDarM7agRuoWeJ\nLthPj1MFaGw21CoVGrUKLw9fKq4cznpzKOYKK10yf6FL5iH0liIUxYKiVjC5a0i7IohzV4QRq7mK\nbp5R1QaKXvj56kor0JVWgMoEGjNgp8JupdA/m2L/fNzLIvEwRlXuJ+sXNDYrpgordjuYVDqOhHUh\nM0yNVWOmxKsQlcaORtEwOGwg942eUO29Xvx9Cu95nmzlGGZrBeoKK+4n8rEVmim3e5Gi70iFzQ17\nsT8Wi7qyc4ECWo0aq63yWFWdMVQqBZ1GhbdH9XFzDX3/6xq/VCV0/NhaZwe51N/dhd8jjZc3pyr8\n2b83i/wCI4UllR1kLqzRN/Z9iNbRbpMU/L6k0FIJpbEncHHpqo5xXmHljBPZ542czChAq1GhUinY\nFAvoMrlKU0gPP086d4omqF9/1AZDtc/HS2unizELS0khZ2yFlHYJwdvbt1pSrGvfRaVm3HVqvBUF\no9FETsU5UspyOJpWirosyNGEqdOqGd0nhP66Qn7++SQ/nirhqCaIEosKRW1B5ZWPm4cVb4MX42MH\nMH5I9zr3d+H36eJEnp3mzg+7syg3WzGWV2C12dHr1LjrNaSdK3HUOKtOFwE+BjRqleOYOWL99ftv\nt9v57LPP6NixI/3798dgqH4vsr7ps+qbvqq5fndVU0v9dCCTQ2nnMWoUqOV9COdr10kKfl9SkITS\ndrjKZ9lQHBcn2AAfNwJ8DM0Sb137Ligu58utKb/uz8ANI6Lw9TI0GOuJEycYOXIkeXl5DBgwgKFD\nhzr+CwkJqVHD8e0Ti9rQ8HyBzf1ZucpnL2rX7pOUEKLl/PLLL4waNYrc3Nxqy7t27epIWMOGDSMm\nJqbGuDkhQJ4nJYRoQVdeeSUbNmyoMSXYyZMnWb58OXPmzGHfvpqDsIWoIklKCNGi+vbty7fffoun\np2eN1+69916mTZsmtShRJ0lSQogWFxcXx1dffVXjaQAvvPACvXv35uuvv+YyvvMgWpAkKSFEqxg1\nahSrV69Gp9Oh0+l44okn0Ol0HDlyhD/84Q+MHTuWM2cubRZ50XZJkhJCtJqxY8fy6aefEhgYyIsv\nvsiRI0eYMmUKUNnJorbH2Yj2zaX7WR47dozFixfj7e1NZGQkt956q7NDEkL8ThMnTnTMXRkZGcmn\nn37Ktm3bKC4uxsPDw8nRCVfj0klq8eLFzJ49m9DQUP76178yZcoUdDqZS0uIy93IkSOr/T18+HAn\nRSJcnUs39+Xl5dGhQ+Wklj4+Pg0+h0kIIUTb4tJJqkOHDmRlZQFQUFBQ7fEWQggh2j6Xbu6bMWMG\nr776Kt7e3lx77bUylkIIIdoZl05SXbt2ZeHChc4OQwghhJO4dHOfEEKI9k2SlBBCCJclSUoIIYTL\nkiQlXIrJZCIpKYl///vfFBYWOjscIYSTuXTHCdH2FRUVkZSUxNatW9m2bRs//vgjFouFNWvW4OPj\n4+zwhBBOJklKtKrMzExHQtq6dSv79+/HZrM5XlcUhY8//pjrr7/eiVEKIVyFJCnRqnJzc1m4cCG7\nd++u9fX//ve/TJs2rZWjEkK4KrknJVpVbGwsL7/8cq0Tib7yyivcddddTohKCOGqpCYlWk1aWhqP\nP/44K1asqPHaP/7xDx555BEnRCWcyWo0UrD/AJaiIjTe3vj2jkV90YMRRfsmSUq0uNLSUhYuXMjC\nhQspLy8H4Nprr8Vms5GYmMgDDzzAvHnznBylaG3nEjeSvfF7rCazY9nZNWsJTriGkNEJToxMuBJp\n7hMtxm638/HHH9OzZ0/mz59PeXk53bp148svv+Tbb7+ld+/eTJ8+nVdffVXmZWxnziVuJHPd+moJ\nCsBqMpO5bj3nEjc6KTLhaqQmJVrE7t27efDBB9m+fTsAXl5ePP3009x///3o9XoApk6dylVXXYVK\nJddK7YnVaCR74/f1lsne+D2Bw4ehNhhaKSrhquTsIJpVVlYWM2bMYODAgWzfvh1FUbjzzjs5fvw4\njz76qCNBAcTFxaHRyHVSe1Ow/0CNGtTFrCYzBfsOtFJEwpXJGUI0q2+++YYPPvgAqHza6uuvv07/\n/v2dHJVwJZaiosaVK25cOdG2SZISzeqOO+7gyy+/ZNq0aUydOlXuNYkaNN7ejSvn1bhyom2TJCWa\nlUqlYtWqVc4OQ7gw396xnF2ztt4mP7Veh2+f2FaMSrgquSclhGhVajc3ghOuqbdMcMI10mlCAFKT\nEkI4QdU4qIvHSan1OhknJaqRJCWEcIqQ0QkEDh9Gwb4DWIqL0Hh549snVmpQohpJUkIIp1EbDATE\nDXR2GMKFyT0pIYQQLkuSlBBCCJclSUoIIYTLkiQlhBDCZUmSEkII4bIkSQkhhHBZkqSEEEK4LElS\nQgghXJYkKSGEEC5LkpQQQgiXJUlKiDZi3bp1bN68Gbvd7uxQhGg2kqREu1ReXs7WrVtJS0tzdijN\nZvjw4dxxxx3ExMTw5ptvUtTIJ+AK4cokSYl2oaCggK+//po5c+YwfPhwfHx8WLp0KeHh4c4Ordl4\ne3vz/vvv88svv3DffffRsWNHZs2axcGDB50dmhCXTLFfxm0DGRkZJCQksHHjRjp16uTscIQLycjI\nYNu2bWzdupVt27Zx4MCBas1gU6dO5X//+x9qtdqJUbaMWbNm8dZbb1VbNnLkSGbNmsWkSZPQarVO\nikyIppMkJdqcQ4cOMWnSJI4fP17r6wMHDuSHH37Azc2tlSNrHSUlJfTp04eUlJQar4WGhrJ8+XIS\nEuShguLyIM19os3p1asX27Ztq7UpLy4ujvz8/DaRoMorytl1Zh+bUraz68w+yivKAfD09OSDDz5A\nUZRq5SOmoP9hAAAgAElEQVQjI1m3bp0kKHFZkSQl2hSbzcaSJUvo06cP6enp1V7r378/AFar1Rmh\nNatNKUks2PIfPj+0ju9ObOHzQ+tYsOU/bEpJAiqb9x588MFq66SmpvL2229jNptr26QQLkmSlGgz\nduzYQVxcHH/5y1/IysrC19eXkSNHApW1q8GDB5OcnIyfn5+TI/19NqUksfDfr/D90m85vfckFeWV\nScdsreC7Ez84EtWCBQuIjo4GKu/BASxatIjRo0eTnZ3tnOCFaCJJUuKyl5GRwW233cbQoUP56aef\nUKlUzJw5k+PHj3PrrbcSFRXFI4884uhM4O/v7+SIL115RTmbU5OIHdOPM4fT+eSJ93jtj/NZ9uDb\nbFq0juNJv7Bu7wbKLSbc3NxYunQpkZGRfPLJJ/znP/9BrVazdetWBg4cyN69e539doRokHScEJct\no9HIK6+8wosvvkhZWRkA8fHxvPbaa8TGxgKwZcsW7HY7EydOpLCwEIDJkyezcuVKp8X9e+w6s4/P\nD60DoKLczKdPLSHj4Kka5SK7RzH66gRGjRpFREQEI0aMAGDTpk1MmTKF/Px83N3dWbJkCVOmTGnN\ntyBEk0hNSlx27HY7K1eu5IorruDpp5+mrKyMqKgoVq9eTWJioiNBAQwaNIiHH37YkaCAy7q5r9hU\n4vi31qBj8vzb6dC9Y41yqcdT2LJlC5GRkY4EBZVJfNeuXfTq1YuysjKmTp3K008/jc1ma5X4hWgq\nSVLisrNp0yamTp3K6dOn8fDw4IUXXuDQoUPceOONNXq0zZ49m59//rnassu5uc9L71ntb72HgSnP\nTycgIrhG2b59+9KlS5cay6OiotixYwcTJ04E4Nlnn+Wmm26iuLi4RWIW4veQJCUuO/Hx8cTHxzN9\n+nSOHz/Ok08+icFgqFHuf//7H++8806N5ZdzTSo2uAc6dfXBuO4+Hkx7cQa+odWT74oVK4iOjuaF\nF16gvLy82mteXl6sWrWKv//97wCsWbOGhIQEqVEJlyNJSlx2FEXhm2++4YMPPiA0NLTOcgMHDiQn\nJ8dxz8Xd3R24vGtSBq2BqyOH1ljuFeDNtBdn4BngjUarYeXKlXTu3JnS0lLmzp3LlVdeyerVq6vN\nuqFSqXj22WdZsWIF7u7uzJ49G5VKTgnCtcg3UlyWdDpdg2W6d++OWq3miy++AGDZsmVMmDDhsq5J\nAcRHDeXabqNq1KiCO4bw+kdvEd09msmTJ3P48GGee+453N3dSU1N5Y9//CMJCQns37+/2npTp07l\n5MmT/OlPf2rNtyFEo0jvPtGmvfXWW8yaNYugoCDOnDmD1WolJyenTUwsW24xcfDcEYpMpXjrPYgJ\n6YlBo+f06dN07tzZUe7MmTPMmTOH5cuXA5U1qHvuuYf58+cTGBjorPCFaBSpSYk2bcmSJQDcdttt\naLVaDAZDm0hQAAaNngEd+xAfNZQBHftg0OgBqiUogI4dO7Js2TJ27NjBoEGDsNlsvP3223Tv3p1/\n//vf8vwp4dIkSYk269ChQ+zatQuA6dOnOzcYFzB48GB27NjB0qVL6dChAwUFBezcubNGj8iG1DVn\noBAtQePsAIRoKUuXLgWgX79+9O7d28nRuAaVSsXtt9/OpEmTePHFF/nb3/7WpPU3pSSxOTUJs7XC\nsezLIxu4OnIo8VE1O3QI8XtJkhJtksVicdyDkVpUTV5eXjz//PNNWmdTShLfnfihxvKqOQMBSVSi\n2Ulzn2iT1q9fT1ZWFlqtlltuucXZ4Vz2quYMrM/m1CTKLaZWiki0F5KkRJtU1WHihhtukB5szeBA\n9tFqTXy1MVsrOHjuSCtFJNoLSVKizcnLy2Pt2rWANPU1lwvnDKxPkam0hSMR7Y0kKdHm/PjjjwAE\nBwdz3XXXOTmatuHiOQPr4q33aOFIRHsjSUq0OePGjePs2bOsWrUKrVbb8AqiQbXNGXgxnVpLTEjP\nVopItBeSpESbFBAQwLBhw5wdRouzGo3kJf/IuQ2J5CX/iNVobJH91DVn4IWujhzqGFAsRHNplS7o\nH374IceOHaO0tJTJkycTFRXFK6+8QmBgIG5ubjzwwAMsWbKEjIwMiouLueeee/Dw8KhRRgjxm3OJ\nG8ne+D1Wk9mx7OyatQQnXEPI6IRm319V9/KLx0np1FoZJyVaTKskqW7dunHbbbdx4sQJVq5cyY8/\n/si0adMYMGAATz75JJmZmWzfvp13332X9PR03nnnHYKDg2uUqW/GayHak3OJG8lct77GcqvJ7Fje\nUolqaMRVtc4ZKERLaJEktWTJEnbs2OH4e/bs2eTn5/Pee+/xwAMP8NZbbxESEgJU3txOT093PD4h\nJCSEnJwcVCpVtTI5OTmSpISgsokve+P39ZbJ3vg9gcOHoa7lOVu/V9WcgUK0hhZJUtOnT6/W9ffI\nkSO8+uqrPPnkk/j4+BAaGsq5c+cIDw8nMzOTqKgozp8/D8DZs2fp2LEjQUFB1cqEhYW1RKhCXHYK\n9h+o1sRXG6vJTMG+AwTEDWylqIRoGa3S3Pf4448TFxfHO++8Q1RUFFOnTmXhwoUkJibSpUsXAgMD\nGTFiBM899xzFxcXMnDkTDw+PGmWEEGApKmpcueLGlRPClcnzpIS4zOQl/0j6is8aLBc+bYrUpMRl\nT7qgC3GZ8e0di1pf/5OJ1Xodvn1iWykiIVqOJCkhLjNqNzeCE66pt0xwwjUt0mlCiNYmj+oQ4jJU\n1b384nFSar2uxcZJCeEMkqSEuEyFjE4gcPgwCvYdwFJchMbLG98+sVKDEm2KJCkhLmNqg0E6R4g2\nTe5JCSGEcFmSpIQQQrgsSVJCCCFcliQpIYQQLks6TgghXIbdbmfp0qXo9XoSEhIIDg52dkjCyaQm\nJYRwGYqiMGHCBJ599llCQkLo06cPjzzyCN988w0lJSXODk84gczdJ4RoEqvRSMH+A1iKitB4e1dO\n0+Tm1qz7SEtLY8iQIZw9e9axTKvVMmTIEEaPHs3o0aMZOHAgGo00BrV1kqSEEI1W29OAW2qWi/37\n9zNixAiKapn1Xa1WM2/ePObOnYuiKM26X+FapLlPCNEoVU8DvvhZVlVPAz6XuLFZ99e7d2/WrFmD\nVqut8dq7777LU089JQmqHZAkJdoUq9FIXvKPnNuQSF7yj1iNRmeH1CY09mnA1vLyZt3vNddcw7Jl\ny2osv/POO7nvvvsoLCxs1v0J1yNJSrQZ5xI38sv850lf8RmZ33xH+orP+GX+881+hd8eNeVpwM3t\n5ptv5l//+hcAERERDBs2DLvdzptvvknPnj35+OOPuYzvWogGSJISLsNssnD0YBZ7ktM4ejALs8nS\n6HVbuymqvXH204Bnz57Nww8/TM+ePdmyZQuLFy/G39+frKwsbrnlFq699lqOHTvWIvsWziVJSriE\nPclpfPRuMls2HOOnpFNs2XCMj95NZk9yWoPrOqspqj3ReHs3rpxX48pdipdffpmZM2eiUqm48847\nOXr0KDNmzAAgMTGR2NhYnnnmGcrlc25TJEkJp9uTnMZPSaewVFirLbdUWPkp6VSDicqZTVHthSs8\nDVilUjFx4kTH34GBgbz33nts3bqVXr16YTabmT9/PjExMaxfv77F4hCtS5KUcCqzycLeXen1ltm7\nK73epj9nN0W1B678NODhw4ezZ88eFi5ciLu7OydPnmT16tWtHodoGZKkhFOlHs+tUYO6mKXCSurx\n3Dpfd4WmqPYgZHQCoePH1qhRqfU6QsePderTgLVaLY899hiHDx/m9ttvZ8GCBU6LRTQvGa4tnKqs\ntP5musaU8+0dy9k1a+tt8mvppqj2wtWfBhwREcHSpUudHYZoRpKkhFO5e9R/n6Mx5aqaojLX1X0f\nwllNUW2RPA1YtCZp7hNOFdk9EI1WXW8ZjVZNZPfAesu4clOUEOLSSU1KOJVOr6HvwHB+SjpVZ5m+\nA8PR6Rv+qrp6U5QrMZsspB7PpazUjLuHjsjugY06xkK0NvlWCqfrFxcBVPbiu7AThUarpu/AcMfr\njSFNUQ3bk5xW41gnbT7Z5GMtRGuQJCVcQr+4CHr1DZOr+xZWNSbtYlVj0gBJVMKlyBlAuAydXkOP\nmA7ODqPNauyYtF59w+TiQLgM6TghRDvRHGPShGhtkqSEaCeaY0yaEK1NkpQQ7URzjEkTorVJkhKi\nnWiuMWlCtCZJUkJchhITE5u8TtWYtPo0dkyaEK1FkpQQl6F3332XL774osnr9YuLYMDQLjVqVBqt\nmgFDu0j3c+Fy5JJJiMtQYGAgt99+O7t27SI6OrpJ68qYNHE5kW+lEJehoKAgioqK+OMf/8jOnTvx\n9PRs0voyJk1cLqS5TwgXsHfvXnbt2tXo8kFBQQAcOnSIv/71r9jt9pYKTQinkiQlhJNYLBY+//xz\nRo0axdSpU4mJiWn0ulVJCmDFihW8/vrrLRGiEE4nzX1CtLK8vDwWL17Mm2++SXp6OoqisGXLFtzc\n3Bq9jeDg4Gp/P/roo/Tv35+RI0c2d7hCOJUkKSFayYEDB/j3v//Nhx9+SHl5uWP5fffdx/Dhw5u0\nrQtrUgBWq5WpU6fy888/ExYW1izxCuEKpLlPiBaWmJhIfHw8vXv3ZvHixdUSVJcuXViwYEGTt3lx\nklIUhfz8fCZPnozZLNMaibZDalJCtLBBgwYxceJEDhw4QG5u9clbFy9e3OSeeQABAQEoioK7uzul\npaWoVCpOnTqFwWDAZrM1V+hCOJ3UpIRoYV5eXgDk5+dXW3733XeTkFD7Y+3tdjvHjh2rc5tqtZqu\nXbuyfft2IiIisFqtLFu2DH9/fwzyJGLRhkiSEqIFVVRUcO+99/LQQw9hs9kcPfg6derEwoUL61zv\n/Pnz/P3vf693219//TV9+vRhxowZQGWtTLqitx9bt27l4Ycfpqys7Hdv64UXXgBg2bJlv3tbzU2a\n+4RoIefPn2fKlCls3LgRgHvuuYfXXnuN4OBg/vvf/+Lj41PnuikpKXz22WccPnyYK664otYyVTNN\n/OUvf2HevHmcPHmSzZs3c8011zT/mxEuZ+3atRQUFPDss88SFhaGSqVi4sSJzJkzhz/84Q/s2bOH\n2NhY9u7dS3R0NCaTCV9fXxISEnjrrbfw9/cnODiY2NhYtm3bxqFDh9i5cydTp07lueeeIzQ0lIKC\nAp566ikmTZrEtGnTSEpK4qWXXmpST9TfS2pSQrSA48ePM3jwYDZu3IhKpeK1117j7bffxmAw8Pbb\nbzN+/Ph61z958iR2u50XX3yxwX1FRERw3XXXAZW1KdE+DBkyhP79+3PjjTdy//33c/ToUaxWKxER\nEUybNg2AsWPHMmbMGLy8vJg5cybJyckYDAZ8fHxwd3dn27Zt9O7dm86dO9OrVy8Atm3bRv/+/Zk1\naxYAZ86cwd/fn5tvvpmuXbuSlpbWqu9TkpQQzWzTpk3ExcVx7NgxvLy8+Oqrr3jwwQdRFAWAW2+9\ntcFtpKSkAPDRRx85/l2fv/71rwB8/vnnNe59ibbtwibeqs40VXQ6HSqVyvF/q9XKmjVrGDlyJHfc\ncQcWi6VR29Tr9QCoVKpW75gjSUqIZrRo0SLGjh3L+fPn6dKlCzt27GDcuHFN3k5VYrJarbz00ksN\nlr/hhhsIDg7GZDLx4YcfNnl/4vIUFRXFV199xRtvvEHv3r1RqRo+pV955ZWsWLGCZcuWoSgKKSkp\nmM1mx7RcI0aMYN++fbz11lsYDAbnj7uzX8bS09Pt0dHR9vT0dGeHIto5i8Vif+ihh+yAHbAPGzbM\nnp2dfcnbi4+Pd2xLp9M16jv++OOP2wF7TEyM3WazXfK+hXAlUpMSohmcPn2a9957D4A///nPbNy4\nscaA26a4sInPbDbzyiuvNLhOVZPfwYMH+fHHHy9530K4EklSQjSDqKgoVqxYwYIFC1i6dKmjDf9S\nmM3mGjenFy1aRHZ2dr3rde/enVGjRgHSgUK0HZKkhGgm48aNY86cOY4OEpcqLS2t2s1pjUaD0Wjk\n1VdfbXDdu+66C4Dc3FwZM3UZMZos7DyYyXfJp9l5MBOjqfYODe2RjJMSwsWkpKQQHBzM2LFjWb58\nOT169ODLL79k06ZN2Gy2em+O33TTTQwZMoSoqKhWjFj8HhuST5O4Kw1zhdWxbPXmE4weGMGYuM5O\njMw1KPbL+HIrIyODhIQENm7cSKdOnZwdjhDNYu/evYSGhnL8+HFGjBiBVqulrKwMjUauKduaDcmn\nWZeUWufr44dGtvtE1eRvfWlpKR4eHi0RixAC6Nu3L4AjKVVUVJCSkuKYYUK0DUaThcRd9Q+MTdyV\nxoi+HTHoGz5Vr1q1iq+//tpRiw4LC+Mvf/lLrWUzMjJ4++23ef7552u8tmjRIjIyMtBoNOTn5/O3\nv/2Nnj17NuId/WblypVce+21/N///R/z589v0roXq/ed33jjjTz33HPVnhj6yCOP8M477/yunQoh\nGhYQEEBQUBA5OTkcPnxYklQbs+94TrUmvtqYK6zsO55DXExoo7Y5YcIEJk6cCIDRaGTu3LmEhYWR\nmZnJ/PnzmTVrFt26dWPKlClA5fn88ccfJyQkhLvuuotbbrkFk8nkSCylpaWUl5eTkZHBo48+ytVX\nX018fDzLli0jKCgIrVbLzJkz+cMf/sCtt97Krl27uP/++9m9ezcjRoxg586drF27lvDwcDZu3Iib\nmxsajYZ77rmn0cep3iTl5eXF+++/z5AhQxxv6lKdPn2am2++mdWrV6MoCq+88gqBgYG4ubnxwAMP\nsGTJEjIyMiguLuaee+7Bw8OjRhkh2psrrrjCkaSqTj6ibSgqbdxzvxpbDirn8zt48CAA/fr1Iyws\nDHd3d06ePEl2djalpaXMmjWLvLw8oLIi8umnnzJmzBi6d+/OyZMn6d+/PwDr1q1jz5492O12pk+f\nTkhICH/729/IzMwkICAAHx8fvvnmG2bOnImnpyd/+tOfcHNzY9++fY54wsLCmDBhAg888ABdu3bF\nZrPxyy+/NPr9QAO9+zw8PPjXv/5Fbm4u//znP6moqGjSxquUl5fz7rvvMnToUAA++eQTpk2bxhNP\nPMHZs2fJzMxk+/bt/P3vf+e+++7jvffeq7WMEO1N1eSyhw8fdnIkorl5e+iatRxU1qSeeuopnnrq\nKXQ6HV5eXsyYMQM/Pz+sVitqtbra5LDDhg1j+/btfPbZZ9x0001ceeWVJCcnAzB+/HiefPJJMjIy\nABzPPXv//fcZN24ct956qyMnVD0eRlGUOqdNuvXWW7n//vuZM2dOo98PNPKe1L333suWLVu49957\nKSwsbLD8kiVL2LFjx2870WiYO3cub7zxBlDZPTYkJASA4OBg0tPT8ff3ByAkJIScnBxUKlW1Mjk5\nOYSGNq7KK0RbIUmq7erTPYjVm0/U2+Sn06rp073xg8IvrEmdPXsWs9lMWVkZQUFBrF+/vkZ5lUrF\nsGHDSEpKomvXrnTt2pXDhw/z9NNPo9PpKCgoqDHXZN++fXnvvfeIiooiJCSEn376qc54oqOjWbx4\nMTNmzODFF18kICAAf3//JjX31du7b9OmTcTHxzv+TktLY8GCBU26J3X27FleffVVIiIiSExMZPTo\n0ajVagYNGsSAAQN47LHHeOKJJ5g7dy6LFi3i1KlTLF26lKCgoBplAgMDq21beveJtu67775j7Nix\neHp6UlRU9LvHYAnX4gq9+zZv3kx2djZTp05t0f1cqlbtgv7kk0/y0EMPodFoWLhwIf7+/nh5eTFr\n1iyWL1/O6dOnKS4uZubMmXh4eNQoczFJUqKtS09PJyIiAqi8SAwPD3dyRKK51TZOSqdVt8o4qQ0b\nNpCYmMizzz6LTtf4ZsXWJOOkhHBhdrsdb29vSkpKWL9+Pddee62zQxItoNxkYd/xHIpKzXh76OjT\nPahR3c7bA5kWSQgXpiiKY4yK3Jdquwx6DXExoYyJ60xcTKgkqAvIkRDCxT3xxBOYTCaGDRvm7FBE\nCymvKOdA9lGKTSV46T2JDe6BQWtwdlguQZKUEC5u8uTJzg5BtKBNKUlsTk3CbP1tiM+XRzZwdeRQ\n4qOGOjEy1yBJSgghnGRTShLfnfihxnKztcKxvL0nKklSQgjhBOUV5WxOTaq3zObUJIZGXIVBU//z\nybKzs3nxxRcJCwujuLiYjh07EhAQwE033dToeJKTk9m9ezczZ86s8drYsWN54oknHEOS5syZQ1hY\nGPfff3+Nsm+88QZlZWVERUVx7bXX4uPjU+31Dz/8kGPHjlFaWsrkyZMZMmRIvXFJkhJCCCc4kH20\nWhNfbczWCg6eO8KAjn3qLXfo0CFCQkJ4+OGHUalUpKSk8N5776EoCnv37iUkJISMjAxiY2PZvn07\nL7zwArNnz2bYsGEUFBQQHh7uGOpw6tQpPvnkE3x9fSksLOSJJ54gJiaG9evXEx8fT35+PiaTCahM\njq+++iodO3akuLiYGTNm8P333xMfH++Yv++LL74gIyOD7OxsHnvsMbp168Ztt93GiRMnWLlyZYNJ\nSnr3CSGEExSbShpVrshU2mCZq6++mujoaObNm8fTTz/Nzp07Ha/FxsYyc+ZMTp06xS233EJ0dDQn\nT57EbDYzdepUHn74Yb755htH+RUrVmCxWKioqODs2bMUFxej1WoJDw/n9OnTrFq1ikmTJgGV0yGF\nhITg4eHBtm3bCAkJITo6utp91M2bNzN37lzmzZuHp6cngwcPJj8/n/fee4/p06c3+N6kJiWEEE7g\npfdsVDlvfcOPRjp27Bjx8fFMmjQJq9XK0KFDGT16NAA6nQ5FUdDrK5sMVSoVVutvA4ftdnuNpzjf\ncMMN9OnTh6ysLLy8vACYOnUqy5Yto7i4mOuuu469e/eyatUqYmJiGD16NKtXr641tqptW61WKioq\nyMzM5KOPPuLJJ5+s0RRYG0lSQgjhBLHBPfjyyIZ6m/x0ai0xIQ0/y8lut/Pss88SFBSE2Wzm7rvv\nJiUlpd51VCoVn332Genp6UyYMMGx/Oabb+b111+nQ4cOWK1Wx4SwwcHB5OXlOZIfQK9evVi6dCnH\njh2jW7dufPvttzX2k5CQwIIFC8jLy+Ohhx7i8ccfJy4ujnfeeYeoqKgGn7AhM04IIYST1NW7r8q1\n3Ua1WO++6dOns2TJkhbZdnOSmpQQQjhJVQK6eJyUTq2VcVK/kiQlhBBOFB81lKERV3Hw3BGKTKV4\n6z2ICenZYLfz3+tyqEWBJCkhhHA6g0bfYDfz9kqSlBC/g8y5JpqD1WikYP8BLEVFaLy98e0di/qC\nJ+i2Z5KkhLhEMueaaA7nEjeSvfF7rCazY9nZNWsJTriGkNEJTozMNchgXiEuwaaUJP638VNMFnO1\n5VVzrm1KqX+6GyGgMkFlrltfLUEBWE1mMtet51ziRidF5jqkJiVEE1XNuZa8civD/5yAT4hfjTKN\nnXNNtF9Wo5Hsjd/XWyZ74/cEDh+G2lB/E3JLz903ceJEBg0ahNVqxWKxMH/+fN544w2GDBlCWloa\nX3/9NVFRUQAMHTqUa665ptH7bYgkKSGa6ED2UUpKSji69SBd+nerNUk1ds410X4V7D9QowZ1MavJ\nTMG+AwTEDay3XEvP3efn58dTTz0FwJ133onFYqm2/wkTJjBx4sTfcTTqJklKiCYqNpVwbNshKsrN\nnDl0ml7xfWst15g510T7ZSkqaly54obLXX311RQUFDBv3jzsdrvjac5QOXff5MmTueWWW3jhhRfI\nycmpNnefp6cnM2bM4J577gFqn7uvoKCAV155hfz8fEJDQ1EUpdr+165dy8GDBwEYM2YMgwYNauxh\naJAkKSGayEvvycHEPQBkHDpdZ7nGzLkm2i+Nt3fjynk1XK6l5+7z9fXl0UcfBWDx4sX88EP1WTKk\nJiWEC/E1uXN6X+W8aDmnzlFeYsTgWb27cGPnXBPtl2/vWM6uWVtvk59ar8O3T2yD22rpufsKCgp4\n6aWXAMjMzGTSpEkcOnSoke/095G5+4RoogULFjja5wEmP3cHXQf2qFamJedcE21HVe++uoSOH9ti\n3dBl7j4h2iC73c7SpUurLTtz6LQjSTVmzjWzyULq8VzKSs24e+iI7B6ITi8/xfaoKgFdPE5KrdfJ\nOKlfyS9DiCZITk7m2LFj1ZaZUou5ttuoRs25tic5jb270rFU/HZPIGnzSfoODKdfXESLxS1cV8jo\nBAKHD6Ng3wEsxUVovLzx7RPbYLfz3+tyqEWBJCkhmmTZsmU1lv2y9xDDOw1Ap9PVu+6e5DR+SjpV\nY7mlwupYLomqfbIqGnI9wimjsnbtpWhQOzsoFyEzTgjRSGazmby8PNasWUNAQAAAs2fPplOnTuzZ\ns6f+dU0W9u5Kr7fM3l3pmE2WesuItmdPchofvZvMlg3H+CnpFFs2HOOjd5PZk5zm7NBcgiQpIRpJ\nq9WyYsUKJkyYQNGvY1ymTJnCwYMH6dixY73rph7PrdbEdzGbzUZZaRmpx3ObNWbh2qpq1xd/N6pq\n15KoJEkJ0WhVAxiNRiMVFZWTyvr6+qLX6xvsXVpWWv/MAqfSDvPLsV0NlhNtR3PWrrOzs3n44Yd5\n5ZVXeOaZZ1i0aBGff/55k+JJTk7mrbfeqvW1sWPHsmnTJsffc+bM4Y033qi17BtvvMFLL73EypUr\nKSwsrPH6zz//zPTp0/niiy8aFZfckxKiiQoKChz/9vX1bdQ67h7136/a90sSxcXncfeY8btiE5eP\nhmrXUFmjSj2eS4+YDvWWa+lpkWJiYli/fj3x8fHk5+djMpmAyuT46quv0rFjR4qLi5kxYwbff/89\n8fHx7N69mxEjRvDFF1+QkZFBdnY2jz32GEFBQdXGZTVEalJCNNGlJKnI7oFotHXfCj9waAf7f9lB\n567+vzs+cXlobK25MeWuvvpqoqOjmTdvHk8//TQ7d+50vBYbG8vMmTM5deoUt9xyC9HR0dWmRXr4\n4dEPke8AAA7eSURBVIf55ptvHOVrmxZJq9USHh7O6dOnWbVqFZMmTQLAYDAQEhKCh4cH27ZtIyQk\nhOjoaCZPnuzY3ubNm5k7dy7z5s3D09OT8PDwxh4iQGpSQjRZVZLS6XQYGtlNWKfX0HdgeK29+/Ly\ns8jIPAnAgYP7GDiw/slERdvQUO26KeVaelokgKlTp7Js2TKKi4u57rrr2Lt3L6tWrSImJobRo0ez\nevXqWmOr2rbVaqWiogIfH59Gve8qkqSEaKKqJOXn51djos36VHUvv3ic1MGjv131fvnll5Kk2onI\n7oEkbT5Zb5OfRqsmsntgg9tq6WmRAIKDg8nLy3MkP4BevXqxdOlSjh07Rrdu3fj2229r7CchIYEF\nCxaQl5fHQw89xPLly/nhhx9QqVRUVFRUq3XVRqZFEqKJPvroI2677TZ69OjBkSNHmrz+xTNO3P/w\nHWzY8B0Affv2bbA7u2g76ho79//t3W1sVNW+x/EfnWkLDBTFaacPDNIeUQMC06SI4RzUY09AeWHi\nOcbGkEijBhQSjDGEBl4oCbmJBBWNBBQ1BU1I1JDjBa4aGelBS0OK0nMlvpCQGqaZ0jPTh3QoZYZp\n931x29rSB9rSzqzp/n5edXZXJ//V7ulvr/2wVq+SlQsm7dk5pkUCpqjekdRor0fdLCPT2Xch/OrV\nq/rXv6r6vldXV6dAIDDm8/ZITcONrp3pDmYh6UFIAWN0uyHV33fffadYbOCF8ePHj+vll1++7fdG\naiheMV+LffkJn88xFUZREiEFjNmzzz6rZcuWjfkC8FCOHTs25DZCyl76j64xECEFjFFRUZGKiopu\n+326u7t14sQJPfzwwzp9+rQkadWqVTp16pQ6OjrkcrFoIkBIAUnS0tKio0ePaunSpcrqWaX1448/\nVkdHh8LhMCEFiJACksbtdsvtdisSiQzY7vP5klQRYB5CCrgNXZ2davvfXxRvb5czK0t3LF0ix4wZ\nt/5Bm2GhR4wXewkwTk0n/YNWVA3+879ZUfUmLPSI20FIAUO41ZF/00m/Gv/n20E/1xWN9W0nqFjo\nEbePkAJucqsj/67OTv3Hf2rE9/iP/5Tcf/nzpC8BbrLRLkWx2JfPqT8Miz0D6Gc0R/7zdWXAKb6h\ndEVjavv3L7prxcB5+IYaoU1VE7kUBeyLkAJ6jPbIP/tP10b1fvFI+4DXw43QFi66c+zFpoCJXIoC\n9kVIAT1Ge+R/pTNzVB8c5+ysvq9HGqFN1SXCJ3IpCtgXix4CPUZ7RD8tJ1+OzJH/sToyM3THsiWS\nbj1CS0tzqMT3qP7+939o1qxZoy/YcLda6FEa/VIUsC9CCugx2iP6WXfMUk7pX0dsk1P6176bJm41\nQsvIyNSG53bqv3a+r7y8vNEXbLjehR5H4lvu5aYJjIiQAnqM5cjf87dS5a1dM2hE5cjMUN7aNQNu\nP7fztZniFfNVsnLBoN+rM90xqWslYergEAboMdIS7736H/l7/lYq91/+rLZ//6J4pF3O2Vm6Y9mS\nQbed2/3aTLKWosDUwF4C9DPWRegc06cPus38ZhO5THiqYikKjBchBdxkoo/8xzpCA/AHPhXAECb6\nyJ9lwoHxIaSABOHaDDB2fDqABOLaDDA2CQmpH374QT/++KOuX7+usrIy3XXXXdqzZ4/cbrdmzJih\nLVu2qLKyUg0NDYpEItq4caNcLtegNgAAe0lISH3xxRdasmSJmpub5Xa7deTIEZWVlamkpEQVFRVq\nbGxUdXW1Dh48qEAgoAMHDignJ2dQm6n0oCMA4NYmJaQqKytVU1PT97qqqkpvv/22rly5og8//FDR\naFQej0eSlJOTo0AgoLlz50qSPB6PQqGQ0tLSBrQJhUKEFADYzKSEVHl5ucrLy/tel5WVKS0tTXPm\nzFFnZ6cKCgrU1NQkr9erxsZGFRUVqbW1VZIUDAZVUFCg7OzsAW3y8/Mno1QAgMEScrpv/fr12rFj\nR9/Xbrdbu3fv1smTJ7VgwQK53W6tWrVKu3btUiQS0aZNm+RyuQa1AQDYyzTLsqxkFzFeDQ0NKi0t\nld/v17x585JdDgBggjHBLADAWIQUAMBYhBQAwFiEFADAWIQUAMBYhBQAwFiEFADAWIQUAMBYhBQA\nwFiEFADAWIQUAMBYhBQAwFiEFADAWAlZqgPA0GLRuOovhnWtI6aZrgwVLnQrI5OPJdCLTwOQJOfP\nXlZdbUDxG119285UXZJvuVfFK+YnsTLAHIQUkATnz17WuTO/D9oev9HVt52gArgmBSRcLBpXXW1g\nxDZ1tQHFovEEVQSYi5ACEqz+YnjAKb6hxG90qf5iOEEVAeYipIAEu9YRm9B2wFRGSAEJNtOVMaHt\ngKmMkAISrHChW850x4htnOkOFS50J6giwFyEFJBgGZlO+ZZ7R2zjW+7leSlA3IIOJEXv7eU3Pyfl\nTHfwnBTQDyEFJEnxivla7MtnxglgBHwagCTKyHTqvgdyk10GYCyuSQEAjEVIAQCMRUgBAIxFSAEA\njEVIAQCMRUgBAIxFSAEAjEVIAQCMxcO8gI3FonFmvIDR2BsBmzp/9vKguQPPVF1i7kAYhZACbOj8\n2cs6d+b3QdvjN7r6thNUMAHXpACbiUXjqqsNjNimrjagWDSeoIqA4RFSgM3UXwwPOMU3lPiNLtVf\nDCeoImB4hBRgM9c6YhPaDphMhBRgMzNdGRPaDphMhBRgM4UL3XKmO0Zs40x3qHChO0EVAcMjpACb\nych0yrfcO2Ib33Ivz0vBCOyFgA313l5+83NSznQHz0nBKIQUYFPFK+ZrsS+fGSdgNPZGwMYyMp26\n74HcZJcBDItrUgAAYxFSAABjEVIAAGMRUgAAYxFSAABjEVIAAGMRUgAAYxFSAABjEVIAAGMRUgAA\nYxFSAABjJWTuvsOHD6ulpUXhcFhPPvmk7r77bu3Zs0dut1szZszQli1bVFlZqYaGBkUiEW3cuFEu\nl2tQGwCAvSQkpGpqarR//379+uuvOnbsmKZPn66ysjKVlJSooqJCjY2Nqq6u1sGDBxUIBHTgwAHl\n5OQMapOXl5eIcgEAhpiUkKqsrFRNTU3f65UrV+q1115TKBTS1q1b9fnnn8vj8UiScnJyFAgENHfu\nXEmSx+NRKBRSWlragDahUIiQAgCbmZSQKi8vV3l5ed/rl156SQcOHFB7e7u2b9+uRYsWqampSV6v\nV42NjSoqKlJra6skKRgMqqCgQNnZ2QPa5OfnT0apAACDJeR03/333693331X7e3tWrt2rR588EHt\n3r1bJ0+e1IIFC+R2u7Vq1Srt2rVLkUhEmzZtksvlGtQGAGAv0yzLspJdxHg1NDSotLRUfr9f8+bN\nS3Y5AIAJxi3oAABjEVIAAGMRUgAAYxFSAABjEVIAAGMl5BZ0AGMTi8ZVfzGsax0xzXRlqHChWxmZ\nfFxhP+z1gGHOn72sutqA4je6+radqbok33KvilfMT2JlQOIRUoBBzp+9rHNnfh+0PX6jq287QQU7\n4ZoUYIhYNK662sCIbepqA4pF4wmqCEg+QgowRP3F8IBTfEOJ3+hS/cVwgioCko+QAgxxrSM2oe2A\nqYCQAgwx05Uxoe2AqYCQAgxRuNAtZ7pjxDbOdIcKF7IiAOyDkAIMkZHplG+5d8Q2vuVenpeCrbC3\nAwbpvb385ueknOkOnpOCLRFSgGGKV8zXYl8+M04AIqQAI2VkOnXfA7nJLgNIOq5JAQCMRUgBAIxF\nSAEAjEVIAQCMRUgBAIxFSAEAjEVIAQCMRUgBAIxFSAEAjJXSM050df3/3GZXrlxJciUA8Ifc3Fw5\nnSn979UYKf1bDIVCkqR169YluRIA+IPf79e8efOSXcaUMM2yLCvZRYzX9evXdeHCBWVnZ8vhGHkd\nHgBIFEZSEyelQwoAMLVx4wQAwFiEFADAWIQUAMBYhBQAwFjcfjKEr776ShcuXJAknT59WqtXr1Y8\nHldzc7MqKioUDof10UcfKSsrS4WFhVq3bp3eeuutMbeZO3du0vq0bds27d+/X9u2bVNJSYl+++23\nlO/TE088oatXryoUCmnTpk2aNm1ayvfpscceUzQaVVtbmzZv3qyurq6U69NQ/fr222/1888/a8OG\nDTp37lxK7n9IEAvDOn78uPXBBx9Y27dvtyzLsmpqaqx9+/ZZW7dutYLBoGVZlvXCCy9Yly5dGleb\nZDh+/Lh14sQJ6/Lly9Z7771n1dbWWpZlTYk+1dTUWJZlWVVVVdahQ4emRJ96/z4nTpywjhw5ktJ9\nsqw/+tXc3Gy98cYb1vr16y3LSu39D5OL033DiEaj8vv9KikpkcfjkSR5PB6FQiE1NzcrNzdXkjRn\nzhyFw+FxtUlWn9auXSuv1zvge1OhTw899JACgYC+/vprPfXUU1OiTyUlJdq3b58+/fRTlZaWpmyf\n+vfr8ccf1969e/Xqq6/2fS+V+4XJRUgNw+/369FHH1VeXp6ampokScFgUAUFBcrNze2biqmtrU1e\nr3dcbZLVp6FMhT7V1NTos88+086dOzV79uyU71NnZ6fq6uq0efNmVVRU6JNPPknZPvXv14ULF9Td\n3a1Dhw4pEAjoyy+/TOl+YXJxTWoYdXV1euaZZ5SXl6fs7Gy9+eabamlpUUVFhVpaWvTOO+8oKytL\nq1evHnebZPVJkt5//31VV1ervr5ebW1tev7551O6T7FYTDt27NCaNWu0d+9e+Xy+lO9TZmamjh49\nqm+++UbBYFDPPfec7rzzzpTsU/9+3XPPPVq6dKkk6aefftLTTz+t4uLilO0XJhczTgAAjMXpPgCA\nsQgpAICxCCkAgLEIKQCAsQgpAICxCCkAgLEIKdhadXV13/M13d3devHFF3Xp0qUkVwWgF89Jwfb2\n7NmjRYsWKRgMavbs2SorK0t2SQB6EFKwvXg8rg0bNigrK0t79+5NdjkA+uF0H2wvEokoLS1Nra2t\nun79erLLAdAPIynY3iuvvKLNmzersbFRVVVVev3115NdEoAejKRga4cPH1ZxcbHuvfdePfLII+ru\n7tb333+f7LIA9GAkBQAwFiMpAICxCCkAgLEIKQCAsQgpAICxCCkAgLEIKQCAsQgpAICx/g/zIn/k\nCvDn1gAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "GeMpy.plot_data(geo_data)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array(['SimpleMafic2', 'SimpleBIF', 'SimpleMafic1', 'EarlyGranite'], dtype=object)" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "geo_data.formations" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "I am in the setting\n", "float32\n", "I am here\n", "[2, 2]\n" ] } ], "source": [ "di = GeMpy.InterpolatorInput(geo_data)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'DefaultBasement': 0,\n", " 'EarlyGranite': 1,\n", " 'SimpleBIF': 3,\n", " 'SimpleMafic1': 4,\n", " 'SimpleMafic2': 2}" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "di.data.get_formation_number()" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "geo_data_s = GeMpy.select_series(geo_data, ['EarlyGranite_Series'])" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "scrolled": false }, "outputs": [ { "ename": "AttributeError", "evalue": "module 'gempy' has no attribute 'set_interpolator'", "output_type": "error", "traceback": [ "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[1;31mAttributeError\u001b[0m Traceback (most recent call last)", "\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m()\u001b[0m\n\u001b[0;32m 4\u001b[0m \u001b[1;31m# use the GPU. Verbose is huge. There is a large list of strings that select what you want to print while\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 5\u001b[0m \u001b[1;31m# the computation.\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 6\u001b[1;33m data_interp = GeMpy.set_interpolator(geo_data,\n\u001b[0m\u001b[0;32m 7\u001b[0m \u001b[0mdtype\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m\"float32\"\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 8\u001b[0m verbose=[])\n", "\u001b[1;31mAttributeError\u001b[0m: module 'gempy' has no attribute 'set_interpolator'" ] } ], "source": [ "# Preprocessing data to interpolate: This rescales the coordinates between 0 and 1 for stability issues.\n", "# Here we can choose also the drift degree (in new updates I will change it to be possible to change the\n", "# grade after compilation). From here we can set also the data type of the operations in case you want to\n", "# use the GPU. Verbose is huge. There is a large list of strings that select what you want to print while\n", "# the computation.\n", "data_interp = GeMpy.set_interpolator(geo_data,\n", " dtype=\"float32\",\n", " verbose=[])" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "scrolled": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[9 9 9]\n" ] } ], "source": [ "# This cell will go to the backend\n", "\n", "# Set all the theano shared parameters and return the symbolic variables (the input of the theano function)\n", "input_data_T = data_interp.interpolator.tg.input_parameters_list()\n", "\n", "# Prepare the input data (interfaces, foliations data) to call the theano function.\n", "#Also set a few theano shared variables with the len of formations series and so on\n", "input_data_P = data_interp.interpolator.data_prep() \n", "\n", "# Compile the theano function.\n", "debugging = theano.function(input_data_T, data_interp.interpolator.tg.whole_block_model(), on_unused_input='ignore',\n", " allow_input_downcast=True, profile=True)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "scrolled": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1 loop, best of 3: 5.34 s per loop\n" ] } ], "source": [ "%%timeit\n", "# Solve model calling the theano function\n", "sol = debugging(input_data_P[0], input_data_P[1], input_data_P[2], input_data_P[3],input_data_P[4], input_data_P[5])\n" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "lith = sol[-1,0,:]" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "np.save('sandstone_lith', lith)" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "run_control": { "marked": false } }, "outputs": [], "source": [ "a = geo_data.grid.grid[:,0].astype(bool)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "run_control": { "marked": false } }, "outputs": [], "source": [ "a2 = a.reshape(50,50,50)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ True, True, True, ..., True, True, True],\n", " [ True, True, True, ..., True, True, True],\n", " [ True, True, True, ..., True, True, True],\n", " ..., \n", " [ True, True, True, ..., True, True, True],\n", " [ True, True, True, ..., True, True, True],\n", " [ True, True, True, ..., True, True, True]], dtype=bool)" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a2[:,:,0]" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 696000. , 6863000. , -20000. ],\n", " [ 696000. , 6863000. , -19551.019531],\n", " [ 696000. , 6863000. , -19102.041016],\n", " ..., \n", " [ 747000. , 6950000. , 1102.040771],\n", " [ 747000. , 6950000. , 1551.020386],\n", " [ 747000. , 6950000. , 2000. ]], dtype=float32)" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "geo_data.grid.grid" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2500" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "50*50" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "geo_data.data_to_pickle('sandstone')" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "ename": "NameError", "evalue": "name 'a1' is not defined", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0ma2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0ma1\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;36m2500\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mNameError\u001b[0m: name 'a1' is not defined" ] } ], "source": [ "a2 = a1[:2500]\n", "\n" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [], "source": [ "g = geo_data.grid.grid\n", "h = geo_data.grid.grid[:2500]" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1 loop, best of 3: 1.73 s per loop\n" ] } ], "source": [ "%%timeit\n", "eu(g,h)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "def squared_euclidean_distances(x_1, x_2):\n", " \"\"\"\n", " Compute the euclidian distances in 3D between all the points in x_1 and x_2\n", " Args:\n", " x_1 (theano.tensor.matrix): shape n_points x number dimension\n", " x_2 (theano.tensor.matrix): shape n_points x number dimension\n", "\n", " Returns:\n", " theano.tensor.matrix: Distancse matrix. shape n_points x n_points\n", " \"\"\"\n", "\n", " # T.maximum avoid negative numbers increasing stability\n", "\n", "\n", " return sqd" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [], "source": [ "x_1 = T.matrix()\n", "x_2 = T.matrix()\n", "\n", "sqd = T.sqrt(T.maximum(\n", " (x_1**2).sum(1).reshape((x_1.shape[0], 1)) +\n", " (x_2**2).sum(1).reshape((1, x_2.shape[0])) -\n", " 2 * x_1.dot(x_2.T), 0\n", "))\n", "eu = theano.function([x_1, x_2], sqd)" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "run_control": { "marked": true } }, "outputs": [ { "ename": "MemoryError", "evalue": "\nApply node that caused the error: Elemwise{add,no_inplace}(Reshape{2}.0, Reshape{2}.0)\nToposort index: 15\nInputs types: [TensorType(float32, col), TensorType(float32, row)]\nInputs shapes: [(125000, 1), (1, 125000)]\nInputs strides: [(4, 4), (500000, 4)]\nInputs values: ['not shown', 'not shown']\nInputs type_num: [11, 11]\nOutputs clients: [[Elemwise{sub,no_inplace}(Elemwise{add,no_inplace}.0, Elemwise{mul,no_inplace}.0)]]\n\nBacktrace when the node is created(use Theano flag traceback.limit=N to make it longer):\n File \"/home/miguel/anaconda3/lib/python3.6/site-packages/ipykernel/kernelbase.py\", line 228, in dispatch_shell\n handler(stream, idents, msg)\n File \"/home/miguel/anaconda3/lib/python3.6/site-packages/ipykernel/kernelbase.py\", line 390, in execute_request\n user_expressions, allow_stdin)\n File \"/home/miguel/anaconda3/lib/python3.6/site-packages/ipykernel/ipkernel.py\", line 196, in do_execute\n res = shell.run_cell(code, store_history=store_history, silent=silent)\n File \"/home/miguel/anaconda3/lib/python3.6/site-packages/ipykernel/zmqshell.py\", line 501, in run_cell\n return super(ZMQInteractiveShell, self).run_cell(*args, **kwargs)\n File \"/home/miguel/anaconda3/lib/python3.6/site-packages/IPython/core/interactiveshell.py\", line 2717, in run_cell\n interactivity=interactivity, compiler=compiler, result=result)\n File \"/home/miguel/anaconda3/lib/python3.6/site-packages/IPython/core/interactiveshell.py\", line 2821, in run_ast_nodes\n if self.run_code(code, result):\n File \"/home/miguel/anaconda3/lib/python3.6/site-packages/IPython/core/interactiveshell.py\", line 2881, in run_code\n exec(code_obj, self.user_global_ns, self.user_ns)\n File \"\", line 6, in \n (x_2**2).sum(1).reshape((1, x_2.shape[0])) -\n\nDebugprint of the apply node: \nElemwise{add,no_inplace} [id A] '' \n |Reshape{2} [id B] '' \n | |Sum{axis=[1], acc_dtype=float64} [id C] '' \n | | |Elemwise{pow,no_inplace} [id D] '' \n | | | [id E] \n | | |TensorConstant{(1, 1) of 2} [id F] \n | |MakeVector{dtype='int64'} [id G] '' \n | |Subtensor{int64} [id H] '' \n | | |Shape [id I] '' \n | | | | [id E] \n | | |Constant{0} [id J] \n | |TensorConstant{1} [id K] \n |Reshape{2} [id L] '' \n |Sum{axis=[1], acc_dtype=float64} [id M] '' \n | |Elemwise{pow,no_inplace} [id N] '' \n | | [id O] \n | |TensorConstant{(1, 1) of 2} [id F] \n |MakeVector{dtype='int64'} [id P] '' \n |TensorConstant{1} [id K] \n |Subtensor{int64} [id Q] '' \n |Shape [id R] '' \n | | [id O] \n |Constant{0} [id J] \n\nStorage map footprint:\n - , Input, Shape: (125000, 3), ElemSize: 4 Byte(s), TotalSize: 1500000 Byte(s)\n - , Input, Shape: (125000, 3), ElemSize: 4 Byte(s), TotalSize: 1500000 Byte(s)\n - Reshape{2}.0, Shape: (1, 125000), ElemSize: 4 Byte(s), TotalSize: 500000 Byte(s)\n - Reshape{2}.0, Shape: (125000, 1), ElemSize: 4 Byte(s), TotalSize: 500000 Byte(s)\n - Constant{0}, Shape: (), ElemSize: 8 Byte(s), TotalSize: 8.0 Byte(s)\n - TensorConstant{1}, Shape: (), ElemSize: 8 Byte(s), TotalSize: 8.0 Byte(s)\n - TensorConstant{(1, 1) of 2}, Shape: (1, 1), ElemSize: 1 Byte(s), TotalSize: 1 Byte(s)\n - TensorConstant{(1, 1) of 0}, Shape: (1, 1), ElemSize: 1 Byte(s), TotalSize: 1 Byte(s)\n TotalSize: 4000018.0 Byte(s) 0.004 GB\n TotalSize inputs: 3000018.0 Byte(s) 0.003 GB\n\n", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mMemoryError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m/home/miguel/anaconda3/lib/python3.6/site-packages/Theano-0.9.0rc4-py3.6.egg/theano/compile/function_module.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m 883\u001b[0m \u001b[0moutputs\u001b[0m \u001b[0;34m=\u001b[0m\u001b[0;31m\\\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 884\u001b[0;31m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0moutput_subset\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;31m\\\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 885\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0moutput_subset\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0moutput_subset\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mMemoryError\u001b[0m: ", "\nDuring handling of the above exception, another exception occurred:\n", "\u001b[0;31mMemoryError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0meu\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0ma2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", "\u001b[0;32m/home/miguel/anaconda3/lib/python3.6/site-packages/Theano-0.9.0rc4-py3.6.egg/theano/compile/function_module.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m 896\u001b[0m \u001b[0mnode\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfn\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnodes\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfn\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mposition_of_error\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 897\u001b[0m \u001b[0mthunk\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mthunk\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 898\u001b[0;31m storage_map=getattr(self.fn, 'storage_map', None))\n\u001b[0m\u001b[1;32m 899\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 900\u001b[0m \u001b[0;31m# old-style linkers raise their own exceptions\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/home/miguel/anaconda3/lib/python3.6/site-packages/Theano-0.9.0rc4-py3.6.egg/theano/gof/link.py\u001b[0m in \u001b[0;36mraise_with_op\u001b[0;34m(node, thunk, exc_info, storage_map)\u001b[0m\n\u001b[1;32m 323\u001b[0m \u001b[0;31m# extra long error message in that case.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 324\u001b[0m \u001b[0;32mpass\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 325\u001b[0;31m \u001b[0mreraise\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mexc_type\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mexc_value\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mexc_trace\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 326\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 327\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/home/miguel/anaconda3/lib/python3.6/site-packages/six.py\u001b[0m in \u001b[0;36mreraise\u001b[0;34m(tp, value, tb)\u001b[0m\n\u001b[1;32m 683\u001b[0m \u001b[0mvalue\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtp\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 684\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mvalue\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__traceback__\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mtb\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 685\u001b[0;31m \u001b[0;32mraise\u001b[0m \u001b[0mvalue\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwith_traceback\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtb\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 686\u001b[0m \u001b[0;32mraise\u001b[0m \u001b[0mvalue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 687\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/home/miguel/anaconda3/lib/python3.6/site-packages/Theano-0.9.0rc4-py3.6.egg/theano/compile/function_module.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m 882\u001b[0m \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 883\u001b[0m \u001b[0moutputs\u001b[0m \u001b[0;34m=\u001b[0m\u001b[0;31m\\\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 884\u001b[0;31m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0moutput_subset\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;31m\\\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 885\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0moutput_subset\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0moutput_subset\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 886\u001b[0m \u001b[0;32mexcept\u001b[0m \u001b[0mException\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mMemoryError\u001b[0m: \nApply node that caused the error: Elemwise{add,no_inplace}(Reshape{2}.0, Reshape{2}.0)\nToposort index: 15\nInputs types: [TensorType(float32, col), TensorType(float32, row)]\nInputs shapes: [(125000, 1), (1, 125000)]\nInputs strides: [(4, 4), (500000, 4)]\nInputs values: ['not shown', 'not shown']\nInputs type_num: [11, 11]\nOutputs clients: [[Elemwise{sub,no_inplace}(Elemwise{add,no_inplace}.0, Elemwise{mul,no_inplace}.0)]]\n\nBacktrace when the node is created(use Theano flag traceback.limit=N to make it longer):\n File \"/home/miguel/anaconda3/lib/python3.6/site-packages/ipykernel/kernelbase.py\", line 228, in dispatch_shell\n handler(stream, idents, msg)\n File \"/home/miguel/anaconda3/lib/python3.6/site-packages/ipykernel/kernelbase.py\", line 390, in execute_request\n user_expressions, allow_stdin)\n File \"/home/miguel/anaconda3/lib/python3.6/site-packages/ipykernel/ipkernel.py\", line 196, in do_execute\n res = shell.run_cell(code, store_history=store_history, silent=silent)\n File \"/home/miguel/anaconda3/lib/python3.6/site-packages/ipykernel/zmqshell.py\", line 501, in run_cell\n return super(ZMQInteractiveShell, self).run_cell(*args, **kwargs)\n File \"/home/miguel/anaconda3/lib/python3.6/site-packages/IPython/core/interactiveshell.py\", line 2717, in run_cell\n interactivity=interactivity, compiler=compiler, result=result)\n File \"/home/miguel/anaconda3/lib/python3.6/site-packages/IPython/core/interactiveshell.py\", line 2821, in run_ast_nodes\n if self.run_code(code, result):\n File \"/home/miguel/anaconda3/lib/python3.6/site-packages/IPython/core/interactiveshell.py\", line 2881, in run_code\n exec(code_obj, self.user_global_ns, self.user_ns)\n File \"\", line 6, in \n (x_2**2).sum(1).reshape((1, x_2.shape[0])) -\n\nDebugprint of the apply node: \nElemwise{add,no_inplace} [id A] '' \n |Reshape{2} [id B] '' \n | |Sum{axis=[1], acc_dtype=float64} [id C] '' \n | | |Elemwise{pow,no_inplace} [id D] '' \n | | | [id E] \n | | |TensorConstant{(1, 1) of 2} [id F] \n | |MakeVector{dtype='int64'} [id G] '' \n | |Subtensor{int64} [id H] '' \n | | |Shape [id I] '' \n | | | | [id E] \n | | |Constant{0} [id J] \n | |TensorConstant{1} [id K] \n |Reshape{2} [id L] '' \n |Sum{axis=[1], acc_dtype=float64} [id M] '' \n | |Elemwise{pow,no_inplace} [id N] '' \n | | [id O] \n | |TensorConstant{(1, 1) of 2} [id F] \n |MakeVector{dtype='int64'} [id P] '' \n |TensorConstant{1} [id K] \n |Subtensor{int64} [id Q] '' \n |Shape [id R] '' \n | | [id O] \n |Constant{0} [id J] \n\nStorage map footprint:\n - , Input, Shape: (125000, 3), ElemSize: 4 Byte(s), TotalSize: 1500000 Byte(s)\n - , Input, Shape: (125000, 3), ElemSize: 4 Byte(s), TotalSize: 1500000 Byte(s)\n - Reshape{2}.0, Shape: (1, 125000), ElemSize: 4 Byte(s), TotalSize: 500000 Byte(s)\n - Reshape{2}.0, Shape: (125000, 1), ElemSize: 4 Byte(s), TotalSize: 500000 Byte(s)\n - Constant{0}, Shape: (), ElemSize: 8 Byte(s), TotalSize: 8.0 Byte(s)\n - TensorConstant{1}, Shape: (), ElemSize: 8 Byte(s), TotalSize: 8.0 Byte(s)\n - TensorConstant{(1, 1) of 2}, Shape: (1, 1), ElemSize: 1 Byte(s), TotalSize: 1 Byte(s)\n - TensorConstant{(1, 1) of 0}, Shape: (1, 1), ElemSize: 1 Byte(s), TotalSize: 1 Byte(s)\n TotalSize: 4000018.0 Byte(s) 0.004 GB\n TotalSize inputs: 3000018.0 Byte(s) 0.003 GB\n\n" ] } ], "source": [] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'/home/miguel/PycharmProjects/GeMpy/Prototype Notebook/sandstone.vtr'" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from evtk.hl import gridToVTK \n", "\n", "import numpy as np \n", "\n", "# Dimensions \n", "\n", "\n", "nx, ny, nz = 50, 50, 50 \n", "\n", "lx = geo_data.extent[0]-geo_data.extent[1]\n", "ly = geo_data.extent[2]-geo_data.extent[3]\n", "lz = geo_data.extent[4]-geo_data.extent[5]\n", "\n", "dx, dy, dz = lx/nx, ly/ny, lz/nz \n", "\n", "ncells = nx * ny * nz \n", "\n", "npoints = (nx + 1) * (ny + 1) * (nz + 1) \n", "\n", "# Coordinates \n", "\n", "x = np.arange(0, lx + 0.1*dx, dx, dtype='float64') \n", "\n", "y = np.arange(0, ly + 0.1*dy, dy, dtype='float64') \n", "\n", "z = np.arange(0, lz + 0.1*dz, dz, dtype='float64') \n", "\n", "x += geo_data.extent[0]\n", "y +=geo_data.extent[2]\n", "z +=geo_data.extent[5]\n", "\n", "# Variables \n", " \n", "litho = sol[-1,0,:].reshape( (nx, ny, nz))\n", "\n", "\n", "\n", "gridToVTK(\"./sandstone\", x, y, z, cellData = {\"lithology\" : litho},) " ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "-20000" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "geo_data.extent[4]" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa0AAABxCAYAAABm3fzXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VNX9+P/XvbNm3zPZgWwgWQBFVkFM2FwKUhHc+inF\nBQviT1wq6udji1Rt+7F1e1SsRauo/WmpUvnUhQJlqWxiICD7GpKQPWGyTDLbnfv9Y8gkk50lG5zn\n48GjnZlz733PnTjvOeee876SqqoqgiAIgtAPyL0dgCAIgiB0lUhagiAIQr8hkpYgCILQb4ikJQiC\nIPQbV1zScjqdFBYW4nQ6ezsUQRAE4TLT9nYAl1tJSQnZ2dmEpN6FRh9w6TuUJHRGX1QJbHVV1FeX\n4Wioxm6txW6twW6tISBsEL7+4Vjrz1FfU4qtvgqrpRJr3Tl0Rn9kWcZhqyMwYhCqqmKznCM4MpWE\nodNQVRXFZUej09BgqaDeXEpNZR4NVWdxWCqw2xtwuVy4VBUkwDPZU0KSZYJNg4kYMJIBmVOQtTpk\njQatYiPFWkJcQykuxUmeFMTpoIFYg42oshM1tARN4Dn0TivXnWwg0uCDVW8gbyBUhVaiqC5ciouC\nd75h17/OtnlaDBjxlQIZEXMT0SGDwG4lpL4Yf62D/79sC2fqK7lr1FTq1TDqf3CArc2doEvXIwVI\nqLUqjgN2r3ZGo4aUodH4+hvwUa1YzhZTY1E4XQF2RUKvUdEFV/JFxeZWux5FNic5QCWlJJHGIOma\nTj9qs1rB92xmODcQLkV5npd1WoLjo9H5GHA02KgpKScwKsLz2FxQjMvR/o+kSrWUvfzH6zkJCV/8\nsVDb5jZhhmhGxk4hSArrdP9C37fetbq3Q7hiSFfalPfCwsKLSlrm0qOUnN4FqgsVFcVhxVpXgaLY\ncSkOVJdy2WIMCBtAxuRFHN76F84VHwEu/iPQGQNJHjWbuqoC/EJjibOfw1h2FKOkUO90cn1kFDfG\nxFPggI8sVdQN8yV2dDyTd50jwBxPmX8yLknj2Z9N7+J0YjHmiEJcTpUf/rqLY1/ta3VcWdIw5fqf\nkJEyEV8d+OpBKykEVO7n4b+/iQrcfc01zLw2gwMhevYcVLFvbspI+kkG9DcakfRN+1TtYN9ixb7Z\nxvUTkghOTUbWapCl88dUFUIqj6MrOkGVBbQDk1l19DscLieV1WWcLjlGgikVs7kSky2WUgqo4Rwy\nMhO4DV3zg7XhuPoDZzhKHEkMkUYAYEpLxpSWhKx1nyOt0YDOaMBpteOwWgFwORVKD56k9OCJNvfr\nUO1sYa3Xc4MHjMKlKBwvzOkwpqGDxnFDxiwspyrb3b/Q94mkdflccT2tixVsGowkaTixdzU2S9UF\nbClxoUmntvIM3/1jGSHRg9Ea/HDa6i5o++Yc1hoOb33P8/hIs6geHJpJvcPBU9u3sLe8lPTICGZk\n3UzcrnMEnYujKHBwq/0Z7DIDT8WRB5wLLyTzvjEYg33Z/9cdXu1cqsKJs7nEhscSEZ+ELIEka9hi\ntnrOxu6iIu5LS2N4pQ1pmIEcDNg329DdqEeXpadZrnTHrAfDFCMjIhIIDElBlvAkLACXpKEyfAgm\nGcJkOBs0hJtDBmDQ+7Bq3ZsAJMemMnbaDHK2befsyVPu7XBxmD1kMqbDc1lBEQDlFDFYHU5UegrR\nw1I9r2uNBvS+RgB0vobz59+KrNV42rWVWOLShxKavxtzbRmSJKG4nMRFpJCScJ1X0tLKeiZdN5fK\n6iLKzuVztvw4h05v53hBDmPSbiP2mkFUHM7r8D0IwpVOJK1mgiKTGZ79BAWH/0XRia3NhuI6cnG9\nJKfNQnnenovatitU4J1D+z2Pg/UGnh85Ck41EFgncdCU0u62RrtKWGkc5pBiVI1C6m2ZGPyNfL9y\nM6qr6f2ePJvLybO5JMUMIfvaW0lNGMrBvFzP6yeqq6mrrUMbEkBauZ0jN/hjz7FRaDzFAEcKsqZ1\nz8fglAmKTEJ1eCes5spDU1FVFVkCg94HAJvD3YvTaw0EGLWkjRzPhpN/9WxTRiFl6lkipdg291mv\n1nmG6mw0UK+tw5SW1NRAAp2P0WsbrY8Bp9WGev5vwJSWRPnRPFzNrqfKOi2mtCSiLIO4ZuBo7E4b\nuw99zfYf1jJ44CgARg6ZxvdH1uF02QkLimF46k2oqsqxghy27PmUGksl/9n3GUF+4QyUhxCumJCk\ndk6OIFzhrriJGJdKo9WTMmouY+/4Nf5hCZ22jxxwPSHRQ9Ea/Hsguotnttt4aOMGdh88TY0h2mtI\nsC1+DRL+NWGexwMmpDLuoeloNTrCgmK5M/tJEiLdvYuTRUd455+/57XVL3DodK7XfnYWl6BXVLQu\nlQENDmyjGzi4YTdava7N4w6sDkVFg9TeX6YEdo0Rh9adrBq/u+3nk5ZBZ0SSQKNx4MTutekhcrCr\n1jZ3W0Gx1+PawFrPkCCAVq+nZZ6QJNAYmt6HrNUQnBDl1SY4PhpZq2FY8o2MTruV0Wm3YtT70WCr\n5diZ75k4/E4mXTeX2GD3j4hDp3ec37fE4ISR/Oy2FxmfOQutRk+1pYJ9yrec5XQ7J0cQrnz9qqd1\n7NgxVq5cSWBgIIMGDeLee+/tluNIGg1BpiTG3fkSp/Z8wcnv16C6nCReNwtJkig5uQvLOfcEhZjk\niQTHpmLwD8RmqaLsTC4nv/s7tvpz7m+1dnprMUNupKb8FHWVBd3yHtpSZbPhq9Hg0Bg7bSu7VLRO\n756QKS2OaTf+hNwf/sOAqKFUmfMYnjKSk2cPczBvL4XlZ1rtZ0dJMTcnJWHTSsj1dvbs2oosy0jt\ndKP0ivuYUgejrqoke15sbGZ3nu9p6dzDdiWVea22c2J3DxOqY1v1VFomreIG78TQXs9GapFdW/bG\ndD7ueGIj3UnJqPcla+Q96LVGkuKGe/Y7MmUqI/TZJMeP8N5eq2dsxo9ITxrP1r1/J6/wAFHO+DZj\nEYSrQb9KWitXrmTJkiVER0fzwAMPcOedd6LXd3xx/WKoinvShazRknz9HUQljebApj8hSRIpo+eQ\nMnoONouZvNyvcbkcqIqCJEkY/cNISMtmQMYUKvL3cXDzSkb9+FeYS45ScWYfxcd3Imt1DMicxuBx\nd6O6XJTn7yMwfCAb//wAPkY/Hhs5HgmosTXgVFW25B3D4rBjdyqUWGoA8DEGEDU0i4r8/dRWnMYv\nJIbbs37CxJqjnCgv5JzVisXp4IaYWP52/Cg3DxyEydeHM7U1pEVEgNJ2b6M5lyzh1Hr3VFAgSBfG\nzImLADCFDuDTja9zz+SHuHn0LFatW0GZ2fvLP6eiHMXuwChr+b9N32OpqcUY4Nvuce0aO37Ojkdm\nJdXV9Fk1budwvye91p0kCstPtdpOg45azJRQQDRNvWin6uAc5V5tKxuKsTRU4+cT5D5OOwGpzWIB\ncDRYWzxuPW1y6KCxrZ4L10cTMjCmzWMABPiGcuv4hzix7TvqzlS0204QrnT9KmlVVlYSFeUefgkK\nCqKuro7Q0NDLfhyHtR6jGuL5FewfGsfoWcuoKW/69a33DSI6aQKoaqv2qstF5KDrCIvPBEkldsiN\nxA65kcSRs8j9+lUGj70bVXVPXIgceC311WXEXnMjwyf9jMzKHWjVpmsi2YMGk3eukte/2wRASpgJ\n38wZxA2/BYu5lGM7/kpG9sNYtFqCis5xs68eUFFcLjSyzC9HjQNAkl0MCAqkKkAmsK4YWU3vcIjQ\n4qNSF1jpeawqKq5qFxXbCoi61T0sGBGajKqqfLR+BbMn/oQA30DO1VbgUBye7QYGBLCnohyzw8bh\nw+7z197QIEBeUBWhZQqqS2578FoFvWJFVVVUSeNJbnanO8HqdQZUFfJLW0+IcOFkDLcit9hxJaWe\n61LNnSjYy7DUSQA47XZ0vj5eQ4SqCoqt6b26nArm/BKvfZgLiokbOdRrqLFVXE6Fs3uPEBRn6rRd\n/Vlzu68LwtWgX13TioqKoqTE/aVgNpsJCQnpngOpKra6Gq+nJFkmyNR0Yd5WV9PUHWijvUtR0Gh1\naBonG6gqAaFxjLljmfsLV1Hc0+hVFaN/KMOn/39gDGR/QBItfXXiIIU17i+r+KSRxA2/BQAf/1CG\nT38Mrd4Hh6xjX2AK9Rp3T0Mjt/hoVYkGvcz+ZD9OxWqJqj3e7tu36iUqIwtRNYq7K6OAWqti32zF\nZXFSevAkADqtgZjza89Wb1lFeGCkV8ICiPL1I8LHhzf2NU2b19jb/2K2aV1Ul7n372qntxVRdYyo\n6uNer9ua9bQsNpWqymKiSEDT7HeZjIY6zMgthvSs1ONHoNdzRvwoy2823Km27kU5G2xeya704Emv\nSRgALkfT+WpP6cGTOBusXWrXcv+CcLXpVz2t+fPn8+qrrxIYGMjUqVO7dQaV/XwSMvgHeh1HPZ+g\n7C2SVMv2qkvBBUgaGRQX0vkkojX6oyoKLsWJrboGrcGI3s/fnQAlif2ByQBk1p5EqzrZcuY4m/KO\neY5THpjgHthS3V+XkrkWggO8th177gD+SgNNg2cStbIvu6JiORLXwKFEO5N3FRJlltpcp5U3qJhz\nEYXuzZ3gqlSwrm3wrLNqnNZtSksiNmIwhWXux7uOeC+gBfiutIQzljpsStM6N7lWg229td11Wrs2\nH+f6CS6CU5OhjXVaUtEJKi3gPxBqI1JQJAnH+Z6Wy6XjxPdHybSOQYOW7ayjAfeSAgUn7qtg3mIY\nyHH2ez1nxUJAqR/F+4551mk5re73f6HrtJqfr+Y9qZbbdbWdIFzNxOLizpyviCFpNKiKgsNa38kF\nlxbtbQ3oDD5IGg2yVovL6Wy9H1nGGBiMrNEiSTIul4LO5SS4/ACfff2W5wsZIGPyI0QOHImzoQGH\n1YLv5zvB3wftxOuQAv3Bz4jeaSPRWUFCXRFqg40zNRpOm1JxBgSg4kQ2FqLxrcJormHEWQd6KZIG\nycBJjUJV9DmIcFffcJW7cB504Mi1t1nRQtZqsYbb2Fq8xut5vSxjd7mYnZTMvwsLqLJ5bxxBDMOk\ncZ1WxDAYNKSmRePrr8dHtVF3tphai8KpCnAoEjqNSqJJiyE6jD98/xYA47W34KM0XTPbrf6baqow\n4ouVegaQSoqU6RVPiVrAAXa1en9DGEGclISs1RKcEIXOx4ijwUpNUTmBMRGex+b8ki71gFrup73t\nutpO6D/E4uLLp1/1tHqFquJosFxS+063d7mwmr0XNNcpDrZu+cQrYQE46xtoqPKeNEBdA86vvm06\nHvDD+X9NSlsdth7YCEDbpZo643I6kYtVd8+yWSK3u9yTE2zaAPQaLbF+Ws5ams6BZ8jOBo6cFpM9\nmrHZFH7YU9ji2aaekkOROFqkYDub3/SqU/LqTOlwD5f64o+V+lYTLgAc2IkkhrLzi4sBwolBQfG8\nz6pT3nG0fNwVbe3nUtoJwtWoX13Tuprk7fuC+poSWg5naXWdT1fvSVpJS0h4uNdzQxLSuXnULBLT\nbmXpT19jeLy77p9fQAAaNGhpfyLGxXAP+52Pp8XvMB3u8cfAgHCyB97F1MSfIOu828SRSC3VXs/V\nUkUC7S/AFgShd4ieVh+kqi5iUm9kYOYM9v7rt9itNWi0RhSnFY3Op7fDayUpIoKq8nIMOiM2hxVL\nQx3Z192KLMk4XS625rsnfYxKHUx1ji9mLu+UbQNGRjABBafXxAuA4HATxRVnMIT6M2L8VKD1NaJ6\namnAuzdsw0otZgLppsk+giBcFNHT6oMkScbHP4Laqjzs1hqQZDJvWoxfUAyaPtbTAkjRRxIdGsvC\nmb9AlmQKyvP47rB7uPL42cPU1lcjSzLX+sQTKkUygNY1Dy+FRtISJpmIlGK9Js2Y0pIJj3Ovx2qw\nNlVTb6wVaEpzT1wppxgDPp5emRYdRnxbLTgWBKH3iaTVh5Xnu4upBkem4hMQydAbHkSnb39hbm8Z\nHn4Nc2+aT3zkQMal3wTA19+twak4yDnmLks0OD6deD93AtF0UkLqcmis+ed7vrxWg631LUDcs/S0\nBBPGeKZjwN2L1aFnHNMJIaLb4xQE4cKI4cE+SnHaqTx7AIDIhOvwXdN6dltfodEHkmAaBMDNo2ZR\nYzEzfdTtuFwufjjlTrzXDR6DIht6LKbGmn8+BvcM0oY2Kuk31gp0nTp/TazZpFBZkkXSEoQ+SPS0\n+qiqoh9wKXY0WgMhMWm9HU6HGuqabkfiY/Bl3vRFRIXGcjAvF5vDhkFnIG3gCBosbd0Jsns01vzz\n8fS06tosxdSyVqAgCH2bSFp9VHm++7YlYbGZaDSXd7bd5ZazpxjZpbQqhNQ4NJiReB1GjZbvc3ru\nGlFjzT8foztpKS4nDmfrmostq1wIgtC3iaTVB9kbajCXuWfcRSRc18vRdK7W4kI5445XPf+vtqGG\nI/nu4c1rU8agnDlOXb2r/Z1cZuaCYlxOBR9DAEH+4USFDWq15q1lrUCD5L6m5ecbTGhifKup8YIg\n9D7xX2UfVF64F1DR+wQTGD6ot8Ppkn9bzzLFpQV5EC5JQ+7J3bhUFwG+QaTG6Pn38Z67BQs01fyL\nHpbKgzN/12ab5rX8TGnJJNjOUXm8BFPUQBLGZBA3cqgonyQIfYxIWn1Q46zBiPhrW92vqU8ygP5G\nI1v1pQTXnWNI9TXkHHUPDSakD2DbwHL0MUbsO2xtloO63GSdluD4aABqisrxN4Uha5rOY8t1Wqa0\nZKKHpRLmHMC4YTPRyO7/LBqnxgMicQlCHyGSVh9jqS6ivtp97Sci4dpejqZrdOl6JD1klmWiJZYi\nSyn5pe77WU1Mvo/IMgP7I/ejS9d3WLbpcjClJbcuOKso1BSVYyk/16qWX+PUeDh/Ly5t6xmOprQk\nyo/mifp/gtAH9IOf8VeXxgkYfsFx+AaaejmarpECJDLLMtEQj4rM/iNbAQgPiSUqIgkN8WSWZSIF\ndF9VfmjqMbW8J5Ws0RAY456+XnWq0Cv5NE6N70jj1HhBEHqfSFp9iKq6qDzuXo8Vaw7Dd80uz7++\nzGjVoSXW89jXGEBwQATDrrnRU6FCSyxGa/d17Jv3mNrTuJi4ucap8Z0RU+OF/uA///kPjz/+OPX1\n9Ze8r5dffhmAVatWXfK+LicxPNiHWOsqceFCQiKK+B4/fuO1IJ2PAUeDzT0Dz9H5kFhScDJqs98/\no4ffwvXDpqMoTduqyCQFpbCX3G6J/UJ6TM0rqDdOje+MmBov9Adr167FbDazfPlyYmJikGWZmTNn\n8swzz3Dbbbexd+9eMjIyyM3NJTU1FZvNRnBwMNnZ2bz11luEhoYSGRlJRkYG3377LQcPHmTnzp3M\nmTOHX//610RHR2M2m3nuueeYNWsWc+fOZfv27fz2t7/Fx6dn6qKKnlYf4hMQwQRuYyST0Es9+8ve\nlJZM+qxsEsZkED0slYQxGaTPyvbU5+uIwccHWqzSkiUZnbbZHR5Rz7frHhfbY2qcGt+RllPjBaGv\nGjt2LNdeey233347ixcv5ujRoyiKQkJCAnPnzgVg2rRpTJkyhYCAABYuXMiuXbswGo0EBQXh6+vL\nt99+S2ZmJgMGDCAtzV3Y4Ntvv+Xaa69l0aJFAJw9e5bQ0FDuuusukpKSyM/Pbzemy00krT5GlmSC\npLAePWa714JaFJZtj62u4XzOUs+v05Ka/XM/jwo2S0P3vAEuvsfUODW+I+I290J/07z6iyRJ+Po2\n1SzV6/XIsuz5X0VR+Mc//sHEiRP56U9/irOdv/WW+zQY3D8UZVnG5eq5NZhiePAq19VrQR3Nnju6\n4zAjhwxBlWRoNUVfQlVdSKqLI9sPXaaoWzMXFBM3cmiHQ4Tt9ZjEbe6FK0liYiL//Oc/2b17N5mZ\nmchy532ToUOH8umnn5Kbm4skSZw6dQq73c7u3bsBmDBhAi+++CKlpaUYjUZiYmK6+220S1LbKsjW\njxUWFpKdnU1I6l1o9AG9Hc4F6+lJF6GJ8SSMyei0Xf7O/R3eTXfsfdPQhrdfYNZZUc6Oj9ZdVIxd\n1dhjbE/xvmMdJiBxm3uhu6x3re7tEK4YoqfVS/rKjMDLMXtO1mmx2CHAbkfW6/G+27KKy27HYncn\nhe5MApfaYxK3uReEvk8kravc5Zg91zhzz251opM0yFoZd+JScTldOKzONmfudYfSgycoP5onekyC\ncIUSSesqdynXghrpfAzojEZ0vu5em0tpdlFWwvN8T611Ej0mQbhyiaR1lWteWLY9nc2eUxwOtJ0M\nM2p9DCh2x0XHKQhXkwabk33Hy6mx2An00zMsJQIfg/i6BpG0ul1fuXbVkUuePad2rTyT2uqOW4Ig\ntLR+1xk27M7H7mhaP7hm8wkmX5/AlNEDejGyvkEkLQG4tGtBGr0WR4MVvW/7w3+OBitavb7d1wVB\ncCesr7afbvW83aF4nr/aE9cFJy2LxYKfn193xCL0sou9FuRosOG0uid06HyMSM06XqrqTlhOq02U\nQhKEDjTYnGzY3XFliQ2785kwPBZjF4YKP//8c7788ksSExMBiImJ4Wc/+1mbbQsLC1mxYgUvvvhi\nq9feeecdCgsL0Wq1VFVV8fDDDzNkyJAuvKMmq1evZurUqfz+97/nhRdeuKBtW+rwnd9+++38+te/\nJj093fPcE088wdtvv31JBxWuLI2TOZxWd/LSGPRIkoSqqig2961IRCkkQejYvuPlXkOCbbE7FPYd\nL2d0enSX9jljxgxmzpwJQENDA88++ywxMTEUFxfzwgsvsGjRIpKTk7nzzjsB9/f7L37xC0wmEw8+\n+CD33HMPNpvNk2gsFgtWq5XCwkKefPJJJk2aRFZWFqtWrSIiIgKdTsfChQu57bbbuPfee9m9ezeL\nFy8mJyeHCRMmsHPnTtauXUt8fDwbN27Ex8cHrVbLggULunyeOkxaAQEBvPfee4wdO9bzpoS29Ydr\nV92l5WSOxkTVnCiFJAgdq7F07V5zXW0H7gK6Bw4cAGDEiBHExMTg6+vLyZMnKSsrw2KxsGjRIior\nKwF3R+Vvf/sbU6ZMISUlhZMnT3Ltte77+n311Vfs3bsXVVWZN28eJpOJhx9+mOLiYsLCwggKCuLr\nr79m4cKF+Pv7c/fdd+Pj48O+ffs88cTExDBjxgweffRRkpKScLlcHDp0YZVyOkxafn5+/OEPf2DF\nihX86le/4rnnnrugnTfas2cPb7zxBrNmzWLmzJmUlpbyyiuvEB4ejo+PD48++ijvv/8+hYWF1NbW\nsmDBAvz8/Dpt09jtFXqfKIUkCJcm0K9r13y72g68e1obNmwgICCAn/70p+zevRtFUdBoNF7V2ceP\nH88f//hHqqurufvuuyktLWXnzp2MHz+eW265hWnTpnmK5vr7+wPw3nvvcccdd5CUlMTatWsBMBrd\n17clSUJR2u493nvvvYSHh1NScmEjMF26pvXzn/+crVu38vOf/5zq6uoLOgBAREQEM2bM8Dz+5JNP\nmDt3LiNHjmTp0qUUFxezbds2/vznP1NQUMDbb79NZGRkp23aGn8Veo9Y2CsIF29YSgRrNp/ocIhQ\nr9MwLKX9cmktNe9pFRUVYbfbqa+vJyIignXrWpdVk2WZ8ePHs337dpKSkkhKSuLw4cM8//zz6PV6\nzGYz9957r9c2w4cP59133yUxMRGTycT333/fbjypqamsXLmS+fPn85vf/IawsDBCQ0Mv3/DgnDlz\nPP9/4sSJDBw4kJdeeqnTnb7//vvs2LHD83jJkiVer1dUVGAyue/KGxkZSUFBAaGhoQCYTCbKy8uR\nZbnTNkLfIxb2CsLF8TFomXx9QpuzBxtNvj6hS5MwAH784x/z4x//uMM2999/PwBxcXGeTkBGRobn\nu7d5m5Ya2996663ceuutXq+9//77AJ5eXmMczz77rKfN8OHDu/Q+Wurw3WdlZXk9TkhI6NIkjHnz\n5jFv3jyv55qPW0ZHR1NaWkp8fDzFxcUkJiZy7tw5wP1rIDY2loiIiE7bCIIgXEkap7O3XKel12l6\nZJ3W+vXr2bBhA8uXL+/W41yKHqny/uGHH7JlyxZkWWbq1KlMmjSJ3/3ud4SGhhIQEMCiRYv48MMP\nOXPmDLW1tSxcuBA/P79O2wwY0PoD7O4q71fzhAtBEC7OhVZ5t7ZREaOrPawrnbg1yQUSSUsQhAt1\nwUnLYeWHsqPU2uoIMPiTETkYo65n72beV4nULQiC0If8+9R2Np/ejl1pqtX5f0fWM2nQOLISx/Vi\nZH2DSFodEL0qQRB60r9PbedfJ7a0et6uODzPX+2JSyQtQRCEPsDqsLL59PYO22w+vZ1xCddh1HZ+\n89aysjJ+85vfEBMTQ21tLbGxsYSFhXHHHXd0OaZdu3aRk5PDwoULW702bdo0nn76ac+EvWeeeYaY\nmBgWL17cqu2bb75JfX09iYmJTJ06laCgIK/XP/roI44dO4bFYmH27NmMHTu23ZhE0hIEQegDfig7\n6jUk2Ba74uBA6RFGxg7rdH8HDx7EZDLx+OOPI8syp06d4t1330WSJHJzczGZTBQWFpKRkcG2bdt4\n+eWXWbJkCePHj8dsNhMfH09CQgIAeXl5fPLJJwQHB1NdXc3TTz9Neno669atIysri6qqKmw2d/3R\nsrIyXn31VWJjY6mtrWX+/Pls2rSJrKwsTzmnL774gsLCQsrKynjqqadITk7mvvvu48SJE6xevbrD\npCVfwDkVBEEQukmtra5L7Wpsli61mzRpEqmpqSxbtoznn3+enTt3el7LyMhg4cKF5OXlcc8995Ca\nmsrJkyex2+3MmTOHxx9/nK+//trT/tNPP8XpdOJwOCgqKqK2thadTkd8fDxnzpzh888/Z9asWYC7\nGobJZMLPz49vv/0Wk8lEamoqs2fP9uxv8+bNPPvssyxbtgx/f3/GjBlDVVUV7777bqvlUi2Jnhbi\n2pUgCL0vwODfpXaBhq7dZePYsWNkZWUxa9YsFEVh3LhxTJ48GQC93l3U2mBwDzPKsuxVbklVVVpO\nLP/Rj37EsGHDKCkpISDAPTN7zpw5rFq1itraWqZPn05ubi6ff/456enpTJ48mTVr1rQZW+O+FUXB\n4XBQXFzAQAWtAAAMWElEQVTMxx9/zNKlS1sNHbYkkpYgCEIfkBE5mP87sr7DIUK9Rke6qWu3BVFV\nleXLlxMREYHdbuehhx7i1KlTHW4jyzJ///vfKSgo8Cq9d9ddd/H6668TFRWFoig888wzgLtaUWVl\npScZAqSlpfHBBx9w7NgxkpOT+eabb1odJzs7m5deeonKykoee+wxfvGLXzB69GjefvttEhMTOyzQ\nLtZpIXpagiB0r66u02pv9mCjqck3duvswXnz5nlKMPVVoqclCILQRzQmpJbrtPQanVindZ5IWoIg\nCH1IVuI4xiVcx4HSI9TYLAQa/Eg3DenSNPdL1dd7WXCVJS0xDCgIQn9g1Bq6NK39anRVJS1BEIT+\nQGlowLz/B5w1NWgDAwnOzEDT7GaNVzORtARBEPqQ0g0bKdu4CcVm9zxX9I+1RGbfhGlydi9G1jdc\nsUnLZ10uWmfXb0stCILQ20o3bKT4q9Z3FFZsds/zV3viumKTliAIQn+iNDRQtnFTh23KNm4i/Ibx\naIyd36aku2sPzpw5k1GjRqEoCk6nkxdeeIE333yTsWPHkp+fz5dffkliYiIA48aN46abburycTsi\nkpYgCEIfYN7/g9eQYFsUmx3zvh8IG319p/vr7tqDISEhPPfccwDcf//9OJ1Or+PPmDGDmTNnXuTZ\naJ9IWoIgCH2As6ama+1qu9Zu0qRJmM1mli1bhqqqDBnSVEkjIyOD2bNnc8899/Dyyy9TXl7uVXvQ\n39+f+fPns2DBAqDt2oNms5lXXnmFqqoqoqOjkSTJ6/hr167lwIEDAEyZMoVRo0Z1Ke7OiKQlCILQ\nB2gDA7vWLqBr7bq79mBwcDBPPvkkACtXrmTLFu9KHqKnJQiCcAULzsyg6B9rOxwi1Bj0BA/L6NL+\nurv2oNls5re//S0AxcXFzJo1i4MHD3YptktxxdYejDk9RMweFAShT+hq7cH2Zg82ir5lWrfOHhS1\nBwVBEIQua0xILddpaQx6sU7rPNHTEgRB6GZd7Wk1UqxWzPt+wFlbgzYgkOBhGV2a5n41ED0tQRCE\nPkZjNHZpWvvVSCQtQRCEPsZuc3L6eAX1Fju+fnoGpYSjN4ivaxBJSxAEoU/Zuyuf3N0FOB1NU9C3\nbz7J8OvjGTE6oRcj6xtE0hIEQegj9u7K5/vtea2edzoUz/NXe+ISSUsQBKEPsNuc5O4u6LBN7u4C\n0obHdGmosLtrD06bNo2nn36arKwsAJ555hliYmJYvHhxq7Zvvvkm9fX1JCYmMnXqVIKCgrxe37Nn\nD2+88QazZs3qdEGySFqCIAh9wOnjFV5Dgm1xOhROH69gcHpUp/vr7tqD6enprFu3jqysLKqqqrDZ\nbIA7Wb766qvExsZSW1vL/Pnz2bRpE1lZWeTk5DBhwgS++OILCgsLKSsr46mnniIiIsJrMXNHeiRp\nffTRRxw7dgyLxcLs2bNJTEzklVdeITw8HB8fHx599FHef/99CgsLqa2tZcGCBfj5+XXaprGCcHON\npUic2o4LTwqCIPQUp9OJVtvx1229pWvfWV1t1921B3U6HXFxcZw5c4b169cza9YscnNzMRqNmEwm\n/Pz8+Prrr3nmmWdITU1l9uzZvPbaawBs3ryZ9957j+rqagCCgoLYvXt3l95XjySt5ORk7rvvPk6c\nOMHq1av57rvvmDt3LiNHjmTp0qUUFxezbds2/vznP1NQUMDbb79NZGRkp21efPHFVscqLy8HoCy+\n43IlgiAIPaWkpIS4uLgO2/j6dW1daVfbdXftQYA5c+awatUqamtrmT59Orm5uXz++eekp6czefJk\n1qxZ02ZsjftWFAWHw9FquLAj3ZK03n//fXbs2OF5vGTJEqqqqnj33Xd59NFHeeuttzCZTABERkZS\nUFBAaGgoACaTifLycmRZ7rRNW9LT0/n444+JiIhAo9F0x9sTBEG4IFFRnQ/nDUoJZ/vmkx0OEWp1\nGgalhHfpmN1dexDc382VlZWeZAiQlpbGBx98wLFjx0hOTuabb75pdZzs7GxeeuklKisreeyxx/jw\nww/ZsmULsizjcDiYPXt2uzH2SEWMI0eO8PHHH/Pkk08SFBTEW2+9xahRoxg5ciRPPfUUTz/9NM8+\n+yzvvPMOeXl5fPDBB0RERHTa5pe//GV3hy4IgtBj2ps92GjkuIHdOnuwP9Qe7JGkNWPGDEaPHo1W\nqyUxMZGbbrqJ3/3ud4SGhhIQEMCiRYv48MMPOXPmDLW1tSxcuBA/P79O2wwYMKC7QxcEQehRba3T\n0uo0PbJOSyQtQRAE4YKJihjtE0lLEARB6Dfk3g5AEARBELrqquxvVldX8+abb6LX6zGZTFRXVyPL\nMoWFhTz44IMcOHCAAwcOALB161bWrVvH73//e5xOJ5WVlSxdupSKigpWrlxJYGAggwYN4t577+1S\nmwuJw+Vy8dFHHxESEoLL5eLxxx+/7HF0FkNVVRWfffYZwcHBJCUlceedd/bIuTh16hT+/v7U1NTw\n+OOPU15eflHHuNxxaLVa3nnnHQ4cOMBf/vIXgG7/TFrGUFVVxR//+EdCQ0PR6XQ8/fTTvXIuSktL\nefvttz1rKZ944okePxchISGoqsrixYsZOnQoCxcu7JVzsWnTJr788ksSExMJCgrikUce6ZY4hKs0\naa1evZqgoCCcTidxcXHs3LmTFStWsH79enbs2MF9993HzJkz+fLLLxkxYgT5+flUVVXx4osvsnPn\nTj755BPy8vJYsmQJ0dHRPPDAA4wdO7bTNnfeeSd6vb7LcTidTm6++WbGjBnDf/3Xf3VLHJ3FkJeX\nxwMPPEBycjIPPfQQo0eP7vZzERQU5FkUuW3bNr744gsOHTp0wcfojjhmzJjBggULeOSRRwC6/TNp\nK4YJEybw7LPPEh4ezv33398jf59txXHDDTfw/PPPExISwrx583rlXMybN4+//OUvZGZm4nQ6e+1c\nBAYG4ufnh1arJTo6utviEK7SpJWfn8/kyZOZOHEiCxcuJDk5mf/5n//h1KlT/O///i8ANpuNjRs3\n8oc//IE9e/Z41ow1rhGrrKz0rL0ICgqioqKi0zZ1dXWetWZdjWPp0qV88cUXDBs2rEvHuNA4OotB\nURRWrVpFaGgo586d65YY2opj9OjRvPrqqwC4XK6LOkZ3xNG8LdAjn0nLGJKSklBVlffee48f/ehH\nvXYukpOTOXToEM899xzjx4/vlXOxc+dOjEYjSUlJ5OTk9Nq5uP3228nKyiI4OJgnnniChISEbolD\nuEqvaYWHNy3Os9lsFBUVsXz5cv77v/+bDz74AICNGzcyadIkAKKjoyktLQWgqKiI2NhYoqKiKCkp\nAfDU6eqsTUhIyAXF8cEHH7B8+XJefvllTp8+3S1xdBaDy+ViwYIFPPzwwxiNxh45FwaDgWuuuYYl\nS5YwaNAgEhISLuoY3RFHS939mbQVg91uZ9myZWRmZnL77bf32meyf/9+Bg4cyIoVK8jJycFkMvX4\nudiwYQOVlZWsWbOGnTt34nQ6e+Vc5OfneypK+Pj4EBMT0y1xCFfp7MHS0lJefvllTCYTJpOJs2fP\nEhwcTFlZGbfddhujR4/mpZdeYs6cOSQnJwPw2muvYbPZqKqqYunSpVRVVfGnP/2JwMBAUlJSmDt3\nbpfaXEgcAN988w0hISFUVlbyq1/9itdff/2yxtFZDOHh4bzyyisEBweTnZ3N5MmTu/1cREVFUVZW\nhs1mo76+nmXLllFYWHhRx7jccRw+fJh169bxzTffMH36dB599FH+9Kc/ddtn0lYMq1atYteuXaSk\npACwaNEiVq5c2ePnIjc3l88++wxfX18URWH58uWX/TPpLIbGMkTNq5H3xt/F8ePHeeedd4iNjcXP\nz49HHnmkW+IQrtKkJQiCIPRPV+XwoCAIgtA/iaQlCIIg9BsiaQmCIAj9hkhagiAIQr8hkpYgCILQ\nb4ikJQiCIPQbImkJQjNvvfUWf/vb3zyPly1bxubNm3svIEEQvIikJQjNPPTQQ/zzn/+krKyMnJwc\n6uvrPZVRBEHofWJxsSC0cPjwYVasWIHZbOaNN94gODi4t0MSBOG8q7JgriB05JprriE0NJTx48eL\nhCUIfYwYHhSENoSFhREWFtbbYQiC0IJIWoIgCEK/Ia5pCYIgCP2G6GkJgiAI/YZIWoIgCEK/IZKW\nIAiC0G+IpCUIgiD0GyJpCYIgCP2GSFqCIAhCvyGSliAIgtBv/D9ucyZVbYywQQAAAABJRU5ErkJg\ngg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Plot the block model. \n", "GeMpy.plot_section(geo_data, 13, block = sol[-1,0,:], direction='x', plot_data = True)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "geo_res = pn.read_csv('olaqases.vox')" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "geo_res = geo_res.iloc[9:]\n" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(array(['Murchison', 'EarlyGranite', 'SimpleMafic', 'SimpleBIF', 'SimpleMafic2', 'out'], dtype=object),\n", " array(['SimpleMafic2', 'SimpleBIF', 'SimpleMafic1', 'EarlyGranite'], dtype=object))" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "geo_res['nx 50'].unique(), geo_data.formations" ] }, { "cell_type": "code", "execution_count": 46, "metadata": {}, "outputs": [], "source": [ "ip_addresses = geo_data.interfaces[\"formation\"].unique()\n", "ip_dict = dict(zip(ip_addresses, range(1, len(ip_addresses) + 1)))\n", "ip_dict['Murchison'] = 0\n", "ip_dict['out'] = 0\n", "ip_dict['SimpleMafic'] = 4\n", "geo_res_num = geo_res['nx 50'].replace(ip_dict)\n" ] }, { "cell_type": "code", "execution_count": 47, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "9 0\n", "10 0\n", "11 0\n", "12 0\n", "13 0\n", "14 1\n", "15 1\n", "16 1\n", "17 1\n", "18 1\n", "19 1\n", "20 1\n", "21 1\n", "22 1\n", "23 1\n", "24 1\n", "25 1\n", "26 1\n", "27 1\n", "28 1\n", "29 1\n", "30 1\n", "31 1\n", "32 1\n", "33 1\n", "34 1\n", "35 1\n", "36 1\n", "37 1\n", "38 1\n", " ..\n", "124979 0\n", "124980 0\n", "124981 0\n", "124982 0\n", "124983 0\n", "124984 0\n", "124985 0\n", "124986 0\n", "124987 0\n", "124988 0\n", "124989 0\n", "124990 0\n", "124991 0\n", "124992 0\n", "124993 0\n", "124994 0\n", "124995 0\n", "124996 0\n", "124997 0\n", "124998 0\n", "124999 0\n", "125000 0\n", "125001 0\n", "125002 0\n", "125003 0\n", "125004 0\n", "125005 0\n", "125006 0\n", "125007 0\n", "125008 0\n", "Name: nx 50, dtype: int64" ] }, "execution_count": 47, "metadata": {}, "output_type": "execute_result" } ], "source": [ "geo_res_num" ] }, { "cell_type": "code", "execution_count": 48, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'EarlyGranite': 1,\n", " 'Murchison': 0,\n", " 'SimpleBIF': 3,\n", " 'SimpleMafic': 4,\n", " 'SimpleMafic1': 4,\n", " 'SimpleMafic2': 2,\n", " 'out': 0}" ] }, "execution_count": 48, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ip_dict" ] }, { "cell_type": "code", "execution_count": 49, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(125000, 125000)" ] }, "execution_count": 49, "metadata": {}, "output_type": "execute_result" } ], "source": [ "(geo_res_num.shape[0]), sol[-1,0,:].shape[0]" ] }, { "cell_type": "code", "execution_count": 50, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.0" ] }, "execution_count": 50, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sol[-1,0, :][7]" ] }, { "cell_type": "code", "execution_count": 52, "metadata": {}, "outputs": [ { "ename": "SyntaxError", "evalue": "invalid syntax (, line 1)", "output_type": "error", "traceback": [ "\u001b[0;36m File \u001b[0;32m\"\"\u001b[0;36m, line \u001b[0;32m1\u001b[0m\n\u001b[0;31m geo_res_num:\u001b[0m\n\u001b[0m ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n" ] } ], "source": [ "geo_res_num:" ] }, { "cell_type": "code", "execution_count": 53, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([0, 0, 0, ..., 0, 0, 0])" ] }, "execution_count": 53, "metadata": {}, "output_type": "execute_result" } ], "source": [ "geo_res_num.as_matrix().astype(int)" ] }, { "cell_type": "code", "execution_count": 76, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 76, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQwAAAEKCAYAAADn1WuOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAACe5JREFUeJzt3c9rVOcex/FPbsYSrUQ6Se7oHVr1RjdBYsDC/Ae9XcdF\nFwotLVIQXLiQOWgodZmCdFGCP2oNLRSyKd02YOmyBiuEMIvQNgV7dUx6PAchQnSc5OmiKLd6r/NJ\nyZmc8b5fG2Eyzvm6mDfPPOeZ2BVCCAIAw982ewAAnYNgALARDAA2ggHARjAA2ApZvOiDBw9Uq9U0\nMDCg7u7uLC4BIAOrq6uK41gHDhxQT0/PMz/PJBi1Wk1HjhzJ4qUBtMGXX36p119//ZnHMwnGwMCA\nJOnv//6nCs2XsrgEgAw0Cw399uovT97DT8skGI8/hhSaLxEMoAP9r60ENj0B2DJZYaBzTddnN3uE\nXPvXP0Y2e4RNxQoDgI1gALARDAA2ggHARjAA2AgGABvBAGAjGABsHNzCn/y3g0kc5sJjrDAA2AgG\nABvBAGAjGABsbHoCz/H//u3Up7HCAGAjGABsBAOAjWAAsLHpiZY4/YnHWGEAsBEMADaCAcBGMADY\nCAYAG8EAYCMYAGwEA4CNYACwcdITf8mLePqTr7K3xgoDgI1gALARDAA2ggHARjAA2AgGABvBAGAj\nGABsHNzChnkRD3Phz1oGI4SgEydOaGhoSCsrK2o2m0qSRFEUqVgstmNGADnR8iPJ5OSkhoeHtba2\npjRNVa1WNTo6qqmpqXbMByBHnhuMa9euqaenRwcPHpQklUqlJ3/GcZz9dABy5bkfSa5evaodO3Zo\nbm5Ot2/fVldXlySpXq+rXC63ZUAA+fHcYIyNjUmSZmZmdOPGDTUaDY2PjytNU0VR1JYBgSzwzdS/\nxrpLUqlUVKlUsp4FQM5xDgOAjWAAsBEMADaCAcBGMADYCAYAG8EAYCMYAGx8vR2ZevpE5WZ83Z1T\nnRuHFQYAG8EAYCMYAGzsYeCFwn5FtlhhALARDAA2ggHARjAA2Nj0RFtt9P9dwiZne7HCAGAjGABs\nBAOAjWAAsLHpiU3HxmXnYIUBwEYwANgIBgAbwQBgIxgAbAQDgI1gALARDAA2ggHARjAA2AgGABvB\nAGAjGABsBAOAjWAAsBEMADaCAcBGMADYCAYAW8vf6Tk/P68LFy6ov79fW7dulSQ1m00lSaIoilQs\nFjMfEkA+tFxhFAoFffDBBzpz5oxmZ2eVpqmq1apGR0c1NTXVjhkB5ETLYOzbt0+Li4s6fvy4KpWK\nSqWSJKlUKimO48wHBJAfLYMxNzenPXv26Pz587p+/bru3LkjSarX6yqXy5kPCCA/Wu5hrKys6MMP\nP9S2bdv02muvqa+vT+Pj40rTVFEUtWNGADnRMhiVSkWVSqUdswDIOW6rArARDAA2ggHARjAA2AgG\nABvBAGAjGABsBAOAjWAAsBEMADaCAcBGMADYCAYAG8EAYCMYAGwEA4CNYACwEQwANoIBwEYwANgI\nBgAbwQBgIxgAbAQDgI1gALARDAA2ggHARjAA2AgGABvBAGAjGABsBAOAjWAAsBEMADaCAcBGMADY\nCAYAG8EAYCMYAGwEA4CNYACwEQwAtkKrJywsLGhiYkLFYlFbtmxRoVBQs9lUkiSKokjFYrEdcwLI\nAWuFcfr0aY2NjWl+fl5pmqparWp0dFRTU1NZzwcgR1quMAYHBxVC0JUrV3To0CGtra1JkkqlkuI4\nznxAAPnRcoXRaDR09uxZDQ8P6/Dhw1paWpIk1et1lcvlzAcEkB8tVxiff/65bt26pW+//VaS9PLL\nL2t8fFxpmiqKoswHBJAfLYNx7NgxHTt2rB2zAMg5bqsCsBEMADaCAcBGMADYCAYAG8EAYCMYAGwE\nA4CNYACwEQwANoIBwEYwANgIBgAbwQBgIxgAbAQDgI1gALARDAA2ggHARjAA2AgGABvBAGAjGABs\nBAOAjWAAsBEMADaCAcBGMADYCAYAG8EAYCMYAGwEA4CNYACwEQwANoIBwEYwANgIBgAbwQBgIxgA\nbAQDgI1gALARDAC2QqsnLC8v69KlS6rVapqcnNS5c+fUbDaVJImiKFKxWGzHnAByoOUK49GjR3r/\n/fcVQtCvv/6qNE1VrVY1OjqqqampdswIICdarjD+cwVx9+5dlUolSVKpVFIcx9lNBiB31rWHsWvX\nLi0tLUmS6vW6yuVyJkMByKeWK4zZ2VlNT0/r5s2b+uKLL9Tb26vx8XGlaaooitoxI4CcaBmMkZER\njYyMqFqttmMeADnGbVUANoIBwEYwANgIBgAbwQBgIxgAbAQDgI1gALARDAA2ggHARjAA2AgGABvB\nAGAjGABsBAOAjWAAsBEMADaCAcBGMADYCAYAG8EAYCMYAGwEA4CNYACwEQwANoIBwEYwANgIBgAb\nwQBgIxgAbAQDgI1gALARDAA2ggHARjAA2AgGABvBAGAjGABsBAOAjWAAsBEMALbCev/Cjz/+qMuX\nL6u3t1d79+7VkSNHspgLQA6te4Vx+fJlnTx5UmNjY/ruu+/UaDSymAtADq07GEmSaOfOnZKkHTt2\n6P79+xs+FIB8Wncwdu7cqcXFRUnSvXv39Morr2z4UADyad17GO+++64+/vhj9fb26o033lBXV9cz\nz1ldXZUkNQt8XAE6yeP37OP38NPWHYzBwUF99NFHz31OHMeSpN9e/WW9Lw8gB+I41u7du595vCuE\nEDb6Yg8ePFCtVtPAwIC6u7s3+uUBZGR1dVVxHOvAgQPq6el55ueZBAPAi4mDWwBsBAOAjWAAsK37\nLomrU4+QLy8v69KlS6rVapqcnNS5c+fUbDaVJImiKFKxWNzsEZ9rYWFBExMTKhaL2rJliwqFQkfN\nPz8/rwsXLqi/v19bt26VpI6aP4SgEydOaGhoSCsrKx01uyVk5NSpU6Fer4cQQnjvvffCw4cPs7rU\nhkqSJCwvL4e333473Lx5M5w+fTqEEML3338fJiYmNnm61n7++ecQx3EIIYR33nmn4+b/6aefQpIk\nYW1tLRw9erTj5v/ss8/CxYsXwyeffNJxszsy+0jSqUfIi8Witm/fLkm6e/euSqWSJKlUKj05X5Jn\ng4OD6uvr05UrV3To0KGOm3/fvn1aXFzU8ePHValUOmr+a9euqaenRwcPHpSkjprdlVkwXoQj5Lt2\n7dLS0pIkqV6vq1wub/JErTUaDZ09e1bDw8M6fPhwx80/NzenPXv26Pz587p+/bru3LkjqTPmv3r1\nqpIk0ddff62ZmRn98MMPkjpjdldm5zAWFhZ08eJF9fb2av/+/XrrrbeyuMyGm52d1fT0tL755hu9\n+eabTx5P01RRFOU+fJ9++qlmZma0f/9+SX8cxOnu7u6Y+WdmZvTVV19p27ZtWl1dVV9fnx4+fNgx\n80t//Btu3LihRqPRcbO3wsEtADZuqwKwEQwANoIBwEYwANgIBgAbwQBgIxgAbAQDgO13DNtegcvc\nzsEAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.imshow( geo_res_num.as_matrix().reshape(50, 50, 50)[:, 23, :], origin=\"bottom\", cmap=\"viridis\" )" ] }, { "cell_type": "code", "execution_count": 77, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 77, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQwAAAEKCAYAAADn1WuOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAChFJREFUeJzt3M9r3HUex/FXNlNJq6Q4aXaaHdS6aS+hpIEK8x+4ntOD\nh8oqShEKPXiQ+VKD2GMKpQcJ/WFtsCDkIl4NVDzaoRZCGJagRqir08Rvv1+EFNKOk3z24LbstruZ\nVyTfyXfa5+MSmExn3gmdJ5/5fD+TnhBCEAAY/rTdAwDoHgQDgI1gALARDAA2ggHAVsjiQe/evat6\nva7BwUH19vZm8RQAMrC2tqY4jnXw4EH19fU98v1MglGv13X06NEsHhpAB3z66ad66aWXHrk9k2AM\nDg5Kkv78z7+q0Hoqi6fAH3Cl9o/tHuGx9PfKyHaPsGVahaZ+ee6HB6/hh2USjPtvQwqtpwhGjpSH\ntnuCx9Pj+H/8/20lsOkJwEYwANgIBgAbwQBgIxgAbAQDgI1gALARDAA2ggHARjAA2DI5Go58+ttf\nxh65bbYxtw2TdK//9Tt8krDCAGAjGABsBAOAjWAAsBEMADaCAcBGMADYCAYAG8EAYOOk5xPu4ZOL\nnPzERlhhALARDAA2ggHARjAA2Nj0BDbwpH+c/WGsMADYCAYAG8EAYCMYAGxseuK/8Hc/sRFWGABs\nBAOAjWAAsLGHgbbY18B9rDAA2AgGABvBAGAjGABsbHoC/8YnU9trG4wQgk6cOKGRkRGtrq6q1Wop\nSRJFUaRisdiJGQHkRNu3JNPT0xodHdX6+rrSNFW1WtX4+LhmZmY6MR+AHNkwGNeuXVNfX58OHTok\nSSqVSg++xnGc/XQAcmXDtyRXr17V7t27NT8/r59//lk9PT2SpEajoXK53JEBAeTHhsGYmJiQJNVq\nNd24cUPNZlOTk5NK01RRFHVkQOQTpz+fTNZVkkqlokqlkvUsAHKOcxgAbAQDgI1gALARDAA2ggHA\nRjAA2AgGABvBAGDj4+3YMt10+pOPsv8xrDAA2AgGABvBAGAjGABsBAOAjWAAsBEMADaCAcDGwS1k\n6uEDUttxkItDWluHFQYAG8EAYCMYAGwEA4CNTU88VtjgzBYrDAA2ggHARjAA2AgGABubnuiorf4z\nfmxydhYrDAA2ggHARjAA2AgGABubnth2bFx2D1YYAGwEA4CNYACwEQwANoIBwEYwANgIBgAbwQBg\nIxgAbG1Pei4sLOj8+fPas2ePdu7cKUlqtVpKkkRRFKlYLGY+JIB8aLvCKBQKev/99/Xee+9pbm5O\naZqqWq1qfHxcMzMznZgRQE60Dcb+/fu1tLSk48ePq1KpqFQqSZJKpZLiOM58QAD50TYY8/Pz2rdv\nn86dO6fr16/r1q1bkqRGo6FyuZz5gADyo+0exurqqj744APt2rVLzz//vAYGBjQ5Oak0TRVFUSdm\nBJATbYNRqVRUqVQ6MQuAnOOyKgAbwQBgIxgAbAQDgI1gALARDAA2ggHARjAA2AgGABvBAGAjGABs\nBAOAjWAAsBEMADaCAcBGMADYCAYAG8EAYCMYAGwEA4CNYACwEQwANoIBwEYwANgIBgAbwQBgIxgA\nbAQDgI1gALARDAA2ggHARjAA2AgGABvBAGAjGABsBAOAjWAAsBEMADaCAcBGMADYCAYAG8EAYCu0\nu8Pi4qKmpqZULBa1Y8cOFQoFtVotJUmiKIpULBY7MSeAHLBWGCdPntTExIQWFhaUpqmq1arGx8c1\nMzOT9XwAcqTtCmN4eFghBF2+fFmHDx/W+vq6JKlUKimO48wHBJAfbVcYzWZTp06d0ujoqI4cOaLl\n5WVJUqPRULlcznxAAPnRdoXxySef6KefftKXX34pSXr66ac1OTmpNE0VRVHmAwLIj7bBOHbsmI4d\nO9aJWQDkHJdVAdgIBgAbwQBgIxgAbAQDgI1gALARDAA2ggHARjAA2AgGABvBAGAjGABsBAOAjWAA\nsBEMADaCAcBGMADYCAYAG8EAYCMYAGwEA4CNYACwEQwANoIBwEYwANgIBgAbwQBgIxgAbAQDgI1g\nALARDAA2ggHARjAA2AgGABvBAGAjGABsBAOAjWAAsBEMADaCAcBGMADYCAYAW6HdHVZWVnTx4kXV\n63VNT0/rzJkzarVaSpJEURSpWCx2Yk4AOdB2hfHbb7/p7bffVghBP/74o9I0VbVa1fj4uGZmZjox\nI4CcaLvC+M8VxO3bt1UqlSRJpVJJcRxnNxmA3NnUHsbQ0JCWl5clSY1GQ+VyOZOhAORT2xXG3Nyc\nZmdndfPmTV25ckX9/f2anJxUmqaKoqgTMwLIibbBGBsb09jYmKrVaifmAZBjXFYFYCMYAGwEA4CN\nYACwEQwANoIBwEYwANgIBgAbwQBgIxgAbAQDgI1gALARDAA2ggHARjAA2AgGABvBAGAjGABsBAOA\njWAAsBEMADaCAcBGMADYCAYAG8EAYCMYAGwEA4CNYACwEQwANoIBwEYwANgIBgAbwQBgIxgAbAQD\ngI1gALARDAA2ggHARjAA2AgGABvBAGArbPYffPvtt7p06ZL6+/v14osv6ujRo1nMBSCHNr3CuHTp\nkt555x1NTEzoq6++UrPZzGIuADm06WAkSaK9e/dKknbv3q07d+5s+VAA8mnTwdi7d6+WlpYkSb/+\n+queffbZLR8KQD5teg/jzTff1NmzZ9Xf36+XX35ZPT09j9xnbW1NktQq8HYF6Cb3X7P3X8MP23Qw\nhoeHdfr06Q3vE8exJOmX537Y7MMDyIE4jvXCCy88cntPCCFs9ZPdvXtX9Xpdg4OD6u3t3eqHB5CR\ntbU1xXGsgwcPqq+v75HvZxIMAI8nDm4BsBEMADaCAcC26askrm49Qr6ysqKLFy+qXq9renpaZ86c\nUavVUpIkiqJIxWJxu0fc0OLioqamplQsFrVjxw4VCoWumn9hYUHnz5/Xnj17tHPnTknqqvlDCDpx\n4oRGRka0urraVbNbQkbefffd0Gg0QgghvPXWW+HevXtZPdWWSpIkrKyshNdffz3cvHkznDx5MoQQ\nwtdffx2mpqa2ebr2vv/++xDHcQghhDfeeKPr5v/uu+9CkiRhfX09vPbaa103/8cffxwuXLgQPvzw\nw66b3ZHZW5JuPUJeLBb1zDPPSJJu376tUqkkSSqVSg/Ol+TZ8PCwBgYGdPnyZR0+fLjr5t+/f7+W\nlpZ0/PhxVSqVrpr/2rVr6uvr06FDhySpq2Z3ZRaMx+EI+dDQkJaXlyVJjUZD5XJ5mydqr9ls6tSp\nUxodHdWRI0e6bv75+Xnt27dP586d0/Xr13Xr1i1J3TH/1atXlSSJPv/8c9VqNX3zzTeSumN2V2bn\nMBYXF3XhwgX19/frwIEDevXVV7N4mi03Nzen2dlZffHFF3rllVce3J6mqaIoyn34PvroI9VqNR04\ncEDS7wdxent7u2b+Wq2mzz77TLt27dLa2poGBgZ07969rplf+v1nuHHjhprNZtfN3g4HtwDYuKwK\nwEYwANgIBgAbwQBgIxgAbAQDgI1gALARDAC2fwEeNmZOEszfygAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.imshow( sol[-1,0,:].reshape(50, 50, 50)[:, 23, :].T, origin=\"bottom\", cmap=\"viridis\" )" ] }, { "cell_type": "code", "execution_count": 59, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 59, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa0AAACmCAYAAABp5tBAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl8XFX9+P/XvXfuLJlM9snarE3TJelKaSmlULpQVKRU\nkfoR5FNRP2j5wlcUtf3w/X3cABFR/Ig/5fFhsYJ+Pyh+QFEEpOx0o9CF7nubpNm3ySSZ9c75/jHJ\nJJNMNpqkSXuePnjYuXPm3jt3MvO+55z3OUcRQggkSZIkaQJQz/cJSJIkSdJQyaAlSZIkTRgyaEmS\nJEkThgxakiRJ0oQhg5YkSZI0YVxwQSsYDFJZWUkwGDzfpyJJkiSNMNP5PoGRVlNTw/Lly8k+NQ1T\n0Hy+T0eSaBb17OZdQoSitiso2EmgDRfTmMckpYiQMDjNEZqoxUUTgr4jUoqYQZEyY6xOXxoBr4We\nO9+ncMG44GpakjTeJCtOyljYZ7tA0IEbgBYaAFAVjXRyCCFiBqxciilk+uiesCSNYxdcTUuShqpd\nuLErjjE5VrqSwzQxl8PsjtreFZYaqUUIgaIoBAkgetXKADLJo4TZKIoyBmcsDeT4I5fF3F589/Yx\nPpOLjwxaFzBDGPjowIsHLx4EBtkUxvzRU3UTSblZ6DYLgUCAVgHCbEJ4fIQqaiAwtD5CsyYoTIM4\nM3T44VQD+I3u40Udx+OjpaKa0BD3PVyaGiAt2YNNVzCH7Bxu8BIwFCyKjXbRyikOU8aCYe1zoPMf\n7L1PUibjE15OcSiyrSs4BfBRxRlaRD3VnOlz3FQymcH88xKwgiKAFw8mTFiVuBHfv094MWPp972N\n5d+MNP5dHEFLN6HmZqLYLMP+EW4SdVSrFZgdyZgt8ZhCJkwtPvSgho658z8LmqINuJ+gCHCSg6ho\nmNAxoaOjYzJZ0TMy0eMcmPwKalUzakAM+cfJEAb1nKWdNqo4jUGQ8E+h0acPJVXJokU0ERIGIQwC\n+AngB7PCDdd8i+yMYtA1MJkIAbXuVs66GtDnTEY5XoVo8wx4DefmCubkCvTOS1He7iIlLcCJepVj\ntQopRbmkTs5F03VUXUdXk7CWOmk6cprWI9WEMGihgVaayaOEBJO1OwgIqMg0E4xXEW6B97AfX5IN\nj9qOt72F9uYqPCE3Xtrx4cXAj2GEUBog2WLF5fehqxq3ZC3jRLXOR2wnQU1FLZwUfk/+ICig6CZE\nIAgCFLMJYXixJjVQGG8QH18I8UX4Q1q4hqQoTLpiLvXV9VD+Djur3qPJsJBkWEiwm3GmWGk2bFSc\n1LEKO2k5uSy0fwoqNE6595NODnWcjVy/Q3zQ42oqdNXDLJY45i+7AWoSEK3+yLWPBEmbiicllYqk\nfOISEjBpELK10uqvIXimkcAeL/j6/u2ETAr+rHjaaaWi/gNaPFVAuNlSKAIhQpEmSrNiYsq0y9BT\nBCKk4T/rJ0MtxFTlAkAtyEHNcYKqhK+fqx3R1kGopgE1M63P380ZcZRqyrFgBcBFI1YRh4MkEkkj\nBSdW4sgqKyGzdDI2q4aqQEiAxzuDmgMnqD1wfMDvRvZ2ByKgEahzInxmFIsfPb2edz+Y1u9rLEGD\n0mYXjkAAt65zIDkRn0nDFFLI9sZhCWn4VANz0IPFb8MUUgmqIdqtHYRUOSPeWLjgg5ZWOhmtdDKY\negSV+TMwDpzAOHBi0NenlV7G2ZZWzpzdNWC5BJHMdC7BoSTFfN6k6MSLJA6yM/qJIPT43cKkx1GS\nfQVZVY4hBS5N0YgXidRxFh8dA5ZtFNUxtxdmzmXfqXfYcfDvdHhcdHS46PC04vd7mDPvM2TNWY26\nJBPh8YG389ev1zWcmyu4tCD6SxsSgp8d2kaHEQhvONn5Xy9JiVn4LR10+FwoKMxmMTOyDYg7TWvI\nT5Xixy38tNb4qQ0FaDaC+BQv/iovIiQgRt9PFwE0+bwArJ5azNJcG1tdu2jvaCW9eCGmS2eC1YJi\ns4RfYBighf9WtOQzzA4epbQqhNtUSJOtJLxPQ9Bh6HjMcWgKFDl0CtNnEtpxhD8cO9z/B3Am3Gel\nmyyAEhWwYp85qKrG3FUrifuEDSUUIliRSfDUDMrKdzFbrUOPt6DE2WhOyCLN5OCD07swrBqzSq4i\npEyieW41jZ85g+/lFvxv+To/F4Mz6XWcathBqLyfm7del9Qvghw49F7kcdFVc3FcsQKjKhPMOopu\nAlUN/4cIR5dAAFQV4fX3+bvJ3R+igWoaqYnssw0XbbiiapqpFVksTbmW6fmziLMlcOjMHtISs5g8\nP/xZDBS4fCfz8J3KQxjd333v4WKSQkFa4t19yl9ZVctV1XXooRBCCPZVVzMn3k7t1CsJ6TmYRDgF\nQDNUlJBGrRqkVjPC19SdRLO9td9zkUaOcqFNmFtZWRnJHrRMnY42u6Tfssbeo/0GrqAI4M630Jro\no7HmEE11R2OWM2OhiFKyKUBVYue1+ISXFhpooaFHbagXRSF/ytVMmbUasyV+wHPrIoTAQzsuGnHR\nRBO1dNA24GuGymKO41Or7qaw8BKOis46mxCI5ugvprH3KNrh49y8sLuG1aUt4OfP5Yf4x9nBbw7C\nFGaykGvzcpibZ/D/H/mALfWVsc8vIY4FX/kU5o4CDv9zC2dPbYNh/imnZkyncPanSMyZgsWegtL1\noxsKoTkrmR04yOxTLgxMlCcuJaRo4R9zBRAKbkPjTONZ3DUHqKw5Tk3Vfirb+v4YAiQ4Urh09qfY\n/dFmmlpj3zwMRNU0HFkpJBdkMQ0nlysW8uIcaBYzb7d42FJ+gqPlH5GalMFXbvkBRnx3M15j4lma\n9DN4/xYOXFrpZJSZk2moPsCxfX/D1RjjTmKIpq24joJpt6JpJuh5k9X1TyMEQhB0u2mrPUVr0xlc\nTWdwVR7B7akjhDHkY8VZ40lLTGf96ntRFIV2j8H2/36dUIzhLVrpZKxrYyestLVbaYx3RQWuK6tq\nWXE2HEBPNDTwpz278AaCfPOGr9NmL6ZV12kz62iGiinU/Yde1SNwAZj+vjvm91ZmD46cC7empWvh\nGtYAtNLJGEdOQ9DAEEFqqMBFEy4aaaeVGF0LEaqmUzh1JblHrZiMvsGqTpylnmpcNAwaSFIzpzN9\n3udJSJ4U89x6qxJnqKMSF00EYrX7xJCUNhmfpwVPe+OQypstcRw+9h6FBXNJUBRahAj/KFl08AWi\nzrOw+SR6jy/uSXczvz3xEUdbGweoA/W18PLbSN3lZk5uEJOqUprsZFvDWUK9gtG0dCdT/vcN6KkO\nMASzO9YxufRTHN/3IlWn32egmldPjbWHaPxnuH9J060kZRTjSCsgNW8axVNVSne6AYV2PSMcsID9\nJz7kZNVhymtOUlVfjhEauJnZarZx5exVeDXBweNbBghY3c2BsYQMA1dlPa7Kek4Dr3Rut2gmfEb4\nHEyazi2fuIMEdFoEiM7AkdyaSXNmDfub38erevBX7aL1UCUhwz+k6zSQw5v/zuHNL2FLcJKcNRVH\naj65ZSvwdTTRdPYQrrqTtNaewN1UgQgNHqAs2AgoHkIxLkWHt43Vn16PNS5cu7baFR551ctP7MlR\n5URAw/1OCe5Wa7/HSW5PwBXXhlAFlqDBVdV1eAMBntqxne1nTgPw3RWraI8rBMARCNKum6ICFkBm\nyES9ahDqvNYDfW+lkXHBBi0lyxndJBiLSUPNyyJ0shJQOMyuqDRjVTOTmJJPsnMydWc/os0VbvPP\nKVxEyezPYLOnEPTt63x9tHqqopo5TOgkkYaGiVoqAIiLdzJ93lrSJ83p2xQYdW7RWmmige4fPws2\nEkkhkVQCSWZOt4T7RpLTS0CEaK4/TmbuPIpmXEsw4OPEgZc4fXgzhtF/wHO7GzhbdRjNZMJkhIj8\noPauTZo07JmJQFNkk91k5khrY+e/dcy6lWZP+K42weHE7/fg9UUH8mVLvkRuyTIcHW+ja+H073YC\nfQLW9dOn8YW5s/ggqIVbGrUQWlor8UYmcxb/G9ddMoOP3n+N7RWVaIpCktVKo8eDWVMxqRodgQCx\nGAEvjZX7aazcjz9QzbKmhZiM8LENzRIpt+PA2xwt3xd5rKCQkZJNXnohMxPtbDnyPodbmtFVlVUl\ns1mw+EuYTRb+46m78Af6Xm+zOQ6ns4jpxYuAEO/t+CNeb/e10TQTScmZNDVWIUTfjMKugAXgTM4k\nIyUbhMAcDOHTw5+VGtJw+FIRmkKdUdnzoyLOkU5iSiHJBemoKS3sf+HdyHOqovS5/rEJPK11eFrr\ngHeZVLqUU7tfovLA673KKdgTMkhIzicxJQ9zVQcf1b4EgA07Rcwgkzy+8Mi7fPkbdX2OkhSfTHpS\nRvdRQ9B4KhXKoq9LoM4Z1SQYiyoU4r023HEdlDa70EMhdF2nLCuLHWdOMzM7m4K8ObR03qwoCOwB\nQaDX11QFkoRGk9IZpAb43kojY0IFraNHj/LEE0+QkJBAYWEhN998c79lFevQBhZ39WVoioZTZKOi\nkkgqyZNnkbRgMapqQogQFcffJSV9KtPn3URiakGf1/eWRhYCQRKpJJGGnQQUReGkOISmmime9WkK\npq1E0/RBz603J9koKCSRSiKpURldR9QTJKYUMP2StaSkh5tGm+qO0lgTrlGYdAtT53yGwszLqarf\nxulj75LgcFJTd4LZZdfgdBaQN6mMqurDkQAe7FkDiPHD6dGizzPDZueWwjJKElKZkpDMQ8f20err\nwBGfSktr3x+jxQs/zyVzruNsewd2W3dQvCo3H8Wm8Pu9+7CZTNxx2QIuK8gDwObvvpNVzN0/3PGZ\n6dyzZDGnmpr4474DrJxciK5q/Hzrdtr80TWL4jk34PW58LY1kZhZjKbptDacIa04L2r/Wqg72Ewr\nmIlZN5ObUURuehHZziIcuhkQTGraw592vcHySbl8eXopZM6m2hoPwMJZyzAsFiZnzubZFx8kEPQx\nu3QlS6/4VwzFjK2zbXXWnKXs2f8qW7e8QHHxPNzuJlbfcCcWq50a20k4tgMOnuVYfTNJFjvZWTPY\nWV1BRe0J5k5dFDnPqJwABUyGTu68EpqO1dPe2hx+XyYLzqwy8kuuJqlMp751GwkOO4uzs/moqobp\nGU7OtrZS627jfy1ewF/2H+ajmvDnZ7VaUOLMlFy7kKr3K0nLXkSHq4YOVy0WewpJmVNw1Z4gMb2I\nBGchzpx8cvPTSE3qvq7x7mepf1zj3rtT+PK/JGA2NwPNQAIHd+Tzs/8O9/86kzKob6mlpa2Z+353\nL8suvY4r5lyDbtJ5uyWFd7Zl0FNym4OU9v5rWd2fa/iaO3rcyFxRNBmbrpPhSCCkRP+GaKJv0ALQ\ne8X1/r630siYUEHriSee4O677yYrK4uvfOUrfO5zn8Nsjh2chHdoTR/C0/2DNEvp8aUnBVUNXx6/\n182sy74Us0bU8/U9ZSiTyGBSn+1W4rgi91+xl8Ye5zGUfacqGaSSEfO5TOsUypauROlRI0pJL4kE\nsC6aN0TRrDUsnvMpXt78K9au+T75+bMjz6el5oICIaC1625biKimwS4nmzUuSyeqT+v63PDxWgM+\naly1hEQoErCmlyzBarGze98rzJ97PYsuvZGQgFafn3ZPCDrjeLJuZUa6k7zEBO654nKyExMi+/eY\nuw8m/N1/xqf1FBaopylMTmLDlYvxBQ3+efxEn4AFcHzPX8grXcmC1f+nuz8L0JLr8Kh7I+Xs/loa\nbdMJKRqLZ69k8eyVnQdW8AsdhB9VGJhay/np4iVMTUoGBA2GN7KPq69ai8+k0lJRT5zVwbVLvkVe\n8SUAuHx+bITfj0mzcOmCTzBr1lV4fG3YrA4sFlv4vOJNOEvyWKgkIPxmCJlx2zKZXPYpBAa6Hr4O\nTa56KjrqSc8vwWKJAwFBLUDGtFyWLP4Xtm7/G63NFRhBH2eOvsGZo2+Quncy2Zfl8q//ay1XVXRg\nBIMca24iIAzKMtMJGAbfy1nK9tNneeL9D5m6sIzsW5YBkD/1cxgt6d2tmwrkz1pF/qxV4cdGCKvV\nwGp1QY8+rJLJJo5uzScurm/z+re+oHGi8jP85d3nuePGezheeYS/vfc/uNqaeXnLn8lOncT0glL8\nqXXQEv1dCKp9b6xiMdTwubj17htHVVGYn5sX7jMT0X8zRj+JUb0DWX/fW2lkTKig1djYSGZmJgCJ\niYm0tbWRkpISs6yorg+3Kw/URBg0CJXH7mMIVdTA/Blg0rDYEsnInTus1/cnW8mHavc5ndtAEhtU\nFEMM/MkGDYzdh9EmZeCPT2DNp76D1RqPCAkUtcc3UEB9yIgkzsf8MgYNfKdr2eNTorIHg6EQ/6w+\nyXNnDtEeDAe6zPTJLLvyNnKyprF73yvMKl3B0sW3oigK9R1eQkGDk4dbWHRJOAAKjyA10cYD1yzH\najJFEi2CmkJ5auedtKFiNHQHM29dAvtTdWbWeSl3uThQV8eWM+XEm3Xa/NEB15kzm7Krvhy+ERFd\nWYgKhiuZ8snxXKK1YDIEmgiS5D1Jk21KdyIGIIRKm2LBIgJktZ4gxayRYu7uX0n2VFOZWEZQ0/Br\n4TsALWDiXz/5Q8xmWzjJDqjp6MBhNqMq4f4YxRLAYo3DYouLBIKQatAW14zHauWSEypaQAUV7N56\ndMdkhGqhKyAkJTr5575X+e//eZj09Hxycqdgnm0lKTme4LuCWYu/xJbX7o/qY2qsOEFjxQmqCg5T\ntmARqZqJaSmpKObwm7WYTIDCwrxcSifl8GRysOvEMFy9voM9rhFCgBAoisBm6f77sZgC3LisHqs5\ndvJSWlkTt116Gd7A9VjNNi4pWURp3iW8vedVKurPMH3STLJdRwi0972Jard2EHInoYr+s29DiqDN\n6gHgQHIinyo/ix4K/6V33Zja/LW44koQioZAoV1XMPfqwgwBLUqP/quP+b2Vhm5CBa3MzExqamrI\nysqipaWF5OTk/gsHDIwDJwbOHjxwov8O00Dw3F4/kPGwb48P48AJ6meXgNWKWYBKCBFSUVSFEFDv\n9VEfCETG2ETSlmOc5+6K8Bd9Tq7ggKuW3538iLMd4X4su2olP1TGvNlrycgsAqAofx6zS1ciCAes\neq83nELvC7GnojMACkgxrCimHj9sIcGBfAfBzm3BijTo0X9hfHSaXU6BqzDIG2/t5UBt3+bIMIWZ\nC28FXxDiTJF9oyoQ1PA05XAgz8XsU+FxSMnekyAELbYiQphAKHgUnSAKrg43M9oqECgoSo8+UREk\npeMEp1KmIRSFkDtEnKVzBg4hwOOjAUFQQIPHS3pcOBALvync5BkSkYy85oQaQmqIkKqyz5bFnNYm\nhM+HatZJai+nOb4QUEEJ4bOYWHntOnx+DwcPbKW29jRdQ8Bs2EneUUh69ixqK6Nn58ibNZ95ty+m\n8nQrqZXhfjVVEYS6fvyFgoLCycnJZOQnIBAYjdmoaChqCNEjSIjOqGXSQmjmEA5bO2qPa7N85k6s\n5titIV8qXwLA1Kn13OiZTotqJmQomHULKy+9HsUIkuE+hs1XTs5WnStzanknu7u2FVIFzfZWUtsS\n+/nsodneiuhsQ/WZNN7OSo9kD3Z/fgYOz2la4ybj1k0IRSGoGlHJGDVqMJKEAefwvZWGTPv+97//\n/fN9EkNVUFDAr371K3bs2MGCBQsoKyvrU6a1tZWnn34aR0saSm0rhARqWlKk6QcI1zT2HR88pby+\n+ZxeP9733VXOkxRPk9+PPxTCGwjS0uGh4vgZ3IdOEjpegahpQE2IH/Q8T7k6+L+VH/BK3UHcAT8q\nKoVKCTPEZSQqqbSdPEuTz0cgOQFhsuEKBKlqa6fd64vaV02rQkhAugPUzhtpxawQ1BT2pZs5kGUH\nRSdYnk6wwtnnfIzTBq1JNtKuKaT+SAWeXjXEZJxMy72GpLzpKCIcHBWTFg4kgQAoCqLdQZ1Vg/g2\n0tr9qIbAFmgk0XMG1d1GTWOQk2o8VR0dVId0fIoJh+FFEMTQwacpuKwqh9M6qIzTsDTF0zOJNBQ0\nqNl9mLqztahpSXSEQiDAZjKhGCoEAihqiJAaoimpmqakGjBUguVOqmpyCVbV4/S1oCoCW6gdFIHP\nnIhLUfBaVVRVpahkDmcbj+FqqO/+2Ajg9tbR3lqDqulRyR2uumqUKp1gWSqYFNI8IWx6iHC1ScXQ\nVPYVJHAwPwGEis2dSqI3oTMYKYASDlwK6JqB3eJFUQQOWwcJtnYgXMNaNXsHK2b2Gq/Yw19d+QA0\nZmrEGU3M238ca6Adm7+V5I6zTHLtQxEtkWa93PYOdmSkYfT4+/Sa/QgE1oAFhe6oElIETfGtfcZp\nnXHEYygKue0daD2ST9RQK0cSHXjMSagokUAnhEpNz3Fanfu1/fEDYrn1ezf1+36l4bmgx2lFZnnv\nzOiJjMovrx7e3dC5vn6873so5YZQpk242M5rQDhZZAqziFPiP/Z56ZqgqJ8ZMQKHDdSM2PsIiRCV\nnOAkBwkS3Xxkxc5lrMSkmKLPo7PpUDHr0f8OerEmNlDoMLD5QrhPG5yshoCh9Hkf2tkqipIM4uLA\nl2XijFnD74bAfj+qYSIpLxPdZiXg8dJSXtM9vqjHfhR/gARAN5sJBL20JzQiHIA1CaPKhr+5FbWy\nCcUIdV8fm4onOZWK5Dw0i8aR45uZs+BqWgO1+E7UsXPXZppD9VHXoSDXxH/9PIfrbq6goDCJhHiV\nDz4MZ21+/jMO1v/7dPxBK/ZyhRw8CKvCu1ou9QEHfp+NkvhWMhJamZ0XHty78+Q0Dp0tJBhUibP4\nyEhsJsXhYmrWGY5U5+P22HHY2pmddzxSw+qqUfX2zrbSyL/n1Tex5lQ5NsNAEwJDUfBoGqJX/9Jf\nCnLZ5ezbVaCEwlmCWkjDUA3arJ5I4InFYhiUNrmIDwRo03UOpCTi08IzYmR547CGNLyqQY3Zg9Vv\n7bPf/uYelOO0Rs7FEbSkMXNaHCGBJFKU2Ikio61R1HCUvbR3zp5uxooVG62EM+au4JOjMn/eWKkX\n4WEXTiU75vNCCN7geeZMvoKsvDyEW+DZ18Fu3zu4eua6A6nJKksuszGtWOcH30nlnu838OiT4ebQ\nBXMt/M9TWWRnjl4PwlCC1lVVtSzv1WwXy+aczKgmwvNFBq3RJ5cmkUZUgTL1vASsDuFmj9jCbt6j\nHTcKKgVM5XJWReZgzGDShA5YAB20cYpD9HevabnaitluYV/ldvyzvVg/YyNpYyoLrlhGPOE+npVX\n2XCmajQ2h/jLy+08/4923tnu4Rf3OXn85+noOry/28eCayvYscsb8zhjpWdmXyy+YJBQKETbIOWk\nC4cMWtKEFhQBjomP2MY/IwOunWSziGsoVmYSJEAb4drDJAaeIWUi8NBOK8000TfBxLzUgmWlFbPd\nStAXYPf/vEZ+cStT53uYt8HMb1+bzpQinc+vcXBkax7f/FoSJhMcPRFg5eeq+Oxt1Sy93Mbrf84h\nPU2jutbg6s+c5ek/DW9OPa+h8UGDkzeqs/mgwYl3kIG+AzmQnEhA7f9nSgEeeedt/vnii7S+t5VA\nfX2/AV26MEyo7EFJ6iKEoIrTnGA//s6prOwkMJXZUTW9rllJbNhJIu28nOtI8hBOaDjFoeixehYw\nXxXOPrTYbbhpobGilbfv28YtX72cxkkae9V8/vbftTQ1BEhM0Pjp99L4ys0J3PODBv6xuYO/vNzO\nP15v55u3J/PG8zncekctu/b5+NL/rmPfIT8P/p9UNG3gSZzfrM7mrZosAj0y7P5WkcfSzGquzqoa\n9vvtL7Ovi9lkYupnbuDgH/9E2779AGjJSdhKSrCWTMFWUozmGJs106SxIWta0oTUhotDfIgfHzpm\npjKXhayIClhCiEjQyiJ/Qi2e2C7ctAlXn+1dQauFBppFd3KFXmamawIHp6m7qWz3R5UcePoQC1/2\nk3lQ4IrLYOG87tkiphab+dsz2bz0hyymFev4/fDgo82s/NxZvnprAmtXhxNpfv5YCw/8Z/OA5/xm\ndTavVU2KClgAgZDGa1WTeLM6dj/cYN7JzmBzTmafGldAVdmck8n7kwvI+NpX0BLC4/WM5hbadrxP\nwzN/oOL/+wFnH/oZTX95EX+1HD91IZA1LWlCcihJ5IhCVDSKmIGu9E26cdEYmaw4i/yxPsVhE0LQ\nSC0VHMNDBwtZ0ed5b2fQAjjFYZIJp/wrjnBALqvz4UHvscwkbPpgH8VpScwKptOR5IBPNfQ59rXL\n7CxfEsdvNrn4/sNNVNcafP3b9bz5fDazZpj5/Z/d3PWV/sc9eQ2Nt2qyBnx/b9VksSi9Bqs2tBkr\nenonO4MdGWkxM/sATCkppP/bl6l59NcIX/QQh2B9PercOehO57CPK40/sqYlTVjTmMdUZU7MgAVQ\nTTkQHpdlU+yjfj4fty8lKAJUiONs41X28B6N1DKDS/osLOrDG7WwZxO1uEQ4I1C4BbohKK33k2CJ\nnvsuJOAnb75PQ3sH5r1mDF/sGqeuK9z11SSObMnn9lsTWH2tnSsXxbHhrhQ++GcuiQn9903tb07p\nU8PqLRDS2N8cewabofBpGrucKbyTncEuZ0okYHWxTMoh/cvrIuuhdVF0Hc3Ra5yhNGHJmpY0YQ3U\n3GfWBJ9MKaMimIwurBhugd8YveZBIQRnOcUkiob8mg7RRgXH+6yxlksxSUrf/jdPjCVuTnOY2VxO\nYL+f3Cs0TCFBgjU6aJWkJTM7O50PKmuY6czEddhByuz+kyucaRq//kk6htEdhK3WgX/wWwPR2Xta\nQJB2NoTZK/BbFRpyVAxdwR0Y3WEotpIppH3h8zQ884fItlCHh8Znn6Ntx05SP/dZzNkD1wil8U0G\nLWncUHUTSblZ6DYLAY+PlopqQoGB16uKZW6uYE6uQNdM0NksGDAEeyqITDc10ty0cJaTQwpaTaKO\nco5FLS/TxYadYrpnehFCRIKzp0fTYJd6qnCLFhy+JMwH/ZAKCRYLZk1lYV42756qpLHDyy3zZpBs\nCaISINg2tGy+wZIuekrQuwdx5x4OknvUQOvx0U3+CCpKNBz5576G12DiL5mL0eqi+a9/J/mG6wk2\nNODesg1zPYacAAAgAElEQVTfqdNUPfwICUuvJGnVSlSLnI19IpJBSxoXMkqLySidjNpjEuFJ82dQ\ne+DEgEuq9zY3V0RN3NtF14hsH43AVc0Z3LQQFAFMysBjhuJwEE8izdT3WcV6OpegKeGvZUiEaKCG\ndMIJDB7a0TFjxkI7brLIJ5EUXDThIInWj4KISyDJYuXbSxcwPT2VrafP0tjh4VDNWVYVh5vmTPED\nz7ji95s5dWYKHZ544mxtFOYfw9zPPIFdypKb+FtFHpkHBQUH++5fC0LhIYOsvBDEHlMcceWiAzG3\n9xx0PJiEpVdhtLjQ050kLr2S+AXzaXzuefwVlbS+8Rbtu/aQ+tkbiJvZdyo4aXyTQUs67zJKi8mK\nMcGvatIi24cSuMxauIY1kDm5gv1VndMwjZCQCFHTubCni6Z+l43pYlVsWIWtz1LzORSRoqRHHrfQ\nQDVnIkHL1jkFVQXHaecIAfxMUrrHnp1qgMtbBFO1ZCbNCGAIhSWFGbx1soZXjpxmVXEKqlmQOC16\n3r2e9uxbwN59CwgGuwPv9p1XM3vm+8yZ+X7/70kzWJpcTfvR/pMdHCY/zduTyVzQjGYZ3bFUiqKQ\nvPrTkaQMS14eWXffhXvLVppfegWjpYW6JzdhK51BymduQE/9+H1t0tiSPZPSeaXqJjJKBx70G66B\nDX5/VZgWvaZXLLoGRSM8XKuRGgKdY8Va6JuZ11NIhDgi9nCY3QgENsIJIhZsTGFmVNkGqmmilpAI\nB7dspQCLYiOOcGp3O9H9Un5DYU9FeCb2OFMQhx7gc6U5ALx7up76di/plzf1GzD27FvAh7sXRwUs\ngGBQ58Pdi9mzb8GA7212SxvJip/eXY2KAgm6H4ceIORXcB0em3FTiqqi2mxRjxOWXEHOxu9gnzsH\nAM+Bg1Q9+FM8R4+NyTlJ507WtKTzKik3K6pJMBbVpJGUl0nTIEuYxw2xj3+o5YaqqnMsGITT7PsT\nFAH2sZ1GagHIooASZvEOf2M6l/RpVmygGoMgzdSTSmZku53wj76HdgxhRGUZdjV9Bow4dA1SRR5p\nloM0+Dz86Pl63nog9gKFfr+ZvYMEpb37FjBj6p5+mwoDbRoOPYBdD+INahhCQVMEVpOB2mP1664+\ntd/mvRtzP/3NSTgSzYYApsQEnP96C/GXLaDxz88jAgEs+XnD2od0/sigJZ1X+hCXJtdtgy+f3jHE\nPv6hlhsKv/DRQPdMDy4aCYkQqhLdiNEh2tjLlshEvsXMJJ8SFEVhqphDmpIZVb5duCNjzOqpjhm0\nADpw4yAp6rW7KxT2V9E5Q75CoV5Ag+8QH7hOEwploqp9m0ZPnZnSp4bVWzCoc7p8CiXFsYOH3tlX\npiKIM/WfQDNYn9pYsU0tIfs73yLY2CiTMiYQGbSk8yowxKXJA57BJ2491QCXTx64iTBgwMmBW/CG\npZZKRI9ahIFBGy4S6F6gtFnU8xHbCOBHQ6OMhVGztPfsl+rSMxA2UI0QcyJZhCZFxyJs+PDQHiNo\nQbjP7ki4QodFFAKH8NDOwpycqH6zLuF+RVuf7SWzPFGPOzz9j3dLnNZG1WtOQv7++wsH61OD/mtg\n/fnSsEpH18xUXcecmTlAaWm8kX1a0nnVUlFNaJD1w0JBg5bywZen6OrTGcieCmVEkzCqOd1nW89+\nrSpxml28QwA/FmzM5+p+lxXpqb5HOryXjj79V121rd7bY7EqcZGa2llOxSwz1JuHOFvftPsumjVE\n+uVN/T4PDNinJklDIYOWdF6FAkFqB1mlufbAie4FEwexu0Jh52mFQK84GDBg52llRNPd20UrrTQT\nR3h+PqXzfy00IoTgmPiIg3yAQJBACgtYjkPpWyvqLSD8ffrG6omebNbemYzRwcC1li45FAJQx1n8\nom+AGsrNg8kUoCBv4ISF9MXNZC5tRDVHBybVLMhc2kj64oHnL5SkwcjmQem860pn7z1OKxQ0hj1O\nC3r36YT7sE42jGyaO4SbAi9nFS6aOcD72IhnNpdTzlHcuCgn/AOfQS4zmN9nWqb+NFIT1eQI4SbC\nQqZHHnfXtIYWtNLIwowFPz5qKCePKVHPd9089B56cPSj7ibD6r0V/N97ZgDwatWefo+VvriZ1Pkt\nuA47CLZpmOINEqe5R62GNVIJHbEMN8lDGn0yaEnjQu2B49QfOd3/kvTD1LNPZ7QkKOF+q0YRPpAJ\nHbviYDqXADBdXIKXDgqZPqwZ5ruaBhUUBAIFFRdN+IUXsxJOSInrDFoduKNmzeiPqqjki6kYBEln\nUswyI3nzoFnEgFNFSdLHJYOWNG6EgsFB09rHoyDhKYx0dLB0LhPiUMh3lxDY74ehdRcB4WmbdMws\n4hp28S4+PJRyKX58uGmJ9E11NQ+GCOGhPdJEOZB8pe8A7t5G+uZBkkaaDFqSdI66gpY53Ur8+kR6\nTjpvuc6G/20v/rcGjlw9511M9eTRUlENnan5Kip5SnFUeTMWMsnDShwqH39l4FiGcvOwKntOzO0D\nNRuOteE2G8YyUmPDpJEjg5YknaNAZ9CyFtjovUqKYgbLynCTXn+Bq795F997/h8Q8MR8jaIolDHw\nYGBJuhDJoCVJ5yio+iEEuq3/qTbMV1nxb/P1aSocaN5FTTdBoM9T0jkazRpYVcyt0kg6L0Hr1ltv\nZerUqQAsWbKEqVOn8vDDD5OWlobNZuOuu+5i06ZNVFZW4na7uf3227Hb7X3KSNJ4YNiD4AbdakY3\nBHmuALagwGNSKE/UCWgKijnc1xX4sHs6jqHMuwigahoMf7FfSbognZegpSgKiYmJ1NfXM2nSJJ59\n9lnWrl3L/Pnz2bBhA9XV1WzZsoXHH3+ciooKHnvsMdLT0/uUycqSi7lJo2+wdb66mgfzAyqfOdyG\nKdSd2j2/2scBp5n96RYUR3SGX2TeRQVMZjOKoiCEIOj30zPj3e5MglHOhBwvPs6yKNLFZUyC1qZN\nm9i2bVvk8caNG5k6dSoul4vvfe97JCQkkJERXs4hPT2diooKUlLCSwVkZGRQX1+PqqpRZerr62XQ\nkkbdUNb5CvrDP6pTvWpUwAIwhQSza8Ntgh+4o5/TbRZMVgu6zRo1M7oeZ4uatkozj+5qv+PFx10W\nZaTEajYcTpOhNDbGJGitW7eOdevWRR6/8cYbTJs2jbi4OAKBAFlZWdTW1pKbm0t1dTVFRUU0N4dH\nzldVVZGTk4PT6Ywqk509+FQ4knQuhrrOlzlgRrdacVj6Dy6ltX4+OGQA3dHJ7kzGHNd3ImBFAXOc\nNTL2yvD7gcEnDJ7IupZF6a1rWRRgTAKXNP6dl+bBw4cP89Zbb2EYBrfddhsFBQU89NBDbN68mYKC\nAtLS0liyZAn33Xcfbreb9evXY7fb+5SRpNEy1HW+6o+c5ubUq1l6CSgJ/c+KprUaFDngSEf3/uMz\nUhGCPutPRXQ+0V7fgo2JManrx0mFH4llUUaCrFVNDOclaK1fv77Ptoceeijq8Re/+MVBy0jSaBnO\nOl9xgQqEWwAhFIfaszIFAoQ7hHAL4szdTyTlZqFqKkGPDz1u4GUxQoYRvc8LzEgsiyJdPGTKuyTF\nMJx1vjo6Jz4XboFoM1BsCmiAAcIjIkkVPdfx6tp/wBvuuzLZLFE1LiHgjrW/pHzHR8Oee3Gi6fAM\nPptHuFz/y6JIFw8ZtCQphuGs8xW1jpcA0dF3Ytje63j13H/A6yXo9aFZdBRFRYjw6xVFwZoYj6qb\norIVLzRxtrYhlut/WRTp4iGXJpGkGIazztfHWcer9/4FgqAvXBXT46yY7VZMNgvJBdmUrVlORmlx\nn31eKArzj2EyDTyKeijLokgXB1nTkqQY+luqo6ee63x1rdM1J1dErZwcMMIBq/c6XrH2r1utUf1b\nXWnvvbMVLzRms5/ZM9+PmT3YZfbM90csCUMmXExsMmhJUj+Gu1THcNfx6rl/TTdh6uznEiIcsILe\n6CbKrmzFC3HG9a509t7jtEymwJiN05ImBhm0JGkAw12qY7jreHXtP3dBGc5pZoQQGL7YNYqubMWx\nWr5lsJlARtqcme8zY+oeTpdPocNjJ87WTkGenBFjON59911eeOEF7rvvPuLi4s5pXz/+8Y/ZuHEj\nTz/9NLfeeusIneG5k0FLkgYx2ut8hYJBvK62PjWrWHTb2AwyHspMIMMRa/xWrLFbZrNfprWfgxdf\nfJGWlhZ+9KMfkZ2djaqqrF69mo0bN3Ldddexe/duZs6cyZ49eygpKcHn85GUlMTy5cv59a9/TUpK\nCunp6cycOZP33nuPAwcOsH37dm666Sbuu+8+srKyaGlp4d5772XNmjWsXbuWrVu38pOf/ASbzTb4\nCY4AmYghSePAcLIVR1vXTCC9x6l19a1dyEkhE92iRYuYN28eN9xwA3feeSdHjhzBMAzy8vJYu3Yt\nAKtWrWLlypU4HA7Wr1/Pjh07sFqtJCYmEhcXx3vvvcesWbPIz8+ntDS8bth7773HvHnzuOOOOwA4\ne/YsKSkpfP7zn2fy5MmUl5eP2XuUNS1JGgdaKqqZNH/GgAOau7IVR9NwZgKZqH1rI7E0yXjXNWwC\nwkMnejYVms1mVFWN/L9hGPzlL3/hyiuvZNasWXzwwQdD2qfFEu6DVVWVUGjsliGQQUuSxoHhZiuO\nluHMBDJWfWvS8BQVFfH3v/+dnTt3MmvWLFR18Aa1GTNm8Mc//pE9e/agKAonT57E7/ezc+dOILyE\n1P33309tbS1Wq/W8zv2qiJ7h8wJQWVnJ8uXLyT41DVPw4pgdW7pwxOpL6i9bcbSOP1Dg7FK99+g5\nn89A8xGeDyNR06q6zB1z+2uh585531KYrGlJ0jgyULbiWGTzjae+NUmKRQYtSRpnYmUrjnQ2X3/G\nS9+aJPVHBi1JGueGuq7XSBgvfWuj6UJKuLgYyZR3SRrHhprNp5pG7v6z9sBxqvce7TP3YihojEhf\nliSdC1nTkqRx7Hxl8w13JhBpZHl8QfYeq6e13U+C3czsKU5sFvlzDTJoSdK4Npx1vUbaaM8EMtpG\nohnwnW2lwypfzPZzPuZrO86weWc5/kB3TfeFt46z4tI8Vi7MP+f9T3TDDlrt7e3Y7XIxNkkaCzKb\n7+Ly2o4z/GPrqT7b/QEjsn2ogev555/npZdeoqioCIDs7Gy+9KUvxSxbWVnJb37zG+6///4+z/3X\nf/0XlZWVmEwmmpqa+NrXvsa0adOG+pYAeO6557jmmmv42c9+xg9/+MNhvba3AYPWDTfcwH333UdZ\nWVlk27e+9S0ee+yxczqoJElDI7P5Lh4eX5DNOweeDmnzznKWzMnBOsSmwuuvv57Vq1eH9+/x8O//\n/u9kZ2dTXV3ND3/4Q+644w6Ki4v53Oc+B4R/37/zne+QkZHBV7/6Vb7whS/g8/kigaa9vR2v10tl\nZSX33HMPS5cuZdmyZTz99NM4nU50XWf9+vVcd9113HzzzezcuZM777yTDz/8kCVLlrB9+3ZefPFF\ncnNzef3117HZbJhMJm6//fYhX6cBEzEcDgdPPfUUzz0nB8ZJ0vnQlc03kImezSeF7T1WH9UkGIs/\nYLD3WP2Q9/niiy9y//33c//99/Pmm2+SnZ1NXFwcJ06coK6ujvb2du644w40LXxTdMMNN/CnP/2J\nw4cPM2XKFE6cOMG8efMA+Mc//sEvfvELfvOb3wCQkZHB1772NRwOB6mpqSQmJvL2228DEB8fz7/8\ny79w5ZVXsnfv3sj5ZGdnc/311/Pb3/4WXdcJhUIcPHhwWNdpwKBlt9v5+c9/TkNDA9///vcJBAZe\nXVSSpJEns/kuDq3tQ1uCZajlIFzTuvfee7n33nsxm804HA5uu+02kpOTMQwDTdOiZmdfvHgxW7Zs\n4c9//jOf/exnmTFjBjt27ADgk5/8JBs2bKCyMtzPGR8fD8BTTz3FJz7xCW6++eZIjLBaw32siqL0\nOy/hzTffzJ133snGjRuH/H5giH1aX//613nnnXf4+te/jsvlGtYBJEk6dzKbb/gm2sS4CfahTTs3\n1HIQrmnt378fgKqqKvx+Px0dHTidTl599dU+5VVVZfHixWzdupXJkyczefJkDh06xH/8x39gNptp\naWnh5ptvjnrNnDlzePLJJykqKiIjI6PfCXcBSkpKeOKJJ7jtttt48MEHSU1NJSUlZVjNgwPOPfjG\nG2+wbNmyyOPy8nIeeOCBcd2nJecelKSJ4XzNPTicoDXs7MG7Y2cPDmXuQY8vyPcf3zZgE6FZ1/jB\nVxcNuU/r43jrrbeoq6vjpptuGrVjnIsB33nPgAWQl5f3sQLWrl27+OUvf8maNWtYvXo1tbW1PPzw\nw6SlpWGz2bjrrrvYtGkTlZWVuN1ubr/9dux2+6BlurJiJEmSJjqbxcSKS/NiZg92WXFp3qgGrNde\ne43Nmzfzox/9aNSOca7GZJyW0+nk+uuvjzx+9tlnWbt2LfPnz2fDhg1UV1ezZcsWHn/8cSoqKnjs\nscdIT08ftEys9ExJksaf8TSj+3huNuxKZ+89Tsusa2MyTmvlypWsXLlyVI9xrkYlaG3atIlt27ZF\nHt99991Rzzc0NJCRkQFAeno6FRUVpKSkAOGMlPr6elRVHbSMJEnShWblwnyWzMnpMyPGaNawJpJR\nuQrr1q1j3bp1Udt6pjVmZWVRW1tLbm4u1dXVFBUV0dzcDIQ7C3NycnA6nYOWkSRJGk1XLjoQc/tw\n+7qGTQ2iJtehxbWhWuJBTUZOYBQ2JlfhmWee4e2330ZVVQKBADfddBMPPfQQmzdvpqCggLS0NJYs\nWcJ9992H2+1m/fr12O32QctIkiRdaN44uZW3Tm3Fb3QPMfrb4ddYWng5y4ouP49nNj7IlYslSRp1\n46lPqz8jkVV4LtmDEA5Y/zz+dr/PX1N81UUfuGR9U5IkidgJGv0Fsv6aDavO4fjegJe3Tm0dsMxb\np7Zyed4lWE2DT6RcV1fHgw8+SHZ2Nm63m5ycHFJTU/nsZz875HPasWMHH374YcyWrVWrVvHd7343\nkmW+ceNGsrOzufPOO/uUffTRR+no6KCoqIhrrrmGxMTEqOd///vfc/ToUdrb27nxxhtZtGhRv+ck\ng5YkSdI4sK/uSFSTYCx+I8D+2sPMz5k96P4OHDhARkYG3/zmN1FVlZMnT/Lkk0+iKAp79uwhIyOD\nyspKZs6cyZYtW/jxj3/M3XffzeLFi2lpaSE3N5e8vDwATp8+zbPPPktSUhIul4vvfve7lJWV8eqr\nr7Js2TKamprw+cKTO9fV1fHII4+Qk5OD2+3mtttu480332TZsmWROQj/+te/UllZSV1dHd/+9rcp\nLi7mlltu4fjx4zz33HMyaEmSJH0cY5ke7/a1Dalcq699SOWWLl1KS0sLP/jBDxBCRM3MPnPmTG68\n8Ua+8IUv8OMf/5j6+npOnDiB3+/npptuIj4+nttuuy0yU8Uf//hHgsEggUCAqqoq3G43uq4zadIk\nzpw5w2uvvcaaNWvYs2cPVquVjIwM7HY7L7/8Mhs3bqSkpIQbb7yRX/ziF0B4APNTTz0VmWHpsssu\no6mpiSeffJK77rprwPclg5YkSdI44LDED6lcgmVoS0MdPXqUZcuWsWbNGgzD4PLLL2fFihUAmM1m\nFEXBYgk3M6qqimF0jwsTQtA73eHTn/40s2fPpqamBofDAcBNN93E008/jdvt5tprr2XPnj08//zz\nlJWVsWLFCl544YWY59a1b8MwCAQCVFdX84c//IENGzb0aTrsTQYtSZKkcWBm+lT+dvi1AZsIzZpO\nWcbQ1rISQvCjH/0Ip9OJ3+/n3/7t3zh58uSAr1FVlT//+c9UVFRETQjx+c9/nv/8z/8kMzMTwzAi\nk9ymp6fT2NgYCYYApaWl/O53v+Po0aMUFxfzyiuv9DnO8uXLeeCBB2hsbOQb3/gG3/nOd1i4cCGP\nPfYYRUVFkaVSYpHZg5IkjbqJkD04ElZlz4m5faJkD65bt45NmzaN2v5HgqxpSZIkjRNdAan3OC2z\npstxWp1kTUuSpPPiQqx9nWtNq4s36GN/7WFafe0kWOyUZUwbUpr7xUDWtCRJksYZq8kypLT2i5EM\nWpIkSeOM4fHQ8tE+gq2tmBISSJo1E63HCsMXMxm0JEk6L/prSuvPhdicGEvt5tepe/1NDJ8/sq3q\nLy+SvvxqMlYsP49nNj7IoCVJkjRO1G5+nep/vNpnu+HzR7YPNXCN9jROq1evZsGCBRiGQTAY5Ic/\n/CGPPvooixYtory8nJdeeimyUO/ll1/O1VdfPeTjDkQGLUmSJoQLvWZmeDzUvf7mgGXqXn+TtCsW\no1mtg+5vtKdxSk5O5t577wXgy1/+MsFgMOr4119/PatXr/6YV6N/MmhJkiSNAy0f7YtqEozF8Plp\n2buP1IWXDrq/0Z7GqaWlhYcffpimpiaysrJQFCXq+C+++CL79+8HwisiL1iwYLiXJCYZtCRJksaB\nYGvr0Mq5h1ZutKdxSkpK4p577gHgiSee4O23owdFy5qWJEnSMEy05kRTQsLQyjmGVm60p3FqaWnh\nJz/5CQDV1dWsWbOGAwdiL9kykuTgYkmSJEYmaJ3L4GLD4+HgD+8fsIlQs5iZ8b3/M6Q+rY9DTuMk\nSZI0QQy3ZjbSNJuN9OVXx8we7JK+/OpRC1gThQxakiRJ40RXOnvvcVqaxTwm47TGey0LZNCSJEka\nVzJWLCftisW07N1H0N2KyZFA0uyZF30Nq4sMWpIkSeOMZrUOKa39YiSDliRJ0jjj9wU5dayBjnY/\ncXYzhVPSMFvkzzWMUdDatWsXv/zlL1mzZg2rV69mx44dPPLII8ycOROAe++9l02bNlFZWYnb7eb2\n22/Hbrfz8MMPk5aWhs1m46677upTpmuKEEmSpAvF7h3l7NlZQTDQPW5q61snmHNpLnMX5p3HMxsf\nxiRoOZ3OqJx/AIvFgtVqJSEhAa/Xy5YtW3j88cepqKjgscceIz09nbVr1zJ//nw2bNhAdXV1nzL3\n33//WJy+JEnSmNi9o5wPtp7usz0YMCLbhxq4RnvuwVWrVvHd736XZcuWAbBx40ays7O58847+5R9\n9NFH6ejooKioiGuuuYbExMSo53tXbAYyKkFr06ZNbNu2LfL47rvvjnq+rKyMhx9+GKfTyU9/+lMO\nHjxISkoKABkZGdTX16OqKhkZGQCkp6dTUVHRp4wkSdKFwu8LsmdnxYBl9uysoHRO9pCaCkd77sGy\nsjJeffVVli1bRlNTEz6fDwgHy0ceeYScnBzcbje33XYbb775JsuWLePDDz9kyZIl/PWvf6WyspK6\nujq+/e1vx6zY9GdUgta6detYt25d1LaDBw9G/l1TU4Ou6wDY7XaCwSDNzc0AVFVVkZOTg9PppLa2\nltzcXKqrqykqKupTRpIk6UJx6lhDVJNgLMGAwaljDUwtyxx0f6M996Cu60yaNIkzZ87w2muvsWbN\nGvbs2YPVaiUjIwO73c7LL7/Mxo0bKSkp4cYbb+QXv/gFAG+99RZPPfUULpcLgMTERHbu3Dmk6zQm\nzYPPPPMMb7/9NqqqEggEuOyyy3jwwQfJz8/H5/Nx6aWXcuTIEe677z7cbjfr16/Hbrfz0EMPsXnz\nZgoKCkhLS2PJkiVRZSRJki4UHe0DT5Y73HKjPfcgwE033cTTTz+N2+3m2muvZc+ePTz//POUlZWx\nYsUKXnjhhZjn1rVvwzAIBAJ9mgsHMiZB64tf/CJf/OIXo7b96le/6lOmt4ceemjQMpIkSReCOPvQ\npp0barnRnnsQwl03jY2NkWAIUFpayu9+9zuOHj1KcXExr7zySp/jLF++nAceeIDGxka+8Y1v9KnY\n3Hjjjf2eo5x7UJIkaZQNZe5Bvy/IHx7fMWAToUnXuPmrC0ct/X0izD2onu8TkCRJksBsMTHn0twB\ny8y5NPeiH691cb97SZKkcaQrnb33OC2Tro3JOK3xXssCGbQkSZLGlbkL8yidky1nxOiHvAqSJEnj\njNliGlJa+8VI9mlJkiRJE4YMWpIkSdKEIYOWJEmSNGHIoCVJkiRNGDJoSZIkSROGDFqSJEnShCGD\nliRJkjRhyKAlSZIkTRgyaEmSJEkThgxakiRJ0oQhg5YkSZI0YcigJUmSJE0YMmhJkiRJE4YMWpIk\nSdKEIYOWJEmSNGHIoCVJkiRNGDJoSZIkSROGDFqSJEnShGEai4P8/ve/5+jRo7S3t3PjjTdSVFTE\nww8/TFpaGjabjbvuuotNmzZRWVmJ2+3m9ttvx263D1qmqKhoLE5fkiRJGifGJGgVFxdzyy23cPz4\ncZ577jnef/991q5dy/z589mwYQPV1dVs2bKFxx9/nIqKCh577DHS09MHLXP//ff3OZZhGAAETf6x\neGuSJEmDCgaDmExj8nN7wRuVq7hp0ya2bdsWeXz33XfT1NTEk08+yV133cWvf/1rMjIyAEhPT6ei\nooKUlBQAMjIyqK+vR1XVQcvE0rW9LvfkaLw1SZKkYaupqWHSpEnn+zQuCKMStNatW8e6desijw8f\nPswjjzzChg0bSExMJCsri9raWnJzc6murqaoqIjm5mYAqqqqyMnJwel0DlomlrKyMv7whz/gdDrR\nNG003p4kSdKwZGZmnu9TuGAoQggx2ge5/vrrWbhwISaTiaKiIq6++moeeughUlJScDgc3HHHHTzz\nzDOcOXMGt9vN+vXrsdvtg5bJz88f7VOXJEmSxpExCVqSJEmSNBJkyrskSZI0YcigJUmSJE0YF3UO\n5l//+lf2798PwDvvvMM111xDMBiksbGRDRs20NDQwBNPPEFCQgKFhYX8v/btHqR1KAzj+KOLUysO\ngg7dxMFBLGRwUxT82gTRQVBQ6WAGcZAWHdRRUAxiERehKriIU4s6VLqUDH7QwUkogoUoaGtAQSua\n9y6mestdbu/Qm/j81pzh/GngJWnO8PAwVlZW/nqN/dXj/9AYDAaxsbGBYDAIRVFwdXXlusbe3l48\nP0QOnd4AAALaSURBVD/j/v4ek5OTqKiocHzjnzo7OjqQz+dhmiZUVcXHx4fjO4sbj4+PcXFxgUAg\ngLOzM1fcr/SPhCQajcrm5qbMzs6KiIiu6xIOh2VmZkYMwxARkfHxcUmn0yWt+R9Eo1GJxWJyc3Mj\na2trcnp6KiLiykZd10VEJJFISCQScVWjyFen/RvGYjHZ29tzVafdmM1mZWFhQUZHR0XEXfcrlebH\nvx7M5/OIx+NQFKVwLsw+B5bNZgufqlZXV+Ph4aGkNeVmN/b19cHn8/12zY2Nra2tyGQyODw8RH9/\nv2sagd87FUVBOBzGzs4OOjs7XdNpN/b09EDTNExPTxeuuaWRSvfjh1Y8Hkd7e3vh7BjwdQ6srq4O\nd3d3AADTNOHz+UpaU25245+4sVHXdezu7mJxcREej8c1jcBX58vLC1KpFFRVRSgUwtbWlms67cbL\ny0tYloVIJIJMJoP9/X3XNFLpfvR/WgCQSqUwODiI+vp61NbWYmlpCblcDqFQCLlcDqurq/B6vejq\n6ip5TbnZjQCwvr6OZDKJ6+trmKaJsbExVzW+vb1hbm4O3d3d0DQNLS0trmkEvjqrqqpwcHCAo6Mj\nGIaBkZER1NTUuKLTbmxoaEBzczMA4Pz8HAMDA/D7/a5opNLxnBYRETnGj389SEREzsGhRUREjsGh\nRUREjsGhRUREjsGhRUREjsGhRUREjsGhRfRNMpksnOOxLAsTExNIp9Nl3hUR2XhOi6jI8vIympqa\nYBgGPB4PhoaGyr0lIvrEoUVU5P39HYFAAF6vF5qmlXs7RPQNXw8SFXl6ekJlZSUeHx/x+vpa7u0Q\n0Td80iIqMjU1BVVVcXt7i0Qigfn5+XJviYg+8UmL6Jvt7W34/X40Njaira0NlmXh5OSk3Nsiok98\n0iIiIsfgkxYRETkGhxYRETkGhxYRETkGhxYRETkGhxYRETkGhxYRETkGhxYRETkGhxYRETnGL995\ng3V4is87AAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Plot the block model. \n", "GeMpy.plot_section(geo_data, 13, block = geo_res_num.as_matrix(), direction='y', plot_data = True)" ] }, { "cell_type": "code", "execution_count": 57, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "125000" ] }, "execution_count": 57, "metadata": {}, "output_type": "execute_result" } ], "source": [ "50*50*50" ] }, { "cell_type": "code", "execution_count": 78, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 0., 1., 2., 3., 4.])" ] }, "execution_count": 78, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.unique(sol[-1,0,:])" ] }, { "cell_type": "code", "execution_count": 79, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "formation number\n", "1 [EarlyGranite]\n", "2 [SimpleMafic2]\n", "3 [SimpleBIF]\n", "4 [SimpleMafic1]\n", "Name: formation, dtype: object" ] }, "execution_count": 79, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Formation number and formation\n", "data_interp.interfaces.groupby('formation number').formation.unique()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "data_interp.interpolator.tg.u_grade_T.get_value()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "np.unique(sol)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#np.save('SandstoneSol', sol)\n", "np.count_nonzero(np.load('SandstoneSol.npy') == sol)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "sol.shape" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", ">> Welcome to the Python client library for Steno3D!\n", "\n", "Credentials file found: /home/miguel/.steno3d_client/credentials\n", "Accessing API developer key for @leguark\n", "Welcome to Steno3D! You are logged in as @leguark\n", "Verifying your quota for public projects...\n", "This PUBLIC project will be viewable by everyone.\n", "Total progress: 100% - Uploading: project Sandstone\n", "Complete!\n", "https://steno3d.com/resource/volume/BWCZxkJUmmdHurZg3mGs\n" ] }, { "data": { "text/html": [ "\n", " \n", " " ], "text/plain": [ "" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "GeMpy.PlotData(geo_data).plot3D_steno(sol[-1,0,:], 'Sandstone', description='The sandstone model')" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(array([ 696000. , 697307.692308, 698615.384615, 699923.076923, 701230.769231, 702538.461538, 703846.153846,\n", " 705153.846154, 706461.538462, 707769.230769, 709076.923077, 710384.615385, 711692.307692, 713000. ,\n", " 714307.692308, 715615.384615, 716923.076923, 718230.769231, 719538.461538, 720846.153846, 722153.846154,\n", " 723461.538462, 724769.230769, 726076.923077, 727384.615385, 728692.307692, 730000. , 731307.692308,\n", " 732615.384615, 733923.076923, 735230.769231, 736538.461538, 737846.153846, 739153.846154, 740461.538462,\n", " 741769.230769, 743076.923077, 744384.615385, 745692.307692, 747000. ]),\n", " 1307.6923076923076)" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.linspace(geo_data.extent[0], geo_data.extent[1], geo_data.resolution[0], retstep=True)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(39,)" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.diff(np.linspace(geo_data.extent[0], geo_data.extent[1], geo_data.resolution[0], retstep=False)).shape" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1271.0" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "(geo_data.extent[1]- geo_data.extent[0])/ geo_data.resolution[0]-4" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "1307.6923076923076" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "(geo_data.extent[1]- geo_data.extent[0])/39" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "scrolled": false }, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "2fa1aa6439df4ed4afcf0b129ac2660c" } }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# So far this is a simple 3D visualization. I have to adapt it into GeMpy \n", "\n", "lith0 = sol == 0\n", "lith1 = sol == 1\n", "lith2 = sol == 2\n", "lith3 = sol == 3\n", "lith4 = sol == 4\n", "np.unique(sol)\n", "\n", "import ipyvolume.pylab as p3\n", "\n", "p3.figure(width=800)\n", "\n", "p3.scatter(geo_data.grid.grid[:,0][lith0],\n", " geo_data.grid.grid[:,1][lith0],\n", " geo_data.grid.grid[:,2][lith0], marker='box', color = 'blue', size = 0.1 )\n", "\n", "p3.scatter(geo_data.grid.grid[:,0][lith1],\n", " geo_data.grid.grid[:,1][lith1],\n", " geo_data.grid.grid[:,2][lith1], marker='box', color = 'yellow', size = 1 )\n", "\n", "p3.scatter(geo_data.grid.grid[:,0][lith2],\n", " geo_data.grid.grid[:,1][lith2],\n", " geo_data.grid.grid[:,2][lith2], marker='box', color = 'green', size = 1 )\n", "\n", "p3.scatter(geo_data.grid.grid[:,0][lith3],\n", " geo_data.grid.grid[:,1][lith3],\n", " geo_data.grid.grid[:,2][lith3], marker='box', color = 'pink', size = 1 )\n", "\n", "p3.scatter(geo_data.grid.grid[:,0][lith4],\n", " geo_data.grid.grid[:,1][lith4],\n", " geo_data.grid.grid[:,2][lith4], marker='box', color = 'red', size = 1 )\n", "\n", "p3.xlim(np.min(geo_data.grid.grid[:,0]),np.min(geo_data.grid.grid[:,0])+2175.0*40)\n", "p3.ylim(np.min(geo_data.grid.grid[:,1]),np.max(geo_data.grid.grid[:,1]))\n", "p3.zlim(np.min(geo_data.grid.grid[:,2]),np.min(geo_data.grid.grid[:,2])+2175.0*40)#np.max(geo_data.grid.grid[:,2]))\n", "\n", "p3.show()" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "scrolled": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Function profiling\n", "==================\n", " Message: :3\n", " Time in 5 calls to Function.__call__: 1.357155e+01s\n", " Time in Function.fn.__call__: 1.357096e+01s (99.996%)\n", " Time in thunks: 1.357014e+01s (99.990%)\n", " Total compile time: 2.592983e+01s\n", " Number of Apply nodes: 95\n", " Theano Optimizer time: 1.642699e+01s\n", " Theano validate time: 3.617525e-02s\n", " Theano Linker time (includes C, CUDA code generation/compiling): 9.462233e+00s\n", " Import time 1.913705e-01s\n", " Node make_thunk time 9.450990e+00s\n", " Node forall_inplace,cpu,scan_fn}(Elemwise{Maximum}[(0, 0)].0, Subtensor{int64:int64:int8}.0, Subtensor{int64:int64:int8}.0, Subtensor{int64:int64:int8}.0, Subtensor{int64:int64:int8}.0, Subtensor{int64:int64:int8}.0, Subtensor{int64:int64:int8}.0, IncSubtensor{InplaceSet;:int64:}.0, grade of the universal drift, , , Value of the formation, Position of the dips, Rest of the points of the layers, Reference points for every layer, Angle of every dip, Azimuth, Polarity, InplaceDimShuffle{x,x}.0, InplaceDimShuffle{x,x}.0, Elemwise{Composite{((sqr(sqr(i0)) * sqr(i0)) * i0)}}.0, Elemwise{Composite{(sqr(sqr(i0)) * i0)}}.0, Elemwise{Composite{(sqr(i0) * i0)}}.0, Elemwise{mul,no_inplace}.0, Elemwise{neg,no_inplace}.0, Elemwise{mul,no_inplace}.0, Elemwise{true_div,no_inplace}.0, Elemwise{Mul}[(0, 1)].0, Elemwise{mul,no_inplace}.0, Elemwise{Composite{(i0 * Composite{sqr(sqr(i0))}(i1))}}.0, Elemwise{Composite{(((i0 * i1) / sqr(i2)) + i3)}}.0, Reshape{2}.0) time 9.275085e+00s\n", " Node Elemwise{Composite{Switch(LT((i0 - Composite{Switch(LT(i0, i1), i0, i1)}(i1, i2)), i3), (i4 - i0), Switch(GE((i0 - Composite{Switch(LT(i0, i1), i0, i1)}(i1, i2)), (i2 - Composite{Switch(LT(i0, i1), i0, i1)}(i1, i2))), (i5 + i0), Switch(LE((i2 - Composite{Switch(LT(i0, i1), i0, i1)}(i1, i2)), i3), (i5 + i0), i0)))}}(Elemwise{Composite{minimum(minimum(minimum(minimum(minimum(i0, i1), i2), i3), i4), i5)}}.0, TensorConstant{1}, Elemwise{add,no_inplace}.0, TensorConstant{0}, TensorConstant{-2}, TensorConstant{2}) time 3.851414e-03s\n", " Node Elemwise{Composite{Switch(i0, Switch(LT((i1 + i2), i3), i3, (i1 + i2)), Switch(LT(i1, i2), i1, i2))}}(Elemwise{lt,no_inplace}.0, Elemwise{Composite{minimum(minimum(minimum(minimum(minimum(i0, i1), i2), i3), i4), i5)}}.0, Elemwise{Composite{Switch(LT((i0 + i1), i2), i2, (i0 + i1))}}.0, TensorConstant{0}) time 3.796577e-03s\n", " Node Elemwise{Composite{minimum(minimum(minimum(minimum(minimum(i0, i1), i2), i3), i4), i5)}}(Elemwise{Composite{Switch(LT((i0 + i1), i2), i2, (i0 + i1))}}.0, Elemwise{sub,no_inplace}.0, Elemwise{Composite{Switch(LT((i0 + i1), i2), i2, (i0 + i1))}}.0, Elemwise{sub,no_inplace}.0, Elemwise{Composite{Switch(LT((i0 + i1), i2), i2, (i0 + i1))}}.0, Elemwise{sub,no_inplace}.0) time 3.589630e-03s\n", " Node Elemwise{Composite{(((i0 - maximum(i1, i2)) - i3) + maximum(i4, i5))}}[(0, 0)](Elemwise{Composite{Switch(LT(i0, i1), (i0 + i2), i0)}}.0, Elemwise{Composite{minimum(((i0 + i1) - i2), i3)}}.0, TensorConstant{1}, TensorConstant{1}, Elemwise{Composite{((maximum(i0, i1) - Switch(LT(i2, i3), (i2 + i4), i2)) + i1)}}[(0, 2)].0, TensorConstant{2}) time 3.567696e-03s\n", "\n", "Time in all call to theano.grad() 0.000000e+00s\n", "Time since theano import 74.908s\n", "Class\n", "---\n", "<% time>