""" Laser_data.py: python Laser_data.py data_dir Script that makes data for the following laser figures: In Chapter 1: LaserLP5.pdf Figure 1.1 LaserLogLike.pdf Figure 1.2 LaserStates.pdf Figure 1.3 LaserForecast.pdf Figure 1.4 In Chapter 3: LaserHist.pdf Figure 3.1 Requires: /LP5.DAT Tang's real laser data Produces: / LaserLP5 Simulated laser data for figure LaserLP5 LaserLogLike Dependence of LL on s and b for figure LaserLogLike LaserStates Sequence of 250 states from EKF for figure LaserStates LaserForecast Sequence of 400 scalar values for figure LaserForecast LaserHist Histogram 600 values for figure LaserHist in Chap 3 Copyright (c) 2005, 2008 Andrew Fraser This file is part of HMM_DS_Code. HMM_DS_Code is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version. HMM_DS_Code is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. """ import sys, os, math, numpy, numpy.linalg as LA, EKF, copy, lorenz import pickle, scipy.optimize as O assert len(sys.argv) == 2,"Try: python Laser_data.py data" data_dir = sys.argv[1] def ParamCalc(rad, Tr, r, s, b, ts, offset, scale, dev_epsilon, dev_eta): """ Translates specification of initial condition (rad,Tr) with radius rad from one of the unstable complex fixed points and relaxation time Tr to (ic[0],ic[1],ic[2]) thus creating a parameter vector suitable for EKF.LogLike(). """ root = math.sqrt(b * (r - 1)) # A frequently used intermediate term f = numpy.array([root, root, (r - 1)]) # One of three fixed points DF = numpy.array([[-s, s, 0],[1, -1, -root], [root, root, -b]]) # DF is the derivative of F at the fixed point f [val, vec] = LA.eig(DF) v = vec.T[1].real u = vec.T[1].imag w = (u - (u[2] / v[2]) * v) x0 = f - rad * (f[0] / w[0]) * w # x0 is in the plane of unstable complex pair of eigenvectors at # radius rad from f ic = lorenz.Lsteps(x0, s, b, r, Tr, 2)[-1,:] # Now ic is result of trajectory starting at x0 after relaxation # time Tr. The idea is avoid funny looking transients. return([ic[0], ic[1], ic[2], r, s, b, ts, offset, scale, dev_epsilon, dev_eta]) def StepNt(P, Nt): """ Simulate a sequence of Nt measurements by: (1) Calculate initial condidtions from radius and time. (2) Integrate Lorenz. (3) Square first component, scale, and add offset. P = [rad, Tr, r, s, b, ts, offset, scale, noise] """ params = ParamCalc(*(list(P[0:8])+[0,0])) # Stick 0, 0 on end of P ic = params[0:3] r,s,b,ts,offset,scale = params[3:9] X = lorenz.Lsteps(ic,s,b,r,ts,Nt) ys = X[:, 0] * X[:, 0] * scale + offset return ys def error250(P): """Calculate an error for the difference between the first 250 points of actual laser data and the simulation that the parameters in P specify. I use this with O.fmin to improve the match between the simulation and the data for the first 250 points """ global Y, best rad = P[0] if rad*rad > 0.25: return 1e20 try: sim = StepNt(P, 250) except ValueError, err: #print 'call StepNt failed at relax with ', err return 1e20 e = 0.0 for t in xrange(250): e += ((Y[t] - sim[t]) ** 2) / Y[t] if e < best-1: best = e print "best=%f"%best return e def NegLogLike(P): """Use EKF.LogLike to calculate the log likelihood of the first Nt values in the Y sequence. """ global Y, Nt,best rad,Tr = P[0:2] devEp,devEta = P[8:10] if rad * rad > 0.25 or Tr < 0.5 or devEp < 0.0 or devEta < 0.0: return 1e20 # If parameters are bad, return big number try: L = ParamCalc(*P[:10]) cost = -EKF.LogLike(L, Y[0:Nt]) + (devEta/1e-2)**2 except (ValueError, RuntimeError ), err: #print 'NegLogLike returning 1e20 because ', err return 1e20 if cost < best-1: best = cost print "best=%f"%best return cost def NLL_Fixed_Noise(P): """Use EKF.LogLike to calculate the log likelihood of the first Nnll values in the Y sequence. Like NegLogLike() but noise is fixed. """ global Y, best Nnll = 250 rad = P[0] Tr = P[1] if rad ** 2 > 0.25 or Tr < 0.5: return 1e20 try: L = ParamCalc(*(list(P[0:8])+[4.0/P[7], 1e-7])) # P[7] is scale. The above line sticks dev_eps = 4.0 channels # and dev_eta = 1e-7 on to the end of the parameter list and # converts from (rad,Tr) to ic. LL = EKF.LogLike(L, Y[0:Nnll]) except ValueError, err: #print 'NLL_Fixed_noise returning 1e20 because ', err return 1e20 if -LL < best-1: best = -LL print "best=%f"%best return(-LL) def P_print(O_Par): print """ ic=[%7.4f, %7.4f, %7.4f], r=%-7.4f, s=%-7.4f, b=%-7.4f, ts=%7.4f, offset=%-7.4f, scale=%-7.4f, sigma_epsilon=%-7.4f, sigma_eta=%-7.4g """%tuple(O_Par) # Read the input data. The first and last line are not data, skip them Y = [int(x.split()[1]) for x in (open( data_dir+'/LP5.DAT','r').readlines())[1:-1]] ################## Begin Making LaserLP5 ############################ """ LaserLP5: Contains the first 250 points of laser data file and 250 simulated points. To get parameters for the simulation, apply fmin to go to a local minimum starting from a point that I found by broadly and coarsely surveying parameter space. I did the broad survey with an old script called survey4.py that is not part of this distribution. The orbits go around each lobe 5 times before swithcing. The optimization finds an unstable period 5 orbit that closely matches the laser data. Note that the simulated orbit is unstable, while the laser data certainly comes from a stable orbit.""" # From Survey4.py #ic = [ 4.09750891, 4.70477676, 22.17708015 ] r = 21.16 s = 1.792 b = 0.366 ts = 0.143 offset = 15.95 scale = 7.5 rad = 0.32915 Tr = 1.645 Params = [rad, Tr, r, s, b, ts, offset, scale] try: # This try/except speeds debugging by saving results of long calculations R = pickle.load(open(data_dir+'/opt_250_results','r')) except: best = 1e20 print 'Before optimization of error250()' R = O.fmin_powell(error250, Params, xtol = 1e-6, ftol = 1e-6) pickle.dump(R,open(data_dir+'/opt_250_results','w')) print 'After optimization of error250(), error=%6.2f'%error250(R) Nt = 250 sim = StepNt(R,Nt) f = open(data_dir + '/LaserLP5','w') for t in xrange(Nt): print >>f, t, Y[t], sim[t] f.close() ############# Finished making LaserLP5 ################################ ################## Begin Optimization ###################################### # print """ I want to find parameters that yield a trajectory with some or all of the following properties: 1. Has a stable period 5 orbit 2. Makes LaserStates do 5 orbits about each unstable fixed point 3. Give a forecast that is close out to at least 500 Ideas to try: 1. Use smoothing to get better orbits about which to linearize 2. Search for period 5 orbits 3. Search for stable period 5 orbits With a procedure that takes r,s,b and finds a period 5 orbit, I could penalize for instability and also solve/optimize other parameters in the following order: 1. Ts. From known orbit time for data 2. IC. Best fit is just 1-d shift to line up simulated peaks to data 3. Scale and offset. Call optimizer 4. Sigma_eta and Sigma_epsilon. Call optimizer on Kalman filter function. Setup affine maps using periodic orbit and do filtering without changing the maps. """ # Params = [rad, Tr, r, s, b, ts, offset, scale, dev_Ep, dev_eta] # Params = [ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9] try: R = pickle.load(open(data_dir+'/opt_Fixed_results','r')) except: best = 1e20 print 'Before optimization of NLL_Fixed_Noise()' R=O.fmin_powell(NLL_Fixed_Noise,Params,xtol=1e-5,ftol=1e-5) pickle.dump(R,open(data_dir+'/opt_Fixed_results','w')) print 'After optimization NLL_Fixed_Noise()=%f'%NLL_Fixed_Noise(R) Params = list(R) + [4.0/R[7], 1e-7] P_print(ParamCalc(*Params[0:10])) try: P = pickle.load(open(data_dir+'/opt_General_results','r')) except: best = 1e20 print 'Before general optimization' P = O.fmin_powell(NegLogLike,Params) pickle.dump(P,open(data_dir+'/opt_General_results','w')) print 'After general optimization\nNegLogLike(P)=',NegLogLike(P) O_Par = ParamCalc(*P[0:10]) P_print(O_Par) ############# End of Optimization ##################################### ### Begin Making LaserLogLike (log likelihood dependence on parameters) ##### s = O_Par[4] b = O_Par[5] Nstep = 5 # Number of steps away from peak in each direction Frange = numpy.arange(-Nstep,Nstep+1,1.0) s_range = (Frange * 5e-3/Nstep + 1.0)*s # Centered at s_max with given step size b_range = (Frange * 5e-3/Nstep + 1.0)*b Par = copy.copy(O_Par) f = open(data_dir+'/LaserLogLike', 'w') for ss in s_range: for bb in (b_range): Par[4] = ss Par[5] = bb ll = EKF.LogLike(Par, Y[0:Nt]) print >>f, ss, bb, ll print >>f, ' ' f.close() ################## Finished Making LaserLogLike ###################### #### Begin Making LaserStates (a trajectory in 3-d state space) ###### Nt = 250 aM = Nt*[None] ll = EKF.LogLike(O_Par,Y[0:Nt],aM=aM) f = open(data_dir+'/LaserStates', 'w') for x in aM: print >>f,x.A[0,0],x.A[2,0] f.close() ################## Finished Making LaserStates ##################### ################## Begin Making LaserForecast ###################### Nf = 400 # Number to forecast ic = aM[-1] X = lorenz.Lsteps(ic,s,b,r,ts,Nf+1) f = open(data_dir+'/LaserForecast','w') for t in xrange(Nf): x0 = float(X[t+1,0]) print >>f, t+Nt,Y[Nt+t],(x0*x0*scale + offset) f.close() ################## Finished Making LaserForecast ################## ### Begin making LaserHist (histogram of the first 600 samples) ### Count = numpy.zeros(256,numpy.int32) for y in Y[0:600]: Count[y] += 1 f = open(data_dir + '/LaserHist', 'w') for t in xrange(0,100): print >>f,t,Count[t] f.close() ############ Finished making LaserHist ################################# # Local Variables: # mode: python # End: