In [27]:

    

import importlib
import theano.tensor as T
import sys, os
sys.path.append("/home/bl3/PycharmProjects/GeMpy/GeMpy")
sys.path.append("/home/bl3/PycharmProjects/GeMpy/pygeomod/pygeomod")
sys.path.append("/home/miguel/PycharmProjects/GeMpy/GeMpy")
import GeoMig
#import geogrid
#importlib.reload(GeoMig)
importlib.reload(GeoMig)
import pandas as pn

import numpy as np

os.environ['CUDA_LAUNCH_BLOCKING'] = '1'
np.set_printoptions(precision = 4, linewidth= 300, suppress =  True)
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
from matplotlib import cm


    

In [ ]:

    




    

In [28]:

    

# Setting extend, grid and compile
# Setting the extent
test = GeoMig.Interpolator(0,10,0,10,0,10,
                                range_var = np.float32(17),
                               u_grade = 3) # Range used in geomodeller

# Setting resolution of the grid
test.set_resolutions(40,40,60)
test.create_regular_grid_3D()

# Compiling
test.theano_compilation_3D()
#test.theano_set_3D_nugget_degree0()


    

    




/home/bl3/anaconda3/lib/python3.5/site-packages/ipykernel/__main__.py:12: DeprecationWarning: stack(*tensors) interface is deprecated, use stack(tensors, axis=0) instead.
WARNING (theano.gof.cmodule): Removing key file /home/bl3/.theano/compiledir_Linux-3.13--generic-x86_64-with-debian-jessie-sid-x86_64-3.5.2-64/tmp84w1cw07/key.pkl because the corresponding module is gone from the file system.
WARNING (theano.gof.cmodule): Removing key file /home/bl3/.theano/compiledir_Linux-3.13--generic-x86_64-with-debian-jessie-sid-x86_64-3.5.2-64/tmp0os76vum/key.pkl because the corresponding module is gone from the file system.
WARNING (theano.gof.cmodule): Removing key file /home/bl3/.theano/compiledir_Linux-3.13--generic-x86_64-with-debian-jessie-sid-x86_64-3.5.2-64/tmpaazz4wrf/key.pkl because the corresponding module is gone from the file system.
WARNING (theano.gof.cmodule): Removing key file /home/bl3/.theano/compiledir_Linux-3.13--generic-x86_64-with-debian-jessie-sid-x86_64-3.5.2-64/tmp9jmlp1tr/key.pkl because the corresponding module is gone from the file system.
WARNING (theano.gof.cmodule): Removing key file /home/bl3/.theano/compiledir_Linux-3.13--generic-x86_64-with-debian-jessie-sid-x86_64-3.5.2-64/tmpaibd3kv6/key.pkl because the corresponding module is gone from the file system.
WARNING (theano.gof.cmodule): Removing key file /home/bl3/.theano/compiledir_Linux-3.13--generic-x86_64-with-debian-jessie-sid-x86_64-3.5.2-64/tmpg26046o7/key.pkl because the corresponding module is gone from the file system.
WARNING (theano.gof.cmodule): Removing key file /home/bl3/.theano/compiledir_Linux-3.13--generic-x86_64-with-debian-jessie-sid-x86_64-3.5.2-64/tmp_komhyop/key.pkl because the corresponding module is gone from the file system.
WARNING (theano.gof.cmodule): Removing key file /home/bl3/.theano/compiledir_Linux-3.13--generic-x86_64-with-debian-jessie-sid-x86_64-3.5.2-64/tmp83ew0tyo/key.pkl because the corresponding module is gone from the file system.
WARNING (theano.gof.cmodule): Removing key file /home/bl3/.theano/compiledir_Linux-3.13--generic-x86_64-with-debian-jessie-sid-x86_64-3.5.2-64/tmp7v_5ct36/key.pkl because the corresponding module is gone from the file system.
WARNING (theano.gof.cmodule): Removing key file /home/bl3/.theano/compiledir_Linux-3.13--generic-x86_64-with-debian-jessie-sid-x86_64-3.5.2-64/tmp3agwgj_t/key.pkl because the corresponding module is gone from the file system.
WARNING (theano.gof.cmodule): Removing key file /home/bl3/.theano/compiledir_Linux-3.13--generic-x86_64-with-debian-jessie-sid-x86_64-3.5.2-64/tmpecl50v8n/key.pkl because the corresponding module is gone from the file system.
WARNING (theano.gof.cmodule): Removing key file /home/bl3/.theano/compiledir_Linux-3.13--generic-x86_64-with-debian-jessie-sid-x86_64-3.5.2-64/tmpdbgcb2cy/key.pkl because the corresponding module is gone from the file system.
WARNING (theano.gof.cmodule): Removing key file /home/bl3/.theano/compiledir_Linux-3.13--generic-x86_64-with-debian-jessie-sid-x86_64-3.5.2-64/tmp5g693w2z/key.pkl because the corresponding module is gone from the file system.
WARNING (theano.gof.cmodule): Removing key file /home/bl3/.theano/compiledir_Linux-3.13--generic-x86_64-with-debian-jessie-sid-x86_64-3.5.2-64/tmpg2otpfgc/key.pkl because the corresponding module is gone from the file system.
WARNING (theano.gof.cmodule): Removing key file /home/bl3/.theano/compiledir_Linux-3.13--generic-x86_64-with-debian-jessie-sid-x86_64-3.5.2-64/tmpdft27hal/key.pkl because the corresponding module is gone from the file system.
WARNING (theano.gof.cmodule): Removing key file /home/bl3/.theano/compiledir_Linux-3.13--generic-x86_64-with-debian-jessie-sid-x86_64-3.5.2-64/tmp9hqbudxq/key.pkl because the corresponding module is gone from the file system.
WARNING (theano.gof.cmodule): Removing key file /home/bl3/.theano/compiledir_Linux-3.13--generic-x86_64-with-debian-jessie-sid-x86_64-3.5.2-64/tmpf78eos80/key.pkl because the corresponding module is gone from the file system.
WARNING (theano.gof.cmodule): Removing key file /home/bl3/.theano/compiledir_Linux-3.13--generic-x86_64-with-debian-jessie-sid-x86_64-3.5.2-64/tmp8l9h5fp2/key.pkl because the corresponding module is gone from the file system.
WARNING (theano.gof.cmodule): Removing key file /home/bl3/.theano/compiledir_Linux-3.13--generic-x86_64-with-debian-jessie-sid-x86_64-3.5.2-64/tmp2g53skc5/key.pkl because the corresponding module is gone from the file system.
WARNING (theano.gof.cmodule): Removing key file /home/bl3/.theano/compiledir_Linux-3.13--generic-x86_64-with-debian-jessie-sid-x86_64-3.5.2-64/tmpgcj5ki7p/key.pkl because the corresponding module is gone from the file system.

In [41]:

    

layer_1 = np.array([[0.5,4,7], [2,4,6], [4,4,7], [5,4,5], [8,4,7], [7,4,8], [1,5,7]])#-np.array([5,5,4]))/8+0.5
layer_2 = np.array([[3,5,2], [6,5,2]])
layer_3 = np.array([[1,5,1],[1,5,2],[6,5,3]])#- np.array([5,5,4]))/8+0.5

dip_pos_1 = np.array([5,5,7])#- np.array([5,5,4]))/8+0.5
dip_pos_2 = np.array([7.,5,2])
dip_pos_3 = np.array([8,4,5])
dip_angle_1 = float(350)
dip_angle_2 = float(10)


layers = np.asarray([layer_1,layer_2])#,np.array([[2,5,11],[3,1,2],[6,2,1]])])
test.layers = np.asarray([layer_1,layer_2])
dips = np.asarray([dip_pos_1, dip_pos_2])#, dip_pos_3])
dips_angles = np.asarray([dip_angle_1, dip_angle_2], dtype="float64")
azimuths = np.asarray([90,90], dtype="float64")
polarity = np.asarray([1,1], dtype="float64")
#print (dips_angles)
np.append(np.tile("Layer_1", len(layer_1)), np.tile("Layer_2", len(layer_2)))
np.append(np.tile("Layer_1", len(dip_pos_1)), np.tile("Layer_2", len(dip_pos_2)))


    

    
Out[41]:





array(['Layer_1', 'Layer_1', 'Layer_1', 'Layer_2', 'Layer_2', 'Layer_2'], 
      dtype='<U7')

In [42]:

    

test.Interfaces = pn.DataFrame(data = {"X" :np.append(layer_1[:, 0],layer_2[:,0]),
                                      "Y" :np.append(layer_1[:, 1],layer_2[:,1]),
                                      "Z" :np.append(layer_1[:, 2],layer_2[:,2]),
"formation" : np.append(np.tile("Layer_1", len(layer_1)), np.tile("Layer_2", len(layer_2)))})
test.Foliations =   pn.DataFrame(data = {"X" :np.append(dip_pos_1[0],dip_pos_2[0]),
                                          "Y" :np.append(dip_pos_1[ 1],dip_pos_2[1]),
                                          "Z" :np.append(dip_pos_1[ 2],dip_pos_2[2]),
                                         "azimuth" : azimuths,
                                         "dip" : dips_angles,
                                         "polarity" : polarity,
"formation" : ["Layer_1", "Layer_2"]})
test.Foliations, test.Interfaces
test.formations = test.Interfaces["formation"].unique()
test.set_series()
test.Interfaces


    

    
Out[42]:






  

    

      

      
X
      
Y
      
Z
      
formation
    
  
  

    

      
0
      
0.5
      
4.0
      
7.0
      
Layer_1
    
    

      
1
      
2.0
      
4.0
      
6.0
      
Layer_1
    
    

      
2
      
4.0
      
4.0
      
7.0
      
Layer_1
    
    

      
3
      
5.0
      
4.0
      
5.0
      
Layer_1
    
    

      
4
      
8.0
      
4.0
      
7.0
      
Layer_1
    
    

      
5
      
7.0
      
4.0
      
8.0
      
Layer_1
    
    

      
6
      
1.0
      
5.0
      
7.0
      
Layer_1
    
    

      
7
      
3.0
      
5.0
      
2.0
      
Layer_2
    
    

      
8
      
6.0
      
5.0
      
2.0
      
Layer_2
    
  

In [43]:

    

% matplotlib notebook
test.compute_potential_field( verbose=1)
test.plot_potential_field_2D(direction = "y", cell_pos = 12, figsize=(7,6), contour_lines = 20)


    

    




The serie formations are Layer_1|Layer_2






    















    

In [33]:

    

%%timeit
# Reset the block
test.block.set_value(np.zeros_like(test.grid[:,0]))

# Compute the block
test.compute_block_model([0], verbose = 0)
test.grad


    

    




1 loop, best of 3: 815 ms per loop

In [10]:

    

%matplotlib qt4
sandstone = test
plot_block =  sandstone.block.get_value().reshape(40,40,60)
plt.imshow(plot_block[:,13,:].T, origin = "bottom", aspect = "equal",
           extent = (sandstone.xmin, sandstone.xmax, sandstone.zmin, sandstone.zmax), interpolation = "none")


    

    




Warning: Cannot change to a different GUI toolkit: qt4. Using notebook instead.






    















    












    
Out[10]:





<matplotlib.image.AxesImage at 0x7ff5bc2d3dd8>

In [3]:

    

rest = np.vstack((i[1:] for i in layers))
ref = np.vstack((np.tile(i[0],(np.shape(i)[0]-1,1)) for i in layers))
dips_angles.dtype
rest = rest.astype("float64")
ref = ref.astype("float64")

test.azimuths = azimuths
test.dips_angles = dips_angles
test.polarity = polarity
test.dips = dips.astype("float32")
test.azimuths = azimuths
test.dips_angles = dips_angles
test.polarity = polarity
test.dips = dips.astype("float32")
dips_angles = dips_angles.astype("float32")
type(dips_angles)
#test.G_x,test.G_y,test.G_z = test.interpolate(test.dips,dips_angles,
#                                              azimuths,polarity, rest, ref)[-3:]
rest, ref;


    

In [26]:

    

par2 = 1/(7.5**2)
w = 1/7.5
test.potential_field = test.interpolate(dips,dips_angles,azimuths,polarity, rest, ref,
par2, w)[0].reshape(40,10,40)
test.potential_field = np.swapaxes(test.potential_field,0,1)

test.azimuths = azimuths
test.dips_angles = dips_angles
test.polarity = polarity
test.dips = dips.astype("float32")
test.calculate_gradient()
test.plot_potential_field_2D(direction="y", cell_pos = 4, figsize = (6.2,5), colorbar = True)


    

    




---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-26-43b7eb110886> in <module>()
      1 par2 = 1/(7.5**2)
      2 w = 1/7.5
----> 3 test.potential_field = test.interpolate(dips,dips_angles,azimuths,polarity, rest, ref,
      4 par2, w)[0].reshape(40,10,40)
      5 test.potential_field = np.swapaxes(test.potential_field,0,1)

NameError: name 'rest' is not defined

In [30]:

    

plot_this_crap("y")


    

In [121]:

    

"""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})


    

    




---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-121-8308869b7d0c> in <module>()
     19 dely = 0.2
     20 delz = 0.2
---> 21 from pyevtk.hl import gridToVTK
     22 # Coordinates
     23 x = np.arange(0, extent_x + 0.1*delx, delx, dtype='float64')

ImportError: No module named 'pyevtk'

In [23]:

    

sys.path.append("/home/bl3/anaconda3/lib/python3.5/site-packages/PyEVTK-1.0.0-py3.5.egg_FILES/pyevtk")
nx = 50
ny = 50
nz = 50

xmin = 1
ymin = 1
zmin = 1
grid =  sol
var_name = "Geology"
#from evtk.hl import gridToVTK
import pyevtk
from pyevtk.hl import gridToVTK

# define coordinates
x = np.zeros(nx + 1)
y = np.zeros(ny + 1)
z = np.zeros(nz + 1)
x[1:] = np.cumsum(delx)
y[1:] = np.cumsum(dely)
z[1:] = np.cumsum(delz)



# plot in coordinates
x += xmin
y += ymin
z += zmin

print (len(x), x)
gridToVTK("GeoMigueller", x, y, z,
          cellData = {var_name: grid})


    

    




51 [ 1.   1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2
  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2]






    
Out[23]:





'/home/bl3/PycharmProjects/GeMpy/GeoMigueller.vtr'

In [6]:

    

%%timeit
test.interpolate(dips,dips_angles,azimuths,polarity, rest, ref,
par2, w);


    

    




100 loops, best of 3: 6.41 ms per loop

In [7]:

    

test.interpolate.profile.summary()


    

    




Function profiling
==================
  Message: /home/miguel/PycharmProjects/GeMpy/GeMpy/GeoMig.py:477
  Time in 412 calls to Function.__call__: 2.742617e+00s
  Time in Function.fn.__call__: 2.696156e+00s (98.306%)
  Time in thunks: 2.541325e+00s (92.661%)
  Total compile time: 1.843985e+01s
    Number of Apply nodes: 254
    Theano Optimizer time: 2.322626e+00s
       Theano validate time: 1.088021e-01s
    Theano Linker time (includes C, CUDA code generation/compiling): 1.604251e+01s
       Import time 1.168010e-01s

Time in all call to theano.grad() 0.000000e+00s
Time since theano import 61.358s
Class
---
<% time> <sum %> <apply time> <time per call> <type> <#call> <#apply> <Class name>
  68.9%    68.9%       1.751s       6.25e-05s     C    28016      68   theano.tensor.elemwise.Elemwise
  11.1%    80.0%       0.282s       6.85e-05s     C     4120      10   theano.tensor.blas.Dot22Scalar
   9.7%    89.7%       0.246s       5.43e-05s     C     4532      11   theano.tensor.basic.Alloc
   4.9%    94.6%       0.125s       4.34e-05s     C     2884       7   theano.tensor.elemwise.Sum
   1.6%    96.2%       0.040s       9.76e-05s     Py     412       1   theano.tensor.nlinalg.MatrixInverse
   1.0%    97.2%       0.026s       6.34e-05s     C      412       1   theano.tensor.blas_c.CGemv
   0.7%    98.0%       0.018s       7.33e-06s     C     2472       6   theano.tensor.basic.Join
   0.6%    98.6%       0.016s       3.83e-05s     Py     412       1   theano.tensor.extra_ops.FillDiagonal
   0.5%    99.1%       0.012s       1.18e-06s     C    10300      25   theano.tensor.subtensor.IncSubtensor
   0.3%    99.4%       0.009s       4.13e-07s     C    21424      52   theano.tensor.elemwise.DimShuffle
   0.2%    99.6%       0.005s       4.87e-07s     C    11124      27   theano.tensor.basic.Reshape
   0.2%    99.8%       0.005s       7.00e-07s     C     7004      17   theano.tensor.subtensor.Subtensor
   0.1%    99.9%       0.002s       4.21e-07s     C     5356      13   theano.tensor.opt.MakeVector
   0.0%    99.9%       0.001s       4.01e-07s     C     2884       7   theano.tensor.basic.ScalarFromTensor
   0.0%   100.0%       0.001s       4.06e-07s     C     2472       6   theano.compile.ops.Shape_i
   0.0%   100.0%       0.000s       5.78e-07s     C      412       1   theano.tensor.basic.AllocEmpty
   0.0%   100.0%       0.000s       5.57e-07s     C      412       1   theano.compile.ops.Rebroadcast
   ... (remaining 0 Classes account for   0.00%(0.00s) of the runtime)

Ops
---
<% time> <sum %> <apply time> <time per call> <type> <#call> <#apply> <Op name>
  30.4%    30.4%       0.771s       1.87e-03s     C      412        1   Elemwise{Composite{(i0 * i1 * LT(Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i2, i3, i4), i5) * (((i6 + ((i7 * Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i2, i3, i4)) / i8)) - ((i9 * Composite{(sqr(i0) * i0)}(Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i2, i3, i4))) / i10)) + ((i11 * Composite{(sqr(sqr(i0)) * i0)}(Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i2, i3, i4))) / i12)))}}[(0, 4)]
  29.3%    59.6%       0.745s       1.81e-03s     C      412        1   Elemwise{Composite{(i0 * ((LT(Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i1, i2, i3), i4) * ((i5 + (i6 * Composite{(sqr(i0) * i0)}((Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i1, i2, i3) / i4))) + (i7 * Composite{((sqr(sqr(i0)) * sqr(i0)) * i0)}((Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i1, i2, i3) / i4)))) - ((i8 * sqr((Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i1, i2, i3) / i4))) + (i9 * Composite{(sqr(sqr(i0)
  11.1%    70.8%       0.282s       6.85e-05s     C     4120       10   Dot22Scalar
   9.7%    80.4%       0.246s       5.43e-05s     C     4532       11   Alloc
   4.1%    84.5%       0.104s       1.27e-04s     C      824        2   Elemwise{Mul}[(0, 1)]
   2.3%    86.8%       0.058s       7.03e-05s     C      824        2   Sum{axis=[0], acc_dtype=float64}
   2.0%    88.8%       0.051s       1.24e-04s     C      412        1   Sum{axis=[0], acc_dtype=float64}
   1.8%    90.7%       0.047s       2.84e-05s     C     1648        4   Elemwise{Cast{float64}}
   1.6%    92.3%       0.040s       9.76e-05s     Py     412        1   MatrixInverse
   1.4%    93.7%       0.035s       7.17e-06s     C     4944       12   Elemwise{sub,no_inplace}
   1.0%    94.7%       0.026s       6.34e-05s     C      412        1   CGemv{inplace}
   0.7%    95.4%       0.018s       7.33e-06s     C     2472        6   Join
   0.6%    96.0%       0.016s       9.86e-06s     C     1648        4   Sum{axis=[1], acc_dtype=float64}
   0.6%    96.7%       0.016s       3.83e-05s     Py     412        1   FillDiagonal
   0.6%    97.3%       0.015s       9.33e-06s     C     1648        4   Elemwise{Sqr}[(0, 0)]
   0.4%    97.6%       0.010s       2.37e-05s     C      412        1   Elemwise{Composite{(i0 + (i1 * i2 * i3) + (i4 * i5 * i6 * i7))}}[(0, 3)]
   0.2%    97.8%       0.005s       4.64e-07s     C     10300       25   Reshape{2}
   0.2%    98.0%       0.004s       1.49e-06s     C     2884        7   IncSubtensor{InplaceSet;int64:int64:, int64:int64:}
   0.1%    98.1%       0.003s       5.95e-07s     C     4944       12   Subtensor{::, int64}
   0.1%    98.2%       0.003s       6.53e-06s     C      412        1   Elemwise{Composite{Switch(EQ(Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i0, i1, i2), i3), i3, ((((i4 * i5) / sqr(Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i0, i1, i2))) * ((i6 * LT(Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i0, i1, i2), i7) * i8 * Composite{(((i0 + ((i1 * i2) / i3)) - ((i4 * i5) / i6)) + ((i7 * i8) / i9))}(i9, i10, Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i0, i1, i2), i11, i12, Composite{(sqr(i0
   ... (remaining 53 Ops account for   1.78%(0.05s) of the runtime)

Apply
------
<% time> <sum %> <apply time> <time per call> <#call> <id> <Apply name>
  30.4%    30.4%       0.771s       1.87e-03s    412   248   Elemwise{Composite{(i0 * i1 * LT(Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i2, i3, i4), i5) * (((i6 + ((i7 * Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i2, i3, i4)) / i8)) - ((i9 * Composite{(sqr(i0) * i0)}(Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i2, i3, i4))) / i10)) + ((i11 * Composite{(sqr(sqr(i0)) * i0)}(Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i2, i3, i4))) / i12)))}}[(0, 4)](Subtensor{:int64:}.0, Join.0, Reshape{2
  29.3%    59.6%       0.745s       1.81e-03s    412   247   Elemwise{Composite{(i0 * ((LT(Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i1, i2, i3), i4) * ((i5 + (i6 * Composite{(sqr(i0) * i0)}((Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i1, i2, i3) / i4))) + (i7 * Composite{((sqr(sqr(i0)) * sqr(i0)) * i0)}((Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i1, i2, i3) / i4)))) - ((i8 * sqr((Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i1, i2, i3) / i4))) + (i9 * Composite{(sqr(sqr(i0)) * i0)}((C
   6.3%    65.9%       0.160s       3.87e-04s    412   240   Alloc(CGemv{inplace}.0, Shape_i{0}.0, TensorConstant{1}, TensorConstant{1}, Elemwise{Add}[(0, 1)].0)
   4.9%    70.8%       0.125s       3.03e-04s    412   166   Dot22Scalar(Elemwise{Cast{float64}}.0, InplaceDimShuffle{1,0}.0, TensorConstant{2.0})
   4.1%    74.9%       0.104s       2.53e-04s    412   249   Elemwise{Mul}[(0, 1)](Subtensor{:int64:}.0, InplaceDimShuffle{1,0}.0, Subtensor{int64::}.0, InplaceDimShuffle{x,x}.0)
   3.7%    78.7%       0.095s       2.30e-04s    412   112   Dot22Scalar(Elemwise{Cast{float64}}.0, InplaceDimShuffle{1,0}.0, TensorConstant{2.0})
   3.2%    81.8%       0.081s       1.96e-04s    412   195   Alloc(Subtensor{:int64:}.0, Shape_i{0}.0, TensorConstant{1}, TensorConstant{1}, Elemwise{Composite{Switch(LT(i0, i1), Switch(LT((i2 + i0), i1), i1, (i2 + i0)), Switch(LT(i0, i2), i0, i2))}}.0)
   2.3%    84.1%       0.058s       1.42e-04s    412   111   Dot22Scalar(Elemwise{Cast{float64}}.0, InplaceDimShuffle{1,0}.0, TensorConstant{2.0})
   2.0%    86.2%       0.051s       1.24e-04s    412   252   Sum{axis=[0], acc_dtype=float64}(Elemwise{Mul}[(0, 1)].0)
   1.8%    88.0%       0.046s       1.11e-04s    412    26   Elemwise{Cast{float64}}(Positions of the points to interpolate)
   1.6%    89.5%       0.040s       9.76e-05s    412   238   MatrixInverse(IncSubtensor{InplaceSet;int64::, int64:int64:}.0)
   1.2%    90.7%       0.030s       7.26e-05s    412   251   Sum{axis=[0], acc_dtype=float64}(Elemwise{Composite{(i0 * i1 * LT(Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i2, i3, i4), i5) * (((i6 + ((i7 * Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i2, i3, i4)) / i8)) - ((i9 * Composite{(sqr(i0) * i0)}(Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i2, i3, i4))) / i10)) + ((i11 * Composite{(sqr(sqr(i0)) * i0)}(Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i2, i3, i4))) / i12)))}}[(0, 4)].0)
   1.1%    91.8%       0.028s       6.79e-05s    412   250   Sum{axis=[0], acc_dtype=float64}(Elemwise{Composite{(i0 * ((LT(Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i1, i2, i3), i4) * ((i5 + (i6 * Composite{(sqr(i0) * i0)}((Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i1, i2, i3) / i4))) + (i7 * Composite{((sqr(sqr(i0)) * sqr(i0)) * i0)}((Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i1, i2, i3) / i4)))) - ((i8 * sqr((Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i1, i2, i3) / i4))) + (i9 * 
   1.0%    92.8%       0.026s       6.34e-05s    412   239   CGemv{inplace}(AllocEmpty{dtype='float32'}.0, TensorConstant{1.0}, MatrixInverse.0, IncSubtensor{InplaceSet;int64:int64:}.0, TensorConstant{0.0})
   0.6%    93.5%       0.016s       3.83e-05s    412   220   FillDiagonal(Elemwise{Composite{Switch(EQ(Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i0, i1, i2), i3), i3, ((((i4 * i5) / sqr(Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i0, i1, i2))) * ((i6 * LT(Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i0, i1, i2), i7) * i8 * Composite{(((i0 + ((i1 * i2) / i3)) - ((i4 * i5) / i6)) + ((i7 * i8) / i9))}(i9, i10, Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i0, i1, i2), i11, i12, Composite{(sqr(
   0.6%    94.1%       0.015s       3.71e-05s    412   191   Sum{axis=[1], acc_dtype=float64}(Elemwise{Sqr}[(0, 0)].0)
   0.6%    94.6%       0.015s       3.60e-05s    412   175   Elemwise{Sqr}[(0, 0)](Elemwise{Cast{float64}}.0)
   0.5%    95.1%       0.012s       2.89e-05s    412   173   Join(TensorConstant{0}, Elemwise{sub,no_inplace}.0, Elemwise{sub,no_inplace}.0, Elemwise{sub,no_inplace}.0)
   0.4%    95.6%       0.011s       2.74e-05s    412   158   Elemwise{sub,no_inplace}(InplaceDimShuffle{0,x}.0, InplaceDimShuffle{1,0}.0)
   0.4%    96.0%       0.011s       2.70e-05s    412   159   Elemwise{sub,no_inplace}(InplaceDimShuffle{0,x}.0, InplaceDimShuffle{1,0}.0)
   ... (remaining 234 Apply instances account for 4.00%(0.10s) of the runtime)

Here are tips to potentially make your code run faster
                 (if you think of new ones, suggest them on the mailing list).
                 Test them first, as they are not guaranteed to always provide a speedup.
We don't know if amdlibm will accelerate this scalar op. deg2rad
We don't know if amdlibm will accelerate this scalar op. deg2rad
  - Try installing amdlibm and set the Theano flag lib.amdlibm=True. This speeds up only some Elemwise operation.

In [54]:

    

%%timeit
test.interpolate(dips,dips_angles,azimuths,polarity, rest, ref,
par2, w);


    

    




100 loops, best of 3: 8.83 ms per loop

In [55]:

    

test.interpolate.profile.summary()


    

    




Function profiling
==================
  Message: /home/miguel/PycharmProjects/GeMpy/GeMpy/GeoMig.py:477
  Time in 412 calls to Function.__call__: 3.664509e+00s
  Time in Function.fn.__call__: 3.621449e+00s (98.825%)
  Time in thunks: 3.093245e+00s (84.411%)
  Total compile time: 2.934427e+00s
    Number of Apply nodes: 291
    Theano Optimizer time: 2.653788e+00s
       Theano validate time: 2.586505e-01s
    Theano Linker time (includes C, CUDA code generation/compiling): 1.999319e-01s
       Import time 0.000000e+00s

Time in all call to theano.grad() 0.000000e+00s
Time since theano import 1682.560s
Class
---
<% time> <sum %> <apply time> <time per call> <type> <#call> <#apply> <Class name>
  51.1%    51.1%       1.582s       1.07e-04s     C    14832      36   theano.tensor.elemwise.Elemwise
   9.2%    60.3%       0.284s       6.89e-05s     C     4120      10   theano.sandbox.cuda.basic_ops.GpuAlloc
   8.9%    69.2%       0.275s       3.71e-05s     C     7416      18   theano.sandbox.cuda.basic_ops.HostFromGpu
   6.0%    75.3%       0.187s       4.54e-05s     C     4120      10   theano.tensor.blas.Dot22Scalar
   5.5%    80.8%       0.170s       1.88e-05s     C     9064      22   theano.sandbox.cuda.basic_ops.GpuIncSubtensor
   4.9%    85.6%       0.150s       9.87e-06s     C    15244      37   theano.sandbox.cuda.basic_ops.GpuElemwise
   3.5%    89.1%       0.107s       1.85e-05s     C     5768      14   theano.sandbox.cuda.basic_ops.GpuFromHost
   3.0%    92.0%       0.091s       3.69e-05s     C     2472       6   theano.sandbox.cuda.basic_ops.GpuJoin
   2.5%    94.6%       0.079s       3.18e-05s     C     2472       6   theano.tensor.elemwise.Sum
   2.4%    97.0%       0.074s       8.98e-06s     C     8240      20   theano.sandbox.cuda.basic_ops.GpuReshape
   1.3%    98.2%       0.040s       9.66e-05s     Py     412       1   theano.tensor.nlinalg.MatrixInverse
   0.3%    98.6%       0.010s       6.03e-07s     C    16480      40   theano.sandbox.cuda.basic_ops.GpuDimShuffle
   0.2%    98.8%       0.007s       9.57e-07s     C     7004      17   theano.sandbox.cuda.basic_ops.GpuSubtensor
   0.2%    99.0%       0.007s       1.62e-05s     Py     412       1   theano.tensor.extra_ops.FillDiagonal
   0.2%    99.2%       0.006s       1.51e-05s     C      412       1   theano.sandbox.cuda.basic_ops.GpuCAReduce
   0.2%    99.4%       0.005s       1.33e-05s     Py     412       1   theano.compile.ops.Rebroadcast
   0.1%    99.5%       0.003s       8.11e-06s     C      412       1   theano.sandbox.cuda.basic_ops.GpuAllocEmpty
   0.1%    99.6%       0.003s       5.75e-07s     C     5356      13   theano.tensor.opt.MakeVector
   0.1%    99.7%       0.003s       5.76e-07s     C     4944      12   theano.tensor.elemwise.DimShuffle
   0.1%    99.8%       0.003s       6.50e-06s     C      412       1   theano.sandbox.cuda.blas.GpuGemv
   ... (remaining 5 Classes account for   0.24%(0.01s) of the runtime)

Ops
---
<% time> <sum %> <apply time> <time per call> <type> <#call> <#apply> <Op name>
  25.4%    25.4%       0.785s       1.90e-03s     C      412        1   Elemwise{Composite{(i0 * i1 * LT(Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i2, i3, i4), i5) * (((i6 + ((i7 * Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i2, i3, i4)) / i8)) - ((i9 * Composite{(sqr(i0) * i0)}(Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i2, i3, i4))) / i10)) + ((i11 * Composite{(sqr(sqr(i0)) * i0)}(Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i2, i3, i4))) / i12)))}}[(0, 4)]
  24.0%    49.3%       0.741s       1.80e-03s     C      412        1   Elemwise{Composite{(i0 * ((LT(Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i1, i2, i3), i4) * ((i5 + (i6 * Composite{(sqr(i0) * i0)}((Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i1, i2, i3) / i4))) + (i7 * Composite{((sqr(sqr(i0)) * sqr(i0)) * i0)}((Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i1, i2, i3) / i4)))) - ((i8 * sqr((Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i1, i2, i3) / i4))) + (i9 * Composite{(sqr(sqr(i0)
   8.9%    58.2%       0.275s       3.71e-05s     C     7416       18   HostFromGpu
   7.0%    65.2%       0.215s       7.46e-05s     C     2884        7   GpuAlloc
   6.0%    71.2%       0.187s       4.54e-05s     C     4120       10   Dot22Scalar
   3.5%    74.7%       0.107s       1.85e-05s     C     5768       14   GpuFromHost
   3.0%    77.6%       0.091s       3.69e-05s     C     2472        6   GpuJoin
   2.4%    80.0%       0.073s       9.81e-06s     C     7416       18   GpuReshape{2}
   2.2%    82.2%       0.069s       5.56e-05s     C     1236        3   GpuAlloc{memset_0=True}
   2.0%    84.2%       0.062s       7.50e-05s     C      824        2   Sum{axis=[0], acc_dtype=float64}
   1.9%    86.1%       0.059s       3.58e-05s     C     1648        4   GpuIncSubtensor{InplaceSet;int64:int64:, int64:int64:}
   1.9%    88.0%       0.058s       1.18e-05s     C     4944       12   GpuElemwise{sub,no_inplace}
   1.3%    89.3%       0.040s       9.66e-05s     Py     412        1   MatrixInverse
   1.0%    90.2%       0.030s       3.61e-05s     C      824        2   GpuIncSubtensor{InplaceSet;int64:int64:, int64::}
   0.9%    91.1%       0.027s       8.25e-06s     C     3296        8   GpuIncSubtensor{InplaceSet;int64:int64:, int64}
   0.9%    92.0%       0.027s       3.30e-05s     C      824        2   GpuIncSubtensor{InplaceSet;int64::, int64:int64:}
   0.6%    92.6%       0.020s       9.60e-06s     C     2060        5   GpuElemwise{mul,no_inplace}
   0.5%    93.2%       0.017s       1.02e-05s     C     1648        4   Sum{axis=[1], acc_dtype=float64}
   0.5%    93.7%       0.015s       9.10e-06s     C     1648        4   Elemwise{Sqr}[(0, 0)]
   0.5%    94.1%       0.015s       8.95e-06s     C     1648        4   GpuIncSubtensor{InplaceSet;:int64:, int64}
   ... (remaining 71 Ops account for   5.86%(0.18s) of the runtime)

Apply
------
<% time> <sum %> <apply time> <time per call> <#call> <id> <Apply name>
  25.4%    25.4%       0.785s       1.90e-03s    412   285   Elemwise{Composite{(i0 * i1 * LT(Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i2, i3, i4), i5) * (((i6 + ((i7 * Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i2, i3, i4)) / i8)) - ((i9 * Composite{(sqr(i0) * i0)}(Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i2, i3, i4))) / i10)) + ((i11 * Composite{(sqr(sqr(i0)) * i0)}(Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i2, i3, i4))) / i12)))}}[(0, 4)](HostFromGpu.0, HostFromGpu.0, Reshape{2
  24.0%    49.3%       0.741s       1.80e-03s    412   284   Elemwise{Composite{(i0 * ((LT(Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i1, i2, i3), i4) * ((i5 + (i6 * Composite{(sqr(i0) * i0)}((Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i1, i2, i3) / i4))) + (i7 * Composite{((sqr(sqr(i0)) * sqr(i0)) * i0)}((Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i1, i2, i3) / i4)))) - ((i8 * sqr((Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i1, i2, i3) / i4))) + (i9 * Composite{(sqr(sqr(i0)) * i0)}((C
   2.5%    51.8%       0.077s       1.87e-04s    412   190   Dot22Scalar(Elemwise{Cast{float64}}.0, InplaceDimShuffle{1,0}.0, TensorConstant{2.0})
   2.3%    54.1%       0.070s       1.70e-04s    412   274   GpuAlloc(GpuGemv{inplace}.0, Shape_i{0}.0, TensorConstant{1}, TensorConstant{1}, Elemwise{Add}[(0, 1)].0)
   2.2%    56.3%       0.068s       1.64e-04s    412   217   GpuAlloc(GpuSubtensor{:int64:}.0, Shape_i{0}.0, TensorConstant{1}, TensorConstant{1}, Elemwise{Composite{Switch(LT(i0, i1), Switch(LT((i2 + i0), i1), i1, (i2 + i0)), Switch(LT(i0, i2), i0, i2))}}.0)
   2.0%    58.3%       0.062s       1.51e-04s    412   146   Dot22Scalar(Elemwise{Cast{float64}}.0, InplaceDimShuffle{1,0}.0, TensorConstant{2.0})
   2.0%    60.3%       0.062s       1.50e-04s    412   281   HostFromGpu(GpuSubtensor{int64:int64:}.0)
   1.9%    62.2%       0.060s       1.46e-04s    412   282   HostFromGpu(GpuSubtensor{:int64:}.0)
   1.7%    63.9%       0.052s       1.26e-04s    412   153   GpuAlloc{memset_0=True}(CudaNdarrayConstant{0.0}, Elemwise{Add}[(0, 1)].0, Elemwise{Add}[(0, 1)].0)
   1.4%    65.3%       0.043s       1.05e-04s    412   271   GpuFromHost(MatrixInverse.0)
   1.4%    66.7%       0.042s       1.02e-04s    412   145   Dot22Scalar(Elemwise{Cast{float64}}.0, InplaceDimShuffle{1,0}.0, TensorConstant{2.0})
   1.3%    68.0%       0.040s       9.66e-05s    412   270   MatrixInverse(HostFromGpu.0)
   1.2%    69.2%       0.037s       8.90e-05s    412   187   GpuJoin(TensorConstant{0}, GpuElemwise{Composite{(i0 * (i1 - i2))},no_inplace}.0, GpuElemwise{Composite{(i0 * (i1 - i2))},no_inplace}.0, GpuElemwise{Composite{(i0 * (i1 - i2))}}[(0, 1)].0, GpuElemwise{Composite{(i0 * (sqr(i1) - sqr(i2)))},no_inplace}.0, GpuElemwise{Composite{(i0 * (sqr(i1) - sqr(i2)))},no_inplace}.0, GpuElemwise{Composite{(i0 * (sqr(i1) - sqr(i2)))},no_inplace}.0, GpuElemwise{Composite{(i0 * ((i1 * i2) - (i3 * i4)))},no_inplace}.0, 
   1.1%    70.3%       0.035s       8.39e-05s    412   197   HostFromGpu(GpuJoin.0)
   1.1%    71.3%       0.033s       7.96e-05s    412   287   Sum{axis=[0], acc_dtype=float64}(Elemwise{Composite{(i0 * ((LT(Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i1, i2, i3), i4) * ((i5 + (i6 * Composite{(sqr(i0) * i0)}((Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i1, i2, i3) / i4))) + (i7 * Composite{((sqr(sqr(i0)) * sqr(i0)) * i0)}((Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i1, i2, i3) / i4)))) - ((i8 * sqr((Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i1, i2, i3) / i4))) + (i9 * 
   0.9%    72.3%       0.029s       7.03e-05s    412   288   Sum{axis=[0], acc_dtype=float64}(Elemwise{Composite{(i0 * i1 * LT(Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i2, i3, i4), i5) * (((i6 + ((i7 * Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i2, i3, i4)) / i8)) - ((i9 * Composite{(sqr(i0) * i0)}(Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i2, i3, i4))) / i10)) + ((i11 * Composite{(sqr(sqr(i0)) * i0)}(Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}(i2, i3, i4))) / i12)))}}[(0, 4)].0)
   0.7%    73.0%       0.022s       5.26e-05s    412   206   HostFromGpu(GpuDimShuffle{1,0}.0)
   0.7%    73.7%       0.022s       5.25e-05s    412   207   HostFromGpu(GpuDimShuffle{1,0}.0)
   0.7%    74.3%       0.021s       5.07e-05s    412    20   HostFromGpu(Positions of the points to interpolate)
   0.6%    74.9%       0.018s       4.43e-05s    412    34   GpuAlloc(GpuFromHost.0, TensorConstant{3}, TensorConstant{1}, Shape_i{0}.0, Shape_i{1}.0)
   ... (remaining 271 Apply instances account for 25.07%(0.78s) of the runtime)

Here are tips to potentially make your code run faster
                 (if you think of new ones, suggest them on the mailing list).
                 Test them first, as they are not guaranteed to always provide a speedup.
  - Try installing amdlibm and set the Theano flag lib.amdlibm=True. This speeds up only some Elemwise operation.