In [1]:

    

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

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


    

In [2]:

    

test = GeoMig.GeoMigSim_pro2(c_o = np.float32(-0.01),range = 17)

test.create_regular_grid_3D(0,10,0,10,0,10,200,200,200)
test.theano_set_3D()


    

    




/home/miguel/anaconda3/lib/python3.5/site-packages/ipykernel/__main__.py:4: DeprecationWarning: stack(*tensors) interface is deprecated, use stack(tensors, axis=0) instead.
WARNING (theano.tensor.blas): We did not found a dynamic library into the library_dir of the library we use for blas. If you use ATLAS, make sure to compile it with dynamics library.

In [81]:

    

layer_1 = np.array([[1,5,7], [9,5,6.8]], dtype = "float32")

layer_2 = np.array([[5,5,1],[7,5,1]], dtype = "float32")

dip_pos_1 = np.array([2,5,6], dtype = "float32")
dip_pos_2 = np.array([6.,4,6], dtype = "float32")
dip_pos_3 = np.array([8,4,5], dtype = "float32")
dip_angle_1 = float(-45)
dip_angle_2 = float(45)


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="float32")
azimuths = np.asarray([2, 20], dtype="float32")
polarity = np.asarray([1, 1], dtype="float32")
#print (dips_angles)
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("float32")
ref = ref.astype("float32")
dips = dips.astype("float32")
dips_angles = dips_angles.astype("float32")
type(dips_angles)


    

    
Out[81]:





numpy.ndarray

In [82]:

    

G_x = np.sin(np.deg2rad(dips_angles)) * np.sin(np.deg2rad(azimuths)) * polarity
G_y = np.sin(np.deg2rad(dips_angles)) * np.cos(np.deg2rad(azimuths)) * polarity
G_z = np.cos(np.deg2rad(dips_angles)) * polarity

G_x, G_y, G_z


    

    
Out[82]:





(array([-0.02,  0.24], dtype=float32),
 array([-0.71,  0.66], dtype=float32),
 array([ 0.71,  0.71], dtype=float32))

In [83]:

    

test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[0]


    

    
Out[83]:





array([ 0.  ,  0.2 ,  0.39, ..., -1.27, -1.34, -1.41], dtype=float32)

In [84]:

    

_,h1 = np.argmin((abs(test.grid - ref[0])).sum(1)), test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[0][np.argmin((abs(test.grid - ref[0])).sum(1))]
_, h2 =np.argmin((abs(test.grid - ref[1])).sum(1)), test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[0][np.argmin((abs(test.grid - ref[1])).sum(1))]


    

In [85]:

    

test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[0][np.argmin((abs(test.grid - ref[0])).sum(1))]


    

    
Out[85]:





41.538002

In [86]:

    

test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[0][np.argmin((abs(test.grid - rest[0])).sum(1))]


    

    
Out[86]:





41.581768

In [87]:

    

sol = test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[0].reshape(200,200,200, order = "C")[:,:,::-1].transpose()
#sol = np.swapaxes(sol,0,1)


    

In [88]:

    

G_x, G_y, G_z = test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[-3:]
G_x, G_y, G_z


    

    
Out[88]:





(array([-0.02,  0.24], dtype=float32),
 array([-0.71,  0.66], dtype=float32),
 array([ 0.71,  0.71], dtype=float32))

In [89]:

    

import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import matplotlib
import matplotlib.cm as cmx

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
h =  np.array([h1,h2])
cm = plt.get_cmap("jet")
cNorm = matplotlib.colors.Normalize(vmin=h.min(), vmax=h.max())
scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=cm)


sol = test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[0].reshape(200,200,200,
    order = "C")[:,:,:]
#sol = np.swapaxes(sol,0,1)


from skimage import measure
isolines = np.linspace(h1,h2,2)
#vertices = measure.marching_cubes(sol, isolines[0], spacing = (0.2,0.2,0.2),
#    gradient_direction = "descent")[0]
for i in isolines[0:10]:
    vertices = measure.marching_cubes(sol, i, spacing = (0.05,0.05,0.05),
    gradient_direction = "ascent")[0]

    ax.scatter(vertices[::40,0],vertices[::40,1],vertices[::40,2],color=scalarMap.to_rgba(i),
     alpha = 0.2) #color=scalarMap.to_rgba(vertices[::10,2])
ax.scatter(layers[0][:,0],layers[0][:,1],layers[0][:,2], s = 50, c = "r" )
ax.scatter(layers[1][:,0],layers[1][:,1],layers[1][:,2], s = 50, c = "g" )
ax.quiver3D(dips[:,0],dips[:,1],dips[:,2], G_x,G_y,G_z, pivot = "tail", linewidths = 2)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")
ax.set_xlim(0,10)
ax.set_ylim(0,10)
ax.set_zlim(0,10)
#ax.scatter(simplices[:,0],simplices[:,1],simplices[:,2])


    

    
Out[89]:





(0, 10)

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [56]:

    




    

    
Out[56]:





array([[ 0.  ,  0.2 ,  8.8 ],
       [ 0.2 ,  0.2 ,  8.8 ],
       [ 0.2 ,  0.16,  8.8 ],
       ..., 
       [ 9.8 ,  9.4 ,  8.65],
       [ 9.8 ,  9.6 ,  8.64],
       [ 9.8 ,  9.8 ,  8.64]])

In [10]:

    

import matplotlib.pyplot as plt

def plot_this_crap(direction):
    fig = plt.figure()
    ax = fig.add_subplot(111)
    dip_pos_1_v = np.array([G_x[0], G_y[0],
                           G_z[0]]) + dip_pos_1

    dip_pos_2_v = np.array([G_x[1], G_y[1],
                           G_z[1]]) + dip_pos_2
    if direction == "x":

        plt.arrow(dip_pos_1[1],dip_pos_1[2], dip_pos_1_v[1]-dip_pos_1[1],
                  dip_pos_1_v[2]-dip_pos_1[2], head_width = 0.2)
        plt.arrow(dip_pos_2[1],dip_pos_2[2],dip_pos_2_v[1]-dip_pos_2[1], 
                  dip_pos_2_v[2]-dip_pos_2[2], head_width = 0.2)

        plt.plot(layer_1[:,1],layer_1[:,2], "o")
        plt.plot(layer_2[:,1],layer_2[:,2], "o")

        plt.plot(layer_1[:,1],layer_1[:,2], )
        plt.plot(layer_2[:,1],layer_2[:,2], )
        plt.contour( sol[25,:,:] ,30,extent = (0,10,0,10) )

    if direction == "y":

        plt.arrow(dip_pos_1[0],dip_pos_1[2], dip_pos_1_v[0]-dip_pos_1[0],
                  dip_pos_1_v[2]-dip_pos_1[2], head_width = 0.2)
        plt.arrow(dip_pos_2[0],dip_pos_2[2],dip_pos_2_v[0]-dip_pos_2[0], 
                  dip_pos_2_v[2]-dip_pos_2[2], head_width = 0.2)

        plt.plot(layer_1[:,0],layer_1[:,2], "o")
        plt.plot(layer_2[:,0],layer_2[:,2], "o")

        plt.plot(layer_1[:,0],layer_1[:,2], )
        plt.plot(layer_2[:,0],layer_2[:,2], )
        plt.contour( sol[:,25,:].T ,30,extent = (0,10,0,10) )

    if direction == "z":

        plt.arrow(dip_pos_1[0],dip_pos_1[1], dip_pos_1_v[0]-dip_pos_1[0],
                  dip_pos_1_v[1]-dip_pos_1[1], head_width = 0.2)
        plt.arrow(dip_pos_2[0],dip_pos_2[1],dip_pos_2_v[0]-dip_pos_2[0], 
                  dip_pos_2_v[1]-dip_pos_2[1], head_width = 0.2)

        plt.plot(layer_1[:,0],layer_1[:,1], "o")
        plt.plot(layer_2[:,0],layer_2[:,1], "o")

        plt.plot(layer_1[:,0],layer_1[:,1], )
        plt.plot(layer_2[:,0],layer_2[:,1], )
        plt.contour( sol[:,:,25] ,30,extent = (0,10,0,10) )
    
    
    
#plt.colorbar()
#plt.xlim(0,10)
#plt.ylim(0,10)
    plt.colorbar()
    plt.title("GeoBulleter v 0.1")


    

In [13]:

    

plot_this_crap("z")


    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [50]:

    

h,j,k =sol[5,10,35], sol[25,5,5], sol[30,15,-25]
layer_1 = np.array([[1,5,7],[5,5,7],[6,5,7], [9,5,7]], dtype = "float32")

layer_2 = np.array([[1,5,1],[5,5,1],[9,5,1]], dtype = "float32")

print(sol[5,25,35], sol[25,25,35],  sol[30,25,35],  sol[45,25,35]) 
print(sol[5,25,5], sol[25,25,5],  sol[45,25,5])


    

    




6.05973 3.17981 2.4532 0.29113
6.14472 3.30846 0.500162

In [51]:

    

list(layer_1[0]*5)


    

    
Out[51]:





[5.0, 25.0, 35.0]

In [52]:

    

interfaces_aux = test.geoMigueller(dips,dips_angles,azimuths,polarity,
                rest, ref)[0]


h = sol[10,20,30]# interfaces_aux[np.argmin(abs((test.grid - ref[0]).sum(1)))]
k = sol[30,15,25]# interfaces_aux[np.argmin(abs((test.grid - dips[0]).sum(1)))]
j = sol[45,25,5]#interfaces_aux[np.argmin(abs((test.grid - dips[-1]).sum(1)))]
h,k,j


    

    
Out[52]:





(5.2851176, 2.330806, 0.50016159)

In [53]:

    

dips[-1], ref[0]


    

    
Out[53]:





(array([ 6.,  3.,  5.], dtype=float32), array([ 1.,  5.,  7.], dtype=float32))

In [54]:

    

sol[30,15,25], sol[30,15,25]


    

    
Out[54]:





(2.330806, 2.330806)

In [59]:

    

sol = test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[0].reshape(50,50,50, order = "C")
sol = np.swapaxes(sol,0,1)
plt.contour(sol[:,25,:].transpose())


    

    
Out[59]:





<matplotlib.contour.QuadContourSet at 0x7fa0222424e0>






    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

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 [201]:

    

len(x)


    

    
Out[201]:





51

In [202]:

    

surf_eq.min()


    

    
Out[202]:





-63.0

In [203]:

    

np.min(z)


    

    
Out[203]:





0.0

In [204]:

    

layers[0][:,0]


    

    
Out[204]:





array([ 1.,  5.,  6.,  9.], dtype=float32)

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [168]:

    

G_x = np.sin(np.deg2rad(dips_angles)) * np.sin(np.deg2rad(azimuths)) * polarity
G_y = np.sin(np.deg2rad(dips_angles)) * np.cos(np.deg2rad(azimuths)) * polarity
G_z = np.cos(np.deg2rad(dips_angles)) * polarity


    

In [183]:

    

a


    

    
Out[183]:





array([[ 2.  ,  6.  ],
       [ 2.03,  6.15]], dtype=float32)

In [ ]:

    

data = [trace1, trace2]
layout = go.Layout(
    xaxis=dict(
        range=[2, 5]
    ),
    yaxis=dict(
        range=[2, 5]
    )
)
fig = go.Figure(data=data, layout=layout)


    

In [ ]:

    




    

In [53]:

    




    

In [ ]:

    

import lxml
lxml??


    

In [11]:

    

# Random Box

#layers = [np.random.uniform(0,10,(10,2)) for i in range(100)]
#dips = np.random.uniform(0,10, (60,2))
#dips_angles = np.random.normal(90,10, 60)
#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))
#rest;


    

In [ ]:

    




    

In [63]:

    

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
X, Y, Z = axes3d.get_test_data(0.05)
cset = ax.contour(X, Y, Z, cmap=cm.coolwarm)
ax.clabel(cset, fontsize=9, inline=1)
print(X)
plt.show()


    

    




[[-30.  -29.5 -29.  ...,  28.5  29.   29.5]
 [-30.  -29.5 -29.  ...,  28.5  29.   29.5]
 [-30.  -29.5 -29.  ...,  28.5  29.   29.5]
 ..., 
 [-30.  -29.5 -29.  ...,  28.5  29.   29.5]
 [-30.  -29.5 -29.  ...,  28.5  29.   29.5]
 [-30.  -29.5 -29.  ...,  28.5  29.   29.5]]






    

In [11]:

    

import matplotlib.pyplot as plt
% matplotlib inline
plt.contour( sol.reshape(100,100) ,30,extent = (0,10,0,10) )


    

    
Out[11]:





<matplotlib.contour.QuadContourSet at 0x7ff55335e7f0>






    

In [131]:

    

import matplotlib.pyplot as plt
% matplotlib inline
dip_pos_1_v = np.array([np.cos(np.deg2rad(dip_angle_1))*1,
                        np.sin(np.deg2rad(dip_angle_1))]) + dip_pos_1

dip_pos_2_v = np.array([np.cos(np.deg2rad(dip_angle_2))*1, 
                        np.sin(np.deg2rad(dip_angle_2))]) + dip_pos_2

plt.arrow(dip_pos_1[0],dip_pos_1[1], dip_pos_1_v[0]-dip_pos_1[0],
          dip_pos_1_v[1]-dip_pos_1[1], head_width = 0.2)
plt.arrow(dip_pos_2[0],dip_pos_2[1],dip_pos_2_v[0]-dip_pos_2[0], 
          dip_pos_2_v[1]-dip_pos_2[1], head_width = 0.2)

plt.plot(layer_1[:,0],layer_1[:,1], "o")
plt.plot(layer_2[:,0],layer_2[:,1], "o")

plt.plot(layer_1[:,0],layer_1[:,1], )
plt.plot(layer_2[:,0],layer_2[:,1], )

plt.contour( sol.reshape(100,100) ,30,extent = (0,10,0,10) )
#plt.colorbar()
#plt.xlim(0,10)
#plt.ylim(0,10)
plt.title("GeoBulleter v 0.1")
print (dip_pos_1_v, dip_pos_2_v, layer_1)


    

    




[ 2.  5.] [ 6.34  3.94] [[ 1.  7.]
 [ 5.  7.]
 [ 6.  7.]
 [ 9.  8.]]






    

In [17]:

    

%%timeit
sol = test.geoMigueller(dips,dips_angles,rest, ref)[0]


    

    




1 loop, best of 3: 5.81 s per loop

In [18]:

    

test.geoMigueller.profile.summary()


    

    




Function profiling
==================
  Message: /home/miguel/PycharmProjects/GeMpy/GeoMig.py:562
  Time in 5 calls to Function.__call__: 2.937774e+01s
  Time in Function.fn.__call__: 2.937736e+01s (99.999%)
  Time in thunks: 2.934835e+01s (99.900%)
  Total compile time: 1.559712e+00s
    Number of Apply nodes: 171
    Theano Optimizer time: 1.410723e+00s
       Theano validate time: 4.508591e-02s
    Theano Linker time (includes C, CUDA code generation/compiling): 9.198022e-02s
       Import time 0.000000e+00s

Time in all call to theano.grad() 0.000000e+00s
Time since theano import 132.105s
Class
---
<% time> <sum %> <apply time> <time per call> <type> <#call> <#apply> <Class name>
  83.6%    83.6%      24.529s       1.29e-01s     C      190      38   theano.tensor.elemwise.Elemwise
   6.4%    90.0%       1.889s       5.40e-02s     C       35       7   theano.tensor.elemwise.Sum
   5.1%    95.1%       1.496s       2.99e-02s     C       50      10   theano.tensor.blas.Dot22Scalar
   2.5%    97.6%       0.743s       1.86e-02s     C       40       8   theano.tensor.basic.Alloc
   2.3%   100.0%       0.682s       2.27e-02s     C       30       6   theano.tensor.basic.Join
   0.0%   100.0%       0.008s       1.55e-03s     Py       5       1   theano.tensor.nlinalg.MatrixInverse
   0.0%   100.0%       0.000s       2.76e-06s     C       95      19   theano.tensor.basic.Reshape
   0.0%   100.0%       0.000s       2.37e-06s     C       65      13   theano.tensor.subtensor.IncSubtensor
   0.0%   100.0%       0.000s       9.48e-07s     C      155      31   theano.tensor.elemwise.DimShuffle
   0.0%   100.0%       0.000s       2.77e-05s     Py       5       1   theano.tensor.extra_ops.FillDiagonal
   0.0%   100.0%       0.000s       2.34e-06s     C       55      11   theano.tensor.subtensor.Subtensor
   0.0%   100.0%       0.000s       1.37e-06s     C       60      12   theano.tensor.opt.MakeVector
   0.0%   100.0%       0.000s       1.14e-06s     C       30       6   theano.compile.ops.Shape_i
   0.0%   100.0%       0.000s       5.77e-06s     C        5       1   theano.tensor.blas_c.CGemv
   0.0%   100.0%       0.000s       7.71e-07s     C       30       6   theano.tensor.basic.ScalarFromTensor
   0.0%   100.0%       0.000s       1.19e-06s     C        5       1   theano.tensor.basic.AllocEmpty
   ... (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>
  47.3%    47.3%      13.894s       2.78e+00s     C        5        1   Elemwise{Composite{(i0 * i1 * LT(Composite{sqrt(((i0 + i1) - i2))}(i2, i3, i4), i5) * Composite{(sqr(i0) * i0)}((i5 - Composite{sqrt(((i0 + i1) - i2))}(i2, i3, i4))) * ((i6 * sqr(Composite{sqrt(((i0 + i1) - i2))}(i2, i3, i4))) + i7 + (i8 * i5 * Composite{sqrt(((i0 + i1) - i2))}(i2, i3, i4))))}}[(0, 1)]
  28.1%    75.5%       8.253s       1.65e+00s     C        5        1   Elemwise{Composite{(i0 * ((LT(Composite{sqrt(((i0 + i1) - i2))}(i1, i2, i3), i4) * ((i5 + (i6 * Composite{(sqr(i0) * i0)}((Composite{sqrt(((i0 + i1) - i2))}(i1, i2, i3) / i4))) + (i7 * Composite{((sqr(sqr(i0)) * sqr(i0)) * i0)}((Composite{sqrt(((i0 + i1) - i2))}(i1, i2, i3) / i4)))) - ((i8 * sqr((Composite{sqrt(((i0 + i1) - i2))}(i1, i2, i3) / i4))) + (i9 * Composite{(sqr(sqr(i0)) * i0)}((Composite{sqrt(((i0 + i1) - i2))}(i1, i2, i3) / i4
   6.1%    81.6%       1.799s       1.20e-01s     C       15        3   Sum{axis=[0], acc_dtype=float64}
   5.1%    86.7%       1.496s       2.99e-02s     C       50       10   Dot22Scalar
   3.6%    90.3%       1.054s       2.11e-01s     C        5        1   Elemwise{Composite{sqrt(((i0 + i1) - i2))}}[(0, 2)]
   2.8%    93.0%       0.810s       1.62e-02s     C       50       10   Elemwise{sub,no_inplace}
   2.5%    95.6%       0.743s       1.86e-02s     C       40        8   Alloc
   2.3%    97.9%       0.682s       2.27e-02s     C       30        6   Join
   1.5%    99.4%       0.431s       2.87e-02s     C       15        3   Elemwise{mul,no_inplace}
   0.3%    99.7%       0.090s       4.51e-03s     C       20        4   Sum{axis=[1], acc_dtype=float64}
   0.2%    99.9%       0.054s       1.09e-02s     C        5        1   Elemwise{Composite{(i0 + ((i1 * i2) / i3) + i4)}}[(0, 0)]
   0.1%   100.0%       0.034s       1.70e-03s     C       20        4   Elemwise{sqr,no_inplace}
   0.0%   100.0%       0.008s       1.55e-03s     Py       5        1   MatrixInverse
   0.0%   100.0%       0.000s       2.76e-06s     C       95       19   Reshape{2}
   0.0%   100.0%       0.000s       2.77e-05s     Py       5        1   FillDiagonal
   0.0%   100.0%       0.000s       2.18e-05s     C        5        1   Elemwise{Composite{Switch(EQ(Composite{sqrt(((i0 + i1) - i2))}(i0, i1, i2), i3), i3, ((i4 * (((i5 * i6 * i7 * LT(Composite{sqrt(((i0 + i1) - i2))}(i0, i1, i2), i8) * Composite{(sqr((i0 - i1)) * (i0 - i1))}(i8, Composite{sqrt(((i0 + i1) - i2))}(i0, i1, i2)) * Composite{((i0 * i1) + i2 + (i3 * i4 * i5))}(i9, sqr(Composite{sqrt(((i0 + i1) - i2))}(i0, i1, i2)), i10, i11, i8, Composite{sqrt(((i0 + i1) - i2))}(i0, i1, i2))) / (sqr(Composite{sqr
   0.0%   100.0%       0.000s       1.37e-06s     C       60       12   MakeVector{dtype='int64'}
   0.0%   100.0%       0.000s       1.54e-05s     C        5        1   Elemwise{Composite{((((LT(Composite{sqrt(((i0 + i1) - i2))}(i0, i1, i2), i3) * ((i4 + (i5 * Composite{(sqr(i0) * i0)}((Composite{sqrt(((i0 + i1) - i2))}(i0, i1, i2) / i3))) + (i6 * Composite{((sqr(sqr(i0)) * sqr(i0)) * i0)}((Composite{sqrt(((i0 + i1) - i2))}(i0, i1, i2) / i3)))) - ((i7 * sqr((Composite{sqrt(((i0 + i1) - i2))}(i0, i1, i2) / i3))) + (i8 * Composite{(sqr(sqr(i0)) * i0)}((Composite{sqrt(((i0 + i1) - i2))}(i0, i1, i2) / i3))))
   0.0%   100.0%       0.000s       1.73e-06s     C       40        8   Subtensor{::, int64}
   0.0%   100.0%       0.000s       2.33e-06s     C       25        5   IncSubtensor{InplaceSet;int64:int64:, int64:int64:}
   ... (remaining 37 Ops account for   0.00%(0.00s) of the runtime)

Apply
------
<% time> <sum %> <apply time> <time per call> <#call> <id> <Apply name>
  47.3%    47.3%      13.894s       2.78e+00s      5   165   Elemwise{Composite{(i0 * i1 * LT(Composite{sqrt(((i0 + i1) - i2))}(i2, i3, i4), i5) * Composite{(sqr(i0) * i0)}((i5 - Composite{sqrt(((i0 + i1) - i2))}(i2, i3, i4))) * ((i6 * sqr(Composite{sqrt(((i0 + i1) - i2))}(i2, i3, i4))) + i7 + (i8 * i5 * Composite{sqrt(((i0 + i1) - i2))}(i2, i3, i4))))}}[(0, 1)](Subtensor{:int64:}.0, Join.0, Reshape{2}.0, Reshape{2}.0, Dot22Scalar.0, InplaceDimShuffle{x,x}.0, TensorConstant{(1, 1) of 3.0}, Elemwise{mul,no_inp
  28.1%    75.5%       8.253s       1.65e+00s      5   164   Elemwise{Composite{(i0 * ((LT(Composite{sqrt(((i0 + i1) - i2))}(i1, i2, i3), i4) * ((i5 + (i6 * Composite{(sqr(i0) * i0)}((Composite{sqrt(((i0 + i1) - i2))}(i1, i2, i3) / i4))) + (i7 * Composite{((sqr(sqr(i0)) * sqr(i0)) * i0)}((Composite{sqrt(((i0 + i1) - i2))}(i1, i2, i3) / i4)))) - ((i8 * sqr((Composite{sqrt(((i0 + i1) - i2))}(i1, i2, i3) / i4))) + (i9 * Composite{(sqr(sqr(i0)) * i0)}((Composite{sqrt(((i0 + i1) - i2))}(i1, i2, i3) / i4)))))) - (L
   4.6%    80.0%       1.346s       2.69e-01s      5   168   Sum{axis=[0], acc_dtype=float64}(Elemwise{Composite{(i0 * i1 * LT(Composite{sqrt(((i0 + i1) - i2))}(i2, i3, i4), i5) * Composite{(sqr(i0) * i0)}((i5 - Composite{sqrt(((i0 + i1) - i2))}(i2, i3, i4))) * ((i6 * sqr(Composite{sqrt(((i0 + i1) - i2))}(i2, i3, i4))) + i7 + (i8 * i5 * Composite{sqrt(((i0 + i1) - i2))}(i2, i3, i4))))}}[(0, 1)].0)
   3.6%    83.6%       1.054s       2.11e-01s      5    98   Elemwise{Composite{sqrt(((i0 + i1) - i2))}}[(0, 2)](Reshape{2}.0, Reshape{2}.0, Dot22Scalar.0)
   2.7%    86.3%       0.780s       1.56e-01s      5   105   Dot22Scalar(Reshape{2}.0, Positions of the points to interpolate.T, TensorConstant{2.0})
   2.5%    88.8%       0.742s       1.48e-01s      5   157   Alloc(CGemv{inplace}.0, Shape_i{0}.0, TensorConstant{1}, TensorConstant{1}, Elemwise{Add}[(0, 1)].0)
   2.3%    91.1%       0.682s       1.36e-01s      5   121   Join(TensorConstant{0}, Elemwise{sub,no_inplace}.0, Elemwise{sub,no_inplace}.0)
   1.5%    92.6%       0.430s       8.61e-02s      5   166   Elemwise{mul,no_inplace}(Subtensor{int64::}.0, Positions of the points to interpolate.T)
   1.4%    94.0%       0.410s       8.19e-02s      5   109   Elemwise{sub,no_inplace}(InplaceDimShuffle{0,x}.0, InplaceDimShuffle{1,0}.0)
   1.4%    95.4%       0.400s       8.00e-02s      5   108   Elemwise{sub,no_inplace}(InplaceDimShuffle{0,x}.0, InplaceDimShuffle{1,0}.0)
   1.2%    96.6%       0.361s       7.22e-02s      5    43   Dot22Scalar(Reference points for every layer, Positions of the points to interpolate.T, TensorConstant{2.0})
   1.2%    97.8%       0.360s       7.20e-02s      5   167   Sum{axis=[0], acc_dtype=float64}(Elemwise{Composite{(i0 * ((LT(Composite{sqrt(((i0 + i1) - i2))}(i1, i2, i3), i4) * ((i5 + (i6 * Composite{(sqr(i0) * i0)}((Composite{sqrt(((i0 + i1) - i2))}(i1, i2, i3) / i4))) + (i7 * Composite{((sqr(sqr(i0)) * sqr(i0)) * i0)}((Composite{sqrt(((i0 + i1) - i2))}(i1, i2, i3) / i4)))) - ((i8 * sqr((Composite{sqrt(((i0 + i1) - i2))}(i1, i2, i3) / i4))) + (i9 * Composite{(sqr(sqr(i0)) * i0)}((Composite{sqrt(((i0 + i1) - 
   1.2%    99.0%       0.355s       7.10e-02s      5    44   Dot22Scalar(Rest of the points of the layers, Positions of the points to interpolate.T, TensorConstant{2.0})
   0.3%    99.4%       0.094s       1.87e-02s      5   169   Sum{axis=[0], acc_dtype=float64}(Elemwise{mul,no_inplace}.0)
   0.3%    99.7%       0.090s       1.80e-02s      5    45   Sum{axis=[1], acc_dtype=float64}(Elemwise{sqr,no_inplace}.0)
   0.2%    99.8%       0.054s       1.09e-02s      5   170   Elemwise{Composite{(i0 + ((i1 * i2) / i3) + i4)}}[(0, 0)](Sum{axis=[0], acc_dtype=float64}.0, TensorConstant{(1,) of -1.75}, Sum{axis=[0], acc_dtype=float64}.0, InplaceDimShuffle{x}.0, Sum{axis=[0], acc_dtype=float64}.0)
   0.1%   100.0%       0.034s       6.79e-03s      5    11   Elemwise{sqr,no_inplace}(Positions of the points to interpolate)
   0.0%   100.0%       0.008s       1.55e-03s      5   155   MatrixInverse(IncSubtensor{InplaceSet;int64::, int64:int64:}.0)
   0.0%   100.0%       0.000s       9.94e-05s      5   134   Sum{axis=[1], acc_dtype=float64}(Elemwise{Sqr}[(0, 0)].0)
   0.0%   100.0%       0.000s       9.12e-05s      5   100   Alloc(Elemwise{sub,no_inplace}.0, TensorConstant{1}, TensorConstant{2}, Elemwise{Composite{Switch(EQ(i0, i1), (i0 // (-i0)), i0)}}.0, Shape_i{0}.0)
   ... (remaining 151 Apply instances account for 0.01%(0.00s) 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
  - Try installing amdlibm and set the Theano flag lib.amdlibm=True. This speeds up only some Elemwise operation.

In [ ]:

    




    

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 [32]:

    

%%timeit
sol = test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref);


    

    




10 loops, best of 3: 19.3 ms per loop

In [33]:

    

test.geoMigueller.profile.summary()


    

    




Function profiling
==================
  Message: /home/miguel/PycharmProjects/GeMpy/GeoMig.py:1346
  Time in 48 calls to Function.__call__: 1.300695e+00s
  Time in Function.fn.__call__: 1.296489e+00s (99.677%)
  Time in thunks: 1.084233e+00s (83.358%)
  Total compile time: 2.702078e+00s
    Number of Apply nodes: 295
    Theano Optimizer time: 2.239654e+00s
       Theano validate time: 1.687617e-01s
    Theano Linker time (includes C, CUDA code generation/compiling): 3.944879e-01s
       Import time 6.902671e-02s

Time in all call to theano.grad() 0.000000e+00s
Time since theano import 533.471s
Class
---
<% time> <sum %> <apply time> <time per call> <type> <#call> <#apply> <Class name>
  37.5%    37.5%       0.407s       3.14e-04s     C     1296      27   theano.tensor.elemwise.Elemwise
  26.3%    63.8%       0.285s       5.93e-04s     C      480      10   theano.tensor.blas.Dot22Scalar
   9.7%    73.5%       0.105s       9.96e-05s     C     1056      22   theano.sandbox.cuda.basic_ops.GpuFromHost
   9.0%    82.5%       0.098s       2.84e-05s     C     3456      72   theano.sandbox.cuda.basic_ops.GpuElemwise
   3.4%    85.9%       0.037s       1.10e-04s     C      336       7   theano.sandbox.cuda.basic_ops.GpuAlloc
   3.3%    89.2%       0.035s       5.26e-05s     C      672      14   theano.sandbox.cuda.basic_ops.HostFromGpu
   2.5%    91.6%       0.027s       3.72e-05s     C      720      15   theano.sandbox.cuda.basic_ops.GpuIncSubtensor
   2.3%    93.9%       0.025s       3.03e-05s     C      816      17   theano.sandbox.cuda.basic_ops.GpuReshape
   2.2%    96.1%       0.024s       1.26e-04s     C      192       4   theano.tensor.elemwise.Sum
   1.4%    97.5%       0.015s       5.22e-05s     C      288       6   theano.sandbox.cuda.basic_ops.GpuJoin
   1.0%    98.5%       0.011s       7.33e-05s     C      144       3   theano.sandbox.cuda.basic_ops.GpuCAReduce
   0.6%    99.1%       0.007s       1.39e-04s     Py      48       1   theano.tensor.nlinalg.MatrixInverse
   0.1%    99.3%       0.002s       3.30e-05s     Py      48       1   theano.tensor.extra_ops.FillDiagonal
   0.1%    99.4%       0.002s       9.38e-07s     C     1680      35   theano.sandbox.cuda.basic_ops.GpuDimShuffle
   0.1%    99.5%       0.001s       2.91e-05s     C       48       1   theano.sandbox.cuda.basic_ops.GpuAllocEmpty
   0.1%    99.7%       0.001s       1.70e-06s     C      720      15   theano.sandbox.cuda.basic_ops.GpuSubtensor
   0.1%    99.7%       0.001s       2.79e-06s     C      288       6   theano.tensor.subtensor.IncSubtensor
   0.1%    99.8%       0.001s       1.03e-06s     C      624      13   theano.tensor.opt.MakeVector
   0.0%    99.8%       0.000s       2.53e-06s     C      192       4   theano.tensor.elemwise.DimShuffle
   0.0%    99.9%       0.000s       4.54e-06s     C       96       2   theano.tensor.basic.Alloc
   ... (remaining 5 Classes account for   0.13%(0.00s) of the runtime)

Ops
---
<% time> <sum %> <apply time> <time per call> <type> <#call> <#apply> <Op name>
  32.6%    32.6%       0.353s       7.36e-04s     C      480       10   Elemwise{Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}}
  26.3%    58.9%       0.285s       5.93e-04s     C      480       10   Dot22Scalar
   9.7%    68.6%       0.105s       9.96e-05s     C     1056       22   GpuFromHost
   3.4%    72.0%       0.037s       1.94e-04s     C      192        4   Elemwise{Cast{float64}}
   3.3%    75.3%       0.035s       5.26e-05s     C      672       14   HostFromGpu
   2.8%    78.0%       0.030s       1.26e-04s     C      240        5   GpuAlloc
   2.3%    80.3%       0.025s       3.03e-05s     C      816       17   GpuReshape{2}
   2.2%    82.5%       0.024s       1.26e-04s     C      192        4   Sum{axis=[1], acc_dtype=float64}
   2.1%    84.7%       0.023s       5.97e-05s     C      384        8   GpuElemwise{Composite{(i0 * (i1 ** i2))},no_inplace}
   2.0%    86.7%       0.022s       4.60e-05s     C      480       10   GpuElemwise{Composite{Cast{float32}(LT(i0, i1))},no_inplace}
   1.9%    88.6%       0.021s       3.12e-05s     C      672       14   GpuElemwise{sub,no_inplace}
   1.4%    90.0%       0.015s       7.92e-05s     C      192        4   Elemwise{Sqr}[(0, 0)]
   1.4%    91.4%       0.015s       4.51e-05s     C      336        7   GpuIncSubtensor{InplaceSet;int64:int64:, int64:int64:}
   1.4%    92.8%       0.015s       5.22e-05s     C      288        6   GpuJoin
   1.0%    93.8%       0.011s       7.33e-05s     C      144        3   GpuCAReduce{add}{1,0}
   0.8%    94.6%       0.009s       9.44e-05s     C       96        2   GpuElemwise{Composite{(i0 * sqr(i1))},no_inplace}
   0.6%    95.3%       0.007s       1.43e-04s     C       48        1   GpuElemwise{Composite{(i0 * i1 * i2 * ((i3 - i4) ** i5) * ((i6 + (i7 * i3 * i4)) + (i5 * sqr(i4))))}}[(0, 1)]
   0.6%    95.9%       0.007s       6.98e-05s     C       96        2   GpuAlloc{memset_0=True}
   0.6%    96.5%       0.007s       1.39e-04s     Py      48        1   MatrixInverse
   0.4%    96.9%       0.004s       4.11e-05s     C       96        2   GpuIncSubtensor{InplaceSet;int64:int64:, int64::}
   ... (remaining 67 Ops account for   3.14%(0.03s) of the runtime)

Apply
------
<% time> <sum %> <apply time> <time per call> <#call> <id> <Apply name>
  20.5%    20.5%       0.223s       4.64e-03s     48   207   Elemwise{Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}}(Reshape{2}.0, Reshape{2}.0, Dot22Scalar.0)
  12.6%    33.1%       0.137s       2.85e-03s     48   160   Dot22Scalar(Elemwise{Cast{float64}}.0, InplaceDimShuffle{1,0}.0, TensorConstant{2.0})
   7.0%    40.1%       0.076s       1.57e-03s     48   108   Dot22Scalar(Elemwise{Cast{float64}}.0, InplaceDimShuffle{1,0}.0, TensorConstant{2.0})
   6.6%    46.7%       0.071s       1.48e-03s     48   199   Elemwise{Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}}(Reshape{2}.0, Reshape{2}.0, Dot22Scalar.0)
   6.1%    52.8%       0.066s       1.38e-03s     48   107   Dot22Scalar(Elemwise{Cast{float64}}.0, InplaceDimShuffle{1,0}.0, TensorConstant{2.0})
   5.4%    58.2%       0.059s       1.23e-03s     48   203   Elemwise{Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}}(Reshape{2}.0, Reshape{2}.0, Dot22Scalar.0)
   3.9%    62.1%       0.042s       8.76e-04s     48   217   GpuFromHost(Elemwise{Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}}.0)
   3.4%    65.5%       0.037s       7.70e-04s     48    38   Elemwise{Cast{float64}}(HostFromGpu.0)
   2.2%    67.7%       0.024s       4.98e-04s     48   179   Sum{axis=[1], acc_dtype=float64}(Elemwise{Sqr}[(0, 0)].0)
   2.0%    69.7%       0.021s       4.46e-04s     48   280   GpuAlloc(GpuGemv{inplace}.0, Shape_i{0}.0, TensorConstant{1}, TensorConstant{1}, Elemwise{add,no_inplace}.0)
   1.8%    71.5%       0.020s       4.08e-04s     48   213   GpuFromHost(Elemwise{Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}}.0)
   1.8%    73.2%       0.019s       3.98e-04s     48    10   HostFromGpu(Positions of the points to interpolate)
   1.4%    74.7%       0.015s       3.18e-04s     48     9   GpuFromHost(Position of the dips)
   1.4%    76.0%       0.015s       3.14e-04s     48   168   Elemwise{Sqr}[(0, 0)](Elemwise{Cast{float64}}.0)
   1.2%    77.2%       0.013s       2.68e-04s     48   209   GpuFromHost(Elemwise{Composite{Cast{float32}(sqrt(((i0 + i1) - i2)))}}.0)
   0.7%    77.9%       0.007s       1.48e-04s     48   235   GpuElemwise{Composite{Cast{float32}(LT(i0, i1))},no_inplace}(GpuFromHost.0, GpuDimShuffle{x,x}.0)
   0.6%    78.5%       0.007s       1.43e-04s     48   287   GpuElemwise{Composite{(i0 * i1 * i2 * ((i3 - i4) ** i5) * ((i6 + (i7 * i3 * i4)) + (i5 * sqr(i4))))}}[(0, 1)](GpuSubtensor{:int64:}.0, GpuJoin.0, GpuElemwise{Composite{Cast{float32}(LT(i0, i1))},no_inplace}.0, GpuDimShuffle{x,x}.0, GpuFromHost.0, CudaNdarrayConstant{[[ 3.]]}, GpuElemwise{Mul}[(0, 1)].0, CudaNdarrayConstant{[[ 9.]]})
   0.6%    79.1%       0.007s       1.42e-04s     48   247   GpuElemwise{Composite{(i0 * (i1 ** i2))},no_inplace}(CudaNdarrayConstant{[[ 0.75]]}, GpuElemwise{TrueDiv}[(0, 0)].0, CudaNdarrayConstant{[[ 7.]]})
   0.6%    79.8%       0.007s       1.39e-04s     48   276   MatrixInverse(HostFromGpu.0)
   0.5%    80.3%       0.006s       1.21e-04s     48   294   HostFromGpu(GpuElemwise{Composite{((((i0 * i1) / i2) + i3) + i4)}}[(0, 1)].0)
   ... (remaining 275 Apply instances account for 19.71%(0.21s) 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.

In [30]:

    

importlib.reload(GeoMig)
test = GeoMig.GeoMigSim_pro2()


    

In [42]:

    




    

    
Out[42]:





array([[ -5.88888884e-01,  -0.00000000e+00,  -0.00000000e+00, ...,
          0.00000000e+00,   1.00000000e+00,   0.00000000e+00],
       [ -0.00000000e+00,  -5.88888884e-01,   4.42373231e-02, ...,
          0.00000000e+00,   1.00000000e+00,   0.00000000e+00],
       [ -0.00000000e+00,   2.12696299e-01,  -5.88888884e-01, ...,
          0.00000000e+00,   1.00000000e+00,   0.00000000e+00],
       ..., 
       [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00, ...,
          2.00000000e+00,  -6.06459351e+02,  -6.13501053e+01],
       [  1.00000000e+00,   1.00000000e+00,   1.00000000e+00, ...,
         -6.06459351e+02,   0.00000000e+00,   0.00000000e+00],
       [  0.00000000e+00,   0.00000000e+00,   0.00000000e+00, ...,
         -6.13501053e+01,   0.00000000e+00,   0.00000000e+00]], dtype=float32)

In [17]:

    




    

    




---------------------------------------------------------------------------
Exception                                 Traceback (most recent call last)
<ipython-input-17-b0516f7436af> in <module>()
----> 1 theano.config.device="cpu"

/home/miguel/anaconda3/lib/python3.5/site-packages/theano/configparser.py in __set__(self, cls, val)
    329         if not self.allow_override and hasattr(self, 'val'):
    330             raise Exception(
--> 331                 "Can't change the value of this config parameter "
    332                 "after initialization!")
    333         # print "SETTING PARAM", self.fullname,(cls), val

Exception: Can't change the value of this config parameter after initialization!

In [1]:

    

from theano import function, config, shared, sandbox
import theano.tensor as T
import numpy
import time

vlen = 10 * 30 * 768  # 10 x #cores x # threads per core
iters = 1000

rng = numpy.random.RandomState(22)
x = shared(numpy.asarray(rng.rand(vlen), config.floatX))
f = function([], T.exp(x))
print(f.maker.fgraph.toposort())
t0 = time.time()
for i in range(iters):
    r = f()
t1 = time.time()
print("Looping %d times took %f seconds" % (iters, t1 - t0))
print("Result is %s" % (r,))
if numpy.any([isinstance(x.op, T.Elemwise) for x in f.maker.fgraph.toposort()]):
    print('Used the cpu')
else:
    print('Used the gpu')


    

    




[Elemwise{exp,no_inplace}(<TensorType(float32, vector)>)]
Looping 1000 times took 2.271379 seconds
Result is [ 1.23178029  1.61879337  1.52278066 ...,  2.20771813  2.29967761
  1.62323284]
Used the cpu

In [1]:

    

from theano import function, config, shared, sandbox
import theano.tensor as T
import numpy
import time

vlen = 10 * 30 * 768  # 10 x #cores x # threads per core
iters = 1000

rng = numpy.random.RandomState(22)
x = shared(numpy.asarray(rng.rand(vlen), config.floatX))
f = function([], T.exp(x))
print(f.maker.fgraph.toposort())
t0 = time.time()
for i in range(iters):
    r = f()
t1 = time.time()
print("Looping %d times took %f seconds" % (iters, t1 - t0))
print("Result is %s" % (r,))
if numpy.any([isinstance(x.op, T.Elemwise) for x in f.maker.fgraph.toposort()]):
    print('Used the cpu')
else:
    print('Used the gpu')


    

    




Using gpu device 0: GeForce GTX 970 (CNMeM is disabled, cuDNN not available)






    




[GpuElemwise{exp,no_inplace}(<CudaNdarrayType(float32, vector)>), HostFromGpu(GpuElemwise{exp,no_inplace}.0)]
Looping 1000 times took 0.353415 seconds
Result is [ 1.23178029  1.61879349  1.52278066 ...,  2.20771813  2.29967761
  1.62323296]
Used the gpu

In [18]:

    

from theano import function, config, shared, sandbox
import theano.tensor as T
import numpy
import time

vlen = 10 * 30 * 768  # 10 x #cores x # threads per core
iters = 1000

rng = numpy.random.RandomState(22)
x = shared(numpy.asarray(rng.rand(vlen), config.floatX))
f = function([], T.exp(x))
print(f.maker.fgraph.toposort())
t0 = time.time()
for i in range(iters):
    r = f()
t1 = time.time()
print("Looping %d times took %f seconds" % (iters, t1 - t0))
print("Result is %s" % (r,))
if numpy.any([isinstance(x.op, T.Elemwise) for x in f.maker.fgraph.toposort()]):
    print('Used the cpu')
else:
    print('Used the gpu')


    

    




[Elemwise{exp,no_inplace}(<TensorType(float32, vector)>)]
Looping 1000 times took 2.291412 seconds
Result is [ 1.23178029  1.61879337  1.52278066 ...,  2.20771813  2.29967761
  1.62323284]
Used the cpu

In [ ]:

    




    

In [759]:

    

np.set_printoptions(precision=2)
test.geoMigueller(dips,dips_angles,rest, ref)[1]


    

    
Out[759]:





array([[-0.59,  0.08,  0.  ,  0.07],
       [ 0.08, -0.59,  0.07,  0.  ],
       [ 0.  ,  0.12, -0.59,  0.13],
       [ 0.07,  0.  ,  0.13, -0.59]])

In [751]:

    

T.fill_diagonal?


    

In [758]:

    

import matplotlib.pyplot as plt
% matplotlib inline
dip_pos_1_v = np.array([np.cos(np.deg2rad(dip_angle_1))*1,
                        np.sin(np.deg2rad(dip_angle_1))]) + dip_pos_1

dip_pos_2_v = np.array([np.cos(np.deg2rad(dip_angle_2))*1, 
                        np.sin(np.deg2rad(dip_angle_2))]) + dip_pos_2

plt.arrow(dip_pos_1[0],dip_pos_1[1], dip_pos_1_v[0]-dip_pos_1[0],
          dip_pos_1_v[1]-dip_pos_1[1], head_width = 0.2)
plt.arrow(dip_pos_2[0],dip_pos_2[1],dip_pos_2_v[0]-dip_pos_2[0], 
          dip_pos_2_v[1]-dip_pos_2[1], head_width = 0.2)

plt.plot(layer_1[:,0],layer_1[:,1], "o")
plt.plot(layer_2[:,0],layer_2[:,1], "o")

plt.plot(layer_1[:,0],layer_1[:,1], )
plt.plot(layer_2[:,0],layer_2[:,1], )

plt.contour( sol.reshape(50,50) ,30,extent = (0,10,0,10) )
#plt.colorbar()
#plt.xlim(0,10)
#plt.ylim(0,10)
plt.title("GeoBulleter v 0.1")
print (dip_pos_1_v, dip_pos_2_v, layer_1)


    

    




[ 2.5   4.87] [ 6.34  3.94] [[ 1.  7.]
 [ 5.  7.]
 [ 6.  7.]
 [ 9.  8.]]






    

In [443]:

    

n = 10
#a = T.horizontal_stack(T.vertical_stack(T.ones(n),T.zeros(n)), T.vertical_stack(T.zeros(n), T.ones(n)))
a = T.zeros(n)

print (a.eval())
#U_G = T.horizontal_stack(([T.ones(n),T.zeros(n)],[T.zeros(n),T.ones(n)]))


    

    




[ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]

In [6]:

    

T.stack?ö+aeg


    

In [13]:

    

x_min = 0
x_max = 10
y_min = 0
y_max = 10
z_min = 0
z_max = 10
nx = 2
ny = 2
nz = 2

g = np.meshgrid(
    np.linspace(x_min, x_max, nx, dtype="float32"),
    np.linspace(y_min, y_max, ny, dtype="float32"),
    np.linspace(z_min, z_max, nz, dtype="float32"), indexing="ij"
)

np.vstack(map(np.ravel, g)).T.astype("float32")


    

    
Out[13]:





array([[  0.,   0.,   0.],
       [  0.,   0.,  10.],
       [  0.,  10.,   0.],
       [  0.,  10.,  10.],
       [ 10.,   0.,   0.],
       [ 10.,   0.,  10.],
       [ 10.,  10.,   0.],
       [ 10.,  10.,  10.]], dtype=float32)

In [18]:

    

map(np.ravel, g)


    

    
Out[18]:





<map at 0x7fa0265314e0>

In [22]:

    

np.ravel(g, order = "F")


    

    
Out[22]:





array([  0.,   0.,   0.,  10.,   0.,   0.,   0.,  10.,   0.,  10.,  10.,   0.,   0.,   0.,  10.,  10.,   0.,  10.,   0.,  10.,
        10.,  10.,  10.,  10.])

In [23]:

    

g


    

    
Out[23]:





[array([[[  0.,   0.],
         [  0.,   0.]],
 
        [[ 10.,  10.],
         [ 10.,  10.]]]), array([[[  0.,   0.],
         [ 10.,  10.]],
 
        [[  0.,   0.],
         [ 10.,  10.]]]), array([[[  0.,  10.],
         [  0.,  10.]],
 
        [[  0.,  10.],
         [  0.,  10.]]])]

In [25]:

    

np.transpose?


    

In [ ]: