#+TITLE: LFP Detection: Common Lisp Version
#+AUTHOR: Christophe Pouzat
#+EMAIL: christophe.pouzat@gmail.com
#+DATE: 2011-10-28 ven.
#+DESCRIPTION:
#+KEYWORDS:
#+LANGUAGE: en
#+OPTIONS: H:3 num:t toc:t \n:nil @:t ::t |:t ^:t -:t f:t *:t <:t
#+OPTIONS: TeX:t LaTeX:t skip:nil d:nil todo:t pri:nil tags:not-in-toc
#+INFOJS_OPT: view:nil toc:nil ltoc:t mouse:underline buttons:0 path:http://orgmode.org/org-info.js
#+EXPORT_SELECT_TAGS: export
#+EXPORT_EXCLUDE_TAGS: noexport
#+LINK_UP:
#+LINK_HOME:
#+XSLT:
#+STYLE:
#+STYLE:
#+STYLE:
#+STYLE:
* Define a small =CLisp= utility function
As an =R= (or =matlab=/=octave=) user one quickly gets used to functions like =ls= (or =who=) to visualize the variables and functions defined in a session. Since =CLisp= does not have such a function, we start by defining it:
#+begin_src lisp
(defun own-symbols ())
(setf (symbol-function 'own-symbols)
(let ((initial-lst ()))
(do-symbols (sym 'cl-user)
(when (or (boundp sym) (fboundp sym))
(push sym initial-lst)))
#'(lambda ()
(let ((new-lst ()))
(do-symbols (sym 'cl-user)
(when (or (boundp sym) (fboundp sym)) (push sym new-lst)))
(set-difference new-lst initial-lst)))))
#+end_src
#+RESULTS:
: #
* Download data with =wget=
The data recorded from 4 electrodes result from a preprocessing briefly described in the main text and are stored as signed integer coded on 4 Bytes. They must be mutiplied by 0.12715626 on channels 1 and 4 and by 0.01271439 on channels 2 and 3 to get voltage in \mu V. They are sampled at 500 Hz and 600 second are stored.
We are going to download the data and save them on our hard disk with the command line program [[http://www.gnu.org/s/wget/][Wget]]. The data are located in a repository whose URL address is: =http://www.biomedicale.univ-paris5.fr/physcerv/C_Pouzat/Data_folder/=. There are four files named =SpinalCordMouseEmbryo_CHX.dat=, where =X= equals 1, 2, 3 and 4.
We start by creating some variables in our Common Lisp environment containing parts of the names of the data files.
#+name: repository-name
#+begin_src lisp
(defparameter repository-name
"http://xtof.disque.math.cnrs.fr/data/"
"URL where data can be found")
#+end_src
#+RESULTS: repository-name
: REPOSITORY-NAME
#+name: data-name-prefix
#+begin_src lisp
(defparameter data-name-prefix
"SpinalCordMouseEmbryo_CH"
"Prefix of the data file names")
#+end_src
#+RESULTS: data-name-prefix
: DATA-NAME-PREFIX
#+name: data-name-suffix
#+begin_src lisp
(defparameter data-name-suffix
".dat"
"Suffix of the data file names")
#+end_src
#+RESULTS: data-name-suffix
: DATA-NAME-SUFFIX
We build the data =URL= with:
#+name: data-urls
#+begin_src lisp
(defparameter data-urls
(mapcar #'(lambda (i) (concatenate
'string
repository-name
data-name-prefix
i
data-name-suffix))
'("1" "2" "3" "4"))
"A list with the URLs of the four data files")
#+end_src
#+RESULTS: data-urls
: DATA-URLS
And we will soon also need the data file names that we build with:
#+name: data-file-names
#+begin_src lisp
(defparameter data-file-names
(mapcar #'(lambda (i) (concatenate
'string
data-name-prefix
i
data-name-suffix))
'("1" "2" "3" "4"))
"A list with the names of the four data files")
#+end_src
#+RESULTS: data-file-names
: DATA-FILE-NAMES
We can use our =own-symbols= function in order to see the "bound symbols" (variables and functions) we have created so far:
#+begin_src lisp :exports both
(own-symbols)
#+end_src
#+RESULTS:
| REPOSITORY-NAME | DATA-NAME-PREFIX | DATA-URLS | DATA-FILE-NAMES | DATA-NAME-SUFFIX |
Fine we can now download the individual data files and quickly check them with [[http://www.gnuplot.info/][Gnuplot]]. Do that we define two utility functions:
#+name: single-data-url(i=0)
#+begin_src lisp
(nth i data-urls)
#+end_src
#+RESULTS: single-data-url
: http://xtof.disque.math.cnrs.fr/data/SpinalCordMouseEmbryo_CH1.dat
#+name: single-data-file(i=0)
#+begin_src lisp
(nth i data-file-names)
#+end_src
#+RESULTS: single-data-file
: SpinalCordMouseEmbryo_CH1.dat
We can download the first data file:
#+name: raw-data-1
#+headers: :exports none :var url=single-data-url(i=0)
#+begin_src sh :cache yes :file SpinalCordMouseEmbryo_CH1.dat
wget $url
#+end_src
#+RESULTS[c9ea9ab637e893ec2b6e8210e75a727934705218]: raw-data-1
[[file:SpinalCordMouseEmbryo_CH1.dat]]
A quick check using Gnuplot is readily obtained:
#+begin_src gnuplot :exports code :eval never
plot 'SpinalCordMouseEmbryo_CH1.dat' binary format='%int32' using 1 with lines
#+end_src
We keep going with the second trace:
#+name: raw-data-2
#+headers: :exports none :var url=single-data-url(i=1)
#+begin_src sh :cache yes :file SpinalCordMouseEmbryo_CH2.dat
wget $url
#+end_src
#+RESULTS[d23fba7eec4961f402004d9d0449237baa902aaf]: raw-data-2
[[file:SpinalCordMouseEmbryo_CH2.dat]]
#+begin_src gnuplot :exports code :eval never
plot 'SpinalCordMouseEmbryo_CH2.dat' binary format='%int32' using 1 with lines
#+end_src
Now the third:
#+name: raw-data-3
#+headers: :exports none :var url=single-data-url(i=2)
#+begin_src sh :cache yes :file SpinalCordMouseEmbryo_CH3.dat
wget $url
#+end_src
#+RESULTS[1f3cd8a4703225a14195495501e1dfd811b07188]: raw-data-3
[[file:SpinalCordMouseEmbryo_CH3.dat]]
#+begin_src gnuplot :exports code :eval never
plot 'SpinalCordMouseEmbryo_CH4.dat' binary format='%int32' using 1 with lines
#+end_src
And the fourth and last one:
#+name: raw-data-4
#+headers: :exports none :var url=single-data-url(i=3)
#+begin_src sh :cache yes :file SpinalCordMouseEmbryo_CH4.dat
wget $url
#+end_src
#+RESULTS[3788e18f9949824788b5b79dc4f8ca8f71e96507]: raw-data-4
[[file:SpinalCordMouseEmbryo_CH4.dat]]
#+begin_src gnuplot :exports code :eval never
plot 'SpinalCordMouseEmbryo_CH4.dat' binary format='%int32' using 1 with lines
#+end_src
* Read data into Lisp session
#+name: data-list
#+begin_src lisp
(defparameter data-list
(mapcar #'(lambda (x)
(declare (optimize (speed 3)))
(with-open-file (str x :direction :input
:element-type '(signed-byte 32))
(let* ((n (file-length str))
(data-v (make-array n :element-type 'fixnum)))
(dotimes (i n) (setf (aref data-v i) (read-byte str)))
data-v)))
data-file-names)
"A list of 1 dimensional arrays, each with the raw or derived data of a single electrode")
#+end_src
#+RESULTS: data-list
: DATA-LIST
* Compute the derivatives
Since the "derivatives" will be normalized by dividing them by their =SD=, we are just going to compute the second order difference of each trace (and not divide these differences by $2 \, \delta$).
#+name: datad-list
#+begin_src lisp :results silent
(defparameter datad-list
(mapcar #'(lambda (v)
(declare (optimize (speed 3)))
(let* ((n (length v))
(res (make-array n
:element-type 'fixnum
:initial-element 0)))
(dotimes (i (- n 2))
(setf (aref res (1+ i))
(- (aref v (+ 2 i))
(aref v i))))
res))
data-list)
"The list of second order differences of the raw data")
#+end_src
* Data derivative normalization
#+name: mean
#+begin_src lisp
(defun mean (v)
"Return mean value of vector v"
(declare (optimize (speed 3)))
(let ((n (length v)))
(do ((i 0 (1+ i))
(ac 0 (+ ac (aref v i))))
((>= i n) (/ ac n))
(declare (fixnum i) (fixnum ac)))))
#+end_src
#+RESULTS: mean
: MEAN
We can easily check that the mean value on each derived trace is essentially zero:
#+name: mean-derivative-values
#+begin_src lisp
(defparameter mean-derivative-values (mapcar #'mean datad-list))
#+end_src
#+RESULTS: mean-derivative-values
: MEAN-DERIVATIVE-VALUES
#+begin_src lisp
mean-derivative-values
#+end_src
#+RESULTS:
| 1/6250 | 159/50000 | 13/75000 | 19/75000 |
We define next a function returning the variance. This function takes an optional argument =mu= which defaults to =nil=. If unspecified the mean value is estimated first. The specified or computed mean value is then subtracted from each obervation.
#+name: variance
#+begin_src lisp
(defun variance (v &optional (mu nil))
"Return variance of vector v"
(declare (optimize (speed 3)))
(labels ((sq (x) (* x x)))
(let* ((n (length v))
(divisor (if (eq mu nil) (1- n) n))
(mu (if (eq mu nil)
(mean v)
mu)))
(do ((i 0 (1+ i))
(ac 0 (+ ac (sq (- (aref v i) mu)))))
((>= i (1- n)) (/ ac divisor))
(declare (fixnum i) (fixnum ac))))))
#+end_src
#+RESULTS: variance
: VARIANCE
We get the variance on each site:
#+name: variance-derivative-values
#+begin_src lisp
(defparameter variance-derivative-values
(mapcar #'(lambda (x) (variance x 0)) datad-list)
"The variance of the derivative on each recording site.")
#+end_src
#+RESULTS: variance-derivative-values
: VARIANCE-DERIVATIVE-VALUES
Out of =VARIANCE-DERIVATIVE-VALUES= we easily get the standard daviations:
#+name: sd-derivative-values
#+begin_src lisp
(defparameter sd-derivative-values
(mapcar #'sqrt variance-derivative-values)
"The standard deviation of the derivative on each recording site.")
#+end_src
#+RESULTS: sd-derivative-values
: SD-DERIVATIVE-VALUES
We define the threshold parameter:
#+name: factor
#+begin_src lisp
(defparameter factor 4 "The detection threshold on the normalized derivatives")
#+end_src
#+RESULTS: factor
: FACTOR
We define next the number of active electrodes required to detect a genuine event:
#+name: active-electrodes-number
#+begin_src lisp
(defparameter active-electrodes-number 3 "Number of detection threshold crossings for genuine events")
#+end_src
#+RESULTS: active-electrodes-number
: ACTIVE-ELECTRODES-NUMBER
* Events detection
We start by "rectifying" the four traces; more explicitely we force positive values to 0 and store the result in a new list of 1D array, =datadr-list=:
#+name: datadr-list
#+begin_src lisp
(defparameter datadr-list
(mapcar #'(lambda (v)
(declare (optimize (speed 3)))
(map-into (make-array (length v) :element-type 'fixnum)
#'(lambda (x)
(if (> x 0) 0 x))
v)
)
datad-list)
"Rectified version of datad-list")
#+end_src
#+RESULTS: datadr-list
: DATADR-LIST
We now take the square of each element and divide it by the corresponding variance and compare the ratio with the square of =factor=. The result is stored in a new list of 1D array, =datadrn-list=:
#+name:datadrn-list
#+begin_src lisp
(defparameter datadrn-list
(mapcar #'(lambda (v d)
(declare (optimize (speed 3)))
(let ((f-sq (* factor factor)))
(map-into
(make-array (length v)
:element-type 'fixnum
:initial-element 0)
#'(lambda (x)
(cond ((and (< x 0) (> (/ (* x x) d) f-sq)) 1)
(t 0)
))
v)))
datadr-list
variance-derivative-values))
#+end_src
#+RESULTS: datadrn-list
: DATADRN-LIST
#+name: detected-events
#+begin_src lisp
(defparameter detected-events
(map-into (make-array (length (car data-list))
:element-type 'fixnum
:initial-element 0)
#'+
(nth 0 datadrn-list)
(nth 1 datadrn-list)
(nth 2 datadrn-list)
(nth 3 datadrn-list)))
#+end_src
#+RESULTS: detected-events
: DETECTED-EVENTS
Create a list, =events-pos= with the indexes at which the value of =detected-events= is larger than or equal to =active-electrodes-number=:
#+name: events-pos
#+begin_src lisp
(defparameter events-pos
(let ((res nil))
(do ((i 0 (1+ i)))
((>= i (length detected-events)) (reverse res))
(if (>= (aref detected-events i) active-electrodes-number)
(push i res)))))
#+end_src
#+RESULTS: events-pos
: EVENTS-POS
We can see how many elements =events-pos= contains:
#+begin_src lisp
(length events-pos)
#+end_src
#+RESULTS:
: 25
Now we keep only events' times which are more than 100 ms, that is, 50 sampling points appart and store the result in a new list, =lfp-pos=:
#+name: minimal-interval
#+begin_src lisp
(defparameter minimal-interval 50
"The minimal number of sampling points between 2 different LFPs.")
#+end_src
#+RESULTS: minimal-interval
: MINIMAL-INTERVAL
#+name: lfp-pos
#+begin_src lisp
(defparameter lfp-pos
(let ((res (list (car events-pos)))
(n (length events-pos))
(last (car events-pos)))
(do* ((i 1 (1+ i))
(new (nth i events-pos) (nth i events-pos)))
((>= i n) (reverse res))
(if (>= (- new last) minimal-interval)
(progn
(push new res)
(setf last new))))))
#+end_src
#+RESULTS: lfp-pos
: LFP-POS
We can convert the result in seconds in order to compare it with what we obtained with =R=:
#+begin_src lisp
(defparameter lfp-pos2 (copy-list lfp-pos))
#+end_src
#+RESULTS:
: LFP-POS2
#+begin_src lisp
(map-into lfp-pos2 #'(lambda (x) (float (/ (1+ x) 500))) lfp-pos)
#+end_src
#+RESULTS:
| 28.164 | 73.984 | 126.06 | 192.304 | 257.194 | 293.352 | 342.952 | 465.13 | 470.484 | 477.416 |
* Figure illustrating the result
We are going to use =Gnuplot= to build our summary figure. We need first to write to disk the normalized derivatives we want to plot.
We are going to write the derivatives as floats and for that we need the =IEEE-Floats= package which is easily installed (if necessary) and loaded with [[http://www.quicklisp.org/][quicklisp]] as follows:
#+begin_src lisp
(ql:quickload "ieee-floats")
#+end_src
#+RESULTS:
| ieee-floats |
We define a function writing to disk the normalized derivative of a channel using a 32 bits encoding:
#+name: write-normalized-derivatives
#+begin_src lisp :cache yes
(defun write-normalized-derivatives (derivative sd file-name)
(with-open-file (str file-name :direction :output
:element-type '(unsigned-byte 32)
:if-exists :overwrite
:if-does-not-exist :create)
(let ((n (length derivative)))
(do* ((i 0 (1+ i))
(x (ieee-floats:encode-float32
(/ (aref derivative i) sd))
(ieee-floats:encode-float32
(/ (aref derivative i) sd))))
((>= i (1- n)) 'done)
(write-byte x str)))))
#+end_src
#+RESULTS[35ef01356040fab33877e1170f8e36a33b0b1a59]: write-normalized-derivatives
: WRITE-NORMALIZED-DERIVATIVES
We now write to disk our four normalized derivative traces:
#+name: processed-data-1
#+headers: :exports none
#+begin_src lisp :cache yes
(write-normalized-derivatives (nth 0 datad-list)
(nth 0 sd-derivative-values)
"SpinalCordMouseEmbryo_dn_CH1.dat")
#+end_src
#+RESULTS[1051febd2a04b88961843e96a667fd87ed9d0bd7]: processed-data-1
: DONE
#+name: processed-data-2
#+headers: :exports none
#+begin_src lisp :cache yes
(write-normalized-derivatives (nth 1 datad-list)
(nth 1 sd-derivative-values)
"SpinalCordMouseEmbryo_dn_CH2.dat")
#+end_src
#+RESULTS[6e66b8ca5d505d84e9aef2fd881312b875dfc913]: processed-data-2
: DONE
#+name: processed-data-3
#+headers: :exports none
#+begin_src lisp :cache yes
(write-normalized-derivatives (nth 2 datad-list)
(nth 2 sd-derivative-values)
"SpinalCordMouseEmbryo_dn_CH3.dat")
#+end_src
#+RESULTS[2a016774a11052f14323f6c7a307ee0d6246e0f1]: processed-data-3
: DONE
#+name: processed-data-4
#+headers: :exports none
#+begin_src lisp :cache yes
(write-normalized-derivatives (nth 3 datad-list)
(nth 3 sd-derivative-values)
"SpinalCordMouseEmbryo_dn_CH4.dat")
#+end_src
#+RESULTS[242fe86f2b9bd1f3d773db81434501328204bf3a]: processed-data-4
: DONE
We extract the global minimal and maximal values of the normalized derivatives in order to impose a common scale on all of our plot:
#+name: global-min
#+begin_src lisp :cache yes
(1- (floor
(apply #'min (mapcar #'(lambda (l s) (/ (reduce #'min l) s))
datad-list
sd-derivative-values))))
#+end_src
#+RESULTS[630217b07c83d03c5d8f8dc5849ee793604c804b]: global-min
: -14
#+name: global-max
#+begin_src lisp :cache yes
(ceiling (apply #'max (mapcar #'(lambda (l s) (/ (reduce #'max l) s))
datad-list
sd-derivative-values)))
#+end_src
#+RESULTS[e79f039622f45dd4ef2e3be650d62b84bb8ba2b9]: global-max
: 7
#+name: figure-lfp-detection
#+header: :var minimum=global-min :var maximum=global-max
#+begin_src gnuplot :exports code :file lfdDetection.png :exports results
set border 0
unset xtics
unset ytics
set tmargin 0
set bmargin 0
set lmargin 2
set rmargin 2
set yrange [minimum:maximum]
unset key
set multiplot layout 4,1
set arrow from 0,-4 to 300000,-4 nohead lt 3 lc rgb "red"
set label 1 "Channel 1" at 250000, 5
plot 'SpinalCordMouseEmbryo_dn_CH1.dat' binary format='%float32' using 1 with lines lc rgb "black"
set label 1 "Channel 2" at 250000, 5
plot 'SpinalCordMouseEmbryo_dn_CH2.dat' binary format='%float32' using 1 with lines lc rgb "black"
set label 1 "Channel 3" at 250000, 5
plot 'SpinalCordMouseEmbryo_dn_CH3.dat' binary format='%float32' using 1 with lines lc rgb "black"
set label 1 "Channel 4" at 250000, 5
plot 'SpinalCordMouseEmbryo_dn_CH4.dat' binary format='%float32' using 1 with lines lc rgb "black"
unset multiplot
#+end_src
#+RESULTS: figure-lfp-detection
[[file:lfdDetection.png]]