Merge pull request #131 from SMAC-Group/rc1

RC1

.Rbuildignore

@@ -8,3 +8,5 @@ CONTRIBUTING\.md
 README\.md
 cran-comments\.md
 LICENSE
+^README\.Rmd$
+^README-.*\.png$

DESCRIPTION

@@ -2,7 +2,7 @@ Package: gmwm
 Type: Package
 Title: Generalized Method of Wavelet Moments
 Version: 2.0.0
-Date: 2016-01-31
+Date: 2016-02-09
 Authors@R: c(
     person("James", "Balamuta", , "[email protected]", c("aut", "cph")),
     person("Stephane", "Guerrier", , "[email protected]", c("ctb","cre","cph")),

NEWS.md

@@ -1,3 +1,93 @@
+# gmwm 2.0.0
+
+This is the second general release of the Generalized Method of Wavelet Moments (GMWM) R Package. Considerable amount of work has been done to improve the user experience across all disciplines. 
+
+## Highlights
+
+* IMU Data Object Improvements
+* New time series modeling term: GM (Gauss-Markov)
+* More flexible MODWT, DWT, and AVAR Wrappers
+* Improved IMU Parameter Search Algorithm
+* Faster Batch WVAR Processing
+* Rewritten visual model comparison
+* Auto-version check on package load
+
+##  Time Series Processed Model Syntax
+
+* New model term that enables users to specify a guass-markov process (GM).
+  * The GM is a reparameterization of an AR1 process using the freqeuncy of the data.
+* Note: The conversion to and fro GM is done outside the C++ backend. Internally, the term is treated as an AR1 with frequency of 1.
+
+## Decomposition Methods
+
+* New flexible R wrappers for
+    * Discrete Wavelet Transform (DWT)
+    * Maximal Overlap Discrete Wavelet Transform (MODWT)
+    * Allan Variance (AVAR)
+* Support is now available for taking up to a specific scale.
+* DWT and MODWT now have built-in `brick.wall` support.
+* New `brick.wall` R wrapper to later handle boundary coefficients. 
+
+## Wavelet Variance
+
+* Calculate the Wavelet Variance for DWT composition in addition to an MODWT decomposition.
+* Change the level of decomposition for the WV.
+* New batch data set processing function in C++
+
+## Model Comparison
+
+* Reimplemented `compare.models()` function that removes the large spaces between plots.
+* Labels are now around the top and right of the plot.
+
+## Data Objects
+
+* `imu` and `lts` objects now have a frequency, unit, and name component.
+* All data objects were rewritten so that data is placed first and foremost. 
+   * Attributes of the data are now accessed via `attr(x,"name")`
+* Smart printing has been enabled for data objects. No longer will console be flooded with text!
+   * Objects print the first and last 10 observations. 
+* `update.*` function
+   * Added a means to update parameters of an object without recreating it.
+
+## IMU
+
+* New grid search algorithm for IMU parameters that relies upon dominance of the first few scales.
+* Improved loading experience with the addition of `read.imu` that returns an R IMU object type.
+* The `wvar.imu` function now calls a new batch wrapper that processes multiple signals quickly.
+* Added ability to subset IMU without losing attributes. (Limited to pairs of data.)
+
+## Package Load
+
+* Auto-version check on package load
+   * Connects to CRAN to see if a new package exists and alerts the user.
+* More specific startup messages to help users access the resources of the package.
+
+## External Data Package Install
+
+* Established the [SMAC Group Data Repo](http://smac-group.com/datarepo)
+* Rewrote install_imudata to make use of the data repo
+* Removed `dl.file` and `install_` functions.
+
+## Internal
+
+* Added quick `create_imu()` and `create_wvar()` objects to avoid long checks.
+* Added `ar1_to_gm()` and `gm_to_ar1()` to convert to and for the GM object.
+* Added `is.*()` functions for gmwm-based objects to avoid constantly calling `inherit(x,"*")`.
+
+## Bug Fixes
+
+* `freq` now affects the `x`-axis.
+* The transformation with `DR()` is no longer the identity. Instead, we opt to use a logarithmic approach.
+   * In transform_data, DR now has the absolute value taken and then a log.
+   * In untransform_data, DR now is exponentiated.
+* Capitalized axis when displaying `imu` graphs.
+
+## Known Bugs
+* The `compare.wvar` does not work with type `wvar.imu`.
+
+
+------------------------------------------------------------
+
 # gmwm 1.0.0
 
 This is the first general release of the Generalized Method of Wavelet Moments (GMWM) R package. The package is meant to act as a platform for using the GMWM estimator.

R/GMWM.R

@@ -127,9 +127,8 @@
 #'  
 #' # ARMA case
 #' set.seed(1336)
-#' data = arima.sim(n = 200, 
-#'               list(ar = c(0.8897, -0.4858), ma = c(-0.2279, 0.2488)),
-#'               sd = sqrt(0.1796))
+#' data = gen.gts(ARMA(ar = c(0.8897, -0.4858), ma = c(-0.2279, 0.2488),
+#'               sigma2 = 0.1796), 200)
 #' #guided.arma = gmwm(ARMA(2,2), data, model.type="ssm")
 #' adv.arma = gmwm(ARMA(ar=c(0.8897, -0.4858), ma = c(-0.2279, 0.2488), sigma2=0.1796),
 #'                 data, model.type="ssm")
@@ -235,28 +234,25 @@ gmwm = function(model, data, model.type="ssm", compute.v="auto",
 
   # Convert from GM to AR1
   if(!starting && any(model$desc == "GM")){
-    idx = model$process.desc %in% c("BETA","SIGMA2_GM")
-    theta[idx,] = gm_to_ar1(theta[idx,], freq)
+    theta = conv.gm.to.ar1(theta, model$process.desc, freq)
   }
 
-
   out = .Call('gmwm_gmwm_master_cpp', PACKAGE = 'gmwm', data, theta, desc, obj, model.type, starting = model$starting,
                                                          p = alpha, compute_v = compute.v, K = K, H = H, G = G,
                                                          robust=robust, eff = eff)
 
   estimate = out[[1]]
   rownames(estimate) = model$process.desc
   colnames(estimate) = "Estimates" 
-  
+
   init.guess = out[[2]]
   rownames(init.guess) = model$process.desc
   colnames(init.guess) = "Starting" 
 
   # Convert from AR1 to GM
   if(any(model$desc == "GM")){
-    idx = model$process.desc %in% c("BETA","SIGMA2_GM")
-    estimate[idx,] = ar1_to_gm(estimate[idx,],freq)
-    init.guess[idx,] = ar1_to_gm(init.guess[idx,],freq)
+    estimate[,1] = conv.ar1.to.gm(estimate[,1], model$process.desc, freq)
+    init.guess[,1] = conv.ar1.to.gm(init.guess[,1], model$process.desc, freq)
   }
 
   # Wrap this into the C++ Lib
@@ -403,9 +399,8 @@ update.gmwm = function(object, model, ...){
   }
 
   # Convert from GM to AR1
-  if(!model$starting && any(model$desc == "GM")){
-    idx = model$process.desc %in% c("BETA","SIGMA2_GM")
-    model$theta[idx,] = gm_to_ar1(model$theta[idx,], object$freq)
+  if(!object$starting && any(model$desc == "GM")){
+    model$theta = conv.gm.to.ar1(model$theta, model$process.desc, object$freq)
   }
 
   out = .Call('gmwm_gmwm_update_cpp', PACKAGE = 'gmwm',
@@ -434,9 +429,8 @@ update.gmwm = function(object, model, ...){
 
   # Convert from AR1 to GM
   if(any(model$desc == "GM")){
-    idx = model$process.desc %in% c("BETA","SIGMA2_GM")
-    estimate[idx,] = ar1_to_gm(estimate[idx,], object$freq)
-    init.guess[idx,] = ar1_to_gm(init.guess[idx,], object$freq)
+    estimate[,1] = conv.ar1.to.gm(estimate[,1], model$process.desc, object$freq)
+    init.guess[,1] = conv.ar1.to.gm(init.guess[,1], model$process.desc, object$freq)
   }
 
   object$estimate = estimate
@@ -633,8 +627,7 @@ summary.gmwm = function(object, inference = NULL,
 
     # Convert from GM to AR1
     if(any(object$model$desc == "GM")){
-      idx = object$model$process.desc %in% c("BETA","SIGMA2_GM")
-      object$estimate[idx,] = gm_to_ar1(object$estimate[idx,], object$freq)
+      object$estimate[,1] = conv.gm.to.ar1(object$estimate[,1], object$model$process.desc, object$freq)
     }
 
     mm = .Call('gmwm_get_summary', PACKAGE = 'gmwm',object$estimate,
@@ -647,12 +640,12 @@ summary.gmwm = function(object, inference = NULL,
                                                     inference, F, # fullV is always false. Need same logic updates.
                                                     bs.gof, bs.gof.p.ci, bs.theta.est, bs.ci,
                                                     B)
-    # Convert back from AR1 to GM
-    if(any(object$model$process.desc == "BETA")){
-      idx = object$model$process.desc %in% c("BETA","SIGMA2_GM")
-      object$estimate[idx,] = ar1_to_gm(object$estimate[idx,], object$freq)
+    # Convert from AR1 to GM
+    if(any(object$model$desc == "GM")){
+      object$estimate[,1] = conv.ar1.to.gm(object$estimate[,1], object$model$process.desc, object$freq)
     }
 
+
   }else{
     mm = vector('list',3)
     mm[1:3] = NA
@@ -829,12 +822,11 @@ plot.gmwm = function(x, process.decomp = FALSE, background = 'white', CI = T, tr
 #' # AR
 #' set.seed(1336)
 #' n = 200
-#' x = gen_ar1(n, phi=.1, sigma2 = 1) + gen_ar1(n,phi=0.95, sigma2 = .1)
+#' x = gen.gts(AR1(phi = .1, sigma2 = 1) + AR1(phi = 0.95, sigma2 = .1), n)
 #' mod = gmwm(AR1(), data=x, model.type="imu")
 #' autoplot(mod)
 #' 
-#' y = gen.gts(AR1(phi = .1, sigma2 = 1) + AR1(phi = 0.95, sigma2 = .1), n)
-#' mod = gmwm(2*AR1(), data = y)
+#' mod = gmwm(2*AR1(), data = x)
 #' autoplot(mod)
 autoplot.gmwm = function(object, process.decomp = FALSE, background = 'white', CI = T, transparence = 0.1, bw = F, 
                      CI.color = "#003C7D", line.type = NULL, line.color = NULL,

R/dwt.R

@@ -27,6 +27,21 @@
 #' x = rnorm(2^8)
 #' dwt(x)
 dwt = function(x, nlevels = floor(log2(length(x))), filter = "haar", boundary="periodic", bw = TRUE) {
+
+  if(is.vector(x) && length(x) %% nlevels != 0){
+    warning("The data has been truncated so that it is divisible by `nlevels` (e.g. 2^*)")
+    x = x[1:2^nlevels]
+  }else if(is.matrix(x) || is.data.frame(x)){
+    if(ncol(x) != 1){
+      stop("Only one column is allowed to be decomposed at a time.")
+    }
+
+    if(nrow(x) %% nlevels !=0){
+      warning("The data has been truncated so that it is divisible by `nlevels` (e.g. 2^*)")
+      idx = 1:2^nlevels
+      x[idx,1] = x[idx,1]
+    }
+  }
   out = .Call('gmwm_dwt_cpp', PACKAGE = 'gmwm', x, filter_name = filter, nlevels, boundary = boundary, brickwall = bw)
   names(out) = paste0("S",1:nlevels)
   mostattributes(out) = list(J=nlevels, filter = filter, boundary = boundary, brick.wall = bw, class="dwt")

R/gts.R

@@ -139,10 +139,17 @@ gts = function(data, start = 0, end = NULL, freq = 1, unit = NULL, name = NULL){
 #' 
 #' set.seed(1336)
 #' # GM + WN
-#' model = GM(beta = 0.06931472, sigma2_gm = 0.1333333) + WN(sigma2=1)
-#' x = gen.gts(model, n, freq = 10, unit = 'sec')
-#' x
-#' plot(x, to.unit = 'min')
+#' # Convert from AR1 to GM values
+#' m = ar1_to_gm(c(.5,.1),10)
+#' # Beta = 6.9314718, Sigma2_gm = 0.1333333
+#' model = GM(beta = m[1], sigma2_gm = m[2]) + WN(sigma2=1)
+#' x2 = gen.gts(model, n, freq = 10, unit = 'sec')
+#' x2
+#' 
+#' plot(x2, to.unit = 'min')
+#' 
+#' # Same time series
+#' all.equal(x, x2, check.attributes = FALSE)
 gen.gts = function(model, N = 1000, start = 0, end = NULL, freq = 1, unit = NULL, name = NULL){
 
   # 1. Do we have a valid model?
@@ -194,8 +201,7 @@ gen.gts = function(model, N = 1000, start = 0, end = NULL, freq = 1, unit = NULL
 
     # Convert from AR1 to GM
     if(any(model$desc == "GM")){
-      idx = model$process.desc %in% c("BETA","SIGMA2_GM")
-      theta[idx] = ar1_to_gm(theta[idx],freq)
+      theta = conv.gm.to.ar1(theta, model$process.desc, freq)
     }
 
     out = .Call('gmwm_gen_model', PACKAGE = 'gmwm', N, theta, desc, obj)
@@ -207,10 +213,10 @@ gen.gts = function(model, N = 1000, start = 0, end = NULL, freq = 1, unit = NULL
 
   out = structure(.Data = out, 
                   start = start, 
-                  end= end, # start and end will not be null now
-                  freq = freq,
-                  unit = unit,
-                  name = name, 
+                  end   = end, # start and end will not be null now
+                  freq  = freq,
+                  unit  = unit,
+                  name  = name, 
                   class = c("gts","matrix"))
 
   out

R/model_selection_methods.R

@@ -94,9 +94,8 @@ output.format = function(out, model.names, scales, N, alpha, robust, eff, B, G,
 
 
   if(any(model.hat$desc == "GM")){
-    idx = model.hat$process.desc %in% c("BETA","SIGMA2_GM")
-    estimate[idx,] = ar1_to_gm(estimate[idx,],freq)
-    init.guess[idx,] = ar1_to_gm(init.guess[idx,],freq)
+    estimate[,1] = conv.ar1.to.gm(estimate[,1], model.hat$process.desc, freq)
+    init.guess[,1] = conv.ar1.to.gm(init.guess[,1], model.hat$process.desc, freq) 
   }
 
   model.hat$theta = as.numeric(estimate)

R/utilities.R

@@ -55,6 +55,31 @@ packageVersionCRAN = function(pkg, cran_url="http://cran.r-project.org/web/packa
 }
 
 
+#' GM Conversion
+#' 
+#' Convert from AR1 to GM and vice-versa
+#' @param theta        A \code{numeric vector} containing the theta values
+#' @param process.desc A \code{character vector} containing the names of parameters.
+#' @param freq         A \code{double} indicating the frequency of the data.
+#' @keywords internal
+#' @author JJB
+#' @rdname gm_conv
+conv.ar1.to.gm = function(theta, process.desc, freq){
+  idx = process.desc %in% c("BETA","SIGMA2_GM")
+  theta[idx] = ar1_to_gm(theta[idx],freq)
+
+  theta
+}
+
+#' @rdname gm_conv
+conv.gm.to.ar1 = function(theta, process.desc, freq){
+  idx = process.desc %in% c("BETA","SIGMA2_GM")
+  theta[idx] = gm_to_ar1(theta[idx],freq)
+
+  theta
+}
+
+
 #' @title Print GMWM Data Object
 #' @description 
 #' Pretty formatting for \code{gts}, \code{imu}, and \code{lts} objects.