Music Dynamics Functions

DESCRIPTION

@@ -18,7 +18,7 @@ Imports:
     ggplot2,
     stats
 LinkingTo: Rcpp, RcppArmadillo
-RoxygenNote: 6.1.1
+RoxygenNote: 7.0.2
 Suggests: 
     knitr,
     rmarkdown,

NAMESPACE

@@ -1,13 +1,19 @@
 # Generated by roxygen2: do not edit by hand
 
+export()
 export(beamSearch)
 export(convert10to4)
+export(convertdynamicstates)
 export(ddirichlet)
 export(getloglike)
 export(kalman)
 export(musicModel)
+export(musicModeldynamics)
 export(plotStates)
+export(plotStates2)
 export(rdirichlet)
+export(resampleOptimal)
+export(resampleSubOptimal)
 import(ggplot2)
 importFrom(Rcpp,evalCpp)
 importFrom(Rcpp,sourceCpp)

R/musicfuns.R

@@ -18,6 +18,27 @@ convert10to4 <- function(path){
   path10
 }
 
+#' Converts from the 0,1,2,3 to 1,2,3,4 notation
+#' 
+#'
+#' @param path path in four-state form
+#'
+#' @return path in two-state form
+#' 
+#' @examples
+#' path = sample(1:10, 25, replace=TRUE)
+#' npath = convert10to4(path)
+#' 
+#' @export
+convertdynamicstates <- function(path){
+  t1 = c(1,2,3,4)
+  path10 = t1[path+1]
+  path10
+}
+
+
+
+
 
 # Deprecated --------------------------------------------------------------
 

R/plots.R

@@ -52,3 +52,128 @@ plotStates <- function(performance, params, y, onset,
                      legend.title = ggplot2::element_blank()) +
       ggplot2::ggtitle(performance)
 }
+
+
+#' Plots the Music Model states and Kalman smoothed estimates
+#'
+#' @param theta parameters in same order
+#' 
+#' @param states vector of states at each time period (not converted)
+#' 
+#' @param y vector of tempos or dynamics of a single performance
+#' 
+#' @param onset vector of onsets.
+#' 
+#' @param model If using the dynamics model put "dynamics".  If using the tempo model put "tempo".
+#' 
+#' @param priormean Prior Mean array
+#' 
+#' @param priorvar Prior Variance array
+#' 
+#' @param title Change the title of the plot
+#'
+#' @return A plot showing the estimated values of the kalman smoother along with the states.
+#' 
+#' @examples
+#'theta = c(426.69980736, 136.33213703, -11.84256691, -34.82234559, 
+#'          439.37886221, 1, 1, 0.84916635, 0.04611644, 0.74119571, 
+#'          0.43966082, 0.02116317, 0.24513563, 0.17253254)
+#'y = matrix(tempos[,'Richter_1976'], 1)
+#'lt = diff(c(tempos$note_onset, 61))
+#'pmats = musicModel(lt, theta[1], theta[2:4], theta[5:7], theta[8:14], 
+#'                   c(132,0), c(400,10)) # prior means and variances on X_1
+#'beam = with(pmats, beamSearch(a0, P0, c(1,0,0,0,0,0,0,0,0,0), dt, ct, Tt, Zt,
+#'                              HHt, GGt, y, transMat, 200))
+#'states = beam$paths[which.max(beam$weights),]
+#'
+#'plotStates2(theta,states,y,tempos$note_onset,model="tempo",c(132,0),c(400,10),title="Richter 1976")
+#'
+#'
+#'theta = c(20,10,
+#'          12,13,14,
+#'          21,22,23,
+#'          .3,.2,.1,.5)
+#'y = matrix(dynamics[,'Richter_1976'], 1)
+#'lt = diff(c(dynamics$note_onset, 61))
+#'
+#'pmats = musicModeldynamics(lt, theta[1], theta[2], theta[3:5], theta[6:8], theta[9:12], c(20,0,0), c(25,10,0)) 
+#'beam = with(pmats, beamSearch(a0, P0, c(1,0,0,0), dt, ct, Tt, Zt,
+#'                              HHt, GGt, y, transMat, 200))
+#'states = beam$paths[which.max(beam$weights),]
+#'kal <- kalman(pmats,states,y)
+#'plotStates2(theta,states,y,dynamics$note_onset,model="dynamics",c(20,0,0), c(25,10,0),title="Richter 1976")
+#' 
+#' @export
+plotStates2 <- function(theta, states, y, onset, model, priormean, priorvar, title="Music Plot"){
+  lt = diff(c(onset,61))
+
+  if(model=="dynamics"){
+    mats = musicModeldynamics(lt, theta[1], theta[2], theta[3:5], theta[6:8], theta[9:12], theta[13],
+                              priormean, priorvar)
+    kal <- kalman(mats,states,y)
+    df = data.frame(
+      measure = onset, 
+      dynamics = c(y), 
+      inferred = c(kal$ests), 
+      state = factor(
+        states,
+        levels=c(0,1,2,3),
+        labels=c('New Value', 'Smooth Progression','Loud Deviation','Soft Deviation')
+      )
+    )
+    df2 <- data.frame(
+      start = c(0,9,17,21,33,45,53),
+      end = c(8,16,20,32,44,52,60)+1,
+      direction = factor(c(1,0,2,0,4,1,0),
+                         levels=c(0,1,2,4),
+                         labels=c("p","f","ff","poco piu vivo"))
+    )
+    myplot <- ggplot2::ggplot(df) + 
+      geom_rect(data=df2, aes(NULL,NULL,xmin=start,xmax=end,fill=direction),
+                ymin=-Inf,ymax=Inf,color="white",size=0.5,alpha=0.3)+
+      scale_fill_manual("Composer Direction",values=c("p"="grey95",
+                                                      "f"="grey90",
+                                                      "ff"="grey75",
+                                                      "sf"="grey50",
+                                                      "poco piu vivo"="lightblue"))+
+      ggplot2::geom_line(ggplot2::aes(x=measure, y=dynamics), color='grey10') +
+      ggplot2::geom_point(ggplot2::aes(x=measure, y=inferred, color=state)) +
+      scale_color_manual("Performer Intention",values=c('New Value'="blue",
+                                                       'Smooth Progression'="darkcyan",
+                                                       'Loud Deviation'="darkorange1",
+                                                       'Soft Deviation'="red"))+
+      ggplot2::theme(legend.position = 'right') +
+      ggplot2::ggtitle(title) +
+      cowplot::theme_cowplot() +
+      ggplot2::labs(x="Measure", y="Dynamics")
+  }
+
+  if(model=="tempo"){
+    mats = musicModel(lt, theta[1], theta[2:4], theta[5:7], theta[8:14], 
+                      priormean, priorvar)
+    kal <- kalman(pmats,states,y)
+    df = data.frame(
+      measure = onset, 
+      tempo = c(y), 
+      inferred = c(kal$ests), 
+      state = factor(
+        convert10to4(states),
+        levels=c(1,2,3,4),
+        labels=c('constant tempo', 'decelerating','accelerating','stress')
+      )
+    )
+    myplot <- ggplot2::ggplot(df) + 
+      ggplot2::geom_rect(
+        data=data.frame(xmin = 33, xmax = 45, ymin = -Inf, ymax = Inf),
+        mapping=ggplot2::aes(xmin=xmin,xmax=xmax,ymin=ymin,ymax=ymax),
+        fill = 'gray90', color = 'gray90') +
+      ggplot2::geom_line(ggplot2::aes(x=measure, y=tempo), color='black') +
+      ggplot2::geom_point(ggplot2::aes(x=measure, y=inferred, color=state)) +
+      ggplot2::scale_color_brewer(palette = 'Paired') +
+      ggplot2::theme_minimal() +
+      ggplot2::theme(legend.position = 'bottom', 
+                     legend.title = ggplot2::element_blank()) +
+      ggplot2::ggtitle(title)
+  }
+  return(myplot)
+}

dpf-master.Rproj

@@ -0,0 +1,17 @@
+Version: 1.0
+
+RestoreWorkspace: Default
+SaveWorkspace: Default
+AlwaysSaveHistory: Default
+
+EnableCodeIndexing: Yes
+UseSpacesForTab: Yes
+NumSpacesForTab: 2
+Encoding: UTF-8
+
+RnwWeave: Sweave
+LaTeX: pdfLaTeX
+
+BuildType: Package
+PackageUseDevtools: Yes
+PackageInstallArgs: --no-multiarch --with-keep.source

extras/cluster_analysis_my_model.R

@@ -1,22 +1,22 @@
-# 1. This file estimates parameters for all performances in parallel using
-# a cluster. See https://mllg.github.io/batchtools/ for setup instructions
-# to run on your own machine or an appropriate managed HPC system.
-# Note that, as specified, each performance takes around 5 hours on my laptop.
-
-library(dpf)
-library(batchtools)
-library(dplyr)
-library(optimr) # cran version
-# library(optimrx) # dev version (installed on my laptop)
-# Note: optimrx correctly handles inability to find an optimum within multistart(), optimr does not
-# fix is to add badval=1e8 to the control list, see optimizer() in djm_optimization.R
-
-data(tempos)
-lt = diff(c(tempos$note_onset,61))
-source("my_model.R")
-makeRegistry("my-mazurka", packages=c('dpf','optimr'), 
-             source = "my_model.R",
-             seed = 20181109)
-batchMap(optimizer, as.list(select(tempos, -meas_num, -note_onset, -beat)),
-         more.args = list(lt=lt))
-#submitJobs(resources = list(ppn=1, nodes=1, memory='16gb', walltime='24:00:00'))
+# 1. This file estimates parameters for all performances in parallel using
+# a cluster. See https://mllg.github.io/batchtools/ for setup instructions
+# to run on your own machine or an appropriate managed HPC system.
+# Note that, as specified, each performance takes around 5 hours on my laptop.
+
+library(dpf)
+library(batchtools)
+library(dplyr)
+library(optimr) # cran version
+# library(optimrx) # dev version (installed on my laptop)
+# Note: optimrx correctly handles inability to find an optimum within multistart(), optimr does not
+# fix is to add badval=1e8 to the control list, see optimizer() in djm_optimization.R
+
+data(tempos)
+lt = diff(c(tempos$note_onset,61))
+source("my_model.R")
+makeRegistry("my-mazurka", packages=c('dpf','optimr'), 
+             source = "my_model.R",
+             seed = 20181109)
+batchMap(optimizer, as.list(select(tempos, -meas_num, -note_onset, -beat)),
+         more.args = list(lt=lt))
+#submitJobs(resources = list(ppn=1, nodes=1, memory='16gb', walltime='24:00:00'))

extras/cluster_analysis_my_model_dynamics.R

@@ -0,0 +1,22 @@
+# 1. This file estimates parameters for all performances in parallel using
+# a cluster. See https://mllg.github.io/batchtools/ for setup instructions
+# to run on your own machine or an appropriate managed HPC system.
+# Note that, as specified, each performance takes around 5 hours on my laptop.
+
+library(dpf)
+library(batchtools)
+library(dplyr)
+library(optimr) # cran version
+# library(optimrx) # dev version (installed on my laptop)
+# Note: optimrx correctly handles inability to find an optimum within multistart(), optimr does not
+# fix is to add badval=1e8 to the control list, see optimizer() in djm_optimization.R
+
+data(dynamics)
+lt = diff(c(tempos$note_onset,61))
+source("my_model_dynamics.R")
+makeRegistry("my-mazurka-dynamics", packages=c('dpf','optimr'), 
+             source = "my_model_dynamics.R",
+             seed = 20200405)
+batchMap(optimizer, as.list(select(dynamics, -meas_num, -note_onset, -beat)),
+         more.args = list(lt=lt))
+submitJobs(resources = list(ppn=1, nodes=1, memory='16gb', walltime='48:00:00'))

extras/collect_cluster_results_dynamics.R

@@ -0,0 +1,15 @@
+#2. This short script collects all parameter estimates from the cluster 
+# (or local) when run via cluster_analysis_my_model.R
+
+library(batchtools)
+loadRegistry('my-mazurka-dynamics', writeable=TRUE)
+
+getStatus()
+
+outL = reduceResultsList()
+fun = function(x) x[which.min(x$value),]
+pvec_ml = t(sapply(outL,fun))
+rownames(pvec_ml) = names(dynamics[-c(1:3)])
+pvec_ml = as.data.frame(pvec_ml)
+
+save(pvec_ml, file='mazurkaDynamicsResults.Rdata')

extras/my_model_dynamics.R

@@ -0,0 +1,100 @@
+# Source file for prior and likelihood functions. Not needed separately.
+
+logprior <- function(theta,yt){
+  p2s = c(theta[9],1-sum(theta[c(9,10,11)]),theta[c(10,11)])
+  p2 = ddirichlet(p2s, alpha=c(5,85,5,5),log=TRUE)
+  p4 = dbeta(theta[12],2,8, log=TRUE)
+  muc = dgamma(theta[1], shape=100, scale=.1, log = TRUE)
+  sig2eps = dgamma(theta[2], shape=10, scale=.5, log = TRUE)
+  mu0 = dnorm(theta[3], mean(yt), sd=10, log = TRUE)
+  mu1 = dnorm(theta[4], 0, sd=.25, log = TRUE)
+  mu2 = dnorm(theta[5], 0, sd=.25, log = TRUE)
+  sig0 = dgamma(theta[6], shape=10, scale=1, log=TRUE)
+  sig1 = dgamma(theta[7], shape=3, scale=1, log=TRUE)
+  sig2 = dgamma(theta[8], shape=3, scale=1, log=TRUE)
+  mue = dgamma(-theta[13], shape=100, scale=.1, log = TRUE)
+  lp = sum(muc, sig2eps, mu0, mu1, mu2, sig0, sig1, sig2, 
+           p2, p4, mue)
+  lp
+}
+
+prior_means <- function(yt){
+  c(10, 5,
+    mean(yt), 0, 0,
+    10, 3, 3,
+    0.05,0.05,.05,0.2,
+    -10)
+}
+
+rprior <- function(n,yt){
+  muc = rgamma(n, shape=100, scale=.1)
+  sig2eps = rgamma(n, shape=10, scale=.5)
+  mu0 = rnorm(n, mean=mean(yt), sd = 10)
+  mu1 = rnorm(n, mean=0, sd = .25)
+  mu2 = rnorm(n, mean=0, sd = .25)
+  sig0 = rgamma(n, shape=10, scale=1)
+  sig1 = rgamma(n, shape=3, scale=1)
+  sig2 = rgamma(n, shape=3, scale=1)
+  p2 = rdirichlet(n, alpha=c(5,85,5,5))
+  p21 = p2[1]
+  p23 = p2[3]
+  p24 = p2[4]
+  p41 = rbeta(n,2,8)
+  mue = -rgamma(n, shape=100, scale=.1)
+  cbind(muc, sig2eps, mu0, mu1, mu2, sig0, sig1, sig2, 
+        p21, p23, p24, p41, mue)
+}
+
+
+
+logStatesGivenParams <- function(states,transProbs){
+  ind = cbind(states[1:(length(states)-1)], states[2:length(states)]) + 1
+  return(sum(log(transProbs[ind])))
+}
+
+
+toOptimize <- function(theta, yt, lt, Npart, badvals=Inf){
+  if(any(theta[c(2,6:8)]<0)){
+    cat('do not allow negative variance')
+    return(badvals)
+  }
+  pmats = musicModeldynamics(lt, theta[1], theta[2], theta[3:5], theta[6:8], theta[9:12], theta[13],
+                             c(18,0,0), 
+                             c(10,2,2))
+  beam = beamSearch(pmats$a0, pmats$P0, c(1,0,0,0), 
+                    pmats$dt, pmats$ct, pmats$Tt, pmats$Zt,
+                    pmats$HHt, pmats$GGt, yt, pmats$transMat, Npart, samplemethod = 1)
+  if(beam$LastStep < length(lt)){
+    cat('beam$LastStep < length(lt)\n')
+    return(badvals)
+  }
+  if(all(is.na(beam$weights))){
+    cat('all weights are NA\n')
+    return(badvals)
+  }
+  states = beam$paths[which.max(beam$weights),]
+  negllike = getloglike(pmats, states, yt)# -log(P(y|params, states))
+  sgp = -1 * logStatesGivenParams(states, pmats$transMat)
+  logp = -1 * logprior(theta,yt)
+  obj = negllike + logp + sgp
+  obj
+}
+
+
+# Cluster funs ------------------------------------------------------------
+
+optimizer <- function(perf, lt, Npart=500, ntries = 10, badvals=1e8){
+  yt = matrix(perf, nrow=1)
+  randos = NULL
+  if(ntries > 1) randos = rprior(ntries-1,yt)
+  init_vals = rbind(prior_means(yt), randos)
+  out1 = multistart(init_vals, toOptimize, yt=yt, lt=lt, Npart=Npart, 
+                    badvals=badvals,
+                    method='Nelder-Mead',
+                    control=list(trace=0, maxit=5000, badval=badvals))
+  out2 = multistart(init_vals, toOptimize, yt=yt, lt=lt, Npart=Npart, 
+                    badvals=badvals,
+                    method='SANN',
+                    control=list(trace=0, maxit=5000,badval=badvals))
+  out = rbind.data.frame(out1, out2)
+}