import numpy as np
import scipy.signal as si
from PIL import Image
def gauss_kern():
h1 = 15
h2 = 15
x, y = np.mgrid[0:h2, 0:h1]
x = x-h2/2
y = y-h1/2
sigma = 1.5
g = np.exp( -( x**2 + y**2 ) / (2*sigma**2) );
return g / g.sum()
def deriv(im1, im2):
g = gauss_kern()
Img_smooth = si.convolve(im1,g,mode='same')
fx,fy=np.gradient(Img_smooth)
ft = si.convolve2d(im1, 0.25 * np.ones((2,2))) + \
si.convolve2d(im2, -0.25 * np.ones((2,2)))
fx = fx[0:fx.shape[0]-1, 0:fx.shape[1]-1]
fy = fy[0:fy.shape[0]-1, 0:fy.shape[1]-1];
ft = ft[0:ft.shape[0]-1, 0:ft.shape[1]-1];
return fx, fy, ft
import matplotlib.pyplot as plt
import numpy as np
import scipy.signal as si
from PIL import Image
import deriv
import numpy.linalg as lin
def lk(im1, im2, i, j, window_size) :
fx, fy, ft = deriv.deriv(im1, im2)
halfWindow = np.floor(window_size/2)
curFx = fx[i-halfWindow-1:i+halfWindow,
j-halfWindow-1:j+halfWindow]
curFy = fy[i-halfWindow-1:i+halfWindow,
j-halfWindow-1:j+halfWindow]
curFt = ft[i-halfWindow-1:i+halfWindow,
j-halfWindow-1:j+halfWindow]
curFx = curFx.T
curFy = curFy.T
curFt = curFt.T
curFx = curFx.flatten(order='F')
curFy = curFy.flatten(order='F')
curFt = -curFt.flatten(order='F')
A = np.vstack((curFx, curFy)).T
U = np.dot(np.dot(lin.pinv(np.dot(A.T,A)),A.T),curFt)
return U[0], U[1]
from pylab import *
from matplotlib.widgets import Slider
ax = subplot(121)
subplots_adjust(left=0.1, bottom=0.25)
l1, = plot(None,None, lw=2, color='red')
axis([-1, 1, -8, 8])
title ('$(a + bi)t$', color='blue')
grid()
ax = subplot(122)
subplots_adjust(left=0.1, bottom=0.25)
l2, = plot(None,None, lw=2, color='red')
axis([-3, 3, -3, 3])
title ('$e^{(a + bi)t}$', color='blue')
grid()
axcolor = 'lightgoldenrodyellow'
axa = axes([0.15, 0.1, 0.65, 0.03], axisbg=axcolor)
axb = axes([0.15, 0.15, 0.65, 0.03], axisbg=axcolor)
slidera = Slider(axa, 'a', -1.0, 1.0, valinit=0)
sliderb = Slider(axb, 'b', -8.0, 8.0, valinit=0)
def update(val):
a = slidera.val
b = sliderb.val
t = arange(-1.0, 1.0, 0.001)
l1.set_xdata(t)
l1.set_ydata((b/a)*t)
t = arange(-3.0, 3.0, 0.001)
l2.set_xdata(exp(a*t)*cos(b*t))
l2.set_ydata(exp(a*t)*sin(b*t))
draw()
slidera.on_changed(update)
sliderb.on_changed(update)
show()
import matplotlib.pyplot as plt
import numpy as np
Img = plt.imread("twoObj.jpg")
n = Img.shape[0]
Img2 = Img.flatten(order='C')
nn = Img2.shape[0]
A = np.zeros((nn,nn))
for i in range(nn):
for j in range(nn):
A[i,j]=np.exp(-((Img2[i]-Img2[j])**2))
V,D = np.linalg.eig(A)
V = np.real(V)
a = np.real(D[0])
print a
threshold = 0 # filter
a = np.reshape(a, (n,n))
Img[a<threshold] = 255
plt.imshow(Img)
plt.show()
import rpy2.robjects as R
import rpy2.rinterface as rinterface
from rpy2.robjects.packages import importr
import rpy2.robjects.numpy2ri
from rpy2.robjects.vectors import FloatVector, StrVector
import datetime
import numpy as np
forecast = importr("forecast")
R.r.source("arima-fn.R")
x = np.loadtxt("garan.dat", skiprows=1)
n = len(x)
m = 10
terms = 4
x = np.insert(x, 10, np.nan) # put in a NaN to see what happens
n += 1
print x
R.r.assign('x',x)
R.r.assign('n',n)
R.r.assign('m',m)
R.r.assign('terms',terms)
R.r('fit = Arima(x=x, order=c(2,0,1), xreg=fourier(1:n,terms,m))')
R.r('f = forecast(fit, h=2*m, xreg=fourier(n+1:(2*m),terms,m))')
res = R.r('forecast_results(f)')
for x in res:
print x
Main thing to do as applied mathematicians it to be willing to listen to the other person's language, try to understand and translate it into something that makes sense mathematically. Engineers are very clever, when you talk to engineers in many different fields, they have an algorithm that works, they don't know why it works, it's not rigorous, but it's great. Then your job is to figure out how to generalize it and make mathematics out of it. That happened [to us] many, many times. Lots of fantastically successful algorithms came from the intuition of engineers, then they were made precise later by mathematicians.
import matplotlib.pyplot as plt
import numpy as np
import scipy.signal as signal
import scipy.ndimage as image
import time
from scipy import ndimage
def bwdist(a):
"""
this is an intermediary function, 'a' has only True, False vals,
so we convert them into 0, 1 values -- in reverse. True is 0,
False is 1, distance_transform_edt wants it that way.
"""
b = np.ones(a.shape)
b[a==True] = 0.
return ndimage.distance_transform_edt(b)
def gauss_kern():
""" Returns a normalized 2D gauss kernel array for convolutions """
h1 = 15
h2 = 15
x, y = np.mgrid[0:h2, 0:h1]
x = x-h2/2
y = y-h1/2
sigma = 1.5
g = np.exp( -( x**2 + y**2 ) / (2*sigma**2) );
return g / g.sum()
Img = plt.imread("twoObj.bmp")
Img = Img[::-1]
g = gauss_kern()
Img_smooth = signal.convolve(Img,g,mode='same')
Iy,Ix=np.gradient(Img_smooth)
absGradI=np.sqrt(Ix**2+Iy**2);
rows, cols = Img.shape
# initial function phi - level set is a square 4 pixels
# away from borders on each side, in 3D it looks like an empty
# box
c0=4
w=4
nrow, ncol=Img.shape
phi=c0*np.ones((nrow,ncol))
phi[w+1:-w-1, w+1:-w-1]=-c0
# edge-stopping function
g = 1 / (1+absGradI**2)
# gradient of edge-stopping function
gy,gx = np.gradient(g)
# gradient descent step size
dt=.4
# number of iterations after which we reinitialize the surface
num_reinit=10
phiOld=np.zeros((rows,cols))
# number of iterations after which we reinitialize the surface
iter=0
plt.ion()
while np.sum(np.sum(np.abs(phi-phiOld))) != 0:
# gradient of phi
gradPhiY, gradPhiX = np.gradient(phi)
# magnitude of gradient of phi
absGradPhi=np.sqrt(gradPhiX**2+gradPhiY**2)
# normalized gradient of phi - eliminating singularities
normGradPhiX=gradPhiX/(absGradPhi+(absGradPhi==0))
normGradPhiY=gradPhiY/(absGradPhi+(absGradPhi==0))
divYnormGradPhiX, divXnormGradPhiX=np.gradient(normGradPhiX)
divYnormGradPhiY, divXnormGradPhiY=np.gradient(normGradPhiY)
# curvature is the divergence of normalized gradient of phi
K = divXnormGradPhiX + divYnormGradPhiY
tmp1 = g * K * absGradPhi
tmp2 = g * absGradPhi
tmp3 = gx * gradPhiX + gy*gradPhiY
dPhiBydT =tmp1 + tmp2 + tmp3
phiOld=phi
# level set evolution equation
phi = phi + ( dt * dPhiBydT )
iter=iter+1
if np.mod(iter,num_reinit)==0:
# reinitialize the embedding function after num_reinit iterations
phi=np.sign(phi)
phi = (phi > 0) * (bwdist(phi <> 0))
if np.mod(iter,10)==0:
time.sleep(0.6)
plt.imshow(Img, cmap='gray')
plt.hold(True)
CS = plt.contour(phi,0, colors='r')
plt.draw()
plt.hold(False)
“Applied math” has come to mean math that could be applied someday. I use the term “very applied math” to denote that subset of applied mathematics that actually is applied to real problems.My graduate advisor used to say “Applied mathematics is not a subject classification. It's an attitude.” By that token, I'd say that very applied math is also an attitude, an interest in the grubby work required to see the math actually used and a willingness to carry it out. This involves not just math but also computing, consulting, managing, marketing, etc.
import numpy as npOnce the code is finished, the last u will have the final values. More detailed code along with the sparse matrix version is here.
import scipy.linalg
# Number of internal points
N = 200
# Calculate Spatial Step-Size
h = 1/(N+1.0)
k = h/2
x = np.linspace(0,1,N+2)
x = x[1:-1] # get rid of the '0' and '1' at each end
# Initial Conditions
u = np.transpose(np.mat(10*np.sin(np.pi*x)))
# second derivative matrix
I2 = -2*np.eye(N)
E = np.diag(np.ones((N-1)), k=1)
D2 = (I2 + E + E.T)/(h**2)
I = np.eye(N)
TFinal = 1
NumOfTimeSteps = int(TFinal/k)
for i in range(NumOfTimeSteps):
# Solve the System: (I - k/2*D2) u_new = (I + k/2*D2)*u_old
A = (I - k/2*D2)
b = np.dot((I + k/2*D2), u)
u = scipy.linalg.solve(A, b)
import numpy as np
def del2(M):
dx = 1
dy = 1
rows, cols = M.shape
dx = dx * np.ones ((1, cols - 1))
dy = dy * np.ones ((rows-1, 1))
mr, mc = M.shape
D = np.zeros ((mr, mc))
if (mr >= 3):
## x direction
## left and right boundary
D[:, 0] = (M[:, 0] - 2 * M[:, 1] + M[:, 2]) / (dx[:,0] * dx[:,1])
D[:, mc-1] = (M[:, mc - 3] - 2 * M[:, mc - 2] + M[:, mc-1]) \
/ (dx[:,mc - 3] * dx[:,mc - 2])
## interior points
tmp1 = D[:, 1:mc - 1]
tmp2 = (M[:, 2:mc] - 2 * M[:, 1:mc - 1] + M[:, 0:mc - 2])
tmp3 = np.kron (dx[:,0:mc -2] * dx[:,1:mc - 1], np.ones ((mr, 1)))
D[:, 1:mc - 1] = tmp1 + tmp2 / tmp3
if (mr >= 3):
## y direction
## top and bottom boundary
D[0, :] = D[0,:] + \
(M[0, :] - 2 * M[1, :] + M[2, :] ) / (dy[0,:] * dy[1,:])
D[mr-1, :] = D[mr-1, :] \
+ (M[mr-3,:] - 2 * M[mr-2, :] + M[mr-1, :]) \
/ (dy[mr-3,:] * dx[:,mr-2])
## interior points
tmp1 = D[1:mr-1, :]
tmp2 = (M[2:mr, :] - 2 * M[1:mr - 1, :] + M[0:mr-2, :])
tmp3 = np.kron (dy[0:mr-2,:] * dy[1:mr-1,:], np.ones ((1, mc)))
D[1:mr-1, :] = tmp1 + tmp2 / tmp3
return D / 4
import random
items = [(10, "low"),
(100, "mid"),
(890, "large")]
w = [10, 100, 890]
v = ["low", "mid", "large"]
def weighted_sample(ws, vs, n):
total = float(sum(w for w in ws))
i = 0
w = ws[0]
v = vs[0]
while n:
x = total * (1 - random.random() ** (1.0 / n))
total -= x
while x > w:
x -= w
i += 1
w = ws[i]
v = vs[i]
w -= x
yield v
n -= 1
for i in weighted_sample(w, v, 500):
print i