In [315]:
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
import matplotlib
import numpy as np
from scipy.spatial import distance_matrix
from scipy import sparse
What does a SLX model do
One interpretation of the SLX model is that it computes the average of all observations within all cells not more then three away from the farm
If we assume we define for each farm a 9x9 cell
In [316]:
N = 9
# make an empty data set
data = np.ones((N, N)) * np.nan
# fill in some fake data
for j in range(3)[::-1]:
data[N//2 - j : N//2 + j +1, N//2 - j : N//2 + j +1] = j
data[data==2] = 1
# make a figure + axes
fig, ax = plt.subplots(1, 1)
# make color map
my_cmap = matplotlib.colors.ListedColormap(['r', 'g', 'b'])
# set the 'bad' values (nan) to be white and transparent
my_cmap.set_bad(color='w', alpha=0)
# draw the grid
for x in range(N + 1):
ax.axhline(x, lw=2, color='k', zorder=5)
ax.axvline(x, lw=2, color='k', zorder=5)
# draw the boxes
ax.imshow(data, interpolation='none', cmap=my_cmap, extent=[0, N, 0, N], zorder=0)
# turn off the axis labels
ax.axis('off')
Out[316]:
This could also be implemented using a 10x10 convolution filter defined as follows
In [317]:
slxFilt = np.zeros((9,9))
slxFilt[2:7,2:7]=1
slxFilt[4,4]=0
slxFilt
Out[317]:
In [318]:
# =============================================================================
# Gernerate data
# =============================================================================
sizeGrid = 9 # size of the grid
nObs = 6 # number of observations
K = 1 # number of features
# set coefficients
b0 = np.random.normal(5, 2, 1)
b1 = np.random.normal(5, 2, K).reshape(K,1) # coef for own characteristics
b2 = np.random.normal(5, 2, K).reshape(K,1) # coef for neighbor characteristics
errStd = 0 # error added to Y
cutW = distance_matrix([[4,4]], [[4+2,4+2]]) # cut distance in W all cells not more then 2 away
In [319]:
def generate_landscape(sizeGrid,nObs,K,b0,b1,b2,errStd,cutW):
"""
Generate a landscape with a fixed number of observations (nObs)
Inputs
nObs = Number of farms
K = Number of characteristics X
b0 = constant
b1 = coef of own effect
b2 = coef of WX
errStd = std of error terms
cutW = disctance cut off
"""
# Create location of observations
locOwn = np.array([[int(sizeGrid/2),int(sizeGrid/2)]])
while True:
loc = np.random.randint(0,sizeGrid, size=(nObs,2))
locAll = np.concatenate((locOwn, loc), axis=0)
# Make sure that farm are not one same location
locAll = np.unique(locAll, axis=0)
if nObs+1 == locAll.shape[0]:
break
# Create own characteristics
X = np.random.randint(0,100, size=(nObs,K))
Xown = np.random.randint(0,100, size=(1,K))
# Create spatial weigthing matrix W
W = distance_matrix(locOwn, loc)<=cutW
row_sum = W.sum(axis=1,keepdims=True) # calc row sum
if row_sum!=0:
W = W/row_sum # row standardize
# Create error
err = np.random.normal(0, errStd, 1)
# Calcualte Y
Y = b0 + np.matmul(Xown,b1)+ np.matmul(np.matmul(W,X),b2)+err
assert(Y.shape==(1,1))
maps = np.zeros((sizeGrid,sizeGrid,K))
#
for k in range(0,K):
I = np.concatenate((locOwn[:,0],loc[:,0]),axis=0)
J = np.concatenate((locOwn[:,1],loc[:,1]),axis=0)
V = np.concatenate((Xown[:,k],X[:,k]),axis=0)
A = sparse.coo_matrix((V,(I,J)),shape=(sizeGrid,sizeGrid))
maps[:,:,k] = A.todense()
#
return maps,Y,X,Xown,W,loc,locOwn
In [320]:
def generate_landscape2(sizeGrid,nobsMax,K,b0,b1,b2,errStd,cutW):
"""
Generate a landscape with a random number of observations, only the
maximum number of neighboring Observations is set (nObsMax)
Inputs
nObs = Number of farms
K = Number of characteristics X
b0 = constant
b1 = coef of own effect
b2 = coef of WX
errStd = std of error terms
cutW = disctance cut off
"""
# Draw nObs
nObs = int(np.random.randint(0,nobsMax, size=(1,1)))
# Create location of observations
locOwn = np.array([[int(sizeGrid/2),int(sizeGrid/2)]])
while True:
loc = np.random.randint(0,sizeGrid, size=(nObs,2))
locAll = np.concatenate((locOwn, loc), axis=0)
# Make sure that farm are not one same location
locAll = np.unique(locAll, axis=0)
if nObs+1 == locAll.shape[0]:
break
# Create own characteristics
X = np.random.randint(0,100, size=(nObs,K))
Xown = np.random.randint(0,100, size=(1,K))
# Create spatial weigthing matrix W
W = distance_matrix(locOwn, loc)<=cutW
row_sum = W.sum(axis=1,keepdims=True) # calc row sum
if row_sum!=0:
W = W/row_sum # row standardize
# Create error
err = np.random.normal(0, errStd, 1)
# Calcualte Y
Y = b0 + np.matmul(Xown,b1)+ np.matmul(np.matmul(W,X),b2)+err
assert(Y.shape==(1,1))
maps = np.zeros((sizeGrid,sizeGrid,K))
#
for k in range(0,K):
I = np.concatenate((locOwn[:,0],loc[:,0]),axis=0)
J = np.concatenate((locOwn[:,1],loc[:,1]),axis=0)
V = np.concatenate((Xown[:,k],X[:,k]),axis=0)
A = sparse.coo_matrix((V,(I,J)),shape=(sizeGrid,sizeGrid))
maps[:,:,k] = A.todense()
#
return maps,Y,X,Xown,W,loc,locOwn
In [322]:
In [324]:
maps,y,x,xown,w,loc,locOwn = generate_landscape3(sizeGrid,nObs,K,b0,b1,b2,errStd,cutW)
print("y: \n",y)
print("x: \n",x)
print("xown: \n",xown)
print("map k=0: \n",maps[:,:,0])
#print("map k=1: \n",maps[:,:,1])
#print("map k=2: \n",maps[:,:,2])
print("loc: \n",loc)
print("w: \n",w)
print("locOwn: \n",locOwn)
In [325]:
from __future__ import print_function
import keras
from keras.datasets import mnist
from keras.models import Sequential
from keras.layers import Dense, Dropout, Flatten
from keras.layers import Conv2D, MaxPooling2D
from keras import backend as K
from scipy.spatial import distance_matrix
from scipy import sparse
from sklearn.preprocessing import MinMaxScaler
from sklearn.model_selection import train_test_split
In [9]:
def mapToX2Canal(N,nObsMax):
# %%
#X = np.zeros((N,sizeGrid,sizeGrid,K))
X = np.zeros((N,sizeGrid,sizeGrid,2))
Y = np.zeros((N,1))
for i in range(0,N):
#
maps,y,x,xown,w,loc,locOwn = generate_landscape2(sizeGrid,nObsMax,K,b0,b1,b2,errStd,cutW)
Y[i,:] = y
#X[i,:,:,:] = maps
#X[i,:,:,:] = maps
X[i,:,:,:] = maps
X[i,4,4,0] = 0
X[i,:,:,1] = np.zeros((sizeGrid,sizeGrid))
X[i,4,4,1] = maps[4,4,0]
return X,Y
In [10]:
def mapToX(N,nObsMax):
# %%
#X = np.zeros((N,sizeGrid,sizeGrid,K))
X = np.zeros((N,sizeGrid,sizeGrid,1))
Y = np.zeros((N,1))
for i in range(0,N):
#
maps,y,x,xown,w,loc,locOwn = generate_landscape2(sizeGrid,nObsMax,K,b0,b1,b2,errStd,cutW)
Y[i,:] = y
X[i,:,:,:] = maps
return X,Y
In [11]:
def standardSplit(X,Y,test_size=0.1,random_state=42):
# %% standardized features and targets
x_min_max = MinMaxScaler()
y_min_max = MinMaxScaler()
X_minmax = x_min_max.fit_transform(X.reshape(X.shape[0],-1)).reshape(X.shape)
Y_minmax = y_min_max.fit_transform(Y)
# %% Split sample in test and training_set
x_train, x_test, y_train, y_test = train_test_split(X_minmax, Y_minmax, test_size=test_size, random_state=random_state)
return x_min_max, y_min_max, x_train, x_test, y_train, y_test
In [12]:
#model 0
def model0(input_shape,num_classes):
model = Sequential()
model.add(Conv2D(32, kernel_size=(3, 3),
activation='relu',
input_shape=input_shape))
#model.add(Conv2D(64, (3, 3), activation='relu'))
#model.add(MaxPooling2D(pool_size=(2, 2)))
#model.add(Dropout(0.25))
model.add(Flatten())
#model.add(Dense(128, activation='relu'))
model.add(Dense(128, activation='relu'))
#model.add(Dropout(0.5))
model.add(Dense(num_classes, activation='relu'))
model.compile(loss='mean_squared_error',
optimizer=keras.optimizers.Adadelta(),
metrics=['mae'])
return model
In [13]:
#model 1
def model1(input_shape,num_classes):
model = Sequential()
model.add(Conv2D(32, kernel_size=(3, 3),
activation='relu',
input_shape=input_shape))
model.add(Conv2D(64, (3, 3), activation='relu'))
model.add(MaxPooling2D(pool_size=(2, 2)))
#model.add(Dropout(0.25))
model.add(Flatten())
model.add(Dense(128, activation='relu'))
model.add(Dense(128, activation='relu'))
#model.add(Dropout(0.5))
model.add(Dense(num_classes, activation='relu'))
model.compile(loss='mean_squared_error',
optimizer=keras.optimizers.Adadelta(),
metrics=['mae'])
return model
In [14]:
#model 2
def model2(input_shape,num_classes):
model = Sequential()
model.add(Conv2D(32, kernel_size=(3, 3),
activation='relu',padding='same',
input_shape=input_shape))
model.add(Conv2D(64, (3, 3), activation='relu',padding='same',))
#model.add(MaxPooling2D(pool_size=(2, 2)))
#model.add(Dropout(0.25))
model.add(Flatten())
model.add(Dense(128, activation='relu'))
model.add(Dense(128, activation='relu'))
#model.add(Dropout(0.5))
model.add(Dense(num_classes, activation='relu'))
model.compile(loss='mean_squared_error',
optimizer=keras.optimizers.Adadelta(),
metrics=['mae'])
return model
In [15]:
#model 3
def model3(input_shape,num_classes):
model = Sequential()
model.add(Conv2D(16, kernel_size=(3, 3),
activation='relu',
input_shape=input_shape, name='conv1'))
model.add(MaxPooling2D(pool_size=(2, 2)))
model.add(Conv2D(32, kernel_size=(3, 3),
activation='relu',
name='conv2'))
#model.add(Conv2D(2, kernel_size=(3, 3),
# activation='relu',
# name='conv3'))
#model.add(Conv2D(64, (3, 3), activation='relu'))
#model.add(Dropout(0.25))
model.add(Flatten())
model.add(Dense(64, activation='relu'))
#model.add(Dense(128, activation='relu'))
#model.add(Dropout(0.5))
model.add(Dense(num_classes, activation='relu'))
model.compile(loss='mean_squared_error',
optimizer=keras.optimizers.Adadelta(),
metrics=['mae'])
return model
In [238]:
#model 4
def model4(input_shape,num_classes):
model = Sequential()
model.add(Conv2D(2, kernel_size=(4, 4),
activation='relu',
input_shape=input_shape, name='conv1'))
#model.add(MaxPooling2D(pool_size=(2, 2)))
model.add(Conv2D(4, kernel_size=(3, 3),
activation='relu',
name='conv2'))
model.add(Conv2D(8, kernel_size=(2, 2),
activation='relu',
name='conv3'))
#model.add(Conv2D(2, kernel_size=(3, 3),
# activation='relu',
# name='conv3'))
#model.add(Conv2D(64, (3, 3), activation='relu'))
#model.add(Dropout(0.25))
model.add(Flatten())
model.add(Dense(12, activation='relu',name='dens1'))
#model.add(Dense(128, activation='relu'))
#model.add(Dropout(0.5))
model.add(Dense(num_classes, activation='relu',name='dens2'))
model.compile(loss='mean_squared_error',
optimizer=keras.optimizers.Adadelta(),
metrics=['mae'])
return model
In [16]:
sizeGrid = 9 # size of the grid
nObs = 6 # number of observations
K = 1 # number of features
N = 60000
nObsMax = 20
# set coefficients
b0 = np.random.normal(5, 2, 1)
b1 = np.random.normal(5, 2, K).reshape(K,1) # coef for own characteristics
b2 = np.random.normal(5, 2, K).reshape(K,1) # coef for neighbor characteristics
errStd = 0 # error added to Y
cutW = distance_matrix([[4,4]], [[4+2,4+2]]) # cut distance in W all cells not more then 2 away
# Generate data
X,Y = mapToX(N,nObsMax)
#X,Y = mapToX2Canal(N,nObsMax)
# Standardize data and split in training and test
x_min_max, y_min_max, x_train, x_test, y_train, y_test = standardSplit(X,Y)
In [17]:
# %%
def r2(y_true, y_pred):
"""Calcualte and return R2.
y_true -- the observed values
y_pred -- the prediced values
"""
SS_res = np.sum(np.square(y_true - y_pred))
SS_tot = np.sum(np.square(y_true - np.mean(y_true)))
return ( 1 - SS_res/SS_tot )
In [269]:
input_shape = X[0,:,:,:].shape
num_classes = 1
#model = model0(input_shape,num_classes)
#model = model1(input_shape,num_classes)
#model = model2(input_shape,num_classes)
model = model3(input_shape,num_classes)
#model = model4(input_shape,num_classes)
In [240]:
tensorboardCall = keras.callbacks.TensorBoard(log_dir='./nn_spatial/run',
histogram_freq=10,
write_graph=True,
write_grads=False,
write_images=False,
embeddings_freq=0,
embeddings_layer_names=None,
embeddings_metadata=None)
In [271]:
batch_size = 128
epochs = 10
hist= model.fit(x_train, y_train,
batch_size=batch_size,
epochs=epochs,
verbose=1,
validation_data=(x_test, y_test))
In [275]:
model.summary()
In [200]:
l1 = model.get_layer('conv1')
l1.get_weights()[0][:,:,0,0]
Out[200]:
In [93]:
print(hist.history)
In [35]:
score = model.evaluate(x_test, y_test, verbose=0)
print('Test loss:', score[0])
print('Test accuracy:', score[1])
In [273]:
# %% Model
Yhat_test = model.predict(x_test,batch_size=32)
oY_test = y_min_max.inverse_transform(y_test)
oY_hat = y_min_max.inverse_transform(Yhat_test)
#oY_test = y_test
#oY_hat = Yhat_test
# %%
fig, ax = plt.subplots()
ax.scatter(oY_test, oY_hat, edgecolors=(0, 0, 0))
ax.plot([oY_test.min(), oY_test.max()], [oY_test.min(), oY_test.max()], 'k--', lw=4)
ax.set_xlabel('Measured')
ax.set_ylabel('Predicted')
plt.show()
r2Model = r2(oY_test,oY_hat)
print("R2 Model: ",r2Model)
In [ ]:
In [39]:
# %% Model that only knows own effect but this one perfectly
oX_test = x_min_max.inverse_transform(x_test.reshape(x_test.shape[0],-1)).reshape(x_test.shape)
oY_test = y_min_max.inverse_transform(y_test)
oY_hat = b0+b1*oX_test[:,4,4,0]
#oY_test = y_test
#oY_hat = Yhat_test
# %%
fig, ax = plt.subplots()
ax.scatter(oY_test, oY_hat, edgecolors=(0, 0, 0))
ax.plot([oY_test.min(), oY_test.max()], [oY_test.min(), oY_test.max()], 'k--', lw=4)
ax.set_xlabel('Measured')
ax.set_ylabel('Predicted')
plt.show()
r2Model = r2(oY_test,oY_hat)
print("R2 Model: ",r2Model)
Based on https://blog.keras.io/how-convolutional-neural-networks-see-the-world.html
In [312]:
'''Visualization of the filters of VGG16, via gradient ascent in input space.
This script can run on CPU in a few minutes.
Results example: http://i.imgur.com/4nj4KjN.jpg
'''
from __future__ import print_function
from scipy.misc import imsave
import numpy as np
import time
from keras.applications import vgg16
from keras import backend as K
# dimensions of the generated pictures for each filter.
img_width = 9
img_height = 9
input_channel = 1
# the name of the layer we want to visualize
# (see model definition at keras/applications/vgg16.py)
layer_name = 'dense_15'
# util function to convert a tensor into a valid image
def deprocess_image(x):
#xs = x.reshape(1,x.shape[0],x.shape[1],x.shape[2])
#xo = x_min_max.inverse_transform(xs.reshape(xs.shape[0],-1)).reshape(xs.shape)
#xe = xo.reshape(x.shape[0],x.shape[1],x.shape[2])
# normalize tensor: center on 0., ensure std is 0.1
#x -= x.mean()
x -= 0.5
x /= (x.std() + 1e-5)
x *= 0.1
# clip to [0, 1]
x += 0.5
x = np.clip(x, 0, 1)
# convert to RGB array
x *= 255
if K.image_data_format() == 'channels_first':
x = x.transpose((1, 2, 0))
x = np.clip(x, 0, 255).astype('uint8')
return x
model.summary()
# this is the placeholder for the input images
input_img = model.input
# get the symbolic outputs of each "key" layer (we gave them unique names).
layer_dict = dict([(layer.name, layer) for layer in model.layers[1:]])
def normalize(x):
# utility function to normalize a tensor by its L2 norm
return x / (K.sqrt(K.mean(K.square(x))) + 1e-5)
kept_filters = []
for filter_index in range(1):
# we only scan through the first 200 filters,
# but there are actually 512 of them
# Hugo: Changed to 12
print('Processing filter %d' % filter_index)
start_time = time.time()
# we build a loss function that maximizes the activation
# of the nth filter of the layer considered
layer_output = layer_dict[layer_name].output
loss = K.mean(layer_output)
# we compute the gradient of the input picture wrt this loss
grads = K.gradients(loss, input_img)[0]
# normalization trick: we normalize the gradient
grads = normalize(grads)
# this function returns the loss and grads given the input picture
iterate = K.function([input_img], [loss, grads])
# step size for gradient ascent
step = 0.1
# we start from a gray image with some random noise
if K.image_data_format() == 'channels_first':
input_img_data = np.random.random((1, input_channel, img_width, img_height))
else:
input_img_data = np.random.random((1, img_width, img_height, input_channel))
#input_img_data = (input_img_data - 0.5) * 20 + 128
input_img_data = np.zeros((1, img_width, img_height, input_channel))
# we run gradient ascent for 20 steps
for i in range(20):
loss_value, grads_value = iterate([input_img_data])
input_img_data += grads_value * step
print('Current loss value:', loss_value)
if loss_value <= 0.:
# some filters get stuck to 0, we can skip them
break
# decode the resulting input image
if loss_value > 0:
img = deprocess_image(input_img_data[0])
kept_filters.append((img, loss_value))
end_time = time.time()
print('Filter %d processed in %ds' % (filter_index, end_time - start_time))
# the filters that have the highest loss are assumed to be better-looking.
# we will only keep the top 64 filters.
kept_filters.sort(key=lambda x: x[1], reverse=True)
#kept_filters = kept_filters[:n * n]
# build a black picture with enough space for
# our 8 x 8 filters of size 128 x 128, with a 5px margin in between
n = len(kept_filters)
zoomFact = 50
img_width =9*zoomFact
img_height = 9*zoomFact
margin = 5
width = img_width + margin
height = n * img_height + (n - 1) * margin
stitched_filters = np.zeros(( height,width, 3))
l = 0
# fill the picture with our saved filters
for i in range(0,n):
#kept_filters[i][0][4,4] = 0
img, loss = kept_filters[i]
#img += abs(img.min())
#img *= 255.0/img.max()
stitched_filters[(img_width + margin) * i: (img_width + margin) * i + img_width,:-margin, :] = zoom(img,[zoomFact,zoomFact,1],order=0)
l +=1
# save the result to disk
imsave('stitched_filters_%dx%d.png' % (n, n), stitched_filters)
In [314]:
img[:,:,0]
Out[314]:
In [294]:
'''Visualization of the filters of VGG16, via gradient ascent in input space.
This script can run on CPU in a few minutes.
Results example: http://i.imgur.com/4nj4KjN.jpg
'''
from __future__ import print_function
from scipy.misc import imsave
import numpy as np
import time
from keras.applications import vgg16
from keras import backend as K
# dimensions of the generated pictures for each filter.
img_width = 9
img_height = 9
input_channel = 1
# the name of the layer we want to visualize
# (see model definition at keras/applications/vgg16.py)
layer_name = 'dense_15'
# util function to convert a tensor into a valid image
def deprocess_image(x):
#xs = x.reshape(1,x.shape[0],x.shape[1],x.shape[2])
#xo = x_min_max.inverse_transform(xs.reshape(xs.shape[0],-1)).reshape(xs.shape)
#xe = xo.reshape(x.shape[0],x.shape[1],x.shape[2])
"""
# normalize tensor: center on 0., ensure std is 0.1
x -= x.mean()
x /= (x.std() + 1e-5)
x *= 0.1
# clip to [0, 1]
x += 0.5
x = np.clip(x, 0, 1)
# convert to RGB array
x *= 255
if K.image_data_format() == 'channels_first':
x = x.transpose((1, 2, 0))
#x = np.clip(x, 0, 255).astype('uint8')
"""
return x
model.summary()
# this is the placeholder for the input images
input_img = model.input
# get the symbolic outputs of each "key" layer (we gave them unique names).
layer_dict = dict([(layer.name, layer) for layer in model.layers[1:]])
def normalize(x):
# utility function to normalize a tensor by its L2 norm
return x / (K.sqrt(K.mean(K.square(x))) + 1e-5)
kept_filters = []
for filter_index in range(1):
# we only scan through the first 200 filters,
# but there are actually 512 of them
# Hugo: Changed to 12
print('Processing filter %d' % filter_index)
start_time = time.time()
# we build a loss function that maximizes the activation
# of the nth filter of the layer considered
layer_output = layer_dict[layer_name].output
loss = K.mean(layer_output)
# we compute the gradient of the input picture wrt this loss
grads = K.gradients(loss, input_img)[0]
# normalization trick: we normalize the gradient
grads = normalize(grads)
# this function returns the loss and grads given the input picture
iterate = K.function([input_img], [loss, grads])
# step size for gradient ascent
step = 1.
# we start from a gray image with some random noise
if K.image_data_format() == 'channels_first':
input_img_data = np.random.random((1, input_channel, img_width, img_height))
else:
input_img_data = np.random.random((1, img_width, img_height, input_channel))
#input_img_data = (input_img_data - 0.5) * 20 + 128
# we run gradient ascent for 20 steps
for i in range(2000):
loss_value, grads_value = iterate([input_img_data])
input_img_data += grads_value * step
print('Current loss value:', loss_value)
if loss_value <= 0.:
# some filters get stuck to 0, we can skip them
break
# decode the resulting input image
if loss_value > 0:
img = deprocess_image(input_img_data[0])
kept_filters.append((img, loss_value))
end_time = time.time()
print('Filter %d processed in %ds' % (filter_index, end_time - start_time))
# the filters that have the highest loss are assumed to be better-looking.
# we will only keep the top 64 filters.
kept_filters.sort(key=lambda x: x[1], reverse=True)
#kept_filters = kept_filters[:n * n]
# build a black picture with enough space for
# our 8 x 8 filters of size 128 x 128, with a 5px margin in between
n = len(kept_filters)
zoomFact = 50
img_width =9*zoomFact
img_height = 9*zoomFact
margin = 5
width = img_width + margin
height = n * img_height + (n - 1) * margin
stitched_filters = np.zeros(( height,width, 3))
l = 0
# fill the picture with our saved filters
for i in range(0,n):
#kept_filters[i][0][4,4] = 0
img, loss = kept_filters[i]
#img += abs(img.min())
#img *= 255.0/img.max()
stitched_filters[(img_width + margin) * i: (img_width + margin) * i + img_width,:-margin, :] = zoom(img,[zoomFact,zoomFact,1],order=0)
l +=1
# save the result to disk
imsave('stitched_filters_%dx%d.png' % (n, n), stitched_filters)
In [295]:
Out[295]:
In [292]:
l = 0
# fill the picture with our saved filters
for i in range(0,n):
#kept_filters[i][0][4,4] = 0
img, loss = kept_filters[i]
img += abs(img.min())
img *= 255.0/img.max()
stitched_filters[(img_width + margin) * i: (img_width + margin) * i + img_width,:-margin, :] = zoom(img,[zoomFact,zoomFact,1],order=0)
l +=1
# save the result to disk
imsave('stitched_filters_%dx%d.png' % (n, n), stitched_filters)
In [243]:
'''Visualization of the filters of VGG16, via gradient ascent in input space.
This script can run on CPU in a few minutes.
Results example: http://i.imgur.com/4nj4KjN.jpg
'''
from __future__ import print_function
from scipy.misc import imsave
import numpy as np
import time
from keras.applications import vgg16
from keras import backend as K
# dimensions of the generated pictures for each filter.
img_width = 9
img_height = 9
input_channel = 1
# the name of the layer we want to visualize
# (see model definition at keras/applications/vgg16.py)
layer_name = 'dens2'
# util function to convert a tensor into a valid image
def deprocess_image(x):
xs = x.reshape(1,x.shape[0],x.shape[1],x.shape[2])
xo = x_min_max.inverse_transform(xs.reshape(xs.shape[0],-1)).reshape(xs.shape)
xe = xo.reshape(x.shape[0],x.shape[1],x.shape[2])
"""
# normalize tensor: center on 0., ensure std is 0.1
x -= x.mean()
x /= (x.std() + 1e-5)
x *= 0.1
# clip to [0, 1]
x += 0.5
x = np.clip(x, 0, 1)
# convert to RGB array
x *= 255
if K.image_data_format() == 'channels_first':
x = x.transpose((1, 2, 0))
#x = np.clip(x, 0, 255).astype('uint8')
"""
return xe
model.summary()
# this is the placeholder for the input images
input_img = model.input
# get the symbolic outputs of each "key" layer (we gave them unique names).
layer_dict = dict([(layer.name, layer) for layer in model.layers[1:]])
def normalize(x):
# utility function to normalize a tensor by its L2 norm
return x / (K.sqrt(K.mean(K.square(x))) + 1e-5)
kept_filters = []
for filter_index in range(8):
# we only scan through the first 200 filters,
# but there are actually 512 of them
# Hugo: Changed to 12
print('Processing filter %d' % filter_index)
start_time = time.time()
# we build a loss function that maximizes the activation
# of the nth filter of the layer considered
layer_output = layer_dict[layer_name].output
if K.image_data_format() == 'channels_first':
loss = K.mean(layer_output[:, filter_index, :, :])
else:
loss = K.mean(layer_output[:, :, :, filter_index])
# we compute the gradient of the input picture wrt this loss
grads = K.gradients(loss, input_img)[0]
# normalization trick: we normalize the gradient
grads = normalize(grads)
# this function returns the loss and grads given the input picture
iterate = K.function([input_img], [loss, grads])
# step size for gradient ascent
step = 1.
# we start from a gray image with some random noise
if K.image_data_format() == 'channels_first':
input_img_data = np.random.random((1, input_channel, img_width, img_height))
else:
input_img_data = np.random.random((1, img_width, img_height, input_channel))
#input_img_data = (input_img_data - 0.5) * 20 + 128
# we run gradient ascent for 20 steps
for i in range(20):
loss_value, grads_value = iterate([input_img_data])
input_img_data += grads_value * step
print('Current loss value:', loss_value)
if loss_value <= 0.:
# some filters get stuck to 0, we can skip them
break
# decode the resulting input image
if loss_value > 0:
img = deprocess_image(input_img_data[0])
kept_filters.append((img, loss_value))
end_time = time.time()
print('Filter %d processed in %ds' % (filter_index, end_time - start_time))
# the filters that have the highest loss are assumed to be better-looking.
# we will only keep the top 64 filters.
kept_filters.sort(key=lambda x: x[1], reverse=True)
#kept_filters = kept_filters[:n * n]
# build a black picture with enough space for
# our 8 x 8 filters of size 128 x 128, with a 5px margin in between
n = len(kept_filters)
zoomFact = 50
img_width =9*zoomFact
img_height = 9*zoomFact
margin = 5
width = img_width + margin
height = n * img_height + (n - 1) * margin
stitched_filters = np.zeros(( height,width, 3))
l = 0
# fill the picture with our saved filters
for i in range(0,n):
img, loss = kept_filters[i]
img += abs(img.min())
img *= 255.0/img.max()
stitched_filters[(img_width + margin) * i: (img_width + margin) * i + img_width,:-margin, :] = zoom(img,[zoomFact,zoomFact,1],order=0)
l +=1
# save the result to disk
imsave('stitched_filters_%dx%d.png' % (n, n), stitched_filters)
In [205]:
img += abs(img.min())
img *= 255.0/img.max()
img
Out[205]:
In [200]:
# build a black picture with enough space for
# our 8 x 8 filters of size 128 x 128, with a 5px margin in between
n = len(kept_filters)
zoomFact = 50
img_width =9*zoomFact
img_height = 9*zoomFact
margin = 5
width = img_width + margin
height = n * img_height + (n - 1) * margin
stitched_filters = np.zeros(( height,width, 3))
l = 0
# fill the picture with our saved filters
for i in range(0,n):
img, loss = kept_filters[i]
stitched_filters[(img_width + margin) * i: (img_width + margin) * i + img_width,:-margin, :] = zoom(img,[zoomFact,zoomFact,1])
l +=1
# save the result to disk
imsave('stitched_filters_%dx%d.png' % (n, n), stitched_filters)
In [194]:
stitched_filters[(img_width + margin) * i: (img_width + margin) * i + img_width,:, :].shape
Out[194]:
In [197]:
(img_width + margin) * i + img_width
Out[197]:
In [131]:
# Pick a landscape
i = 500
marEff = np.zeros((x_test.shape[0],1))
for i in range(0,x_test.shape[0]):
xi = x_test[i:i+1,:,:,:]
yi = y_test[i,:]
yi_hat = model.predict(xi)
oyi_hat = y_min_max.inverse_transform(yi_hat.reshape(-1, 1))
oyi = y_min_max.inverse_transform(yi.reshape(-1, 1))
#print("oyi:",oyi," oyi_hat: ",oyi_hat )
#Make a change to own characteristic
oxi = x_min_max.inverse_transform(xi.reshape(xi.shape[0],-1)).reshape(xi.shape)
oxi[:,4,4,0] +=1
xi_prime = x_min_max.transform(oxi.reshape(oxi.shape[0],-1)).reshape(oxi.shape)
# Make new prediction
yi_hat_prime = model.predict(xi_prime)
oyi_hat_prime = y_min_max.inverse_transform(yi_hat_prime.reshape(-1, 1))
#print('oyi_hat_prime: ',oyi_hat_prime )
#print('del pred: ', oyi_hat_prime-oyi_hat )
#print('b1',b1)
marEff[i,:] = oyi_hat_prime-oyi_hat
In [41]:
print("average marginal effect own characteristics: \n",np.mean(marEff))
print('b1: \n', b1)
In [44]:
nSize = 80
ownSize = 30
res = np.zeros((sizeGrid,sizeGrid))
for k in range(0,sizeGrid):
for l in range(0,sizeGrid):
imap = np.zeros((1,sizeGrid,sizeGrid,1))
imap[0,4,4,0] = ownSize
imap[0,k,l,0] = nSize
#imap[0,4,4,0] = 0
imap_prime = x_min_max.transform(imap.reshape(imap.shape[0],-1)).reshape(imap.shape)
# Make new prediction
yi_hat_prime = model.predict(imap_prime)
oyi_hat_prime = y_min_max.inverse_transform(yi_hat_prime.reshape(-1, 1))
res[k,l] = oyi_hat_prime
# calculate effect of neighbor by substracting own effect (b1*ownSize) and devide by neighboring size nSize
print((res-b1*ownSize)/nSize)
print(b2)
In [48]:
res = np.zeros((sizeGrid,sizeGrid))
fig = plt.subplots(9, 9, sharex=True, sharey=True)
i = 0
for k in range(0,sizeGrid):
for l in range(0,sizeGrid):
i +=1
seli = (x_test[:,k,l,0]>0)
imap = np.copy(x_test[seli,:,:,:])
iY = y_test[seli,:]
oimap = x_min_max.inverse_transform(imap.reshape(imap.shape[0],-1)).reshape(imap.shape)
oiY = y_min_max.inverse_transform(iY)
yi_hat = model.predict(imap)
oimap_prime = np.copy(oimap)
oimap_prime[:,k,l,:] += 10
imap_prime = x_min_max.transform(oimap_prime.reshape(oimap_prime.shape[0],-1)).reshape(oimap_prime.shape)
yi_hat_prime = model.predict(imap_prime)
oyi_hat = y_min_max.inverse_transform(yi_hat.reshape(-1, 1))
oyi_hat_prime = y_min_max.inverse_transform(yi_hat_prime.reshape(-1, 1))
#res[k,l] = np.mean(oyi_hat-oyi_hat_prime)
res[k,l] = np.mean(oyi_hat_prime-oyi_hat)
# Plot hist
if k==0 & l==0:
ax = plt.subplot(9, 9, i)
ax0 = ax
else:
ax = plt.subplot(9, 9, i, sharex=ax0)
ax.hist(oyi_hat_prime-oyi_hat, bins=100)
#ax.axvline(x=0,color='black')
for label in ax.get_xticklabels() + ax.get_yticklabels():
label.set_visible(False)
In [26]:
%matplotlib notebook
# use gridspec to partition the figure into subplots
import matplotlib.gridspec as gridspec
In [35]:
plt.figure()
gspec = gridspec.GridSpec(9, 9)
for k in range(0,sizeGrid):
for l in range(0,sizeGrid):
ax = plt.subplot(gspec[k, l])
ax.hist(oyi_hat_prime-oyi_hat, bins=100)
for label in ax.get_xticklabels() + ax.get_yticklabels():
label.set_visible(False)
In [36]:
oyi_hat_prime-oyi_hat
Out[36]:
In [49]:
np.set_printoptions(precision=1,suppress=True)
print(res)
print("b1: ",str(b1))
print("b2: ",str(b2))
In [118]:
imap.shape
Out[118]:
In [117]:
imap_prime.shape
Out[117]:
In [92]:
np.mean(oyi_hat_prime-oyi_hat)
Out[92]:
In [ ]:
In [ ]:
In [ ]:
In [123]:
xx = np.pad(imap,pad_width=1,mode='constant')
In [136]:
nSize = 80
ownSize = 30
res = np.zeros((sizeGrid,sizeGrid))
for k in range(0,sizeGrid):
for l in range(0,sizeGrid):
imap = np.zeros((1,sizeGrid,sizeGrid,2))
imap[0,4,4,1] = ownSize
imap[0,k-1:k+2,l-1:l+2,0] = nSize
imap[0,4,4,0] = 0
#print(imap[0,:,:,0])
imap_prime = x_min_max.transform(imap.reshape(imap.shape[0],-1)).reshape(imap.shape)
# Make new prediction
yi_hat_prime = model.predict(imap_prime)
oyi_hat_prime = y_min_max.inverse_transform(yi_hat_prime.reshape(-1, 1))
res[k,l] = oyi_hat_prime
# calculate effect of neighbor by substracting own effect (b1*ownSize) and devide by neighboring size nSize
print((res-b1*ownSize)/nSize)
print(b2)
In [126]:
imap[0,:,:,0]
Out[126]:
In [98]:
# Pick a landscape
pick_i = np.random.randint(0,x_test.shape[0], size=(200,1))
sensAvg = np.zeros((pick_i.shape[0]*sizeGrid*sizeGrid,3))
rr = 0
marEff = np.zeros((pick_i.shape[0],sizeGrid,sizeGrid))
#marEff = np.zeros((10,sizeGrid,sizeGrid))
for ii in range(0,pick_i.shape[0]):
i = int(pick_i[ii])
xi = x_test[i:i+1,:,:,:]
yi = y_test[i,:]
yi_hat = model.predict(xi)
oyi_hat = y_min_max.inverse_transform(yi_hat.reshape(-1, 1))
oyi = y_min_max.inverse_transform(yi.reshape(-1, 1))
for k in range(0,sizeGrid):
for l in range(0,sizeGrid):
#Make a change to neigh. characteristic
oxi = x_min_max.inverse_transform(xi.reshape(xi.shape[0],-1)).reshape(xi.shape)
#Save neighboring average before
xx = np.copy(oxi[:,2:7,2:7,0])
xx[xx==0] = np.nan
sensAvg[rr,0] = np.nanmean(xx)
numNeigh = np.count_nonzero(~np.isnan(xx))
# Check if k and l are withing cutoff
if (k>=2) & (k<7) & (l>=2) & (l<7):
oldMean = np.nanmean(xx)
if np.isnan(oldMean):
oldMean=0
oldSum = np.nansum(xx)
if np.isnan(oldSum):
oldSum=0
if oxi[:,k,l,0]>0:
addPart = (oldMean+1)*(numNeigh)-oldSum
else:
addPart = (oldMean+1)*(numNeigh+1)-oldSum
else:
addPart = 1
if np.isnan(addPart):
addPart = 1
# Add old average +2 such that average is increased by 1
oxi[:,k,l,0] += addPart
#Save neighboring average after
xx = np.copy(oxi[:,2:7,2:7,0])
xx[xx==0] = np.nan
sensAvg[rr,1] = np.nanmean(xx)
xi_prime = x_min_max.transform(oxi.reshape(oxi.shape[0],-1)).reshape(oxi.shape)
# Make new prediction
yi_hat_prime = model.predict(xi_prime)
oyi_hat_prime = y_min_max.inverse_transform(yi_hat_prime.reshape(-1, 1))
marEff[ii,k,l] = oyi_hat_prime-oyi_hat
sensAvg[rr,2] = oyi_hat_prime-oyi_hat
rr +=1
In [99]:
np.set_printoptions(precision=1,suppress=True)
meanMarEff = np.mean(marEff,axis=0)
meanMarEff[4,4] = 0
print('Mean marginal effects of neighboring characteristics \n',meanMarEff)
print('b2 \n',b2)
In [537]:
np.set_printoptions(precision=3)
aa = sensAvg[sensAvg[:,0]!=sensAvg[:,1]]
aa = aa[~np.isnan(np.sum(aa,axis=1))]
plt.scatter((aa[:,1]-aa[:,0]),aa[:,2]/(aa[:,1]-aa[:,0]))
print('mean',np.mean(aa[:,2]/(aa[:,1]-aa[:,0])))
print('median',np.median(aa[:,2]/(aa[:,1]-aa[:,0])))
print('b2',b2)
In [538]:
np.set_printoptions(precision=3)
aa = sensAvg[sensAvg[:,0]!=sensAvg[:,1]]
aa = aa[~np.isnan(np.sum(aa,axis=1))]
aa = aa[(aa[:,1]-aa[:,0])>0,:]
plt.scatter((aa[:,1]-aa[:,0]),aa[:,2]/(aa[:,1]-aa[:,0]))
print('mean',np.mean(aa[:,2]/(aa[:,1]-aa[:,0])))
print('median',np.median(aa[:,2]/(aa[:,1]-aa[:,0])))
print('b2',b2)
In [494]:
ll = 25
k = 6
l = 2
xx = np.copy(oxi[:,2:7,2:7,0])
xx[xx==0] = np.nan
print(xx)
#print(np.nanmean(xx))
print(sensAvg[rr-ll,0])
print(sensAvg[rr-ll,1])
print((sensAvg[rr-ll,1]-sensAvg[rr-ll,0]))
print(marEff[ii,k,l])
print(b2)
print(marEff[ii,k,l]/(sensAvg[rr-ll,1]-sensAvg[rr-ll,0]))
In [266]:
k = 6
l = 1
oxi = x_min_max.inverse_transform(xi.reshape(xi.shape[0],-1)).reshape(xi.shape)
oxi[:,k,l,0] +=1
xi_prime = x_min_max.transform(oxi.reshape(oxi.shape[0],-1)).reshape(oxi.shape)
# Make new prediction
yi_hat_prime = model.predict(xi_prime)
oyi_hat_prime = y_min_max.inverse_transform(yi_hat_prime.reshape(-1, 1))
marEff[ii,k,l] = oyi_hat_prime-oyi_hat
In [294]:
xx = oxi[:,2:7,2:7]
xx[xx==0] = np.nan
Out[294]:
In [268]:
x_min_max.inverse_transform(xi.reshape(xi.shape[0],-1)).reshape(xi.shape)[:,:,:,0]
Out[268]:
In [269]:
print(oyi_hat_prime)
print(oyi_hat)
print(oyi_hat_prime-oyi_hat)
print(b2)
In [253]:
((39+1)/2)*b2
Out[253]:
In [255]:
39*b2
Out[255]:
In [272]:
oxi[:,2:7,2:7].shape
Out[272]:
In [45]:
model.summary()
In [502]:
#import pydot
from IPython.display import SVG
from keras.utils.vis_utils import model_to_dot
from keras.utils import plot_model
plot_model(model, to_file='SLX_CNN_Model.png')
#SVG(model_to_dot(model).create(prog='dot', format='svg'))
In [ ]: