Running examples in ‘MOSAlloc-Ex.R’ failed
The error most likely occurred in:

> ### Name: mosalloc
> ### Title: Multiobjective sample allocation for constraint multivariate and
> ###   multidomain optimal allocation in survey sampling
> ### Aliases: mosalloc
> 
> ### ** Examples
> 
> # Artificial population of 50 568 business establishments and 5 business
> # sectors (data from Valliant, R., Dever, J. A., & Kreuter, F. (2013).
> # Practical tools for designing and weighting survey samples. Springer.
> # https://doi.org/10.1007/978-1-4614-6449-5, Example 5.2 pages 133-9)
> 
> # See also <https://umd.app.box.com/s/9yvvibu4nz4q6rlw98ac/file/297813512360>
> # file: Code 5.3 constrOptim.example.R
> 
> Nh <- c(6221, 11738, 4333, 22809, 5467) # stratum sizes
> ch <- c(120, 80, 80, 90, 150) # stratum-specific cost of surveying
> 
> # Revenues
> mh.rev <- c(85, 11, 23, 17, 126) # mean revenue
> Sh.rev <- c(170.0, 8.8, 23.0, 25.5, 315.0) # standard deviation revenue
> 
> # Employees
> mh.emp <- c(511, 21, 70, 32, 157) # mean number of employees
> Sh.emp <- c(255.50, 5.25, 35.00, 32.00, 471.00) # std. dev. employees
> 
> # Proportion of estabs claiming research credit
> ph.rsch <- c(0.8, 0.2, 0.5, 0.3, 0.9)
> 
> # Proportion of estabs with offshore affiliates
> ph.offsh <- c(0.06, 0.03, 0.03, 0.21, 0.77)
> 
> budget <- 300000 # overall available budget
> n.min  <- 100 # minimum stratum-specific sample size
> 
> # Examples
> #----------------------------------------------------------------------------
> # Example 1: Minimization of the variation of estimates for revenue subject
> # to cost restrictions and precision restrictions to the coefficient of
> # variation of estimates for the proportion of businesses with offshore
> # affiliates.
> 
> l <- rep(n.min, 5) # minimum sample size per stratum
> u <- Nh            # maximum sample size per stratum
> C <- rbind(ch,
+            ch * c(-1, -1, -1, 0, 0))
> c <- c(budget,                        # Maximum overall survey budget
+        - 0.5 * budget)                # Minimum overall budget for strata 1-3
> 
> # We require at maximum 5 % relative standard error for estimates of
> # proportion of businesses with offshore affiliates
> A <- matrix(ph.offsh * (1 - ph.offsh) * Nh**3/(Nh - 1)/sum(Nh * ph.offsh)**2,
+  nrow = 1)
> a <- sum(ph.offsh * (1 - ph.offsh) * Nh**2/(Nh - 1)
+ )/sum(Nh * ph.offsh)**2 + 0.05**2
> 
> D <- matrix(Sh.rev**2 * Nh**2, nrow = 1) # objective variance components
> d <- sum(Sh.rev**2 * Nh) # finite population correction
> 
> opts = list(sense = "max_precision",
+             f = NULL, df = NULL, Hf = NULL,
+             init_w = 1,
+             mc_cores = 1L, pm_tol = 1e-05,
+             max_iters = 100L, print_pm = FALSE)
> 
> sol <- mosalloc(D = D, d = d, A = A, a = a, C = C, c = c, l = l, u = u,
+                 opts = opts)
> 
> # Check solution statement of the internal solver to verify feasibility
> sol$Ecosolver$Ecoinfostring # [1] "Optimal solution found"
[1] "Optimal solution found"
> 
> # Check constraints
> c(C[1, ] %*% sol$n) # [1] 3e+05
[1] 3e+05
> c(C[2, ] %*% sol$n) # [1] -150000
[1] -150000
> c(sqrt(A %*% (1 / sol$n) - A %*% (1 / Nh))) # 5 % rel. std. err.
[1] 0.05
> 
> #----------------------------------------------------------------------------
> # Example 2: Minimization of the maximum relative variation of estimates for
> # the total revenue, the number of employee, the number of businesses claimed
> # research credit, and the number of businesses with offshore affiliates
> # subject to cost restrictions
> 
> l <- rep(n.min, 5) # minimum sample size ber stratum
> u <- Nh            # maximum sample size per stratum
> C <- rbind(ch, ch * c(-1, -1, -1, 0, 0))
> c <- c(budget, - 0.5 * budget)
> A <- NULL # no precision constraint
> a <- NULL # no precision constraint
> 
> # Precision components (Variance / Totals^2) for multidimensional objective
> D <- rbind(Sh.rev**2 * Nh**2/sum(Nh * mh.rev)**2,
+            Sh.emp**2 * Nh**2/sum(Nh * mh.emp)**2,
+            ph.rsch * (1 - ph.rsch) * Nh**3/(Nh - 1)/sum(Nh * ph.rsch)**2,
+            ph.offsh * (1 - ph.offsh) * Nh**3/(Nh - 1)/sum(Nh * ph.offsh)**2)
> 
> d <- as.vector(D %*% (1 / Nh)) # finite population correction
> 
> opts = list(sense = "max_precision",
+             f = NULL, df = NULL, Hf = NULL,
+             init_w = 1,
+             mc_cores = 1L, pm_tol = 1e-05,
+             max_iters = 100L, print_pm = FALSE)
> 
> sol <- mosalloc(D = D, d = d, C = C, c = c, l = l, u = u, opts = opts)
> 
> # Obtain optimal objective value
> sol$J # [1] 0.0017058896 0.0004396972 0.0006428475 0.0017058896
[1] 0.0017058896 0.0004396972 0.0006428475 0.0017058896
> 
> # Obtain corresponding normal vector
> sol$Normal # [1] 6.983113e-01 1.337310e-11 1.596167e-11 3.016887e-01
[1] 6.983113e-01 1.337427e-11 1.596306e-11 3.016887e-01
> 
> # => Revenue and offshore affiliates are dominating the solution with a
> #    ratio of approximately 2:1 (sol$Normal[1] / sol$Normal[4])
> 
> #----------------------------------------------------------------------------
> # Example 3: Example 2 with preference weighting
> 
> w <- c(1, 3.85, 3.8, 1.3) # preference weighting
> l <- rep(n.min, 5) # minimum sample size ber stratum
> u <- Nh            # maximum sample size per stratum
> C <- rbind(ch, ch * c(-1, -1, -1, 0, 0))
> c <- c(budget, - 0.5 * budget)
> A <- NULL # no precision constraint
> a <- NULL # no precision constraint
> 
> D <- rbind(Sh.rev**2 * Nh**2/sum(Nh * mh.rev)**2,
+            Sh.emp**2 * Nh**2/sum(Nh * mh.emp)**2,
+            ph.rsch * (1 - ph.rsch) * Nh**3/(Nh - 1)/sum(Nh * ph.rsch)**2,
+            ph.offsh * (1 - ph.offsh) * Nh**3/(Nh - 1)/sum(Nh * ph.offsh)**2)
> 
> d <- as.vector(D %*% (1 / Nh))
> 
> opts = list(sense = "max_precision",
+             f = NULL, df = NULL, Hf = NULL,
+             init_w = w,
+             mc_cores = 1L, pm_tol = 1e-05,
+             max_iters = 100L, print_pm = FALSE)
> 
> mosalloc(D = D, d = d, C = C, c = c, l = l, u = u, opts = opts)
$w
[1] 1.00 3.85 3.80 1.30

$n
[1] 755.8117 499.7609 241.5215 691.6864 584.9882

$J
[1] 0.0018907289 0.0004917598 0.0004982304 0.0014563656

$Objective
NULL

$Utiopian
NULL

$Normal
[1] 2.364034e-08 2.237152e+00 1.186254e+00 1.387742e-01

$dfJ
NULL

$Sensitivity
$Sensitivity$D
[1] 2.364034e-08 5.810784e-01 3.121722e-01 1.067494e-01

$Sensitivity$A
numeric(0)

$Sensitivity$C
[1]  1.050676e-08 -7.307475e-09

$Sensitivity$lbox
[1] 9.094066e-17 1.547199e-16 4.498679e-16 1.020057e-16 1.253216e-16

$Sensitivity$ubox
[1] 1.129214e-17 5.617962e-18 1.503178e-17 2.807553e-18 1.252124e-17


$Qbounds
[1] 1 1 1 1 1

$Dbounds
[1] 0.001890729 0.001893275 0.001893275 0.001893275

$Scalingfactor
[1] 0.0016654338 0.0020770702 0.0022015450 0.0002711488 0.0005464200

$Ecosolver
$Ecosolver$Ecoinfostring
[1] "Optimal solution found"

$Ecosolver$Ecoredcodes
exitFlag     iter  mi_iter   numerr 
       0       12       -1        0 

$Ecosolver$Ecosummary
       pcost        dcost         pres         dres         pinf         dinf 
1.069071e+00 1.069071e+00 1.187801e-12 4.864704e-12 0.000000e+00 0.000000e+00 
     pinfres      dinfres          gap       relgap           r0 
7.107497e-02          NaN 2.309588e-09 2.160369e-09 1.000000e-10 


$Timing
     TotalTime InnerTime ECOS_runtime ECOS_tsetup ECOS_tsolve
[1,]     0.001         0   8.7458e-05   7.375e-06  8.0083e-05

$Iteration
NULL

> 
> #----------------------------------------------------------------------------
> # Example 4: Example 2 with multiple preference weightings for simultaneous
> # evaluation
> 
> w <- matrix(c(1.0, 1.0, 1.0, 1.0,       # matrix of preference weightings
+               1.0, 3.9, 3.9, 1.3,
+               0.8, 4.2, 4.8, 1.5,
+               1.2, 3.5, 4.8, 2.0,
+               2.0, 1.0, 1.0, 2.0), 5, 4, byrow = TRUE)
> w <- w / w[,1]     # rescale w (ensure the first weighting to be one)
> l <- rep(n.min, 5) # minimum sample size ber stratum
> u <- Nh            # maximum sample size per stratum
> C <- rbind(ch, ch * c(-1, -1, -1, 0, 0))
> c <- c(budget, - 0.5 * budget)
> A <- NULL # no precision constraint
> a <- NULL # no precision constraint
> 
> D <- rbind(Sh.rev**2 * Nh**2/sum(Nh * mh.rev)**2,
+            Sh.emp**2 * Nh**2/sum(Nh * mh.emp)**2,
+            ph.rsch * (1 - ph.rsch) * Nh**3/(Nh - 1)/sum(Nh * ph.rsch)**2,
+            ph.offsh * (1 - ph.offsh) * Nh**3/(Nh - 1)/sum(Nh * ph.offsh)**2)
> 
> d <- as.vector(D %*% (1 / Nh))
> 
> opts = list(sense = "max_precision",
+             f = NULL, df = NULL, Hf = NULL,
+             init_w = w,
+             mc_cores = 1L, pm_tol = 1e-05,
+             max_iters = 100L, print_pm = FALSE)
> 
> sols <- mosalloc(D = D, d = d, C = C, c = c, l = l, u = u, opts = opts)
> lapply(sols, function(sol){sol$Qbounds})
[[1]]
[1] 1 1 1 1 1

[[2]]
[1] 1 1 1 1 1

[[3]]
[1] 1 1 1 1 1

[[4]]
[1] 1 1 1 1 1

[[5]]
[1] 1 1 1 1 1

> 
> #----------------------------------------------------------------------------
> # Example 5: Example 2 where a weighted sum scalarization of the objective
> # components is minimized
> 
> l <- rep(n.min, 5) # minimum sample size ber stratum
> u <- Nh            # maximum sample size per stratum
> C <- matrix(ch, nrow = 1)
> c <- budget
> A <- NULL # no precision constraint
> a <- NULL # no precision constraint
> 
> # Objective variance components
> D <- rbind(Sh.rev**2 * Nh**2/sum(Nh * mh.rev)**2,
+            Sh.emp**2 * Nh**2/sum(Nh * mh.emp)**2,
+            ph.rsch * (1 - ph.rsch) * Nh**3/(Nh - 1)/sum(Nh * ph.rsch)**2,
+            ph.offsh * (1 - ph.offsh) * Nh**3/(Nh - 1)/sum(Nh * ph.offsh)**2)
> 
> d <- as.vector(D %*% (1 / Nh)) # finite population correction
> 
> # Simple weighted sum as decision functional
> wss <- c(1, 1, 0.5, 0.5) # preference weighting (weighted sum scalarization)
> 
> Dw <- wss %*% D
> dw <- as.vector(wss %*% d)
> 
> opts = list(sense = "max_precision",
+             f = NULL, df = NULL, Hf = NULL,
+             init_w = 1,
+             mc_cores = 1L, pm_tol = 1e-05,
+             max_iters = 1000L, print_pm = FALSE)
> 
> # Solve weighted sum scalarization (WSS) via mosalloc
> sol_wss <- mosalloc(D = Dw, d = dw, C = C, c = c, l = l, u = u, opts = opts)
> 
> # Obtain optimal objective values
> J <- D %*% (1 / sol_wss$n) - d
> 
> # Reconstruct solution via a weighted Chebyshev minimization
> wcm <- J[1] / J
> opts = list(sense = "max_precision",
+             f = NULL, df = NULL, Hf = NULL,
+             init_w = matrix(wcm, 1),
+             mc_cores = 1L, pm_tol = 1e-05,
+             max_iters = 1000L, print_pm = FALSE)
> 
> sol_wcm <- mosalloc(D = D, d = d, C = C, c = c, l = l, u = u, opts = opts)
> 
> # Compare solutions
> rbind(t(J), sol_wcm$J)
           [,1]         [,2]         [,3]        [,4]
[1,] 0.00155645 0.0004037429 0.0005934474 0.001327165
[2,] 0.00155645 0.0004037429 0.0005934474 0.001327165
> #            [,1]         [,2]         [,3]        [,4]
> # [1,] 0.00155645 0.0004037429 0.0005934474 0.001327165
> # [2,] 0.00155645 0.0004037429 0.0005934474 0.001327165
> 
> rbind(sol_wss$n, sol_wcm$n)
         [,1]     [,2]     [,3]     [,4]     [,5]
[1,] 582.8247 236.6479 116.7866 839.5988 841.4825
[2,] 582.8226 236.6475 116.7871 839.5989 841.4841
> #          [,1]     [,2]     [,3]     [,4]     [,5]
> # [1,] 582.8247 236.6479 116.7866 839.5988 841.4825
> # [2,] 582.8226 236.6475 116.7871 839.5989 841.4841
> 
> rbind(wss, sol_wcm$Normal / sol_wcm$Normal[1])
    [,1]      [,2]      [,3]      [,4]
wss    1 1.0000000 0.5000000 0.5000000
       1 0.9976724 0.4997552 0.4997463
> #    [,1]      [,2]      [,3]      [,4]
> #wss    1 1.0000000 0.5000000 0.5000000
> #       1 0.9976722 0.4997552 0.4997462
> 
> #----------------------------------------------------------------------------
> # Example 6: Example 1 with two subpopulations and a p-norm as decision
> # functional
> 
> l <- rep(n.min, 5) # minimum sample size per stratum
> u <- Nh            # maximum sample size per stratum
> C <- rbind(ch, ch * c(-1, -1, -1, 0, 0))
> c <- c(budget, - 0.5 * budget)
> 
> # At maximum 5 % relative standard error for estimates of proportion of
> # businesses with offshore affiliates
> A <- matrix(ph.offsh * (1 - ph.offsh) * Nh**3/(Nh - 1)/sum(Nh * ph.offsh)**2,
+  nrow = 1)
> a <- sum(ph.offsh * (1 - ph.offsh) * Nh**2/(Nh - 1)
+ )/sum(Nh * ph.offsh)**2 + 0.05**2
> 
> D <- rbind((Sh.rev**2 * Nh**2)*c(0,0,1,1,0),
+            (Sh.rev**2 * Nh**2)*c(1,1,0,0,1))# objective variance components
> d <- as.vector(D %*% (1 / Nh)) # finite population correction
> 
> # p-norm solution
> p <- 5 # p-norm
> opts = list(sense = "max_precision",
+             f = function(x) sum(x**p),
+             df = function(x) p * x**(p - 1),
+             Hf = function(x) diag(p * (p - 1) * x**(p - 2)),
+             init_w = 1,
+             mc_cores = 1L, pm_tol = 1e-05,
+             max_iters = 1000L, print_pm = TRUE)
> 
> sol <- mosalloc(D = D, d = d, C = C, c = c, l = l, u = u, opts = opts)

 Projection method:
k =  0 ; stop_crit =  1 ; angle = 88.6434 ; ||J^k-J^{k-1}|| =  NA 
k =  1 ; stop_crit =  1 ; angle = 78.178 ; ||J^k-J^{k-1}|| =  3673247824 
k =  2 ; stop_crit =  0.9953 ; angle = 43.2535 ; ||J^k-J^{k-1}|| =  1138623818 
k =  3 ; stop_crit =  0.8264 ; angle = 12.7467 ; ||J^k-J^{k-1}|| =  782643345 
k =  4 ; stop_crit =  0.0029 ; angle = 0.0244 ; ||J^k-J^{k-1}|| =  529932729 
k =  5 ; stop_crit =  0 ; angle = 0 ; ||J^k-J^{k-1}|| =  988913.9 
> 
> c(sol$Normal/sol$dfJ)/mean(c(sol$Normal/sol$dfJ))
[1] 0.9999971 1.0000029
> # [1] 0.9999972 1.0000028
> 
> #----------------------------------------------------------------------------
> # Example 7: Example 2 with p-norm as decision functional and only one
> # overall cost constraint
> 
> l <- rep(n.min, 5) # minimum sample size ber stratum
> u <- Nh            # maximum sample size per stratum
> C <- matrix(ch, nrow = 1)
> c <- budget
> A <- NULL # no precision constraint
> a <- NULL # no precision constraint
> 
> # Objective precision components
> D <- rbind(Sh.rev**2 * Nh**2/sum(Nh * mh.rev)**2,
+            Sh.emp**2 * Nh**2/sum(Nh * mh.emp)**2,
+            ph.rsch * (1 - ph.rsch) * Nh**3/(Nh - 1)/sum(Nh * ph.rsch)**2,
+            ph.offsh * (1 - ph.offsh) * Nh**3/(Nh - 1)/sum(Nh * ph.offsh)**2)
> 
> d <- as.vector(D %*% (1 / Nh)) # finite population correction
> 
> # p-norm solution
> p <- 5 # p-norm
> opts = list(sense = "max_precision",
+             f = function(x) sum(x**p),
+             df = function(x) p * x**(p - 1),
+             Hf = function(x) diag(p * (p - 1) * x**(p - 2)),
+             init_w = 1,
+             mc_cores = 1L, pm_tol = 1e-05,
+             max_iters = 1000L, print_pm = TRUE)
> 
> sol <- mosalloc(D = D, d = d, C = C, c = c, l = l, u = u, opts = opts)

 Projection method:
k =  0 ; stop_crit =  0.9999 ; angle = 84.6214 ; ||J^k-J^{k-1}|| =  NA 
k =  1 ; stop_crit =  0.9914 ; angle = 41.3652 ; ||J^k-J^{k-1}|| =  0.0005733031 
k =  2 ; stop_crit =  0.9968 ; angle = 33.8186 ; ||J^k-J^{k-1}|| =  0.0002696888 
k =  3 ; stop_crit =  0.8617 ; angle = 5.5978 ; ||J^k-J^{k-1}|| =  6.351207e-05 
k =  4 ; stop_crit =  0.9511 ; angle = 4.106 ; ||J^k-J^{k-1}|| =  3.079688e-05 
k =  5 ; stop_crit =  1 ; angle = 85.6952 ; ||J^k-J^{k-1}|| =  1.100972e-05 
k =  6 ; stop_crit =  0.8643 ; angle = 1.4772 ; ||J^k-J^{k-1}|| =  5.016101e-06 
k =  7 ; stop_crit =  0.998 ; angle = 54.8637 ; ||J^k-J^{k-1}|| =  8.804734e-08 
k =  8 ; stop_crit =  0.9973 ; angle = 31.039 ; ||J^k-J^{k-1}|| =  1.189502e-09 
k =  9 ; stop_crit =  0.9959 ; angle = 14.6729 ; ||J^k-J^{k-1}|| =  5.344541e-10 
k =  10 ; stop_crit =  0.6858 ; angle = 1.0945 ; ||J^k-J^{k-1}|| =  1.373068e-06 
k =  11 ; stop_crit =  0.9941 ; angle = 16.7977 ; ||J^k-J^{k-1}|| =  8.056878e-09 
k =  12 ; stop_crit =  0.9903 ; angle = 9.0123 ; ||J^k-J^{k-1}|| =  2.421504e-10 
k =  13 ; stop_crit =  0.8457 ; angle = 0.8525 ; ||J^k-J^{k-1}|| =  7.227261e-07 
k =  14 ; stop_crit =  0.9747 ; angle = 4.3227 ; ||J^k-J^{k-1}|| =  7.364256e-10 
k =  15 ; stop_crit =  0.9538 ; angle = 2.3117 ; ||J^k-J^{k-1}|| =  6.755616e-11 
k =  16 ; stop_crit =  0.7225 ; angle = 0.3259 ; ||J^k-J^{k-1}|| =  1.09439e-06 
k =  17 ; stop_crit =  0.7087 ; angle = 0.4905 ; ||J^k-J^{k-1}|| =  3.81818e-11 
k =  18 ; stop_crit =  0.7156 ; angle = 0.2248 ; ||J^k-J^{k-1}|| =  3.765222e-07 
k =  19 ; stop_crit =  0.3662 ; angle = 0.1662 ; ||J^k-J^{k-1}|| =  1.224304e-11 
k =  20 ; stop_crit =  0.7106 ; angle = 0.206 ; ||J^k-J^{k-1}|| =  1.98921e-07 
k =  21 ; stop_crit =  0.2835 ; angle = 0.1257 ; ||J^k-J^{k-1}|| =  1.004688e-11 
k =  22 ; stop_crit =  0.1722 ; angle = 0.1185 ; ||J^k-J^{k-1}|| =  9.162814e-12 
k =  23 ; stop_crit =  0.5073 ; angle = 0.2342 ; ||J^k-J^{k-1}|| =  2.000299e-11 
k =  24 ; stop_crit =  0.1373 ; angle = 0.1152 ; ||J^k-J^{k-1}|| =  1.014242e-11 
k =  25 ; stop_crit =  0.5472 ; angle = 0.2705 ; ||J^k-J^{k-1}|| =  2.676948e-11 
k =  26 ; stop_crit =  0.0794 ; angle = 0.1107 ; ||J^k-J^{k-1}|| =  1.073186e-11 
k =  27 ; stop_crit =  0.624 ; angle = 0.3657 ; ||J^k-J^{k-1}|| =  4.571775e-11 
k =  28 ; stop_crit =  0.0384 ; angle = 0.1084 ; ||J^k-J^{k-1}|| =  1.267456e-11 
k =  29 ; stop_crit =  0.705 ; angle = 0.5248 ; ||J^k-J^{k-1}|| =  8.816326e-11 
k =  30 ; stop_crit =  0.052 ; angle = 0.1065 ; ||J^k-J^{k-1}|| =  1.771919e-11 
k =  31 ; stop_crit =  0.8526 ; angle = 1.2888 ; ||J^k-J^{k-1}|| =  5.118877e-10 
k =  32 ; stop_crit =  0.2859 ; angle = 0.1261 ; ||J^k-J^{k-1}|| =  4.031493e-11 
k =  33 ; stop_crit =  0.2415 ; angle = 0.1263 ; ||J^k-J^{k-1}|| =  9.908977e-12 
k =  34 ; stop_crit =  0.4298 ; angle = 0.1814 ; ||J^k-J^{k-1}|| =  1.343929e-11 
k =  35 ; stop_crit =  0.218 ; angle = 0.1235 ; ||J^k-J^{k-1}|| =  9.602227e-12 
k =  36 ; stop_crit =  0.4749 ; angle = 0.2096 ; ||J^k-J^{k-1}|| =  1.781596e-11 
k =  37 ; stop_crit =  0.0938 ; angle = 0.1116 ; ||J^k-J^{k-1}|| =  8.859693e-12 
k =  38 ; stop_crit =  0.5624 ; angle = 0.2864 ; ||J^k-J^{k-1}|| =  3.305285e-11 
k =  39 ; stop_crit =  0.0978 ; angle = 0.1119 ; ||J^k-J^{k-1}|| =  1.130016e-11 
k =  40 ; stop_crit =  0.5674 ; angle = 0.2919 ; ||J^k-J^{k-1}|| =  3.289195e-11 
k =  41 ; stop_crit =  0.0636 ; angle = 0.1063 ; ||J^k-J^{k-1}|| =  1.000995e-11 
k =  42 ; stop_crit =  0.8856 ; angle = 1.7365 ; ||J^k-J^{k-1}|| =  1.032579e-09 
k =  43 ; stop_crit =  0.4428 ; angle = 0.1888 ; ||J^k-J^{k-1}|| =  5.126245e-11 
k =  44 ; stop_crit =  0.2431 ; angle = 0.1264 ; ||J^k-J^{k-1}|| =  1.051365e-11 
k =  45 ; stop_crit =  0.4869 ; angle = 0.2182 ; ||J^k-J^{k-1}|| =  1.8235e-11 
k =  46 ; stop_crit =  0.0351 ; angle = 0.1069 ; ||J^k-J^{k-1}|| =  8.483267e-12 
k =  47 ; stop_crit =  0.775 ; angle = 0.7596 ; ||J^k-J^{k-1}|| =  1.808586e-10 
k =  48 ; stop_crit =  0.0754 ; angle = 0.106 ; ||J^k-J^{k-1}|| =  2.48278e-11 
k =  49 ; stop_crit =  0.9439 ; angle = 3.8402 ; ||J^k-J^{k-1}|| =  4.699882e-09 
k =  50 ; stop_crit =  0.6613 ; angle = 0.4291 ; ||J^k-J^{k-1}|| =  1.101886e-10 
k =  51 ; stop_crit =  0.0992 ; angle = 0.1054 ; ||J^k-J^{k-1}|| =  1.412875e-11 
k =  52 ; stop_crit =  0.7203 ; angle = 0.166 ; ||J^k-J^{k-1}|| =  7.924969e-07 
k =  53 ; stop_crit =  0.0873 ; angle = 0.0202 ; ||J^k-J^{k-1}|| =  5.594265e-12 
k =  54 ; stop_crit =  0.0034 ; angle = 0.0018 ; ||J^k-J^{k-1}|| =  7.588138e-13 
k =  55 ; stop_crit =  0.0034 ; angle = 0.0057 ; ||J^k-J^{k-1}|| =  0 
k =  56 ; stop_crit =  0.0195 ; angle = 0.0047 ; ||J^k-J^{k-1}|| =  4.37123e-13 
k =  57 ; stop_crit =  0.0181 ; angle = 0.0044 ; ||J^k-J^{k-1}|| =  6.756965e-13 
k =  58 ; stop_crit =  0.0021 ; angle = 0.0017 ; ||J^k-J^{k-1}|| =  8.626247e-13 
k =  59 ; stop_crit =  0.037 ; angle = 0.0086 ; ||J^k-J^{k-1}|| =  1.825068e-12 
k =  60 ; stop_crit =  0.0126 ; angle = 0.0033 ; ||J^k-J^{k-1}|| =  1.094704e-12 
k =  61 ; stop_crit =  0.0021 ; angle = 0.0017 ; ||J^k-J^{k-1}|| =  2.24088e-13 
k =  62 ; stop_crit =  0.0321 ; angle = 0.0075 ; ||J^k-J^{k-1}|| =  3.063634e-12 
k =  63 ; stop_crit =  0.0099 ; angle = 0.0027 ; ||J^k-J^{k-1}|| =  1.929717e-12 
k =  64 ; stop_crit =  0.0022 ; angle = 0.0017 ; ||J^k-J^{k-1}|| =  1.545681e-12 
k =  65 ; stop_crit =  0.0022 ; angle = 0.0055 ; ||J^k-J^{k-1}|| =  0 
k =  66 ; stop_crit =  0.0022 ; angle = 0.0053 ; ||J^k-J^{k-1}|| =  0 
k =  67 ; stop_crit =  0.0022 ; angle = 0.0054 ; ||J^k-J^{k-1}|| =  0 
k =  68 ; stop_crit =  0.0022 ; angle = 0.0032 ; ||J^k-J^{k-1}|| =  0 
Procedur terminated due to numerical accuracy! Verify result!
> 
> c(sol$Normal/sol$dfJ)/mean(c(sol$Normal/sol$dfJ))
[1] 0.9992928 1.0013168 1.0001418 0.9992486
> # [1] 1.0014362 0.9780042 1.0197807 1.0007789
> 
> #----------------------------------------------------------------------------
> # Example 8: Minimization of sample sizes subject to precision constraints
> 
> l <- rep(n.min, 5) # minimum sample size ber stratum
> u <- Nh            # maximum sample size per stratum
> 
> # We require at maximum 4.66 % relative standard error for the estimate of 
> # total revenuee, 5 % for the number of employees, 3 % for the proportion of
> # businesses claiming research credit, and 3 % for the proportion of
> # businesses with offshore affiliates
> A <- rbind(Sh.rev**2 * Nh**2/sum(Nh * mh.rev)**2,
+            Sh.emp**2 * Nh**2/sum(Nh * mh.emp)**2,
+            ph.rsch * (1 - ph.rsch) * Nh**3/(Nh - 1)/sum(Nh * ph.rsch)**2,
+            ph.offsh * (1 - ph.offsh) * Nh**3/(Nh - 1)/sum(Nh * ph.offsh)**2)
> a <- as.vector(A%*%(1 / Nh) + c(0.0466, 0.05, 0.03, 0.03)**2)
> 
> # We do not consider any additional sample size or cost constraints
> C <- NULL # no cost constraint
> c <- NULL # no cost constraint
> 
> # Since we minimize the sample size, we define D and d as follows:
> D <- matrix(1, nrow = 1, ncol = length(Nh)) # objective cost components
> d <- as.vector(0)                           # vector of possible fixed cost
> 
> opts = list(sense = "min_cost", # Sense of optimization is survey cost
+             f = NULL,
+             df = NULL,
+             Hf = NULL,
+             init_w = 1,
+             mc_cores = 1L, pm_tol = 1e-05,
+             max_iters = 100L, print_pm = TRUE)
> 
> sol <- mosalloc(D = D, d = d, A = A, a = a, l = l, u = u, opts = opts)
> 
> sum(sol$n) # [1] 2843.219
[1] 2843.219
> sol$J # [1] 2843.219
[1] 2843.219
> 
> #----------------------------------------------------------------------------
> #----------------------------------------------------------------------------
> # Note: Sample size optimization for two-stage cluster sampling can be
> #       reduced to the structure of optimal stratified random samplin when
> #       considering expected costs. Therefore, mosalloc() can handle such
> #       designs. A benefit is that mosalloc() allows relatively complex
> #       sample size restrictions such as box constraints for subsampling.
> #       Optimal sample sizes at secondary stages have to be reconstructed
> #       from sol$n.
> #
> # Example 9: Optimal number of primary sampling units (PSU) and secondary
> # sampling units (SSU) in 2-stage cluster sampling.
> 
> set.seed(1234)
> pop <- data.frame(value = rnorm(100, 100, 35),
+                   cluster = sample(1:4, 100, replace = TRUE))
> 
> CI <- 36  # Sampling cost per PSU/cluster
> CII <- 10 # Average sampling cost per SSU
> 
> NI <- 4                   # Number of PSUs/clusters
> NII <- table(pop$cluster) # PSU/cluster sizes
> 
> S2I <- var(by(pop$value, pop$cluster, sum)) # between cluster variance
> S2II <- by(pop$value, pop$cluster, var)     # within cluster variances
> 
> D <- matrix(c(NI**2 * S2I - NI * sum(NII * S2II), NI * NII**2 * S2II), 1)
> d <- as.vector(NI * S2I)
> 
> C <- cbind(c(CI, rep(2, NI), -NII),
+            rbind(rep(CII / NI, 4), -diag(4), diag(4)))
> c <- as.vector(c(500, rep(0, 8)))
> 
> l <- c(2, rep(4, 4))
> u <- c(NI, NI * NII)
> 
> opts = list(sense = "max_precision",
+             f = NULL,
+             df = NULL,
+             Hf = NULL,
+             init_w = 1,
+             mc_cores = 1L, pm_tol = 1e-05,
+             max_iters = 100L, print_pm = TRUE)
> 
> sol <- mosalloc(D = D, d = d, C = C, c = c, l = l, u = u, opts = opts)
Error in mosalloc(D = D, d = d, C = C, c = c, l = l, u = u, opts = opts) : 
  d is not an utopian vector! d is too large.
Execution halted