mle2 start parameters don't match lower and upper bounds (using L-BFGS-B)

When you pass lower and upper bounds to mle2 for use with L-BFGS-B, the bounds do not automatically line up with the formula and start values.

Here is some code to illustrate the problem.
```
# Create test data frame with 14 species, 2 canopy levels, and Height measurements every 0.1 meter
spp_levels <- c("abba","acru","acsa","acsp","alin","amel","beal","bepa","coco","cose","frni","pist","poba","potr")
canopy_levels <- c("Open","Shaded")
test_data <- expand.grid(Spp = spp_levels , CanopyClass = canopy_levels , Height = seq(0.1, 20, length.out = 20/0.1) , stringsAsFactors = FALSE)
test_data$Spp <- as.factor(test_data$Spp)
test_data$CanopyClass <- as.factor(test_data$CanopyClass )

# Set test parameters
Tmax_t <- 20
Hmax_t <- 2.5
b_t <- 0.6
overdisp_t <- 0.6

# Create complex power Ricker function with a negative binomial distribution:
# power ricker function:
test_data$twigs <- Tmax_t * Hmax_t^(-b_t*Hmax_t) * exp(b_t*Hmax_t) * test_data$Height^(b_t*Hmax_t) * exp(-b_t*test_data$Height)
# Add negative binomial random noise:
test_data$twigs_r <- rnbinom( n = length(test_data$twigs), mu = test_data$twigs, size = overdisp_t )

plot(twigs_r ~ Height, data = test_data)

# Use mle2 to fit a model:
mod_test <- mle2( twigs_r ~ dnbinom(mu = Tmax * Hmax^(-b*Hmax) * exp(b*Hmax) * Height^(b*Hmax) * exp(-b*Height),
size = overdisp ), data = test_data,
parameters = list(Tmax ~ 1, Hmax ~ 1 , b ~ 1 , overdisp ~ 1 ),
optimizer = "optimx", method="L-BFGS-B", start = list( Tmax = 17 , Hmax = 0.9, b = 0.5, overdisp = 2.2),
lower = c(Tmax=0.1, Hmax = 0.5, b = 0.01, overdisp = 0.01),
upper = c(Tmax=200, Hmax = 10, b = 3.0 , overdisp = 10.0),
control=list( trace = 3, maxit = 100, factr = 1e12), itnmax =100 )

# mod_test converges just fine.

```
Now, if I try to use mle2 to fit more complex model:
```
mod_test2 <- mle2( twigs_r ~ dnbinom(mu = Tmax * Hmax^(-b*Hmax) * exp(b*Hmax) * Height^(b*Hmax) * exp(-b*Height),
size = overdisp ), data = test_data,
parameters = list(Tmax ~ Spp:CanopyClass -1, Hmax ~ Spp -1 , b ~ Spp -1 , overdisp ~ Spp -1 ),
optimizer = "optimx", method="L-BFGS-B", start = list( Tmax = 17 , Hmax = 0.9, b = 0.5, overdisp = 2.2),
lower = c(Tmax=0.1, Hmax = 0.5, b = 0.01, overdisp = 0.01),
upper = c(Tmax=200, Hmax = 10, b = 3.0 , overdisp = 10.0),
control=list( trace = 3, maxit = 100, factr = 1e12), itnmax =100 )
```
If throws an error:
```
fn is fn1
gr NULL, npar > 50, kkt set FALSE
Looking for method = L-BFGS-B
Methods to be used:[1] "L-BFGS-B"
Function has 70 arguments
par[ 1 ]: 0.1

Start by reproducing issue 38 with the mle2 examples, especially the models using expanded parameters and L-BFGS-B bounds. Trace how parameter formulas, start values, lower bounds, and upper bounds are expanded and ordered; done means the reported models accept correctly aligned values without the bounds error.

Written by the indexing model from the issue text.