In [1]:
%matplotlib inline
import matplotlib as mpl
import matplotlib.pyplot as plt
mpl.rc('font', size=15)
mpl.rc('figure', figsize=(8, 5))
import numpy as np
import scipy.signal as sig
import keras
from keras.layers import Input, Dense, Activation, Dropout
from keras.models import Model
from keras.layers.advanced_activations import ELU
from keras import regularizers
from keras.models import load_model
from keras.initializers import glorot_normal, glorot_uniform
from keras.optimizers import Adam
from mpl_toolkits.basemap import Basemap
from matplotlib.path import Path
In [2]:
# fix random seed for reproducibility
np.random.seed(7)
# Load data and exclude nan value
data = np.genfromtxt('IRISwest.txt')
In [3]:
#####maybe if we cut out the IRIS P waves we can get better results##
#going to chooose 6000
eqgpstime = data[:,1]
peakgpstime = data[:,17]
arrivaltime = np.subtract(peakgpstime,eqgpstime)
distance = data[:,13]
Velocity = np.divide(distance, arrivaltime)
pwaveomit = 6000
Velocity1 = Velocity[Velocity<6000]
data = data[Velocity<6000]
print (len(data))
In [4]:
########cutting to only use western hemi data, not needed for IRISwest.txt##########
######side note for sky, refernces w/ [:,] instead of np.array
#print(len(data1))
#eq_lon1 =data1[:,11]
#print(eq_lon1)
#data = data1[(eq_lon1>=-180) & (eq_lon1<=-30)]
#print(len(data))
#########cutting out ocean points ############
eq_lat2 = data[:,11]
eq_lon2 = data[:,12]
map1 = Basemap(projection='aeqd', lon_0 = 10, lat_0 = 50, resolution='h')
lats = eq_lat2 #[:100] \
lons = eq_lon2
x, y = map1(lons, lats)
locations = np.c_[x, y]
polygons = [Path(p.boundary) for p in map1.landpolygons]
result = np.zeros(len(locations), dtype=bool)
for polygon in polygons:
result += np.array(polygon.contains_points(locations))
eq_lat1=lats[result]
eq_lon1=lons[result]
print (len(data))
print (result)
data =data[result]
print (len(data))
In [5]:
# Extract X and y and divide into train, val, and test set
X = data[:, [2, 11, 12, 13, 14, 15]] #side note nikils has log10distnace maybe we should try that
# y = (data[:, 24] - data[:, 0])
y = np.log10(data[:, 18]) #(tri's org)
#y = data[:,25]
# Data preprocessing
# Exclude bad data
#z = np.log10(1e-6)
mask = y > -6.0 #-6.5 #(tri's orig)
#mask = y > 1e-6
y = y[mask]
X = X[mask]
###if i were to try and test on all data (good or bad idea?)
print(y.shape)
# Normalizing
X -= np.mean(X, axis=0) #these standard deviations need to be changed if im not doing log?
X /= np.std(X, axis=0)
mean_y = np.mean(y, axis=0)
stdv_y = np.std(y, axis=0)
y = (y-mean_y)/stdv_y
# Shuffle and divide into train and val set
mask = np.random.permutation(X.shape[0]) #(does this work with seed)
X = X[mask]
y = y[mask]
tfrac = int(0.8*y.size)
X_train = X[:tfrac]
y_train = y[:tfrac]
X_val = X[tfrac:]
y_val = y[tfrac:]
#trying to test against all of itself
#tfrac = int(1*y.size)
#X_train = X[:tfrac]
#y_train = y[:tfrac]
#X_val = X[:tfrac]
#y_val = y[:tfrac]
print('')
print('X_train shape: {}'.format(X_train.shape))
print('y_train shape: {}'.format(y_train.shape))
print('X_val shape: {}'.format(X_val.shape))
print('y_val shape: {}'.format(y_val.shape))
In [6]:
def QuakeNet(input_shape, lr=1e-3, reg=0.00, dropout=0.2):
#orig (input_shape, lr=1e-3, reg=0.00, dropout=0.0)
X_input = Input(input_shape)
X = Dense(512, kernel_regularizer=regularizers.l2(reg),
activation='relu')(X_input)
X = Dense(256, kernel_regularizer=regularizers.l2(reg),
activation='relu')(X)
X = Dense(128, kernel_regularizer=regularizers.l2(reg),
activation='relu')(X)
####X = Dense(16, kernel_regularizer=regularizers.l2(reg),
### # activation='sigmoid')(X)
X = Dense(64, kernel_regularizer=regularizers.l2(reg),
activation='relu')(X)
X = Dense(32, kernel_regularizer=regularizers.l2(reg),
activation='relu')(X)
X = Dense(16, kernel_regularizer=regularizers.l2(reg),
activation='relu')(X)
X = Dense(8, kernel_regularizer=regularizers.l2(reg),
activation='relu')(X)
# X = Dense(64, kernel_regularizer=regularizers.l2(reg),
# activation='relu')(X)
X = Dropout(rate=dropout)(X)
X = Dense(1, kernel_regularizer=regularizers.l2(reg))(X)
model = Model(inputs=X_input, outputs=X, name='QuakeNet')
model.compile(optimizer=Adam(lr=lr), loss='mse')
return model
In [7]:
input_shape = (X_train.shape[1], )
model = QuakeNet(input_shape=input_shape)
model.summary()
In [8]:
stats = model.fit(X_train, y_train, epochs=200, batch_size=32, validation_data=(X_val, y_val))
In [9]:
model.save('iwgv6.hdf5')
#from keras.model import load_model
model = load_model('iwgv6.hdf5')
In [10]:
y_pred = model.predict(X_val) # X_val could be new data too?
# Inverse-normalize
y_val = y_val*stdv_y + mean_y
y_pred = y_pred*stdv_y + mean_y
In [ ]:
In [11]:
print(y_val.shape)
y_pred = y_pred.flatten()
#print(y_pred.shape)
#fig, ax = plt.subplots()
#v_min = min(np.min(y_val), np.min(y_pred))
#v_max = max(np.max(y_val), np.max(y_pred))
#x = np.linspace(v_min, v_max, 1000)
#ax.plot(y_val, y_pred, '.')
#ax.plot(x, x)
#ax.set(xlabel='Prediction', ylabel='True')
#fig.tight_layout()
In [12]:
absy_val = abs(y_val)
absy_pred = abs(y_pred)
#taking antilog so we can see what we are used too
antiy_val = 10**y_val
antiy_pred = 10**y_pred
v_min = min(np.min(antiy_val), np.min(antiy_pred))
v_max = max(np.max(antiy_val), np.max(antiy_pred))
x = np.linspace(v_min, v_max, 1000)
fig, ax = plt.subplots()
ax.loglog(antiy_val, antiy_pred, '.')
ax.plot(x, x)
ax.set(ylabel='Predicted ground velocity [m/s]', xlabel='Actual ground velocity [m/s]',title = 'Actual vs. predicted ground velocities' )
ax.set(adjustable='box-forced', aspect='equal')
#ax.axis([yedges.min(),yedges.max(),yedges.min(),yedges.max()])
fig.tight_layout()
plt.savefig('iriswestomit1.png', dpi =300,bbox_inches='tight')
#plt.show()
In [13]:
x = np.linspace(v_min, v_max, 1000)
fig2, ax, = plt.subplots()
x_bins = np.logspace(np.log10(antiy_val.min()), np.log10(antiy_val.max()),np.sqrt(5000)) #12279
y_bins = np.logspace(np.log10(antiy_pred.min()), np.log10(antiy_pred.max()),np.sqrt(5000))
H, xedges, yedges = np.histogram2d(antiy_val, antiy_pred, bins=[x_bins,y_bins])
#ax2 = fig.add_subplot(212)
h = ax.pcolormesh(xedges, yedges, H.T)
#ax.set_aspect('equal')
#ax.set(adjustable='box-forced', aspect='equal')
#a2.imshow(img, origin='lower', extent=extent, aspect='auto')
ax.set_xscale('log')
ax.set_yscale('log')
ax.axis([yedges.min(),yedges.max(),yedges.min(),yedges.max()])
ax.set(ylabel='Predicted ground velocity [m/s]', xlabel='Actual ground velocity [m/s]',title = 'Actual vs. predicted ground velocities')
cbar = plt.colorbar(h, ax=ax)
ax.plot(x, x, c='r',linewidth=.5)
#ax.set_ylim([0, 10e-2])
#ax.set_xlim([0, 10e-2])
#ax.set_aspect('equal')
#cbar =plt.colorbar()
#cbar.ax.set_ylabel('Counts')
cbar.set_label('Counts', rotation=270,labelpad=9)
fig.tight_layout()
ax.set(adjustable='box', aspect='equal')
plt.savefig('iriswestdensityomit1.png', dpi =300,bbox_inches='tight')
plt.show()
In [ ]:
In [ ]:
In [ ]:
In [14]:
z = np.array(abs((antiy_val -antiy_pred)/antiy_val))
print(z)
print(z.shape)
print(np.min(z))
print(np.max(z))
print (np.average(z))
#x_bins = np.logspace(np.log10(antiy_val.min()), np.log10(antiy_val.max()),np.sqrt(12279))
#y_bins = np.logspace(np.log10(antiy_pred.min()), np.log10(antiy_pred.max()),np.sqrt(12279))
plt.hist(z, bins=30,range =[0,3.5], facecolor='blue', alpha=0.5)
plt.xlabel('(Predicted-Actual)/Actual Error')
plt.ylabel('Counts')
plt.title('Predicted ground velocity amount falling within error')
plt.savefig('iriswesthistomit1.png', dpi =300,bbox_inches='tight')
plt.show()
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]: