{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Mineral Resource estimation with Python (Quick Demo)\n",
    "\n",
    "Welcome to the **world of Python**. \n",
    "\n",
    "IN this demo we: \n",
    " - import drillhole tables\n",
    " - create drillholes\n",
    " - tag drillholes with domain\n",
    " - composite\n",
    " - create a block model\n",
    " - do statistical analysis (no variography)\n",
    " - interpolate grade and validate interpolations\n",
    " - Report resources\n",
    " \n",
    " <img src ='fig1.JPG'>\n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import drillhole tables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Import libraries \n",
    "import pandas as pd \n",
    "import pygslib\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt    \n",
    "%matplotlib inline "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Import drillhole data located in folder data/ into pandas dataframes\n",
    "collar = pd.read_csv('data/collar.csv')\n",
    "survey = pd.read_csv('data/survey.csv')\n",
    "assay  = pd.read_csv('data/assay.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create a drillholes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create drillhole object\n",
    "mydholes = pygslib.drillhole.Drillhole(collar, survey)   \n",
    "\n",
    "# add intervall table\n",
    "mydholes.addtable(table = assay, table_name = 'assay')  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# validate\n",
    "mydholes.validate()\n",
    "mydholes.validate_table('assay')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Desurvey\n",
    "mydholes.desurvey('assay', endpoints=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Export drillhole as vtk file\n",
    "mydholes.intervals2vtk(table_name= 'assay', filename= 'assay')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Tag drillholes with domain (wireframe solid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import wireframe\n",
    "dom1_poly= pygslib.vtktools.loadVTP('data/dom1.vtp')\n",
    "\n",
    "# creating array to tag samples in domain1\n",
    "inside1=pygslib.vtktools.pointinsolid(dom1_poly, \n",
    "                       x=mydholes.table['assay']['xm'].values, # .values this extracts numpy array from pandas\n",
    "                       y=mydholes.table['assay']['ym'].values, \n",
    "                       z=mydholes.table['assay']['zm'].values)\n",
    "\n",
    "# creating a new domain field \n",
    "mydholes.table['assay']['Domain1']=inside1.astype(int)\n",
    "\n",
    "# Export drillhole as vtk file\n",
    "mydholes.intervals2vtk(table_name= 'assay', filename= 'assay')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Composite"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# extract samples from domain 1 and addit as a new table\n",
    "mydholes.addtable(mydholes.table['assay'].loc[mydholes.table['assay']['Domain1']==1, :], 'assay_d1', overwrite = True)\n",
    "\n",
    "# composite only using samples in Dom1\n",
    "mydholes.downh_composite('assay_d1','Au', 'cmp', cint = 1, overwrite = True)\n",
    "\n",
    "# need to desurvey cmp\n",
    "mydholes.desurvey('cmp')\n",
    "\n",
    "# Export drillhole as vtk file\n",
    "mydholes.intervals2vtk(table_name= 'cmp', filename= 'cmp')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Create a block model with definition: \n",
    "```\n",
    "xorg = 0\n",
    "yorg = 0\n",
    "zorg = 0\n",
    "dx = 10\n",
    "dy = 10\n",
    "dz = 10\n",
    "nx = 18\n",
    "ny = 20\n",
    "nz = 15\n",
    "\n",
    "```\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#create an empty model\n",
    "mymodel = pygslib.blockmodel.Blockmodel(xorg = 0,\n",
    "                                        yorg = 0,\n",
    "                                        zorg = 0,\n",
    "                                        dx = 5,\n",
    "                                        dy = 5,\n",
    "                                        dz = 5,\n",
    "                                        nx = 36,\n",
    "                                        ny = 40,\n",
    "                                        nz = 30)\n",
    "\n",
    "#generate blocks and calculate percent in domain 1\n",
    "_ = mymodel.fillwireframe(dom1_poly)\n",
    "\n",
    "# remove blocks out of domain 1\n",
    "mymodel.set_blocks (mymodel.bmtable[mymodel.bmtable['__in']>0])\n",
    "\n",
    "# and save the model as vtkunestructured grid\n",
    "mymodel.blocks2vtkUnstructuredGrid('bmodel') # this will have extension vtu"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Statistical analysis (no variography)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Declustering"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#declustering parameters \n",
    "parameters_declus = { \n",
    "        'x'      :  mydholes.table['cmp'].loc[:, 'xm'], \n",
    "        'y'      :  mydholes.table['cmp'].loc[:, 'ym'],  \n",
    "        'z'      :  mydholes.table['cmp'].loc[:, 'zm'], \n",
    "        'vr'     :  mydholes.table['cmp'].loc[:, 'Au'],   \n",
    "        'anisy'  :  1.,       \n",
    "        'anisz'  :  .05,              \n",
    "        'minmax' :  0,                 \n",
    "        'ncell'  :  100,                  \n",
    "        'cmin'   :  10., \n",
    "        'cmax'   :  100.,                 \n",
    "        'noff'   :  8,                    \n",
    "        'maxcel' :  -1}               \n",
    "\n",
    "# declustering \n",
    "wtopt,vrop,wtmin,wtmax,error, \\\n",
    "xinc,yinc,zinc,rxcs,rycs,rzcs,rvrcr = pygslib.gslib.declus(parameters_declus)\n",
    "\n",
    "#Plotting declustering optimization results\n",
    "plt.plot (rxcs, rvrcr, '-o')\n",
    "plt.xlabel('X cell size')\n",
    "plt.ylabel('declustered mean')\n",
    "plt.show()\n",
    "plt.plot (rycs, rvrcr, '-o')\n",
    "plt.xlabel('Y cell size')\n",
    "plt.ylabel('declustered mean')\n",
    "plt.show()\n",
    "plt.plot (rzcs, rvrcr, '-o')\n",
    "plt.xlabel('Z cell size')\n",
    "plt.ylabel('declustered mean')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we fix the cell size to 60 x 60 x 3 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#declustering parameters \n",
    "parameters_declus['ncell'] = 1\n",
    "parameters_declus['cmin'] = 60.\n",
    "parameters_declus['cmax'] = 60.\n",
    "   \n",
    "# declustering \n",
    "wtopt,vrop,wtmin,wtmax,error, \\\n",
    "xinc,yinc,zinc,rxcs,rycs,rzcs,rvrcr = pygslib.gslib.declus(parameters_declus)\n",
    "\n",
    "# Adding declustering weight to a drillhole interval table\n",
    "mydholes.table['cmp']['declustwt'] = 1\n",
    "mydholes.table['cmp'].loc[:, 'declustwt'] = wtopt\n",
    "\n",
    "# calculating declustered mean\n",
    "decl_mean = rvrcr[0]\n",
    "\n",
    "print ('declustered mean:', decl_mean)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Variography\n",
    "\n",
    "For now use this variogram model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In Pygslib: \n",
    "```\n",
    "vario_model = {\n",
    "            # Variogram parameters Pygslib\n",
    "            # ----------\n",
    "            'c0'         : 0.1,   \n",
    "            'it'         : [1],    # \n",
    "            'cc'         : [.9],     \n",
    "            'aa'         : [100],   \n",
    "            'aa1'        : [100],  \n",
    "            'aa2'        : [20],   \n",
    "            'ang1'       : [0],   \n",
    "            'ang2'       : [0],  \n",
    "            'ang3'       : [-15]}  \n",
    "\n",
    "```\n",
    "\n",
    "\n",
    "In gslib: \n",
    "\n",
    "```\n",
    "1    0.1                      -nst, nugget effect\n",
    "1    0.9  0.0   0.0  -15.0     -it,cc,ang1,ang2,ang3\n",
    "       100.0  100.0  20.0     -a_hmax, a_hmin, a_vert\n",
    "\n",
    "\n",
    "nst and c0: the number of variogram structures and the nugget\n",
    "it:  the type of structure\n",
    "cc:  the c parameter \"sill\"\n",
    "ang1,ang2,ang3: the angles defining the geometric anisotropy\n",
    "aa: also aa_hmax, the maximum horizontal range\n",
    "aa1: also aa_hmin, the minimum horizontal range\n",
    "aa2: also aa_vert, the vertical range\n",
    "\n",
    "\n",
    "it is \n",
    "\n",
    "1. Spherical (use actual range)\n",
    "2. Exponential (use practical range)\n",
    "3. Gaussian (use practical range)\n",
    "4. Power law variogram\n",
    "5. Cosine hole effect model\n",
    "\n",
    "\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interpolation and validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# creating BHID of type integer, this is to be able to use drillhole id in Fortran!\n",
    "mydholes.txt2intID('cmp')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# creating parameter dictionary for estimation in one block\n",
    "kt3d_Parameters = {\n",
    "            # Input Data (Only using intervals in the mineralized domain)\n",
    "            # ----------\n",
    "            'x' : mydholes.table[\"cmp\"]['xm'].values, \n",
    "            'y' : mydholes.table[\"cmp\"]['ym'].values,\n",
    "            'z' : mydholes.table[\"cmp\"]['zm'].values,\n",
    "            'vr' : mydholes.table[\"cmp\"]['Au'].values,\n",
    "            'bhidint' : mydholes.table[\"cmp\"]['BHIDint'].values, # an interger BHID\n",
    "            # Output (Target) \n",
    "            # ----------\n",
    "            'nx' : 100,  # these parameters are only used to define supperblock search\n",
    "            'ny' : 100,  \n",
    "            'nz' : 100, \n",
    "            'xmn' : 0,  \n",
    "            'ymn' : 0,  \n",
    "            'zmn' : 0,  \n",
    "            'xsiz' : 5,  \n",
    "            'ysiz' : 5,   \n",
    "            'zsiz' : 5, \n",
    "            'nxdis' : 5,  \n",
    "            'nydis' : 5,  \n",
    "            'nzdis' : 3,  \n",
    "            'outx' : mymodel.bmtable['XC'].values,  \n",
    "            'outy' : mymodel.bmtable['YC'].values,\n",
    "            'outz' : mymodel.bmtable['ZC'].values,\n",
    "            # Search parameters \n",
    "            # ----------\n",
    "            'radius'     : 60,   \n",
    "            'radius1'    : 60,   \n",
    "            'radius2'    : 8,   \n",
    "            'sang1'      : 0,  \n",
    "            'sang2'      : 0,   \n",
    "            'sang3'      : -15,   \n",
    "            'ndmax'      : 20,    \n",
    "            'ndmin'      : 7,  \n",
    "            'noct'       : 0,\n",
    "            'nbhid'      : 5,   \n",
    "            # Kriging parameters and options \n",
    "            # ----------\n",
    "            'ktype'      : 1,   # 1 Ordinary kriging \n",
    "            'idbg'       : 0,   # 0 no debug \n",
    "            # Variogram parameters Pygslib\n",
    "            # ----------\n",
    "            'c0'         : 0.1,   \n",
    "            'it'         : [1],    \n",
    "            'cc'         : [.9],     \n",
    "            'aa'         : [100],   \n",
    "            'aa1'        : [100],  \n",
    "            'aa2'        : [20],   \n",
    "            'ang1'       : [0],   \n",
    "            'ang2'       : [0],  \n",
    "            'ang3'       : [-15]}   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# estimating in all blocks\n",
    "estimate, debug, summary = pygslib.gslib.kt3d(kt3d_Parameters)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# adding the estimate into the model\n",
    "mymodel.bmtable['Au_OK'] = estimate['outest']\n",
    "mymodel.bmtable['Au_ID2'] = estimate['outidpower']\n",
    "mymodel.bmtable['Au_NN'] = estimate['outnn']\n",
    "mymodel.bmtable['Au_Lagrange'] = estimate['outlagrange']\n",
    "mymodel.bmtable['Au_KVar']= estimate['outkvar']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# exporting block model to VTK (unstructured grid) \n",
    "mymodel.blocks2vtkUnstructuredGrid(path='bmodel.vtu')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Validating results\n",
    "\n",
    "Basic validations are:\n",
    "\n",
    " - visual validation\n",
    " - comparison of mean grades\n",
    " - swath plots \n",
    " - global change of support (GCOS)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Visual validations \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Mean comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print (\"Mean in model OK   :\",  mymodel.bmtable['Au_OK'].mean())\n",
    "print (\"Mean in model ID2   :\",  mymodel.bmtable['Au_ID2'].mean())\n",
    "print (\"Mean in model NN   :\",  mymodel.bmtable['Au_NN'].mean())\n",
    "print (\"Mean in data    :\", mydholes.table[\"cmp\"]['Au'].mean())\n",
    "print (\"Declustered mean:\", decl_mean)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Swath plots\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mymodel.bmtable.groupby('XC')[['Au_OK','Au_ID2','Au_NN']].mean().plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mymodel.bmtable.groupby('YC')[['Au_OK','Au_ID2','Au_NN']].mean().plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mymodel.bmtable.groupby('ZC')[['Au_OK','Au_ID2','Au_NN']].mean().plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Global change of support  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fit anamorphosis by changing, zmax, zmin, and extrapolation function\n",
    "PCI, H, raw, zana, gauss, z, P, raw_var, PCI_var, fig1 = pygslib.nonlinear.anamor(\n",
    "                         z = mydholes.table[\"cmp\"]['Au'], \n",
    "                         w = mydholes.table[\"cmp\"]['declustwt'], \n",
    "                         zmin = mydholes.table[\"cmp\"]['Au'].min(), \n",
    "                         zmax = mydholes.table[\"cmp\"]['Au'].max(),\n",
    "                         zpmin = None, zpmax = None,\n",
    "                         ymin=-5, ymax=5,\n",
    "                         ndisc = 5000,\n",
    "                         ltail=1, utail=4, ltpar=1, utpar=1.5, K=40)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calc block covariance\n",
    "parameters_cov = {\n",
    "           # Coordinates of the discretization points (using 500 points located at random in the block)\n",
    "            'xdb'  :  np.random.rand(1000)*5, # array('f')\n",
    "            'ydb'  :  np.random.rand(1000)*5, # array('f')\n",
    "            'zdb'  :  np.random.rand(1000)*5, # array('f')\n",
    "           # Variogram model\n",
    "            'c0'         : 0.1,   \n",
    "            'it'         : [1],    \n",
    "            'cc'         : [.9],     \n",
    "            'aa'         : [100],   \n",
    "            'aa1'        : [100],  \n",
    "            'aa2'        : [20],   \n",
    "            'ang1'       : [0],   \n",
    "            'ang2'       : [0],  \n",
    "            'ang3'       : [-15]} \n",
    "\n",
    "\n",
    "cbb = pygslib.gslib.block_covariance(parameters_cov)\n",
    "\n",
    "# calculate the support correction coefficient r\n",
    "r = pygslib.nonlinear.get_r(Var_Zv = cbb, PCI = PCI)\n",
    "\n",
    "print ('cbb :', cbb)\n",
    "print ('r   :', r)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fit block anamorphosis\n",
    "ZV, PV, fig2 = pygslib.nonlinear.anamor_blk( PCI, H, r = r, gauss = gauss, Z = z,\n",
    "                  ltail=1, utail=1, ltpar=1, utpar=1,\n",
    "                  raw=raw, zana=zana)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate grade tonnage courve\n",
    "\n",
    "cutoff = np.arange(0, 3, 0.01)\n",
    "tt = []\n",
    "gg = []\n",
    "label = []\n",
    "\n",
    "# calculate GTC from gaussian in block support \n",
    "t,ga,gb = pygslib.nonlinear.gtcurve (cutoff = cutoff, z=ZV, p=PV, varred = 1, ivtyp = 0, zmin = 0, zmax = None,\n",
    "             ltail = 1, ltpar = 1, middle = 1, mpar = 1, utail = 1, utpar = 1,maxdis = 1000)\n",
    "tt.append(t)\n",
    "gg.append(ga)\n",
    "label.append('DGM with block support')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# see how it look with GCOS. This is also know as global estimation\n",
    "fig = pygslib.nonlinear.plotgt(cutoff = cutoff, t = tt, g = gg, label = label)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# to compare global resources with the one estimated we calculate the CDF of the blocks\n",
    "\n",
    "# cdf of kriging estimate\n",
    "parameters_probplt = {\n",
    "        'iwt'  : 0,                             #int, 1 use declustering weight\n",
    "        'va'   : mymodel.bmtable['Au_OK'][mymodel.bmtable['Au_OK'].notnull()].values,             # array('d') with bounds (nd)\n",
    "        'wt'   : np.ones(mymodel.bmtable['Au_OK'][mymodel.bmtable['Au_OK'].notnull()].shape[0])} # array('d') with bounds (nd), wight variable (obtained with declust?)\n",
    "\n",
    "\n",
    "binval_ok,cl_ok,xpt025,xlqt,xmed,xuqt,xpt975,xmin,xmax, \\\n",
    "xcvr,xmen,xvar,error = pygslib.gslib.__plot.probplt(**parameters_probplt)\n",
    "\n",
    "# cdf of id2\n",
    "parameters_probplt = {\n",
    "        'iwt'  : 0,                             #int, 1 use declustering weight\n",
    "        'va'   : mymodel.bmtable['Au_ID2'][mymodel.bmtable['Au_OK'].notnull()].values,             # array('d') with bounds (nd)\n",
    "        'wt'   : np.ones(mymodel.bmtable['Au_OK'][mymodel.bmtable['Au_OK'].notnull()].shape[0])} # array('d') with bounds (nd), wight variable (obtained with declust?)\n",
    "\n",
    "binval_id2,cl_id2,xpt025,xlqt,xmed,xuqt,xpt975,xmin,xmax, \\\n",
    "xcvr,xmen,xvar,error = pygslib.gslib.__plot.probplt(**parameters_probplt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate GTC ok \n",
    "t,ga,gb = pygslib.nonlinear.gtcurve (cutoff = cutoff, z=cl_ok, p=binval_ok, varred = 1, ivtyp = 2, zmin = 0, zmax = None,\n",
    "             ltail = 1, ltpar = 1, middle = 1, mpar = 1, utail = 1, utpar = 1,maxdis = 1000)\n",
    "tt.append(t)\n",
    "gg.append(ga)\n",
    "label.append('Ordinary Kriging')\n",
    "\n",
    "# calculate GTC in block support\n",
    "t,ga,gb = pygslib.nonlinear.gtcurve (cutoff = cutoff, z=cl_id2, p=binval_id2, varred = 1, ivtyp = 2, zmin = 0, zmax = None,\n",
    "             ltail = 1, ltpar = 1, middle = 1, mpar = 1, utail = 1, utpar = 1,maxdis = 1000)\n",
    "tt.append(t)\n",
    "gg.append(ga)\n",
    "label.append('Inverse of the Distance 2)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = pygslib.nonlinear.plotgt(cutoff = cutoff, t = tt, g = gg, label = label)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we can plot diferences (relative error in grade)\n",
    "plt.plot (cutoff, gg[0]-gg[1], label = 'DGM - OK')\n",
    "plt.plot (cutoff, gg[0]-gg[2], label = 'DGM - ID2')\n",
    "plt.plot (cutoff, np.zeros(cutoff.shape[0]),'--k', label = 'Zero error')\n",
    "plt.title('relative error in grade')\n",
    "plt.xlabel ('cutoff')\n",
    "plt.ylabel ('relative error')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we can plot diferences (relative error in tonnage)\n",
    "plt.plot (cutoff, tt[0]-tt[1], label = 'DGM - OK')\n",
    "plt.plot (cutoff, tt[0]-tt[2], label = 'DGM - ID2')\n",
    "plt.plot (cutoff, np.zeros(cutoff.shape[0]),'--k', label = 'Zero error')\n",
    "plt.legend()\n",
    "plt.xlabel ('cutoff')\n",
    "plt.ylabel ('relative error')\n",
    "plt.title('relative error in tonnage')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
