monogenea

# Check first five columns in matrix.csv
#!cut -d, -f-5 matrix.csv | head

# Import data with Bash command discarding first column
matrix = dt.fread(cmd='cut -d, -f2- matrix.csv',
                  header=True, sep=',', columns=dt.int32) # ~7 GB (76533, 50281)
# Import metadata
metadata = pd.read_csv('metadata.csv')

monogenea

# Define UMAP
brain_umap = umap.UMAP(random_state=999, n_neighbors=30, min_dist=.25)

# Fit UMAP and extract latent vars 1-2
embedding = pd.DataFrame(brain_umap.fit_transform(matrix), columns = ['UMAP1','UMAP2'])

# Produce sns.scatterplot and pass metadata.subclasses as color
sns_plot = sns.scatterplot(x='UMAP1', y='UMAP2', data=embedding,
                hue=metadata.subclass_label.to_list(),
                alpha=.1, linewidth=0, s=1)

monogenea

# Imports
#!pip install datatable
import os, umap
import numpy as np
import pandas as pd
import datatable as dt
import seaborn as sns

os.chdir('PATH/TO/WDIR')

monogenea

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("pcaMethods")
library(pcaMethods)

monogenea

# read, downsample, clip, mel spec, normalize and remove noise
melspec <- function(x, start, end){
        mp3 <- readMP3(filename = x) %>% 
                extractWave(xunit = "time",
                            from = start, to = end)
        
        # return log-spectrogram with 256 Mel bands and compression
        sp <- melfcc(mp3, nbands = 256, usecmp = T,
                     spec_out = T,
                     hoptime = (end-start) / 256)$aspectrum

monogenea

# Encode species from fnames regex
species <- str_extract(fnames, patt = "[A-Za-z]+-[a-z]+") %>%
      gsub(patt = "-", rep = " ") %>% factor()

# Stratified sampling: train (80%), val (10%) and test (10%)
set.seed(100)
idx <- createFolds(species, k = 10)
valIdx <- idx$Fold01
testIdx <- idx$Fold02
# Define samples for train, val and test

monogenea

# Define targets and augment data
target <- model.matrix(~0+species)

targetTrain <- do.call("rbind", lapply(1:(dim(Xtrain)[1]/length(fnamesTrain)),
                                       function(x) target[-c(valIdx, testIdx),]))
targetVal <- do.call("rbind", lapply(1:(dim(Xval)[1]/length(fnamesVal)),
                                     function(x) target[valIdx,]))
targetTest <- do.call("rbind", lapply(1:(dim(Xtest)[1]/length(fnamesTest)),
                                      function(x) target[testIdx,]))
# Assemble Xs and Ys

monogenea

# Plot spectrogram from random training sample - range is 0-22.05 kHz
image(train$X[sample(dim(train$X)[1], 1),,,],
      xlab = "Time (s)",
      ylab = "Frequency (kHz)",
      axes = F)
# Generate mel sequence from Hz points, standardize to plot
freqs <- c(0, 1, 5, 15, 22.05)
mels <- 2595 * log10(1 + (freqs*1e3) / 700) # https://en.wikipedia.org/wiki/Mel_scale
mels <- mels - min(mels)
mels <- mels / max(mels)

monogenea

# Fri Feb  7 15:49:46 2020 ------------------------------
setwd("~/Documents/Tutorials/birdsong")
library(keras)
use_condaenv("plaidml")
use_backend("plaidml")
k_backend() # plaidml
library(tidyverse)
library(caret)
library(e1071)
library(pheatmap)

monogenea

# Build model
model <- keras_model_sequential() %>% 
        layer_conv_2d(input_shape = dim(train$X)[2:4], 
                      filters = 16, kernel_size = c(3, 3),
                      activation = "relu") %>% 
        layer_max_pooling_2d(pool_size = c(2, 2)) %>% 
        layer_dropout(rate = .2) %>% 
        
        layer_conv_2d(filters = 32, kernel_size = c(3, 3),
                      activation = "relu") %>%