Machine Learning Python: Wavelet Spectrum examples

Another test on my MacBook Pro:  The wavelet transform example in Megies et al (2011).

You may need to install the Machine Learning Python package.  The command

conda install mlpy

did not work for me, so I search the internet for a method to install mlpy

The mlpy package requires the GNU Scientific Library (GSL), which was not installed on my laptop, and may not be installed on yours either.  This package can be found on the web.  My method was taken from the book "Practical Data Analysis Cookbook" by Tomasz Drabas.  First, direct your web browser to

http://gnu.mirro.iweb.com/gsl

and download the latest version of GSL, typically called gsl-latest.tar-gz. Currently the latest is V2.4. Move the downloaded file into your Applications directory (on a Mac) -- I dont know the equivalent place on a Windows machine.  Then unpack this gzipped tarfile by double-clicking on it ot else use the command

gzip -c gsl-latest.tar-gz | tar -xv

The unpacked file will produce a folder with the GSL software in it.  In your terminal-app go to this folder and prepare and install GSL

./configure 

make

make install

Hopefully this compilation and installation proceeds without error.  Now download a compressed tarfile with the lastest version of mlpy from this website into your anaconda dirctory, which on my laptop is /Users/jjpark/anaconda.  The current latest version is V3.5.

http://sourceforge.net/projects/mlpy/files

Double-click the tarfile to extract its contents into a subdirectory, or else

gzip -c mlpy-.tar.gz | tar -xv

or, if already uncompressed,

tar -xvf mlpy-.tar

This action produces a subdirectory mlpy-.  cd to this directory in your terminal app and type

python setup.py install

and mlpy will be installed if all goes well.  Python modules that one obtains from sourceforge and other websites typically include a setup.py scriptfile to instruct python how to set the module up for later use.

I discovered coding bugs in the mlpy module while trying to execute the coding examples.  The continuous wavelet transform relies on a script file ./mlpy-3.5.0/mlpy/wavelet/continuous.py that had a mistake in its Line 162:

161     J = floor(dj**-1 * log2((N * dt) / s0))
162     s = empty(J + 1)
   
The python function "floor" from the numpy module returns a floating-point number, but the python function "empty" expects an integer to define how long an array to generate. This combination of mismatched variable types leads to an error message:

'numpy.float64' object cannot be interpreted as an integer

Fixing the code to read, e.g.,

     J = floor(dj**-1 * log2((N * dt) / s0))
     JJ = int(J)
     s = empty(JJ + 1)

*should* solve the problem, but it doesn't.  Why? Because the original continuous.py is loaded into the python libraries during the install-operation, and changing the code after installation doesnt change the installed code.  You must make changes to the script continuous.py in the mlpy-3.5.0 subdirectory, and then re-install the module by executing

python setup.py install

while the terminal app is "in" the ./mlpy-3.5.0 directory.  Note that python will tell you that the execution error resides in a copy of continuous.py that lies deep within a subdirectory tree.  This is NOT the copy that needs to be debugged.  You must debug the copy of continuous.py that sits close within the subdirectory anaconda/mlpy-3.5.0/mlpy/wavelet, and then re-install the mlpy module.

Here is a simple example of a wavelet transform of a randomly-generated array of numbers

import numpy as np
import matplotlib.pyplot as plt
import mlpy.wavelet as wave
x = np.random.sample(512)
scales = wave.autoscales(N=x.shape[0], dt=1, dj=0.025, wf='morlet', p=2)
X = wave.cwt(x=x, dt=1, scales=scales, wf='dog', p=2)
fig = plt.figure(1)
ax1 = plt.subplot(2,1,1)
p1 = ax1.plot(x)
ax1.autoscale_view(tight=True)
ax2 = plt.subplot(2,1,2)
p2 = ax2.imshow(np.abs(X), interpolation='nearest')
plt.show()


Here is the wavelet transform of the data series we downloaded in the first ObsPy exercise

import numpy as np, matplotlib.pyplot as plt
import mlpy 
import mlpy.wavelet as wave
from obspy.core import UTCDateTime
from obspy.clients.arclink import Client
from obspy.signal.filter import lowpass
from obspy.signal.invsim import corn_freq_2_paz, simulate_seismometer
client = Client("webdc.eu", port=18001)
t = UTCDateTime("2009-08-24 00:20:03")
st = client.get_waveforms("BW", "RJOB", "", "EHZ", t, t+30)
omega0 = 2.0
trace = st[0].data
npts = len(trace)
dt = st[0].stats.delta
dj = 0.05
wf = 'morlet'
p=omega0
scl = mlpy.wavelet.autoscales(npts, dt, dj , wf, p) 
spec = mlpy.wavelet.cwt(trace, dt, scl , wf='morlet', p=omega0) 
freq = (omega0 + np.sqrt(2.0 + omega0**2)) / (4*np.pi * scl[1:])
t = np.arange(st[0].stats.npts) / st[0].stats.sampling_rate 
plt.imshow(np.abs(spec), extent=(t[0], t[-1], freq[-1], freq[0])) 
plt.xlabel('Time [s]') 
plt.ylabel('Frequency [Hz]') 
plt.show()

Try these scripts to be sure that your installation of Python and ObsPy is working properly.