# tests/extras.R In weiyaw/blackbox: A SMC-SA Algorithm for Constrained Optimisation

```#### Dugong ####
scaled_dugong <- as.data.frame(scale(dugong, FALSE, apply(abs(dugong), 2, max)))
plot(Length ~ Age, data = scaled_dugong, main = "Dugong",
cex = 1.3, cex.lab = 1.3, cex.axis = 1.3)

Bmat <- get_bs_design(scaled_dugong\$Age, 2, 2)\$design

rtnorm_rw_dg <- function(states, k) {
sd <- 2
fraction <- 0.99
sigma <- (sd * fraction^k)^2 * diag(NROW(states))
F <- diag(NROW(states)) * -1
F[row(F) + 1 == col(F)] <- 1
F <- F[1:(NROW(states) - 1), ]
lower <- rep(0, len = NROW(F))
initial <- seq(1, NROW(states))

res <- matrix(NA, NROW(states), NCOL(states))
for (i in 1:NCOL(states)) {
res[, i] <- tnorm::rmvtnorm(1, states[, i], cov = sigma, g = lower, F = F,
initial = initial, burn = 10)
}
res
}

## set sd = 100 and k = 0, and use rtnorm_rw_dg
set.seed(100)
starting <- rtnorm_rw_dg(matrix(rep(1:NCOL(Bmat), times = 1000), ncol = 1000), 1)

## SMCSA, log schedule
dg1_crss <- SMCSA(loss, rtnorm_rw_dg, starting, log_schedule, 2000, 1000, TRUE, TRUE)

## SMCSA, log schedule
dg2_crss <- SMCSA(loss, rtnorm_rw_dg, starting, recip_schedule, 2000, 1000, TRUE, TRUE)

dg_pilot <- solve(crossprod(Bmat)) %*% crossprod(Bmat, scaled_dugong\$Length)

Bmat_plot <- with(scaled_dugong, get_bs_design(seq(min(Age), max(Age), len = 100), 2, 2)\$design)
lines((Bmat_plot %*% dg1_crss\$state) ~ with(scaled_dugong, seq(min(Age), max(Age), len = 100)), col = 'black')
lines((Bmat_plot %*% dg_pilot) ~ with(scaled_dugong, seq(min(Age), max(Age), len = 100)), col = 'red')

#### Rabbit ####
scaled_rabbit <- as.data.frame(scale(rabbit, FALSE, apply(abs(rabbit), 2, max)))
plot(Lens ~ Days, data = scaled_rabbit, main = "Rabbit",
cex = 1.3, cex.lab = 1.3, cex.axis = 1.3)
Bmat <- get_bs_design(scaled_rabbit\$Days, 2, 2)\$design

pilot <- solve(crossprod(Bmat)) %*% crossprod(Bmat, scaled_rabbit\$Lens)
lines((Bmat %*% pilot) ~ scaled_rabbit\$Days, col = 'blue')

## ## get starting values using the rearrange and regress method
## rar_rabbit <- coef(lm(Lens ~ Days + I(-Lens*Days), scaled_rabbit))
## names(rar_rabbit) <- letters[1:3]

## randomly generate 1000 starting values
set.seed(2)
starting <- matrix(rep(rar_rabbit, times = 1000), length(rar_rabbit), 1000) +
rcauchy(1000 * length(rar_rabbit), scale = 1)

## SMCSA, log schedule
rnorm_rw <- rnorm_rw_fac(sd = 2, fraction = 0.99)
rb1_rss <- SMCSA(loss, rnorm_rw, starting, log_schedule, 1000, 500, FALSE, TRUE)