gin commit from pteromys

New files: 11
Modified files: 1

README.md

@@ -1,3 +1,120 @@
 # p_discolor_auditory_feedback
 
-ABR recordings, acoustic call parameters and analysis programs, as reported in Lattenkamp et al., “Learning how to call: the vocal development of the pale spear-nosed bat is dependent on auditory feedback”.
+Acoustic parameters of P. discolor vocalisations, auditory brainstem
+responses, and data analysis programs for Lattenkamp et al., “Learning
+how to call: the vocal development of the pale spear-nosed bat is
+dependent on auditory feedback”.
+
+## Vocalisations
+
+The `vocalisations` directory provides data in four CSV files:
+
+* `pup_calls.csv`, `adult_calls.csv`: Acoustic call parameters of the
+  vocalisations recorded from six bats, three of which were deafened
+  after the first recording session.  In `pup_calls.csv`, these bats
+  were up to half a year old; in `adult_calls.csv`, the same individuals
+  were 2 years old.  The juvenile recordings were acquired with pairs of
+  bats: (deafened or normal-hearing) juveniles and their (always
+  normal-hearing) respective mothers, as well as with pairs of one
+  deafened and one normal-hearing juvenile.  The adult recordings were
+  acquired from groups of two or three bats such that within each
+  recording session, all were either deafened or of normal hearing.
+
+  `pup_calls.csv` includes the vocalisations from the juvenile bats’
+  respective mothers which were placed in the recording chamber with
+  their offspring in some of the sessions.  These data have not been
+  analysed further by the authors.
+
+* `pup_sessions.csv`, `adult_sessions.csv`: Metadata for the recording
+  sessions.
+
+### Acoustic call parameters
+
+Each row of `*_calls.csv` corresponds to one extracted vocalisation from
+one bat.  The columns are:
+
+* `session_id`: ID of the recording session in which this vocalisation
+  occurred, for lookup in the corresponding `*_sessions.csv` file.
+* `call_id`: Sequence number of the vocalisation within the recording
+  session.  The pairs (`session_id`, `call_id`) are unique within each
+  CSV file.
+* `start_sample`: Sample index in the microphone buffer where the
+  extracted vocalisation started.
+* `start_time`: Absolute time in Central European (Summer) Time at which
+  the extracted vocalisation started. Can be calculated from
+  `start_sample` and the metadata in `*_sessions.csv`, based on the
+  sampling rate of 192 kHz.
+* `calling_bat` (`pup_calls.csv`): Identifier of the animal which
+  emitted the vocalisation.  This can be determined by comparing levels
+  in the microphone array in the recording chamber, as the position of
+  the bats within the chamber was known and fixed.  This column is empty
+  in `adult_calls.csv`, where the recording setup did not facilitate
+  this analysis.
+* `other_bat` (`pup_calls.csv`): Identifier of the animal which was in
+  the recording chamber with the calling bat, but did not emit the
+  vocalisation.  Empty in `adult_calls.csv`, where this information is
+  unknown.
+* `calling_bat_deaf`: 1 if the calling bat was deafened, 0 otherwise.
+  Note that in `pup_calls.csv`, the first session of each
+  juvenile–mother pairing was recorded before deafening, yet
+  `calling_bat_deaf` is 1 even in those sessions.
+* `before_deafening`: 1 if this vocalisation occurred in the first
+  session of a juvenile (before any deafening), 0 otherwise.
+* `calling_bat_mother`: 1 if the calling bat was a mother, 0 otherwise.
+* `level_difference` (`pup_calls.csv`; in dB): Difference between the
+  root-mean-square vocalisation levels of the microphones closer to one
+  bat vs. the microphones closer to the other bat.  Used to determine
+  `calling_bat` and `other_bat` in sessions with juveniles.
+* `call_rms` (in dB): Root-mean-square sound pressure of the
+  vocalisation at the microphone were it was loudest, relative to an
+  arbitrary reference which was constant across sessions in
+  `pup_calls.csv` and `adult_calls.csv`, but not the same in both.
+* `call_duration` (in s): Duration of the vocalisation.
+* `mean_aperiodicity`: Aperiodicity value extracted with the YIN
+  algorithm for fundamental frequency estimation.
+* `f0_mean`, `f0_min`, `f0_max`, `f0_start`, `f0_end`, `f0_slope` (in
+  Hz): Mean, minimum, maximum, initial and final fundamental frequencies
+  of the vocalisation, respectively, as estimated with the YIN
+  algorithm.
+* `spectral_centroid` (in Hz): The center of mass of the frequency
+  spectrum of the vocalisation.
+
+### Session files
+
+`pup_sessions.csv` contains the session IDs (`session_id`), start times
+in Central European (Summer) Time (`start_time`), the IDs of the two
+individuals which were placed in the recording chamber together in each
+session (`animal1` and `animal2`), and `before_deafening` indicating
+whether a recording was the first of a juvenile.  `adult_sessions.csv`
+is structured similarly, but replaces the two `animalX` columns with
+`group`, which is either `deaf` or `hearing`, and does not contain
+`before_deafening`.
+
+### Animals
+
+The individuals referred to in the data files were:
+
+* `b2`, born on 2017/01/26, never deafened.
+* `b3`, born on 2017/01/29, deafened after the first recording session.
+* `b4`, born on 2017/01/30, never deafened.
+* `b5`, born on 2017/02/02, deafened after the first recording session.
+* `b7`, born on 2017/02/08, never deafened.
+* `b8`, born on 2017/02/08, deafened after the first recording session.
+* `mh`, `mt`, `mb`, `my` and `mp`, the mothers of these six bats.
+
+### Scripts
+
+`scripts` contains the Jupyter notebook which was used to create the
+figures and tables included in the publication.
+
+`scripts_for_reference` is of documentary nature and contains the MATLAB
+and Python code which was used to extract the acoustic parameters and
+generate the CSV files from the raw microphone data (which is not
+included in this repository for reasons of space).  The YIN
+implementation used by some of the MATLAB scripts is available from
+http://audition.ens.fr/adc/sw/yin.zip (de Cheveigné & Kawahara, JASA
+111:1917:1930, 2002, doi:10.1121/1.1458024).
+
+## Auditory brainstem responses
+
+TODO

datacite.yml

@@ -0,0 +1,53 @@
+authors:
+  -
+    firstname: "Ella Z."
+    lastname: "Lattenkamp"
+    affiliation: "Department Biology II, LMU München, Martinsried, Germany; Neurogenetics of Vocal Communication Group, MPI for Psycholinguistics, Nijmegen"
+    id: "ORCID:0000-0002-8928-8770"
+  -
+    firstname: "Meike"
+    lastname: "Linnenschmidt"
+    affiliation: "Department Biology II, LMU München, Martinsried, Germany"
+    id: "ORCID:0000-0001-7090-8610"
+  -
+    firstname: "Eva"
+    lastname: "Mardus"
+    affiliation: "Department Biology II, LMU München, Martinsried, Germany"
+  -
+    firstname: "Sonja C."
+    lastname: "Vernes"
+    affiliation: "Neurogenetics of Vocal Communication Group, MPI for Psycholinguistics, Nijmegen"
+    id: "ORCID:0000-0003-0305-4584"
+  -
+    firstname: "Lutz"
+    lastname: "Wiegrebe"
+    affiliation: "Department Biology II, LMU München, Martinsried, Germany"
+    id: "ORCID:0000-0002-9289-6187"
+  -
+    firstname: "Michael"
+    lastname: "Schutte"
+    affiliation: "Department Biology II, LMU München, Martinsried, Germany"
+    id: "ORCID:0000-0002-2017-3465"
+
+title: "Data for “Learning how to call: the vocal development of the pale spear-nosed bat is dependent on auditory feedback”"
+
+description: |
+  Acoustic call parameters, auditory brainstem responses, and data analysis
+  programs for Lattenkamp et al., “Learning how to call: the vocal development
+  of the pale spear-nosed bat is dependent on auditory feedback”.
+
+keywords:
+  - Phyllostomus discolor
+  - vocal learning
+  - deafening
+  - vocal development
+  - auditory feedback
+  - hearing impairment
+
+license:
+  name: "CC-BY-NC-SA"
+  url: "https://creativecommons.org/licenses/by-nc-sa/4.0/"
+
+funding:
+  - "HFSP, RGP0058/2016"
+  - "DFG, wi1518/17"

vocalisations/adult_sessions.csv

@@ -0,0 +1,219 @@
+session_id,group,start_time
+0,deaf,2018-11-22 13:35:39
+1,deaf,2018-11-22 13:37:27
+2,deaf,2018-11-22 13:42:36
+3,deaf,2018-11-22 13:43:57
+4,deaf,2018-11-22 13:44:39
+5,deaf,2018-11-22 13:45:12
+6,deaf,2018-11-22 13:46:17
+7,deaf,2018-11-22 13:46:55
+8,deaf,2018-11-22 13:47:36
+9,deaf,2018-11-22 13:49:35
+10,deaf,2018-11-22 13:50:14
+11,deaf,2018-11-22 13:56:27
+12,deaf,2018-11-22 13:59:20
+13,deaf,2018-11-22 14:01:33
+14,deaf,2018-11-22 14:04:11
+15,deaf,2018-11-22 14:05:51
+16,deaf,2018-11-22 14:08:23
+17,deaf,2018-11-22 14:11:30
+18,deaf,2018-11-22 14:12:28
+19,deaf,2018-11-22 14:15:22
+20,deaf,2018-11-22 14:16:18
+21,deaf,2018-11-22 14:18:21
+22,deaf,2018-11-22 14:21:24
+23,deaf,2018-11-22 14:22:41
+24,deaf,2018-11-22 14:25:27
+25,deaf,2018-11-22 14:27:42
+26,deaf,2018-11-22 14:29:27
+27,deaf,2018-11-22 14:30:33
+28,deaf,2018-11-22 14:31:01
+29,deaf,2018-11-22 14:32:35
+30,deaf,2018-11-22 14:33:24
+31,deaf,2018-11-22 14:34:22
+32,deaf,2018-11-22 14:37:36
+33,deaf,2018-11-22 14:42:50
+34,deaf,2018-11-22 14:43:57
+35,deaf,2018-11-22 14:45:11
+36,deaf,2018-11-22 14:47:07
+37,deaf,2018-11-22 14:47:36
+38,deaf,2018-11-22 14:48:59
+39,deaf,2018-11-22 14:51:24
+40,deaf,2018-11-22 14:55:23
+41,deaf,2018-11-22 14:59:12
+42,deaf,2018-11-22 15:00:05
+43,deaf,2018-11-22 15:00:41
+44,deaf,2018-11-22 15:03:25
+45,deaf,2018-11-22 15:06:43
+46,deaf,2018-11-22 15:07:25
+47,deaf,2018-11-22 15:09:17
+48,deaf,2018-11-22 15:17:50
+49,deaf,2018-11-22 15:18:34
+50,deaf,2018-11-22 15:20:32
+51,deaf,2018-11-27 13:43:38
+52,deaf,2018-11-27 13:50:25
+53,deaf,2018-11-27 13:51:46
+54,deaf,2018-11-27 13:54:22
+55,deaf,2018-11-27 13:56:06
+56,deaf,2018-11-27 14:00:30
+57,deaf,2018-11-27 14:03:24
+58,deaf,2018-11-27 14:04:22
+59,deaf,2018-11-27 14:05:02
+60,deaf,2018-11-27 14:05:32
+61,deaf,2018-11-27 14:06:41
+62,deaf,2018-11-27 14:08:49
+63,deaf,2018-11-27 14:12:19
+64,deaf,2018-11-27 14:13:04
+65,deaf,2018-11-27 14:16:01
+66,deaf,2018-11-27 14:16:55
+67,deaf,2018-11-27 14:18:48
+68,deaf,2018-11-27 14:19:33
+69,deaf,2018-11-27 14:20:43
+70,deaf,2018-11-27 14:22:47
+71,deaf,2018-11-27 14:23:16
+72,deaf,2018-11-27 14:23:43
+73,deaf,2018-11-27 14:30:02
+74,deaf,2018-11-27 14:32:18
+75,deaf,2018-11-27 14:35:09
+76,deaf,2018-11-27 14:35:49
+77,deaf,2018-11-27 14:36:23
+78,deaf,2018-11-27 14:37:14
+79,deaf,2018-11-27 14:39:55
+80,deaf,2018-11-27 14:44:24
+81,deaf,2018-11-27 14:45:16
+82,deaf,2018-11-27 14:51:59
+83,deaf,2018-11-27 14:52:57
+84,deaf,2018-11-27 14:54:24
+85,deaf,2018-11-27 14:56:27
+86,deaf,2018-11-27 14:57:35
+87,deaf,2018-11-27 14:58:57
+88,deaf,2018-11-27 15:13:55
+89,deaf,2018-11-27 15:14:54
+90,deaf,2018-11-27 15:15:24
+91,deaf,2018-11-27 15:16:30
+92,deaf,2018-11-27 15:18:21
+93,deaf,2018-11-27 15:19:15
+94,deaf,2018-11-27 15:20:26
+95,deaf,2018-11-27 15:21:20
+96,deaf,2018-11-27 15:24:15
+97,deaf,2018-11-27 15:26:32
+98,deaf,2018-11-27 15:27:35
+99,deaf,2018-11-27 15:29:50
+100,deaf,2018-11-27 15:32:36
+101,deaf,2018-11-27 15:34:01
+102,deaf,2018-11-27 15:34:56
+103,deaf,2018-11-28 13:44:31
+104,deaf,2018-11-28 13:48:21
+105,deaf,2018-11-28 13:59:17
+106,deaf,2018-11-28 14:01:08
+107,deaf,2018-11-28 14:10:06
+108,deaf,2018-11-28 14:11:39
+109,deaf,2018-11-28 14:15:19
+110,deaf,2018-11-28 14:28:00
+111,deaf,2018-11-28 14:32:01
+112,deaf,2018-11-28 14:37:03
+113,deaf,2018-11-28 14:56:11
+114,deaf,2018-11-28 14:57:15
+115,deaf,2018-11-28 14:59:03
+116,deaf,2018-11-28 15:00:25
+117,deaf,2018-11-28 15:09:11
+118,deaf,2018-11-29 12:51:33
+119,deaf,2018-11-29 12:52:02
+120,deaf,2018-11-29 12:55:16
+121,deaf,2018-11-29 12:56:38
+122,deaf,2018-11-29 12:57:43
+123,deaf,2018-11-29 13:00:08
+124,deaf,2018-11-29 13:01:30
+125,deaf,2018-11-29 13:03:25
+126,deaf,2018-11-29 13:04:33
+127,deaf,2018-11-29 13:14:33
+128,deaf,2018-11-29 13:22:07
+129,deaf,2018-11-29 13:22:55
+130,deaf,2018-11-29 13:23:31
+131,deaf,2018-11-29 13:24:00
+132,deaf,2018-11-29 13:24:33
+133,deaf,2018-11-29 13:25:22
+134,deaf,2018-11-29 13:27:36
+135,deaf,2018-11-29 13:28:44
+136,deaf,2018-11-29 13:35:29
+137,deaf,2018-11-29 13:36:27
+138,deaf,2018-11-29 13:37:58
+139,deaf,2018-11-29 13:40:02
+140,deaf,2018-11-29 13:41:29
+141,deaf,2018-11-29 13:50:35
+142,deaf,2018-11-29 14:10:29
+143,deaf,2018-11-29 14:13:06
+144,deaf,2018-11-29 14:23:44
+145,deaf,2018-11-29 14:28:26
+146,deaf,2018-11-29 14:43:11
+147,deaf,2018-11-30 10:27:25
+148,deaf,2018-11-30 10:28:34
+149,deaf,2018-11-30 10:31:01
+150,deaf,2018-11-30 10:34:33
+151,deaf,2018-11-30 10:36:44
+152,deaf,2018-11-30 10:37:18
+153,deaf,2018-11-30 10:38:57
+154,deaf,2018-11-30 10:40:39
+155,deaf,2018-11-30 10:41:35
+156,deaf,2018-11-30 10:42:55
+157,deaf,2018-11-30 10:44:19
+158,deaf,2018-11-30 10:45:53
+159,deaf,2018-11-30 10:47:10
+160,deaf,2018-11-30 10:49:24
+161,deaf,2018-11-30 10:52:42
+162,deaf,2018-11-30 10:54:21
+163,deaf,2018-11-30 10:54:56
+164,deaf,2018-11-30 10:56:27
+165,deaf,2018-11-30 10:57:28
+166,deaf,2018-11-30 10:58:05
+167,deaf,2018-11-30 10:59:02
+168,deaf,2018-11-30 11:04:03
+169,deaf,2018-11-30 11:04:29
+170,deaf,2018-11-30 11:09:16
+171,deaf,2018-11-30 11:12:49
+172,deaf,2018-11-30 11:16:20
+173,deaf,2018-11-30 11:19:00
+174,deaf,2018-11-30 11:19:36
+175,deaf,2018-11-30 11:23:05
+176,deaf,2018-11-30 11:26:23
+177,deaf,2018-11-30 11:27:42
+178,deaf,2018-11-30 11:28:38
+179,deaf,2018-11-30 11:29:47
+180,deaf,2018-11-30 11:30:33
+181,deaf,2018-11-30 11:34:39
+182,deaf,2018-11-30 11:35:27
+183,deaf,2018-11-30 11:37:12
+184,deaf,2018-11-30 11:40:42
+185,deaf,2018-11-30 11:42:54
+186,deaf,2018-11-30 11:47:34
+187,deaf,2018-11-30 11:50:52
+188,deaf,2018-11-30 11:52:13
+189,deaf,2018-11-30 11:52:57
+190,deaf,2018-11-30 11:55:00
+191,deaf,2018-11-30 11:58:52
+192,deaf,2018-11-30 12:00:51
+193,deaf,2018-11-30 12:03:10
+194,deaf,2018-11-30 12:04:07
+195,deaf,2018-11-30 12:05:21
+196,deaf,2018-11-30 12:07:11
+197,deaf,2018-11-30 12:10:56
+198,deaf,2018-11-30 12:12:08
+199,deaf,2018-11-30 12:16:38
+200,deaf,2018-11-30 12:21:27
+201,deaf,2018-11-30 12:27:09
+202,hearing,2019-01-09 10:18:25
+203,hearing,2019-01-11 09:33:58
+204,hearing,2019-01-11 09:37:00
+205,hearing,2019-01-21 09:36:24
+206,hearing,2019-01-21 10:00:12
+207,hearing,2019-01-21 10:58:52
+208,hearing,2019-01-21 11:05:51
+209,hearing,2019-01-23 09:42:37
+210,hearing,2019-01-23 10:26:29
+211,hearing,2019-01-24 09:12:27
+212,hearing,2019-01-24 09:14:19
+213,hearing,2019-01-24 09:22:37
+214,hearing,2019-01-24 09:25:26
+215,hearing,2019-01-24 09:26:20
+216,hearing,2019-01-24 09:57:26
+217,hearing,2019-01-24 10:44:06

vocalisations/pup_calls.csv

@@ -0,0 +1 @@
+/annex/objects/MD5-s308030911--99d2bed8ac778abed4536d01ec99b63c

vocalisations/pup_sessions.csv

@@ -0,0 +1,353 @@
+session_id,animal1,animal2,start_time,before_deafening
+0,b2,b3,2017-09-04 10:36:36,0
+1,b2,b3,2017-08-21 10:07:58,0
+2,b2,b3,2017-08-28 10:24:38,0
+3,b2,b5,2017-09-06 10:06:57,0
+4,b2,b5,2017-08-24 11:20:51,0
+5,b2,b5,2017-08-30 09:48:17,0
+6,b2,b8,2017-09-05 10:03:12,0
+7,b2,b8,2017-08-23 10:15:24,0
+8,b2,b8,2017-08-29 09:58:12,0
+9,b2,mh,2017-08-01 14:21:14,0
+10,b2,mh,2017-03-01 10:41:11,0
+11,b2,mh,2017-03-02 13:08:42,0
+12,b2,mh,2017-05-02 14:26:42,0
+13,b2,mh,2017-07-03 10:02:11,0
+14,b2,mh,2017-05-03 13:24:20,0
+15,b2,mh,2017-04-04 11:19:40,0
+16,b2,mh,2017-05-04 13:37:07,0
+17,b2,mh,2017-04-05 12:14:24,0
+18,b2,mh,2017-04-06 11:52:53,0
+19,b2,mh,2017-02-06 14:20:49,1
+20,b2,mh,2017-02-07 14:16:02,0
+21,b2,mh,2017-06-07 10:09:40,0
+22,b2,mh,2017-03-07 14:44:40,0
+23,b2,mh,2017-08-08 17:00:29,0
+24,b2,mh,2017-02-08 13:57:45,0
+25,b2,mh,2017-03-08 13:28:58,0
+26,b2,mh,2017-05-08 10:09:27,0
+27,b2,mh,2017-03-09 14:28:14,0
+28,b2,mh,2017-04-10 10:32:00,0
+29,b2,mh,2017-07-10 10:57:52,0
+30,b2,mh,2017-04-11 12:38:19,0
+31,b2,mh,2017-04-12 10:26:44,0
+32,b2,mh,2017-06-12 10:13:12,0
+33,b2,mh,2017-03-13 16:29:53,0
+34,b2,mh,2017-02-14 15:10:10,0
+35,b2,mh,2017-03-14 13:48:49,0
+36,b2,mh,2017-02-15 09:43:24,0
+37,b2,mh,2017-03-15 11:29:15,0
+38,b2,mh,2017-02-17 14:43:01,0
+39,b2,mh,2017-07-17 09:35:03,0
+40,b2,mh,2017-05-17 09:10:50,0
+41,b2,mh,2017-04-18 13:55:12,0
+42,b2,mh,2017-04-19 11:50:13,0
+43,b2,mh,2017-02-20 13:20:28,0
+44,b2,mh,2017-06-20 14:10:31,0
+45,b2,mh,2017-04-21 09:42:40,0
+46,b2,mh,2017-02-21 09:40:37,0
+47,b2,mh,2017-03-21 11:50:53,0
+48,b2,mh,2017-03-22 12:27:34,0
+49,b2,mh,2017-05-22 10:31:50,0
+50,b2,mh,2017-04-24 13:22:13,0
+51,b2,mh,2017-02-24 13:26:59,0
+52,b2,mh,2017-07-24 09:40:41,0
+53,b2,mh,2017-03-24 10:43:28,0
+54,b2,mh,2017-04-25 09:49:14,0
+55,b2,mh,2017-04-26 10:26:26,0
+56,b2,mh,2017-06-27 14:55:01,0
+57,b2,mh,2017-03-27 13:39:31,0
+58,b2,mh,2017-02-28 13:51:55,0
+59,b2,mh,2017-03-28 09:46:56,0
+60,b2,mh,2017-03-29 10:14:09,0
+61,b2,mh,2017-05-29 10:01:31,0
+62,b3,mt,2017-08-01 14:42:58,0
+63,b3,mt,2017-03-01 11:02:36,0
+64,b3,mt,2017-04-03 14:33:30,0
+65,b3,mt,2017-07-03 10:28:13,0
+66,b3,mt,2017-03-03 13:47:14,0
+67,b3,mt,2017-05-03 13:46:39,0
+68,b3,mt,2017-04-04 12:42:14,0
+69,b3,mt,2017-05-05 09:24:39,0
+70,b3,mt,2017-04-06 12:14:02,0
+71,b3,mt,2017-02-07 07:46:44,1
+72,b3,mt,2017-06-07 10:32:22,0
+73,b3,mt,2017-03-07 15:05:36,0
+74,b3,mt,2017-03-08 13:02:04,0
+75,b3,mt,2017-05-08 10:32:26,0
+76,b3,mt,2017-05-09 13:06:51,0
+77,b3,mt,2017-04-10 10:54:47,0
+78,b3,mt,2017-08-10 09:18:50,0
+79,b3,mt,2017-07-10 12:42:39,0
+80,b3,mt,2017-03-10 11:58:02,0
+81,b3,mt,2017-05-10 10:13:34,0
+82,b3,mt,2017-04-12 10:52:28,0
+83,b3,mt,2017-06-12 10:35:40,0
+84,b3,mt,2017-04-13 14:47:05,0
+85,b3,mt,2017-03-13 16:50:40,0
+86,b3,mt,2017-02-14 15:32:30,0
+87,b3,mt,2017-03-14 13:26:25,0
+88,b3,mt,2017-02-15 10:10:12,0
+89,b3,mt,2017-03-15 12:48:18,0
+90,b3,mt,2017-02-17 13:25:08,0
+91,b3,mt,2017-07-17 10:18:15,0
+92,b3,mt,2017-05-17 09:34:08,0
+93,b3,mt,2017-04-18 14:16:36,0
+94,b3,mt,2017-04-19 13:01:48,0
+95,b3,mt,2017-05-19 09:41:45,0
+96,b3,mt,2017-02-20 13:47:14,0
+97,b3,mt,2017-06-20 13:49:03,0
+98,b3,mt,2017-04-21 10:05:30,0
+99,b3,mt,2017-02-21 10:06:18,0
+100,b3,mt,2017-03-21 12:14:29,0
+101,b3,mt,2017-03-22 12:50:30,0
+102,b3,mt,2017-05-22 10:54:27,0
+103,b3,mt,2017-04-24 13:47:16,0
+104,b3,mt,2017-02-24 12:36:33,0
+105,b3,mt,2017-07-24 10:02:15,0
+106,b3,mt,2017-03-24 11:04:38,0
+107,b3,mt,2017-04-25 10:10:56,0
+108,b3,mt,2017-04-26 10:49:47,0
+109,b3,mt,2017-06-27 14:33:47,0
+110,b3,mt,2017-03-27 14:00:18,0
+111,b3,mt,2017-02-28 13:25:52,0
+112,b3,mt,2017-03-28 10:08:59,0
+113,b3,mt,2017-03-29 10:35:52,0
+114,b3,mt,2017-05-29 10:30:44,0
+115,b4,b3,2017-09-05 10:27:13,0
+116,b4,b3,2017-08-23 10:36:33,0
+117,b4,b3,2017-08-29 10:25:52,0
+118,b4,b5,2017-09-04 11:08:05,0
+119,b4,b5,2017-08-21 10:29:11,0
+120,b4,b5,2017-08-28 10:51:35,0
+121,b4,b8,2017-09-06 10:31:07,0
+122,b4,b8,2017-08-24 12:36:04,0
+123,b4,b8,2017-08-30 10:12:18,0
+124,b4,mb,2017-08-01 15:06:08,0
+125,b4,mb,2017-03-01 11:23:21,0
+126,b4,mb,2017-03-02 13:32:08,0
+127,b4,mb,2017-07-03 10:50:17,0
+128,b4,mb,2017-05-03 14:07:20,0
+129,b4,mb,2017-04-04 13:04:14,0
+130,b4,mb,2017-05-04 14:02:07,0
+131,b4,mb,2017-04-05 12:35:53,0
+132,b4,mb,2017-05-05 10:57:30,0
+133,b4,mb,2017-04-06 12:36:07,0
+134,b4,mb,2017-03-07 15:29:39,0
+135,b4,mb,2017-02-08 08:14:36,1
+136,b4,mb,2017-06-08 10:22:51,0
+137,b4,mb,2017-03-08 11:20:17,0
+138,b4,mb,2017-05-08 10:54:24,0
+139,b4,mb,2017-05-09 13:30:09,0
+140,b4,mb,2017-04-10 11:20:43,0
+141,b4,mb,2017-08-10 15:34:50,0
+142,b4,mb,2017-07-10 10:36:40,0
+143,b4,mb,2017-03-10 12:20:13,0
+144,b4,mb,2017-05-10 10:38:24,0
+145,b4,mb,2017-04-11 13:07:17,0
+146,b4,mb,2017-04-12 11:28:33,0
+147,b4,mb,2017-06-13 09:56:21,0
+148,b4,mb,2017-03-13 17:11:36,0
+149,b4,mb,2017-02-14 15:58:02,0
+150,b4,mb,2017-03-14 10:53:41,0
+151,b4,mb,2017-02-15 10:35:26,0
+152,b4,mb,2017-03-15 13:11:45,0
+153,b4,mb,2017-02-17 14:16:54,0
+154,b4,mb,2017-07-17 09:56:38,0
+155,b4,mb,2017-05-17 09:57:19,0
+156,b4,mb,2017-04-18 14:40:42,0
+157,b4,mb,2017-05-18 10:36:52,0
+158,b4,mb,2017-04-19 13:22:50,0
+159,b4,mb,2017-02-20 14:13:10,0
+160,b4,mb,2017-04-21 10:28:10,0
+161,b4,mb,2017-02-21 10:30:02,0
+162,b4,mb,2017-06-21 10:20:13,0
+163,b4,mb,2017-03-21 12:36:21,0
+164,b4,mb,2017-03-22 13:12:12,0
+165,b4,mb,2017-05-22 11:16:36,0
+166,b4,mb,2017-04-24 14:11:10,0
+167,b4,mb,2017-07-24 15:36:44,0
+168,b4,mb,2017-03-24 11:28:24,0
+169,b4,mb,2017-04-25 10:32:29,0
+170,b4,mb,2017-04-26 11:12:29,0
+171,b4,mb,2017-06-27 14:12:24,0
+172,b4,mb,2017-03-27 14:22:00,0
+173,b4,mb,2017-02-28 14:15:09,0
+174,b4,mb,2017-03-28 10:31:15,0
+175,b4,mb,2017-03-29 10:56:26,0
+176,b4,mb,2017-05-29 10:58:47,0
+177,b5,mt,2017-08-01 11:11:55,0
+178,b5,mt,2017-03-01 12:54:12,0
+179,b5,mt,2017-04-03 14:54:18,0
+180,b5,mt,2017-07-03 11:11:12,0
+181,b5,mt,2017-03-03 13:18:51,0
+182,b5,mt,2017-05-03 14:28:17,0
+183,b5,mt,2017-04-04 13:26:07,0
+184,b5,mt,2017-05-04 14:26:06,0
+185,b5,mt,2017-04-05 12:58:32,0
+186,b5,mt,2017-05-05 09:47:41,0
+187,b5,mt,2017-03-07 15:50:41,0
+188,b5,mt,2017-08-08 11:18:14,0
+189,b5,mt,2017-06-08 09:12:53,0
+190,b5,mt,2017-03-08 10:58:07,0
+191,b5,mt,2017-05-08 11:17:43,0
+192,b5,mt,2017-05-09 14:09:49,0
+193,b5,mt,2017-04-10 14:35:35,0
+194,b5,mt,2017-07-10 11:19:28,0
+195,b5,mt,2017-03-10 12:42:32,0
+196,b5,mt,2017-05-10 11:01:03,0
+197,b5,mt,2017-04-11 13:35:11,0
+198,b5,mt,2017-04-12 12:53:32,0
+199,b5,mt,2017-06-12 09:50:33,0
+200,b5,mt,2017-02-13 09:32:54,1
+201,b5,mt,2017-03-13 17:32:25,0
+202,b5,mt,2017-08-14 09:54:44,0
+203,b5,mt,2017-03-14 10:32:31,0
+204,b5,mt,2017-03-15 13:32:15,0
+205,b5,mt,2017-07-17 10:39:39,0
+206,b5,mt,2017-05-17 10:20:46,0
+207,b5,mt,2017-04-18 15:01:30,0
+208,b5,mt,2017-05-18 10:59:48,0
+209,b5,mt,2017-04-19 12:40:30,0
+210,b5,mt,2017-05-19 10:04:35,0
+211,b5,mt,2017-02-20 14:41:28,0
+212,b5,mt,2017-04-21 10:51:38,0
+213,b5,mt,2017-02-21 10:58:02,0
+214,b5,mt,2017-06-21 10:41:49,0
+215,b5,mt,2017-03-21 13:02:13,0
+216,b5,mt,2017-03-22 13:34:02,0
+217,b5,mt,2017-05-22 11:39:26,0
+218,b5,mt,2017-05-23 09:04:59,0
+219,b5,mt,2017-04-24 14:35:19,0
+220,b5,mt,2017-02-24 12:06:41,0
+221,b5,mt,2017-07-24 13:26:29,0
+222,b5,mt,2017-03-24 12:03:11,0
+223,b5,mt,2017-05-24 09:04:00,0
+224,b5,mt,2017-04-25 10:54:00,0
+225,b5,mt,2017-04-26 11:37:10,0
+226,b5,mt,2017-06-26 10:20:06,0
+227,b5,mt,2017-03-27 14:44:38,0
+228,b5,mt,2017-02-28 14:36:31,0
+229,b5,mt,2017-03-28 10:58:48,0
+230,b5,mt,2017-03-29 11:17:57,0
+231,b5,mt,2017-05-29 12:38:32,0
+232,b7,b3,2017-09-06 10:55:05,0
+233,b7,b3,2017-09-06 11:07:41,0
+234,b7,b3,2017-08-24 12:57:24,0
+235,b7,b3,2017-08-30 10:35:41,0
+236,b7,b5,2017-09-05 10:50:27,0
+237,b7,b5,2017-08-23 10:58:00,0
+238,b7,b5,2017-08-29 10:49:54,0
+239,b7,b8,2017-09-04 13:11:03,0
+240,b7,b8,2017-08-21 10:50:58,0
+241,b7,b8,2017-08-28 11:14:30,0
+242,b7,my,2017-08-01 15:59:41,0
+243,b7,my,2017-06-01 09:06:19,0
+244,b7,my,2017-03-01 13:15:36,0
+245,b7,my,2017-03-02 14:16:22,0
+246,b7,my,2017-04-03 15:15:52,0
+247,b7,my,2017-07-03 12:44:16,0
+248,b7,my,2017-05-03 14:49:58,0
+249,b7,my,2017-04-04 13:47:41,0
+250,b7,my,2017-05-04 14:48:43,0
+251,b7,my,2017-04-05 13:19:36,0
+252,b7,my,2017-05-05 10:10:46,0
+253,b7,my,2017-03-07 16:11:32,0
+254,b7,my,2017-06-08 09:35:01,0
+255,b7,my,2017-03-08 10:36:15,0
+256,b7,my,2017-05-08 12:45:30,0
+257,b7,my,2017-08-09 09:56:35,0
+258,b7,my,2017-03-09 15:16:07,0
+259,b7,my,2017-07-10 13:03:18,0
+260,b7,my,2017-05-10 11:24:42,0
+261,b7,my,2017-04-11 16:12:14,0
+262,b7,my,2017-05-11 10:23:10,0
+263,b7,my,2017-04-12 13:15:45,0
+264,b7,my,2017-06-12 15:17:55,0
+265,b7,my,2017-08-14 12:33:00,0
+266,b7,my,2017-03-14 10:11:07,0
+267,b7,my,2017-03-15 13:53:38,0
+268,b7,my,2017-03-16 10:36:41,0
+269,b7,my,2017-02-17 13:49:39,1
+270,b7,my,2017-07-17 11:00:50,0
+271,b7,my,2017-05-17 10:43:41,0
+272,b7,my,2017-04-18 15:22:33,0
+273,b7,my,2017-05-18 11:22:41,0
+274,b7,my,2017-04-19 13:43:46,0
+275,b7,my,2017-05-19 10:28:18,0
+276,b7,my,2017-02-20 15:08:20,0
+277,b7,my,2017-04-21 11:13:26,0
+278,b7,my,2017-02-21 11:24:00,0
+279,b7,my,2017-06-21 11:02:33,0
+280,b7,my,2017-03-21 14:02:26,0
+281,b7,my,2017-03-22 13:56:12,0
+282,b7,my,2017-05-22 12:01:51,0
+283,b7,my,2017-05-23 09:27:50,0
+284,b7,my,2017-07-24 16:03:10,0
+285,b7,my,2017-03-24 12:25:15,0
+286,b7,my,2017-05-24 09:26:20,0
+287,b7,my,2017-04-25 11:18:10,0
+288,b7,my,2017-04-26 13:07:22,0
+289,b7,my,2017-06-26 10:40:37,0
+290,b7,my,2017-04-27 11:09:33,0
+291,b7,my,2017-02-28 10:55:24,0
+292,b7,my,2017-03-28 11:20:23,0
+293,b7,my,2017-03-29 12:58:26,0
+294,b7,my,2017-05-29 13:08:24,0
+295,b7,my,2017-03-30 15:41:15,0
+296,b7,my,2017-05-30 09:15:21,0
+297,b8,mp,2017-08-01 16:21:54,0
+298,b8,mp,2017-06-01 09:34:19,0
+299,b8,mp,2017-03-01 13:36:38,0
+300,b8,mp,2017-03-02 13:54:18,0
+301,b8,mp,2017-07-03 13:05:29,0
+302,b8,mp,2017-05-03 15:14:06,0
+303,b8,mp,2017-04-04 14:08:45,0
+304,b8,mp,2017-04-05 13:41:41,0
+305,b8,mp,2017-05-05 10:34:16,0
+306,b8,mp,2017-04-06 13:50:58,0
+307,b8,mp,2017-03-07 16:32:28,0
+308,b8,mp,2017-06-08 09:57:29,0
+309,b8,mp,2017-03-08 10:15:49,0
+310,b8,mp,2017-05-08 13:14:55,0
+311,b8,mp,2017-08-09 12:59:38,0
+312,b8,mp,2017-03-09 14:54:46,0
+313,b8,mp,2017-04-10 15:00:07,0
+314,b8,mp,2017-07-10 13:24:24,0
+315,b8,mp,2017-05-10 11:47:37,0
+316,b8,mp,2017-05-11 10:46:03,0
+317,b8,mp,2017-04-12 13:39:42,0
+318,b8,mp,2017-06-12 15:39:46,0
+319,b8,mp,2017-04-13 15:29:25,0
+320,b8,mp,2017-08-14 12:55:14,0
+321,b8,mp,2017-03-14 09:49:14,0
+322,b8,mp,2017-03-15 14:15:38,0
+323,b8,mp,2017-03-16 10:12:30,0
+324,b8,mp,2017-02-17 08:40:52,1
+325,b8,mp,2017-07-17 11:21:27,0
+326,b8,mp,2017-05-17 11:06:27,0
+327,b8,mp,2017-04-18 15:43:41,0
+328,b8,mp,2017-05-18 11:44:33,0
+329,b8,mp,2017-04-19 12:15:07,0
+330,b8,mp,2017-05-19 10:52:19,0
+331,b8,mp,2017-04-21 11:36:15,0
+332,b8,mp,2017-02-21 11:48:47,0
+333,b8,mp,2017-06-21 11:23:28,0
+334,b8,mp,2017-03-21 14:27:55,0
+335,b8,mp,2017-03-22 14:17:47,0
+336,b8,mp,2017-05-22 12:23:42,0
+337,b8,mp,2017-05-23 09:50:06,0
+338,b8,mp,2017-02-24 13:52:45,0
+339,b8,mp,2017-07-24 16:25:30,0
+340,b8,mp,2017-03-24 12:50:30,0
+341,b8,mp,2017-05-24 09:49:06,0
+342,b8,mp,2017-04-25 11:41:45,0
+343,b8,mp,2017-04-26 13:30:52,0
+344,b8,mp,2017-06-26 11:01:26,0
+345,b8,mp,2017-04-27 11:30:33,0
+346,b8,mp,2017-03-27 15:07:07,0
+347,b8,mp,2017-02-28 14:57:24,0
+348,b8,mp,2017-03-28 13:00:08,0
+349,b8,mp,2017-03-29 13:19:26,0
+350,b8,mp,2017-05-29 13:33:43,0
+351,b8,mp,2017-05-30 10:55:03,0

vocalisations/scripts/.ipynb_checkpoints/figures-checkpoint.ipynb

@@ -0,0 +1,714 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib inline"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import IPython.display as ipd\n",
+    "import sys\n",
+    "import io\n",
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "import matplotlib as mpl\n",
+    "import matplotlib.pyplot as plt\n",
+    "import scipy.stats as sps"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Global information\n",
+    "==="
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mpl.rcParams[\"font.family\"] = \"TeX Gyre Pagella\"\n",
+    "mpl.rcParams[\"font.size\"] = 8\n",
+    "rst_column_width = 3.3\n",
+    "rst_total_width = 7"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "deaf_pups = {\"b3\", \"b5\", \"b8\"}\n",
+    "hearing_pups = {\"b2\", \"b4\", \"b7\"}\n",
+    "all_pups = deaf_pups.union(hearing_pups)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "group_colors = dict(hearing=\"#7bb992\", deaf=\"#5f0f00\")\n",
+    "bat_colors = {**{bat: group_colors[\"deaf\"] for bat in deaf_pups},\n",
+    "              **{bat: group_colors[\"hearing\"] for bat in hearing_pups}}\n",
+    "bat_markers = dict(zip(sorted(all_pups), [\"P\", \"o\", \"*\", \"s\", \"X\", \"D\"]))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Load data for juveniles\n",
+    "==="
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pup_sessions = pd.read_csv(\"../pup_sessions.csv\", index_col=0, parse_dates=[\"start_time\"],\n",
+    "                           dtype=dict(before_deafening=np.bool))\n",
+    "pup_calls = pd.read_csv(\"../pup_calls.csv\", index_col=[0, 1], parse_dates=[\"start_time\"], na_values=[\"--\"],\n",
+    "                        dtype=dict(calling_bat_deaf=np.bool, calling_bat_mother=np.bool, before_deafening=np.bool))\n",
+    "\n",
+    "# filter out mother vocalisations\n",
+    "pup_calls = pup_calls[~pup_calls[\"calling_bat_mother\"]].reset_index()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "bat_birthdates = pd.DataFrame([(\"b2\", \"2017/01/26\"),\n",
+    "                               (\"b4\", \"2017/01/30\"),\n",
+    "                               (\"b7\", \"2017/02/08\"),\n",
+    "                               (\"b3\", \"2017/01/29\"),\n",
+    "                               (\"b5\", \"2017/02/02\"),\n",
+    "                               (\"b8\", \"2017/02/08\")], columns=[\"bat\", \"birthdate\"])\n",
+    "bat_birthdates[\"birthdate\"] = pd.to_datetime(bat_birthdates[\"birthdate\"])\n",
+    "bat_birthdates.set_index(\"bat\", inplace=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# determine ages (in full days) of the calling bats\n",
+    "\n",
+    "birth_dates = bat_birthdates.loc[pup_calls[\"calling_bat\"]][\"birthdate\"].reset_index(drop=True)\n",
+    "ages = pup_calls[\"start_time\"] - birth_dates\n",
+    "pup_calls[\"calling_bat_age_in_days\"] = ages.dt.days\n",
+    "pup_calls[\"calling_bat_age_in_weeks\"] = ages.dt.days // 7\n",
+    "pup_calls[\"call_rms\"] = pup_calls[\"call_rms\"] + 55"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Load data for adults\n",
+    "==="
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def single_value(values):\n",
+    "    value_set = set(values)\n",
+    "    if len(value_set) != 1:\n",
+    "        raise ValueError(\"non-unity set\")\n",
+    "    return next(iter(value_set))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "adult_calls = pd.read_csv(\"../adult_calls.csv\",\n",
+    "                          parse_dates=[\"start_time\"],\n",
+    "                          na_values=[\"--\"],\n",
+    "                          dtype=dict(calling_bat_deaf=np.bool, session_id=np.int),\n",
+    "                          index_col=[\"session_id\", \"call_id\"])\n",
+    "adult_calls[\"recording_day\"] = adult_calls[\"start_time\"].dt.date\n",
+    "adult_calls[\"call_rms\"] -= adult_calls[\"call_rms\"].min()\n",
+    "adult_calls[\"group\"] = np.where(adult_calls[\"calling_bat_deaf\"], \"deaf\", \"hearing\")\n",
+    "\n",
+    "adult_sessions = pd.read_csv(\"../adult_sessions.csv\",\n",
+    "                             parse_dates=[\"start_time\"],\n",
+    "                             dtype=dict(group=\"category\"),\n",
+    "                             index_col=0)\n",
+    "\n",
+    "adult_recording_days = adult_sessions[[\"group\"]].groupby(adult_sessions[\"start_time\"].dt.date).agg(single_value)\n",
+    "adult_recording_days.index.name = \"recording_day\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Data to be dropped\n",
+    "==="
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Possible echolocation calls (3 ms or shorter)\n",
+    "---"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pup_calls = pup_calls[pup_calls[\"call_duration\"] > 3e-3]\n",
+    "adult_calls = adult_calls[adult_calls[\"call_duration\"] > 3e-3]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Sessions shorter than 20 minutes\n",
+    "---"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "session_gaps = pup_sessions.sort_values(\"start_time\")[\"start_time\"].diff()\n",
+    "new_index = session_gaps.index.insert(0, -1)\n",
+    "new_index = new_index[:-1]\n",
+    "session_gaps.index = new_index\n",
+    "session_gaps.drop(-1, inplace=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ids_to_drop = session_gaps[session_gaps < pd.Timedelta(20, \"m\")].index"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pup_calls = pup_calls[~pup_calls[\"session_id\"].isin(ids_to_drop)]\n",
+    "pup_sessions = pup_sessions.drop(ids_to_drop)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Last two recordings with mother for pups that have more than others\n",
+    "---"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sorted_pup_mother_sessions = pup_sessions[pup_sessions[\"animal2\"].str.startswith(\"m\")].sort_values(\"start_time\")\n",
+    "session_counts_with_mothers = sorted_pup_mother_sessions[\"animal1\"].value_counts()\n",
+    "num_sessions_to_drop = session_counts_with_mothers - session_counts_with_mothers.min()\n",
+    "\n",
+    "ids_to_drop = set()\n",
+    "for bat, drop_count in num_sessions_to_drop.iteritems():\n",
+    "    if drop_count == 0:\n",
+    "        continue\n",
+    "    this_bat_sessions = sorted_pup_mother_sessions[sorted_pup_mother_sessions[\"animal1\"] == bat]\n",
+    "    ids_to_drop = ids_to_drop.union(this_bat_sessions.tail(drop_count).index)\n",
+    "\n",
+    "pup_calls = pup_calls[~pup_calls[\"session_id\"].isin(ids_to_drop)]\n",
+    "pup_sessions = pup_sessions.drop(index=ids_to_drop)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Compute activity statistics (separately for pups before and after deafening)\n",
+    "==="
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "calls_per_pup_and_session = pup_calls.reset_index().groupby([\"calling_bat\", \"session_id\"])[[\"calling_bat\"]].count() / (20*6) # per 10 seconds\n",
+    "calls_per_pup_and_session.columns = [\"calls_per_10s\"]\n",
+    "calls_per_pup_and_session.reset_index(inplace=True)\n",
+    "\n",
+    "session_dates = pup_sessions.loc[calls_per_pup_and_session[\"session_id\"]][\"start_time\"]\n",
+    "birth_dates = bat_birthdates.loc[calls_per_pup_and_session[\"calling_bat\"]][\"birthdate\"]\n",
+    "calls_per_pup_and_session[\"calling_bat_age_in_days\"] = \\\n",
+    "    (session_dates.reset_index(drop=True) - birth_dates.reset_index(drop=True)).dt.days\n",
+    "calls_per_pup_and_session[\"calling_bat_age_in_weeks\"] = \\\n",
+    "    calls_per_pup_and_session[\"calling_bat_age_in_days\"] // 7\n",
+    "calls_per_pup_and_session = pd.merge(calls_per_pup_and_session, pup_sessions[[\"before_deafening\"]],\n",
+    "                                     left_on=\"session_id\", right_index=True,)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dfs = []\n",
+    "\n",
+    "for before_deafening in [False, True]:\n",
+    "    mask = pup_calls[\"before_deafening\"]\n",
+    "    if not before_deafening:\n",
+    "        mask = ~mask\n",
+    "    gr = pup_calls[mask].groupby([\"calling_bat\", \"calling_bat_age_in_weeks\"])\n",
+    "    gr = gr[[\"call_duration\", \"mean_aperiodicity\", \"f0_mean\", \"call_rms\"]]\n",
+    "    params_per_pup_and_week = gr.median()\n",
+    "    params_per_pup_and_week[\"calls_this_week\"] = gr.size()\n",
+    "    \n",
+    "    mask = calls_per_pup_and_session[\"before_deafening\"]\n",
+    "    if not before_deafening:\n",
+    "        mask = ~mask\n",
+    "    gr = calls_per_pup_and_session[mask].groupby([\"calling_bat\", \"calling_bat_age_in_weeks\"])\n",
+    "    params_per_pup_and_week[\"calls_per_10s\"] = gr[\"calls_per_10s\"].mean()\n",
+    "    params_per_pup_and_week.reset_index(inplace=True)\n",
+    "    params_per_pup_and_week[\"before_deafening\"] = before_deafening\n",
+    "    \n",
+    "    dfs.append(params_per_pup_and_week)\n",
+    "    \n",
+    "params_per_pup_and_week = pd.concat(dfs)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "gr = adult_calls.groupby([\"group\", \"session_id\"], observed=True)\n",
+    "calls_per_adult_group_and_session = gr.size().to_frame(name=\"calls_per_10s\").reset_index()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Quantities of extracted data\n",
+    "==="
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(\"Total calls by pups: {}\\nCalls by deaf pups: {} ({:.1f} %)\\nCalls by hearing pups: {} ({:.1f} %)\".format(\n",
+    "        len(pup_calls),\n",
+    "        np.sum(pup_calls[\"calling_bat_deaf\"]),\n",
+    "        np.sum(pup_calls[\"calling_bat_deaf\"]) / len(pup_calls) * 100,\n",
+    "        np.sum(~pup_calls[\"calling_bat_deaf\"]),\n",
+    "        np.sum(~pup_calls[\"calling_bat_deaf\"] / len(pup_calls) * 100)))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(\"Total calls by adults: {}\\nCalls by deaf adults: {} ({:.1f} %)\\nCalls by hearing adults: {} ({:.1f} %)\".format(\n",
+    "        len(adult_calls),\n",
+    "        np.sum(adult_calls[\"calling_bat_deaf\"]),\n",
+    "        np.sum(adult_calls[\"calling_bat_deaf\"]) / len(adult_calls) * 100,\n",
+    "        np.sum(~adult_calls[\"calling_bat_deaf\"]),\n",
+    "        np.sum(~adult_calls[\"calling_bat_deaf\"] / len(adult_calls) * 100)))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Main figure\n",
+    "==="
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def plot_trajectories(ax, values, parameter, scale):\n",
+    "    line_artists = {}\n",
+    "    for bat, df in values.groupby(\"calling_bat\"):\n",
+    "        bat_color, bat_marker = bat_colors[bat], bat_markers[bat]\n",
+    "        line_artist, = ax.plot(df[\"calling_bat_age_in_weeks\"] + 1, df[parameter] * scale,\n",
+    "                               c=bat_color, linewidth=0.8)\n",
+    "        ax.plot(df[\"calling_bat_age_in_weeks\"] + 1, df[parameter] * scale,\n",
+    "                c=bat_color, marker=bat_marker, markersize=5, clip_on=False,\n",
+    "                markerfacecolor=\"none\", linestyle=\"\")\n",
+    "        line_artists[bat] = line_artist\n",
+    "    return line_artists"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def plot_histograms(ax, values, parameter, scale, y_range):\n",
+    "    data = {}\n",
+    "\n",
+    "    if \"group\" in values:\n",
+    "        for group, hearing in [(\"deaf\", False), (\"hearing\", True)]:\n",
+    "            data[hearing] = (values.loc[values[\"group\"] == group, parameter] * scale).values\n",
+    "    else:\n",
+    "        for hearing, df in values.groupby(values[\"calling_bat\"].isin(hearing_pups)):\n",
+    "            data[hearing] = (df[parameter] * scale).values\n",
+    "\n",
+    "    ax.hist([data[True], data[False]], density=True,\n",
+    "            range=y_range, bins=15,\n",
+    "            color=[group_colors[\"hearing\"], group_colors[\"deaf\"]],\n",
+    "            orientation=\"horizontal\", rwidth=0.85)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def plot_boxplots(ax, values, parameter, scale, test_significance=False):\n",
+    "    data = []\n",
+    "    for group in [\"deaf\", \"hearing\"]:\n",
+    "        data.append(scale * values.loc[values[\"group\"] == group, parameter].dropna().values)\n",
+    "\n",
+    "    artists = boxplot_ax.boxplot(data, vert=True, manage_ticks=False, showfliers=False, widths=0.75,\n",
+    "                                 patch_artist=True,\n",
+    "                                 boxprops=dict(linewidth=0.5),\n",
+    "                                 whiskerprops=dict(linewidth=0.5),\n",
+    "                                 medianprops=dict(linewidth=1, color=\"bisque\"),\n",
+    "                                 capprops=dict(linewidth=0.5))\n",
+    "\n",
+    "    artists[\"boxes\"][0].set_facecolor(group_colors[\"deaf\"])\n",
+    "    artists[\"boxes\"][1].set_facecolor(group_colors[\"hearing\"])\n",
+    "    \n",
+    "    if test_significance:\n",
+    "        p = sps.mannwhitneyu(data[0], data[1], alternative=\"two-sided\").pvalue\n",
+    "        if p < 0.05:\n",
+    "            ax.annotate(\"*\", (0.8, 0.8), xycoords=\"axes fraction\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "call_parameters = [(\"calls_per_10s\", 1, \"Vocal activity\\n[calls/10 s]\", (0, 100), 20),\n",
+    "                   (\"call_rms\", 1, \"Amplitude\\n[dB]\", (0, 54), 12),\n",
+    "                   (\"call_duration\", 1e3, \"Duration\\n[ms]\", (0, 80), 20),\n",
+    "                   (\"f0_mean\", 1e-3, \"Fundamental\\nfrequency [kHz]\", (0, 40), 10),\n",
+    "                   (\"mean_aperiodicity\", 1, \"Aperiodicity\\n[1]\", (0, 1), 0.25)]\n",
+    "\n",
+    "fig_width = rst_total_width\n",
+    "fig_height = 0.7 * fig_width\n",
+    "left_margin, right_margin = 0.75, 0.1\n",
+    "bottom_margin, top_margin = 0.45, 0.45\n",
+    "spacing = 0.2\n",
+    "boxplot_extra_spacing = 0.1\n",
+    "\n",
+    "row_height = (fig_height - bottom_margin - top_margin - (len(call_parameters) - 1) * spacing) / len(call_parameters)\n",
+    "trajectory_width = 3.2\n",
+    "boxplot_width = 0.2\n",
+    "histogram_width = (fig_width - left_margin - right_margin - trajectory_width - boxplot_width - 4 * spacing - boxplot_extra_spacing) / 3\n",
+    "\n",
+    "fig = plt.figure(figsize=(fig_width, fig_height))\n",
+    "\n",
+    "roman_numerals = [\"I\", \"II\", \"III\", \"IV\", \"V\"]\n",
+    "\n",
+    "for i, (parameter, scale, y_label, y_range, y_tick_interval) in enumerate(reversed(call_parameters)):\n",
+    "    early_ax = \\\n",
+    "        fig.add_axes([left_margin / fig_width,\n",
+    "                      (bottom_margin + i * (row_height + spacing)) / fig_height,\n",
+    "                      histogram_width / fig_width,\n",
+    "                      row_height / fig_height])\n",
+    "    late_ax = \\\n",
+    "        fig.add_axes([(left_margin + histogram_width + spacing) / fig_width,\n",
+    "                      (bottom_margin + i * (row_height + spacing)) / fig_height,\n",
+    "                      histogram_width / fig_width,\n",
+    "                      row_height / fig_height])\n",
+    "    adult_ax = \\\n",
+    "        fig.add_axes([(left_margin + 2 * (histogram_width + spacing)) / fig_width,\n",
+    "                      (bottom_margin + i * (row_height + spacing)) / fig_height,\n",
+    "                      histogram_width / fig_width,\n",
+    "                      row_height / fig_height])\n",
+    "    trajectory_ax = \\\n",
+    "        fig.add_axes([(left_margin + 3 * (histogram_width + spacing)) / fig_width,\n",
+    "                      (bottom_margin + i * (row_height + spacing)) / fig_height,\n",
+    "                      trajectory_width / fig_width,\n",
+    "                      row_height / fig_height])\n",
+    "    boxplot_ax = \\\n",
+    "        fig.add_axes([(left_margin + 3 * histogram_width + trajectory_width + 4 * spacing + boxplot_extra_spacing) / fig_width,\n",
+    "                      (bottom_margin + i * (row_height + spacing)) / fig_height,\n",
+    "                      boxplot_width / fig_width,\n",
+    "                      row_height / fig_height])\n",
+    "    overall_ax = \\\n",
+    "        fig.add_axes([left_margin / fig_width,\n",
+    "                      (bottom_margin + i * (row_height + spacing)) / fig_height,\n",
+    "                      (fig_width - left_margin - right_margin) / fig_width,\n",
+    "                      row_height / fig_height], frameon=False)\n",
+    "    overall_ax.set_xticks([])\n",
+    "    overall_ax.set_yticks([])\n",
+    "    \n",
+    "    if parameter == \"calls_per_10s\":\n",
+    "        df = calls_per_pup_and_session\n",
+    "    else:\n",
+    "        df = pup_calls\n",
+    "    plot_histograms(early_ax,\n",
+    "                    df.query(\"before_deafening\"),\n",
+    "                    parameter, scale, y_range)\n",
+    "    plot_histograms(late_ax,\n",
+    "                    df.query(\"not before_deafening\"),\n",
+    "                    parameter, scale, y_range)\n",
+    "    \n",
+    "    if parameter == \"calls_per_10s\":\n",
+    "        df = calls_per_adult_group_and_session\n",
+    "    else:\n",
+    "        df = adult_calls\n",
+    "    plot_histograms(adult_ax, df,\n",
+    "                    parameter, scale, y_range if parameter != \"calls_per_10s\" else (0, 400))\n",
+    "    \n",
+    "    plot_trajectories(trajectory_ax,\n",
+    "                      params_per_pup_and_week.query(\"not before_deafening and calling_bat_age_in_weeks <= 24\"),\n",
+    "                      parameter, scale)\n",
+    "    \n",
+    "    trajectory_ax.xaxis.set_major_locator(mpl.ticker.MultipleLocator(4))\n",
+    "    trajectory_ax.set_xlim(2, 25)\n",
+    "    trajectory_ax.spines[\"left\"].set_position((\"outward\", 3))\n",
+    "\n",
+    "    plot_boxplots(boxplot_ax, df, parameter, scale, test_significance=True)\n",
+    "    boxplot_ax.set_xlim(0.25, 2.75)\n",
+    "    \n",
+    "    if i == 0:\n",
+    "        trajectory_ax.set_xlabel(\"Week of age\")\n",
+    "        trajectory_ax.xaxis.set_major_formatter(mpl.ticker.FuncFormatter(lambda x, pos: \"{}th\".format(int(x))))\n",
+    "        #trajectory_ax.xaxis.set_label_coords(0.5, -0.4)\n",
+    "    else:\n",
+    "        for ax in [early_ax, late_ax, trajectory_ax]:\n",
+    "            ax.set_xticklabels([])\n",
+    "            \n",
+    "    if i == len(call_parameters) - 1:\n",
+    "        early_ax.set_title(\"a)\\n\\n< 2 weeks\", fontsize=8)\n",
+    "        late_ax.set_title(\"b)\\n\\n2–25 weeks\", fontsize=8)\n",
+    "        adult_ax.set_title(\"c)\\n\\n2 years\", fontsize=8)\n",
+    "        boxplot_ax.set_title(\"e)\\n\\n2 years\", fontsize=8)\n",
+    "        trajectory_ax.set_title(\"d)\\nDevelopment trajectory of median parameter values\\nat 2–25 weeks of age\", fontsize=8)\n",
+    "        early_ax.annotate(\"hearing\",\n",
+    "                         xy=(0, 1), xycoords=\"axes fraction\",\n",
+    "                         xytext=(5, -10), textcoords=\"offset points\",\n",
+    "                         verticalalignment=\"top\", horizontalalignment=\"left\",\n",
+    "                         color=group_colors[\"hearing\"])\n",
+    "        early_ax.annotate(\"deafened\",\n",
+    "                         xy=(0, 1), xycoords=\"axes fraction\",\n",
+    "                         xytext=(5, 0), textcoords=\"offset points\",\n",
+    "                         verticalalignment=\"top\", horizontalalignment=\"left\",\n",
+    "                         color=group_colors[\"deaf\"])\n",
+    "            \n",
+    "    for ax in [early_ax, late_ax, adult_ax, trajectory_ax, boxplot_ax]:\n",
+    "        if ax is adult_ax and parameter == \"calls_per_10s\":\n",
+    "            ax.set_ylim(0, 400)\n",
+    "            ax.yaxis.set_major_locator(mpl.ticker.MultipleLocator(100))\n",
+    "        else:\n",
+    "            ax.set_ylim(y_range[0], y_range[1])\n",
+    "            ax.yaxis.set_major_locator(mpl.ticker.MultipleLocator(y_tick_interval))\n",
+    "\n",
+    "        for spine in [\"right\", \"top\"]:\n",
+    "            ax.spines[spine].set_visible(False)\n",
+    "        ax.spines[\"bottom\"].set_position((\"outward\", 3))\n",
+    "            \n",
+    "    for ax in [late_ax, adult_ax, trajectory_ax, boxplot_ax]:\n",
+    "        if (ax is not adult_ax or parameter != \"calls_per_10s\") \\\n",
+    "                and (ax is not trajectory_ax or parameter != \"calls_per_10s\"):\n",
+    "            ax.set_yticklabels([])\n",
+    "        \n",
+    "    for ax in [early_ax, late_ax, adult_ax, boxplot_ax]:\n",
+    "        ax.set_xticks([])\n",
+    "        ax.spines[\"bottom\"].set_visible(False)\n",
+    "        \n",
+    "    early_ax.set_ylabel(y_label)\n",
+    "    early_ax.yaxis.set_label_coords(-0.32 / histogram_width, 0.5)\n",
+    "    \n",
+    "    overall_ax.set_ylabel(roman_numerals[-i - 1],\n",
+    "                          rotation=\"0\",\n",
+    "                          verticalalignment=\"center\",\n",
+    "                         labelpad=49)\n",
+    "    \n",
+    "fig.savefig(\"../parameters.pdf\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Parameter table\n",
+    "==="
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "per_session_parameters = [(\"call_rms\", 1, \"Amplitude [dB]\"),\n",
+    "                          (\"call_duration\", 1e3, \"Duration [ms]\"),\n",
+    "                          (\"f0_min\", 1e-3, \"Minimum fundamental frequency [kHz]\"),\n",
+    "                          (\"f0_mean\", 1e-3, \"Mean fundamental frequency [kHz]\"),\n",
+    "                          (\"f0_max\", 1e-3, \"Maximum fundamental frequency [kHz]\"),\n",
+    "                          (\"f0_start\", 1e-3, \"Fundamental frequency at call onset [kHz]\"),\n",
+    "                          (\"f0_end\", 1e-3, \"Fundamental frequency at end of call [kHz]\"),\n",
+    "                          (\"f_min\", 1e-3, \"Minimum frequency [kHz]\"),\n",
+    "                          (\"f_max\", 1e-3, \"Maximum frequency [kHz]\"),\n",
+    "                          (\"bandwidth\", 1e-3, \"Bandwidth [kHz]\"),\n",
+    "                          (\"spectral_centroid\", 1e-3, \"Spectral centroid frequency [kHz]\"),\n",
+    "                          (\"mean_aperiodicity\", 1, \"Aperiodicity [1]\")]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "adult_calls.loc[:, \"bandwidth\"] = adult_calls[\"f_max\"] - adult_calls[\"f_min\"]\n",
+    "def q25(series):\n",
+    "    return series.quantile(0.25)\n",
+    "def q75(series):\n",
+    "    return series.quantile(0.75)\n",
+    "\n",
+    "adult_groups = adult_calls[[p[0] for p in per_session_parameters] + [\"group\"]].groupby(\"group\")\n",
+    "adult_summary = adult_groups.agg([q25, \"median\", q75])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "per_session_parameters.insert(0, (\"calls_per_10s\", 1, \"Vocal activity [calls/10 s]\"))\n",
+    "adult_activity_summary = calls_per_adult_group_and_session.groupby(\"group\")[[\"calls_per_10s\"]].agg([q25, \"median\", q75])\n",
+    "adult_summary = pd.concat([adult_activity_summary, adult_summary], axis=1)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fout = io.StringIO()\n",
+    "print(\"<!DOCTYPE html><body><table>\", file=fout)\n",
+    "print(\"<tr><th></th><th colspan='3'>Hearing (N={})</th><th colspan='3'>Deaf (N={})</th></tr>\".format(*adult_calls[\"group\"].value_counts().loc[[\"hearing\", \"deaf\"]]), file=fout)\n",
+    "print(\"<tr><th></th>{}<th>p &lt; 0.05 (Mann–Whitney)</th></tr>\".format(\"<th>Q25</th><th>Q50</th><th>Q75</th>\" * 2), file=fout)\n",
+    "\n",
+    "for parameter_name, scale, parameter_description in per_session_parameters:\n",
+    "    print(\"<tr><th>{}</th>\".format(parameter_description), file=fout, end=\"\")\n",
+    "    test_data = []\n",
+    "    for group in [\"hearing\", \"deaf\"]:\n",
+    "        for statistic in [\"q25\", \"median\", \"q75\"]:\n",
+    "            print(\"<td>{:.2f}</td>\".format(adult_summary.loc[group, (parameter_name, statistic)] * scale), file=fout, end=\"\")\n",
+    "        if parameter_name != \"calls_per_10s\":\n",
+    "            test_data.append(adult_calls.loc[adult_calls[\"group\"] == group, parameter_name].values)\n",
+    "    if not test_data:\n",
+    "        print(\"<td>(see plot)</td>\", file=fout)\n",
+    "    else:\n",
+    "        p = sps.mannwhitneyu(test_data[0], test_data[1], alternative=\"two-sided\").pvalue\n",
+    "        if p < 0.05:\n",
+    "            print(\"<td>*</td>\".format(p), file=fout)\n",
+    "        else:\n",
+    "            print(\"<td>{:.2f}</td>\".format(p), file=fout)\n",
+    "    print(\"</tr>\", file=fout)\n",
+    "print(\"</table></body>\", file=fout)\n",
+    "\n",
+    "parameters_table = fout.getvalue()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ipd.display(ipd.HTML(parameters_table))\n",
+    "with open(\"../parameters_table.html\", \"wt\") as fout:\n",
+    "    print(parameters_table, file=fout, end=\"\")"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.8.5"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}

vocalisations/scripts/figures.ipynb

@@ -0,0 +1,714 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib inline"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import IPython.display as ipd\n",
+    "import sys\n",
+    "import io\n",
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "import matplotlib as mpl\n",
+    "import matplotlib.pyplot as plt\n",
+    "import scipy.stats as sps"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Global information\n",
+    "==="
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mpl.rcParams[\"font.family\"] = \"TeX Gyre Pagella\"\n",
+    "mpl.rcParams[\"font.size\"] = 8\n",
+    "rst_column_width = 3.3\n",
+    "rst_total_width = 7"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "deaf_pups = {\"b3\", \"b5\", \"b8\"}\n",
+    "hearing_pups = {\"b2\", \"b4\", \"b7\"}\n",
+    "all_pups = deaf_pups.union(hearing_pups)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "group_colors = dict(hearing=\"#7bb992\", deaf=\"#5f0f00\")\n",
+    "bat_colors = {**{bat: group_colors[\"deaf\"] for bat in deaf_pups},\n",
+    "              **{bat: group_colors[\"hearing\"] for bat in hearing_pups}}\n",
+    "bat_markers = dict(zip(sorted(all_pups), [\"P\", \"o\", \"*\", \"s\", \"X\", \"D\"]))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Load data for juveniles\n",
+    "==="
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pup_sessions = pd.read_csv(\"../pup_sessions.csv\", index_col=0, parse_dates=[\"start_time\"],\n",
+    "                           dtype=dict(before_deafening=np.bool))\n",
+    "pup_calls = pd.read_csv(\"../pup_calls.csv\", index_col=[0, 1], parse_dates=[\"start_time\"], na_values=[\"--\"],\n",
+    "                        dtype=dict(calling_bat_deaf=np.bool, calling_bat_mother=np.bool, before_deafening=np.bool))\n",
+    "\n",
+    "# filter out mother vocalisations\n",
+    "pup_calls = pup_calls[~pup_calls[\"calling_bat_mother\"]].reset_index()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "bat_birthdates = pd.DataFrame([(\"b2\", \"2017/01/26\"),\n",
+    "                               (\"b4\", \"2017/01/30\"),\n",
+    "                               (\"b7\", \"2017/02/08\"),\n",
+    "                               (\"b3\", \"2017/01/29\"),\n",
+    "                               (\"b5\", \"2017/02/02\"),\n",
+    "                               (\"b8\", \"2017/02/08\")], columns=[\"bat\", \"birthdate\"])\n",
+    "bat_birthdates[\"birthdate\"] = pd.to_datetime(bat_birthdates[\"birthdate\"])\n",
+    "bat_birthdates.set_index(\"bat\", inplace=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# determine ages (in full days) of the calling bats\n",
+    "\n",
+    "birth_dates = bat_birthdates.loc[pup_calls[\"calling_bat\"]][\"birthdate\"].reset_index(drop=True)\n",
+    "ages = pup_calls[\"start_time\"] - birth_dates\n",
+    "pup_calls[\"calling_bat_age_in_days\"] = ages.dt.days\n",
+    "pup_calls[\"calling_bat_age_in_weeks\"] = ages.dt.days // 7\n",
+    "pup_calls[\"call_rms\"] = pup_calls[\"call_rms\"] + 55"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Load data for adults\n",
+    "==="
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def single_value(values):\n",
+    "    value_set = set(values)\n",
+    "    if len(value_set) != 1:\n",
+    "        raise ValueError(\"non-unity set\")\n",
+    "    return next(iter(value_set))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "adult_calls = pd.read_csv(\"../adult_calls.csv\",\n",
+    "                          parse_dates=[\"start_time\"],\n",
+    "                          na_values=[\"--\"],\n",
+    "                          dtype=dict(calling_bat_deaf=np.bool, session_id=np.int),\n",
+    "                          index_col=[\"session_id\", \"call_id\"])\n",
+    "adult_calls[\"recording_day\"] = adult_calls[\"start_time\"].dt.date\n",
+    "adult_calls[\"call_rms\"] -= adult_calls[\"call_rms\"].min()\n",
+    "adult_calls[\"group\"] = np.where(adult_calls[\"calling_bat_deaf\"], \"deaf\", \"hearing\")\n",
+    "\n",
+    "adult_sessions = pd.read_csv(\"../adult_sessions.csv\",\n",
+    "                             parse_dates=[\"start_time\"],\n",
+    "                             dtype=dict(group=\"category\"),\n",
+    "                             index_col=0)\n",
+    "\n",
+    "adult_recording_days = adult_sessions[[\"group\"]].groupby(adult_sessions[\"start_time\"].dt.date).agg(single_value)\n",
+    "adult_recording_days.index.name = \"recording_day\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Data to be dropped\n",
+    "==="
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Possible echolocation calls (3 ms or shorter)\n",
+    "---"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pup_calls = pup_calls[pup_calls[\"call_duration\"] > 3e-3]\n",
+    "adult_calls = adult_calls[adult_calls[\"call_duration\"] > 3e-3]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Sessions shorter than 20 minutes\n",
+    "---"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "session_gaps = pup_sessions.sort_values(\"start_time\")[\"start_time\"].diff()\n",
+    "new_index = session_gaps.index.insert(0, -1)\n",
+    "new_index = new_index[:-1]\n",
+    "session_gaps.index = new_index\n",
+    "session_gaps.drop(-1, inplace=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ids_to_drop = session_gaps[session_gaps < pd.Timedelta(20, \"m\")].index"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pup_calls = pup_calls[~pup_calls[\"session_id\"].isin(ids_to_drop)]\n",
+    "pup_sessions = pup_sessions.drop(ids_to_drop)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Last two recordings with mother for pups that have more than others\n",
+    "---"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sorted_pup_mother_sessions = pup_sessions[pup_sessions[\"animal2\"].str.startswith(\"m\")].sort_values(\"start_time\")\n",
+    "session_counts_with_mothers = sorted_pup_mother_sessions[\"animal1\"].value_counts()\n",
+    "num_sessions_to_drop = session_counts_with_mothers - session_counts_with_mothers.min()\n",
+    "\n",
+    "ids_to_drop = set()\n",
+    "for bat, drop_count in num_sessions_to_drop.iteritems():\n",
+    "    if drop_count == 0:\n",
+    "        continue\n",
+    "    this_bat_sessions = sorted_pup_mother_sessions[sorted_pup_mother_sessions[\"animal1\"] == bat]\n",
+    "    ids_to_drop = ids_to_drop.union(this_bat_sessions.tail(drop_count).index)\n",
+    "\n",
+    "pup_calls = pup_calls[~pup_calls[\"session_id\"].isin(ids_to_drop)]\n",
+    "pup_sessions = pup_sessions.drop(index=ids_to_drop)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Compute activity statistics (separately for pups before and after deafening)\n",
+    "==="
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "calls_per_pup_and_session = pup_calls.reset_index().groupby([\"calling_bat\", \"session_id\"])[[\"calling_bat\"]].count() / (20*6) # per 10 seconds\n",
+    "calls_per_pup_and_session.columns = [\"calls_per_10s\"]\n",
+    "calls_per_pup_and_session.reset_index(inplace=True)\n",
+    "\n",
+    "session_dates = pup_sessions.loc[calls_per_pup_and_session[\"session_id\"]][\"start_time\"]\n",
+    "birth_dates = bat_birthdates.loc[calls_per_pup_and_session[\"calling_bat\"]][\"birthdate\"]\n",
+    "calls_per_pup_and_session[\"calling_bat_age_in_days\"] = \\\n",
+    "    (session_dates.reset_index(drop=True) - birth_dates.reset_index(drop=True)).dt.days\n",
+    "calls_per_pup_and_session[\"calling_bat_age_in_weeks\"] = \\\n",
+    "    calls_per_pup_and_session[\"calling_bat_age_in_days\"] // 7\n",
+    "calls_per_pup_and_session = pd.merge(calls_per_pup_and_session, pup_sessions[[\"before_deafening\"]],\n",
+    "                                     left_on=\"session_id\", right_index=True,)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dfs = []\n",
+    "\n",
+    "for before_deafening in [False, True]:\n",
+    "    mask = pup_calls[\"before_deafening\"]\n",
+    "    if not before_deafening:\n",
+    "        mask = ~mask\n",
+    "    gr = pup_calls[mask].groupby([\"calling_bat\", \"calling_bat_age_in_weeks\"])\n",
+    "    gr = gr[[\"call_duration\", \"mean_aperiodicity\", \"f0_mean\", \"call_rms\"]]\n",
+    "    params_per_pup_and_week = gr.median()\n",
+    "    params_per_pup_and_week[\"calls_this_week\"] = gr.size()\n",
+    "    \n",
+    "    mask = calls_per_pup_and_session[\"before_deafening\"]\n",
+    "    if not before_deafening:\n",
+    "        mask = ~mask\n",
+    "    gr = calls_per_pup_and_session[mask].groupby([\"calling_bat\", \"calling_bat_age_in_weeks\"])\n",
+    "    params_per_pup_and_week[\"calls_per_10s\"] = gr[\"calls_per_10s\"].mean()\n",
+    "    params_per_pup_and_week.reset_index(inplace=True)\n",
+    "    params_per_pup_and_week[\"before_deafening\"] = before_deafening\n",
+    "    \n",
+    "    dfs.append(params_per_pup_and_week)\n",
+    "    \n",
+    "params_per_pup_and_week = pd.concat(dfs)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "gr = adult_calls.groupby([\"group\", \"session_id\"], observed=True)\n",
+    "calls_per_adult_group_and_session = gr.size().to_frame(name=\"calls_per_10s\").reset_index()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Quantities of extracted data\n",
+    "==="
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(\"Total calls by pups: {}\\nCalls by deaf pups: {} ({:.1f} %)\\nCalls by hearing pups: {} ({:.1f} %)\".format(\n",
+    "        len(pup_calls),\n",
+    "        np.sum(pup_calls[\"calling_bat_deaf\"]),\n",
+    "        np.sum(pup_calls[\"calling_bat_deaf\"]) / len(pup_calls) * 100,\n",
+    "        np.sum(~pup_calls[\"calling_bat_deaf\"]),\n",
+    "        np.sum(~pup_calls[\"calling_bat_deaf\"] / len(pup_calls) * 100)))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(\"Total calls by adults: {}\\nCalls by deaf adults: {} ({:.1f} %)\\nCalls by hearing adults: {} ({:.1f} %)\".format(\n",
+    "        len(adult_calls),\n",
+    "        np.sum(adult_calls[\"calling_bat_deaf\"]),\n",
+    "        np.sum(adult_calls[\"calling_bat_deaf\"]) / len(adult_calls) * 100,\n",
+    "        np.sum(~adult_calls[\"calling_bat_deaf\"]),\n",
+    "        np.sum(~adult_calls[\"calling_bat_deaf\"] / len(adult_calls) * 100)))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Main figure\n",
+    "==="
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def plot_trajectories(ax, values, parameter, scale):\n",
+    "    line_artists = {}\n",
+    "    for bat, df in values.groupby(\"calling_bat\"):\n",
+    "        bat_color, bat_marker = bat_colors[bat], bat_markers[bat]\n",
+    "        line_artist, = ax.plot(df[\"calling_bat_age_in_weeks\"] + 1, df[parameter] * scale,\n",
+    "                               c=bat_color, linewidth=0.8)\n",
+    "        ax.plot(df[\"calling_bat_age_in_weeks\"] + 1, df[parameter] * scale,\n",
+    "                c=bat_color, marker=bat_marker, markersize=5, clip_on=False,\n",
+    "                markerfacecolor=\"none\", linestyle=\"\")\n",
+    "        line_artists[bat] = line_artist\n",
+    "    return line_artists"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def plot_histograms(ax, values, parameter, scale, y_range):\n",
+    "    data = {}\n",
+    "\n",
+    "    if \"group\" in values:\n",
+    "        for group, hearing in [(\"deaf\", False), (\"hearing\", True)]:\n",
+    "            data[hearing] = (values.loc[values[\"group\"] == group, parameter] * scale).values\n",
+    "    else:\n",
+    "        for hearing, df in values.groupby(values[\"calling_bat\"].isin(hearing_pups)):\n",
+    "            data[hearing] = (df[parameter] * scale).values\n",
+    "\n",
+    "    ax.hist([data[True], data[False]], density=True,\n",
+    "            range=y_range, bins=15,\n",
+    "            color=[group_colors[\"hearing\"], group_colors[\"deaf\"]],\n",
+    "            orientation=\"horizontal\", rwidth=0.85)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def plot_boxplots(ax, values, parameter, scale, test_significance=False):\n",
+    "    data = []\n",
+    "    for group in [\"deaf\", \"hearing\"]:\n",
+    "        data.append(scale * values.loc[values[\"group\"] == group, parameter].dropna().values)\n",
+    "\n",
+    "    artists = boxplot_ax.boxplot(data, vert=True, manage_ticks=False, showfliers=False, widths=0.75,\n",
+    "                                 patch_artist=True,\n",
+    "                                 boxprops=dict(linewidth=0.5),\n",
+    "                                 whiskerprops=dict(linewidth=0.5),\n",
+    "                                 medianprops=dict(linewidth=1, color=\"bisque\"),\n",
+    "                                 capprops=dict(linewidth=0.5))\n",
+    "\n",
+    "    artists[\"boxes\"][0].set_facecolor(group_colors[\"deaf\"])\n",
+    "    artists[\"boxes\"][1].set_facecolor(group_colors[\"hearing\"])\n",
+    "    \n",
+    "    if test_significance:\n",
+    "        p = sps.mannwhitneyu(data[0], data[1], alternative=\"two-sided\").pvalue\n",
+    "        if p < 0.05:\n",
+    "            ax.annotate(\"*\", (0.8, 0.8), xycoords=\"axes fraction\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "call_parameters = [(\"calls_per_10s\", 1, \"Vocal activity\\n[calls/10 s]\", (0, 100), 20),\n",
+    "                   (\"call_rms\", 1, \"Amplitude\\n[dB]\", (0, 54), 12),\n",
+    "                   (\"call_duration\", 1e3, \"Duration\\n[ms]\", (0, 80), 20),\n",
+    "                   (\"f0_mean\", 1e-3, \"Fundamental\\nfrequency [kHz]\", (0, 40), 10),\n",
+    "                   (\"mean_aperiodicity\", 1, \"Aperiodicity\\n[1]\", (0, 1), 0.25)]\n",
+    "\n",
+    "fig_width = rst_total_width\n",
+    "fig_height = 0.7 * fig_width\n",
+    "left_margin, right_margin = 0.75, 0.1\n",
+    "bottom_margin, top_margin = 0.45, 0.45\n",
+    "spacing = 0.2\n",
+    "boxplot_extra_spacing = 0.1\n",
+    "\n",
+    "row_height = (fig_height - bottom_margin - top_margin - (len(call_parameters) - 1) * spacing) / len(call_parameters)\n",
+    "trajectory_width = 3.2\n",
+    "boxplot_width = 0.2\n",
+    "histogram_width = (fig_width - left_margin - right_margin - trajectory_width - boxplot_width - 4 * spacing - boxplot_extra_spacing) / 3\n",
+    "\n",
+    "fig = plt.figure(figsize=(fig_width, fig_height))\n",
+    "\n",
+    "roman_numerals = [\"I\", \"II\", \"III\", \"IV\", \"V\"]\n",
+    "\n",
+    "for i, (parameter, scale, y_label, y_range, y_tick_interval) in enumerate(reversed(call_parameters)):\n",
+    "    early_ax = \\\n",
+    "        fig.add_axes([left_margin / fig_width,\n",
+    "                      (bottom_margin + i * (row_height + spacing)) / fig_height,\n",
+    "                      histogram_width / fig_width,\n",
+    "                      row_height / fig_height])\n",
+    "    late_ax = \\\n",
+    "        fig.add_axes([(left_margin + histogram_width + spacing) / fig_width,\n",
+    "                      (bottom_margin + i * (row_height + spacing)) / fig_height,\n",
+    "                      histogram_width / fig_width,\n",
+    "                      row_height / fig_height])\n",
+    "    adult_ax = \\\n",
+    "        fig.add_axes([(left_margin + 2 * (histogram_width + spacing)) / fig_width,\n",
+    "                      (bottom_margin + i * (row_height + spacing)) / fig_height,\n",
+    "                      histogram_width / fig_width,\n",
+    "                      row_height / fig_height])\n",
+    "    trajectory_ax = \\\n",
+    "        fig.add_axes([(left_margin + 3 * (histogram_width + spacing)) / fig_width,\n",
+    "                      (bottom_margin + i * (row_height + spacing)) / fig_height,\n",
+    "                      trajectory_width / fig_width,\n",
+    "                      row_height / fig_height])\n",
+    "    boxplot_ax = \\\n",
+    "        fig.add_axes([(left_margin + 3 * histogram_width + trajectory_width + 4 * spacing + boxplot_extra_spacing) / fig_width,\n",
+    "                      (bottom_margin + i * (row_height + spacing)) / fig_height,\n",
+    "                      boxplot_width / fig_width,\n",
+    "                      row_height / fig_height])\n",
+    "    overall_ax = \\\n",
+    "        fig.add_axes([left_margin / fig_width,\n",
+    "                      (bottom_margin + i * (row_height + spacing)) / fig_height,\n",
+    "                      (fig_width - left_margin - right_margin) / fig_width,\n",
+    "                      row_height / fig_height], frameon=False)\n",
+    "    overall_ax.set_xticks([])\n",
+    "    overall_ax.set_yticks([])\n",
+    "    \n",
+    "    if parameter == \"calls_per_10s\":\n",
+    "        df = calls_per_pup_and_session\n",
+    "    else:\n",
+    "        df = pup_calls\n",
+    "    plot_histograms(early_ax,\n",
+    "                    df.query(\"before_deafening\"),\n",
+    "                    parameter, scale, y_range)\n",
+    "    plot_histograms(late_ax,\n",
+    "                    df.query(\"not before_deafening\"),\n",
+    "                    parameter, scale, y_range)\n",
+    "    \n",
+    "    if parameter == \"calls_per_10s\":\n",
+    "        df = calls_per_adult_group_and_session\n",
+    "    else:\n",
+    "        df = adult_calls\n",
+    "    plot_histograms(adult_ax, df,\n",
+    "                    parameter, scale, y_range if parameter != \"calls_per_10s\" else (0, 400))\n",
+    "    \n",
+    "    plot_trajectories(trajectory_ax,\n",
+    "                      params_per_pup_and_week.query(\"not before_deafening and calling_bat_age_in_weeks <= 24\"),\n",
+    "                      parameter, scale)\n",
+    "    \n",
+    "    trajectory_ax.xaxis.set_major_locator(mpl.ticker.MultipleLocator(4))\n",
+    "    trajectory_ax.set_xlim(2, 25)\n",
+    "    trajectory_ax.spines[\"left\"].set_position((\"outward\", 3))\n",
+    "\n",
+    "    plot_boxplots(boxplot_ax, df, parameter, scale, test_significance=True)\n",
+    "    boxplot_ax.set_xlim(0.25, 2.75)\n",
+    "    \n",
+    "    if i == 0:\n",
+    "        trajectory_ax.set_xlabel(\"Week of age\")\n",
+    "        trajectory_ax.xaxis.set_major_formatter(mpl.ticker.FuncFormatter(lambda x, pos: \"{}th\".format(int(x))))\n",
+    "        #trajectory_ax.xaxis.set_label_coords(0.5, -0.4)\n",
+    "    else:\n",
+    "        for ax in [early_ax, late_ax, trajectory_ax]:\n",
+    "            ax.set_xticklabels([])\n",
+    "            \n",
+    "    if i == len(call_parameters) - 1:\n",
+    "        early_ax.set_title(\"a)\\n\\n< 2 weeks\", fontsize=8)\n",
+    "        late_ax.set_title(\"b)\\n\\n2–25 weeks\", fontsize=8)\n",
+    "        adult_ax.set_title(\"c)\\n\\n2 years\", fontsize=8)\n",
+    "        boxplot_ax.set_title(\"e)\\n\\n2 years\", fontsize=8)\n",
+    "        trajectory_ax.set_title(\"d)\\nDevelopment trajectory of median parameter values\\nat 2–25 weeks of age\", fontsize=8)\n",
+    "        early_ax.annotate(\"hearing\",\n",
+    "                         xy=(0, 1), xycoords=\"axes fraction\",\n",
+    "                         xytext=(5, -10), textcoords=\"offset points\",\n",
+    "                         verticalalignment=\"top\", horizontalalignment=\"left\",\n",
+    "                         color=group_colors[\"hearing\"])\n",
+    "        early_ax.annotate(\"deafened\",\n",
+    "                         xy=(0, 1), xycoords=\"axes fraction\",\n",
+    "                         xytext=(5, 0), textcoords=\"offset points\",\n",
+    "                         verticalalignment=\"top\", horizontalalignment=\"left\",\n",
+    "                         color=group_colors[\"deaf\"])\n",
+    "            \n",
+    "    for ax in [early_ax, late_ax, adult_ax, trajectory_ax, boxplot_ax]:\n",
+    "        if ax is adult_ax and parameter == \"calls_per_10s\":\n",
+    "            ax.set_ylim(0, 400)\n",
+    "            ax.yaxis.set_major_locator(mpl.ticker.MultipleLocator(100))\n",
+    "        else:\n",
+    "            ax.set_ylim(y_range[0], y_range[1])\n",
+    "            ax.yaxis.set_major_locator(mpl.ticker.MultipleLocator(y_tick_interval))\n",
+    "\n",
+    "        for spine in [\"right\", \"top\"]:\n",
+    "            ax.spines[spine].set_visible(False)\n",
+    "        ax.spines[\"bottom\"].set_position((\"outward\", 3))\n",
+    "            \n",
+    "    for ax in [late_ax, adult_ax, trajectory_ax, boxplot_ax]:\n",
+    "        if (ax is not adult_ax or parameter != \"calls_per_10s\") \\\n",
+    "                and (ax is not trajectory_ax or parameter != \"calls_per_10s\"):\n",
+    "            ax.set_yticklabels([])\n",
+    "        \n",
+    "    for ax in [early_ax, late_ax, adult_ax, boxplot_ax]:\n",
+    "        ax.set_xticks([])\n",
+    "        ax.spines[\"bottom\"].set_visible(False)\n",
+    "        \n",
+    "    early_ax.set_ylabel(y_label)\n",
+    "    early_ax.yaxis.set_label_coords(-0.32 / histogram_width, 0.5)\n",
+    "    \n",
+    "    overall_ax.set_ylabel(roman_numerals[-i - 1],\n",
+    "                          rotation=\"0\",\n",
+    "                          verticalalignment=\"center\",\n",
+    "                         labelpad=49)\n",
+    "    \n",
+    "fig.savefig(\"../parameters.pdf\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Parameter table\n",
+    "==="
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "per_session_parameters = [(\"call_rms\", 1, \"Amplitude [dB]\"),\n",
+    "                          (\"call_duration\", 1e3, \"Duration [ms]\"),\n",
+    "                          (\"f0_min\", 1e-3, \"Minimum fundamental frequency [kHz]\"),\n",
+    "                          (\"f0_mean\", 1e-3, \"Mean fundamental frequency [kHz]\"),\n",
+    "                          (\"f0_max\", 1e-3, \"Maximum fundamental frequency [kHz]\"),\n",
+    "                          (\"f0_start\", 1e-3, \"Fundamental frequency at call onset [kHz]\"),\n",
+    "                          (\"f0_end\", 1e-3, \"Fundamental frequency at end of call [kHz]\"),\n",
+    "                          (\"f_min\", 1e-3, \"Minimum frequency [kHz]\"),\n",
+    "                          (\"f_max\", 1e-3, \"Maximum frequency [kHz]\"),\n",
+    "                          (\"bandwidth\", 1e-3, \"Bandwidth [kHz]\"),\n",
+    "                          (\"spectral_centroid\", 1e-3, \"Spectral centroid frequency [kHz]\"),\n",
+    "                          (\"mean_aperiodicity\", 1, \"Aperiodicity [1]\")]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "adult_calls.loc[:, \"bandwidth\"] = adult_calls[\"f_max\"] - adult_calls[\"f_min\"]\n",
+    "def q25(series):\n",
+    "    return series.quantile(0.25)\n",
+    "def q75(series):\n",
+    "    return series.quantile(0.75)\n",
+    "\n",
+    "adult_groups = adult_calls[[p[0] for p in per_session_parameters] + [\"group\"]].groupby(\"group\")\n",
+    "adult_summary = adult_groups.agg([q25, \"median\", q75])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "per_session_parameters.insert(0, (\"calls_per_10s\", 1, \"Vocal activity [calls/10 s]\"))\n",
+    "adult_activity_summary = calls_per_adult_group_and_session.groupby(\"group\")[[\"calls_per_10s\"]].agg([q25, \"median\", q75])\n",
+    "adult_summary = pd.concat([adult_activity_summary, adult_summary], axis=1)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fout = io.StringIO()\n",
+    "print(\"<!DOCTYPE html><body><table>\", file=fout)\n",
+    "print(\"<tr><th></th><th colspan='3'>Hearing (N={})</th><th colspan='3'>Deaf (N={})</th></tr>\".format(*adult_calls[\"group\"].value_counts().loc[[\"hearing\", \"deaf\"]]), file=fout)\n",
+    "print(\"<tr><th></th>{}<th>p &lt; 0.05 (Mann–Whitney)</th></tr>\".format(\"<th>Q25</th><th>Q50</th><th>Q75</th>\" * 2), file=fout)\n",
+    "\n",
+    "for parameter_name, scale, parameter_description in per_session_parameters:\n",
+    "    print(\"<tr><th>{}</th>\".format(parameter_description), file=fout, end=\"\")\n",
+    "    test_data = []\n",
+    "    for group in [\"hearing\", \"deaf\"]:\n",
+    "        for statistic in [\"q25\", \"median\", \"q75\"]:\n",
+    "            print(\"<td>{:.2f}</td>\".format(adult_summary.loc[group, (parameter_name, statistic)] * scale), file=fout, end=\"\")\n",
+    "        if parameter_name != \"calls_per_10s\":\n",
+    "            test_data.append(adult_calls.loc[adult_calls[\"group\"] == group, parameter_name].values)\n",
+    "    if not test_data:\n",
+    "        print(\"<td>(see plot)</td>\", file=fout)\n",
+    "    else:\n",
+    "        p = sps.mannwhitneyu(test_data[0], test_data[1], alternative=\"two-sided\").pvalue\n",
+    "        if p < 0.05:\n",
+    "            print(\"<td>*</td>\".format(p), file=fout)\n",
+    "        else:\n",
+    "            print(\"<td>{:.2f}</td>\".format(p), file=fout)\n",
+    "    print(\"</tr>\", file=fout)\n",
+    "print(\"</table></body>\", file=fout)\n",
+    "\n",
+    "parameters_table = fout.getvalue()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ipd.display(ipd.HTML(parameters_table))\n",
+    "with open(\"../parameters_table.html\", \"wt\") as fout:\n",
+    "    print(parameters_table, file=fout, end=\"\")"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.8.5"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}

vocalisations/scripts_for_reference/extract_adult_call_parameters.m

@@ -0,0 +1,193 @@
+function result = extract_adult_call_parameters(recbuf)
+
+% fixed variables
+fs=192e3;
+env_flow=500;
+
+pk_thrs=-45;    % dB
+pk_dist=0.005;    % in seconds (changed from 0.02 in Meike's extraction script)
+dur_thrs=-20;   % dB
+freq_thrs=-20;  % dB
+
+hpf_display=50e3;
+nfft = 128;
+
+[eb,ea]=butter(2,2*env_flow/fs);
+[hpb,hpa] = butter(8,2*hpf_display/fs,'high'); % 50 (!) kHz high pass
+
+% parameters for envelope (SAM) analysis
+fftpts=2^16;
+env_f=linspace(0,fs,fftpts);
+amf_min=50;
+amf_max=800;
+amf_ind1=min(find(env_f>amf_min));
+amf_ind2=max(find(env_f<amf_max));
+amf_inds=[amf_ind1 amf_ind2];
+amfs=env_f(amf_inds(1):amf_inds(2));
+am_mindur=10e-3;
+am_pk_thrs=-90;    % dB
+am_pk_prom=9;  % dB
+
+recbuf=double(recbuf(:,1:16))/32767;
+% highpass filter
+recbuf=filtfilt(hpb,hpa,recbuf);
+
+aud_time=(0:size(recbuf,1)-1)'/fs;
+
+% search for calls
+sumblockbuf=sum(abs(recbuf),2);
+
+% smooth the envelope
+env_dB=20*log10(max(filtfilt(eb,ea,sumblockbuf),1e-4));
+
+% find calls
+t=[0:length(env_dB)-1]/fs;
+[pk_val pk_pos]=findpeaks(env_dB,'minpeakheight',pk_thrs,'minpeakdistance',round(pk_dist*fs));
+cn=0;
+
+if isempty(pk_pos)
+    return;
+end
+
+% call cutting and saving
+call_num=1;
+while call_num<=length(pk_pos),
+    % determine start and stop of each call as -x dB duration
+    env_dB=env_dB-env_dB(pk_pos(call_num));
+    % find start of call
+    start=max(find(env_dB(1:pk_pos(call_num))<dur_thrs));
+    stop=min(find(env_dB(pk_pos(call_num):end)<dur_thrs))+pk_pos(call_num);
+    
+    if 1-isempty(stop) & 1-isempty(start)
+        cn=cn+1;
+        
+        call_dur=(stop-start)/fs;
+        
+        % SET UP RESULTS STRUCTURE
+        result.call(cn).call_start_sample=pk_pos(cn);
+        result.call(cn).call_dur=call_dur;
+        
+        call=recbuf(start:stop,1:16);
+        call_levels=20*log10(std(call));
+        [call_level_max call_level_max_pos]=max(call_levels);
+        xcall=call(:,call_level_max_pos);
+        xct=(0:length(xcall)-1)'/fs*1e3;
+        result.call(cn).call_levels=call_levels;
+        result.call(cn).loudest_call=xcall;
+        result.call(cn).loudest_mic=call_level_max_pos;
+        
+        % --- call frequency parameter ----
+        % spectrogram method
+        spects = pwelch(xcall,256,250,nfft,fs);
+        spectsdB = 10*log10(spects);
+        [ampl_pf pos_pf]= max(spectsdB);
+        r = find(spectsdB >= ampl_pf+freq_thrs);
+        faxis = (0:1/nfft:1/2).*fs;
+        PF = faxis(pos_pf);    % peak frequency
+        Fmin = faxis(r(1));    % min frequency
+        Fmax= faxis(r(end));  % max frequency
+        BW= Fmax - Fmin;    % bandwidth
+        
+        % --- spectral centroid frequency
+        SCF = sum(spects.*faxis')/sum(spects);
+        SCF_dB = sum(spectsdB.*faxis')/sum(spectsdB);
+        
+        result.call(cn).PF=PF;
+        result.call(cn).Fmin=Fmin;
+        result.call(cn).Fmax=Fmax;
+        result.call(cn).BW=BW;
+        result.call(cn).SCF=SCF;
+        result.call(cn).SCF_dB=SCF_dB;
+        
+        
+        % use yin to calculate f0 as a function of time
+        p.sr = fs;
+        
+        R=yin(xcall,p,0);
+        f0s=R.f0;
+        nframes=length(f0s);
+        f0ts=(1:nframes)/(fs/R.hop);
+        aperiodicity=R.ap0;
+        
+        % correct for f0 jumps
+        fcin=[f0ts' f0s' aperiodicity'];
+        fdiff=abs(diff(fcin(:,2)));
+        finds=find(fdiff>9000);
+        if 1-isempty(finds) % corrections are necessary
+            f0_corrected=1;
+            f0_corr=fcin(:,2);
+            % check start frequency of the irregularity
+            if f0_corr(finds(1))>30e3;
+                % correct first jump
+                f0_corr(1:finds(1))=f0_corr(1:finds(1))/2;
+                if length(finds)>1,
+                    finds=finds(2:end);
+                end
+            end
+            while length(finds)>1
+                xfinds=finds(1:2);
+                f0_corr(xfinds(1)+1:xfinds(2))=f0_corr(xfinds(1)+1:xfinds(2))/2;
+                if length(finds)>2,
+                    finds=finds(3:end);
+                elseif length(finds)==2,
+                    finds=[];
+                    break
+                else
+                    break
+                end
+            end
+            if length(finds)>0,
+                f0_corr(finds(1)+1:end)=f0_corr(finds(1)+1:end)/2;
+            end
+            %                         figure(10), plot(fcin(:,1),fcin(:,2),'b',fcin(:,1),f0_corr,'r');
+            %corr_ok=input('correction (1) ok or (0) not ok');
+            corr_ok=1;
+            if corr_ok==1,
+                fcin(:,2)=f0_corr;
+            end
+        else
+            f0_corrected=0;
+        end
+        
+        result.call(cn).f0=fcin;
+        result.call(cn).f0_corrected=f0_corrected;
+        
+        % fit f0(t)
+        f0inds=find(1-isnan(fcin(:,2)));
+        f0data=fcin(f0inds,1:2);
+        [f0slope,sfm_amp,sfm_frq,sfm_phase]=fit_f0(f0data,0);
+        
+        result.call(cn).f0_slope=f0slope;
+        result.call(cn).sfm_amp=sfm_amp;
+        result.call(cn).sfm_frq=sfm_frq;
+        result.call(cn).sfm_phase=sfm_phase;
+        
+        % quantify AM in call envelope
+        if call_dur>=am_mindur,
+            if length(xcall)<fftpts,
+                lxcall=[xcall;zeros(fftpts-length(xcall),1)];
+            end
+            
+            [modspec,mfs]=pwelch(abs(lxcall),fftpts,0.5*fftpts,fftpts,fs);
+            modspec(1:amf_inds(1))=0;
+            modspec(amf_inds(2):end)=0;
+            modspecdb=10*log10(modspec);
+            [AM_val AM_pos]=findpeaks(modspecdb,'minpeakheight',am_pk_thrs,'MinPeakProminence',am_pk_prom);
+            if isempty(AM_pos),
+                result.call(cn).amf=0;
+            else
+                % select the highest peak
+                [mist,ind]=max(AM_val);
+                maxAMfreq=mfs(AM_pos(ind));
+                result.call(cn).amf=maxAMfreq;
+            end
+            
+        else
+            result.call(cn).amf=0;
+        end
+    end
+    call_num=call_num+1;
+    start_sec = start/fs;
+    stop_sec = stop/fs;
+    disp(sprintf('DONE: call %u, %f-%f s', cn, start_sec, stop_sec))
+end

vocalisations/scripts_for_reference/extract_pup_call_parameters.m

@@ -0,0 +1,222 @@
+clear all
+close all
+
+rootpath='vpl_behav';
+bats = {'vpl_b2b3', 'vpl_b2b5', 'vpl_b2b8', 'vpl_b2mh', 'vpl_b3mt', ...
+	'vpl_b4b3', 'vpl_b4b5', 'vpl_b4b8', 'vpl_b4mb', 'vpl_b5mt', ...
+	'vpl_b7b3', 'vpl_b7b5', 'vpl_b7b8', 'vpl_b7my', 'vpl_b8mp'};
+
+% fixed variables
+fs=192e3;
+env_mad=0.005;   % duration of the boxcar kernel(for moving average) to smooth envelope (s)
+env_flow=500;
+
+pk_thrs=-45;    % dB
+pk_dist=0.02;    % in seconds
+
+dur_thrs=-20;   % dB
+freq_thrs=-20;  % dB
+
+
+
+[hpb,hpa] = butter(8,2*50e3/fs,'high'); % 10 kHz high pass
+
+blockdur=60;
+
+mph=0.002;
+mpd=0.01*fs;
+nfft=128;
+
+for batnum=1:length(bats),
+    filepath=[rootpath 'recordings\' bats{batnum} '\'];
+    d=dir(filepath);
+    wavnames=strvcat(d.name);
+    for wavnum=3:size(wavnames,1);
+        wavname=deblank(wavnames(wavnum,:));
+        if 1-isempty(findstr(wavname,'m1.wav'));
+            wavname
+            cn=0;
+            result=[];
+            savname=[rootpath 'results\' bats{batnum} '\' wavname(1:end-4) '.mat'];
+            fid=fopen(savname);
+            if fid<0,   % this file has not been analysed
+                filesize=wavread([filepath wavname],'size');
+                filename1=wavname(1:end-5);
+                nblocks=ceil(filesize(1)/(blockdur*fs));
+                
+                for bn=1:nblocks,
+                    disp(['working on block ' num2str(bn) ' of ' num2str(nblocks)]);
+                    blockbuf=zeros(round(blockdur*fs),1);
+                    block_start=(bn-1)*round(blockdur*fs)+1;
+                    
+                    if bn < nblocks
+                        block_stop=bn*blockdur*fs;
+                    else
+                        block_stop=filesize(1);
+                    end
+                    
+                    blockbuf=zeros(block_stop-block_start+1,6);
+                    
+                    for micnum=1:6,
+                        micfile=[filename1 num2str(micnum) '.wav'];
+                        blockbuf(:,micnum)=wavread([filepath micfile],[block_start block_stop]);
+                    end
+                    compblockbuf=filtfilt(hpb,hpa,blockbuf);
+                    
+                    % call finder
+                    sumblockbuf=sum(abs(compblockbuf),2);
+                    
+                    % smooth the envelope
+                    [eb,ea]=butter(2,2*env_flow/fs);
+                    env_dB=20*log10(max(filtfilt(eb,ea,sumblockbuf),1e-4));
+                    
+                    %% find calls
+                    t=[0:length(env_dB)-1]/fs;
+                    
+                    
+                    %                     figure(2),
+                    %                     plot(env_dB)
+                    %                     hold on
+                    %                     xs=get(gca,'xlim');
+                    %                     ys=[pk_thrs pk_thrs];
+                    %                     lh=line(xs,ys);
+                    %                     set(lh,'color','r')
+                    
+                    [pk_val pk_pos]=findpeaks(env_dB,'minpeakheight',pk_thrs,'minpeakdistance',round(pk_dist*fs));
+                    
+                    if 1-isempty(pk_pos)
+                        % draw found peaks positions
+                        %                         hold on
+                        %                         plot(pk_pos,pk_val,'r*')
+                        %                         hold off
+                        
+                        % call cutting and saving
+                        for call_num=1:size(pk_pos)
+                            cn=cn+1;
+                            
+                            % determine start and stop of each call as -x dB duration
+                            env_dB=env_dB-env_dB(pk_pos(call_num));
+                            % find start of call
+                            start=max(find(env_dB(1:pk_pos(call_num))<dur_thrs));
+                            stop=min(find(env_dB(pk_pos(call_num):end)<dur_thrs))+pk_pos(call_num);
+                            
+                            if 1-isempty(stop) & 1-isempty(start)
+                                
+                                call_starts=(bn-1)*blockdur*fs+start;
+                                call_dur=(stop-start)/fs;
+                                
+                                result.call(cn).call_start_sample=call_starts;
+                                result.call(cn).call_dur=call_dur;
+                                
+                                call=blockbuf(start:stop,1:6);
+                                call_levels=20*log10(std(call));
+                                [call_level_max call_level_max_pos]=max(call_levels);
+                                xcall=call(:,call_level_max_pos);
+                                %                             figure(3)
+                                %                             subplot(2,1,1)
+                                %                             plot(xcall)
+                                %                             subplot(2,1,2)
+                                %                             spectrogram(xcall,nfft,120,nfft,fs,'yaxis');
+                                
+                                result.call(cn).call_levels=call_levels;
+                                result.call(cn).loudest_call=xcall;
+                                result_call(cn).loudest_mic=call_level_max_pos;
+                                
+                                %% --- call frequency parameter ----
+                                %% spectrogram method
+                                spects = pwelch(xcall,256,250,nfft,fs);
+                                spectsdB = 10*log10(spects);
+                                [ampl_pf pos_pf]= max(spectsdB);
+                                r = find(spectsdB >= ampl_pf+freq_thrs);
+                                faxis = (0:1/nfft:1/2).*fs;
+                                PF = faxis(pos_pf);    % peak frequency
+                                Fmin = faxis(r(1));    % min frequency
+                                Fmax= faxis(r(end));  % max frequency
+                                BW= Fmax - Fmin;    % bandwidth
+                                
+                                %% --- spectral centroid frequency
+                                SCF = sum(spects.*faxis')/sum(spects);
+                                SCF_dB = sum(spectsdB.*faxis')/sum(spectsdB);
+                                
+                                result.call(cn).PF=PF;
+                                result.call(cn).Fmin=Fmin;
+                                result.call(cn).Fmax=Fmax;
+                                result.call(cn).BW=BW;
+                                result.call(cn).SCF=SCF;
+                                result.call(cn).SCF_dB=SCF_dB;
+                                
+                                
+                                % use yin to calculate f0 as a function of time
+                                p.sr = fs;
+                                
+                                R=yin(xcall,p,0);
+                                f0s=R.f0;
+                                nframes=length(f0s);
+                                f0ts=(1:nframes)/(fs/R.hop);
+                                aperiodicity=R.ap0;
+                                
+                                % correct for f0 jumps
+                                fcin=[f0ts' f0s' aperiodicity'];
+                                fdiff=abs(diff(fcin(:,2)));
+                                finds=find(fdiff>9000);
+                                if 1-isempty(finds) % corrections are necessary
+                                    f0_corrected=1;
+                                    f0_corr=fcin(:,2);
+                                    % check start frequency of the irregularity
+                                    if f0_corr(finds(1))>30e3;
+                                        % correct first jump
+                                        f0_corr(1:finds(1))=f0_corr(1:finds(1))/2;
+                                        if length(finds)>1,
+                                            finds=finds(2:end);
+                                        end
+                                    end
+                                    while length(finds)>1
+                                        xfinds=finds(1:2);
+                                        f0_corr(xfinds(1)+1:xfinds(2))=f0_corr(xfinds(1)+1:xfinds(2))/2;
+                                        if length(finds)>2,
+                                            finds=finds(3:end);
+                                        elseif length(finds)==2,
+                                            finds=[];
+                                            break
+                                        else
+                                            break
+                                        end
+                                    end
+                                    if length(finds)>0,
+                                        f0_corr(finds(1)+1:end)=f0_corr(finds(1)+1:end)/2;
+                                    end
+                                    %                                 figure(10), plot(fcin(:,1),fcin(:,2),'b',fcin(:,1),f0_corr,'r');
+                                    %corr_ok=input('correction (1) ok or (0) not ok');
+                                    corr_ok=1;
+                                    if corr_ok==1,
+                                        fcin(:,2)=f0_corr;
+                                    end
+                                else
+                                    f0_corrected=0;
+                                end
+                                
+                                result.call(cn).f0=fcin;
+                                result.call(cn).f0_corrected=f0_corrected;
+                                
+                                % fit f0(t)
+                                f0inds=find(1-isnan(fcin(:,2)));
+                                f0data=fcin(f0inds,1:2);
+                                [f0slope,sfm_amp,sfm_frq,sfm_phase]=fit_f0(f0data,0);
+                                
+                                result.call(cn).f0_slope=f0slope;
+                                result.call(cn).sfm_amp=sfm_amp;
+                                result.call(cn).sfm_frq=sfm_frq;
+                                result.call(cn).sfm_phase=sfm_phase;
+                            end
+                        end
+                        %                         pause;
+                    end
+                    drawnow;
+                end
+                save(savname,'result');
+            else
+                fclose(fid);
+            end
+        end
+    end
+end

vocalisations/scripts_for_reference/fit_f0.m

@@ -0,0 +1,57 @@
+function [f0slope,sfm_amp,sfm_frq,sfm_phase]=fit_f0(f0data, ploton);
+% this function tries to fit f0 as a function of time with a combination of
+% a linear fm plus a sinusoidal fm;
+% input is an x by 2 matrix with the time values in the first column and
+% the corresponding f0 values in the 2nd column
+% the matrix must have at least 5 rows
+% outputs are four values
+% 1st output: linear fm slope (in kHz/ms)
+% 2nd output: amplitude of sfm (in kHz)
+% 3rd output: frequency of SFM component (in Hz)
+% 4th output: phase of SFM component (in rad)
+
+if size(f0data,1)>5,
+    % fit a linear regression
+    [p,s]=polyfit(f0data(:,1),f0data(:,2),1);
+    f0slope=p(1)/1e6;   % in kHz per millisecond
+    linfit=polyval(p,f0data(:,1));
+    if ploton==1,
+        figure(1),
+        subplot(2,1,1);
+        plot(f0data(:,1),f0data(:,2),f0data(:,1),linfit);
+        title(['f0 slope = ' num2str(f0slope) ' kHz/ms']);
+    end
+    
+    % fit a sinusoid to the residual
+    x=f0data(:,1)-min(f0data(:,1)); % reset time to start at 0
+    y=f0data(:,2)-linfit;
+    
+    % estimate for sfm_amp
+    b0(1)=(max(y)-min(y))/2;
+    % estimate for sfm_freq
+    b0(2)= 1/(mean(diff(find(diff(sign(y)))))*2*mean(diff(x)));
+    % estimate for sfm_phase
+    if diff(y(1:2))>0
+        b0(3)=0;
+    else
+        b0(3)=1;
+    end
+    
+    b=nlinfit(x,y,@(b,x) b(1).*sin(2*pi*x*b(2)+b(3)*pi),b0);
+    sinfit=b(1).*sin(2*pi*x*b(2)+b(3)*pi);
+    sfm_amp=b(1)/1e3;   % SFM Amplitude (kHz)
+    sfm_frq=b(2);
+    sfm_phase=b(3);
+    
+    if ploton==1,
+        figure(1),
+        subplot(2,1,2);
+        plot(x,y,'b',  x,sinfit, 'r')
+        title(['SFM FRQ = ' num2str(sfm_frq) ' Hz; SFM AMP = ' num2str(sfm_amp) ' kHz']);
+    end
+else
+    f0slope=nan;
+    sfm_amp=nan;
+    sfm_frq=nan;
+    sfm_phase=nan;
+end

vocalisations/scripts_for_reference/prepare_data.py

@@ -0,0 +1,195 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+import pandas as pd
+import numpy as np
+import scipy.io as sio
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+import os
+import re
+
+FS = 192e3
+DEAF_BATS = {"b3", "b5", "b8"}
+
+def calculate_f0_slope(f0s):
+    unmasked_indices, = np.where(~f0s.mask)
+    if len(unmasked_indices) >= 2:
+        first_index, last_index = unmasked_indices[0], unmasked_indices[-1]
+        first_f0, last_f0 = f0s[first_index], f0s[last_index]
+        return (last_f0 - first_f0) / (last_index - first_index) * 6000  # factor 6000: converts to Hz/seconds
+    else:
+        return np.nan
+
+def process_file(session_id, animals, session_start_time, filename):
+    data = sio.loadmat(filename)
+    if data["result"].size == 0:
+        return pd.DataFrame()
+    data = data["result"][0, 0]
+    start_times = session_start_time + pd.to_timedelta(
+        np.concatenate(data["call"][0, :]["call_start_sample"], axis=1)[0] / FS,
+        unit="s")
+    valid_indices = np.where(list(map(lambda x: x.shape[1] != 0, data["call"][0, :]["call_levels"])))[0]
+    start_samples = np.concatenate(data["call"][0, valid_indices]["call_start_sample"], axis=1)[0]
+    call_durations = np.concatenate(data["call"][0, valid_indices]["call_dur"], axis=1)[0]
+    f_mins = np.concatenate(data["call"][0, valid_indices]["Fmin"], axis=1)[0]
+    f_maxs = np.concatenate(data["call"][0, valid_indices]["Fmax"], axis=1)[0]
+    f0s = [np.ma.masked_invalid(it[:, 1]) for it in data["call"][0, valid_indices]["f0"]]
+    f0s_compressed = [it.compressed() for it in f0s]
+    f0s_start = [it[np.isfinite(it)][0] if len(it) > 0 else np.nan for it in f0s_compressed]
+    f0s_end = [it[np.isfinite(it)][-1] if len(it) > 0 else np.nan for it in f0s_compressed]
+    f0s_slope = [calculate_f0_slope(it) for it in f0s]
+    aperiodicities = [np.ma.masked_invalid(it[:, 2]).compressed() for it in data["call"][0, valid_indices]["f0"]]
+    spectral_centroid = np.concatenate(data["call"][0, valid_indices]["SCF"], axis=1)[0]
+    if animals:
+        call_levels = np.concatenate(data["call"][0, valid_indices]["call_levels"])
+        level_differences = np.mean(call_levels[:, [1, 4, 5]], axis=1) - np.mean(call_levels[:, [0, 2, 3]], axis=1)
+        calling_bat = np.where(level_differences > 0, animals[1], animals[0])
+        other_bat = np.where(level_differences > 0, animals[0], animals[1])
+    else:
+        level_differences = np.nan
+        calling_bat = other_bat = ""
+    call_rms = [10*np.log10(np.mean(it**2)) for it in data["call"][0, valid_indices]["loudest_call"]]
+    return pd.DataFrame(dict(session_id=session_id,
+                             call_id=np.arange(len(start_times)),
+                             start_sample=start_samples,
+                             start_time=start_times,
+                             calling_bat=calling_bat,
+                             other_bat=other_bat,
+                             level_difference=level_differences,
+                             call_rms=call_rms,
+                             call_duration=call_durations,
+                             f_min=f_mins,
+                             f_max=f_maxs,
+                             mean_aperiodicity=[np.mean(it) for it in aperiodicities],
+                             f0_mean=[np.mean(it) for it in f0s],
+                             f0_min=[np.min(it) for it in f0s],
+                             f0_max=[np.max(it) for it in f0s],
+                             f0_start=f0s_start,
+                             f0_end=f0s_end,
+                             f0_slope=f0s_slope,
+                             spectral_centroid=spectral_centroid))
+
+###
+# Pup sessions
+###
+
+PUP_RESULTS_ROOT_DIR = "../raw_data/pups/results"
+
+sessions = []
+calls = []
+
+dn_re = re.compile(r"^vpl_(b\d)(..)$")
+fn_re = re.compile(r"^vpl_...._(\d\d)-(...)-(\d\d\d\d)_(\d\d)x(\d\d)x(\d\d)_m1.mat$")
+
+session_id = 0
+dlist = sorted(os.listdir(PUP_RESULTS_ROOT_DIR))
+for i, dn in enumerate(dlist):
+    dpath = os.path.join(PUP_RESULTS_ROOT_DIR, dn)
+    dn_mt = dn_re.match(dn)
+    if not dn_mt:
+        continue
+    animals = [dn_mt.group(1), dn_mt.group(2)]
+    for fn in sorted(os.listdir(dpath)):
+        fn_mt = fn_re.match(fn)
+        if not fn_mt:
+            continue
+        fpath = os.path.join(dpath, fn)
+        print("\r{}/{}: {}... ".format(i + 1, len(dlist), fpath), end="")
+        start_time = pd.to_datetime("{} {} {} {}:{}:{}".format(*[fn_mt.group(1 + i) for i in range(6)]))
+        df = process_file(session_id, animals, start_time, fpath)
+        calls.append(df)
+        sessions.append((session_id, animals[0], animals[1], start_time))
+        session_id += 1
+
+calls = pd.concat(calls)
+
+sessions = pd.DataFrame(sessions, columns=["session_id", "animal1", "animal2", "start_time"])
+
+sessions.set_index("session_id", inplace=True)
+calls.set_index(["session_id", "call_id"], inplace=True)
+
+calls.insert(calls.columns.get_loc("level_difference"),
+             "calling_bat_deaf",
+             calls["calling_bat"].isin(DEAF_BATS))
+
+sorted_sessions = sessions[sessions["animal2"].str.startswith("m")].sort_values("start_time")
+
+sessions["before_deafening"] = False
+calls.insert(calls.columns.get_loc("level_difference"),
+             "before_deafening",
+             False)
+
+for pup, sessions_per_pup in sorted_sessions.groupby("animal1"):
+    first_id = sessions_per_pup.index[0]
+    sessions.loc[first_id, "before_deafening"] = True
+    calls.loc[calls.index.get_level_values(0) == first_id, "before_deafening"] = True
+
+calls.insert(calls.columns.get_loc("level_difference"),
+             "calling_bat_mother",
+             calls["calling_bat"].str.startswith("m"))
+
+for bool_column in ["before_deafening"]:
+    sessions[bool_column] = sessions[bool_column].astype(np.int)
+for bool_column in ["calling_bat_deaf", "before_deafening", "calling_bat_mother"]:
+    calls[bool_column] = calls[bool_column].astype(np.int)
+
+sessions.to_csv("../pup_sessions.csv")
+calls.to_csv("../pup_calls.csv")
+
+print("\npups done")
+
+###
+# Adults
+###
+
+ADULT_RESULTS_ROOT_DIR = "../raw_data/adults/"
+
+sessions = []
+calls = []
+
+fn_re = re.compile(r"^(\d\d\d\d)_(\d\d)_(\d\d)_(\d\d)_(\d\d)_(\d\d)_call_parameters.mat$")
+
+session_id = 0
+for dn in ["deaf", "hearing"]:
+    dpath = os.path.join(ADULT_RESULTS_ROOT_DIR, dn)
+    for fn in sorted(os.listdir(dpath)):
+        fn_mt = fn_re.match(fn)
+        if not fn_mt:
+            continue
+            
+        print("\r{} {}... ".format(session_id, fn), end="")
+
+        fpath = os.path.join(dpath, fn)
+        groups = list(map(int, fn_mt.groups()))
+        start_time = pd.Timestamp(year=groups[0], month=groups[1], day=groups[2],
+                                  hour=groups[3], minute=groups[4], second=groups[5])
+        df = process_file(session_id, None, start_time, fpath)
+        if len(df) == 0:
+            continue
+
+        df.insert(df.columns.get_loc("level_difference"),
+                  "calling_bat_deaf",
+                  int(dn == "deaf"))
+
+        calls.append(df)
+        sessions.append((session_id, dn, start_time))
+        session_id += 1
+        
+calls = pd.concat(calls)
+calls.set_index(["session_id", "call_id"], inplace=True)
+
+calls.insert(calls.columns.get_loc("level_difference"),
+             "before_deafening", 0)
+calls.insert(calls.columns.get_loc("level_difference"),
+             "calling_bat_mother", 0)
+
+sessions = pd.DataFrame(sessions, columns=["session_id", "group", "start_time"])
+sessions.set_index(["session_id"], inplace=True)
+
+sessions.to_csv("../adult_sessions.csv")
+calls.to_csv("../adult_calls.csv")
+
+print("\nadults done")
+
+# vim:sw=4:sts=4:et: