I recently started developing a BuddyPress Theme/Plugin using the WordPress 3.4.2 and BuddyPress 1.6.1.
To remove the Admin Bar up at the top, you can create a file called bp-custom.php which gets thrown into wp-content/plugins. The nice thing about bp-custom.php is that you don’t need to activate it or add any special code to make it run as long as BuddyPress is activated.
So to remove the Admin Bar create bp-custom.php and add the following code:
<?php
function remove_bp_adminbar(){
show_admin_bar(false);
}
add_action('after_setup_theme','remove_bp_adminbar');
?>
After updating and refreshing the page, your admin bar should be removed. But, be careful because this can cause issues if you don’t have alternative links to log out in place.
Cosine similarity measures the angle between two vectors and can be used to perform similarity between text strings. In the following code, the two input strings are vectorized and the similarity is returned as a floating point value between 0 and 1.
SET ANSI_NULLS ON
GO
SET QUOTED_IDENTIFIER ON
GO
-- =============================================
-- Author: Adam Maus (http://adammaus.com)
-- Create date: 2012-06-20
-- Description: Cosine Similarity: http://en.wikipedia.org/wiki/Cosine_similarity
-- Determines the angle between two vectors
-- =============================================
CREATE FUNCTION [dbo].[cosine_distance]
(
@s nvarchar(4000), @t nvarchar(4000)
)
RETURNS float
AS
BEGIN
-- Create an array of letter frequencies using the unicode of @s and @t
DECLARE @sLet table(letter int, freq int) -- Pretend these are vectors
DECLARE @tLet table(letter int, freq int)
DECLARE @i int, @j int
DECLARE @c int -- The current character
-- Create the arrays
SET @i = 1
WHILE @i <= LEN(@s)
BEGIN
-- Use the UNICODE values for the character but not necessary
SET @c = UNICODE(SUBSTRING(@s, @i, 1)) + 1
-- Determine whether this character is in the vector already
SELECT @j = COUNT(letter) FROM @sLet WHERE letter = @c
IF @j > 0 BEGIN
UPDATE @sLet SET freq += 1 WHERE letter = @c
END ELSE BEGIN
INSERT INTO @sLet VALUES (@c, 1)
END
-- We want to keep @tLet consistent with @sLet so insert a 0 for that letter
SELECT @j = COUNT(letter) FROM @tLet WHERE letter = @c
IF @j = 0 BEGIN
INSERT INTO @tLet VALUES (@c, 0)
END
SELECT @i += 1
END
SET @i = 1
WHILE @i <= LEN(@t)
BEGIN
SET @c = UNICODE(SUBSTRING(@t, @i, 1)) + 1
SELECT @j = COUNT(letter) FROM @tLet WHERE letter = @c
IF @j > 0 BEGIN
UPDATE @tLet SET freq += 1 WHERE letter = @c
END ELSE BEGIN
INSERT INTO @tLet VALUES (@c, 1)
END
-- We want to keep @sLet consistent with @tLet so insert a 0 for that letter
SELECT @j = COUNT(letter) FROM @sLet WHERE letter = @c
IF @j = 0 BEGIN
INSERT INTO @sLet VALUES (@c, 0)
END
SELECT @i += 1
END
-- Compute the similarity
-- Declare the numerator for the similarity
DECLARE @numer float
SET @numer = 0
SELECT @numer += s.freq * t.freq FROM @sLet s LEFT JOIN @tLet t ON s.letter = t.letter
-- Declare the norm values and calculate the denominator for the similarity
DECLARE @sNorm int, @tNorm int
SET @sNorm = 0
SET @tNorm = 0
SELECT @sNorm += freq * freq FROM @sLet
SELECT @tNorm += freq * freq FROM @tLet
DECLARE @denom float
SET @denom = SQRT(@sNorm) * SQRT(@tNorm)
RETURN (@numer) / (@denom+1) -- The +1 eliminates the possibility 0 = @denom
END
To “install” this script using Microsoft SQL Server Management Studio, go to your database, and open Programmability > Functions and right-click on Scalar-valued Functions to add a new function.
Multidimensional Scaling (MDS) is a linear embedding method used when we only know the pairwise distances between data points. For linear systems, MDS works well with as little as 10 points and the system is 2 dimensional. For this program, you choose the window size or embedding for the system and MDS identifies what the true embedding should be. The following program uses Logistic map data when r = 3.2 but you can change r to see how well MDS works as the map bifurcates and eventually becomes chaotic and nonlinear at 4. MDS breaks down at this point and the embedding defaults to just 1 over the embedding you choose.
# Requires Numpy and Scipy libraries
from numpy import *
import scipy.linalg
import random
r=3.2 # parameter on Logistic map
n = 10 # the number of timeseries points to use
d = 40 # your embedding choice
# Create a time series
def timeseries(number,throw, d):
xList = []
dataX0 = []
t = 0
while t < d:
xList.append(.1)
t += 1
while t <= number+throw+d:
x = r * xList[t-1] * (1 - xList[t-1]) # Logistic Map
xList.append(x)
if t > throw:
y = xList[t-1]
dataX0.append(x)
t = t + 1
return dataX0
# Construct an n x n centering matrix
# The form is P = I - (1/n) U where U is a matrix of all ones
def centering_matrix(n):
P = eye(n) - 1/float(n) * ones((n,n))
return P
def create_data(number, d):
dataX0 = timeseries(number, 500, d)
# Create a sliding window of size d
data = []
i = d
while i < len(dataX0):
data.append([])
j = i - d
while j < i:
data[len(data)-1].append(dataX0[j])
j = j + 1
i = i + 1
return data
def main():
print "n =",n, ": number of timeseries points used"
print "d =",d, ": embedding"
P = centering_matrix(n)
# Create the data with the number of dimensions
data = create_data(n, d)
X = pairwise_distances(data)
A = -1/2.0 * P * X * P
# Calculate the eigenvalues/vectors
[vals, vectors] = scipy.linalg.eig(A)
# Sort the values
vals = sort(vals)
embedding = 0
for x in vals:
if x > 10**(-10):
embedding += 1
print "mds embedding =",embedding
# Compute the pairwise distance between vector x and y
def metric(x1, y1):
d = 2
summ = []
i = 0
while i < len(x1):
# in this case use euclidean distance
summ.append((x1[i] - y1[i])**d)
i = i + 1
return sum(summ) ** (1 / float(d))
# Return a matrix of pairwise distances
def pairwise_distances(data):
distances = []
i = 0
while i < len(data):
distances.append([])
x1 = data[i]
j = 0
while j < len(data):
y1 = data[j]
distances[i].append(metric(x1, y1)**2)
j = j + 1
i = i + 1
return distances
main()
Step-by-step instructions for rewiring where your RSS link icon points to within the Platform theme (by PageLines) for WordPress:
Requirements:
1) WordPress
2) Platform theme by PageLines
Instructions
1) Navigate to the directory where the theme is stored on your server, usually wp-content > themes > platform
2) Download functions.php
3) Add the following line of code to the end of the file:
4) Upload functions.php and make sure that you have enabled the RSS link in the PageLines Settings under Header and Nav.
5) As of 2012-06-21, for the latest version of Platform, download includes > core.init.php and add the filter code at the very bottom.
Wolf et al. (1985) outlined an algorithm that estimates the Lyapunov spectra of systems whose equations are known using local Jacobian matrices and Gram-Schmidt orthonormalization. Python code is available for Wolf’s algorithm and discrete maps and their inverted counterparts. I have adapted this code to estimate Lyapunov spectra for continuous-time systems like the Lorenz attractor and Rossler attractor. Additionally, Python code is available to generate time series for ordinary differential equations. Lyapunov spectrum code is also available on Clint Sprott’s website.
Lorenz Attractor
Source Code:
import math, operator, random
h = 0.01
cmax = 1000 # number of iterations to perform
choice = input("Which system would you like?\n (1) Rossler \n (2) Lorenz\n")
while str(choice) != "1" and str(choice) != "2":
print("Please select 1 or 2")
choice = input("Which system would you like?\n (1) Rossler \n (2) Lorenz\n")
print "\n",
a = 0.2
b = 0.2
c = 5.7
def derivs(x, xnew, n):
if choice == 1:
return Rossler(x, xnew, n)
else:
return Lorenz(x, xnew, n)
def Rossler(x, xnew, n):
# Nonlinear Rossler Equations
xnew[1] = -x[2]-x[3]
xnew[2] = x[1] + a * x[2]
xnew[3] = b + x[3] * (x[1] - c)
# Linearized Rossler Equations
xnew[4] = -1*x[7]-x[10]
xnew[5] = -1*x[8]-x[11]
xnew[6] = -1*x[9]-x[12]
xnew[7] = x[4] + a*x[7]
xnew[8] = x[5] + a*x[8]
xnew[9] = x[6] + a*x[9]
xnew[10] = x[3]*x[4] + x[1]*x[10] - c*x[10]
xnew[11] = x[3]*x[5] + x[1]*x[11] - c*x[11]
xnew[12] = x[3]*x[6] + x[1]*x[12] - c*x[12]
return [x, xnew]
def Lorenz(x, xnew, n):
# Nonlinear Lorenz Equations
xnew[1] = 10 * (x[2] - x[1])
xnew[2] = -1*x[1] * x[3] + 28 * x[1] - x[2]
xnew[3] = x[1] * x[2] - 8/3.0 * x[3]
# Linearized Lorenz Equations
xnew[4] = -10 * x[4] + 10 * x[7]
xnew[5] = -10 * x[5] + 10 * x[8]
xnew[6] = -10 * x[6] + 10 * x[9]
xnew[7] = 28*x[4]-x[3]*x[4] - x[7] - x[1]*x[10]
xnew[8] = 28*x[5]-x[3]*x[5] - x[8] - x[1]*x[11]
xnew[9] = 28*x[6]-x[3]*x[6] - x[9] - x[1]*x[12]
xnew[10] = x[2]*x[4] + x[1]*x[7] - 8/3.0 * x[10]
xnew[11] = x[2]*x[5] + x[1]*x[8] - 8/3.0 * x[11]
xnew[12] = x[2]*x[6] + x[1]*x[9] - 8/3.0 * x[12]
return [x, xnew]
def timeseries(cmax):
X0 = []
Y0 = []
Z0 = []
xList = []
yList = []
zList = []
changeInTime = h
# Initial conditions
if choice == 1:
# Rossler
X0.append(0.01)
Y0.append(0.01)
Z0.append(0.01)
else:
# Lorenz
X0.append(0)
Y0.append(1)
Z0.append(0)
t = 0
while len(xList) <= cmax:
[x,y,z] = Rk4o(X0, Y0, Z0, h, len(X0))
X0.append(x)
Y0.append(y)
Z0.append(z)
if 200 < t:
xList.append(x)
yList.append(y)
zList.append(z)
changeInTime += h
t = t + 1
return [xList, yList, zList]
def f(x,y,z):
if choice == 1:
dxdt = -y-z
else:
dxdt = 10 * (y - x)
return dxdt
def g(x,y,z):
if choice == 1:
dydt = x + a * y
else:
dydt = 28 * x - y - x*z
return dydt
def e(x,y,z):
if choice == 1:
dzdt = b + z * (x - c)
else:
dzdt = x * y - 8/3.0 * z
return dzdt
def Rk4o(xList, yList, zList, h, t):
k1x = h*f(xList[t-1],yList[t-1], zList[t-1])
k1y = h*g(xList[t-1],yList[t-1], zList[t-1])
k1z = h*e(xList[t-1],yList[t-1], zList[t-1])
k2x = h*f(xList[t-1] + k1x/2,yList[t-1] + k1y/2, zList[t-1] + k1y/2)
k2y = h*g(xList[t-1] + k1x/2,yList[t-1] + k1y/2, zList[t-1] + k1y/2)
k2z = h*e(xList[t-1] + k1x/2,yList[t-1] + k1y/2, zList[t-1] + k1y/2)
k3x = h*f(xList[t-1] + k2x/2,yList[t-1] + k2y/2, zList[t-1] + k2y/2)
k3y = h*g(xList[t-1] + k2x/2,yList[t-1] + k2y/2, zList[t-1] + k2y/2)
k3z = h*e(xList[t-1] + k2x/2,yList[t-1] + k2y/2, zList[t-1] + k2y/2)
k4x = h*f(xList[t-1] + k3x/2,yList[t-1] + k3y/2, zList[t-1] + k3y/2)
k4y = h*g(xList[t-1] + k3x/2,yList[t-1] + k3y/2, zList[t-1] + k3y/2)
k4z = h*e(xList[t-1] + k3x/2,yList[t-1] + k3y/2, zList[t-1] + k3y/2)
x = xList[t-1] + k1x/6 + k2x/3 + k3x/3 + k4x/6
y = yList[t-1] + k1y/6 + k2y/3 + k3y/3 + k4y/6
z = zList[t-1] + k1z/6 + k2z/3 + k3z/3 + k4z/6
return [x,y,z]
n = 3 # number of variables in nonlinear system
nn=n*(n+1) # total number of variables (nonlinear + linear)
m = 0
x = []
xnew = []
v = []
ltot = []
znorm = []
gsc = []
A = []
B = []
C = []
D = []
i = 0
while i <= nn:
x.append(0)
xnew.append(0)
v.append(0)
A.append(0)
B.append(0)
C.append(0)
D.append(0)
i = i + 1
i = 0
while i <= n:
ltot.append(0)
znorm.append(0)
gsc.append(0)
i = i + 1
irate=10 # integration steps per reorthonormalization
io= 100 # number of iterations between normalization
# initial conditions for nonlinear maps
# must be within the basin of attraction
# Generate a random transient before starting the initial conditions
i = 1
while i <= n:
v[i] = 0.001
i = i + 1
transient = random.randint(n,100000)
# Generate the initial conditions for the system
[tempx,tempy,tempz] = timeseries(transient)
v[1] = tempx[len(tempx)-1]
v[2] = tempy[len(tempy)-1]
v[3] = tempz[len(tempz)-1]
i = n+1
while i <= nn: # initial conditions for linearized maps
v[i]=0 # Don't mess with these; they are problem independent!
i = i + 1
i = 1
while i <= n:
v[(n+1)*i]=1
ltot[i]=0
i = i + 1
#print "v = ",v
t=0
w = 0
while (w < cmax):
j = 1
while j <= irate:
i = 1
while i <= nn:
x[i]=v[i]
i = i + 1
[x, xnew] = derivs(x, xnew, n)
i = 1
while i <= nn:
A[i] = xnew[i]
x[i] = v[i] + (h*A[i]) / 2.0
i = i + 1
[x, xnew] = derivs(x, xnew, n)
i = 1
while i <= nn:
B[i] = xnew[i]
x[i] = v[i] + (h*B[i]) / 2.0
i = i + 1
[x, xnew] = derivs(x, xnew, n)
i = 1
while i <= nn:
C[i] = xnew[i]
x[i] = v[i] + h*C[i]
i = i + 1
[x, xnew] = derivs(x, xnew, n)
i = 1
while i <= nn:
D[i] = xnew[i]
v[i] = v[i] + h*(A[i] + D[i] + 2*(B[i] + C[i]))/6.0
i = i + 1
t = t + h
j = j + 1
#construct new orthonormal basis by gram-schmidt:
znorm[1]=0 #normalize first vector
j = 1
while j <= n:
znorm[1]=znorm[1]+v[n*j+1]**2
j = j + 1
znorm[1] = math.sqrt(znorm[1])
j = 1
while j <= n:
v[n*j+1]=v[n*j+1]/znorm[1]
j = j + 1
#generate new orthonormal set:
j = 2
while j <= n:
k = 1
while k <= j-1:
gsc[k]=0
l = 1
while l <= n:
gsc[k]=gsc[k]+v[n*l+j]*v[n*l+k]
l = l + 1
k = k + 1
k = 1
while k <= n: #construct a new vector
l = 1
while l <= j-1:
v[n*k+j]=v[n*k+j]-gsc[l]*v[n*k+l]
l = l + 1
k = k + 1
znorm[j]=0 #calculate the vector's norm
k = 1
while k <= n: #construct a new vector
znorm[j]=znorm[j]+v[n*k+j]**2
k = k + 1
znorm[j]=math.sqrt(znorm[j])
k = 1
while k <= n: #normalize the new vector
v[n*k+j] = v[n*k+j] / znorm[j]
k = k + 1
j = j + 1
k = 1
while k <= n: #update running vector magnitudes
if znorm[k] > 0:
ltot[k] = ltot[k] + math.log(znorm[k])
k = k + 1
m = m + 1
if m % io == 0 or w == cmax-1: #normalize exponent and print every io iterations
lsum=0
kmax=0
k = 1
while k <= n:
le = ltot[k] / t
lsum = lsum + le
if lsum > 0:
lsum0 = lsum
kmax = k
k = k + 1
w = w + 1
if choice == 1:
print "Rossler:"
else:
print "Lorenz:"
print n, "LEs = "
lsum=0
kmax=0
k = 1
while k <= n:
le = ltot[k] / t
lsum = lsum + le
if lsum > 0:
lsum0 = lsum
kmax = k
print le
k = k + 1
Similar to Neural Networks, Support Vector Machines (SVM) are powerful modelling techniques in machine learning that can be applied to a variety of tasks. SVMs are primarily used to classify data while one of its variants, Support Vector Regression (SVR), can be used for time series analysis. In this post, we will perform SVR on chaotic time series using common Python libraries and a simple wrapper program.
The following libraries will be required to use the script I will demonstrate using Python 2.7:
Numpy
Scipy
Scikits.Learn (a toolbox that hosts a variety of machine learning algorithms and requires Numpy and Scipy)
MatPlotLib (used for graphing and visualization)
Other libraries such as SetupTools may be required depending on your system
With the packages installed, we can begin to build a program that fits an SVR model to chaotic data. First, we define a time series. In this case, I will use the Delayed Henon Map with a delay of 4, for more information please see my posts on that system.
The model will be embedded in a d dimension space.
def DelayedHenon(cmax, transient, d):
temp = []
rF = [] # Features
rT = [] # Targets
dX = []
i = 0
while i < d:
dX.append(0.1)
i = i + 1
while i <= cmax + transient:
x = 1 - 1.6 * dX[i-1] ** 2 + 0.1 * dX[i-4]
if i > transient-d:
temp.append(x)
if i > transient:
rT.append(x)
dX.append(x)
i = i + 1
rF = SlidingWindow(temp, d)
return {"F":rF, "T":rT}
I added a function that formats the data properly so it can be used by the SVR Fit function:
def SlidingWindow(x, d):
i = d
y = []
while i < len(x):
temp = []
j = 1
while j <= d:
temp.append(x[i - d + j])
j = j + 1
y.append(temp)
i = i + 1
return y
We then define the model, note that SVR has a few user defined parameters which can be chosen to fit the data better.
In this case, you can choose the following parameters to produce a model with an error near 0.0035 for 1024 points taken from the time series. In practice, you can create a training algorithm to find these values.
gam = 0.117835847192
c = 1.02352954164
eps = 0.00163811344978
tolerance = 0.00539604526663
d = 5
cmax = 1024
Finally, we create functions to fit the time series and plot the results:
def Fit(cmax, d, tolerance, gam, c, eps):
data = TimeSeries(cmax, random.randint(0,10000), d)
test = TimeSeries(cmax, random.randint(0,10000), d)
x = data["F"]
y = data["T"]
sout = test["T"]
clf = svm.SVR(kernel='rbf', degree=d, tol=tolerance, gamma=gam, C=c, epsilon=eps)
try:
m = clf.fit(x, y)
except:
return [1000, clf, sout, sout]
mout = clf.predict(test["F"])
# Determine the error for the system
err = 0
i = 0
while i < len(mout):
err += (mout[i] - sout[i])**2
i = i + 1
err = math.sqrt(err / (len(mout)+.0))
return [err, clf, mout, sout]
def Plotter(err, mout, sout, SaveGraph):
plt.title("SVR: Err=" + str(err))
plt.scatter(sout[0:len(sout)-2], sout[1:len(sout)-1], s=1, marker='o', edgecolors='none')
plt.scatter(mout[0:len(mout)-2], mout[1:len(mout)-1], s=1, marker='o', edgecolors='none', facecolor='r')
plt.show()
if SaveGraph:
try:
plt.savefig("svr")
except:
# Do Nothing
SaveGraph = SaveGraph
plt.clf()
Here is the final program after putting all of the functions together (and moving some stuff around) with example output:
#!/usr/bin/env python
from scikits.learn import svm
import random
import math
import matplotlib.pyplot as plt
d = 5 # Degree of the SVR model
cmax = 1024 # Number of points to run
def TimeSeries(cmax, transient, d):
return DelayedHenon(cmax, transient, d)
def main():
global cmax, d, perturbation
gam = 0.117835847192
c = 1.02352954164
eps = 0.00163811344978
tolerance = 0.00539604526663
# Create and Fit a model to data, we just need to create support vectors
[err, svr, mout, sout] = Fit(cmax, d, tolerance, gam, c, eps)
Plotter(err, mout, sout, False)
print "error: ", err
def DelayedHenon(cmax, transient, d):
temp = []
rF = [] # Features
rT = [] # Targets
dX = []
i = 0
while i < d:
dX.append(0.1)
i = i + 1
while i <= cmax + transient:
x = 1 - 1.6 * dX[i-1] ** 2 + 0.1 * dX[i-4]
if i > transient-d:
temp.append(x)
if i > transient:
rT.append(x)
dX.append(x)
i = i + 1
rF = SlidingWindow(temp, d)
return {"F":rF, "T":rT}
def Fit(cmax, d, tolerance, gam, c, eps):
data = TimeSeries(cmax, random.randint(0,10000), d)
test = TimeSeries(cmax, random.randint(0,10000), d)
x = data["F"]
y = data["T"]
sout = test["T"]
clf = svm.SVR(kernel='rbf', degree=d, tol=tolerance, gamma=gam, C=c, epsilon=eps)
try:
m = clf.fit(x, y)
except:
return [1000, clf, sout, sout]
mout = clf.predict(test["F"])
# Determine the error for the system
err = 0
i = 0
while i < len(mout):
err += (mout[i] - sout[i])**2
i = i + 1
err = math.sqrt(err / (len(mout)+.0))
return [err, clf, mout, sout]
def Plotter(err, mout, sout, SaveGraph):
plt.title("Delayed Henon Map and SVR Model")
p1 = plt.scatter(sout[0:len(sout)-2], sout[1:len(sout)-1], s=1, marker='o', edgecolors='none')
p2 = plt.scatter(mout[0:len(mout)-2], mout[1:len(mout)-1], s=1, marker='o', edgecolors='none', facecolor='r')
plt.legend([p2, p1], ["Model (Err="+ str(round(err,4))+")", "Map"], loc=8)
plt.draw()
plt.show()
if SaveGraph:
try:
plt.savefig("svr")
except:
# Do Nothing
SaveGraph = SaveGraph
plt.clf()
def SlidingWindow(x, d):
i = d
y = []
while i < len(x):
temp = []
j = 1
while j <= d:
temp.append(x[i - d + j])
j = j + 1
y.append(temp)
i = i + 1
return y
def WriteFile(err, BestTol, BestGam, BestC, BestEps, s):
# Write the model, sensitivities to a file
f = open("sensitivities", "w")
f.write("Model Parameters:\n")
f.write("d=" + str(d) + "\n")
f.write("gam=" + str(BestGam) + "\n")
f.write("c=" + str(BestC) + "\n")
f.write("eps=" + str(BestEps) + "\n")
f.write("tolerance=" + str(BestTol) + "\n")
f.write("Model Data:\n")
f.write("err=" + str(err) + "\n")
f.write("s=" + str(s) + "\n")
f.close()
main()
Example Output:
err: 0.00381722720161
Example Output
To fit data taken from other systems, simply change the Timeseries function.
At work, we use VB.net and one problem that I encountered involved a forum where they wanted to automatically link any words with http:// or https:// within the posts. Aside from security concerns, you can do this using the following script.
However, it is important to note, that it will only link URLS that are one word. You could modify the regular expression so that if a person put quotes around the URL, they could link the entire thing, but I am going to leave that to someone else.
' Function: LinkText()
' Takes a msg and automatically adds HTML to link all words starting
' http:// or https:// in the text
' Parameters
' Msg: Is the Message we want to link the links in
' Returns:
' NewMsg: The linked the links in the text
Function LinkText(ByVal Msg As String)
' Remove all hrefs (so we can replace them)
Msg = Msg.Replace("", "")
Msg = System.Text.RegularExpressions.Regex.Replace(Msg, "]*>?","")
Dim MatchObj as MatchCollection = Regex.Matches(Msg,
"(htt(p|ps)://.+?)($|\s+)")
Dim NewMsg As String = Msg
For Each m As Match In MatchObj
NewMsg = System.Text.RegularExpressions.Regex.Replace(NewMsg,
m.Groups(0).Value, "<a
href=""" & m.Groups(1).Value & """>"
& m.Groups(1).Value & "</a>"
& m.Groups(3).Value)
Next
Return NewMsg
End Function
Using the iterative 4-th order Runge-Kutta method as described here, we can create low dimensional slices of the system’s attractor known as Poincare Sections, Return maps, or Recurrence maps.
with a, b, and c set to 0.2, 0.2, and 5.7, respectively.
A short time series from the Rosseler Attractor
Poincare sections are important for visualizing an attractor that is embedded in greater than 3 spatial dimensions. Additionally, certain dynamical and topological properties are invariant to these transformations such as Lyapunov exponents and bifurcation diagrams (See Sprott 2003). However, the Lyapunov exponent with a value of 0 is lost due to the removal of the direction of the flow.
To build a Poincaré section, imagine that you have selected a slice of phase space. Every time the orbit passes through that slice, we collect a point. We repeat this over and over to create an image of what the attractor looks like in that slice of space.
Return map for the Y and Z variables where the value of X is between 3.5 and 3.75
We can also find local maximums to create a return map. In the following graph, we look for points that are local maximums and plot them versus the previous local maximum.
The local maximums of the variable X versus previous local maximums of the variable X
The following Python code uses the 4-th order Runge-Kutta method to create this return map.
import math, operator, random
dataX0 = []
dataY0 = []
dataZ0 = []
yList = []
xList = []
zList = []
count = 1
t = 1
x = -9
y = 0
z = 0
xList.append(x)
yList.append(y)
zList.append(z)
h = .01
a = .2
b = .2
c = 5.7
def f(x,y,z):
global a,b,c
dxdt = -y-z
return dxdt
def g(x,y,z):
global a,b,c
dydt = x + a * y
return dydt
def e(x,y,z):
global a,b,c
dzdt = b + z * (x - c)
return dzdt
def rk4o(x, y, z):
global h
k1x = h*f(x, y, z)
k1y = h*g(x, y, z)
k1z = h*e(x, y, z)
k2x = h*f(x + k1x/2.0, y + k1y/2.0, z + k1z/2.0)
k2y = h*g(x + k1x/2.0, y + k1y/2.0, z + k1z/2.0)
k2z = h*e(x + k1x/2.0, y + k1y/2.0, z + k1z/2.0)
k3x = h*f(x + k2x/2.0, y + k2y/2.0, z + k2z/2.0)
k3y = h*g(x + k2x/2.0, y + k2y/2.0, z + k2z/2.0)
k3z = h*e(x + k2x/2.0, y + k2y/2.0, z + k2z/2.0)
k4x = h*f(x + k3x, y + k3y, z + k3z)
k4y = h*g(x + k3x, y + k3y, z + k3z)
k4z = h*e(x + k3x, y + k3y, z + k3z)
x = x + k1x/6.0 + k2x/3.0 + k3x/3.0 + k4x/6.0
y = y + k1y/6.0 + k2y/3.0 + k3y/3.0 + k4y/6.0
z = z + k1z/6.0 + k2z/3.0 + k3z/3.0 + k4z/6.0
return [x,y,z]
X0 = []
Y0 = []
Z0 = []
previousLocalMax = x
localMax = x
t = 1
changeInTime = h
avgChange = []
LastPoint = changeInTime
while changeInTime < 20000 and len(dataX0) < 1000:
[x,y,z] = rk4o(xList[t-1], yList[t-1], zList[t-1])
xList.append(x)
yList.append(y)
zList.append(z)
if 201 < changeInTime:
# Print y and z points when x is between 3.5 and 3.75 (as shown in the figure)
if x < 3.75 and x > 3.5:
print y , "," , z
if x < xList[t-1] and xList[t-2] < xList[t-1]:
previousLocalMax = localMax
localMax = xList[t-1]
dataX0.append(localMax)
dataY0.append(previousLocalMax)
# Calculate the change between this and last point
avgChange.append(changeInTime - LastPoint)
LastPoint = changeInTime
t = t + 1
changeInTime += h
# Print the average steps between points
print "Average Number of Steps between points:", sum(avgChange[1:]) / float(len(avgChange[1:]))
print "First Return Map Length: " + str(len(dataX0))
f = open("rosseler.csv", "w")
i = 0
while i < len(X0) and i < 10000:
f.write(str(X0[i]) + "," + str(Y0[i])+ "," + str(Z0[i])+"\n")
i = i + 1
f.close()
f = open("fr-rossler-x", "w")
i = 0
while i < len(dataX0):
f.write(str(dataX0[i])+"\r")
i = i + 1
f.close()
Ordinary Differential Equations (ODEs) can be used to define systems in fields as varied as biology to engineering to mathematics. These equations express a relationship between an unknown function and its derivative. One example, the Rossler attractor,
A short time series from the Rossler Attractor
We can produce a time series from these equations by solving the equations using the iterative 4-th order Runge-Kutta method and plugging each of the solutions back into the equations. Along with the following script which allows you to implement the Runge-Kutta method on ODEs, I have included code to numerically estimate the largest Lyapunov exponent. The largest Lyapunov exponent is used to indicate chaos, or sensitive dependence to initial conditions, within a system. For the Rossler attractor, defined parameters (a=b=.2, c=5.7), and initial conditions (x=-9, y=z= 0), the largest Lyapunov exponent is about 0.0714. Numerous resources are available for more information about ordinary differential equations and other systems that you may want to explore with this script.
import math, random
# Step Size
h = .001
#Initial conditions for Rossler system
x = -9
y = 0
z = 0
# Parameters for the Rossler System
a = .2
b = .2
c = 5.7
# The perturbation used to calculate the Largest Lyapunov exponent
perturb = .000000001
# Functions that define the system
def f(x,y,z):
global a,b,c
dxdt = -y-z
return dxdt
def g(x,y,z):
global a,b,c
dydt = x + a * y
return dydt
def e(x,y,z):
global a,b,c
dzdt = b + z * (x - c)
return dzdt
# randomly perturb the initial conditions to create variable time series
x = x + random.random() / 2.0
y = y + random.random() / 2.0
z = z + random.random() / 2.0
dataX0 = []
dataY0 = []
dataZ0 = []
yList = []
xList = []
zList = []
lamdaList = []
lyapunovList = []
t = 1
xList.append(x)
yList.append(y)
zList.append(z)
# Use the 4th order Runge-Kutta method
def rk4o(x, y, z):
global h
k1x = h*f(x, y, z)
k1y = h*g(x, y, z)
k1z = h*e(x, y, z)
k2x = h*f(x + k1x/2.0, y + k1y/2.0, z + k1z/2.0)
k2y = h*g(x + k1x/2.0, y + k1y/2.0, z + k1z/2.0)
k2z = h*e(x + k1x/2.0, y + k1y/2.0, z + k1z/2.0)
k3x = h*f(x + k2x/2.0, y + k2y/2.0, z + k2z/2.0)
k3y = h*g(x + k2x/2.0, y + k2y/2.0, z + k2z/2.0)
k3z = h*e(x + k2x/2.0, y + k2y/2.0, z + k2z/2.0)
k4x = h*f(x + k3x, y + k3y, z + k3z)
k4y = h*g(x + k3x, y + k3y, z + k3z)
k4z = h*e(x + k3x, y + k3y, z + k3z)
x = x + k1x/6.0 + k2x/3.0 + k3x/3.0 + k4x/6.0
y = y + k1y/6.0 + k2y/3.0 + k3y/3.0 + k4y/6.0
z = z + k1z/6.0 + k2z/3.0 + k3z/3.0 + k4z/6.0
return [x,y,z]
t = 1
changeInTime = h
startLE = True
while changeInTime < 20000: # Perform 20000 / h iterations
[x,y,z] = rk4o(xList[t-1], yList[t-1], zList[t-1])
xList.append(x)
yList.append(y)
zList.append(z)
if 200 < changeInTime: # Remove the transient after 200 / h iterations
if startLE:
cx = xList[t-1] + perturb
cy = yList[t-1]
cz = zList[t-1]
startLE = False
# Calculate the Largest Lyapunov Exponent
[cx, cy, cz] = rk4o(cx, cy, cz)
delx = cx - x
dely = cy - y
delz = cz - z
delR1 = ((delx)**2+(dely)**2+(delz)**2)
df = 1.0 / (perturb**2) * delR1
rs = 1.0 / math.sqrt(df)
cx = x + rs*delx
cy = y + rs*dely
cz = z + rs*delz
lamda = math.log(df)
lamdaList.append(lamda)
#if t % 1000 == 0: # Print the Lyapunov Exponent as you go
# print t, " ", .5*sum(lamdaList) / (len(lamdaList)+.0) / h
t = t + 1
changeInTime += h
lyapunov = .5*sum(lamdaList) / (len(lamdaList)+.0) / h
print lyapunov
# Output the x-component to a file
f = open("rossler-x", "w")
i = 0
while i < len(dataX0):
f.write(str(dataX0[i])+"\r")
i = i + 1
f.close()
A relatively simple and straightforward way to post information to a server using Python.
import urllib
# Tries to open the url with the params through the method specified
def fetch_url(url, params, method):
params = urllib.urlencode(params)
if method=="GET":
f = urllib.urlopen(url+"?"+params)
else:
# Usually a POST
f = urllib.urlopen(url, params)
return (f.read(), f.code)
url = "http://google.com"
method = "POST"
params = {"Param1": "Value1"}
# Fetch the content and response code
[content, response_code] = fetch_url(url, params, method)
# Check if the server responded with a success code (200)
if (response_code == 200):
print content
else:
print response_code
