gempy_single_layer.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import pandas as pd\n",
    "import pyvista as pv\n",
    "import pyvistaqt as pvqt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def external_kriging(G_position,G_orientation,layer_position,range):\n",
    "\n",
    "    G_1 = G_position\n",
    "    G_1_o = G_orientation\n",
    "    a_T = range\n",
    "    layer1 = layer_position\n",
    "    \n",
    "    def squared_euclidean_distance(x_1,x_2):\n",
    "        sqd = np.sqrt(np.reshape(np.sum(x_1**2,1),newshape =(x_1.shape[0],1))+\\\n",
    "        np.reshape(np.sum(x_2**2,1),newshape =(1,x_2.shape[0]))-\\\n",
    "        2*(x_1@x_2.T))\n",
    "        return sqd\n",
    "\n",
    "    def cartesian_dist(x_1,x_2):\n",
    "        return np.concatenate([\n",
    "            np.tile(x_1[:,0] - np.reshape(x_2[:,0],[x_2[:,0].shape[0],1]),[1,3]),\n",
    "            np.tile(x_1[:,1] - np.reshape(x_2[:,1],[x_2[:,1].shape[0],1]),[1,3]),\n",
    "            np.tile(x_1[:,2] - np.reshape(x_2[:,2],[x_2[:,2].shape[0],1]),[1,3])],axis = 0) \n",
    "\n",
    "    def cov_gradients(dist_tiled):\n",
    "        \n",
    "        condition1 = 0\n",
    "        a = (h_u*h_v)\n",
    "        b = dist_tiled**2\n",
    "\n",
    "        t1 =  np.divide(a, b, out=np.zeros_like(a),casting='unsafe', where=b!=0)\n",
    "        t2 = np.where(dist_tiled < a_T,(-c_o_T*((-14/a_T**2)+\n",
    "                    105/4*dist_tiled/a_T**3 -\n",
    "                    35/2 * dist_tiled**3 / a_T **5 +\n",
    "                    21 /4 * dist_tiled**5/a_T**7)),0)+\\\n",
    "            np.where(dist_tiled < a_T,(c_o_T * 7 * (9 * dist_tiled ** 5 -\n",
    "                        20 * a_T ** 2 * dist_tiled ** 3 +\n",
    "                        15 * a_T ** 4 * dist_tiled -\n",
    "                        4 * a_T ** 5) / (2 * a_T ** 7)),0)\n",
    "\n",
    "        t3 = perpendicularity_matrix * \\\n",
    "            np.where(dist_tiled < a_T,(c_o_T * ((-14 / a_T ** 2) + 105 / 4 * dist_tiled / a_T ** 3 -\n",
    "                    35 / 2 * dist_tiled ** 3 / a_T ** 5 +\n",
    "                    21 / 4 * dist_tiled ** 5 / a_T ** 7)),0)\n",
    "        t4 = 1/3*np.eye(dist_tiled.shape[0])\n",
    "\n",
    "        condition2 = t1 * t2 - t3 + t4\n",
    "\n",
    "        C_G = np.where(dist_tiled==0, condition1, condition2) \n",
    "        return C_G\n",
    "\n",
    "    def set_rest_ref_matrix(number_of_points_per_surface):\n",
    "        ref_layer_points = np.repeat(np.stack([layer1[-1]],axis = 0),repeats=number_of_points_per_surface-1,axis = 0)\n",
    "        rest_layer_points = np.concatenate([layer1[0:-1]],axis = 0)\n",
    "        return ref_layer_points,rest_layer_points\n",
    "\n",
    "    def cov_interface(ref_layer_points,rest_layer_points):\n",
    "        sed_rest_rest = squared_euclidean_distance(rest_layer_points,rest_layer_points)\n",
    "        sed_ref_rest = squared_euclidean_distance(ref_layer_points,rest_layer_points)\n",
    "        sed_rest_ref = squared_euclidean_distance(rest_layer_points,ref_layer_points)\n",
    "        sed_ref_ref = squared_euclidean_distance(ref_layer_points,ref_layer_points)\n",
    "        \n",
    "        C_I = c_o_T*(\\\n",
    "                    (np.where(sed_rest_rest < a_T, (1 - 7 * (sed_rest_rest / a_T) ** 2 +\\\n",
    "                    35 / 4 * (sed_rest_rest / a_T) ** 3 -\\\n",
    "                    7 / 2 * (sed_rest_rest / a_T) ** 5 +\\\n",
    "                    3 / 4 * (sed_rest_rest / a_T) ** 7),0)) -\\\n",
    "                    (np.where(sed_ref_rest < a_T, (1 - 7 * (sed_ref_rest / a_T) ** 2 +\\\n",
    "                    35 / 4 * (sed_ref_rest / a_T) ** 3 -\\\n",
    "                    7 / 2 * (sed_ref_rest / a_T) ** 5 +\\\n",
    "                    3 / 4 * (sed_ref_rest / a_T) ** 7),0)) -\\\n",
    "                    (np.where(sed_rest_ref < a_T, (1 - 7 * (sed_rest_ref / a_T) ** 2 +\\\n",
    "                    35 / 4 * (sed_rest_ref / a_T) ** 3 -\\\n",
    "                    7 / 2 * (sed_rest_ref / a_T) ** 5 +\\\n",
    "                    3 / 4 * (sed_rest_ref / a_T) ** 7),0)) +\\\n",
    "                    (np.where(sed_ref_ref < a_T, (1 - 7 * (sed_ref_ref / a_T) ** 2 +\\\n",
    "                    35 / 4 * (sed_ref_ref / a_T) ** 3 -\\\n",
    "                    7 / 2 * (sed_ref_ref / a_T) ** 5 +\\\n",
    "                    3 / 4 * (sed_ref_ref / a_T) ** 7),0)))\n",
    "        \n",
    "        return C_I\n",
    "\n",
    "    def cartesian_dist_no_tile(x_1,x_2):\n",
    "        return np.concatenate([\n",
    "            np.transpose((x_1[:,0] - np.reshape(x_2[:,0],[x_2.shape[0],1]))),\n",
    "            np.transpose((x_1[:,1] - np.reshape(x_2[:,1],[x_2.shape[0],1]))),\n",
    "            np.transpose((x_1[:,2] - np.reshape(x_2[:,2],[x_2.shape[0],1])))],axis = 0) \n",
    "\n",
    "    def cov_interface_gradients(hu_rest,hu_ref):\n",
    "        C_GI = hu_rest*np.where(sed_dips_rest < a_T,(- c_o_T * ((-14 / a_T ** 2) + 105 / 4 * sed_dips_rest / a_T ** 3 -\n",
    "                                    35 / 2 * sed_dips_rest ** 3 / a_T ** 5 +\n",
    "                                    21 / 4 * sed_dips_rest ** 5 / a_T ** 7)),0)-\\\n",
    "        hu_ref*np.where(sed_dips_ref < a_T,(-c_o_T * ((-14 / a_T ** 2) + 105 / 4 * sed_dips_ref / a_T ** 3 -\n",
    "                                    35 / 2 * sed_dips_ref ** 3 / a_T ** 5 +\n",
    "                                    21 / 4 * sed_dips_ref ** 5 / a_T ** 7)),0)\n",
    "        return C_GI\n",
    "\n",
    "    def safe_arange(start, stop, step):\n",
    "        return step * np.arange(start / step, stop / step)\n",
    "\n",
    "    def perpendicularity(G_1):\n",
    "        a = np.concatenate([np.ones([G_1.shape[0],G_1.shape[0]]),np.zeros([G_1.shape[0],G_1.shape[0]]),np.zeros([G_1.shape[0],G_1.shape[0]])],axis = 1)\n",
    "        b = np.concatenate([np.zeros([G_1.shape[0],G_1.shape[0]]),np.ones([G_1.shape[0],G_1.shape[0]]),np.zeros([G_1.shape[0],G_1.shape[0]])],axis = 1)\n",
    "        c = np.concatenate([np.zeros([G_1.shape[0],G_1.shape[0]]),np.zeros([G_1.shape[0],G_1.shape[0]]),np.ones([G_1.shape[0],G_1.shape[0]])],axis = 1)\n",
    "        return np.concatenate([a,b,c],axis = 0)    \n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "    c_o_T = a_T**2/14/3\n",
    "    G_1_tiled = np.tile(G_1,[3,1])\n",
    "\n",
    "    h_u = cartesian_dist(G_1,G_1)\n",
    "    h_v = h_u.T\n",
    "\n",
    "    perpendicularity_matrix = perpendicularity(G_1)\n",
    "\n",
    "    dist_tiled = squared_euclidean_distance(G_1_tiled,G_1_tiled)\n",
    "\n",
    "    dist_tiled = dist_tiled + np.eye(dist_tiled.shape[0])\n",
    "\n",
    "    C_G = cov_gradients(dist_tiled)\n",
    "\n",
    "    number_of_layer = 1\n",
    "\n",
    "    number_of_points_per_surface = np.array([layer1.shape[0]])\n",
    "\n",
    "    ref_layer_points,rest_layer_points = set_rest_ref_matrix(number_of_points_per_surface)\n",
    "\n",
    "    sed_rest_rest = squared_euclidean_distance(rest_layer_points,rest_layer_points)\n",
    "    sed_ref_rest = squared_euclidean_distance(ref_layer_points,rest_layer_points)\n",
    "    sed_rest_ref = squared_euclidean_distance(rest_layer_points,ref_layer_points)\n",
    "    sed_ref_ref = squared_euclidean_distance(ref_layer_points,ref_layer_points)\n",
    "\n",
    "    C_I = cov_interface(ref_layer_points,rest_layer_points)\n",
    "\n",
    "    sed_dips_rest = squared_euclidean_distance(G_1_tiled,rest_layer_points)\n",
    "    sed_dips_ref = squared_euclidean_distance(G_1_tiled,ref_layer_points)\n",
    "\n",
    "    hu_rest = cartesian_dist_no_tile(G_1,rest_layer_points)\n",
    "    hu_ref = cartesian_dist_no_tile(G_1,ref_layer_points)\n",
    "\n",
    "\n",
    "    C_GI = cov_interface_gradients(hu_rest,hu_ref)\n",
    "    C_IG = C_GI.T\n",
    "\n",
    "    xx = safe_arange(0,10,0.5)\n",
    "    yy = safe_arange(0,10,0.5)\n",
    "    zz = safe_arange(0,10,0.5)\n",
    "    XX,YY,ZZ = np.meshgrid(xx,yy,zz)\n",
    "    X = (np.reshape(XX,[-1])).T\n",
    "    Y = (np.reshape(YY,[-1])).T\n",
    "    Z = (np.reshape(ZZ,[-1])).T\n",
    "    grid = np.stack([X,Y,Z],axis = 1)\n",
    "\n",
    "\n",
    "    hu_Simpoints = cartesian_dist_no_tile(G_1,grid)\n",
    "    sed_dips_SimPoint = squared_euclidean_distance(G_1_tiled,grid)\n",
    "    sed_rest_SimPoint = squared_euclidean_distance(rest_layer_points,grid)\n",
    "    sed_ref_SimPoint = squared_euclidean_distance(ref_layer_points,grid)\n",
    "\n",
    "\n",
    "    K = np.concatenate([np.concatenate([C_G,C_GI],axis = 1),np.concatenate([C_IG,C_I],axis = 1)],axis = 0)\n",
    "\n",
    "    n = G_1.shape[0]\n",
    "    sub_x = np.tile(np.array([[1,0,0]]),[n,1])\n",
    "    sub_y = np.tile(np.array([[0,1,0]]),[n,1])\n",
    "    sub_z = np.tile(np.array([[0,0,1]]),[n,1])\n",
    "    sub_block1 = np.concatenate([sub_x,sub_y,sub_z],0)\n",
    "\n",
    "    sub_x_2 = np.reshape(1*G_1[:,0],[n,1])\n",
    "    sub_y_2 = np.reshape(1*G_1[:,1],[n,1])\n",
    "    sub_z_2 = np.reshape(1*G_1[:,2],[n,1])\n",
    "    sub_block2 = np.concatenate([sub_x_2, sub_y_2, sub_z_2], 0)\n",
    "\n",
    "    sub_xy = np.reshape(np.concatenate([G_1[:, 1],G_1[:, 0]], 0), [2 * n, 1])\n",
    "    sub_xy = np.pad(sub_xy, [[0, n], [0, 0]])\n",
    "\n",
    "    sub_xz = np.concatenate([np.pad(np.reshape(G_1[:, 2], [n, 1]), [[0, n], [0, 0]]), np.reshape(G_1[:, 0], [n, 1])], 0)\n",
    "\n",
    "    sub_yz = np.reshape(np.concatenate([G_1[:, 2], G_1[:, 1]], 0), [2 * n, 1])\n",
    "    sub_yz = np.pad(sub_yz, [[n, 0], [0, 0]])\n",
    "\n",
    "    sub_block3 = np.concatenate([sub_xy, sub_xz, sub_yz], 1)\n",
    "\n",
    "    U_G = np.concatenate([sub_block1, sub_block2, sub_block3], 1)\n",
    "\n",
    "    U_I = -np.stack([(rest_layer_points[:, 0] - ref_layer_points[:, 0]), \n",
    "                    (rest_layer_points[:, 1] - ref_layer_points[:, 1]),\n",
    "                    (rest_layer_points[:, 2] - ref_layer_points[:, 2]),\n",
    "                    (rest_layer_points[:, 0] ** 2 - ref_layer_points[:, 0] ** 2),\n",
    "                    (rest_layer_points[:, 1] ** 2 - ref_layer_points[:, 1] ** 2),\n",
    "                    (rest_layer_points[:, 2] ** 2 - ref_layer_points[:, 2] ** 2),\n",
    "                    (rest_layer_points[:, 0] * rest_layer_points[:, 1] - ref_layer_points[:, 0] * ref_layer_points[:, 1]),\n",
    "                    (rest_layer_points[:, 0] * rest_layer_points[:, 2] - ref_layer_points[:, 0] * ref_layer_points[:, 2]),\n",
    "                    (rest_layer_points[:, 1] * rest_layer_points[:, 2] - ref_layer_points[:, 1] * ref_layer_points[:, 2])], 1)\n",
    "\n",
    "    length_of_CG = C_G.shape[1]\n",
    "    length_of_CGI = C_GI.shape[1]\n",
    "\n",
    "    U_G = U_G[:length_of_CG, :3]\n",
    "    U_I = U_I[:length_of_CGI, :3]\n",
    "\n",
    "\n",
    "    U = np.concatenate([U_G,U_I],0)\n",
    "\n",
    "    # build zero matrix\n",
    "    zero_matrix = np.zeros([U.shape[1],U.shape[1]])\n",
    "\n",
    "\n",
    "    # concatenate drift matrix and kriging matrix\n",
    "    K_D = np.concatenate([np.concatenate([K,U],axis = 1),np.concatenate([U.T,zero_matrix],axis = 1)],axis = 0)\n",
    "\n",
    "    # build right side matrix of cokriging system\n",
    "    bk = np.concatenate([G_1_o[:,0],G_1_o[:,1],G_1_o[:,2],np.zeros(K_D.shape[0]-G_1.shape[0]*3)],axis = 0)\n",
    "    bk = np.reshape(bk,newshape = [bk.shape[0],1])\n",
    "\n",
    "    # solve kriging weight\n",
    "    w = np.linalg.lstsq(K_D,bk)[0]\n",
    "\n",
    "    # gradient contribution\n",
    "    sigma_0_grad = w[:G_1.shape[0]*3] * (-hu_Simpoints*(sed_dips_SimPoint < a_T)*(- c_o_T * ((-14 / a_T ** 2) + 105 / 4 * sed_dips_SimPoint / a_T ** 3 - 35 / 2 * sed_dips_SimPoint ** 3 / a_T ** 5 + 21 / 4 * sed_dips_SimPoint ** 5 / a_T ** 7)))\n",
    "\n",
    "    sigma_0_grad = np.sum(sigma_0_grad,axis=0)\n",
    "    # surface point contribution\n",
    "    sigma_0_interf = -w[G_1.shape[0]*3:-3]*(c_o_T  * ((sed_rest_SimPoint < a_T)*(1 - 7 * (sed_rest_SimPoint / a_T) ** 2 +\n",
    "                35 / 4 * (sed_rest_SimPoint / a_T) ** 3 -\n",
    "                7 / 2 * (sed_rest_SimPoint / a_T) ** 5 +\n",
    "                3 / 4 * (sed_rest_SimPoint / a_T) ** 7) -\n",
    "                (sed_ref_SimPoint < a_T)*(1 - 7 * (sed_ref_SimPoint / a_T) ** 2 +\n",
    "                35 / 4 * (sed_ref_SimPoint / a_T) ** 3 -\n",
    "                7 / 2 * (sed_ref_SimPoint / a_T) ** 5 +\n",
    "                3 / 4 * (sed_ref_SimPoint / a_T) ** 7)))\n",
    "    sigma_0_interf = np.sum(sigma_0_interf,axis = 0)\n",
    "\n",
    "    # 2nd order drift contribution\n",
    "    sigma_0_2nd_drift_1 = np.sum(grid* (w[-3:]).T,axis = 1)\n",
    "\n",
    "    sigma_0_2nd_drift = sigma_0_2nd_drift_1\n",
    "\n",
    "\n",
    "    interpolate_result = 1*sigma_0_grad+1*sigma_0_interf + sigma_0_2nd_drift \n",
    "\n",
    "    intp = np.reshape(interpolate_result,[20,20,20]) # reshape the result to matrix shape\n",
    "    intp_c = intp\n",
    "\n",
    "    def plot_3D(grid,value,mode='block',surfaces_nr=10):\n",
    "        x = grid[:,0]\n",
    "        y = grid[:,1]\n",
    "        z = grid[:,2]\n",
    "\n",
    "        grid = pv.UniformGrid(dims=(20,20,20), spacing=(0.5,0.5,0.5),origin=(0,0,0))\n",
    "        p = pvqt.BackgroundPlotter()\n",
    "        layer1_ = layer_position\n",
    "        if mode == 'block':\n",
    "            p.add_mesh(grid,scalars=value.flatten(),show_edges=False,style='surface')\n",
    "            p.add_axes()\n",
    "            p.show()\n",
    "\n",
    "        elif mode == 'surfaces':\n",
    "            grid['vol'] = value.ravel()\n",
    "            lvl_t=value[np.where((XX==layer1[0,0])&(YY==layer1[0,1])&(ZZ==layer1[0,2]))]\n",
    "            contours_1 = grid.contour(np.linspace(value.min(),value.max(),surfaces_nr))\n",
    "            contours_2 = grid.contour(lvl_t)\n",
    "            #p.add_mesh(contours_1, show_scalar_bar=True, label='scalar_field',style='surface')\n",
    "            p.add_mesh(contours_2, show_scalar_bar=True, label='scalar_field',style='surface',color='red')\n",
    "        p.add_points(layer1[:,[2,0,1]], render_points_as_spheres=True,\n",
    "                        point_size=20.0)\n",
    "        \n",
    "        orientation1 = pv.Arrow(start=(G_1[:,[2,0,1]][0]),direction=(G_1_o[:,[2,0,1]][0]))\n",
    "        orientation2 = pv.Arrow(start=(G_1[:,[2,0,1]][1]),direction=(G_1_o[:,[2,0,1]][1]))\n",
    "        p.add_mesh(orientation1)\n",
    "        p.add_mesh(orientation2)\n",
    "        p.add_axes()\n",
    "        p.show()\n",
    "\n",
    "    plot_3D(grid,intp,mode='surfaces',surfaces_nr=5)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "G_position = np.array([[5,6,0],[5,4,0]])\n",
    "G_orientation = np.array([[0,1,0],[0,1,0]])\n",
    "range = 10\n",
    "layer_position = np.array([[8.5,2,0],[7,4,0],[5.5,5.4,0],[4,5.4,0],[2.5,4,0],[1,2,0]])\n",
    "\n",
    "external_kriging(G_position,G_orientation,layer_position,range)"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "32398c9479ff240dde49d7b1f724314064821a1f7902a5fffc4bc8ead0815669"
  },
  "kernelspec": {
   "display_name": "Python 3.9.12 ('gempy_test')",
   "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.9.12"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}