## Software License Agreement (BSD License)
##
## Copyright (c) 2014, Ross Linscott (rossklin@gmail.com)
##
## Redistribution and use in source and binary forms, with or without
## modification, are permitted provided that the following conditions are
## met:
##
## Redistributions of source code must retain the above copyright
## notice, this list of conditions and the following disclaimer.
##
## Redistributions in binary form must reproduce the above copyright
## notice, this list of conditions and the following disclaimer in
## the documentation and/or other materials provided with the
## distribution.
##
## The names of its contributors may not be used to endorse or promote products
## derived from this software without specific prior written permission.
##
## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
## "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
## LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
## A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
## HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
## SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
## LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
## DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
## THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
## (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#' #' LPoly system builder
#' @param mcm model coefficient matrix
#' @param mspec model specification, as created by lpoly_model_spec
lpoly_make_system <- function(mcm, mspec){
lsys <- lpoly_make_system_xptr(mcm, mspec)
attr(lsys, "mcm") <- mcm
attr(lsys, "mspec") <- mspec
attr(lsys, "dimension") <- ncol(mspec)
lsys
}
#' LPoly system modelfun
#' @param sys lpoly_system_type object created with lpoly_make_system
lpoly_system_modelfun <- function(sys){
function(x){
m <- lpoly_model_matrix(sys, x)
colnames(m) <- rownames(attr(sys, "mspec"))
m
}
}
#' LPoly system evalfun
#' @param sys lpoly_system_type object created with lpoly_make_system
lpoly_system_evalfun <- function(sys){
function(x) {
m <- lpoly_model_matrix(sys, x) %*% t(attr(sys, "mcm"))
colnames(m) <- rownames(attr(sys, "mcm"))
m
}
}
#' LPoly system jacobian
#' @param sys lpoly_system_type object created with lpoly_make_system
lpoly_system_jacobian <- function(sys){
function(x) lpoly_compute_jacobian(sys, x)
}
#' LPoly system specs
#' @param sys lpoly_system_type object created with lpoly_make_system
lpoly_system_specs <- function(sys){
res <- NULL
res$mcm <- attr(sys, "mcm")
res$mspec <- attr(sys, "mspec")
res$dimension <- attr(sys, "dimension")
res
}
#' LPoly model spec builder
#' @param d number of measurement variables
#' @param ord array of polynomial orders to include in the model
#' @export
lpoly_model_spec <- function(d, ord){
x <- data.frame(unique(c(0:max(ord), ord)))
s <- expand.grid(rep(x, d))
s <- s[rowSums(abs(s)) %in% abs(ord),,drop=F]
rownames(s) <- alply(s, 1, function(v) paste0(v, collapse = "."))
colnames(s) <- paste0("V", seq_len(d))
as.matrix(s)
}
#' LPoly steplength guesser: sampling method
#'
#' (this function is still under development)
#' @export
lpoly_guess_steplength_sample <- function(sys, goal = 0.1, freq = 400, h0 = 1){
h <- 2 * h0
ml2 <- 1
steps <- 10
f <- lpoly_system_evalfun(sys)
while(ml2 > goal){
h <- h / 2
print.noquote(sprintf("Testing h = %e...", h))
tt <- synthetic.dataset.quick(sys = sys
, num.entities = 10
, steps = steps * as.integer(freq * h + 1)
, tmax = steps * h
, at.times = seq(0, steps*h, h)
, observation.noise.sd = 0
, process.noise.sd = 0)
z <- as.matrix(measurement(step_derivatives(tt)))
x <- as.matrix(measurement(tt))
y <- f(x)
r <- y - z
l2e <- sqrt(rowSums(r^2) / rowSums(y^2))
l2e[is.na(l2e)] <- 1
ml2 <- median(l2e)
print.noquote(sprintf("ml2 = %f", ml2))
}
h
}
#' LPoly steplength guesser: linear method
#'
#' (this function is still under development)
#' @export
lpoly_guess_steplength_linear <- function(sys){
f <- lpoly_system_evalfun(sys)
j <- lpoly_system_jacobian(sys)
d <- nrow(lpoly_system_specs(sys)$mcm)
n <- 100
x <- matrix(rnorm(n * d), n, d)
js <- alply(1:n, 1, function(x) j(matrix(x, 1, d)))
ls <- laply(js, function(m) eigen(m)$values)
ys <- aaply(1:n, 1, function(i) log(2 * f(matrix(x[i,], 1, d)) / min(abs(ls[i,]))^2) / min(abs(ls[i,])))
h <- 0.1 * quantile(ys, 0.1)
h
}
#' LPoly random dynamics generator
#'
#' Generates an lpoly_system_type object representing a polynomial dynamic
#' with randomly generated coefficients.
#' @param d Dimension of the system
#' @param ord Polynomial orders to include: array
#' @param scale.factor Scaling for normally distributed polynomial coefficients
#' @param damp.factor Scaling for the damping terms
#' @param pzero Probability of removing terms from the set of possible terms
#' @param porder.decay Exponential reduction of coefficients w.r.t. order
#' @export
lpoly_random_dynamic <- function(d, ord = 1:2, scale.factor = 1, damp.factor = 0.02, pzero = 0.8, model.size = NULL, porder.decay = 0.5){
next_odd <- function(x) 2 * as.integer(x / 2) + 1
ms <- lpoly_model_spec(d, ord)
if (!(nrow(ms) > 0)) stop("Can not generate dynamic without model spec!")
if (is.null(model.size)){
use.rows <- 0
while (sum(use.rows) < 1) use.rows <- which(runif(nrow(ms), 0, 1) > pzero)
}else {
use.rows <- sample(seq_len(nrow(ms)), model.size)
}
ms <- ms[use.rows,,drop=F]
msB <- diag(rep(next_odd(max(ord) + 1), d))
msC <- diag(rep(-next_odd(max(-ord) + 1), d))
mspec <- ms
msnames <- rownames(ms)
if (max(ord) > 0) {
mspec <- rbind(mspec, msB)
msnames <- c(msnames, paste0("neg.feed.", seq_len(d)))
}
if (min(ord) < 0) {
mspec <- rbind(mspec, msC)
msnames <- c(msnames, paste0("pos.feed.", seq_len(d)))
}
rownames(mspec) <- msnames
colnames(mspec) <- paste0("u.", seq_len(d))
A <- matrix(scale.factor * rnorm(nrow(ms) * d), d, nrow(ms))
## Yes, R is wonderfully concistent, that's why we have to write
## special case clauses for different sizes of arrays...
if (ncol(A) > 1){
A <- t(aaply(seq_len(ncol(A)), 1, function(cidx) porder.decay^(sum(abs(ms[cidx,])) - 1) * A[,cidx]))
}else{
A <- porder.decay^sum(abs(ms)) * A
}
B <- diag(rep(-damp.factor, d))
C <- diag(rep(damp.factor, d))
mcm <- A
if (max(ord) > 0) mcm <- cbind(mcm, B)
if (min(ord) < 0) mcm <- cbind(mcm, C)
rownames(mcm) <- colnames(mspec)
colnames(mcm) <- rownames(mspec)
lsys <- lpoly_make_system(mcm, mspec)
}
#' LPoly example dynamic: lorenz system
#' @export
lpoly_examples_lorenz <- function(s = 16, r = 45.6, b = 5){
cm = t(matrix(c( -s, s, 0, 0, 0, 0,
r, -1, 0, 0, -20, 0,
0, 0, -b, 5, 0, 0), ncol = 3))
trm = t(matrix(c( 1,0,0,
0,1,0,
0,0,1,
1,1,0,
1,0,1,
0,1,1), nrow = 3))
rownames(trm) <- c("x", "y", "z", "xy", "xz", "yz")
colnames(trm) <- c("dx", "dy", "dz")
rownames(cm) <- colnames(trm)
colnames(cm) <- rownames(trm)
lpoly_make_system(cm, trm)
}
examples.gensys.det.lorenz <- function(u,t){
## lorenz system proposed by Chi-Chung Chen and Kung Yao
## in STOCHASTIC CALCULUS NUMERICAL EVALUATION OF CHAOTIC COMMUNICATION SYSTEM PERFORMANCE
## note: behaviour seems stable with respect to each parameter value
sigma <- 16
r <- 45.6
b <- 5
x <- u[,1]
y <- u[,2]
z <- u[,3]
du <- cbind(sigma * (y - x), r * x - y - 20 * x * z, 5 * x * y - b * z)
}
attr(examples.gensys.det.lorenz, "description") <- list("dx/dt = 16 y - 16 x", "dy/dt = 45.6 x - y - 20 x z", "dz/dt = 5 x y - 5 z")
examples.gensys.jacob.lorenz <- function(u_,t){
u <- as.numeric(u_)
sigma <- 16
r <- 45.6
b <- 5
x <- u[1]
y <- u[2]
z <- u[3]
t(matrix(c( -sigma, sigma, 0
, r - 20 * z, -1, -20 * x
, 5 * y, 5 * x, -b), nrow = 3, ncol = 3))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.