---
title: 'When Stars Collide - G2Net EDA'
date: '`r Sys.Date()`'
output:
  html_document:
    number_sections: true
    fig_caption: true
    toc: true
    fig_width: 7
    fig_height: 4.5
    theme: cosmo
    highlight: tango
    code_folding: hide
---

```{r setup, include=FALSE, echo=FALSE}
knitr::opts_chunk$set(echo=TRUE, error=FALSE)
knitr::opts_chunk$set(out.width="100%", fig.height = 4.5, split=FALSE, fig.align = 'default')
options(dplyr.summarise.inform = FALSE)
```


<div 
    style="background-image: url('https://svs.gsfc.nasa.gov/vis/a010000/a013000/a013086/image-072-_000150.png'); 
    width:100%; 
    height:400px; 
    background-position:center;">&nbsp;
</div>

*Image: NASA*



# Introduction

Welcome to an initial Exploratory Data Analysis for the [G2Net Gravitational Wave Detection](https://www.kaggle.com/c/g2net-gravitational-wave-detection/overview) competition with [tidy R](http://tidyverse.org/) and [ggplot2](http://ggplot2.tidyverse.org/).

The aim of this challenge is to detect the miniscule signatures of ripples in space time created by super heavy objects in the distant universe. When dense stellar remnants like [Black Holes](https://en.wikipedia.org/wiki/Black_hole) or [Neutron Stars](https://en.wikipedia.org/wiki/Neutron_star) collide with each other, it shakes the very fabric of the universe. The resulting gravitational waves can be measured millions of light years away on earth with highly sensitive instruments. All of this was predicted by Einstein in the early 1900s, but it took humankind 100 years to build detectors that are sensitive enough to discover the first gravitational waves.

The [data](https://www.kaggle.com/c/g2net-gravitational-wave-detection/data) comes in the shape of a large set (70 GB) of time series, together with a label file that describes whether a gravitational wave signal was detected (1) or not (0) in this data. Thus, we are dealing with a **binary classification problem** and are tasked to predict probabilities. The [evaluation metric](https://www.kaggle.com/c/g2net-gravitational-wave-detection/overview/evaluation) is the ROC AUC.

Each time series file contains simultaneous measurements from the 3 detectors LIGO Hanford, LIGO Livingston, and Virgo. Another important aspect to note is that we are dealing with simulated data, since gravitational wave detections are still pretty rare. **The simulated part is the gravitational wave signal, while the detector noise is always real.** These hyper-sensitive instruments are operating at the edge of what's measurable against a noise that comes from simply being located on earth. The [data description](https://www.kaggle.com/c/g2net-gravitational-wave-detection/data) further tells us that the signal simulation is based on astrophysically motivated distributions of 15 parameters, which are hidden from us. Still, it might be expected that the frequencies of the random combinations of those parameters are reflected in the simulated signals.

(Also, yes: it's not technically stars that are colliding, but rather super dense remnants of formerly massive stars. But that doesn't sound quite as poetic now, does it?)

Let's get started!


# Preparations {.tabset .tabset-fade .tabset-pills}

## Load libraries

We load a range of libraries for general data wrangling and general visualisation together with more specialised tools.

```{r, message = FALSE}
# general visualisation
library('ggplot2') # visualisation
library('scales') # visualisation
library('patchwork') # visualisation
library('RColorBrewer') # visualisation
library('corrplot') # visualisation
library('glue') # visualisation
library('ggthemes') # visualisation
library('gt') # tables

# general data manipulation
library('dplyr') # data manipulation
library('readr') # input/output
library('vroom') # input/output
library('skimr') # overview
library('tibble') # data wrangling
library('tidyr') # data wrangling
library('stringr') # string manipulation
library('forcats') # factor manipulation
library('janitor') # data cleaning

# specific visualisation
library('alluvial') # visualisation
library('ggrepel') # visualisation
library('ggforce') # visualisation
library('ggridges') # visualisation
library('gganimate') # animations
library('GGally') # visualisation
library('wesanderson') # visualisation
library('corrplot') # visualisation

# specific data manipulation
library('lazyeval') # data wrangling
library('broom') # data wrangling
library('purrr') # data wrangling
library('reshape2') # data wrangling
library('rlang') # encoding

# Date plus forecast
library('lubridate') # date and time
library('lomb') # power spectrum
library('tuneR') # spectra
library('seewave') #spectra

# Modelling
#library('xgboost')
#library('rsample')
#library('recipes')
#library('parsnip')
#library('yardstick')
```


## Reticulate support for Python files

The time series are being provided as python numpy files. Nice try, but that doesn't stop us at all from using R ;-)

The file conversion is powered by [reticulate](https://rstudio.github.io/reticulate/).

```{r}
library(reticulate)
np <- import("numpy")
```


## Load data

We use the [vroom](https://cran.r-project.org/web/packages/vroom/vignettes/vroom.html) package for reading data.

```{r echo=FALSE}
if (dir.exists("/kaggle")){
  path <- "/kaggle/input/g2net-gravitational-wave-detection/"
} else {
  path <- ""
}
```

```{r warning=FALSE, results=FALSE}
train <- vroom(str_c(path, 'training_labels.csv'), col_types = cols())
sample_submit <- vroom(str_c(path, 'sample_submission.csv'), col_types = cols())
```


# Overview: File structure and content {.tabset .tabset-fade .tabset-pills}

To start, we'll get a quick overview of the datasets by looking at the raw content of the files and at their shapes.


## Training data set


The train labels are provided in `training_labels.csv`. Here is what the first rows of that file look like:

```{r}
head(train, 50) %>% 
  gt() %>% 
  tab_header(
    title = md("*G2Net Training data labels*")
    ) %>% 
  opt_row_striping() %>% 
  fmt_markdown(columns = everything()) %>% 
  tab_style(
    style = list(
      cell_fill("lightblue"),
      cell_text(color = "white", weight = "bold")
      ),
    locations = cells_title(groups = "title")
  ) %>% 
  tab_source_note("First 50 rows.") %>% 
  tab_options(container.height = px(350))
```


```{r}
glue("Number of rows: { nrow(train) }; Number of columns { ncol(train) }")
```

We find:

- This file only contains the `target` labels, together with the file `id` of the time series data.

- The targets look pretty balanced; let's have a closer look further below.

- There are 560k time series in the training set, which is a sizeable amount.



## Test data - sample submission

The test data IDs are provided in the `sample_submission.csv` file, since there are no additional features other than what's in the time series data. This sample submission shows us what our prediction should look like: a probability for a signal being present for each of the time series IDs:


```{r}
sample_submit %>% 
  head(50) %>% 
  gt() %>% 
  tab_header(
    title = md("*G2Net Test data IDs*")
    ) %>% 
  opt_row_striping() %>% 
  fmt_markdown(columns = everything()) %>% 
  tab_style(
    style = list(
      cell_fill("lightgreen"),
      cell_text(color = "white", weight = "bold")
      ),
    locations = cells_title(groups = "title")
  ) %>% 
  tab_source_note("First 50 rows.") %>% 
  tab_options(container.height = px(350))
```



```{r}
glue("Number of sample submit rows: { nrow(sample_submit) }; Number of sample submit columns { ncol(sample_submit) }")
```

We find:

- With 226k rows, the test data is about 40% as large as the training data. Still a pretty large dataset, though.

- The [leaderboard](https://www.kaggle.com/c/g2net-gravitational-wave-detection/leaderboard) tells us that the public LB is calculated with about 16% of the data, which is about 36k files. The remaining 84% are reserved for the private LB.



## Training folders and file paths


The individual time series are organised in a nested folder structure with 3 additional levels inside the parent `train` and `test` folders. The names of these 3 folders come from the first 3 digits of the file `id`. For instance, `id = 00000e74ad` can be found in folder `train/0/0/0/`. Or folder `train/b/8/2/` contains the file `b820212efb.npy`. You get the idea.

To make things easier, we can add the path and the filename to our train and test label files.

```{r}
train <- train %>% 
  mutate(path = str_c(path, "train/", str_sub(id, 1, 1), "/", str_sub(id, 2, 2), "/", str_sub(id, 3, 3), "/"),
         filename = str_c(path, id, ".npy"))

test <- sample_submit %>% 
  select(id) %>% 
  mutate(path = str_c(path, "test/", str_sub(id, 1, 1), "/", str_sub(id, 2, 2), "/", str_sub(id, 3, 3), "/"),
         filename = str_c(path, id, ".npy"))
```


Here's what the amended table looks like for the training data:

```{r}
train %>% 
  head(50) %>% 
  gt() %>% 
  tab_header(
    title = md("*G2Net Training data labels - augmented*")
    ) %>% 
  opt_row_striping() %>% 
  fmt_markdown(columns = everything()) %>% 
  tab_style(
    style = list(
      cell_fill("lightblue"),
      cell_text(color = "white", weight = "bold")
      ),
    locations = cells_title(groups = "title")
  ) %>% 
  tab_source_note("First 50 rows.") %>% 
  tab_options(container.height = px(350))
```




## Missing values

Just in case, since it's always better to check: there are no missing values in our training labels file:

```{r}
sum(is.na(train))
```


## Example Time Series

Let's also load an example time series, to get an impression of what we will be training on. The time series data is stored in numpy format, but with the fantastic `reticulate` package at our fingertips we have no problem in reading this data into R.


```{r warning=FALSE}
ts <- np$load(train$filename[1]) %>% 
  as.matrix() %>%
  t() %>% 
  as_tibble() %>% 
  rownames_to_column(var = "time") %>% 
  mutate(time = as.numeric(time)) %>% 
  rename(site_1 = V1, site_2 = V2, site_3 = V3)
```


```{r}
ts %>% 
  head(100) %>% 
  gt() %>% 
  tab_header(
    title = md("*G2Net Time Series - ID 00000e74ad*")
    ) %>% 
  opt_row_striping() %>% 
  fmt_markdown(columns = everything()) %>% 
  tab_style(
    style = list(
      cell_fill("orange1"),
      cell_text(color = "white", weight = "bold")
      ),
    locations = cells_title(groups = "title")
  ) %>% 
  tab_source_note("First 100 rows.") %>% 
  tab_options(container.height = px(350))
```


We find:

- Each file contains 3 time series from 3 different detectors at different sites: LIGO Hanford, LIGO Livingston, and Virgo. I'm choosing to label those columns `site_1` - `site_3`.

- The values themselves are notably small, which reflects the difficulties in detecting graviational waves in the first place.


We're turning our code to read in the time series into a short convenience function:

```{r}
read_ts <- function(filename){
  
  ts <- np$load(file = filename) %>% 
    as.matrix() %>%
    t() %>% 
    as_tibble() %>% 
    rownames_to_column(var = "time") %>% 
    mutate(time = as.numeric(time)) %>% 
    rename(site_1 = V1, site_2 = V2, site_3 = V3)
  
  return(ts)
  
}
```



# Individual feature visualisations

Looking at the raw data itself is always a great first step to see immediately what kind of features we're dealing with. This also prepares us for the following visual analysis, where we aim to gain a broad understanding of our data shapes and distributions.


## Target labels


```{r fig.cap ="Fig. 1"}
train %>% 
  count(target) %>% 
  mutate(target = as.factor(target)) %>% 
  mutate(perc = n / nrow(train)) %>% 
  ggplot(aes(target, perc, fill = target)) +
  geom_col() +
  scale_y_continuous(labels = scales::percent) +
  theme_hc() +
  theme(legend.position = "none") + 
  labs(x = "Target", y = "", title = "Target labels are well balanced")
```

We find:

- Now this is the advantage of working with simulated data: the target labels are perfectly balanced at 50/50.

- However, this might not translate well to any practical applications, where the probability of actually detecting gravitational waves is rather small.



## Example Time Series - Target vs No Target {.tabset .tabset-fade .tabset-pills}


Let's also look at one of the time series files.

We will plot the time 3 individual time series together with their distributions for `id = 00000e74ad.npy`, the first ID in the training data, which is supposed to contain a signal (`target = 1`). We also convert the time axis to seconds.


```{r}
ts <- read_ts(train$filename[1])
```


```{r fig.cap ="Fig. 2"}
p1 <- ts %>% 
  pivot_longer(starts_with("site"), names_to = "site", values_to = "val") %>% 
  mutate(site = str_to_sentence(str_replace(site, "_", " "))) %>% 
  mutate(time = time/max(time) * 2) %>% 
  ggplot(aes(time, val, col = site)) +
  geom_line() +
  scale_y_continuous(labels=function(x) sprintf("%.1e", x)) +
  facet_wrap(~site, ncol = 1) +
  theme_hc() +
  theme(legend.position = "none") + 
  labs(x = "Time [s]", y = "values")

p2 <- ts %>% 
  pivot_longer(starts_with("site"), names_to = "site", values_to = "val") %>% 
  mutate(site = str_to_sentence(str_replace(site, "_", " "))) %>% 
  ggplot(aes(val, fill = site)) +
  geom_density() +
  coord_flip() +
  facet_wrap(~site, ncol = 1) +
  theme_hc() +
  theme(legend.position = "none", axis.text = element_blank()) + 
  labs(x = "", y = "Density")

design <- "
1112
"

p1 + p2 + plot_layout(design = design) + plot_annotation(title = "Example Time Series - 00000e74ad (target = 1)")
```

We find:

- Measured on the same y-axis scale, the difference in variance between site 3 vs 1 and 2 is very notable. Site 3 might be Virgo, whereas 1 & 2 could be the 2 LIGO detectors.

- The sponsors weren't kidding when they wrote in the data description that "in nearly all cases [...] these signals are not visible by eye in the time series".


Let's turn this plotting layout into another function and use it to compare the above plot to the second ID in the train data, "00001f4945", which has `target = 0`:


```{r}
plot_ts <- function(file_name){
  
  foo <- train %>% 
    filter(filename == file_name)
  
  file_id <- foo$id[1]
  target_val <- foo$target[1] 

  ts <- read_ts(file_name)
  
  p1 <- ts %>% 
    pivot_longer(starts_with("site"), names_to = "site", values_to = "val") %>% 
    mutate(site = str_to_sentence(str_replace(site, "_", " "))) %>% 
    mutate(time = time/max(time) * 2) %>% 
    ggplot(aes(time, val, col = site)) +
    geom_line() +
    scale_y_continuous(labels=function(x) sprintf("%.1e", x)) +
    facet_wrap(~site, ncol = 1) +
    theme_hc() +
    theme(legend.position = "none") + 
    labs(x = "Time [s]", y = "values")
  
  p2 <- ts %>% 
    pivot_longer(starts_with("site"), names_to = "site", values_to = "val") %>% 
    mutate(site = str_to_sentence(str_replace(site, "_", " "))) %>% 
    ggplot(aes(val, fill = site)) +
    geom_density() +
    coord_flip() +
    facet_wrap(~site, ncol = 1) +
    theme_hc() +
    theme(legend.position = "none", axis.text = element_blank()) + 
    labs(x = "", y = "Density")
  
  design <- "
  1112
  "
  
  print(p1 + p2 + plot_layout(design = design) + plot_annotation(title = glue("Example Time Series - {file_id} (target = {target_val})")))

}
```


```{r fig.cap ="Fig. 3"}
plot_ts(train$filename[2])
```

We find:

- Well, this might be tough. Those two IDs look pretty similar, which means that most of that variance we see are noise fluctuations.

- Still, no reason to despair. We have quite a few tricks left up our sleeves.


We also don't want to draw too many conclusions from looking at only two files, so I've prepared a small, random selection of 6 observations with `target == 1` vs 6 files with `target == 0`. I have arranged those time series in individual tabs so that you can click and compare between them. Here is a table of the 12 files:


```{r}
set.seed(4321)
samp <- train %>% 
  filter(target == 1) %>% 
  sample_n(6) %>% 
  bind_rows(train %>% filter(target == 0) %>% sample_n(6))

samp %>% 
  gt() %>% 
  tab_header(
    title = md("*G2Net Sample Rows*")
    ) %>% 
  opt_row_striping() %>% 
  fmt_markdown(columns = everything()) %>% 
  tab_style(
    style = list(
      cell_fill("purple1"),
      cell_text(color = "white", weight = "bold")
      ),
    locations = cells_title(groups = "title")
  ) %>% 
  tab_options(container.height = px(350))
```


And here are the tabs. We see that the overall shapes of the time series' and their distributions are indeed virtually indistinguishable to the eye for the scenarios of target 1 vs target 0. 


### ID 7addd6b05e - target 1

```{r fig.cap ="Fig. 4a"}
plot_ts(samp$filename[1])
```


### ID 7151c27f11 - target 1

```{r fig.cap ="Fig. 4b"}
plot_ts(samp$filename[2])
```


### ID f5fc77cff1 - target 1

```{r fig.cap ="Fig. 4c"}
plot_ts(samp$filename[3])
```


### ID ca5745412b - target 1

```{r fig.cap ="Fig. 4d"}
plot_ts(samp$filename[4])
```


### ID f38ac3cb63 - target 1

```{r fig.cap ="Fig. 4e"}
plot_ts(samp$filename[5])
```



### ID d234601c58 - target 1

```{r fig.cap ="Fig. 4f"}
plot_ts(samp$filename[6])
```



### ID 467ff7639e - target 0


```{r fig.cap ="Fig. 4g"}
plot_ts(samp$filename[7])
```



### ID fad2fe47ce - target 0

```{r fig.cap ="Fig. 4h"}
plot_ts(samp$filename[8])
```



### ID b8e3092db7 - target 0

```{r fig.cap ="Fig. 4i"}
plot_ts(samp$filename[9])
```


### ID b2363a8b95 - target 0

```{r fig.cap ="Fig. 4j"}
plot_ts(samp$filename[10])
```


### ID 156cc79304 - target 0

```{r fig.cap ="Fig. 4k"}
plot_ts(samp$filename[11])
```


### ID 71b9b975a5 - target 0

```{r fig.cap ="Fig. 4l"}
plot_ts(samp$filename[12])
```



## Periodograms {.tabset .tabset-fade .tabset-pills}

Whenever we're dealing with time series in some shape or form, it's always a good idea to check some periodograms. Those will tells us about the underlying patterns in the data, and about periodicities on different time scales.

We'll start with simple 1d periodograms. Here we'll be using [Lomb-Scargle](https://en.wikipedia.org/wiki/Least-squares_spectral_analysis), i.e. fitting sinusoids via the `lomb` package, instead of Fourier spectra. Partly to get a slightly different take, partly because this is an astronomy competition and Lomb-Scargle is popular in astronomy.


```{r}
ts <- read_ts(train$filename[1])
```


```{r}
lomb_1 <- lsp(ts$site_1, ts$time, type = 'period', alpha = 0.05,
              main = "Lomb-Scargle Periodogram with 95% confidence level",
              xlab = "Period [time]", ylab = "Normalised Power",
              cex.lab = 1.5, cex.axis = 1.5, plot = FALSE)

lomb_2 <- lsp(ts$site_2, ts$time, type = 'period', alpha = 0.05,
              main = "Lomb-Scargle Periodogram with 95% confidence level",
              xlab = "Period [time]", ylab = "Normalised Power",
              cex.lab = 1.5, cex.axis = 1.5, plot = FALSE)

lomb_3 <- lsp(ts$site_3, ts$time, type = 'period', alpha = 0.05,
              main = "Lomb-Scargle Periodogram with 95% confidence level",
              xlab = "Period [time]", ylab = "Normalised Power",
              cex.lab = 1.5, cex.axis = 1.5, plot = FALSE)
```



Here is the time series for `id = 00000e74ad.npy` and site 1, along with the corresponding power spectrum.

In terms of styling, I've removed the values from the time series y-axis, since this is more about the relative variations then the absolute amplitudes. This gives more space to the time series graph itself. The `classic` plot theme does a better job at aligning the left and right panels to the eye than my otherwise preferred `theme_hc` or `theme_minimal`.


```{r fig.cap = "Fig. 5", fig.height=3.5}
lomb_tidy <- tibble(
  period = lomb_1$scanned,
  power = lomb_1$power
)

lomb_max <- lomb_tidy %>% 
  arrange(desc(power)) %>% 
  dplyr::slice(1)

p1 <- ts %>% 
    ggplot(aes(time, site_1)) +
    geom_line(col = "red2") +
    theme_classic() +
    theme(legend.position = "none", title = element_text(size = 8), axis.text.y = element_blank()) + 
    labs(y = "", title = glue("ID {train$id[1]}\nSite 1"))

p2 <- lomb_tidy %>% 
  ggplot(aes(period, power)) +
  geom_line() +
  geom_hline(yintercept = lomb_1$sig.level, linetype = 2) +
  scale_x_log10() +
  theme_classic() +
  theme(title = element_text(size = 8)) +
  labs(title = str_c("Lomb-Scargle Periodogram with 95% SL\nMax Power ",
                     sprintf("%.0f", lomb_max$power),
                     " at period ",
                     sprintf("%.1f", lomb_max$period),
                     " units"))

p1 + p2
```

We find:

- A very strong signal at about 186 time units, which at 2048 Hz corresponds to about 0.09 s or about 11 Hz. (Remember that each time series only covers a total of 2 seconds). This periodicity is 

- The bad news: This is not the signal that we're looking for. The good news: This feature might help in understanding the noise characteristics.



The next step is to turn the spectrogram plot into a function, so that we can compare `target = 0` vs `target = 1` side by side. This will display the 3 sites separately. I'm also including parameters to narrow the x-axis range, to allow us to better focus on the interesting region of the spectrum.


```{r}
plot_lomb <- function(file_name, xmin = 10, xmax = 2000, ymax = 0){
  
  if (ymax == 0){
    ylim <- NULL
  } else {
    ylim <- c(0, ymax)
  }
  
  foo <- train %>% 
    filter(filename == file_name)
  
  file_id <- foo$id[1]
  target_val <- foo$target[1] 
  
  ts <- read_ts(file_name)

  lomb_1 <- lsp(ts$site_1, ts$time, type = 'period', alpha = 0.05,
                main = "Lomb-Scargle Periodogram with 95% confidence level",
                xlab = "Period [time]", ylab = "Normalised Power",
                cex.lab = 1.5, cex.axis = 1.5, plot = FALSE)
  
  lomb_2 <- lsp(ts$site_2, ts$time, type = 'period', alpha = 0.05,
                main = "Lomb-Scargle Periodogram with 95% confidence level",
                xlab = "Period [time]", ylab = "Normalised Power",
                cex.lab = 1.5, cex.axis = 1.5, plot = FALSE)
  
  lomb_3 <- lsp(ts$site_3, ts$time, type = 'period', alpha = 0.05,
                main = "Lomb-Scargle Periodogram with 95% confidence level",
                xlab = "Period [time]", ylab = "Normalised Power",
                cex.lab = 1.5, cex.axis = 1.5, plot = FALSE)
  
  lomb_tidy_1 <- tibble(
    period = lomb_1$scanned,
    power = lomb_1$power
  ) %>% 
    mutate(site = "site_1")
  
  lomb_tidy_2 <- tibble(
    period = lomb_2$scanned,
    power = lomb_2$power
  ) %>% 
    mutate(site = "site_2")
  
  lomb_tidy_3 <- tibble(
    period = lomb_3$scanned,
    power = lomb_3$power
  ) %>% 
    mutate(site = "site_3")
  
  lomb_tidy <- lomb_tidy_1 %>% 
    bind_rows(lomb_tidy_2) %>% 
    bind_rows(lomb_tidy_3)
  
  p <- lomb_tidy %>% 
    ggplot(aes(period, power, col = site)) +
    geom_line() +
    scale_x_log10() +
    coord_cartesian(xlim = c(xmin, xmax), ylim = ylim) +
    facet_wrap(~site, ncol = 1) +
    theme_minimal() +
    theme(legend.position = "none", title = element_text(size = 8)) +
    labs(title = glue("LS Periodogram - {file_id} (target = {target_val})"))
  
  return(p)
  
  
}
```


```{r fig.cap = "Fig. 6"}
p1 <- plot_lomb(train$filename[1], 90, 500)
p2 <- plot_lomb(train$filename[2], 90, 500)

p1 + p2
```

We find:

- See what I mean? Most of the power in this periodogram appears to come from the noise. The shapes and peaks are very similar between `target = 1` vs `target = 0`. Or at least that's true for those 2 IDs. There's nothing that immediately screams "signal", though.

- We could now sample the interesting region of the power spectrum with higher frequency resolution. Maybe we'll run some parametrisations later. For now, let's look at some more IDs for context.


To get a somewhat broader view, we will look at the power spectra of our same 12 sample IDs from above. We could choose a different sample here, but I think there is value in being able to compare a given spectrum to its time series. I'm again arranging the different views in individual tabs to make them easier to navigate.


### ID 7addd6b05e - target 1

```{r fig.cap ="Fig. 7a"}
plot_lomb(samp$filename[1], 90, 500)
```


### ID 7151c27f11 - target 1

```{r fig.cap ="Fig. 7b"}
plot_lomb(samp$filename[2], 90, 500)
```


### ID f5fc77cff1 - target 1

```{r fig.cap ="Fig. 7c"}
plot_lomb(samp$filename[3], 90, 500)
```


### ID ca5745412b - target 1

```{r fig.cap ="Fig. 7d"}
plot_lomb(samp$filename[4], 90, 500)
```


### ID f38ac3cb63 - target 1

```{r fig.cap ="Fig. 7e"}
plot_lomb(samp$filename[5], 90, 500)
```



### ID d234601c58 - target 1

```{r fig.cap ="Fig. 7f"}
plot_lomb(samp$filename[6], 90, 500)
```



### ID 467ff7639e - target 0


```{r fig.cap ="Fig. 7g"}
plot_lomb(samp$filename[7], 90, 500)
```



### ID fad2fe47ce - target 0

```{r fig.cap ="Fig. 7h"}
plot_lomb(samp$filename[8], 90, 500)
```



### ID b8e3092db7 - target 0

```{r fig.cap ="Fig. 7i"}
plot_lomb(samp$filename[9], 90, 500)
```


### ID b2363a8b95 - target 0

```{r fig.cap ="Fig. 7j"}
plot_lomb(samp$filename[10], 90, 500)
```


### ID 156cc79304 - target 0

```{r fig.cap ="Fig. 7k"}
plot_lomb(samp$filename[11], 90, 500)
```


### ID 71b9b975a5 - target 0

```{r fig.cap ="Fig. 7l"}
plot_lomb(samp$filename[12], 90, 500)
```




## Two-dimensional spectrograms

This is likely the way to go in this competition: turn the data into 2d spectrogram images and then use CNN based models. This is essentially what we're doing with audio data and melspectrograms. And if you think about it, our task to find weak signals in strong environmental noise is not that different from detecting bird songs in recordings from microphones in nature. I don't have much experience in this area, but I will give it a try and hopefully learn a thing or two.

Here I'm using the `tuneR` package (which appears to be the R equivalent of the popular Python library `librosa`) together with a tool called `seewave` which allows us to visualise spectrograms quickly. This workflow requires converting the time series into a `Wave` object, from which we then extract and plot the spectrogram.


```{r}
ts <- read_ts(train$filename[1])
bar <- Wave(left = ts$site_1, samp.rate = 2048, bit = 16)
```


Here is the spectrogram for the frequency range between 0 and 25 Hz. The plot is produced by `seewave`, which also defines most of the styling parameters. In the small bottom panel you can see the original time series below the main spectrogram. We're looking again at the first `id = 00000e74ad.npy`, site 1:

```{r fig.cap = "Fig. 8"}
spectro(bar, osc = TRUE, flim = c(0, 0.025), wl = 512, collevels = seq(-30, 0, .5), palette = viridis::inferno, ovlp = 50,
        main = glue("Zoomed Spectrogram - {train$id[1]} - site 1"))
```

We find:

- This spectrogram confirms the view we got in Fig. 5 above: There is a strong (noise) signal around 11 Hz with some secondary peaks around it. We can see the areas in the time series, before 0.5s and 1.5s, where this signal is weaker.

- Here we cut off the amplitude scale at -30 dB (see the colour axis on the right), because there just isn't that much signal in the rest of the spectrogram. This is also the reason for the narrow frequency range. This narrow range in turn is responsible for the rather low frequency resolution, since the tool is apparently first sampling over the entire range and then zooming in afterwards.


If we look at the full spectrogram, which we can compute up to about 1 kHz, then this is what we get (with amplitudes down to -100 dB).

```{r fig.cap = "Fig. 9"}
spectro(bar, osc = TRUE, flim = c(0, 1), wl = 512, collevels = seq(-100, 0, .5), palette = viridis::inferno, ovlp = 50, flog = TRUE,
        main = glue("Full Spectrogram - {train$id[1]}- Site 1"))
```

We find:

- There really isn't that much else in the spectrum. You can see the interesting region squashed into the bottom of the spectrogram near zero.

- Everything above 25 Hz looks pretty similar. Maybe the real signal is still hiding in there. The images shown in the [competition description](https://www.kaggle.com/c/g2net-gravitational-wave-detection/) show power spectra with a strong signal around 100 - 200 Hz. It's possible that we want to search in a similar region.


---

To be continued ...