This week I wrote Matlab code, most of which can be found in /users/masha

First, I wrote a simulation seismicFilter.m which filters noisy seismic noise with a desired signal of non-seismic noise. The signals are purely simulated, so I played around with zero-pole-gain generation of transfer functions to obtain them. The function takes the number of taps, the filter type (Wiener or adaptive nlms) as well as an iteration step size and number of iterations, and generates PSD plots of the witness signal, the desired signal, the estimated (filtered) signal, and the error. I'm not sure that I am properly implementing the Wiener part of the code, and I assume the line "[W, R, P] = miso_firlev(TAPS, noisySeismicSignal1, seismicSignal2); " generates W, a filter with TAPS number of weights, but  then "[y, error] = filter(W, 1, noisySeismicSignal1);" generates an error signal of size TAPS rather than N, the size of the original signal. Perhaps I should calculate error using e(t) = d(t + a) - w(t)*x(t), where "a" is the delay.

I have various screenshots in my directory of what seismicFilter.m generates, and I will take a larger screenshot, as well as generate a learning curve (for error vs. number of taps) when I can use Sasha's computer for a bit, since it both has more computing power and a larger screen.

The funciton filterConvergence.m, meanwhile, is similar, except it takes two file names as real data, and uses realDataFilter.m to run the filtering. Currently, I am working with data from C1:IOO-MC_F_DQ-Online  and C1:PEM-SEIS_GUR1_X_IN1_DQ-Online, and I will include screenshots of these once I get on Sasha's computer.

In order to generate the data, meanwhile, I had to modify the python script, and thus wrote mashaImportingData.py for myself. Likewise, plotSignalFromFile.m visualizes this data, both in the time domain and in the frequency domain.

On the side, I wrote an NLMS filter in adaptiveFilterSimulationNLMS.m, and compared is to Matlab's NLMS filter in NLMStest.m (using generated data) and adaptiveFilterSimulation.m using twn input signals. Right now, it's faster on smaller inputs and smaller tap sizes, but then begins to choke and run slower than the Matlab one when these get to big. In order to improve it, I have to develop a better method of generating the initial weights.

As far as machine learning goes, once I find the number of taps for the convergence of both the Wiener filter and the NLMS filter, I will email Denis for further instructions. At some point, however, I should generate learning sets from the seismometers and the MCL (or the DARM), and thus have to find adequate times at which I can take data (probably not from the DARM, however, because it was rarely on).

Thanks for reading!

Hi everybody,

Last night I (with the help of Jenne and Jenne's advice - not to implicate her in this or anything) changed the filters for GUR1, GUR2, and STS in C1:PEM-RMS, adding a butterworth bandpass filter at each corresponding frequency band as well as a gain to convert from counts to micros/sec, and then adding a low pass filter in case of aliasing upon squaring.

Currently the seismic signals are going crazy, and producing "Nan" output on the strip graph (which leads to the instantaneously sharp spikes - which leads to the entire signal being filled on the visualizer on the wall). I checked the DataViewer output and the tdsdata output using both grep and wc, and it seems that both every single signal point is present and is a real number (also not a small real number, thereby debunking floating-point error). I'm currently not sure why seismic-strip reads out 'Nan' - perhaps because it's taking the log of 0, taking a negative log, taking the root of a negative number, or dividing by zero.

Does anyone know if the seismic-strip Nan issue is a program bug? If it's not (and therefore a filter bug), please let me know as well.

I'll be in lab for the rest of the night changing the butterworth filters to odd-order elliptic filters (at Rana's suggestion), as well as changing the cut-off frequency for the low-pass filters.

I'll E-log about it when I'm done.

Just to be sure that my numbers are correct - The STS, GUR1, and GUR2 channels all have gain 10, right? (I parsed through the e-log, and these seem to be the most recent numbers

Thanks for your help,

Masha

UPDATE: I changed all of the GUR1Z channels to order-5 elliptic filters. I approximated the attenuation for each one by setting the integral from _CutoffFrequency to 10^3 Hz of 10^(-Percent(f)/20) df to 0.01, where Percent(f) is a linear approximation of the relationship between the log of the frequency and the dB level (with the attenuation defining one of the points). Right now the Nan problem continues to persist, even after I loaded the coefficients. In Dataviewer, the channels look relatively normal for the past 10 minutes, as does the data when viewed with tdsdata.

 FIGURED IT OUT - THERE WAS A PROBLEM WITH THE LOW PASS FILTERS (TOO HIGH ORDER). FIXING IT NOW, SHOULD BE GOOD IN AN HOUR. 

Hi everybody,

Sorry for flooding the ELOG about the PEM channels. Today I

- Changed all of the GUR1 and GUR2 filters to elliptic, and lowered the orders of their low-pass filters.

- Lowered the order of the low-pass filters on the STS channels

- Changed the parameters in seismic.strip, which I saved as MashaTemplate2.

Attached is the most recent status of the channels as seen with StripTools:

I'm not currently sure how to apply my template to seismic.strip shown on the wall (I saved it as seismic.strip on Pianossa and copied the old file to seismic.stripOld). I understand the job is being run on Megatron. I'll play around with this later tomorrow. (In other words, the display currently on the wall, while it does not have the Nan spikes like yesterday and this morning does not currently display the template I made).

The RMS signals generated by the updated filtering process are now on the wall. The NaN issue is gone it seems, and the template has been changed. Thanks for your help, Jenne. 

Today I REPOSITIONED THE SEISMOMETERS in order to triangulate noise sources (as Rana suggested).

I re-levelled all of them, locked them, and turned them on. They should be located out of sight, but just in case:

GUR 1 IS DOWN THE X-ARM, behind the interferometer.

GUR 2 IS BETWEEN THE TWO ARMS, BEHIND THE CABLE TRAP THAT RUNS PARALLEL TO THE X-ARM.

STS 1 IS DOWN THE Y-ARM behind the interferometer.

I'll wait a day for them to stabilize (continuing to reset STS-1 every hour or so) and then begin taking data tomorrow morning, depending on the condition of the signal.

Ideally, I'd like a few days' worth of data, so I'll update when I've changed the configuration back to the way it was prior.

This week, my work fell into two categories: Artificial Neural Networks and lab-related projects.

Artificial Neural Networks

- I played around with radial basis functions and k-means classification algorithms for a bit in order to develop an algorithm to pick out various features of seismic signals. However, I soon realized that k-means is an extremely slow algorithm in practice, and that radial basis functions are thus difficult to implement since their centers are chosen by the k-means algorithm in practice.

- Thus, I moved on to artificial neural networks. Specifically, I chose to implement a sigmoidal neural network, where the activation function of each neuron is f(u) = 1/ (1 + e-u/T), T constant, which is nice because it's bounded in [0, 1]. Classification, then, is achieved by generating a final output vector from the output layer of the form [c1, c2, c3, ..., cN] where N is the number of classes, ci = 1 (ideally) if the input is of class i, and ck = 0 otherwise.

- First, I built a network with randomly generated weights, ten neurons in the one hidden layer, and two output neurons - to simply classify [1, 0] (earthquake) and [0, 1] (not an earthquake). I ran this on fake input I generated myself, and it quickly converged to error 0. Thus, I decided to built a network for real data.

- My current network is a 2-layer, 10 neuron / 2 neuron sigmoidal network that also classified earthquake / not an earthquake. It trains in roughly 80 - 100 iterations (it's learning curve on training data it attached). It decimates full data from DataViewer by a factor of 256 in order to run faster.

- Next steps: currently, my greatest limitation is data - I can use US Geological Survey statistics to classify each earthquake (so that N = 10, rather than 2, for example), but I would like definite training data on people, cars, trucks, etc. for supervised learning, in order to develop those classes. Currently, however, the seismometers are being used for mine and Yaakov's triangulation project, so this may have to wait a few days.

Lab-Related Projects

- I apologize for all of the E-logs, but I changed the filters in the RMS system (to elliptic and butterworth filters) and changed the seismic.strip display file.

- I repositioned the seismometers so that Yaakov and I can triangulate signals and determine seismic noise hot-spots (as a side-project).

Right now I'm going to try for more classes based on USGS statistics, and I will also explore other data sources Den suggested.

Thanks for your help, everybody in 40m!

I noticed on DataViewer today that GUR2 was outputting only noise (somewhere around 2 counts). Jenne suggested that GUR 2 might not be plugged in. I turned off the ADC, and tried several times to plug GUR 2 back in. I thought something might be wrong with the cable, but when I plugged the GUR1 cable into GUR2, there was still no readout (although the GUR1 cable works fine when I plug it into GUR1). Perhaps I'm just inept at plugging in GUR2, or perhaps there's another issue. Either way, I'll ask Jenne about it tomorrow and try again.

It seems that the STS-1 ADC channels had the same mismatch issue as the GUR-2 channels. The PEM_MONITOR has STS_1 listed as channels 6, 7, 8 (+1 on the actual ADC) whereas it was plugged into channels 13, 14, 15 (+1 on the actual ADC as well) with nothing in channels 6, 7, 8. Thus, I moved the cables and reset STS_1. The readout, however, is still only a magnitude of ~10 counts (I checked, however, that this is indeed the readout when the seismometer is plugged in vs. when it is unplugged), but hopefully it will stabilize during the day, as did GUR 2.

Den and investigated the STS-1 problem (which is currently plugged into ADC channels 13, 14, and 15, which correspond to the STS-3 channels in dtt). It turns out that I had plugged in the power to the monitor in the host box rather than the remote. The X, Y, and Z readout is currently approaching a mean of zero, and I will let it continue to do so overnight (pressing auto-zero as necessary). Attached is a plot of the coherence with GUR 1, and the time-domain signals.

The Streckeisen is currently plugged into ADC channels 13, 14, and 15 (corresponding to STS-3 in the channels). The X, Y, and Z components are correct. The signals is zeroed (it's been so for at least the past 10 hours), the coherence with GUR1 looks decent, and the signal looks similar to the GUR1 signal.

This week, I continued to work with my Artificial Neural Network. Specifically, I implemented a 3-hidden layer sigmoidal, gradient-descent supervised network, with 3 neurons in the final output layer, since I have introduced a new class, trucks. I have overcome my past data limitation, since I observed that there is a multitude of trucks that comes between 9 and 10 am, and thus I have observed a bunch of trucks after the fact (their seismic patterns are rather distinct, and thus there could prove to be a very large supply of this data - I have gathered data on the past 50 days so far, and have gathered 60+ truck patterns).

With 3 classes, the two-layer network converges in ~200 epochs, while the 3-layer network takes around ~1200 (and more time per iteration). Since the error gradients in the stochastic gradient descent are recursively calculated, the only real time limitation in the algorithm is just lots of multiplication of weight / input vectors, lots of computation of sigmoidal functions, and lots of data I/O (actually, since the sigmoidal function is technically an exponentiation to a decimal power and a division, I would be curious to know if theory or Matlab has any clever ways of computing this faster that can be easily implemented - I will look into this today). Thus, the networks take a long time to train. I'm currently looking at optimizing the number of layers / number of neurons, but this will be a background process for at least several days, if not the next week. In the greater scope of things, however, training time isn't really a problem, since the actual running of the algorithm requires only one pass through the network, and the network should be as well-trained as possible. However, due to the fact that I am only here until the end of August, it would be nice to speed things up.

As far as other classifications, I can simulate signals either by dropping the copper block from the Stacis experiment, or by applying transfer functions to general seismic noise. However, I would like more real data on noise sources, but the only other one plausible to LIGO that I can currently think of (cars don't show up very well) is the LA Metro. Perhaps I will take a day to clock trains as they come in (since the schedule is imprecise) and see if there is any visible seismic pattern.

I also, with the help of Yaakov, Jenne, and Den, now have three working, triangulated seismometers, which can now begin taking triangulation data (the rock tumblers are still working, so there should be opportunities to do this), both to find hot-spots as Rana suggested, and to measure the velocity and test out my algorithm, as Den suggested.

I realized what the ADC channel mismatch was, and apologize for plotting a terrible coherence in log scale. The channels are now properly matched (there is decent coherence between GUR1_X/STS_X, etc.).

Today I worked with the BLRMS channels, re-triangulated the seismometers (the STS is now on the very end of the Y-arm, while the GUR2 is on the X-arm - this GUR2 cable will need to be either extended or replaced - Jenne and I will look at parts tomorrow), and added 0.01 - 0.03 Hz and 0.03 - 0.1 Hz RMS channels (However, the MEDM files for these are not yet complete - I will finish these tomorrow) in order to be able to better see earthquakes. I also did some things for the neural network project, including beginning Simulink tutorials so that I can run my code by applying a force on a damped harmonic oscillator + white noise until it stops.

I will explain the methodology behind the new RMS filters tomorrow morning, when the seismometers have settled down and I can make coherence plots.

I'll post a better E-log tomorrow when it's not 2 in the morning.

As Jenne suggested last night, I changed the C1PEM overview in Epics. Previously, the C1PEM_OVERVIEW.adl screen had two separate visualizations, one for LP and one for BP for each channel. I changed the format so that there is only one frame per RMS channel, showing all of the input and output as before, but continuously, to demonstrate the actual RMS process.

First, the input is bandpass filtered, then multiplied by itself (squared), then lowpass filtered, and then square rooted to yield the final output. I have kept the previous files, but have applied this new format to all of the RMS ACC*, GUR1, GUR2, and STS1 channels. 

As Rana suggested yesterday, I made channels for 0.01 - 0.03 Hz and 0.03 - 0.1 Hz, and created filters for them, which I will work on some more now.

RMS Filters

How Den, Rana, and I chose RMS filters:

- Because filter ring-down generates negative outputs, which then show up as NaN when the log is taken in StripTool, we decided to only use low-pass filters with real poles (using ZPK in Foton).

- The band-pass filters were chosen by looking at the dB drop from the cutoff frequencies to the next (usually aiming for 40 dB, or 99%), and checking that the BP_IN and BP_OUT had a coherence of 1 in the pass-band.

- The low-pass filters were chosen by finding the lowest filter order at which there was coherence of ~1 in the passband between the input signal to the filter and the filter output. The cutoff frequencies were chosen to be lower than the first passband frequency, in order to get rid of the cos(2*\omega*t)/2 terms that arise during the squaring of a signal of the form Asin(\omega*t) and to assure that only terms related to A^2/2 were kept. A plot of the 3-10 region is attached - in the Coherence plot, the coherence in ~1 in the 3 to 10 hZ region. Likewise,

PEM Sitemap

My previous post had digital zeros in two of the BLRMS channels. Jenne figured out that this was because the file Csqrt.c, which performs the square root operation in the root-mean square processing only accepted 5 inputs. I modified and committed the code so that it now accepts 7 inputs (for our 7 frequency bands) and returns 7 outputs. The new PEM sitemap seems to currently work.

StripTool

I have modified the StripTool file in order to show our new 0.01 Hz - 0.03 Hz and 0.03 Hz - 0.1 Hz channels.

This week, I have been mainly working on my new neural network project, and working with the BLRMS channels.

New Neural Network Project:

Rana and I read a paper which demonstrated the use of a neural network to control an inverted pendulum. Since the test masses are, in simplest terms coupled harmonic oscillators with six degrees of freedom, we suspect something similar can be done. The logic behind trying to use a neural network as a control system lies in the fact that the actual motion of the test masses is a complicated function of many variables (including environment variables, such as temperature and seismic disturbance as well), and neural networks, at the core, are function approximators, which can approximate a function as a combination of polynomials (with error 0 in a closed interval) or of sigmoidal functions (which, in the papers I've read, seems to be possible in practice - the authors also reference some theoretical results regarding this). Thus, this network could, in theory take the displacement in the degrees of freedom of the OSEMS and generate the actuation force.

However, there was the question of building a neural network which explicitly took the history of a time series into account, and not only implicitly in the neuron weights. Thus, I found recurrent neural networks, which, for some T (this will ultimately be related to the sampling frequency, or spectrum of the signal),  has T hidden sub-layers which pass in the output of the T previous iterations.

I have written Matlab code that generates, for a given T and neuron vector (the number of neurons in each of the T layers) a recurrent neural network (it still, like my old network, uses gradient descent and sigmoidal activation functions). However, I needed a system to run it on, so I began playing around with Simulink. I have never used this system before, so I took a bit of time doing the tutorials, but I currently have a damped, driven harmonic oscillator which accepts a force at each time step and outputs a displacement (taking into account the earlier forces). This afternoon, I think I will link it to my recurrent neural network, which can accept a displacement and output a force, and see what happens. These networks are notoriously slow to converge, but I could at least check to see that my code (and understanding of the algorithm) is correct. Should this work, I can also try hooking it up to Sasha's simulation, which is more complex than an oscillator with one degree of freedom.

As Far as Seismic Things Go:

I haven't worked with seismic noise itself very much this week, although the signals could be used as inputs to a neural network. I moved the STS1 and GUR1 cables, which I realized were crossed and thus stupidly placed (my fault) and now have STS1 on the very end of the Y arm (GUR1 is about 10 meters down the X arm, however, so Jenne and I will probably have to make a cable). I also found granite and foam, so at least the STS can be given a final position at the end of the arm.

BLRMS:

I added two new channels to the BRLMS, for 0.01 - 0.03 Hz and 0.03 - 0.1 Hz, so that we can better observe distant earthquakes. I also changed the MEDM screen for all of the ACC, GUR, and STS RMS channels (I think Liz is also using this for her MIC channels), combining the low pass and high pass filter screens into a screen demonstrating the complete RMS process, which squaring and square roots. I have also changed the RMS filters, and wrote an Elog about that process.

The paper Rana and I read is Charles W. Anderson, Learning to Control and Inverted Pendulum Using Neural Networks

A good paper of recurrent neural network basics is Mikael Boden, "A guide to recurrent neural networks and backpropagation"

I checked the connections specified in the old Gulap Pin Map and found that they do not correspond to the current values. I mapped out the current connections (in this case, the letter refers to the labeled pin on the mil/spec while the number refers to the pin on the 37 pin DSub, labeled consecutively):

A-1, B-2, C-3, D-4, E-5, F-6, G-7, H-Unused, J-8, K-unused, L-9, M-10, N -11, P-12, S-13, T-Unused, U-14, V-15, W-16, X-17, Y-18, Z-Unused, a-Unused, b-19, c-20, UnlabeledPin-Unused.

There are 20 pins in use of 26 total, which is good because that means Jenne and I can use the ~70m long 24 wire cable to make a new Gurlap 1 cable.

 A lot of my time this week was spent struggling with implementing my neural network code in Simulink in order to experiment with using neural network control systems, as Rana suggested. Perhaps I'm inept at S-Functions, but I decided try to use the Model Reference Controller block in Simulink instead of my own code, and experiment with using that to control a driven, damped harmonic oscillator with noise. The block consists of two neural networks, one which is trained on a model of the plant to simulate the plant, and one that it trained on a reference model (which, in my case, is just input -> gain 0 -> output since my desired signal for now is 0. So far, I have managed to adjust the parameters enough to stop the neural network controller from outputting too much force (this causing the amplitude of the oscillator to increase with each iteration), but outputting enough to keep the plant oscillating with a maximal displacement of 2 m/s (with 30 neurons, 100 reference time delays, and 100 plant output time delays). I will continue to work on this, especially with added noise sources, and see how feasible it is to train such a controller to actually perform well.

As far as classification, I took up that project again, and realized that I have been approaching it the wrong way the whole time - instead of using a neural network to classify an entire time series, I could just classify it online using a recurrent neural network. This, however, meant that my data (which used to be in packets of 900 seconds) had to be parsed through to generate time-dependent desired output vectors. I did this last night, and have been trying various combinations of neuron numbers, time delays, and learning parameters this morning. Below is my current best result for mean square error with time, obtained with 1 neuron in the hidden layer, 50 time delays (so that there are actually 51 neurons feeding into the hidden layer, and a subnetwork of 50 neurons connected to 50 neurons, and learning parameter 0.7). The peak is due to the fact that a large amount of sharp earthquakes occur around that time, essentially giving the neural network a surprise, and causing it to rapidly learn. However, I suspect this sharp rise would decrease if I were to stop decimating my data by a factor of 256, and using all of the inputs as they come in (this, however would be drastically slower). Currently, I have a massive loop running which tries different combinations of neurons and time delays.

In terms of other lab stuff, Jenne and I ordered parts to make a cable for Gurlap-1, and I updated the pin map for Gurlap-1. Also, I wrote my progress report.

The parts Jenne and I ordered arrived today, so we made a long cable for Guralp 1 using a 24 + 1 wire 70 meter long cable, a female 37-pin DSub, and a 26-pin milspec. The pin map is the same as the one I specified in my previous E-log. I soldered both the milspec attachment and the DSub attachment, and used a Multimeter to check the connectivity of the cables. 20 of 20 connections worked (beeped), so I plugged  the cable into the Gurlap 1 seismometer and the Guralp box.

The time series comparison for the two cables

Old cable:

New cable: (I had to move GUR 1, so it's still stabilizing in the X and Y time series)

 New

The current signal spectrum

The BLRMS on the seismic strip also look similar using the two cables - it's more visible on the wall, but I will include a StripTool picture:

New Cable BLRMS (similar to old cable BLRMS)

Den and I moved the Streckeisen, Guralp 2, and Trillium seismometers to the isolation box in order to measure the noise of the Streckeisen while we have the Trillium.

[Sasha, Masha, Liz, Eric]

A bunch of surfs in the lab just noticed that ETMX is going crazy (laser is shifting everywhere) due to a 4.5 EQ that just hit LA. The optic is already shut down according to the watchdogs.

The main thing that I did this week was write a C block that, given static weights, would classify seismic signals (into one of three categories - Earthquake, Truck, Quiet). I have successfully debugged the C block so that it works without segmentation faults, etc, and have made various versions - one that uses a recurrent neural network, and one that uses a time-delayed input vector (rather than keeping recurrent weights). I've timed my code, and it works very fast (to the point where clock() differences in <time.h> are 0.000000 seconds per iteration). This is good news, because it means that operations can be performed in real-time, whether we are sampling at 2048 Hz, or, as Rana suggested, 256 Hz (currently, my weights are for 256 Hz, and I can decimate the incoming signal which is at 2048 Hz right now).

In order to optimize my code, since at the core it involves a lot of matrix multiplications, I considered how the data is actually stored in the computer, and attempted to minimize pointer movement. Suppose you have an array in C of the form A[num_row][num_col] - the way this array is actually stored on the stack or heap is row_1 / row_2 / row_3 / ... / row_num_row, so it makes sense to move across a matrix from left to right and then down (as though reading on a page). Likewise, there's no efficient algorithm for matrix multiplication which is less that O(N^2) (I think), so it's essentially impossible to avoid double for loops (however, the way I process the matricies, as mentioned before, minimizes this time).

The code is also fast because, rather than using an actual e^-u operation for the sigmoidal activation function, it uses a parametrized hyperbola - this arithmetic operations are the only ones that occur, and this is much faster than exponentiation (which I believe is just computer by Taylor series in the C math library..)

The weight vectors for this block are static (since they're made from training data where the signal type is already known). I am not currently satisfied with the performance of the block on data read from a file, so I am retraining my network. I realized that what is currently happening is that, given a time-dependent desired output vector, the network first trains to output a "quiet" vector, then a "disturbance" vector, and then retrains again to output a "quiet vector" and completely forgets how to classify disturbance. Currently, I am trying to get around this problem by shifting my earthquake data time-series, so that when I train in batch (on all of my data), there is probably an earthquake at all time points, so that the network does not only train on "quiet" at certain iterations. Likewise, I realized that I should perform several epochs (not just one) on the data - I tried this last night, and training performance MSE decreased by a factor of 1 per iteration (when on average, it's about 40, and 20 at best).

After I input the static weight vectors (which shouldn't take long since my code is very generalized), the C block can be added to the c1pem frame, and a channel can be made for the seismic disturbance class. I've made sure to keep with all of the C block rules when writing my code (both in terms of function input/output, and in terms of not using any C libraries).

As for neural networks for control, I talked to Denis about the controller block, and he realized that we should, instead of adding noises, at first attempt to use a reference plant with a lower Q pendulum and a real plant with a higher Q pendulum (since we want pendulum motion to be damped). I've tried training the controller block several times, but each time so far the plant pendulum has started oscillating greatly. My current guess at fixing this is training more.

Also, Jenne and I made a cable for Guralp 1 (I soldered, she explained how to do it), and it seems to work well, as mentioned in my previous E-log. Hopefully it can be used to permanently keep the seismometer all the way down the arm.

The c1lsc and c1sus screens are red in the front-end status. I restarted the frame builder, and hit the "global diag reset" button, but to no avail. Yesterday, the only thing Den and I did to c1sus was install a new c1pem model. I got rid of the changes and switched to the old one (I ran rtcds build, install, restart), but the status is still the same.

[Masha, Sasha]

Sorry to spam the e-log, but did someone come knock loudly on the emergency exit door a few moments ago? It gave Sasha and I quite a fright, and we are rather worried.

Den and I decided to try to classify seismic signals in the frequency domain rather than the time domain. We looked at amplitude spectral density plots of all of the data in our set, and noted that there were noticeable differences in the frequency domain for midnight quiet, trucks, and earthquakes.

For example, here is the time series of quiet, midnight seismic noise as compared to the seismic noise at the peak of an earthquake - the earthquake signal is noticeably higher in the 1 - 3 Hz region. Likewise, for the truck signal, there are noticeable bumps that arise at 10 and 30 Hz during the peak of the truck's motion due to the resonant frequency of the truck bouncing on its wheels.

We investigated this potential means of classification further by considering the linear separability of the power of our signals in various frequency bands. Below is a plot of the power of a normalized signal in the 0.1 - 3.0 Hz region vs. the power of the normalized signal in the 3.0 - 30.0 Hz region - calculated by means of fft and separation of the discrete resulting frequencies (in short, an ideal filter).

There is rather clear linear separability of the normalized signals in this case, as two lines could potentially be drawn to separate trucks from quiet and earthquake in this case (with a few misclassified points due to quiet - since the lab isn't actually empty and quiet in the middle of the night, and man-made seismic disturbances to occur). The reason we have to normalize our signals lies in the fact that the data set had different gains for various seismometers at different times. Normalization not only allows us to use our data set for training effectively, but it also assures that the online classification, if the online signals are also normalized, will allow for variable seismometer gains in the future and still be able to classify signals.

I looked at the linear separability of our training set using various combinations of frequency bands, and deduced that the current separation in the BLRMS preforms best (coincidentally, since the BLRMS separations are just decades), which meant that we could use the current BLRMS system we have for online classification of seismic noise.

Thus, I built a neural network which performed classification with the following parameters:

- One hidden layer of 20 neurons

- Gradient descent backpropagation with learning parameter mu = 0.175

- Sigmoidal activation functions for each neuron (computationally achieved by a parametrized hyperbola rather than an actual hyper-tangent in order to save on computation time). 

- 5 inputs - the normalized fft^2 of the signal (since the root of a signal doesn't add linearly to 1) in the following frequency regions: 0.1 - 0.3, 0.3 - 1.0, 1.0 - 3.0, 3.0 - 10.0 and 10.0 - 30.0 Hz. Since this division was done through the (frequency, fft value) return in Matlab, the signal was essentially filtered ideally into these frequency bands.

- 3 output neurons representing an output vector, with desired output vectors of [1, 0, 0] for earthquake, [0, 1, 0] for truck, and [0, 0, 1] for quiet.

- 1,600,000 training epochs (batch backpropagation on all of the data)

Below is the best learning curve for this network, representing the total amount of inputs misclassified out of 224. The best result achieved was 30 misclassified signals out of 224. Obviously this is not ideal, but our data is not totally linearly separable. This could, however, be reduced with further iterations, but given the close to 0 slope of the learning curve between iteration number 1,000,000 and number 1,500,000, this could take a very long time.

Thus, I trained the network, generated the weight vectors and optimal activation function parameters, and was ready to implement a feed-forward neural network (with no online training). My next e-log (Part 2) will be about this system and will be posted shortly.

As promised in previous e-log, this log is all about the current online seismic noise classification system.

While we had the BLRMS system already in place (which I helped make), Den realized that we would need better filters for the BLRMS channels, as we wanted a strong cut-off, but we also wanted a short step-response so that we could quickly classify seismic signals. Likewise, having a step response which oscillates is also undesirable as this could lead to false classifications of post-truck signal as trucks as a filter adjusts and then dips back down. Thus, after experimenting with many different filters, Den chose to use a combination of

chebyl("LowPass", 1, 1, 0.03)*chebyl("LowPass", 1, 1, 0.03)

as our low-pass filter. The step response and bode plot are below.

 The next step was to write C code that would implement the feedforward neural network with my newly generated weights.

Next, I had to implement the code in the c1pem model, and normalize the inputs. Below is an overview of the model, and a close up of the C block section.

The above close-up includes the process of normalization (dividing by the square of the incoming signal), feeding through the neural network, and classifying.

Each seismometer channel set (GUR1X, GUR1Y, GUR1Z, GUR2X, GUR2Y, GUR2Z, STS1X, STS1Y, STS1Z) now has channels (and corresponding DQ channels) of the following form:

SEIS_CLASS : The class of seismic noise 1.0 means Earthquake, 0.5 means Quiet, and 0.0 means Truck. (There are only these 3 digital values).

SEIS_CLASS_EQ, SEIS_CLASS_TRUCK, SEIS_CLASS_QUIET: These channels represent the confidence of the neural network's classification. The class of the current signal will have an output of 1, where the other two channels will have an output between 0 and 1 representing the ratio of the neural network's output in that class neuron to the output in the classification vector neuron. To simply - suppose the neural network classified an earthquake. Ideally, the neural network output neurons would have the value [1, 0, 0], and SEIS_CLASS would equal 1.0 for earthquake. However, the output neurons probably read something along the lines of [0.9, 0.3, 0.5] - SEIS_CLASS is still 1.0, but SEIS_CLASS_EQ would be 1.0, and SEIS_CLASS_TRUCK would be 0.5 / 0.9 and SEIS_CLASS_QUIET would be 0.3 / 0.9. The lower the other two signals are, the better - this means that we are more confident in our classification.

The MEDM screen for this system (in the RMS system) has the following form for all seismometer channels (this one is GUR1X):

These are the screens I edited earlier in the summer, with modifications. The bottom filter banks represent the norm of the seismometer signal, which we use to normalize the inputs to the neural network.

Here a close-up of the most important part:

The orange meter on the right points to the current signal type. Here it reads truck - this is ok because it's the middle of the day, and there are a lot of trucks around. The left side represents our confidence in the signal - the signal is classified as a truck, so the "Truck" bar is saturated. The quiet signal bar is very low, which is good since it means that the neural network thinks that it's definitely not quiet. The earthquake bar has some magnitude, since earthquake signals and trucks have some degree of linear non-separability.

How has this been performing? Firstly, all of the seismometer channels have the same classification readout, which is good. Last night, all of the classes were "quiet", with an "earthquake" which occurred when Den jumped around GUR1 to simulate an EQ. This morning it was on "truck" as expected. The filters are still not fine enough to detect individual trucks, but I will continue to monitor the performance over the coming days.

If anyone has ideas on how better to represent this information, please let me know. This was the first thing that came into my head that would work with my MEDM monitor options, and I'm open to suggestions!

I added a widget to the C1PEM_OVERVIEW MEDM screen. The screen shows the nine seismometer channels (GUR1, GUR2, and STS1 X, Y, and Z), the current signal class in dark red, and the overall confidence in the classification, as Rana suggested. The confidence indication thresholds range from 0.1 to 0.9, in intervals of 0.1. Basically, if a signal class is completely dark red, and the other two classes show only white, or, better yet, nothing at all, this means that we have a clear classification. If, however, the other regions have some yellow, or even red indicators, this means that we are not very confident in our signal classification.

This is a screenshot of the widget. The nine seismometer channels are classifying the signal as quiet, which is good both because it's the middle of the night, and because the nine seismometer signals somehow agree (I'd use the word correspond with one another, but that implies a strong level of coherence..). The confidence is high, seeing as there's little indication in the truck and earthquake regions (none whatsoever in the truck, meaning that the signal, given our classification method, could not possibly be a truck, and some in the earthquake region (below 0.1 of the quiet signal classification strength, however), possibly due to low seismic disturbance).

Since the classification finally works (or seems to work..), I wrote triangulation scripts in Python which triangulate the signals, and a plotting script in Matlab which generates a heat map of seismic noise source locations. I switched the ADC Streckeisen and Trillium connections in order to better triangulate with the current channels, and will return them either tomorrow, or when I come back from Livingston so that we can have weekday data as well.

There was a 5.6 Earthquake that occurred near Tofino, Canada about 30 minutes ago. It showed up rather strongly on the BLRMS.

The neural network classification system also picked up on it, but oscillated from Earthquake (1.0) to Quiet (0.5) perhaps due to the filters we currently have installed. Here is a shot of the GUR1X classification channel at the time of the EQ:

I cleaned up my directory (/users/masha) today. A lot of the files are just code that I experimented with, but the important files for training the classification neural network are in "neural_network_classification". The "EarthquakeData" subdirectory contains my entire dataset. Files of the form "GenerateRNNInput" are used to create input vector sets to the network, while files of the "*NeuralNetworkClassification* actually run the code that generates the neural network vectors for the classification code block in the c1pem model.

Also, the folder "feed_c", which can also be found in Den's directory, contains the neural network controller code we played around with.

I recently realized that I may have over-trained my classification neural network and used too many parameters, so that my weight vectors are too fine-tuned to my particular data set and do not generalize well. I lowered the number of hidden neurons in the network to 15, and the number of epochs to 25000, and regularized based on the deltas (the gradient). Here is the most recent learning curve:

The old weights and code are saved in the c1pem directory in the file "classify_seismic_20neurons.c", while the current 15 neuron network is saved as "classify_seismic.c". I'll monitor the performance of this current network throughout the day, and decide which one we should keep.

 Hi everyone,

I've been working a bit on neural network code for a controller, and thus far I have code that creates a reference plant neural network. This is necessary to perform a gradient-descent learning algorithm on a controller neural network (one that reads an error signal and outputs actuation force). Because the error signal is read after the previous output of the controller neural network has passed through the plant, in order to calculate the gradient, either the inverse of the plant needs to be calculated, or the plant can be simulated through a neural network, and the error signal can thus be back-propagated through the plant neural network to find the gradient with respect to the output (as opposed to with respect to the plant), and thus back-propagated through the controller network in order to learn. 

I have uploaded to my directory a directory neural_plant. The most important file is reference_plant.c, which compiles with the command

 gcc reference_plant.c -o reference_plant -lfann -lm

The code runs on a file called reference_plant.data, which consists of a series of delayed inputs i_1, i_2, i_3 ... i_{n - 1} of the plant signals and then an output that is i_{n}, the subsequent signal. Parallel streams may also be used, if more than one signal is to be read. The top of the file must contain the number of total training packets (input-output groups), followed by the number of inputs, and the number of outputs. reference_plant.c also has constant variables which specify the number of hidden neurons, which can be changed. 

Hi everyone,

I've been on break this week, so in addition to working at my lab here, I've done some NN stuff. In response to Den's response to my last post, I've included learning curve plotting capabilities, 

I've explored all of the currently documented capabilities of FANN (Fast Artificial Neural Network - it's a C library) (most likely, there are additions to the library floating around in open-source communities, but I have yet to look into those). There is extensive FANN documentation on the FANN website (http://leenissen.dk/fann/html/files/fann-h.html), but I'll cut it down to the basics here:

FANN Neural Network Architectures

standard: This creates a fully connected network, useful for small networks, as in the reference plant case 

sparse: This creates a sparsely connected network (not all of the connections between all neurons exist at all times), useful for large networks, but not useful in the reference plant case, since the number of neurons is relatively small

shortcut: This creates some connections in the network which skip over various hidden layers. Not useful in the harmonic oscillator case since there are no hidden layers. Probably won't be useful in a better-modeled referrence plant since this reduces the non-linear capabilities of the model.

FANN Training  

TRAIN_INCREMENTAL: updates the weights after every iteration, rather than after each epoch. This is faster than the other algorithms for the reference plant.

TRAIN_BATCH: updates the weights after training on the whole set. This should not be used on batches of data for the reference plant, seeing as the time history dependence of the plant is smaller than the size of the entire data set. 

TRAIN_RPROP: batch training algorithm which updates the learning parameter.

TRAIN_QUICKPROP: updates the learning parameter, and uses second derivative information, instead of just first derivative, for backpropagation. 

FANN Activation Functions

FANN offers a bunch of activation functions, including a function FANN_ELLIOT, which is essentially the "signmoid like" activation function Den and I used this summer, which runs in the order of multiplication and addition. The function parameters (steepness) can also be set.

FANN Parameters

As usual, the learning parameter can be set. While over the summer we worked with lower learning parameters, in the case of the harmonic oscillator reference plant, since the error is low after the first iteration, higher learning parameters (0.9, for example), work better. However, this is a very isolated case, and, in general, lower parameters, though convergence is slower, produce more optimal results. 

The learning momentum is another parameter that can be set - the momentum factor is a coefficient in the weight adjustment equation which allows for the difference in weights beyond the previous weight to be factored in. In the case of the reference plant, a higher learning momentum (0.9) is optimal, although in most cases, a lower learning momentum is optimal so that the learning curve doesn't oscillate terribly. 

FANN does not explicitly include regularization, but this can be implemented by checking the MSE  at each iteration against the MSE at the n previous iterations, where n is the regularization parameter, and stopping training if there is no significant decrease (also determined by a parameter). The error bound I specified during training was 0.0001

The best result for the reference plant was obtained using FANN_TRAIN_INCREMENTAL, a "standard" architecture, a learning rate of 0.9 (as explained above) and a learning momentum of 0.9 (these values should NOT be used for highly non-linear and more complicated systems). 

I have included plots of the learning curves - each title includes the architecture, the learning algorithm, the learning parameter, and the learning momentum if I modified it explicitly.

All of my code (and more plots!) can be found in /users/masha/neural_plant

On the whole, FANN has rather limited capabilities, especially in terms of learning algorithms, where it only has 4 (+ all of the changes one can make to parameters and rates). It is, however, much more intuitive to code with and faster han the Matlab NN library, although the later has more algorithms. I'll browse around for more open-source packages. 

Best,

Masha

Yaakov and I investigated the GUR 2 problem. It turns out that the ADC channels that GUR 2 was plugged into, ADC channels 6 through 8 (on the actual ADC they are C7 through C9), did not correspond to the channels labelled "GUR 2" in the PEM, ADC channels 3 through 5. We modified them so that GUR 2 is now plugged into ADC channels 3 through 5 (on the ADC it's +1).

Before we discovered that this was the problem, we attempted to take the cover off of GUR 1 to check the gains, and discovered a stripped Allen screw on the side by the "Vertical" pot, which we removed.

Now the GUR 2 readout looks good, and we will give it more time to settle down before we take data.

I mounted the new projector to the pipe where the old projector was attached. The mounting hardware wasn't designed for attaching to a pipe, but with Steve's help I mounted the projector. The projector is angled away from the area of the wall designated as the screen, and I am going to meet with Steve Monday to fix this.

I finished mounting the new projector. The projector and computer monitor now display information.