<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Sandstone Model
First we make a GeMpy instance with most of the parameters default (except range that is given by the project). Then we also fix the extension and the resolution of the domain we want to interpolate. Finally we compile the function, only needed once every time we open the project (the guys of theano they are working on letting loading compiled files, even though in our case it is not a big deal).
General note. So far the reescaling factor is calculated for all series at the same time. GeoModeller does it individually for every potential field. I have to look better what this parameter exactly means
Step1:  All input data is stored in pandas dataframes under, self.Data.Interances and self.Data.Foliations
Step2:  Plotting raw data
Step3:  Class Interpolator
This class will take the data from the class Data and calculate potential fields and block
Step4:  Early granite
Step5:  BIF Series
Step6:  SImple mafic
Step7:  Optimizing the export of lithologies
Here I am going to try to return from the theano interpolate function the internal type of the result (in this case DK I guess) so I can make another function in python I guess to decide which potential field I calculate at every grid_pos
Step9:  Export vtk
Step10:  Performance Analysis
CPU
Step11:  GPU
<ASSISTANT_TASK:>
Python Code:
# Setting extend, grid and compile
# Setting the extent
sandstone = GeMpy_core.GeMpy()
# Create Data class with raw data
sandstone.import_data( 696000,747000,6863000,6950000,-20000, 2000,
                                     path_f = os.pardir+"/input_data/a_Foliations.csv",
                                     path_i = os.pardir+"/input_data/a_Points.csv")
sandstone.Data.set_series({"EarlyGranite_Series":sandstone.Data.formations[-1], 
                      "BIF_Series":(sandstone.Data.formations[0], sandstone.Data.formations[1]),
                      "SimpleMafic_Series":sandstone.Data.formations[2]}, 
                       order = ["EarlyGranite_Series",
                              "BIF_Series",
                              "SimpleMafic_Series"]) 
Explanation: Sandstone Model
First we make a GeMpy instance with most of the parameters default (except range that is given by the project). Then we also fix the extension and the resolution of the domain we want to interpolate. Finally we compile the function, only needed once every time we open the project (the guys of theano they are working on letting loading compiled files, even though in our case it is not a big deal).
General note. So far the reescaling factor is calculated for all series at the same time. GeoModeller does it individually for every potential field. I have to look better what this parameter exactly means
End of explanation
sandstone.Data.Foliations;
# Create a class Grid so far just regular grid
sandstone.create_grid()
sandstone.Grid.grid
Explanation: All input data is stored in pandas dataframes under, self.Data.Interances and self.Data.Foliations:
End of explanation
sandstone.Plot.plot_data(serie = sandstone.Data.series.columns.values[1])
Explanation: Plotting raw data
End of explanation
sandstone.set_interpolator()
Explanation: Class Interpolator
This class will take the data from the class Data and calculate potential fields and block
End of explanation
sandstone.Plot.plot_potential_field(10, n_pf=0)
Explanation: Early granite
End of explanation
sandstone.Plot.plot_potential_field(10, n_pf=1, cmap = "magma",  plot_data = False)
Explanation: BIF Series
End of explanation
sandstone.Plot.plot_potential_field(10, n_pf=2)
Explanation: SImple mafic
End of explanation
# Reset the block
sandstone.Interpolator.block.set_value(np.zeros_like(sandstone.Grid.grid[:,0]))
# Compute the block
sandstone.Interpolator.compute_block_model([0,1,2], verbose = 1)
sandstone.Plot.plot_block_section()
Explanation: Optimizing the export of lithologies
Here I am going to try to return from the theano interpolate function the internal type of the result (in this case DK I guess) so I can make another function in python I guess to decide which potential field I calculate at every grid_pos
End of explanation
Export model to VTK
Export the geology blocks to VTK for visualisation of the entire 3-D model in an
external VTK viewer, e.g. Paraview.
..Note:: Requires pyevtk, available for free on: https://github.com/firedrakeproject/firedrake/tree/master/python/evtk
**Optional keywords**:
    - *vtk_filename* = string : filename of VTK file (default: output_name)
    - *data* = np.array : data array to export to VKT (default: entire block model)
vtk_filename = "noddyFunct2"
extent_x = 10
extent_y = 10
extent_z = 10
delx = 0.2
dely = 0.2
delz = 0.2
from pyevtk.hl import gridToVTK
# Coordinates
x = np.arange(0, extent_x + 0.1*delx, delx, dtype='float64')
y = np.arange(0, extent_y + 0.1*dely, dely, dtype='float64')
z = np.arange(0, extent_z + 0.1*delz, delz, dtype='float64')
# self.block = np.swapaxes(self.block, 0, 2)
gridToVTK(vtk_filename, x, y, z, cellData = {"geology" : sol})
Explanation: Export vtk
End of explanation
%%timeit
sol =  interpolator.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[0]
sandstone.block_export.profile.summary()
Explanation: Performance Analysis
CPU
End of explanation
%%timeit
# Reset the block
sandstone.block.set_value(np.zeros_like(sandstone.grid[:,0]))
# Compute the block
sandstone.compute_block_model([0,1,2], verbose = 0)
sandstone.block_export.profile.summary()
Explanation: GPU
End of explanation
<END_TASK>