In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
import matplotlib
import numpy as np
import h5py
In [2]:
def generateLand(N,sizeGrid,cellLow,cellUp,cellMid,errStd,eff,shareZero=0.5,neigPat='block'):
"""
Generate a landscape with size (sizeGrid x sizeGrid) for N observations
Inputs:
neigEff : Specify the type of neighboring effect "sum", "mean","nanmean" or HHL (a measure
of the concentration in the neighborhood)
neigPat : Specify the type of neighbroing relations "block", "cross", "star"
"""
K = b1.shape[0]
# memory efficient way to generate landscape with characteristics between 0-255 with around also zero observations
land = np.random.randint(0,high=255,size=(N,sizeGrid, sizeGrid,K),dtype='uint8')
landmask = np.random.choice(a=[False, True], size=(N,sizeGrid, sizeGrid,K), p=[1-shareZero, shareZero])
land[landmask] = 0
del landmask
land[:,cellMid,cellMid,:] = 0 # Set mid point to zero
# create onw characteristics
ownX = np.random.randint(maxX,size=(N,K))
# create an error
err = np.random.normal(0, errStd, (N,1))
# Add constand, own effect and error
Y = eff['const']+np.matmul(ownX,eff['own'])+err
# --------------------------------------------------------
# Create a neigboring mask, different option are avaliable
m = np.zeros((sizeGrid,sizeGrid))
if neigPat =='block':
m[cellLow:cellUp,cellLow:cellUp] = 1
elif neigPat == 'cross':
m[cellMid,cellLow:cellUp] = 1
m[cellLow:cellUp,cellMid] = 1
elif neigPat == 'star':
m[cellMid,cellLow:cellUp] = 1
m[cellLow:cellUp,cellMid] = 1
for j in range(0,numNeigCell):
m[cellMid+j,cellMid+j] = 1
m[cellMid-j,cellMid+j] = 1
m[cellMid+j,cellMid-j] = 1
m[cellMid-j,cellMid-j] = 1
else:
raise Warning("neigPat not defined")
maskN = np.ma.make_mask(m)
# --------------------------------------------------------
# Add neighboring effects
if 'sum' in eff:
neigX = np.sum(land[:,maskN,:],axis=(1)).reshape(N,K) # reshape necessary for the case N=1
Y += np.matmul(neigX,eff['sum'])
if 'mean' in eff:
neigX = np.mean(land[:,maskN,:],axis=(1)).reshape(N,K)
Y += np.matmul(neigX,eff['mean'])
if 'nanmean' in eff:
# Calcualte mean without zero
nSum = np.sum(land[:,maskN,:],axis=(1))
nCount = np.sum(land[:,maskN,:]>0,axis=(1))
neigX = nSum/nCount
Y += np.matmul(neigX,eff['nanmean'])
if 'hhl' in eff:
# Calculate the Herfindahl index in the neighborhood
# see: https://en.wikipedia.org/wiki/Herfindahl_index
sumX = np.sum(land[:,maskN,:],axis=(1)).reshape(N,K)
# Repeat the sum for observaton N and channel K to the size of the landscape (sizeGrid x sizeGrid)
zz= sumX.reshape(N,1,1,K)
zz = np.repeat(zz, sizeGrid, axis=1)
zz = np.repeat(zz, sizeGrid, axis=2)
# Calcualte the a_i over a_bar (i.e. the mean in the maskN) to the for the entire landscape
# i.e. select the maskN in the next step
aa = (land/zz)**2
# select only the maskN to calcualte the sum
neigX = np.sum(aa[:,maskN,:],axis=(1))
Y += np.matmul(neigX,eff['hhl'])
return land, ownX,Y
In [44]:
sizeGrid = 51
maxX = 255
K =2
# set coefficients
b0 = np.random.normal(5, 2, 1)
b1 = np.random.normal(4, 1, K).reshape(K,1) # coef for own characteristics
b2 = np.random.normal(1, 0.1, K).reshape(K,1) # coef for neighbor characteristics
b3 = np.random.normal(1, 0.1, K).reshape(K,1) # coef for neighbor characteristics
b4 = np.random.normal(1, 0.1, K).reshape(K,1) # coef for neighbor characteristics
b5 = np.random.normal(1, 0.1, K).reshape(K,1) # coef for neighbor characteristics
errStd = 0 # error added to Y
numNeigCell = 15 # number of cells in the middle that are relevant for neighboring effects
cellLow = int((sizeGrid-2*numNeigCell)/2)
cellUp = cellLow+2*numNeigCell+1
cellMid = int(sizeGrid/2)
"""
eff = {'const':b0,
'own':b1}
"""
eff = {'const':b0,
'own':b1,
'sum':b2,
'mean':b3,
'nanmean':b4,
'hhl':b5
}
shareZero = 0.2
In [19]:
In [47]:
f = h5py.File("nn_spatial.hdf5", "w")
dName = {'test':2000,
'train':20000,
'val':2000
}
numI = 5 # multiplyier of the dataset
for k,v in dName.items():
N = v
dLand = f.create_dataset((k+"/land"), (N*numI,sizeGrid,sizeGrid,K), dtype='i')
dOwnX = f.create_dataset((k+"/ownX"), (N*numI,K), dtype='f')
dY = f.create_dataset((k+"/Y"), (N*numI,1), dtype='f')
for i in range(0,numI):
dLand[((numI-1)*N):((numI)*N),:,:,:], dOwnX[((numI-1)*N):((numI)*N),:],dY[((numI-1)*N):((numI)*N),:] = generateLand(N,sizeGrid,cellLow,cellUp,cellMid,errStd,eff,shareZero=shareZero,neigPat='block')
f.close()
In [46]:
f.close()
In [8]:
def printname(name):
print(name)
f.visit(printname)
In [2]:
from __future__ import print_function
import numpy as np
import keras
from keras import layers
from keras.layers import Input, Dense, Activation, ZeroPadding2D, BatchNormalization, Flatten, Conv2D
from keras.layers import AveragePooling2D, MaxPooling2D, Dropout, GlobalMaxPooling2D, GlobalAveragePooling2D
from keras.models import Model
from keras.preprocessing import image
from keras.utils import layer_utils
from keras.utils.data_utils import get_file
from keras.applications.imagenet_utils import preprocess_input
from IPython.display import SVG
from keras.utils.vis_utils import model_to_dot
from keras.utils import plot_model
import pydot
import keras.backend as K
K.set_image_data_format('channels_last')
import matplotlib.pyplot as plt
from matplotlib.pyplot import imshow
%matplotlib inline
In [42]:
#model 0
def model0(input_shape,num_classes):
model = Sequential()
model.add(Conv2D(8, kernel_size=(7, 7),
activation='relu',
input_shape=input_shape,
name='conv1'))
model.add(MaxPooling2D(pool_size=(2, 2)))
model.add(Conv2D(16, kernel_size=(7, 7),
activation='relu',
name='conv2',padding='same'))
model.add(MaxPooling2D(pool_size=(2, 2)))
model.add(Conv2D(32, 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(128, activation='relu'))
model.add(Dense(9, 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 [3]:
def model1(input_shape,K_own,num_classes):
"""
Implementation of the HappyModel.
Arguments:
input_shape -- shape of the images of the dataset
Returns:
model -- a Model() instance in Keras
"""
X_input = Input(input_shape)
# Zero-Padding: pads the border of X_input with zeroes
#X = ZeroPadding2D((3, 3))(X_input)
# CONV -> BN -> RELU Block applied to X
X = Conv2D(8, (7, 7), strides = (1, 1), name = 'conv0')(X_input)
#X = BatchNormalization(axis = 3, name = 'bn0')(X)
X = Activation('relu')(X)
# MAXPOOL
X = MaxPooling2D((2, 2), name='max_pool1')(X)
X = Conv2D(16, (7, 7), strides = (1, 1), padding='same', name = 'conv1')(X)
X = Activation('relu')(X)
X = MaxPooling2D((2, 2), name='max_pool2')(X)
X = Conv2D(32, (3, 3), strides = (1, 1), padding='valid', name = 'conv2')(X)
X = Activation('relu')(X)
# FLATTEN X (means convert it to a vector) + FULLYCONNECTED
X = Flatten()(X)
xOwn = Input([K_own])
X = keras.layers.concatenate([X,xOwn], axis=1)
X = Dense(9, activation='relu', name='fc1')(X)
X = Dense(num_classes, activation='relu', name='fc2')(X)
# Create model. This creates your Keras model instance, you'll use this instance to train/test the model.
model = Model(inputs = [X_input,xOwn], outputs = X, name='model1')
return model
In [131]:
K.tf.reset_default_graph()
X_input = Input(input_shape)
# Zero-Padding: pads the border of X_input with zeroes
#X = ZeroPadding2D((3, 3))(X_input)
# CONV -> BN -> RELU Block applied to X
X = Conv2D(8, (7, 7), strides = (1, 1), name = 'conv0')(X_input)
#X = BatchNormalization(axis = 3, name = 'bn0')(X)
X = Activation('relu')(X)
# MAXPOOL
X = MaxPooling2D((2, 2), name='max_pool1')(X)
X = Conv2D(16, (7, 7), strides = (1, 1), padding='same', name = 'conv1')(X)
X = Activation('relu')(X)
X = MaxPooling2D((2, 2), name='max_pool2')(X)
X = Conv2D(32, (3, 3), strides = (1, 1), padding='valid', name = 'conv2')(X)
X = Activation('relu')(X)
# FLATTEN X (means convert it to a vector) + FULLYCONNECTED
X = Flatten()(X)
xOwn = Input([1],name='xOwn')
X = keras.layers.concatenate([X,xOwn], axis=1)
In [132]:
X
Out[132]:
In [128]:
xOwn
Out[128]:
In [4]:
# %%
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 [5]:
f = h5py.File("nn_spatial.hdf5", "r")
In [16]:
y_train = f['train/Y']
x_train = f['train/land']
ownX_train = f['train/ownX']
y_test = f['test/Y']
x_test = f['test/land']
ownX_test = f['test/ownX']
In [17]:
input_shape = f['test/land'][0,:,:,:].shape
num_classes = 1
model = model1(input_shape,ownX_train.shape[1],num_classes)
In [18]:
model.summary()
In [92]:
from keras.utils import plot_model
plot_model(model, to_file='model.png')
In [19]:
model.compile(loss='mean_squared_error',
optimizer='Adadelta',
metrics=['mae'])
In [21]:
batch_size = 256
epochs = 1
hist= model.fit([x_train,ownX_train], y_train,
batch_size=batch_size,
epochs=epochs,
verbose=1,
validation_data=([x_test,ownX_test], y_test),
shuffle="batch")
In [ ]: