%matplotlib inline
%load_ext autoreload
%autoreload 2
import numpy as np
import matplotlib.pyplot as plt
from pymks import MKSRegressionModel
from pymks import FastMKSRegressionModel
from pymks import FiPyCHModel
from sklearn import metrics
mse = metrics.mean_squared_error
The autoreload extension is already loaded. To reload it, use: %reload_ext autoreload
from IPython.parallel import Client
rc = Client()
dview = rc[:]
#serial_result = map(lambda x:x**10, range(32))
#parallel_result = dview.map_sync(lambda x: x**10, range(32))
np.random.seed(101)
X = [np.random.random((100, 11, 11)) for x in range(32)]
#fipymodel = FiPyCHModel(dx=dx, dy=dx, dt=dt, a=1.)
def fipymodel(X):
from pymks import FiPyCHModel
dx = 0.25
dt = 1e-3
fipymodel = FiPyCHModel(dx=dx, dy=dx, dt=dt, a=1.)
return fipymodel.predict(X)
dview['fipymodel'] = fipymodel
#y = fipymodel.predict(X)
#y = dview.map_sync(lambda x: x, X)
y = dview.map_sync(lambda x: fipymodel(x), X)
y = np.concatenate(y, axis=0)
X = np.concatenate(X, axis=0)
print y.shape
(3200, 11, 11)
model = MKSRegressionModel(Nbin=5)
model.fit(X, y)
coeff = np.fft.ifftn(model.Fcoeff, axes=(0, 1)).real
Nroll = coeff.shape[0] / 2
coeff_rolled = np.roll(np.roll(coeff[:,:,0], Nroll, axis=0), Nroll, axis=1)
plt.contourf(coeff_rolled, 250)
plt.colorbar()
<matplotlib.colorbar.Colorbar instance at 0x664bf80>
from sklearn.cross_validation import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.5, random_state=3)
from sklearn.grid_search import GridSearchCV
tuned_parameters = [{'Nbin': np.arange(2, 20)}]
def neg_mse(a, b):
return -mse(a, b)
gridSearch = GridSearchCV(MKSRegressionModel(Nbin=10), tuned_parameters, cv=5, score_func=neg_mse)
gridSearch.fit(X_train, y_train)
print(gridSearch.best_estimator_)
for params, mean_score, scores in gridSearch.grid_scores_:
print("%0.5f (+/-%0.5f) for %r"% (mean_score, scores.std() / 2, params))
y_true, y_pred = y_test, gridSearch.predict(X_test)
print(mse(y_true, y_pred))
MKSRegressionModel(Nbin=6) -0.00410 (+/-0.00009) for {'Nbin': 2} -0.00411 (+/-0.00009) for {'Nbin': 3} -0.00400 (+/-0.00008) for {'Nbin': 4} -0.00399 (+/-0.00008) for {'Nbin': 5} -0.00399 (+/-0.00008) for {'Nbin': 6} -0.00399 (+/-0.00008) for {'Nbin': 7} -0.00399 (+/-0.00008) for {'Nbin': 8} -0.00400 (+/-0.00008) for {'Nbin': 9} -0.00400 (+/-0.00008) for {'Nbin': 10} -0.00400 (+/-0.00008) for {'Nbin': 11} -0.00401 (+/-0.00008) for {'Nbin': 12} -0.00401 (+/-0.00008) for {'Nbin': 13} -0.00401 (+/-0.00008) for {'Nbin': 14} -0.00402 (+/-0.00008) for {'Nbin': 15} -0.00402 (+/-0.00008) for {'Nbin': 16} -0.00402 (+/-0.00008) for {'Nbin': 17} -0.00402 (+/-0.00008) for {'Nbin': 18} -0.00403 (+/-0.00008) for {'Nbin': 19} 0.00392022262831
np.random.seed(11)
X0 = np.random.random((1, 21, 21))
dt = 1e-2
dx = 0.25
fipymodel = FiPyCHModel(dx=dx, dy=dx, dt=dt, a=1.)
#a = fipymodel.predict(X0.copy())
np.maximum(np.array((0,1,2)), 1)
array([1, 1, 2])
import numexpr as ne
X1 = X0.copy()
data = []
for step in range(701):
y = fipymodel.predict(X1)
X1 = ne.evaluate('X1 + dt * y')
X1 = np.minimum(X1, 1)
X1 = np.maximum(X1, 0)
if step % 10 == 0:
data.append(X1[0])
print step,
0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200 210 220 230 240 250 260 270 280 290 300 310 320 330 340 350 360 370 380 390 400 410 420 430 440 450 460 470 480 490 500 510 520 530 540 550 560 570 580 590 600 610 620 630 640 650 660 670 680 690 700
# JSAnimation import available at https://github.com/jakevdp/JSAnimation
from JSAnimation import IPython_display
from matplotlib import animation
# create a simple animation
fig = plt.figure()
ax = plt.axes()
im = plt.imshow(X0[0], vmin=0.3, vmax=0.7)
plt.colorbar()
def init():
im.set_array(X0[0])
return im
def animate(i):
im.set_array(data[i])
return im
anim = animation.FuncAnimation(fig, animate, init_func=init,
frames=len(data), interval=1, blit=True)
from IPython.display import HTML
plt.close(anim._fig)
HTML(IPython_display.anim_to_html(anim))
2**8
256
from pymks import FastMKSRegressionModel
import numexpr
np.random.seed(11)
X0 = np.random.random((1, 256, 256))
dt = 1e-3
model = FastMKSRegressionModel(Nbin=6, threads=8)
print X.shape
print y.shape
model.fit(X, y)
model.resize_coeff(X0.shape[1:])
data = []
X1 = X0.copy()
for step in range(15001):
y1 = model.predict(X1)
X1 = numexpr.evaluate('X1 + dt * y1')
X1 = np.minimum(X1, 1)
X1 = np.maximum(X1, 0)
if step % 100 == 0:
data.append(X1[0])
print step,
(3200, 11, 11) (3200, 11, 11) 0 100 200 300 400 500 600 700 800 900 1000 1100 1200 1300 1400 1500 1600 1700 1800 1900 2000 2100 2200 2300 2400 2500 2600 2700 2800 2900 3000 3100 3200 3300 3400 3500 3600 3700 3800 3900 4000 4100 4200 4300 4400 4500 4600 4700 4800 4900 5000 5100 5200 5300 5400 5500 5600 5700 5800 5900 6000 6100 6200 6300 6400 6500 6600 6700 6800 6900 7000 7100 7200 7300 7400 7500 7600 7700 7800 7900 8000 8100 8200 8300 8400 8500 8600 8700 8800 8900 9000 9100 9200 9300 9400 9500 9600 9700 9800 9900 10000 10100 10200 10300 10400 10500 10600 10700 10800 10900 11000 11100 11200 11300 11400 11500 11600 11700 11800 11900 12000 12100 12200 12300 12400 12500 12600 12700 12800 12900 13000 13100 13200 13300 13400 13500 13600 13700 13800 13900 14000 14100 14200 14300 14400 14500 14600 14700 14800 14900 15000
print len(data)
print data[0][0,0]
print data[1][0,0]
print data[150][0,0]
151 0.444784411698 0.48151012531 0.51896227618
for d in data[::10]:
print max(d.flatten())
print min(d.flatten())
0.881927544437 0.114925168086 0.613731183259 0.384969616014 0.669862815132 0.33748519272 0.772712312324 0.251608785198 0.952949536614 0.113355627358 1.0 0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0 0.0
from IPython.display import clear_output
import time
fig = plt.figure()
ax = plt.axes()
a = plt.imshow(X0[0], vmin=0.3, vmax=0.7)
plt.colorbar()
for d in data:
a.set_array(d)
plt.show()
time.sleep(1)
clear_output()
display(fig)
# JSAnimation import available at https://github.com/jakevdp/JSAnimation
from JSAnimation import IPython_display
from matplotlib import animation
# create a simple animation
fig = plt.figure()
ax = plt.axes()
im = plt.imshow(X0[0], vmin=0.0, vmax=1.0)
plt.colorbar()
def init():
im.set_array(X0[0])
return im
def animate(i):
im.set_array(data[i])
return im
anim = animation.FuncAnimation(fig, animate, init_func=init,
frames=len(data), interval=1, blit=True)
from IPython.display import HTML
plt.close(anim._fig)
HTML(IPython_display.anim_to_html(anim))