{"metadata":{"accelerator":"GPU","colab":{"gpuType":"A100","machine_shape":"hm","provenance":[]},"kernelspec":{"display_name":"Python 3","name":"python3"},"language_info":{"name":"python"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":87230,"databundleVersionId":9911932,"sourceType":"competition"}],"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# This notebook generates the synthetic data. It was developed for Google Colab!!\n# Mount Google Drive.\n# Only for Google Colab\nfrom google.colab import drive\ndrive.mount('/content/drive')","metadata":{"id":"ydEZGDKO_gH3","outputId":"342e5de6-d648-4012-fb96-4883bb31ace4","pycharm":{"name":"#%%\n"}},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Install conda in Google colab\n# Only for Google COlab\n!pip install -q condacolab\nimport condacolab\ncondacolab.install()\ncondacolab.check()","metadata":{"id":"cvJx1SR-_IJ-","outputId":"8c61a052-d77d-4c85-bca8-77f499041253","pycharm":{"name":"#%%\n"}},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Download gprMax\n!git clone https://github.com/gprMax/gprMax.git\n\n# !!! Very Important part !!!!\n# Go to gprMax folder, at materials.py, and change the last function at the end of the page to the following function\n# ----------------------------------------------------------------------------------------------------------\n# def calculate_debye_properties(self, nbins, G, fractalboxname):\n\n#         self.startmaterialnum = len(G.materials)\n\n#         mubins = np.linspace(self.mu[0], self.mu[1], nbins)\n#         mumaterials = mubins + (mubins[1] - mubins[0]) / 2\n\n#         muiter = np.nditer(mumaterials, flags=['c_index'])\n#         while not muiter.finished:\n#             er = muiter[0]\n\n#             sig = er/10000\n\n#             digitscount =  len(str(int(nbins)))\n#             materialID = '|{}_{}|'.format(fractalboxname, str(muiter.index + 1).zfill(digitscount))\n#             m = Material(len(G.materials), materialID)\n#             m.type = ''\n#             m.averagable = False\n\n#             m.er = er\n#             m.se = sig\n\n#             G.materials.append(m)\n\n#             muiter.iternext()\n\n# ----------------------------------------------------------------------------------------------------------\n\n#  This changes the peplinski soil and instead of minimum and maximum water fraction, now it takes as inputs minimum and maximum permittivity.\n#  The conductivity is chosen automatically as er/1000 i.e. low conductivity environment (simialr to planetary setups).\n\n","metadata":{"id":"qk54F-fb_dh8","pycharm":{"name":"#%%\n"}},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Install pycude for running gprMax with GPUs'\n!pip install pycuda","metadata":{"id":"RMGrC1uWAD0v","outputId":"5fe9fbe5-cf29-4378-e53d-e5119880e7cc","pycharm":{"name":"#%%\n"},"_kg_hide-output":false,"_kg_hide-input":false,"scrolled":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os\n# This is the directory where gprMax is saved.\n# The directory is for Google Colab, you need to change accordingly.\nos.chdir('/content/drive/MyDrive/gprMax')\n\n\n# Update conda\n!mamba env update -n base -f conda_env.yml\n# Remove gprMax environments\n!conda remove --name gprMax\n# Build and install gprMax\n!conda activate gprMax\n!python setup.py build\n!python setup.py install","metadata":{"id":"vXo4WQrNAzU1","outputId":"2c2adc5f-f74a-4753-deaa-b51cebb32fc4","pycharm":{"name":"#%%\n"},"scrolled":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np\n\n\ndef write_file(fil):\n    # Function that writes the input gprMax file.\n    # fil: Name of the file\n\n    # File will be saved in /content/drive/MyDrive\n    # This is the default for Google Colab.\n    # If you are not using Google Colab you need to change accordingly\n    f = open('/content/drive/MyDrive/' + fil, 'w')\n\n    # Define the domain of the problem\n    # 2D with 50.02x26 dimensions.\n    # z dimension equals with spatial step i.e. 2D geometry\n    f.write('#domain: 50.02 26 0.02 \\n')\n\n    # Choosing the spatial step which is 2 cm\n    f.write('#dx_dy_dz: 0.02 0.02 0.02 \\n')\n\n    # Simulation will run for 15000 iterations\n    f.write('#time_window: 15000 \\n')\n    f.write('#time_step_stability_factor: 0.99 \\n')\n\n    # In every realisation of the function the model will choose a random shape for the pulse\n    wave = {0: 'gaussiandot', 1: 'ricker', 2: 'gaussiandotnorm', 3: 'gaussiandotdot'}\n    # Randomly choose the index of the pulse\n    ind = np.random.randint(4)\n\n    # Randomly choose the central frequency of the pulse from 60 to 100 MHz\n    cf = 60+40*np.random.rand()\n\n    # Define the waveform with the randomly selected shape and central frequency\n    f.write('#waveform: {} 1 {}e6 my_pulse \\n'.format( wave[ind], cf ));\n\n    # Place your transmitters at y = 22.8 meters.\n    # If you use one transmiter it will be one-shot configuration\n    f.write('#hertzian_dipole: z 2 22.8 0 my_pulse \\n');\n    # f.write('#hertzian_dipole: z 48 22.8 0 my_pulse \\n');\n    # f.write('#hertzian_dipole: z 25 22.8 0 my_pulse \\n');\n\n    # Place your receivers from 101*0.02 with step 10*0.02 all the way to 2400*0.02\n    for i in range(101,2400,10):\n        f.write('#rx: {} 22.8 0 \\n'.format(i*0.02) );\n\n    # We use 60 cells for the PML to ensure no artifacts from the boundaries\n    f.write('#pml_cells: 60 60 0 60 60 0 \\n');\n    # You can play with the PML parameters and reduce the number of PML cells\n    # f.write('#pml_cfs: constant reverse 0.015 0 constant forward 1 5 quadratic forward 0 None \\n');\n\n\n\n    ###### First layer #########################\n    # Choose the fractal dimension of the first layer\n    d = 0.5 + 3*np.random.rand()\n\n    # Choose how assymetrical the fractal distribution will be\n    dx = 0.5 + 1*np.random.rand()\n\n    # Choose the lower permittivity of the formation randomly from 1.5 to 9.5 (similar to planetary setups)\n    r1 = 1.5 + 8*np.random.rand()\n\n    # Choose the upper permittivity of the formation randomly from r1 to 10\n    r2 = r1 + (10-r1)*np.random.rand()\n\n    # z1 and z2 and the lower and upper y-axis of the formation. The first formation extents from y = 0-22\n    # Subsequent layers will be written on tom of the first layer.\n    z1 = 0\n    z2 = 22\n\n    # Choose the discretisation bins of the first formation\n    bin = 3 + 37*np.random.rand()\n\n    # pro_1 is the ID of the soil for the first formation\n    pro = \"pro_{}\".format(1)\n\n    # geo_1 is the ID of the fractal box for the first layer\n    geo = \"geo_{}\".format(1)\n\n    # Define the soil properties for the first layer\n    f.write('#soil_peplinski: 0.5 0.5 2 2.66 {} {} {} \\n'.format(r1, r2, pro))\n    f.write('#fractal_box: 0 {} 0 50.02 {} 0.02 {} {} 1 1 {} {} {} \\n'.format(z1, z2, d, dx, int(np.round(bin)), pro, geo))\n\n    # Select the fractal dimension of the surface of the first formation\n    sd = 0.9 + 2*np.random.rand()\n\n    # Select the minimum and maximum amplitude of the first layer.\n    # The first layer is chosen to be relatively flat i.e. have a topography variation of 20 cm\n    t1 = 0.2\n    f.write('#add_surface_roughness: 0 {} 0 50.02 {} 0.02 {} 1 1 {} {} {} \\n'.format(z2, z2, sd, z2 - t1, z2 + t1, geo))\n    #########################\n\n\n    for i in range(2,200000):\n        # This for-lop will generate layers one after the other until it reaches the limits of the model\n\n        # Choose the fractal dimension of the layer\n        d = 0.5 + 3*np.random.rand()\n\n        # Choose the assymetry factor for the layer\n        dx = 0.5 + 1*np.random.rand()\n\n        # Choose the minimum permittivity\n        r1 = 1.5 + 8*np.random.rand()\n\n        # Choose the maximum permittivity\n        r2 = r1 + (10-r1)*np.random.rand()\n\n        # Choose the upper depth of the layer. The lower depth is always zero i.e. the next layer is written on the previous one\n        # The upper depth is updated based on the previous one i.e. iteration by iteration the depth for every layer decreases.\n        z2 = z2 - 10*np.random.rand()\n\n        # Make the depth a multiple of the spatial step i.e. 0.02.\n        z2 = 0.02 * np.round(z2/0.02)\n        z1 = 0\n\n        # Choose the number of discritised bins\n        bin = 3 + 37*np.random.rand()\n\n        # Choose the maximum topographic variation of the topography of hte layer. Randomly from 0.2 - 1.7 i.e. maximum 3.4 m\n        t1 = 0.2 + 1.5*np.random.rand()\n\n        # Make t1 a multiple of 0.02\n        t1 = np.round(t1/0.02)*0.02\n\n        # If the upper depth of the layer is below 1 meter then stop making more layers. Notice that top surface is at 22 meters.\n        if z2 - t1 < 1:\n          break\n\n        # pro_{i} is the ID of the ith soil parameters\n        pro = \"pro_{}\".format(i)\n\n        # geo_{i} is the ID of the ith fractal box\n        geo = \"geo_{}\".format(i)\n\n        # Define the ith layer\n        f.write('#soil_peplinski: 0.5 0.5 2 2.66 {} {} {} \\n'.format(r1, r2, pro))\n        f.write('#fractal_box: 0 {} 0 50.02 {} 0.02 {} {} 1 1 {} {} {} \\n'.format(z1, z2, d, dx, int(np.round(bin)), pro, geo))\n\n        # Choose the fractal dimension of the topography of the ith layer\n        sd = 0.9 + 2*np.random.rand()\n\n        # Build the topographu for the ith layer\n        f.write('#add_surface_roughness: 0 {} 0 50.02 {} 0.02 {} 1 1 {} {} {} \\n'.format(z2, z2, sd, z2 - t1, z2 + t1, geo))\n\n    ##############\n\n    # In this saves the model parameters from x,y,z = (0,0,0) to x,y,z = (50, 22, 0.02) i.e. it saves the permittivity of the model exluding the free space\n    # We need this file in order to get the ground truth for training ML\n    f.write(\"#geometry_objects_write: 0 0 0 50 22 0.02 file_geo \\n\")\n\n\n    # f.write('#geometry_view: 0 0 0 50 25 0.02 0.04 0.04 0.02 geo_ml n \\n');\n\n\n\n\n\n    f.close()\n\n","metadata":{"id":"5QcF86MFwry1","pycharm":{"name":"#%%\n"}},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# gprMax gives error when using Numpy 2\n!pip uninstall numpy -y\n!pip install numpy==1.22.1","metadata":{"id":"6Hcz6w2dyFCD","outputId":"b62dd337-5a4d-4068-d1fe-0e02b0b5587e","pycharm":{"name":"#%%\n"},"scrolled":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Choose the directory where you want the Bscans and Ground-Truth nunpy arrays to be saved\nos.chdir('/content/drive/MyDrive')\n\nfrom matplotlib import pyplot as plt\nfrom skimage.transform import resize\nimport numpy as np\nfrom matplotlib import pyplot as pl\n\nimport argparse\nimport glob\nimport os\nimport copy\nimport h5py\nimport numpy as np\nfrom sklearn import preprocessing as p\n\nfrom scipy.interpolate import NearestNDInterpolator\n\n\n# This is a gprMax function that is used to open gprMax output files\ndef get_output_data(filename, rxnumber, rxcomponent):\n    \"\"\"Gets B-scan output data from a model.\n\n    Args:\n        filename (string): Filename (including path) of output file.\n        rxnumber (int): Receiver output number.\n        rxcomponent (str): Receiver output field/current component.\n\n    Returns:\n        outputdata (array): Array of A-scans, i.e. B-scan data.\n        dt (float): Temporal resolution of the model.\n    \"\"\"\n\n    # Open output file and read some attributes\n    f = h5py.File(filename, 'r')\n    nrx = f.attrs['nrx']\n    dt = f.attrs['dt']\n\n    # Check there are any receivers\n    if nrx == 0:\n        raise CmdInputError('No receivers found in {}'.format(filename))\n\n    path = '/rxs/rx' + str(rxnumber) + '/'\n    availableoutputs = list(f[path].keys())\n\n    # Check if requested output is in file\n    if rxcomponent not in availableoutputs:\n        raise CmdInputError('{} output requested to plot, but the available output for receiver 1 is {}'.format(rxcomponent, ', '.join(availableoutputs)))\n\n    outputdata = f[path + '/' + rxcomponent]\n    outputdata = np.array(outputdata)\n    f.close()\n\n    return outputdata, dt\n\n\n\ndef make_bscan(fil,plot=False):\n    # The function opens .out file, process the BScan (gain) and saves it in a numpy array\n    # plot: When True it plots the BScan after saturated its values for visual inspection\n    # fil: The name of the .out file\n\n\n    # Number of receivers\n    # This is for the case study examined in this script. If you decide to change the numnber of receivers in the write_file, then you need to ...\n    # change the number of receivers here as well.\n    n_receivers = np.shape(np.arange(101,2400,10))[0]\n\n    Bscan=[]\n    for i in range(0,n_receivers):\n        # Go through all the receivers from the fil\n        [fi, t]=get_output_data(fil, i+1, \"Ez\")\n\n        # Apply gain. You can change this, or you can remove it completely\n        gain = np.arange(0,np.shape(fi)[0],40)**3\n\n        Bscan.append(fi[0:np.shape(fi)[0]:40]*gain)\n        # plt.plot(fi[0:np.shape(fi)[0]:40])\n        # plt.show()\n\n    Bs = np.array(Bscan)\n    image = Bs.T\n\n    # Resize your BScan to 230x230 dimensions.\n    B = resize(image, (230, 230))\n    B2 = copy.copy(B)\n    if plot:\n        B[B>np.max(B)*0.05] = np.max(B)*0.05\n        B[B<np.min(B)*0.05] = np.min(B)*0.05\n        plt.imshow(B,aspect='auto',cmap= 'bone')\n        plt.show()\n\n    # B2 is a numpy array with the resized (230x230) and processed BScan\n    return B2\n\n\n\ndef make_ground_truth(file_name, plot=False):\n    # This function opens the file_name to extract the ground truth for every model\n\n\n    # Id is a 2D array with the IDs of the materials. The IDs of the materials are saved in the file\n    # file_geo_materials.txt\n    f = h5py.File(file_name, \"r\")\n    n=np.array(f['data'])\n    Idf = np.array(n[:,:,0])\n    Id = np.zeros(np.shape(Idf))\n\n\n    # Open the materials file, and isolate the permittivity associated with every ID\n    f = open('file_geo_materials.txt')\n    mat = []\n    for c, x in enumerate(f):\n        paok = x.split(\" \")\n        mat.append(float(paok[1]))\n        # Replace the IDs with their permittivitie associated with these IDs\n        Id[Idf == c] = float(paok[1])\n\n\n    # Reduce the discretisation undersample by 8 and 4 the x and y axis respectively\n    km = np.shape(Id)\n    Id2 = Id[0:km[0]:8, 0:km[1]:4]\n\n    Id = np.flipud(Id.T)\n    Id2 = np.flipud(Id2.T)\n\n\n\n\n\n\n\n    # min_max_scaler = p.MinMaxScaler()\n    # normalizedData = min_max_scaler.fit_transform(Id2)\n    # Id2 = mm + normalizedData*(MM-mm)\n\n\n    if plot:\n        mm = np.min(np.array(mat[2:]))\n        MM = np.max(np.array(mat[2:]))\n        plt.imshow(Id,vmin=mm, vmax=MM)\n        plt.colorbar()\n        plt.show()\n\n        plt.imshow(Id2,vmin=mm, vmax=MM)\n        plt.colorbar()\n        plt.show()\n\n\n    # B2 is a numpy array with the resized ground truth\n    return Id2\n\n\n\n\ndef lookupNearest(x0, y0, x, y, data):\n    # Function for reshaping the ground truth based on nearest point\n    xi = np.abs(x-x0).argmin()\n    yi = np.abs(y-y0).argmin()\n    return data[xi,yi]\n\n\n# The dimensions of the ground truth 275x313\nx = np.linspace(0, 1, 275)\ny = np.linspace(0, 1, 313)\n\n# Resizing it to be square i.e. 224x224\nx_new = np.linspace(0, 1, 224)\ny_new = np.linspace(0, 1, 224)\n\n\nfor i in range(0,4581):\n    # Generate 5000 models\n\n    # Generate an input file with random layers and permittivity distributions\n    write_file('test_test.in')\n\n    # Run gprMax for the test_test.in file generated before\n    !python -m gprMax test_test.in -gpu\n\n    # Call make_bscan to read the .out file and save the Bscan in the numpy array B2\n    B2 = make_bscan('test_test.out', plot=False)\n\n    # Save B2 in the folder Bscan_data_2 and the file Inversion_Bscan_2_{i}\n    np.save(\"Kaggle_Bscan/Bscan_{}\".format(i), B2)\n\n    # Call make_ground_truth to open file_geo_10.h5 to save the ground truth permittivity in the numpy array Id2\n    temporary_gd = make_ground_truth(file_name='file_geo.h5', plot=False)\n\n    Id2 = np.zeros((224,224))\n    # Interpolate temporary_gd to resize the ground truth to a 224x224\n    for ax, q in enumerate(x_new):\n        for ay, w in enumerate(y_new):\n          Id2[ax, ay] = lookupNearest(q, w, x, y, temporary_gd)\n\n\n    # Save B2 in the folder Ground_truth_2 and the file Inversion_Model_2_{i}\n    np.save(\"Kaggle_Labels/Model_{}\".format(i), Id2)\n\n\n# This jupyter notebook will generate 4580 models with their resized BScans and corresponding permittivity ground truth\n# The Bscans will be saved at the folder Kaggle_Bscan\n# Ground truth will be saved at the folder Kaggle_Labels\n# All images will be saved as numpy arrays i.e. .npy\n\n","metadata":{"id":"F-3z4-KQ4_9V","outputId":"e63b2b4f-09fb-41d7-ddb0-ac81554e23ca","pycharm":{"name":"#%%\n"},"scrolled":true},"execution_count":null,"outputs":[]}]}