In [3]:
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
In [ ]:
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)
In [ ]:
In [ ]:
In [ ]:
In [5]:
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()
In [1]:
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)
In [53]:
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]
Out[53]:
In [54]:
test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[7]
Out[54]:
In [55]:
test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[2]
Out[55]:
In [ ]:
In [42]:
sol = test.geoMigueller(dips,dips_angles,azimuths,polarity, rest, ref)[0].reshape(50,50,50, order = "F")
#sol = np.swapaxes(sol,0,1)
In [43]:
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;
In [ ]:
In [44]:
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])
In [45]:
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
Out[45]:
In [46]:
dip_pos_1[1]
Out[46]:
In [ ]:
In [ ]:
In [50]:
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,
)
In [56]:
vertices
In [131]:
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)
Out[131]:
In [29]:
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
Out[29]:
In [128]:
fig.layout["shapes"]
Out[128]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [120]:
from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot
init_notebook_mode(connected=False)
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})
In [201]:
len(x)
Out[201]:
In [202]:
surf_eq.min()
Out[202]:
In [203]:
np.min(z)
Out[203]:
In [204]:
layers[0][:,0]
Out[204]:
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]:
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()
In [11]:
import matplotlib.pyplot as plt
% matplotlib inline
plt.contour( sol.reshape(100,100) ,30,extent = (0,10,0,10) )
Out[11]:
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)
In [17]:
%%timeit
sol = test.geoMigueller(dips,dips_angles,rest, ref)[0]
In [18]:
test.geoMigueller.profile.summary()
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})
Out[23]:
In [19]:
%%timeit
sol = test.geoMigueller(dips,dips_angles,rest, ref)[0]
In [23]:
test.geoMigueller.profile.summary()
In [30]:
importlib.reload(GeoMig)
test = GeoMig.GeoMigSim_pro2()
In [42]:
Out[42]:
In [17]:
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')
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')
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')
In [ ]:
In [759]:
np.set_printoptions(precision=2)
test.geoMigueller(dips,dips_angles,rest, ref)[1]
Out[759]:
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)
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)]))
In [6]:
T.stack?ö+aeg
In [ ]:
_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)
In [ ]:
In [ ]:
In [137]:
np.ravel?
In [142]:
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)
Out[142]:
In [145]:
np.vstack(map(np.ravel, g))
Out[145]:
In [ ]: