import importlib
import theano.tensor as T
import sys, os
sys.path.append("/home/bl3/PycharmProjects/GeMpy/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 = 8, linewidth= 130, suppress = True)
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
from matplotlib import cm
%matplotlib inline
os.getcwd()
'/home/bl3/PycharmProjects/GeMpy/Prototype Notebook'
import geogrid
a = geogrid.GeoGrid()
a.xmax = 10
a.xmin = 0
a.ymax = 10
a.ymin = 0
a.zmax = 10
a.zmin = 0
a.define_regular_grid(50,50,50)
test = GeoMig.GeoMigSim_pro2(c_o = np.float32(-0.2),range = np.float32(13))
test.create_regular_grid_3D(0,10,0,10,0,10,50,50,50)
test.theano_set_3D()
/home/bl3/anaconda3/lib/python3.5/site-packages/ipykernel/__main__.py:4: DeprecationWarning: stack(*tensors) interface is deprecated, use stack(tensors, axis=0) instead.
layer_1 = np.array([[1,5,7],[4,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")
dip_pos_1 = np.array([2,5,5], dtype = "float32")
dip_pos_2 = np.array([8.,5.,5], dtype = "float32")
#dip_pos_3 = np.array([8,4,5], dtype = "float32")
dip_angle_1 = float(0)
dip_angle_2 = float(0)
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([90, 90], 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)
numpy.ndarray
dips
array([[ 2., 5., 5.], [ 8., 5., 5.]], dtype=float32)
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[0],G_y[1],G_z[0],G_x[1], G_y[1],G_z[1]
(0.0, -0.0, 1.0, 0.0, -0.0, 1.0)
test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[7]
array([[ 0., 0., 1., 1., 1., 1.], [ 0., 0., 1., 1., 1., 1.], [ 1., 1., 0., 0., 1., 1.], [ 1., 1., 0., 0., 1., 1.], [ 1., 1., 1., 1., 0., 0.], [ 1., 1., 1., 1., 0., 0.]], dtype=float32)
test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[2]
array([[-0.2 , -0.04439974, 0. , -0.12408967, 0. , -0.12408967], [-0.04439974, -0.2 , -0.12408967, 0. , -0.12408967, 0. ], [ 0. , -0.12408967, -0.2 , 0. , 0. , -0.12408967], [-0.12408967, 0. , 0. , -0.2 , -0.12408967, 0. ], [ 0. , -0.12408967, 0. , -0.12408967, -0.2 , 0. ], [-0.12408967, 0. , -0.12408967, 0. , 0. , -0.2 ]])
sol = test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[0].reshape(50,50,50, order = "F")
#sol = np.swapaxes(sol,0,1)
dip_pos_1_v = [test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[-3][0],
test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[-1][0],]+ dip_pos_1[0::2]
dip_pos_2_v = [test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[-3][1],
test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[-1][1],]+ dip_pos_2[0::2]
print (test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[-1])
dip_pos_1_v, dip_pos_2_v, dip_pos_2;
[ 0.70710678 0.70710678]
h,j,k =sol[5,10,35], sol[25,5,5], sol[30,15,-25]
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])
0.733862734467 3.74963487998 4.40451413599 6.48626327603 0.711641729939 3.54708831021 6.61943538563
interfaces_aux = test.geoMigueller(dips,dips_angles,azimuths,polarity,
rest, ref)[0]
h = sol[10,25,30]# interfaces_aux[np.argmin(abs((test.grid - ref[0]).sum(1)))]
k = sol[30,25,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
(1.5142873563459891, 4.3996999023327437, 6.6194353856317658)
dip_pos_1[1]
5.0
a.grid = sol
a.plot_section('z',cell_pos=27,colorbar = True, alpha = 0.8, cmap = 'coolwarm_r',
figsize=(6,6),interpolation= 'nearest' ,ve = 1, geomod_coord= False, contour = False,
)
vertices
--------------------------------------------------------------------------- NameError Traceback (most recent call last) <ipython-input-56-6a2207daec15> in <module>() ----> 1 vertices NameError: name 'vertices' is not defined
import plotly.plotly as py
from plotly.graph_objs import *
from plotly.tools import FigureFactory as FF
import plotly.graph_objs as go
import plotly.plotly as py
py.sign_in('leguiom', 'wphx7fdchi')
import numpy as np
from skimage import measure
# ========
# PLOTING ISOSURFACES
h = sol[35,25,25]
vertices, simplices = measure.marching_cubes(sol, h, spacing = (1,1,1),
gradient_direction = "descent")
x,y,z = zip(*vertices)
#colormap=['rgb(255,105,180)','rgb(255,255,51)','rgb(0,191,255)']
fig = FF.create_trisurf(x=x,
y=y,
z=z,
plot_edges=False,
# colormap=colormap,
simplices=simplices,
title="Isosurface")
"""
vertices, simplices = measure.marching_cubes(sol,k, spacing = (1,1,1),
gradient_direction = "descent")
x,y,z = zip(*vertices)
#colormap=['rgb(255,105,180)','rgb(255,255,51)','rgb(0,191,255)']
fig2 = FF.create_trisurf(x=x,
y=y,
z=z,
plot_edges=False,
# colormap=colormap,
simplices=simplices,
title="Isosurface",
)
vertices, simplices = measure.marching_cubes(sol,j, spacing = (1,1,1),
gradient_direction = "descent")
x,y,z = zip(*vertices)
#colormap=['rgb(255,105,180)','rgb(255,255,51)','rgb(0,191,255)']
fig3 = FF.create_trisurf(x=x,
y=y,
z=z,
plot_edges=False,
# colormap=colormap,
simplices=simplices,
title="Isosurface",
)
"""
# ============
# Plotting interface points
import plotly.graph_objs as go
trace1 = go.Scatter3d(
x=layers[0][:,0]*5,
y=layers[0][:,1]*5,
z=layers[0][:,2]*5,
mode='markers',
marker=dict(
size=12,
line=dict(
color='rgba(217, 217, 217, 0.14)',
width=0.5
),
opacity=0.8
)
)
trace12 = go.Scatter3d(
x=layers[1][:,0]*5,
y=layers[1][:,1]*5,
z=layers[1][:,2]*5,
mode='markers',
marker=dict(
size=12,
line=dict(
color='rgba(217, 217, 217, 0.14)',
width=0.5
),
opacity=0.8
)
)
# ================
# Plotting traces
a = np.vstack((dips[:,0],dips[:,0]+G_x))
b = np.vstack((dips[:,1],dips[:,1]+G_y))
c = np.vstack((dips[:,2],dips[:,2]+G_z))
trace2 = go.Scatter3d(
x=a[:,0]*5, y=b[:,0]*5, z=c[:,0]*5,
marker=dict(
size=4,
color=z,
colorscale='Viridis',
),
line=dict(
color='#1f77b4',
width=1
)
)
trace3 = go.Scatter3d(
x=a[:,1]*5, y=b[:,1]*5, z=c[:,1]*5,
marker=dict(
size=4,
color=z,
colorscale='Viridis',
),
line=dict(
color='#1f77b4',
width=1
)
)
fig.layout['scene']['zaxis']['range'] = [0,50]
#fig.data.append(fig2.data[0])
#fig.data.append(fig3.data[0])
fig.data.append(trace12)
fig.data.append(trace1)
fig.data.append(trace2)
fig.data.append(trace3)
fig.data[1]['visible'] = False
py.iplot(fig)
a = np.vstack((dips[:,0],dips[:,0]+G_x))
b = np.vstack((dips[:,1],dips[:,1]+G_y))
c = np.vstack((dips[:,2],dips[:,2]+G_z))
a,b,c
(array([[ 2. , 6. ], [ 2.14454389, 5.54548073]], dtype=float32), array([[ 4. , 3. ], [ 4.39713144, 2.45832467]], dtype=float32), array([[ 6. , 5. ], [ 6.9063077 , 5.70710659]], dtype=float32))
fig.layout["shapes"]
[]
from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot
init_notebook_mode(connected=False)
"""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'
len(x)
51
surf_eq.min()
-63.0
np.min(z)
0.0
layers[0][:,0]
array([ 1., 5., 6., 9.], dtype=float32)
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
a
array([[ 2. , 6. ], [ 2.03, 6.15]], dtype=float32)
data = [trace1, trace2]
layout = go.Layout(
xaxis=dict(
range=[2, 5]
),
yaxis=dict(
range=[2, 5]
)
)
fig = go.Figure(data=data, layout=layout)
import lxml
lxml??
# 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;
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]]
import matplotlib.pyplot as plt
% matplotlib inline
plt.contour( sol.reshape(100,100) ,30,extent = (0,10,0,10) )
<matplotlib.contour.QuadContourSet at 0x7ff55335e7f0>
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.]]
%%timeit
sol = test.geoMigueller(dips,dips_angles,rest, ref)[0]
1 loop, best of 3: 5.81 s per loop
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.
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]
'/home/bl3/PycharmProjects/GeMpy/GeoMigueller.vtr'
%%timeit
sol = test.geoMigueller(dips,dips_angles,rest, ref)[0]
10 loops, best of 3: 171 ms per loop
test.geoMigueller.profile.summary()
Function profiling ================== Message: /home/miguel/PycharmProjects/GeMpy/GeoMig.py:562 Time in 43 calls to Function.__call__: 7.149101e+00s Time in Function.fn.__call__: 7.146105e+00s (99.958%) Time in thunks: 4.351497e+00s (60.868%) Total compile time: 2.354024e+00s Number of Apply nodes: 227 Theano Optimizer time: 2.021928e+00s Theano validate time: 1.219668e-01s Theano Linker time (includes C, CUDA code generation/compiling): 2.757676e-01s Import time 0.000000e+00s Time in all call to theano.grad() 0.000000e+00s Time since theano import 281.359s Class --- <% time> <sum %> <apply time> <time per call> <type> <#call> <#apply> <Class name> 86.8% 86.8% 3.779s 1.73e-02s C 219 5 theano.sandbox.cuda.basic_ops.HostFromGpu 6.5% 93.3% 0.282s 7.92e-04s C 356 8 theano.sandbox.cuda.basic_ops.GpuAlloc 2.4% 95.7% 0.103s 3.19e-05s C 3235 73 theano.sandbox.cuda.basic_ops.GpuElemwise 1.5% 97.2% 0.066s 2.48e-04s C 266 6 theano.sandbox.cuda.basic_ops.GpuJoin 0.7% 97.9% 0.030s 6.69e-05s C 446 10 theano.sandbox.cuda.blas.GpuDot22Scalar 0.5% 98.4% 0.021s 3.60e-05s C 583 13 theano.sandbox.cuda.basic_ops.GpuIncSubtensor 0.5% 98.8% 0.020s 6.53e-05s C 307 7 theano.sandbox.cuda.basic_ops.GpuCAReduce 0.4% 99.3% 0.019s 4.51e-04s Py 43 1 theano.tensor.nlinalg.MatrixInverse 0.3% 99.6% 0.014s 1.71e-05s C 847 19 theano.sandbox.cuda.basic_ops.GpuReshape 0.3% 99.9% 0.011s 3.20e-05s C 356 8 theano.sandbox.cuda.basic_ops.GpuFromHost 0.0% 99.9% 0.001s 8.66e-07s C 1385 31 theano.sandbox.cuda.basic_ops.GpuDimShuffle 0.0% 99.9% 0.001s 2.48e-05s Py 45 1 theano.tensor.extra_ops.FillDiagonal 0.0% 100.0% 0.001s 1.81e-06s C 485 11 theano.sandbox.cuda.basic_ops.GpuSubtensor 0.0% 100.0% 0.000s 8.73e-07s C 534 12 theano.tensor.opt.MakeVector 0.0% 100.0% 0.000s 1.27e-06s C 358 8 theano.tensor.elemwise.Elemwise 0.0% 100.0% 0.000s 8.67e-06s C 43 1 theano.sandbox.cuda.blas.GpuGemv 0.0% 100.0% 0.000s 8.26e-06s C 45 1 theano.sandbox.cuda.basic_ops.GpuAllocEmpty 0.0% 100.0% 0.000s 9.69e-07s C 266 6 theano.compile.ops.Shape_i 0.0% 100.0% 0.000s 7.56e-07s C 268 6 theano.tensor.basic.ScalarFromTensor ... (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> 86.8% 86.8% 3.779s 1.73e-02s C 219 5 HostFromGpu 6.3% 93.1% 0.274s 1.54e-03s C 178 4 GpuAlloc 1.5% 94.7% 0.066s 2.48e-04s C 266 6 GpuJoin 0.7% 95.3% 0.030s 6.69e-05s C 446 10 GpuDot22Scalar 0.6% 95.9% 0.026s 6.52e-05s C 401 9 GpuElemwise{sub,no_inplace} 0.6% 96.5% 0.025s 7.23e-05s C 352 8 GpuElemwise{Composite{(i0 * (i1 ** i2))},no_inplace} 0.5% 97.1% 0.023s 5.15e-05s C 444 10 GpuElemwise{Composite{Cast{float32}(LT(i0, i1))},no_inplace} 0.4% 97.5% 0.019s 4.51e-04s Py 43 1 MatrixInverse 0.3% 97.8% 0.014s 1.71e-05s C 847 19 GpuReshape{2} 0.3% 98.2% 0.014s 1.09e-04s C 129 3 GpuCAReduce{add}{1,0} 0.3% 98.4% 0.011s 3.20e-05s C 356 8 GpuFromHost 0.2% 98.7% 0.011s 1.25e-04s C 86 2 GpuElemwise{Composite{(i0 * sqr(i1))},no_inplace} 0.2% 98.9% 0.010s 4.40e-05s C 225 5 GpuIncSubtensor{InplaceSet;int64:int64:, int64:int64:} 0.2% 99.1% 0.008s 4.40e-05s C 178 4 GpuAlloc{memset_0=True} 0.1% 99.2% 0.006s 3.35e-05s C 178 4 GpuCAReduce{pre=sqr,red=add}{0,1} 0.1% 99.3% 0.004s 3.97e-05s C 90 2 GpuIncSubtensor{InplaceSet;int64:int64:, int64::} 0.1% 99.4% 0.003s 3.70e-05s C 90 2 GpuIncSubtensor{InplaceSet;int64::, int64:int64:} 0.1% 99.4% 0.003s 6.39e-06s C 444 10 GpuElemwise{Composite{sqrt(((i0 + i1) - i2))}}[(0, 2)] 0.0% 99.5% 0.002s 3.50e-05s C 45 1 GpuIncSubtensor{InplaceSet;int64, :int64:} 0.0% 99.5% 0.001s 1.57e-05s C 90 2 GpuElemwise{sqr,no_inplace} ... (remaining 56 Ops account for 0.50%(0.02s) of the runtime) Apply ------ <% time> <sum %> <apply time> <time per call> <#call> <id> <Apply name> 68.8% 68.8% 2.993s 6.96e-02s 43 215 HostFromGpu(GpuDimShuffle{1,0}.0) 17.4% 86.2% 0.756s 1.76e-02s 43 122 HostFromGpu(GpuElemwise{Composite{sqrt(((i0 + i1) - i2))}}[(0, 2)].0) 6.1% 92.3% 0.267s 6.21e-03s 43 211 GpuAlloc(GpuGemv{inplace}.0, Shape_i{0}.0, TensorConstant{1}, TensorConstant{1}, Elemwise{Add}[(0, 1)].0) 1.2% 93.5% 0.053s 1.18e-03s 45 145 GpuJoin(TensorConstant{0}, GpuElemwise{sub,no_inplace}.0, GpuElemwise{Sub}[(0, 0)].0) 0.7% 94.2% 0.028s 6.62e-04s 43 226 HostFromGpu(GpuElemwise{Composite{((((i0 * i1) / i2) + i3) + i4)}}[(0, 1)].0) 0.4% 94.6% 0.019s 4.51e-04s 43 208 MatrixInverse(HostFromGpu.0) 0.2% 94.8% 0.008s 1.84e-04s 43 112 GpuDot22Scalar(GpuReshape{2}.0, GpuDimShuffle{1,0}.0, TensorConstant{2.0}) 0.2% 95.0% 0.008s 1.82e-04s 43 141 GpuJoin(TensorConstant{0}, GpuElemwise{sub,no_inplace}.0, GpuElemwise{sub,no_inplace}.0) 0.2% 95.2% 0.008s 1.78e-04s 43 184 GpuElemwise{Composite{Cast{float32}(LT(i0, i1))},no_inplace}(GpuElemwise{Composite{sqrt(((i0 + i1) - i2))}}[(0, 2)].0, GpuDimShuffle{x,x}.0) 0.2% 95.3% 0.008s 1.67e-04s 45 32 GpuDot22Scalar(GpuFromHost.0, GpuDimShuffle{1,0}.0, TensorConstant{2.0}) 0.2% 95.5% 0.007s 1.54e-04s 43 118 GpuElemwise{sub,no_inplace}(GpuDimShuffle{0,x}.0, GpuDimShuffle{1,0}.0) 0.1% 95.6% 0.006s 1.46e-04s 43 154 GpuElemwise{Composite{(i0 * (i1 ** i2))},no_inplace}(CudaNdarrayConstant{[[ 8.75]]}, GpuElemwise{TrueDiv}[(0, 0)].0, CudaNdarrayConstant{[[ 3.]]}) 0.1% 95.8% 0.006s 1.40e-04s 43 153 GpuElemwise{Composite{(i0 * (i1 ** i2))},no_inplace}(CudaNdarrayConstant{[[ 0.75]]}, GpuElemwise{TrueDiv}[(0, 0)].0, CudaNdarrayConstant{[[ 7.]]}) 0.1% 95.9% 0.006s 1.39e-04s 43 123 GpuElemwise{Composite{Cast{float32}(LT(i0, i1))},no_inplace}(GpuElemwise{Composite{sqrt(((i0 + i1) - i2))}}[(0, 2)].0, GpuDimShuffle{x,x}.0) 0.1% 96.0% 0.006s 1.38e-04s 43 117 GpuElemwise{sub,no_inplace}(GpuDimShuffle{0,x}.0, GpuDimShuffle{1,0}.0) 0.1% 96.2% 0.006s 1.36e-04s 43 148 GpuElemwise{Composite{(i0 * (i1 ** i2))},no_inplace}(CudaNdarrayConstant{[[ 8.75]]}, GpuElemwise{TrueDiv}[(0, 0)].0, CudaNdarrayConstant{[[ 3.]]}) 0.1% 96.3% 0.006s 1.33e-04s 43 39 GpuDot22Scalar(GpuFromHost.0, GpuDimShuffle{1,0}.0, TensorConstant{2.0}) 0.1% 96.4% 0.006s 1.33e-04s 43 146 GpuElemwise{Composite{(i0 * sqr(i1))},no_inplace}(CudaNdarrayConstant{[[ 7.]]}, GpuElemwise{TrueDiv}[(0, 0)].0) 0.1% 96.6% 0.006s 1.25e-04s 45 86 GpuElemwise{sub,no_inplace}(GpuDimShuffle{x,0}.0, GpuReshape{2}.0) 0.1% 96.7% 0.005s 1.23e-04s 43 116 GpuElemwise{Composite{Cast{float32}(LT(i0, i1))},no_inplace}(GpuElemwise{Composite{sqrt(((i0 + i1) - i2))}}[(0, 2)].0, GpuDimShuffle{x,x}.0) ... (remaining 207 Apply instances account for 3.31%(0.14s) 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. Sorry, no tip for today.
importlib.reload(GeoMig)
test = GeoMig.GeoMigSim_pro2()
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)
--------------------------------------------------------------------------- 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!
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
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
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
np.set_printoptions(precision=2)
test.geoMigueller(dips,dips_angles,rest, ref)[1]
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]])
T.fill_diagonal?
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.]]
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.]
T.stack?รถ+aeg
_squared_euclidean_distances2 = T.sqrt(
(dips ** 2).sum(1).reshape((dips.shape[0], 1)) + (aux_Y ** 2).sum(1).reshape(
(1, aux_Y.shape[0])) - 2 * dips.dot(aux_Y.T))
_squared_euclidean_distances3 = T.sqrt(
(dips ** 2).sum(1).reshape((dips.shape[0], 1)) + (aux_X ** 2).sum(1).reshape(
(1, aux_X.shape[0])) - 2 * dips.dot(aux_X.T))
h3 = T.vertical_stack(
(dips[:, 0] - aux_Y[:, 0].reshape((aux_Y[:, 0].shape[0], 1))).T,
(dips[:, 1] - aux_Y[:, 1].reshape((aux_Y[:, 1].shape[0], 1))).T
)
h4 = T.vertical_stack(
(dips[:, 0] - aux_X[:, 0].reshape((aux_X[:, 0].shape[0], 1))).T,
(dips[:, 1] - aux_X[:, 1].reshape((aux_X[:, 1].shape[0], 1))).T)
r_3 = T.tile(_squared_euclidean_distances2, (2, 1)) # Careful with the number of dimensions
r_4 = T.tile(_squared_euclidean_distances3, (2, 1)) # Careful with the number of dimensions
_ans_d1_3 = (r_3 < self.a) * (
-7 * (self.a - r_3) ** 3 * r_3 * (8 * self.a ** 2 + 9 * self.a * r_3 + 3 * r_3 ** 2) * 1)
/ (4 * self.a ** 7)
_ans_d1_4 = (r_4 < self.a) * (
-7 * (self.a - r_4) ** 3 * r_4 * (8 * self.a ** 2 + 9 * self.a * r_4 + 3 * r_4 ** 2) * 1)
/ (4 * self.a ** 7)
_C_GI = (h3 / r_3 * _ans_d1_3 - h4 / r_4 * _ans_d1_4).T
self._f_CGI = theano.function([dips, aux_X, aux_Y], _C_GI)
np.ravel?
x_min = 0
x_max = 10
y_min =0
y_max = 10
z_min = 0
z_max =10
nx = 50
ny = 50
nz = 50
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")
)
np.ravel(g)
array([ 0. , 0. , 0. , ..., 9.59183693, 9.79591846, 10. ])
np.vstack(map(np.ravel, g))
array([[ 0. , 0. , 0. , ..., 10. , 10. , 10. ], [ 0. , 0. , 0. , ..., 10. , 10. , 10. ], [ 0. , 0.20408164, 0.40816328, ..., 9.59183693, 9.79591846, 10. ]])