diff -Nru adagio-0.7.1/debian/changelog adagio-0.8.4/debian/changelog --- adagio-0.7.1/debian/changelog 2021-05-02 00:18:11.000000000 +0000 +++ adagio-0.8.4/debian/changelog 2021-05-02 00:16:07.000000000 +0000 @@ -1,15 +1,22 @@ -adagio (0.7.1-1cran1.1604.0) xenial; urgency=medium +adagio (0.8.4-1cran1.1604.0) xenial; urgency=medium - * Compilation for Ubuntu 16.04.6 LTS + * Compilation for Ubuntu 16.04.7 LTS * Build for c2d4u for R 4.0.0 plus - -- Michael Rutter Thu, 04 Jun 2020 20:52:37 +0000 + -- Michael Rutter Sun, 02 May 2021 00:16:07 +0000 + +adagio (0.8.4-1cran1) testing; urgency=low + + * cran2deb svn: 362M with DB version 1. + + -- cran2deb4ubuntu Sat, 01 May 2021 13:48:28 -0400 + adagio (0.7.1-1cran1) testing; urgency=low * cran2deb svn: 362M with DB version 1. - -- cran2deb4ubuntu Sun, 27 May 2018 10:39:48 -0400 + -- cran2deb4ubuntu Sun, 27 May 2018 10:39:57 -0400 adagio (0.6.5-1cran2) testing; urgency=low diff -Nru adagio-0.7.1/debian/compat adagio-0.8.4/debian/compat --- adagio-0.7.1/debian/compat 2021-05-02 00:18:11.000000000 +0000 +++ adagio-0.8.4/debian/compat 2021-05-01 17:48:28.000000000 +0000 @@ -1 +1 @@ -7 \ No newline at end of file +9 \ No newline at end of file diff -Nru adagio-0.7.1/debian/control adagio-0.8.4/debian/control --- adagio-0.7.1/debian/control 2021-05-02 00:18:11.000000000 +0000 +++ adagio-0.8.4/debian/control 2021-05-01 17:48:28.000000000 +0000 @@ -2,13 +2,13 @@ Section: gnu-r Priority: optional Maintainer: cran2deb4ubuntu -Build-Depends: r-base-dev, xvfb, xauth, xfonts-base, r-base-core, - debhelper (>> 4.1.0), cdbs -Standards-Version: 3.9.1 +Build-Depends: r-base-dev, r-cran-lpsolve, xvfb, xauth, xfonts-base, + r-base-core, debhelper (>> 4.1.0), cdbs +Standards-Version: 4.1.4 Package: r-cran-adagio -Architecture: any -Depends: r-base-core, ${shlibs:Depends} +Architecture: all +Depends: r-base-core, r-cran-lpsolve Description: GNU R package "Discrete and Global Optimization Routines" . The R package 'adagio' will provide methods and algorithms for discrete @@ -16,6 +16,6 @@ Nelder-Mead and Hooke-Jeeves minimization, and some (evolutionary) global optimization functions. . - Author: Hans Werner Borchers + Author: Hans W. Borchers [aut, cre] . Maintainer: Hans W. Borchers diff -Nru adagio-0.7.1/debian/copyright adagio-0.8.4/debian/copyright --- adagio-0.7.1/debian/copyright 2021-05-02 00:18:11.000000000 +0000 +++ adagio-0.8.4/debian/copyright 2021-05-01 17:48:28.000000000 +0000 @@ -2,8 +2,8 @@ automatically using cran2deb4ubuntu by cran2deb4ubuntu . -The original GNU R package is Copyright (C) 2018 Hans Werner Borchers -and possibly others. +The original GNU R package is Copyright (C) 2021 Hans W. Borchers [aut, +cre] and possibly others. The original GNU R package is maintained by Hans W. Borchers and was obtained from: diff -Nru adagio-0.7.1/debian/source/format adagio-0.8.4/debian/source/format --- adagio-0.7.1/debian/source/format 1970-01-01 00:00:00.000000000 +0000 +++ adagio-0.8.4/debian/source/format 2021-05-01 17:48:28.000000000 +0000 @@ -0,0 +1 @@ +3.0 (quilt) \ No newline at end of file diff -Nru adagio-0.7.1/DESCRIPTION adagio-0.8.4/DESCRIPTION --- adagio-0.7.1/DESCRIPTION 2018-05-17 21:45:48.000000000 +0000 +++ adagio-0.8.4/DESCRIPTION 2021-04-30 19:00:11.000000000 +0000 @@ -1,24 +1,21 @@ Package: adagio Type: Package Title: Discrete and Global Optimization Routines -Version: 0.7.1 -Date: 2018-05-16 -Author: Hans Werner Borchers +Version: 0.8.4 +Date: 2021-04-30 +Authors@R: person("Hans W.", "Borchers", email = "hwborchers@googlemail.com", + role = c("aut", "cre")) Maintainer: Hans W. Borchers Depends: R (>= 3.1.0) -Imports: graphics, stats +Imports: graphics, stats, lpSolve (>= 5.6.15) Description: The R package 'adagio' will provide methods and algorithms for discrete optimization, e.g. knapsack and subset sum procedures, derivative-free Nelder-Mead and Hooke-Jeeves minimization, and some (evolutionary) global optimization functions. License: GPL (>= 3) -LazyLoad: yes -LazyData: yes +NeedsCompilation: no +Packaged: 2021-04-30 17:28:14 UTC; hwb +Author: Hans W. Borchers [aut, cre] Repository: CRAN -Repository/R-Forge/Project: optimist -Repository/R-Forge/Revision: 454 -Repository/R-Forge/DateTimeStamp: 2018-05-16 19:12:35 -Date/Publication: 2018-05-17 21:45:48 UTC -NeedsCompilation: yes -Packaged: 2018-05-16 19:35:09 UTC; rforge +Date/Publication: 2021-04-30 19:00:11 UTC diff -Nru adagio-0.7.1/man/assignment.Rd adagio-0.8.4/man/assignment.Rd --- adagio-0.7.1/man/assignment.Rd 2012-03-17 11:59:53.000000000 +0000 +++ adagio-0.8.4/man/assignment.Rd 2021-04-25 18:39:30.000000000 +0000 @@ -7,19 +7,21 @@ Linear (sum) assignment problem, or LSAP. } \usage{ -assignment(cmat) +assignment(cmat, dir = "min") } \arguments{ - \item{cmat}{quadratic integer matrix, the cost matrix.} + \item{cmat}{quadratic (numeric) matrix, the cost matrix.} + \item{dir}{direction, can be "min" or "max".} } \details{ - Solves the linear (sum) assignment problem for quadratic matrices with - integer entries. + Solves the linear (sum) assignment problem for quadratic matrices. + Uses the \code{lp.assign} function from the \code{lpSolve} package, + that is it solves LSAP as a mixed integer linear programming problem. } \value{ - List with components \code{perm}, the permutation that defines the minimum - solution, \code{min}, the minimum value, and \code{err}, which is \code{-1} - if an integer overflow occured. + List with components \code{perm}, the permutation that defines the + minimum solution, \code{min}, the minimum value, and \code{err} is + always \code{0}, i.e. not used at the moment. } \references{ Burkard, R., M. Dell'Amico, and S. Martello (2009). Assignment Problems. @@ -28,15 +30,8 @@ Martello, S., and P. Toth (1990). Knapsack Problems: Algorithms and Computer Implementations. John Wiley & Sons, Ltd. } -\author{ - Copyright(c) 1993 A. H. Morris, Jr., Naval Surface Warfare Center, using - Fortran routines written by S. Martello and P. Toth, University of Bologna. - Released for free and general use, now under GPL license, wrapped for R by - Hans W Borchers . -} \note{ - Faster than the Hungarian algorithm, but only applicable to quadratic cost - matrices with integer components. + Slower than the Hungarian algorithm in package \code{clue}. } \seealso{ \code{clue::solve_LSAP} @@ -47,9 +42,10 @@ x <- matrix(sample(1:100), nrow = 10) y <- assignment(x) # show permutation and check minimum sum -y$perm; y$t # 4 5 7 2 6 1 3 8 10 9 -z <- cbind(1:10, y$perm) # 156 -x[z] # 5 4 11 8 20 7 38 15 22 26 +y$perm # 7 6 10 5 8 2 1 4 9 3 +y$min # 173 +z <- cbind(1:10, y$perm) +x[z] # 16 9 49 6 17 14 1 44 10 7 y$min == sum(x[z]) # TRUE \dontrun{ @@ -58,12 +54,13 @@ x <- rt(n, df=3) + 1i * rt(n, df=3) y <- runif(n) + 1i * runif(n) cmat <- round(outer(x, y, FUN = function(x,y) Mod(x - y)), 2) -dmat <- round(100*cmat, 2) -system.time(T1 <- assignment(dmat)) # elapsed: 0.003 -T1$min / 100 # 139.32 +system.time(T1 <- assignment(cmat)) # elapsed: 0.003 +T1$min / 100 # 145.75 +## Hungarian algorithm in package 'clue' library("clue") -system.time(T2 <- solve_LSAP(cmat)) # elapsed: 0.034 -sum(cmat[cbind(1:n, T2)]) # 139.32 +system.time(T2 <- solve_LSAP(cmat)) # elapsed: 0.014 +sum(cmat[cbind(1:n, T2)]) # 145.75 +} } -} \ No newline at end of file +\keyword{ discrete-optimization } diff -Nru adagio-0.7.1/man/binpacking.Rd adagio-0.8.4/man/binpacking.Rd --- adagio-0.7.1/man/binpacking.Rd 1970-01-01 00:00:00.000000000 +0000 +++ adagio-0.8.4/man/binpacking.Rd 2021-04-30 13:51:33.000000000 +0000 @@ -0,0 +1,91 @@ +\name{bpp_approx} +\alias{bpp_approx} +\title{ + Approximate Bin Packing +} +\description{ + Solves the Bin Packing problem approximately. +} +\usage{ +bpp_approx(S, cap, method = c("firstfit", "bestfit", "worstfit")) +} +\arguments{ + \item{S}{vector of weights (or sizes) of items.} + \item{cap}{same capacity for all the bins.} + \item{method}{which approximate method to use.} +} +\details{ + Solves approximately the Bin Packing problem for numeric weights and bins, + all having the same volume. + + Possible methods are "firstfit", "bestfit", and "worstfit". "firstfit" tries + to place each item as early as possible, "bestfit" such that the remaining + space in the bin is as small as possible, and "worstfit" such that the + remaining space is as big as possible. + + Best results are achieved with the "bestfit" method. "firstfit" may be a + reasonable alternative. For smaller and medium-sized data the approximate + results will come quite close to the exact solution, see the examples. + + In general, the results are much better if the items in \code{S} are sorted + decreasingly. If they are not, an immediate warning is issued. +} +\value{ + A list of the following components: + \item{nbins}{minimum number of bins.} + \item{xbins}{index of the bin each item is assigned to.} + \item{sbins}{sum of item sizes in each bin.} + \item{filled}{total volume filled in the bins (as percentage).} +} +\references{ + Silvano Martello. "Bin packing problems". In: 23rd Belgian Mathematical + Optimization Workshop, La-Roche-en-Ardennes 2019. +} +\author{ + Hans W. Borchers +} +\note{ + The Bin Packing problem can be solved as a Linear Program. The formulation + is a bit tricky, and it turned out 'lpSolve' does not solve medium-sized + problems in acceptable time. (Tests with 'Rglpk' will follow.) +} +\seealso{ + Function \code{binpacking} in package 'knapsack' (on R-Forge). +} +\examples{ +## (1) +S <- c(50, 3, 48, 53, 53, 4, 3, 41, 23, 20, 52, 49) +cap <- 100 +bpp_approx(S, cap, method = "bestfit") +## exact -- $nbins 4, filled 99.75 % +## firstfit -- $nbins 6, filled 66.5 % +## bestfit -- $nbins 5, filled 79.8 % +## ! when decreasingly sorted, 'bestfit' with nbins = 4 + +## (2) +S <- c(100,99,89,88,87,75,67,65,65,57,57,49,47,31,27,18,13,9,8,1) +cap <- 100 +bpp_approx(S, cap, method = "firstfit") +# firstfit: 12 bins; exact: 12 bins + +\dontrun{ +## (3) +S <- c(99,99,96,96,92,92,91,88,87,86, + 85,76,74,72,69,67,67,62,61,56, + 52,51,49,46,44,42,40,40,33,33, + 30,30,29,28,28,27,25,24,23,22, + 21,20,17,14,13,11,10, 7, 7, 3) +cap <- 100 +bpp_approx(S, cap) +# exact: 25; firstfit: 25; bestfit: 25 nbins + +## (4) +# 20 no.s in 1..100, capacity 100 +set.seed(7013) +S <- sample(1:100, 20, replace = TRUE) +cap <- 100 +bpp_approx(sort(S, decreasing = TRUE), cap, method = "bestfit") +# exact: 12 bins; firstfit and bestfit: 13; worstfit: 14 bins +} +} +\keyword{ discrete-optimization } diff -Nru adagio-0.7.1/man/cmaes.Rd adagio-0.8.4/man/cmaes.Rd --- adagio-0.7.1/man/cmaes.Rd 2013-09-14 08:56:31.000000000 +0000 +++ adagio-0.8.4/man/cmaes.Rd 2021-04-29 07:42:41.000000000 +0000 @@ -38,26 +38,17 @@ for several minutes; avoid problem dimensions of 30 and more! } \references{ - Hansen, N., and A. Ostermeier (2001). Completely Derandomized Self-Adaptation - in Evolution Strategies. Evolutionary Computation 9(2), pp. 159-195.\cr - URL: https://www.lri.fr/~hansen/cmaartic.pdf - - Hansen, N., S.D. Mueller, P. Koumoutsakos (2003). Reducing the Time - Complexity of the Derandomized Evolution Strategy with Covariance Matrix - Adaptation (CMA-ES). Evolutionary Computation 11(1), pp. 1-18.\cr - URL: https://www.lri.fr/~hansen/evco_11_1_1_0.pdf - Hansen, N. (2011). The CMA Evolution Strategy: A Tutorial.\cr - URL: https://www.lri.fr/~hansen/cmatutorial.pdf + \url{https://arxiv.org/abs/1604.00772} Hansen, N., D.V. Arnold, and A. Auger (2013). Evolution Strategies. - To appear in Janusz Kacprzyk and Witold Pedrycz (Eds.): Handbook of - Computational Intelligence, Springer-Verlag (accepted for publication).\cr - URL: https://www.lri.fr/~hansen/es-overview-2014.pdf + J. Kacprzyk and W. Pedrycz (Eds.). Handbook of Computational Intelligence, + Springer-Verlag, 2015. } \author{ Copyright (c) 2003-2010 Nikolas Hansen for Matlab code PURECMAES; converted to R by Hans W Borchers. + (Hansen's homepage: www.cmap.polytechnique.fr/~nikolaus.hansen/) } \note{ There are other implementations of Hansen's CMAES in package `cmaes' diff -Nru adagio-0.7.1/man/hamiltonian.Rd adagio-0.8.4/man/hamiltonian.Rd --- adagio-0.7.1/man/hamiltonian.Rd 2015-12-06 17:39:08.000000000 +0000 +++ adagio-0.8.4/man/hamiltonian.Rd 2021-04-25 18:41:59.000000000 +0000 @@ -105,4 +105,4 @@ hamiltonian(sq_edges) }} } -\keyword{ graphs } +\keyword{ graph-algorithms } diff -Nru adagio-0.7.1/man/knapsack.Rd adagio-0.8.4/man/knapsack.Rd --- adagio-0.7.1/man/knapsack.Rd 2017-02-06 14:46:16.000000000 +0000 +++ adagio-0.8.4/man/knapsack.Rd 2021-04-26 11:58:12.000000000 +0000 @@ -20,13 +20,13 @@ Maximize \code{sum(x*p)} such that \code{sum(x*w) <= cap}, where \code{x} is a vector with \code{x[i] == 0 or 1}. + + Knapsack procedures can even solve subset sum problems, + see the examples 3 and 3' below. } \value{ A list with components \code{capacity}, \code{profit}, and \code{indices}. } -\note{ - Will be replaced by a compiled version. -} \author{ HwB email: } @@ -54,5 +54,26 @@ cap <- 50 (is <- knapsack(w, p, cap)) # [1] 1 4 , capacity 50 and total profit 107 + +\dontrun{ +## Example 3: subset sum +p <- seq(2, 44, by = 2)^2 +w <- p +is <- knapsack(w, p, 2012) +p[is$indices] # 16 36 64 144 196 256 324 400 576 + +## Example 3': maximize number of items +# w <- seq(2, 44, by = 2)^2 +# p <- numeric(22) + 1 +# is <- knapsack(w, p, 2012) + +## Example 4 from Rosetta Code: +w = c( 9, 13, 153, 50, 15, 68, 27, 39, 23, 52, 11, + 32, 24, 48, 73, 42, 43, 22, 7, 18, 4, 30) +p = c(150, 35, 200, 160, 60, 45, 60, 40, 30, 10, 70, + 30, 15, 10, 40, 70, 75, 80, 20, 12, 50, 10) +cap = 400 +system.time(is <- knapsack(w, p, cap)) # 0.001 sec +} } -\keyword{ optimize } +\keyword{ discrete-optimization } diff -Nru adagio-0.7.1/man/maxempty.Rd adagio-0.8.4/man/maxempty.Rd --- adagio-0.7.1/man/maxempty.Rd 2015-07-25 11:53:55.000000000 +0000 +++ adagio-0.8.4/man/maxempty.Rd 2021-04-25 18:39:55.000000000 +0000 @@ -61,3 +61,4 @@ do.call(rect, as.list(R$rect)) grid()} } +\keyword{ discrete-optimization } diff -Nru adagio-0.7.1/man/maxsub.Rd adagio-0.8.4/man/maxsub.Rd --- adagio-0.7.1/man/maxsub.Rd 2013-09-27 13:36:58.000000000 +0000 +++ adagio-0.8.4/man/maxsub.Rd 2021-04-25 18:40:20.000000000 +0000 @@ -8,7 +8,7 @@ Find a subarray with maximal positive sum. } \usage{ -maxsub(x, inds = TRUE, compiled = TRUE) +maxsub(x, inds = TRUE) maxsub2d(A) } @@ -16,21 +16,20 @@ \item{x}{numeric vector.} \item{A}{numeric matrix} \item{inds}{logical; shall the indices be returned?} - \item{compiled}{logical; shall the compiled version be used?} } \details{ \code{maxsub} finds a contiguous subarray whose sum is maximally positive. -This is sometimes called Kadane's algorithm.\cr -\code{maxsub} will use a compiled and very fast version with a running time -of \code{O(n)} where \code{n} is the length of the input vector \code{x}. +This is sometimes called Kadane's algorithm. +\code{maxsub} will use a very fast version with a running time of +\code{O(n)} where \code{n} is the length of the input vector \code{x}. \code{maxsub2d} finds a (contiguous) submatrix whose sum of elements is maximally positive. The approach taken here is to apply the one-dimensional routine to summed arrays between all rows of \code{A}. This has a run-time of \code{O(n^3)}, though a run-time of \code{O(n^2 log n)} seems possible -see the reference below.\cr -\code{maxsub2d} uses a Fortran workhorse and can solve a 1000-by-1000 matrix -in a few seconds---but beware of biggere ones +see the reference below. +\code{maxsub2d} can solve a 100-by-100 matrix in a few seconds -- +but beware of bigger ones. } \value{ Either just a maximal sum, or a list this sum as component \code{sum} plus @@ -47,7 +46,7 @@ In special cases, the matrix \code{A} may be sparse or (as in the example section) only have one nonzero element in each row and column. Expectation is that there may exists a more efficient (say \code{O(n^2)}) algorithm in - this extreme case. + these special cases. } \author{ HwB @@ -56,7 +55,7 @@ ## Find a maximal sum subvector set.seed(8237) x <- rnorm(1e6) -system.time(res <- maxsub(x, inds = TRUE, compiled = FALSE)) +system.time(res <- maxsub(x, inds = TRUE)) res ## Standard example: Find a maximal sum submatrix @@ -89,4 +88,4 @@ rect(X[msr$inds[3]], Y[msr$inds[1]], X[msr$inds[4]+1], Y[msr$inds[2]+1]) } } -\keyword{ optimize } +\keyword{ discrete-optimization } diff -Nru adagio-0.7.1/man/mknapsack.Rd adagio-0.8.4/man/mknapsack.Rd --- adagio-0.7.1/man/mknapsack.Rd 2012-03-17 11:59:53.000000000 +0000 +++ adagio-0.8.4/man/mknapsack.Rd 2021-04-26 11:58:49.000000000 +0000 @@ -7,54 +7,48 @@ Solves the 0-1 (binary) multiple knapsack problem. } \usage{ -mknapsack(p, w, k, bck = -1) +mknapsack(w, p, cap) } \arguments{ - \item{p}{integer vector of profits.} - \item{w}{integer vector of weights.} - \item{k}{integer vector of capacities of the knapsacks.} - \item{bck}{maximal number of backtrackings allowed; default: -1.} + \item{w}{vector of (positive) weights.} + \item{p}{vector of (positive) profits.} + \item{cap}{vector of capacities of different knapsacks.} } \details{ - Solves the 0-1 multiple knapsack problem for integer profits and weights - + Solves the 0-1 multiple knapsack problem for a set of profits and weights.\cr A multiple 0-1 knapsack problem can be formulated as: \code{ maximize vstar = p(1)*(x(1,1) + ... + x(m,1)) + ... ... + p(n)*(x(1,n) + ... + x(m,n)) subject to - w(1)*x(i,1) + ... + w(n)*x(i,n) <= k(i) for i=1,...,m + w(1)*x(i,1) + ... + w(n)*x(i,n) <= cap(i) for i=1,...,m x(1,j) + ... + x(m,j) <= 1 for j=1,...,n x(i,j) = 0 or 1 for i=1,...,m , j=1,...,n , } - The input problem description must satisfy the following conditions: - \itemize{ - \item \code{vs=-1} if n < 2 or m < 1 - \item \code{vs=-2} if some p(j) , w(j) or k(i) are not positive - \item \code{vs=-3} if a knapsack cannot contain any item - \item \code{vs=-4} if an item cannot fit into any knapsack - \item \code{vs=-5} if knapsack m contains all the items - \item \code{vs=-7} if array k is not correctly sorted - \item \code{vs=-8} [should not happen] - } + The multiple knapsack problem is reformulated as a linear program and + solved with the help of package \code{lpSolve}. + + This function can be used for the single knapsack problem as well, + but the 'dynamic programming' version in the \code{knapsack} function + is faster (but: allows only integer values). + + The solution found is most often not unique and may not be the most compact + one. In the future, we will attempt to 'compactify' through backtracking. + The number of backtracks will be returned in list element \code{bs}. } \value{ - A list with compomnents, \code{ksack} the knapsack numbers the items are + A list with components, \code{ksack} the knapsack numbers the items are assigned to, \code{value} the total value/profit of the solution found, and \code{bs} the number of backtracks used. } \note{ - With some care, this function can be used for the bounded and unbounded - single knapsack problem as well. -} -\author{ - The Fortran source code is adapted from the free NSCW Library of - Mathematical Subroutines. + Contrary to earlier versions, the sequence of profits and weights has been + interchanged: first the weights, then profits. - The wrapping code has been written by yours package maintainer,\cr - HwB email: + The compiled version was transferred to the \code{knapsack} package on + R-Forge (see project 'optimist'). } \references{ Kellerer, H., U. Pferschy, and D. Pisinger (2004). Knapsack Problems. @@ -68,35 +62,61 @@ } \examples{ ## Example 1: single knapsack -p <- c(15, 100, 90, 60, 40, 15, 10, 1) w <- c( 2, 20, 20, 30, 40, 30, 60, 10) +p <- c(15, 100, 90, 60, 40, 15, 10, 1) cap <- 102 -(is <- mknapsack(p, w, cap)) +(is <- mknapsack(w, p, cap)) which(is$ksack == 1) # [1] 1 2 3 4 6 , capacity 102 and total profit 280 ## Example 2: multiple knapsack -p <- c(110, 150, 70, 80, 30, 5) w <- c( 40, 60, 30, 40, 20, 5) -k <- c(65, 85) -is <- mknapsack(p, w, k) -# kps 1: 2,6; kps 2: 1,4; value: 345; backtracks: 14 +p <- c(110, 150, 70, 80, 30, 5) +cap <- c(85, 65) +is <- mknapsack(w, p, cap) +# kps 1: 1,4; kps 2: 2,6; value: 345 ## Example 3: multiple knapsack p <- c(78, 35, 89, 36, 94, 75, 74, 79, 80, 16) w <- c(18, 9, 23, 20, 59, 61, 70, 75, 76, 30) -k <- c(103, 156) -is <- mknapsack(p, w, k) -# kps 1: 1,3,6; kps 2: 4,5,9; value: 452; backtracks: 4 - -## Example 4: subset sum -p <- seq(2, 44, by = 2)^2 -w <- p -is <- mknapsack(p, w, 2012) -sum((2 * which(is$ksack == 1))^2) - -## Example 5: maximize number of items -w <- seq(2, 44, by = 2)^2 -p <- numeric(22) + 1 -is <- mknapsack(p, w, 2012) +cap <- c(103, 156) +is <- mknapsack(w, p, cap) +# kps 1: 3,4,5; kps 2: 1,6,9; value: 452 + +\dontrun{ +# How to Cut Your Planks with R +# R-bloggers, Rasmus Baath, 2016-06-12 +# +# This is application of multiple knapsacks to cutting planks into pieces. + +planks_we_have <- c(120, 137, 220, 420, 480) +planks_we_want <- c(19, 19, 19, 19, 79, 79, 79, 103, 103, + 103, 135, 135, 135, 135, 160) +s <- mknapsack(planks_we_want, planks_we_want + 1, planks_we_have) +s$ksack +## [1] 5 5 5 5 3 5 5 4 1 5 4 5 3 2 4 + +# Solution w/o backtracking +# bin 1 : 103 | Rest: 17 +# bin 2 : 135 | Rest: 2 +# bin 3 : 79 + 135 | Rest: 6 +# bin 4 : 103 + 135 + 160 | Rest: 22 +# bin 5 : 4*19 + 2*79 + 103 + 135 | Rest: 8 +# +# Solution with reversing the bins (bigger ones first) +# bin 1 : 103 | Rest: 4 +# bin 2 : 2*19 + 79 | Rest: 20 +# bin 3 : 79 + 135 | Rest: 6 +# bin 4 : 2*19 + 79 + 135 + 160 | Rest: 8 +# bin 5 : 2*103 + 2*135 | Rest: 17 +# +# Solution with backtracking (compactification) +# sol = c(1, 4, 4, 1, 1, 3, 4, 5, 5, 5, 5, 4, 2, 3, 4) +# bin 1 : 2*19 + 79 | Rest: 3 +# bin 2 : 135 | Rest: 2 +# bin 3 : 79 + 135 | Rest: 6 +# bin 4 : 2*19 + 79 + 135 + 160 | Rest: 8 +# bin 5 : 3*103 + 135 | Rest: 36 +} } +\keyword{ discrete-optimization } diff -Nru adagio-0.7.1/man/occurs.Rd adagio-0.8.4/man/occurs.Rd --- adagio-0.7.1/man/occurs.Rd 1970-01-01 00:00:00.000000000 +0000 +++ adagio-0.8.4/man/occurs.Rd 2021-04-28 09:35:45.000000000 +0000 @@ -0,0 +1,42 @@ +\name{occurs} +\alias{occurs} +\title{ + Finding Subsequences +} +\description{ + Find subsequences of (integer) sequences. +} +\usage{ + occurs(subseq, series) +} +\arguments{ + \item{subseq}{vector of integers.} + \item{series}{vector of integers.} +} +\details{ + If \code{m} and \code{n} are the lengths of \code{s} and \code{S} resp., + \code{occurs(s, S)} determines all positions \code{i} such that + \code{s == S[i, ..., i+m-1]}. + + The code is vectorized and relatively fast. It is intended to complement + this with an implementation of Rabin-Karp, and possibly Knuth-Morris-Pratt + and Boyer-Moore algorithms. +} +\value{ + Returns a vector of indices. +} +\examples{ +## Examples +patrn <- c(1,2,3,4) +exmpl <- c(3,3,4,2,3,1,2,3,4,8,8,23,1,2,3,4,4,34,4,3,2,1,1,2,3,4) +occurs(patrn, exmpl) +## [1] 6 13 23 + +\dontrun{ +set.seed(2437) +p = sample(1:20, 1000000, replace=TRUE) +system.time(i <- occurs(c(1,2,3,4,5), p)) #=> [1] 799536 +## user system elapsed +## 0.017 0.000 0.017 [sec] +} +} diff -Nru adagio-0.7.1/man/setcover.Rd adagio-0.8.4/man/setcover.Rd --- adagio-0.7.1/man/setcover.Rd 1970-01-01 00:00:00.000000000 +0000 +++ adagio-0.8.4/man/setcover.Rd 2021-04-25 18:41:25.000000000 +0000 @@ -0,0 +1,53 @@ +\name{setcover} +\alias{setcover} +\title{ + Set cover problem +} +\description{ + Solves the Set Cover problem as an integer linear program. +} +\usage{ + setcover(Sets, weights) +} +\arguments{ + \item{Sets}{matrix of 0s and 1s, each line defining a subset.} + \item{weights}{numerical weights for each subset.} +} +\details{ + The Set Cover problems attempts to find in subsets (of a 'universe') + a minimal set of subsets that still covers the whole set. + + Each line of the matrix \code{Sets} defines a characteristic function of + a subset. It is required that each element of the universe is contained + in at least one of these subsets. + + The problem is treated as an Integer Linear Program (ILP) and solved + with the \code{lp} solver in \code{lpSolve}. +} +\value{ + Returns a list with components \code{sets}, giving the indices of subsets, + and \code{objective}, the sum of weights of subsets present in the solution. +} +\references{ + See the Wikipedia article on the "set cover problem". +} +\seealso{ + \code{\link{knapsack}} +} +\examples{ +# Define 12 subsets of universe {1, ..., 10}. +set.seed(7*11*13) +A <- matrix(sample(c(0,1), prob = c(0.8,0.2), size = 120, replace =TRUE), + nrow = 12, ncol = 10) +sol <- setcover(Sets = A, weights = rep(1, 12)) +sol +## $sets +## [1] 1 2 9 12 +## $no.sets +##[1] 4 + +# all universe elements are covered: +colSums(A[sol$sets, ]) +## [1] 1 1 2 1 1 1 2 1 1 2 +} +\keyword{ discrete-optimization } diff -Nru adagio-0.7.1/man/subsetsum.Rd adagio-0.8.4/man/subsetsum.Rd --- adagio-0.7.1/man/subsetsum.Rd 2013-09-27 13:36:58.000000000 +0000 +++ adagio-0.8.4/man/subsetsum.Rd 2021-04-28 18:18:11.000000000 +0000 @@ -1,5 +1,6 @@ \name{subsetsum} \alias{subsetsum} +\alias{sss_test} \title{ Subset Sum Problem } @@ -8,23 +9,36 @@ } \usage{ subsetsum(S, t, method = "greedy") + +sss_test(S, t) } \arguments{ \item{S}{vector of positive integers.} - \item{t}{target value.} + \item{t}{target value, bigger than all items in \code{S}.} \item{method}{can be ``greedy'' or ``dynamic'', where ``dynamic'' stands for the dynamic programming approach.} } \details{ - Searching for a set of elements in \code{S} that sum up to \code{t} by - continuously adding more elements of \code{S}. + \code{subsetsum} is searching for a set of elements in \code{S} that + sum up to \code{t} by continuously adding more elements of \code{S}. + + It is not required that \code{S} is decreasingly sorted. But for reasons + of efficiency and smaller execution times it is urgently recommended to + sort the item set in decreasing order. See the examples to find out how + to handle your data. The first components will be preferred, i.e., if \code{S} is decreasing, the sum with larger elements will be found, if increasing, the sum with - smaller elements. + smaller elements. Because of timing considerations, the default is to + sort decreasingly before processing. The dynamic method may be faster for large sets, but will also require - much more memory if the target value is large. + much more memory if the target value is large. + + \code{sss_test} will find the biggest number below or equal to \code{t} + that can be expressed as a sum of items in \code{S}. It will not return + any indices. It can be quite fast, though it preprocesses the set \code{S} + to be sorted decreasingly, too. } \value{ List with the target value, if reached, and vector of indices of elements @@ -32,9 +46,14 @@ If no solution is found, the dynamic method will return indices for the largest value below the target, the greedy method witll return NULL. + + \code{sss_test} will simply return maximum sum value found. } \note{ - Will be replaced by a compiled version. + A compiled version -- and much faster, in Fortran -- can be found in + package 'knapsack' (R-Forge, project 'optimist') as \code{subsetsum}. + A recursive version, returning *all* solutions, is much too slow in R, + but is possible in Julia and can be asked from the author. } \author{ HwB email: @@ -47,6 +66,17 @@ \code{\link{maxsub}} } \examples{ +t <- 5842 +S <- c(267, 493, 869, 961, 1000, 1153, 1246, 1598, 1766, 1922) + +# S is not decreasingly sorted, so ... +o <- order(S, decreasing = TRUE) +So <- S[o] # So is decreasingly sorted + +sol <- subsetsum(So, t) # $inds: 2 4 6 7 8 w.r.t. So +is <- o[sol$inds] # is: 9 7 5 4 3 w.r.t. S +sum(S[is]) # 5842 + \dontrun{ amount <- 4748652 products <- @@ -66,13 +96,20 @@ # now find one solution system.time(is <- subsetsum(prods, amount)) # user system elapsed -# 0.320 0.032 0.359 +# 0.030 0.000 0.029 prods[is] # [1] 70740 70740 71190 76950 77382 80000 83800 # [8] 90900 456000 533600 829350 1110000 1198000 sum(prods[is]) == amount -# [1] TRUE} +# [1] TRUE + +# Timings: +# unsorted decr.sorted +# "greedy" 22.930 0.030 (therefore the default settings) +# "dynamic" 2.515 0.860 (overhead for smaller sets) +# sss_test 8.450 0.040 (no indices returned) +} } -\keyword{ optimize } +\keyword{ discrete-optimization } diff -Nru adagio-0.7.1/man/testfunctions.Rd adagio-0.8.4/man/testfunctions.Rd --- adagio-0.7.1/man/testfunctions.Rd 2015-07-25 11:53:55.000000000 +0000 +++ adagio-0.8.4/man/testfunctions.Rd 2021-04-25 11:35:45.000000000 +0000 @@ -59,6 +59,8 @@ Solution: \tab xi = 1, i = 1:n \cr } +\bold{Nesterov1} -- Simlar to \code{Nesterov}, except the terms added are taken with absolute value, which makes this function nonsmooth and painful for gradient-based optimization routines; no gradient prpovided. + \bold{Rastrigin} -- Rastrigin's function is a famous, non-convex example from 1989 for global optimization. It is a typical example of a multimodal function with many local minima: \deqn{10 n + \sum_1^n (x_i^2 - 10 \cos(2 \pi x_i))} \tabular{ll}{ diff -Nru adagio-0.7.1/MD5 adagio-0.8.4/MD5 --- adagio-0.7.1/MD5 2018-05-17 21:45:48.000000000 +0000 +++ adagio-0.8.4/MD5 2021-04-30 19:00:11.000000000 +0000 @@ -1,7 +1,8 @@ -ad6be58fd8522bf69a24f18fa8c822e4 *DESCRIPTION -bb5996d7ab51d07ab67591d28923ce49 *NAMESPACE -39f93ff8cabf4babe6de22c718414dc9 *NEWS -0e70f23f927819f7c8144c44b21cc29e *R/assignment.R +e3f48e5c5924c2c797ddbc77eb2f40f5 *DESCRIPTION +eee5fab879bf100e1fb36d9bf92f150a *NAMESPACE +80ea66caabe9ffd3c837e32fb774fd48 *NEWS +8be7862cf0a4b8b2a8210bb1498f8838 *R/assignment.R +6be00c860625858214957318fcd892c4 *R/binpacking.R 352ae4d29230c581f3caed401a9340c8 *R/cmaes.R ea81f864b25f9cf6deec45909cb88181 *R/fminviz.R 7160b272cd05915b3bbd40a27f59358c *R/gradients.R @@ -10,35 +11,34 @@ a2019b193d89d79701abba21bf8a0c9c *R/knapsack.R f842fe7ab721077d49dd00889d6039df *R/maxempty.R 93f38426b2771a3253768ac46b5beb5c *R/maxquad.R -6ee743bf5f43dac47d951525bc474246 *R/maxsub.R -778f6e936ddfdc4bebcdcea8572752e3 *R/mknapsack.R +1657df87f9cbe0e99f2eab1fdc16af11 *R/maxsub.R +cb70e35a8b6dd09ae8c4662952fbe9fb *R/mknapsack.R f0cebd647f616d6449eac13a9fec986c *R/nelmin.R +56a0e87c23455f44df51d50a78ae212c *R/occurs.R +cb6167ffe5f479392d7abda340c6bf0b *R/setcover.R a51a39815609ab33e9ec0a8cc5103d4e *R/simpleDE.R 2ed28c83c6c3b050c0e09c47df30279e *R/simpleEA.R -eb7a19c4093585e2d79a6a5175f11e7e *R/subsetsum.R +63517b4d342dc0a4f07a9d348f073a6e *R/subsetsum.R 2decf9cde3b507637a7e9ec2eb3f47b0 *R/testfunctions.R 305c9fe88cff3d3506799f8ccdeb22a2 *R/transfinite.R 53e27f9adf2a87445a24610792a26e43 *R/trefethen.R -7bb3d0d7477b53b97393841782c1cfba *man/assignment.Rd -d02552de1d18a85e8030419393af743e *man/cmaes.Rd +f0694675110dfeb527ad7966b560313e *man/assignment.Rd +a8cc39dda6ab7cead5d2945d261bb42c *man/binpacking.Rd +9c91900c41d7c8ed136e5519e884219b *man/cmaes.Rd f02eb8a97850184dbf3d30221d983182 *man/fminviz.Rd -8ce2a49c7384dce354231245ab991547 *man/hamiltonian.Rd +980829d65f2ccb46e6ce2ef402baa371 *man/hamiltonian.Rd 5094b24cc63bcc9c976366583694f9b6 *man/hookejeeves.Rd -f29ed0e0c83f616807d00496e9a23215 *man/knapsack.Rd -4d0edcf5390e437082ddafe30080bba6 *man/maxempty.Rd +b8cd6059a41f362f4c541f4eca82c47b *man/knapsack.Rd +e7f1766b1cc553f231011650596d73ad *man/maxempty.Rd a468c6be9714a530ce3d448766639460 *man/maxquad.Rd -7956c59835fe68cc9f4b16dc082b066b *man/maxsub.Rd -5ab62c06526e327dfefdb2a39401f28a *man/mknapsack.Rd +3d29b5a9f3d3625144c6a0f13e4da777 *man/maxsub.Rd +73daaf5f534a86e2337af6429756ceb6 *man/mknapsack.Rd 32a224e61b1b319d8e926fb06d0d7c2f *man/nelmin.Rd +1547b638eadec6ddf2e2749c4f22e0b7 *man/occurs.Rd +154d5038b93548d8800c8558aaf4c6fd *man/setcover.Rd 4d7972348c182c2e8b52b33c4b09b60b *man/simpleDE.Rd 258c7990b092d8fe0b8a109ecbad792a *man/simpleEA.Rd -a9d918dc27fe8bbfb40a75f8fb33c387 *man/subsetsum.Rd -6be402e6c64b8bb92e960ca594cb1dfa *man/testfunctions.Rd +d5898d5097abd338cc81e7dd8a8bd9f1 *man/subsetsum.Rd +11c603d58d230f456c61642f6ffb3cca *man/testfunctions.Rd df61a392882b1cc7a34124b28406e25b *man/transfinite.Rd 453eac04967f8a02c455ff2461d90cd4 *man/trefethen.Rd -cc8b078de83ee0886536fbfd68cfb5b8 *src/assgn.f -780d3288cfda5d030fe5669630fa20ec *src/maxsubf.f -52693d51172a3f2ac4cef0f35e8b31af *src/mkp.f -21175f956a07bbee3d760385159a989b *src/nscw_util.f -8eb8046e7d8c44fe15039ca5e58c5c81 *src/ratfor/maxempty_rat.for -418188af7e38bead4f17b10212a3fab5 *src/ratfor/maxsub_rat.for diff -Nru adagio-0.7.1/NAMESPACE adagio-0.8.4/NAMESPACE --- adagio-0.7.1/NAMESPACE 2017-02-06 14:46:16.000000000 +0000 +++ adagio-0.8.4/NAMESPACE 2021-04-30 13:55:59.000000000 +0000 @@ -2,11 +2,14 @@ ## NAMESPACE ## -useDynLib(adagio) +# Fortran codes +# useDynLib(adagio) importFrom("graphics", "grid", "lines", "plot", "points") importFrom("stats", "rnorm", "runif") +importFrom(lpSolve, lp, lp.assign) + # Public routines export(fnRosenbrock, grRosenbrock, fnRastrigin, grRastrigin, @@ -15,10 +18,11 @@ export(fnTrefethen, fnWagon) export(maxquad) -export(maxsub, maxsub2d) -export(maxempty) -export(subsetsum, knapsack) -export(mknapsack, assignment) +export(occurs) +export(maxsub, maxsub2d, maxempty) +export(assignment, subsetsum, sss_test, + setcover, knapsack, mknapsack, + bpp_approx) export(neldermead, neldermeadb) export(hookejeeves) diff -Nru adagio-0.7.1/NEWS adagio-0.8.4/NEWS --- adagio-0.7.1/NEWS 2018-05-16 17:35:57.000000000 +0000 +++ adagio-0.8.4/NEWS 2021-04-30 13:57:25.000000000 +0000 @@ -3,6 +3,41 @@ ------------------------------------------------------------------------------ +Changes in Version 0.8.4 (2021-04-30) + + o bpp_approx() solves the Bin Packing Problem approximately. + +Changes in Version 0.8.3 (2021-04-28) + + o Small corrections in subsetsum(); added sss_test() -- no indices. + Changed default values of arguments for methods and sorting. + o occurs() finds index positions of subsequences in integer sequences. + o Removed inactive URL links (Hansen's homepage changed). + +Changes in Version 0.8.2 (2021-04-24) + + o setcover() solves the Set Cover problem in discrete optimization. + o Added a one-line description for testfunction 'fnNesterov1'. + o Corrected help pages for 'knapsack()' and 'mknapsack()'. + +Changes in Version 0.8.1 (2021-04-22) + + o Removed Fortran code for 'maxsub()' and 'maxsub2d()'. + o 'assignment()' uses 'lp.assign' from the 'lpSolve' package; + now numeric cost matrices are allowed, not only integer ones. + o 'mknapsack()': Fortran code replaced by a linear program; + also interchanged the weight and profit arguments. + +Changes in Version 0.8.0 (2021-04-21) + + o Restart 'adagio' by avoiding Fortran code now and forever. + o Imports 'lpSolve' for treating mixed-integer linear programs. + + +Changes in Version 0.7.2 (2019-12-12) + + o Corrected maxsub2d() and Fortran, fix provided by Yogesh Bansal. + Changes in Version 0.7.1 (2018-05-16) o Removed 'adagio-package.Rd' on request of K. Hornick, CRAN. diff -Nru adagio-0.7.1/R/assignment.R adagio-0.8.4/R/assignment.R --- adagio-0.7.1/R/assignment.R 2012-03-17 11:59:53.000000000 +0000 +++ adagio-0.8.4/R/assignment.R 2021-04-24 13:08:03.000000000 +0000 @@ -3,28 +3,47 @@ ## -assignment <- function(cmat) { +assignment <- function(cmat, dir = "min") { stopifnot(is.numeric(cmat), is.matrix(cmat)) + if (! dir %in% c("min", "max")) + stop("Argument 'dir' must be 'min' or 'max'.") n <- nrow(cmat); m <- ncol(cmat) - if (n != m || n <= 1) + if (n != m || n <= 1) stop("Argument 'cmat' must be a quadratic, numeric matrix.") - if (any(floor(cmat) != ceiling(cmat))) { - warning("Matrix 'cmat' not integer; will take floor of it.") - cmat <- floor(cmat) - } - - a <- cbind(t(cmat), rep(0, n)) - b <- rep(0, n) - t <- 0; ierr <- 0 - # dummies - iwk <- vector("integer", 7*n+2) + + smat <- lp.assign(cmat, direction = dir)$solution + inds <- which(round(smat) == 1, arr.ind = TRUE) + + ordr <- order(inds[, 1]) + perm <- inds[, 2][ordr] + mini <- sum(cmat * smat) + + return(list(perm = perm, min = mini, err = 0)) +} - S <- .Fortran("assgn", as.integer(n), as.integer(a), - b = as.integer(b), t = as.integer(t), - as.integer(iwk), err = as.integer(ierr), - PACKAGE = "adagio") - if (S$err != 0) - warning("Integer overflow happened; result may not be correct.") - return(list(perm = S$b, min = S$t, err = S$err)) -} +# assignment <- function(cmat) { +# stopifnot(is.numeric(cmat), is.matrix(cmat)) +# n <- nrow(cmat); m <- ncol(cmat) +# if (n != m || n <= 1) +# stop("Argument 'cmat' must be a quadratic, numeric matrix.") +# if (any(floor(cmat) != ceiling(cmat))) { +# warning("Matrix 'cmat' not integer; will take floor of it.") +# cmat <- floor(cmat) +# } +# +# a <- cbind(t(cmat), rep(0, n)) +# b <- rep(0, n) +# t <- 0; ierr <- 0 +# # dummies +# iwk <- vector("integer", 7*n+2) +# +# S <- .Fortran("assgn", as.integer(n), as.integer(a), +# b = as.integer(b), t = as.integer(t), +# as.integer(iwk), err = as.integer(ierr), +# PACKAGE = "adagio") +# +# if (S$err != 0) +# warning("Integer overflow happened; result may not be correct.") +# return(list(perm = S$b, min = S$t, err = S$err)) +# } diff -Nru adagio-0.7.1/R/binpacking.R adagio-0.8.4/R/binpacking.R --- adagio-0.7.1/R/binpacking.R 1970-01-01 00:00:00.000000000 +0000 +++ adagio-0.8.4/R/binpacking.R 2021-04-30 13:49:04.000000000 +0000 @@ -0,0 +1,55 @@ +## +## b i n p a c k i n g . r Approximate Bin Packing +## + + +## ===================================================================== +## first-fit and best-fit heuristic for bin packing + +# Approximate solution to the bin packing problem +# Method: first_fit, input item size, and bin size +# Returns: no.of bins, assigned bin, bin allocation +# First-Fit, Best-Fit, decreasing = TRUE +bpp_approx <- function(S, cap, + method = c("firstfit", "bestfit", "worstfit")) { + stopifnot(is.numeric(S), is.numeric(cap)) + method <- match.arg(method) + if (length(cap) != 1 || any(S > cap)) + stop("Argument 'cap' must be a scalar and all S <= cap.") + if (is.unsorted(rev(S))) { + warning("Argument 'S' of item weights is not decreasingly sorted.", + call. = FALSE, immediate. = TRUE) + } + + n <- length(S) + b <- rep(cap, n) + p <- rep(0, n) + + if (method == "firstfit") { + for (i in 1:n) { + j <- which(b >= S[i])[1] + p[i] <- j + b[j] <- b[j] - S[i] + } + } else if (method == "bestfit") { + for (i in 1:n) { + inds <- which(b >= S[i]) + j <- inds[which.min(b[inds])] + p[i] <- j + b[j] <- b[j] - S[i] + } + } else if (method == "worstfit") { + for (i in 1:n) { + inds <- which(b >= S[i] & b < cap) + if (length(inds) == 0) inds <- which(b >= S[i]) + j <- inds[which.max(b[inds])] + p[i] <- j + b[j] <- b[j] - S[i] + } + } else + stop("Unknown method; may be implemented in the course of time.") + + m <- max(p) + f <- sum(S) / m / cap + return(list(nbins = m, xbins = p, sbins = cap - b[1:m], filled = f)) +} diff -Nru adagio-0.7.1/R/maxsub.R adagio-0.8.4/R/maxsub.R --- adagio-0.7.1/R/maxsub.R 2012-01-25 18:56:39.000000000 +0000 +++ adagio-0.8.4/R/maxsub.R 2021-04-24 09:49:57.000000000 +0000 @@ -3,85 +3,123 @@ ## m a x s u b 2 d Maximal Sum Subrectangle Problem ## - -maxsub <- function(x, inds = TRUE, compiled = TRUE) { +maxsub <- function(x, inds = TRUE) { if (!is.numeric(x)) stop("Argument 'x' must be a numeric vector.") n <- length(x) i1 <- 0; i2 <- 0 s <- 0.0 - - if (compiled) { - S <- .Fortran("maxsubf", x = as.numeric(x), n = as.integer(n), - s = as.numeric(s), - i1 = as.integer(i1), i2 = as.integer(i2), - PACKAGE = "adagio") - if (inds) - return(list(sum = S$s, inds = c(S$i1, S$i2))) - else - return(S$s) - + + if (!inds) { + m1 <- m2 <- 0.0 + for (i in 1:n) { + m2 <- max(m2 + x[i], 0.0) + m1 <- max(m1, m2) + } + return(m1) } else { - if (!inds) { - m1 <- m2 <- 0.0 - for (i in 1:n) { - m2 <- max(m2 + x[i], 0.0) - m1 <- max(m1, m2) - } - return(m1) - } else { - m1 <- m2 <- 0 - p1 <- p2 <- 0 - q1 <- q2 <- 1 - - for (i in 1:n) { - if (m2 > -x[i]) { - m2 <- m2 + x[i] - q2 <- i - if (m2 > m1) { - m1 <- m2 - p1 <- q1; p2 <- q2 - } - } else { - m2 <- 0 - q1 <- q2 <- i+1 + m1 <- m2 <- 0 + p1 <- p2 <- 0 + q1 <- q2 <- 1 + for (i in 1:n) { + if (m2 > -x[i]) { + m2 <- m2 + x[i] + q2 <- i + if (m2 > m1) { + m1 <- m2 + p1 <- q1; p2 <- q2 } + } else { + m2 <- 0 + q1 <- q2 <- i+1 } - return(list(sum = m1, inds = c(p1, p2))) } + return(list(sum = m1, inds = c(p1, p2))) } } - maxsub2d <- function(A) { stopifnot(is.numeric(A), is.matrix(A)) n <- nrow(A); m <- ncol(A) - + if (all(A <= 0)) stop("At least on element of matrix 'A' must be positive.") if (all(A >= 0)) return(list(sum = sum(A), inds = c(1, n, 1, m), submat = A)) - - mi <- vector("integer", 4) + S <- matrix(0, nrow = n+1, ncol = m) - aa <- numeric(m) - b <- 0.0 - - fm <- 0.0 - R <- .Fortran("maxsub2f", as.numeric(A), as.numeric(S), - as.integer(n), as.integer(m), - fmax = as.numeric(fm), mind = as.integer(mi), - as.numeric(aa), as.numeric(b), - PACKAGE = "adagio") - - fm <- R$fmax - mi <- R$mind - - return(list(sum = fm, inds = mi, + S[1, ] <- 0 + for (i in 2:(n+1)) + S[i, ] <- S[i-1, ] + cumsum(A[i-1, ]) + + mm <- 0 + mi <- c(0, 0, 0, 0) + + for (i in 2:(n+1)) { + for (j in i:(n+1)) { + a <- numeric(n) + a[1] <- S[j, 1] - S[i-1, 1] + for (k in 2:m) ### HwB replaced n with m -- 2019-12-12 + a[k] <- S[j, k] - S[(i-1), k] - sum(a[1:(k-1)]) + ms <- maxsub(a) + if (ms$sum > mm) { + mm <- ms$sum + mi <- c(i-1, j-1, ms$inds[1], ms$inds[2]) + } + } + } + return(list(sum = mm, inds = mi, submat = A[mi[1]:mi[2], mi[3]:mi[4], drop = FALSE])) } - +# maxsub <- function(x, inds = TRUE, compiled = TRUE) { +# if (!is.numeric(x)) +# stop("Argument 'x' must be a numeric vector.") +# n <- length(x) +# i1 <- 0; i2 <- 0 +# s <- 0.0 +# +# if (compiled) { +# S <- .Fortran("maxsubf", x = as.numeric(x), n = as.integer(n), +# s = as.numeric(s), +# i1 = as.integer(i1), i2 = as.integer(i2), +# PACKAGE = "adagio") +# if (inds) +# return(list(sum = S$s, inds = c(S$i1, S$i2))) +# else +# return(S$s) +# +# } else { +# if (!inds) { +# m1 <- m2 <- 0.0 +# for (i in 1:n) { +# m2 <- max(m2 + x[i], 0.0) +# m1 <- max(m1, m2) +# } +# return(m1) +# } else { +# m1 <- m2 <- 0 +# p1 <- p2 <- 0 +# q1 <- q2 <- 1 +# +# for (i in 1:n) { +# if (m2 > -x[i]) { +# m2 <- m2 + x[i] +# q2 <- i +# if (m2 > m1) { +# m1 <- m2 +# p1 <- q1; p2 <- q2 +# } +# } else { +# m2 <- 0 +# q1 <- q2 <- i+1 +# } +# } +# return(list(sum = m1, inds = c(p1, p2))) +# } +# } +# } +# # maxsub2d <- function(A) { # stopifnot(is.numeric(A), is.matrix(A)) # n <- nrow(A); m <- ncol(A) @@ -91,27 +129,21 @@ # if (all(A >= 0)) # return(list(sum = sum(A), inds = c(1, n, 1, m), submat = A)) # +# mi <- vector("integer", 4) # S <- matrix(0, nrow = n+1, ncol = m) -# S[1, ] <- 0 -# for (i in 2:(n+1)) -# S[i, ] <- S[i-1, ] + cumsum(A[i-1, ]) -# -# mm <- 0 -# mi <- c(0, 0, 0, 0) -# -# for (i in 2:(n+1)) { -# for (j in i:(n+1)) { -# a <- numeric(n) -# a[1] <- S[j, 1] - S[i-1, 1] -# for (k in 2:n) -# a[k] <- S[j, k] - S[(i-1), k] - sum(a[1:(k-1)]) -# ms <- maxsub(a) -# if (ms$sum > mm) { -# mm <- ms$sum -# mi <- c(i-1, j-1, ms$inds[1], ms$inds[2]) -# } -# } -# } -# return(list(sum = mm, inds = mi, -# submat = A[mi[1]:mi[2], mi[3]:mi[4], drop = FALSE])) +# aa <- numeric(m) +# b <- 0.0 +# +# fm <- 0.0 +# R <- .Fortran("maxsub2f", as.numeric(A), as.numeric(S), +# as.integer(n), as.integer(m), +# fmax = as.numeric(fm), mind = as.integer(mi), +# as.numeric(aa), as.numeric(b), +# PACKAGE = "adagio") +# +# fm <- R$fmax +# mi <- R$mind +# +# return(list(sum = fm, inds = mi, +# submat = A[mi[1]:mi[2], mi[3]:mi[4], drop = FALSE])) # } diff -Nru adagio-0.7.1/R/mknapsack.R adagio-0.8.4/R/mknapsack.R --- adagio-0.7.1/R/mknapsack.R 2012-03-11 23:04:34.000000000 +0000 +++ adagio-0.8.4/R/mknapsack.R 2021-04-25 08:25:20.000000000 +0000 @@ -1,39 +1,61 @@ -mknapsack <- function(p, w, k, bck = -1) { - stopifnot(is.numeric(p), is.numeric(w), is.numeric(k)) - if (any(w <= 0)) - stop("'weights' must be a vector of positive numbers.") - if (any(p <= 0)) - stop("'profits' must be a vector of positive numbers.") - - if (any(floor(p) != ceiling(p)) || - any(floor(w) != ceiling(w)) || - any(floor(k) != ceiling(k)) || - any(p >= 2^31) || any(w >= 2^31) || any(k >= 2^31)) - stop("All inputs must be positive integers < 2^31 !") +# w <- c( 40, 60, 30, 40, 20, 5) +# p <- c(110, 150, 70, 80, 30, 5) +# cap <- 90 # c(65, 85) + +mknapsack <- function(w, p, cap) { + stopifnot(is.numeric(p), is.numeric(w), is.numeric(cap)) + if (any(w <= 0) || any(p <= 0)) + stop("Profits 'w' and weights 'p' must be vectors of positive numbers.") - n <- length(p); m <- length(k) + m <- length(cap) + n <- length(p) # == length(w) if (length(w) != n) - stop("Profit 'p' and weight 'w' must be vectors of equal length.") - # bck <- -1, i.e., request exact solution, - # else maximal number of allowed backtrackings - - xstar <- vector("integer", n) # no. of knapsack item 'i' is assigned to - vstar <- 0 # value of optimal solution - - # dummy variables - num <- 5*m + 14*n + 4*m*n + 3 - wk <- numeric(n) - iwk <- vector("integer", num) - - S <- .Fortran("mkp", as.integer(n), as.integer(m), - as.integer(p), as.integer(w), as.integer(k), - bs = as.integer(bck), - xs = as.integer(xstar), vs = as.integer(vstar), - as.numeric(wk), as.integer(iwk), as.integer(num), - PACKAGE = 'adagio') + stop("Profits 'p' and weights 'w' must be vectors of equal length.") - if (S$vs < 0) - warning("Error condition raised: check input data ...!") + if (n == 1) { + inds <- which(cap >= w) + if (length(inds) == 0) { + return(list(ksack = 0, val = 0)) + } else { + ind <- inds[1] + return(list(ksack = ind, val = p, bs = 0)) + } + } + + if (m == 1) { + sol <- lp(direction = "max", + objective.in = p, + const.mat = matrix(w, nrow = 1), + const.dir = "<=", + const.rhs = cap, + all.bin = TRUE) + return(list(ksack = sol$solution, val = sol$objval, bs = 0)) + } + + obj <- rep(p, m) + cm1 <- matrix(0, nrow = m, ncol = m * n) + for (k in 1:m) { + cm1[k, ((k-1)*n+1):(k*n)] <- w + } + cm2 <- diag(1, n) + for (k in 2:m) { + cm2 <- cbind(cm2, diag(1, n)) + } + cm <- rbind(cm1, cm2) + + sol <- lp(direction = "max", + objective.in = obj, + const.mat = cm, + const.dir = rep("<=", n + m), + const.rhs = c(cap, rep(1, n)), + all.bin = TRUE) + + sls = sol$solution + my = sls[1:n] + for (k in 2:m) { + my = my + k * sls[((k-1)*n+1):(k*n)] + } - return(list(ksack = S$xs, value = S$vs, btracks = S$bs)) + return(list(ksack = my, val = sol$objval, bs = 0)) } + diff -Nru adagio-0.7.1/R/occurs.R adagio-0.8.4/R/occurs.R --- adagio-0.7.1/R/occurs.R 1970-01-01 00:00:00.000000000 +0000 +++ adagio-0.8.4/R/occurs.R 2021-04-28 09:47:05.000000000 +0000 @@ -0,0 +1,18 @@ +## +## o c c u r s . R Finding subsequences +## + + +occurs <- function(subseq, series){ + # Find all indices i such that series[i, ..., i+m-1] == subseq + stopifnot(is.numeric(subseq), is.numeric(series)) + m <- length(subseq) + n <- length(series) + if (m > n) + return(as.integer(c())) + + inds <- seq.int(length.out = n-m+1) + for (i in seq.int(length.out = m)) + inds <- inds[ subseq[i] == series[inds + i - 1] ] + inds +} diff -Nru adagio-0.7.1/R/setcover.R adagio-0.8.4/R/setcover.R --- adagio-0.7.1/R/setcover.R 1970-01-01 00:00:00.000000000 +0000 +++ adagio-0.8.4/R/setcover.R 2021-04-25 17:53:32.000000000 +0000 @@ -0,0 +1,26 @@ +## +## s e t c o v e r . R Set cover problem +## + +setcover <- function(Sets, weights = rep(1, nrow(Sets))) { + stopifnot(is.numeric(Sets), is.numeric(weights)) + if (!is.matrix(Sets) || !all(Sets %in% c(0, 1))) { + stop("Argument 'Sets' must be a matrix of zeros and ones.") + } + n = nrow(Sets) # number of sets + m = ncol(Sets) # size of universe + if (any(colSums(Sets) == 0)) { + stop("Not all elements covered by sets in 'Sets'.") + } + + sol <- lp(direction = "min", + objective.in = weights, + const.mat = Sets, + const.dir = rep(">=", m), + const.rhs = rep(1, m), + transpose.constraints = FALSE, + all.bin =TRUE) + + sets <- which(sol$solution == 1) + return(list(sets = sets, objective = sol$objval)) +} diff -Nru adagio-0.7.1/R/subsetsum.R adagio-0.8.4/R/subsetsum.R --- adagio-0.7.1/R/subsetsum.R 2012-02-16 19:37:07.000000000 +0000 +++ adagio-0.8.4/R/subsetsum.R 2021-04-28 18:18:46.000000000 +0000 @@ -3,6 +3,21 @@ ## +# A simple test to find the largest summable value +sss_test <- function(S, t) { + stopifnot(is.numeric(S), is.numeric(t)) + S = sort(S, decreasing = TRUE) + n <- length(S) + L <- c(0) + for (i in 1:n) { + L <- sort(unique(c(L, L + S[i]))) + L <- L[L <= t] + if(max(L) == t) break + } + return(max(L)) +} + + # Assume S decreasing, no elements > t, total sum >= t subsetsum <- function(S, t, method = "greedy") { stopifnot(is.numeric(S), is.numeric(t)) @@ -14,11 +29,11 @@ floor(t) != ceiling(t) || t <= 0) stop("Arguments 'S' and 't' must be positive integer vectors.") - if (any(S >= t)) - stop("No element of 'S' shall be greater or equal to 't'.") + if (any(S > t)) + stop("No element of 'S' shall be greater than 't'.") if (sum(S) < t) { warning("Total sum of 'S' is smaller than 't'; no solution possible.") - return(NA) + return(list(val = NA, inds = NA)) } method <- pmatch(method, c("greedy", "dynamic")) diff -Nru adagio-0.7.1/src/assgn.f adagio-0.8.4/src/assgn.f --- adagio-0.7.1/src/assgn.f 2018-05-16 19:35:09.000000000 +0000 +++ adagio-0.8.4/src/assgn.f 1970-01-01 00:00:00.000000000 +0000 @@ -1,268 +0,0 @@ -*** ******************************************************************** - subroutine assgn (n,a,c,t,iwk,ierr) -C ------------------- -C solution of the assignment problem -C ------------------- - integer a(n,*), c(n), t, iwk(*) - - i1 = n + 1 - i2 = i1 + n - i3 = i2 + n - i4 = i3 + n + 1 - i5 = i4 + n - i6 = i5 + n - call assgn1(n,a,c,t,iwk(1),iwk(i1),iwk(i2),iwk(i3),iwk(1), - + iwk(i3),iwk(i4),iwk(i5),iwk(i6),ierr) - return - end - -*** ******************************************************************** - subroutine assgn1(n,a,c,t,ch,lc,lr,lz,nz,rh,slc,slr, - + u,ierr) - integer a(n,*), c(n), ch(n), lc(n), lr(n), lz(n), - + nz(n), rh(*), slc(n), slr(n), u(*) - integer h, q, r, s, t -C ___________________________________________________________________ -C -C This subroutine solves the square assignment problem -C The meaning of the input parameters is -C n = number of rows and columns of the cost matrix -C a(i,j) = element in row i and column j of the cost matrix -C ( at the end of computation the elements of a are changed) -C the meaning of the output parameters is -C c(j) = row assigned to column j (j=1,n) -C t = cost of the optimal assignment -C All parameters are integer -C The meaning of the local variables is -C a(i,j) = element of the cost matrix if a(i,j) is positive, -C column of the unassigned zero following in row i -C (i=1,n) the unassigned zero of column j (j=1,n) -C if a(i,j) is not positive -C a(i,n+1) = column of the first unassigned zero of row i -C (i=1,n) -C ch(i) = column of the next unexplored and unassigned zero -C of row i (i=1,n) -C lc(j) = label of column j (j=1,n) -C lr(i) = label of row i (i=1,n) -C lz(i) = column of the last unassigned zero of row i(i=1,n) -C nz(i) = column of the next unassigned zero of row i(i=1,n) -C rh(i) = unexplored row following the unexplored row i -C (i=1,n) -C rh(n+1) = first unexplored row -C slc(k) = k-th element contained in the set of the labelled -C columns -C slr(k) = k-th element contained in the set of the labelled -C rows -C u(i) = unassigned row following the unassigned row i -C (i=1,n) -C u(n+1) = first unassigned row -C ierr = 0 if the routine terminates successfully. otherwise -C ierr = 1 -C -C The vectors c,ch,lc,lr,lz,nz,slc,slr must be dimensioned -C at least at (n), the vectors rh,u at least at (n+1), -C and the matrix a at least at (n,n+1). To save storage -C lz and rh may use the same storage area, and nz and ch -C may use the same storage area. -C ___________________________________________________________________ - -C initialization - maxnum = 2147483647 - ierr = 0 - np1 = n+1 - do 10 j=1,n - c(j) = 0 - lz(j) = 0 - nz(j) = 0 - u(j) = 0 - 10 continue - u(np1) = 0 - t = 0 -C reduction of the initial cost matrix - do 40 j=1,n - s = a(1,j) - do 15 l=2,n - if ( a(l,j) .lt. s ) s = a(l,j) - 15 continue - if (s) 20,40,30 - 20 mm = maxnum + s - if (t .lt. -mm) go to 400 - t = t + s - do 25 i = 1,n - if (a(i,j) .gt. mm) go to 400 - a(i,j) = a(i,j) - s - 25 continue - go to 40 - 30 mm = maxnum - s - if (t .gt. mm) go to 400 - t = t + s - do 35 i = 1,n - a(i,j) = a(i,j) - s - 35 continue - 40 continue - do 70 i=1,n - q = a(i,1) - do 50 l=2,n - if ( a(i,l) .lt. q ) q = a(i,l) - 50 continue - mm = maxnum - q - if (t .gt. mm) go to 400 - t = t + q - l = np1 - do 60 j=1,n - a(i,j) = a(i,j)-q - if ( a(i,j) .ne. 0 ) go to 60 - a(i,l) = -j - l = j - 60 continue - 70 continue -C choice of the initial solution - k = np1 - do 140 i=1,n - lj = np1 - j = -a(i,np1) - 80 if ( c(j) .eq. 0 ) go to 130 - lj = j - j = -a(i,j) - if ( j .ne. 0 ) go to 80 - lj = np1 - j = -a(i,np1) - 90 r = c(j) - lm = lz(r) - m = nz(r) - 100 if ( m .eq. 0 ) go to 110 - if ( c(m) .eq. 0 ) go to 120 - lm = m - m = -a(r,m) - go to 100 - 110 lj = j - j = -a(i,j) - if ( j .ne. 0 ) go to 90 - u(k) = i - k = i - go to 140 - 120 nz(r) = -a(r,m) - lz(r) = j - a(r,lm) = -j - a(r,j) = a(r,m) - a(r,m) = 0 - c(m) = r - 130 c(j) = i - a(i,lj) = a(i,j) - nz(i) = -a(i,j) - lz(i) = lj - a(i,j) = 0 - 140 continue -C research of a new assignment - 150 if ( u(np1) .eq. 0 ) return - do 160 i=1,n - ch(i) = 0 - lc(i) = 0 - lr(i) = 0 - rh(i) = 0 - 160 continue - rh(np1) = -1 - kslc = 0 - kslr = 1 - r = u(np1) - lr(r) = -1 - slr(1) = r - if ( a(r,np1) .eq. 0 ) go to 220 - 170 l = -a(r,np1) - if ( a(r,l) .eq. 0 ) go to 180 - if ( rh(r) .ne. 0 ) go to 180 - rh(r) = rh(np1) - ch(r) = -a(r,l) - rh(np1) = r - 180 if ( lc(l) .eq. 0 ) go to 200 - if ( rh(r) .eq. 0 ) go to 210 - 190 l = ch(r) - ch(r) = -a(r,l) - if ( a(r,l) .ne. 0 ) go to 180 - rh(np1) = rh(r) - rh(r) = 0 - go to 180 - 200 lc(l) = r - if ( c(l) .eq. 0 ) go to 360 - kslc = kslc+1 - slc(kslc) = l - r = c(l) - lr(r) = l - kslr = kslr+1 - slr(kslr) = r - if ( a(r,np1) .ne. 0 ) go to 170 - 210 continue - if ( rh(np1) .gt. 0 ) go to 350 -C reduction of the current cost matrix - 220 h = maxnum - do 240 j=1,n - if ( lc(j) .ne. 0 ) go to 240 - do 230 k=1,kslr - i = slr(k) - if ( a(i,j) .lt. h ) h = a(i,j) - 230 continue - 240 continue - mm = maxnum - h - if (mm .eq. 0 .or. t .gt. mm) go to 400 - t = t + h - do 290 j=1,n - if ( lc(j) .ne. 0 ) go to 290 - do 280 k=1,kslr - i = slr(k) - a(i,j) = a(i,j)-h - if ( a(i,j) .ne. 0 ) go to 280 - if ( rh(i) .ne. 0 ) go to 250 - rh(i) = rh(np1) - ch(i) = j - rh(np1) = i - 250 l = np1 - 260 nl = -a(i,l) - if ( nl .eq. 0 ) go to 270 - l = nl - go to 260 - 270 a(i,l) = -j - 280 continue - 290 continue - if ( kslc .eq. 0 ) go to 350 - do 340 i=1,n - if ( lr(i) .ne. 0 ) go to 340 - do 330 k=1,kslc - j = slc(k) - if ( a(i,j) .gt. 0 ) go to 320 - l = np1 - 300 nl = - a(i,l) - if ( nl .eq. j ) go to 310 - l = nl - go to 300 - 310 a(i,l) = a(i,j) - a(i,j) = h - go to 330 - 320 mm = maxnum - h - if (a(i,j) .gt. mm) go to 400 - a(i,j) = a(i,j) + h - 330 continue - 340 continue - 350 r = rh(np1) - go to 190 -C assignment of a new row - 360 c(l) = r - m = np1 - 370 nm = -a(r,m) - if ( nm .eq. l ) go to 380 - m = nm - go to 370 - 380 a(r,m) = a(r,l) - a(r,l) = 0 - if ( lr(r) .lt. 0 ) go to 390 - l = lr(r) - a(r,l) = a(r,np1) - a(r,np1) = -l - r = lc(l) - go to 360 - 390 u(np1) = u(r) - u(r) = 0 - go to 150 -C error return - integer overflow occurs - 400 ierr = 1 - return - end diff -Nru adagio-0.7.1/src/maxsubf.f adagio-0.8.4/src/maxsubf.f --- adagio-0.7.1/src/maxsubf.f 2018-05-16 19:35:09.000000000 +0000 +++ adagio-0.8.4/src/maxsubf.f 1970-01-01 00:00:00.000000000 +0000 @@ -1,77 +0,0 @@ -C Maximal Sum Subarray Problem - subroutine maxsubf(x, n, s, i1, i2) - integer n, i1, i2, j1, j2 - double precision x(n) - double precision s, ss - - ss = 0 - j1 = 1 - j2 = 1 - - do 3000 i = 1,n - if (ss .gt. -x(i)) then - ss = ss + x(i) - j2 = i - if (ss .gt. s) then - s = ss - i1 = j1 - i2 = j2 - endif - else - ss = 0 - j1 = i+1 - j2 = i+1 - endif - 3000 continue - - return - end - - -C Maximal Sum Subrectangle Problem - subroutine maxsub2f(a, s, n, m, fmax, mind, aa, b) - integer n, m, mind(4) - integer mi1, mi2 - double precision a(n, m), aa(m), b, s(n+1, m) - double precision fsum, fmax - - do 3000 i = 1, m - s(1, i) = 0 - 3000 continue - - do 3002 i = 2, n+1 - b = a(i-1, 1) - s(i, 1) = s(i-1, 1) + b - do 3004 j = 2, m - b = b + a(i-1, j) - s(i, j) = s(i-1, j) + b - 3004 continue - 3002 continue - - fsum = 0.0 - mi1 = 0 - mi2 = 0 - - do 3006 i = 2, n+1 - do 3008 j = i, n+1 - aa(1) = s(j, 1) - s(i-1, 1) - b = aa(1) - do 3010 k = 2, n - aa(k) = s(j, k) - s(i-1, k) - b - b = b + aa(k) - 3010 continue - - call maxsubf(aa, m, fsum, mi1, mi2) - if (fsum .gt. fmax) then - fmax = fsum - mind(1) = i-1 - mind(2) = j-1 - mind(3) = mi1 - mind(4) = mi2 - endif - - 3008 continue - 3006 continue - - return - end diff -Nru adagio-0.7.1/src/mkp.f adagio-0.8.4/src/mkp.f --- adagio-0.7.1/src/mkp.f 2018-05-16 19:35:09.000000000 +0000 +++ adagio-0.8.4/src/mkp.f 1970-01-01 00:00:00.000000000 +0000 @@ -1,397 +0,0 @@ -*** ******************************************************************** -*** S u b r o u t i n e M K P -*** ******************************************************************** - - subroutine mkp (n,m,p,w,k,bck,xstar,vstar,wk,iwk,num) -* ____________________________________________________________________ -* -* Subroutine to solve a 0-1 multiple knapsack problem of n -* items (with n .ge. 2) and m knapsacks (with m .ge. 1 , -* i.e. also suitable for a 0-1 single knapsack problem). -* The problem to be solved is -* -* maximize vstar = p(1)*(x(1,1) + ... + x(m,1)) + ... -* ... + p(n)*(x(1,n) + ... + x(m,n)) -* subject to -* w(1)*x(i,1) + ... + w(n)*x(i,n) .le. k(i) for i=1,...,m -* x(1,j) + ... + x(m,j) .le. 1 for j=1,...,n -* x(i,j) = 0 or 1 for i=1,...,m , j=1,...,n , -* -* where all p(j), w(j), and k(i) are positive integers. -* Before mkp is called, array k must be sorted so that -* k(1) .le. k(2) .le. ... .le. k(m) . -* -* Meaning of the input parameters ... -* -* n = number of items. -* m = number of knapsacks. -* p(j) = profit of item j (j=1,...,n) . -* w(j) = weight of item j (j=1,...,n) . -* k(i) = capacity of knapsack i (i=1,...,m) . -* bck = -1 if exact solution is required. -* = maximum number of backtrackings to be performed, if -* heuristic solution is required. -* wk = real work space of dimension n. -* iwk = work space of dimension .ge. 5*m + 14*n + 4*m*n + 3 -* num = dimension of iwk -* -* Meaning of the output parameters ... -* -* xstar(j) = 0 if item j is not in the optimal solution -* (i.e. if x(i,j) = 0 for all i ). -* = knapsack where item j is inserted, otherwise -* (i.e. if x(xstar(j),j) = 1 ). -* vstar = value of the optimal solution if vstar .gt. 0. -* = error condition (infeasibility or triviality) -* in the input data if vstar .lt. 0 . -* = -1 if n .lt. 2 or m .lt. 1 . -* = -2 if some p(j) , w(j) or k(i) are not -* positive. -* = -3 if a knapsack cannot contain any item. -* = -4 if an item cannot fit into any knapsack. -* = -5 if knapsack m contains all the items. -* = -7 if array k is not correctly sorted. -* = -8 if num .lt. 5*m + 14*n + 4*m*n + 3. -* -* (in all the above cases array xstar is not -* defined). -* -* All the parameters except wk are of integer type. when mkp -* terminates, all the input parameters are unchanged except -* bck, which gives the number of backtrackings performed. -* ____________________________________________________________________ - - integer p(n),w(n),k(m),xstar(n),bck,vstar,iwk(num) - real wk(n) - integer bb, bl, x, xl - integer b, ubb - integer f, pbl, q, v - integer bs, ps, ws, xs -C -C check the input data -C - if (m .lt. 1 .or. n .lt. 2) go to 100 - mn = m*n - if (num .lt. 5*m + 14*n + 4*mn + 3) go to 160 - - if (p(1) .le. 0 .or. w(1) .le. 0) go to 110 - ap = p(1) - aw = w(1) - wk(1) = -ap/aw - maxw = w(1) - minw = w(1) - isumw = w(1) - do 10 j = 2,n - if (p(j) .le. 0 .or. w(j) .le. 0) go to 110 - ap = p(j) - aw = w(j) - wk(j) = -ap/aw - if (w(j) .gt. maxw) maxw = w(j) - if (w(j) .lt. minw) minw = w(j) - isumw = isumw + w(j) - 10 continue - - if (k(1) .le. 0) go to 110 - if (m .eq. 1) go to 30 - do 20 i = 2,m - if (k(i) .le. 0) go to 110 - if (k(i) .lt. k(i-1)) go to 150 - 20 continue - - 30 if (minw .gt. k(1)) go to 120 - if (maxw .gt. k(m)) go to 130 - if (isumw .le. k(m)) go to 140 - vstar = 0 -C -C reorder the arrays p and w so that -C p(j)/w(j) .ge. p(j+1)/w(j+1) -C - n5 = 5*n - do 40 j = 1,n - jj = n5 + j - iwk(jj) = j - 40 continue - call risort (wk, iwk(n5 + 1), n) - - do 50 j = 1,n - iwk(j) = p(j) - jn = j + n - iwk(jn) = w(j) - 50 continue - - do 60 j = 1,n - jj = n5 + j - l = iwk(jj) - p(j) = iwk(l) - npl = n + l - w(j) = iwk(npl) - 60 continue -C -C partition the work space iwk -C - lx = jj + 1 - lxi = lx + n - bs = lxi + n - xs = bs + n - ubb = xs + n - - np1 = n + 1 - b = ubb + n - ps = b + np1 - ws = ps + np1 - - f = ws + np1 - pbl = f + m - q = pbl + m - v = q + m - - bb = v + m - x = bb + mn - xl = x + mn - - bl = xl + mn -C -C solve the problem -C - call mkp1 (n, m, p, w, k, bck, xstar, vstar, np1, n5, - + iwk(bb), iwk(bl), iwk(x), iwk(xl), - + iwk(b), iwk(ubb), iwk(lx), iwk(lxi), - + iwk(f), iwk(pbl), iwk(q), iwk(v), - + iwk(bs), iwk(ps), iwk(ws), iwk(xs), iwk(1)) -C -C restore the initial ordering to p and w, -C and reorder xstar accordingly -C - do 70 j = 1,n - iwk(j) = p(j) - jn = j + n - iwk(jn) = w(j) - jnn = jn + n - iwk(jnn) = xstar(j) - 70 continue - - do 80 j = 1,n - jj = n5 + j - l = iwk(jj) - p(l) = iwk(j) - jn = j + n - w(l) = iwk(jn) - jnn = jn + n - xstar(l) = iwk(jnn) - 80 continue - return -C -C error return -C - 100 vstar = -1 - return - 110 vstar = -2 - return - 120 vstar = -3 - return - 130 vstar = -4 - return - 140 vstar = -5 - return - 150 vstar = -7 - return - 160 vstar = -8 - return - end - - -*** ******************************************************************** -*** S u b r o u t i n e M K P 1 -*** ******************************************************************** - - subroutine mkp1 (n, m, p, w, k, bck, xstar, vstar, np1, n5, - + bb, bl, x, xl, b, ubb, lx, lxi, - + f, pbl, q, v, bs, ps, ws, xs, iwk) - -* -------------------------------------------------------------------- -* meaning of the main internal variables and arrays ... -* -* i = knapsack currently considered. -* lb = lower bound on the optimal solution. -* ub = upper bound on the optimal solution. -* vb = value of the current solution. -* x(i,j) = 1 if item j is inserted in knapsack i in -* the current solution. -* = 0 otherwise. -* f(i) = pointer to the last item inserted in knapsack i -* ( = -1 if knapsack i is empty). -* bb(i,j) = pointer to the item inserted in knapsack i -* just before item j ( = -1 if j is the first -* item inserted in knapsack i ). -* q(i) = current available capacity of knapsack i . -* b(j) = 1 if item j is not inserted in any knapsack. -* = 0 if item j is inserted in a knapsack. -* pbl(i) = number of the items which can be inserted in -* knapsack i . -* bl(i,s) = pointer to the s-th item which can be inserted -* in knapsack i . -* xl(i,j) = 1 if item j was inserted in knapsack i in -* the last execution of subroutine pi1. -* = 0 otherwise. -* iwk work space for the subroutine sknp. -* -------------------------------------------------------------------- - - integer p(n), w(n), k(m), xstar(n), bck, vstar - integer bb(m,n), bl(m,np1), x(m,n), xl(m,n) - integer b(np1), ubb(n), lx(n), lxi(n) - integer f(m), pbl(m), q(m), v(m) - integer bs(n), ps(np1), ws(np1), xs(n), iwk(n5) - integer s, u, ub, vb - - if (m .eq. 1) go to 250 -C -C step 1 (initialization) -C - jbck = bck - bck = 0 - kub = 0 - n1 = n + 1 - b(n1) = 1 - m1 = m - 1 - do 40 j=1,n - b(j) = 1 - do 30 i=1,m - x(i,j) = 0 - bb(i,j) = 0 - 30 continue - 40 continue - do 50 i=1,m1 - q(i) = k(i) - f(i) = -1 - 50 continue - q(m) = k(m) - vstar = 0 - vb = 0 - i = 1 - call sigma1 (n,m,p,w,k,1,b,kub,ub,np1,n5,lx,lr, - + bs,ps,ws,xs,iwk) - do 60 j=1,n - lxi(j) = lx(j) - 60 continue - lri = lr - lubi = ub - iflag = 0 -C -C step 2 (heuristic) -C - 70 kub = vstar - vb - call pi1 (n,m,p,w,q,i,b,bb,kub,bl,lb,pbl,v,xl, - + np1,n5,bs,ps,ws,xs,iwk) - if ( lb + vb .le. vstar ) go to 140 - vstar = lb + vb - do 90 j=1,n - xstar(j) = 0 - do 80 s=1,i - if ( x(s,j) .eq. 0 ) go to 80 - xstar(j) = s - go to 90 - 80 continue - 90 continue - ip = pbl(i) - if ( ip .eq. 0 ) go to 110 - do 100 j=1,ip - jj = bl(i,j) - if ( xl(i,j) .eq. 1 ) xstar(jj) = i - 100 continue - 110 i1 = i + 1 - do 130 ii=i1,m - ip = pbl(ii) - if ( ip .eq. 0 ) go to 130 - do 120 j=1,ip - jj = bl(ii,j) - if ( xl(ii,j) .eq. 1 ) xstar(jj) = ii - 120 continue - 130 continue - if ( ub .eq. lb ) go to 200 -C -C step 3 (updating) -C - 140 if ( v(i) .eq. 0 ) go to 180 - iuv = ub + vb - u = pbl(i) - ibv = 0 - do 170 s=1,u - if ( xl(i,s) .eq. 0 ) go to 170 - j = bl(i,s) - x(i,j) = 1 - q(i) = q(i) - w(j) - vb = vb + p(j) - b(j) = 0 - bb(i,j) = f(i) - ubb(j) = iuv - if ( iflag .eq. 1 ) go to 150 - lub = iuv - lj = j - li = i - 150 f(i) = j - ibv = ibv + p(j) - if ( ibv .eq. v(i) ) go to 180 - call parc (i,i,ub,iflag,vb,lub,lj,li,f,bb,q,b,n,m,np1, - + lx,lxi,lr,lri,lubi) - if ( iflag .eq. 1 ) go to 160 - kub = vstar - vb - call sigma1 (n,m,p,w,q,i,b,kub,ub,np1,n5,lx,lr, - + bs,ps,ws,xs,iwk) - lj = n1 - 160 iuv = ub + vb - if ( iuv .le. vstar ) go to 200 - 170 continue - 180 if ( i .eq. m - 1 ) go to 200 - ip1 = i + 1 - call parc (ip1,i,ub,iflag,vb,lub,lj,li,f,bb,q,b,n,m,np1, - + lx,lxi,lr,lri,lubi) - if ( iflag .eq. 1 ) go to 190 - kub = vstar - vb - call sigma1 (n,m,p,w,q,ip1,b,kub,ub,np1,n5,lx,lr, - + bs,ps,ws,xs,iwk) - lj = n1 - 190 if ( ub + vb .le. vstar ) go to 200 - i = i + 1 - go to 140 -C -C step 4 (backtracking) -C - 200 if ( i .gt. 0 ) go to 210 - bck = bck - 1 - return - 210 if ( bck .eq. jbck ) return - bck = bck + 1 - if ( f(i) .ne. (-1) ) go to 230 - do 220 j=1,n - bb(i,j) = 0 - 220 continue - i = i - 1 - go to 200 - 230 j = f(i) - x(i,j) = 0 - b(j) = 1 - vb = vb - p(j) - q(i) = q(i) + w(j) - do 240 s=1,n - if ( bb(i,s) .eq. j ) bb(i,s) = 0 - 240 continue - f(i) = bb(i,j) - if ( ubb(j) .le. vstar ) go to 200 - ub = ubb(j) - vb - iflag = 1 - go to 70 -C -C particular case (0-1 single knapsack problem) -C - 250 k1 = k(1) - do 260 j=1,n - ps(j) = p(j) - ws(j) = w(j) - 260 continue - call sknp (n,k1,0,vstar,n,np1,n5,ps,ws,xs,iwk) - do 270 j=1,n - xstar(j) = xs(j) - 270 continue - bck = 0 - return - end diff -Nru adagio-0.7.1/src/nscw_util.f adagio-0.8.4/src/nscw_util.f --- adagio-0.7.1/src/nscw_util.f 2018-05-16 19:35:09.000000000 +0000 +++ adagio-0.8.4/src/nscw_util.f 1970-01-01 00:00:00.000000000 +0000 @@ -1,469 +0,0 @@ -*** ******************************************************************** - subroutine risort (a, m, n) -c----------------------------------------------------------------------- -c the shell sorting procedure is used to reorder the elements of a -c so that a(i).le.a(i+1) for i=1,...,n-1. the same permutations are -c performed on m that are performed on a. it is assumed that n.ge.1. -c----------------------------------------------------------------------- - real a(n) - integer m(n), t - integer k(10) -c------------------------ - data k(1)/1/, k(2)/4/, k(3)/13/, k(4)/40/, k(5)/121/, k(6)/364/, - 1 k(7)/1093/, k(8)/3280/, k(9)/9841/, k(10)/29524/ -c------------------------ -c -c selection of the increments k(i) = (3**i-1)/2 -c - if (n .lt. 2) return - imax = 1 - do 10 i = 3,10 - if (n .le. k(i)) go to 20 - imax = imax + 1 - 10 continue -c -c stepping through the increments k(imax),...,k(1) -c - 20 i = imax - do 40 ii = 1,imax - ki = k(i) -c -c sorting elements that are ki positions apart -c so that a(j).le.a(j+ki) for j=1,...,n-ki -c - jmax = n - ki - do 32 j = 1,jmax - l = j - ll = j + ki - s = a(ll) - t = m(ll) - 30 if (s .ge. a(l)) go to 31 - a(ll) = a(l) - m(ll) = m(l) - ll = l - l = l - ki - if (l .gt. 0) go to 30 - 31 a(ll) = s - m(ll) = t - 32 continue -c - 40 i = i - 1 - return - end - -*** ******************************************************************** - subroutine sigma1 (n,m,p,w,q,i,b,kub,ub,np1,n5,lx,lr, - * bs,ps,ws,xs,iwk) -c -c subroutine to compute an upper bound ub on the best -c final solution which can be obtained from the current -c solution. -c - integer p(n),w(n),q(m),b(np1),ub,iwk(n5) - integer lx(n),bs(n),ps(np1),ws(np1),xs(n) - integer qs,sb -c - ns = 0 - qs = 0 - do 10 j=i,m - qs = qs + q(j) - 10 continue - sb = 0 - do 20 j=1,n - lx(j) = 0 - if ( b(j) .eq. 0 ) go to 20 - ns = ns + 1 - bs(ns) = j - ps(ns) = p(j) - ws(ns) = w(j) - sb = sb + w(j) - 20 continue - if ( sb .gt. qs ) go to 40 - lr = qs - sb - ub = 0 - if ( ns .eq. 0 ) return - do 30 j=1,ns - ub = ub + ps(j) - xs(j) = 1 - 30 continue - go to 50 - 40 call sknp (ns,qs,kub,ub,n,np1,n5,ps,ws,xs,iwk) - lr = qs - 50 do 60 j=1,ns - jj = bs(j) - lx(jj) = xs(j) - 60 continue - return - end - -*** ******************************************************************** - subroutine pi1 (n,m,p,w,q,i,b,bb,kub,bl,lb,pbl,v,xl, - * np1,n5,bs,ps,ws,xs,iwk) -c -c subroutine to compute a feasible solution to the current -c problem. the solution is stored in array xl , the -c corresponding value in lb . -c - integer bb(m,n),bl(m,np1),xl(m,n),iwk(n5) - integer p(n),w(n),q(m),b(np1),pbl(m),v(m) - integer bs(n),ps(np1),ws(np1),xs(n) - integer pb,qs,sb,u -c -c step 1 -c - u = 0 - do 10 j=1,n - if ( b(j) .eq. 0 ) go to 10 - u = u + 1 - bs(u) = j - 10 continue - do 20 j=i,m - pbl(j) = 0 - v(j) = 0 - 20 continue - lb = 0 - ikub = kub - if ( u .eq. 0 ) return - ns = 0 - sb = 0 - do 30 j=1,u - jj = bs(j) - if ( bb(i,jj) .ne. 0 ) go to 30 - if ( w(jj) .gt. q(i) ) go to 30 - ns = ns + 1 - sb = sb + w(jj) - bl(i,ns) = jj - ps(ns) = p(jj) - ws(ns) = w(jj) - 30 continue - ii = i -c -c step 2 -c - 40 pbl(ii) = ns - if ( sb .gt. q(ii) ) go to 60 - pb = 0 - if ( ns .eq. 0 ) go to 80 - do 50 j=1,ns - pb = pb + ps(j) - xl(ii,j) = 1 - 50 continue - go to 80 - 60 qs = q(ii) - kub = 0 - if ( ii .eq. m ) kub = ikub - call sknp (ns,qs,kub,pb,n,np1,n5,ps,ws,xs,iwk) - do 70 j=1,ns - xl(ii,j) = xs(j) - 70 continue - 80 lb = lb + pb - ikub = ikub - pb - v(ii) = pb - bl(ii,ns+1) = n + 1 -c -c step 3 -c - if ( ii .eq. m ) return - jb = 1 - jbs = 0 - do 100 j=1,u - if ( bs(j) .lt. bl(ii,jb) ) go to 90 - jb = jb + 1 - if ( xl(ii,jb-1) .eq. 1 ) go to 100 - 90 jbs = jbs + 1 - bs(jbs) = bs(j) - 100 continue - u = jbs - if ( u .eq. 0 ) return - ns = 0 - sb = 0 - ii = ii + 1 - do 110 j=1,u - jj = bs(j) - if( w(jj) .gt. q(ii) ) go to 110 - ns = ns + 1 - sb = sb + w(jj) - bl(ii,ns) = jj - ps(ns) = p(jj) - ws(ns) = w(jj) - 110 continue - go to 40 - end - -*** ******************************************************************** - subroutine parc (i,ii,ub,iflag,vb,lub,lj,li,f,bb,q,b,n,m,np1, - * lx,lxi,lr,lri,lubi) -c -c subroutine for parametric computation of the upper bounds. -c - integer f(m),bb(m,n),q(m),b(np1),ub,vb,r,s - integer lx(n),lxi(n) -c - iflag = 0 - if ( b(lj) .ne. 0 ) go to 60 - i1 = i - 1 - if ( i1 .lt. li ) go to 20 - iq = 0 - do 10 r=li,i1 - iq = iq + q(r) - 10 continue - if ( iq .gt. lr ) return - 20 r = ii - s = f(r) - 30 if ( s .ne. (-1) ) go to 40 - r = r - 1 - s = f(r) - go to 30 - 40 if ( lx(s) .eq. 0 ) return - if ( s .eq. lj ) go to 50 - s = bb(r,s) - go to 30 - 50 ub = lub - vb - iflag = 1 - return - 60 i1 = i - 1 - if ( i1 .lt. 1 ) go to 80 - iq = 0 - do 70 r=1,i1 - iq = iq + q(r) - 70 continue - if ( iq .gt. lri ) return - 80 do 90 j=1,n - if ( b(j) .eq. 1 ) go to 90 - if ( lxi(j) .eq. 0 ) return - 90 continue - ub = lubi - vb - iflag = 1 - return - end - -*** ******************************************************************** - subroutine sknp (ns,qs,kub,vs,n,np1,n5,ps,ws,xs,iwk) -c -c subroutine to solve the 0-1 single knapsack problem -c -c maximize vs = ps(1)*xs(1) + ... + ps(ns)*xs(ns) -c subject to ws(1)*xs(1) + ... + ws(ns)*xs(ns) .le. qs -c xs(j) = 0 or 1 for j=1,...,ns -c vs .gt. kub -c -c this subroutine is a modified version of subroutine kp01 -c which appeared in computing 21, 81-86(1978). -c - integer qs, vs - integer ps(np1), ws(np1), xs(n), iwk(n5) -c - i1 = 1 - i2 = i1 + n - i3 = i2 + n - i4 = i3 + n - i5 = i4 + n - call sknp1 (ns,qs,kub,vs,n,np1,ps,ws,xs,iwk(i1),iwk(i2), - * iwk(i3),iwk(i4),iwk(i5)) - return - end - -*** ******************************************************************** - subroutine sknp1 (ns,qs,kub,vs,n,np1,ps,ws,xs,d,min, - * pbar,wbar,zbar) -c -c subroutine to solve the 0-1 single knapsack problem -c -c maximize vs = ps(1)*xs(1) + ... + ps(ns)*xs(ns) -c subject to ws(1)*xs(1) + ... + ws(ns)*xs(ns) .le. qs -c xs(j) = 0 or 1 for j=1,...,ns -c vs .gt. kub -c -c this subroutine is a modified version of subroutine kp01 -c which appeared in computing 21, 81-86(1978). -c - integer qs,vs,diff,pr,r,t - integer ps(np1),ws(np1),xs(n) - integer d(n),min(n),pbar(n),wbar(n),zbar(n) -c - vs = kub - ip = 0 - ms = qs - do 10 l=1,ns - ll = l - if ( ws(l) .gt. ms ) go to 20 - ip = ip + ps(l) - ms = ms - ws(l) - 10 continue - 20 ll = ll - 1 - if ( ms .eq. 0 ) go to 50 - ps(ns+1) = 0 - ws(ns+1) = qs + 1 - lim = ip + (ms*ps(ll+2))/ws(ll+2) - a = ip + ps(ll+1) - b = (ws(ll+1) - ms)*ps(ll) - c = ws(ll) - lim1 = a - b/c - if ( lim1 .gt. lim ) lim = lim1 - if ( lim .le. vs ) return - mink = qs + 1 - min(ns) = mink - do 30 j=2,ns - kk = ns + 2 - j - if ( ws(kk) .lt. mink ) mink = ws(kk) - min(kk-1) = mink - 30 continue - do 40 j=1,ns - d(j) = 0 - 40 continue - pr = 0 - lold = ns - ii = 1 - go to 170 - 50 if ( vs .ge. ip ) return - vs = ip - do 60 j=1,ll - xs(j) = 1 - 60 continue - nn = ll + 1 - do 70 j=nn,ns - xs(j) = 0 - 70 continue - qs = 0 - return - 80 if ( ws(ii) .le. qs ) go to 90 - ii1 = ii + 1 - if ( vs .ge. (qs*ps(ii1))/ws(ii1) + pr ) go to 280 - ii = ii1 - go to 80 - 90 ip = pbar(ii) - ms = qs - wbar(ii) - in = zbar(ii) - ll = ns - if ( in .gt. ns) go to 110 - do 100 l=in,ns - ll = l - if ( ws(l) .gt. ms ) go to 160 - ip = ip + ps(l) - ms = ms - ws(l) - 100 continue - 110 if ( vs .ge. ip + pr ) go to 280 - vs = ip + pr - mfirst = ms - nn = ii - 1 - do 120 j=1,nn - xs(j) = d(j) - 120 continue - do 130 j=ii,ll - xs(j) = 1 - 130 continue - if ( ll .eq. ns ) go to 150 - nn = ll + 1 - do 140 j=nn,ns - xs(j) = 0 - 140 continue - 150 if ( vs .ne. lim ) go to 280 - qs = mfirst - return - 160 l = ll - ll = ll - 1 - if ( ms .eq. 0 ) go to 110 - if ( vs .ge. pr + ip + (ms*ps(l))/ws(l) ) go to 280 - 170 wbar(ii) = qs - ms - pbar(ii) = ip - zbar(ii) = ll + 1 - d(ii) = 1 - nn = ll - 1 - if ( nn .lt. ii ) go to 190 - do 180 j=ii,nn - wbar(j+1) = wbar(j) - ws(j) - pbar(j+1) = pbar(j) - ps(j) - zbar(j+1) = ll + 1 - d(j+1) = 1 - 180 continue - 190 j1 = ll + 1 - do 200 j=j1,lold - wbar(j) = 0 - pbar(j) = 0 - zbar(j) = j - 200 continue - lold = ll - qs = ms - pr = pr + ip - if ( ll - (ns - 2) ) 240, 220, 210 - 210 ii = ns - go to 250 - 220 if ( qs .lt. ws(ns) ) go to 230 - qs = qs - ws(ns) - pr = pr + ps(ns) - d(ns) = 1 - 230 ii = ns - 1 - go to 250 - 240 ii = ll + 2 - if ( qs .ge. min(ii-1) ) go to 80 - 250 if ( vs .ge. pr ) go to 270 - vs = pr - do 260 j=1,ns - xs(j) = d(j) - 260 continue - mfirst = qs - if ( vs .eq. lim ) return - 270 if ( d(ns) .eq. 0 ) go to 280 - d(ns) = 0 - qs = qs + ws(ns) - pr = pr - ps(ns) - 280 nn = ii - 1 - if ( nn .eq. 0 ) go to 300 - do 290 j=1,nn - kk = ii - j - if ( d(kk) .eq. 1 ) go to 310 - 290 continue - 300 qs = mfirst - return - 310 r = qs - qs = qs + ws(kk) - pr = pr - ps(kk) - d(kk) = 0 - if ( r .lt. min(kk) ) go to 320 - ii = kk + 1 - go to 80 - 320 nn = kk + 1 - ii = kk - 330 if ( vs .ge. pr + (qs*ps(nn))/ws(nn) ) go to 280 - diff = ws(nn) - ws(kk) - if ( diff ) 390, 340, 350 - 340 nn = nn + 1 - go to 330 - 350 if ( diff .gt. r ) go to 340 - if ( vs .ge. pr + ps(nn) ) go to 340 - vs = pr + ps(nn) - do 360 j=1,kk - xs(j) = d(j) - 360 continue - jj = kk + 1 - do 370 j=jj,ns - xs(j) = 0 - 370 continue - xs(nn) = 1 - mfirst = qs - ws(nn) - if ( vs .ne. lim ) go to 380 - qs = mfirst - return - 380 r = r - diff - kk = nn - nn = nn + 1 - go to 330 - 390 t = r - diff - if ( t .lt. min(nn) ) go to 340 - n1 = nn + 1 - if ( vs .ge. pr + ps(nn) + (t*ps(n1))/ws(n1) ) go to 280 - qs = qs - ws(nn) - pr = pr + ps(nn) - d(nn) = 1 - ii = nn + 1 - wbar(nn) = ws(nn) - pbar(nn) = ps(nn) - zbar(nn) = ii - do 400 j=n1,lold - wbar(j) = 0 - pbar(j) = 0 - zbar(j) = j - 400 continue - lold = nn - go to 80 - end diff -Nru adagio-0.7.1/src/ratfor/maxempty_rat.for adagio-0.8.4/src/ratfor/maxempty_rat.for --- adagio-0.7.1/src/ratfor/maxempty_rat.for 2012-03-14 12:59:47.000000000 +0000 +++ adagio-0.8.4/src/ratfor/maxempty_rat.for 1970-01-01 00:00:00.000000000 +0000 @@ -1,77 +0,0 @@ -# Converted from R code provided by Hans Werner Borchers -# ax = x-limits for region of interest -# ay = y-limits " " -# x, y = coordinates of points to avoid -# Assume x, y are sorted in y-order, e.g -# o = order(y); x <- x[o]; y <- y[o] -# n = length(x) = length(y) -# z = c(D[m], d[m], d[m+1]), d=sort(c(ax,x)), D=diff(d), m=which.max(D) -# Output: area, rect[4] -# To convert to Fortran: -# sudo apt-get install ratfor -# ratfor -o ../maxempr.f maxempr.r -SUBROUTINE maxempr(ax, ay, x, y, n, w, h, z, area, rect) - -IMPLICIT DOUBLE PRECISION (a-h,o-z) -INTEGER n -DOUBLE PRECISION ax(2), ay(2), x(n), y(n), z(3), rect(4), maxr, li -# check vertical slices -maxr = z(1) * dabs(ay(2) - ay(1)) -rect(1) = z(2) -rect(2) = ay(1) -rect(3) = z(3) -rect(4) = ay(2) - -do i=1,n { - tl = ax(1); tr = ax(2) - if (i < n) { - do j=(i+1),n { - if (x(j) > tl & x(j) < tr) { - ## check horizontal slices (j == i+1) - ## and (all) rectangles above (x(i), y(i)) - area = (tr - tl) * (y(j) - y(i)) - if (area > maxr & ((tr - tl) > w) & ((y(j) - y(i)) > h)) { - - maxr = area - rect(1) = tl - rect(2) = y(i) - rect(3) = tr - rect(4) = y(j) - } - if (x(j) > x(i)) tr = x(j) - else tl = x(j) - } - } - } - ## check open rectangles above (x(i), y(i)) - area = (tr - tl) * (ay(2) - y(i)) - if (area > maxr & ((tr - tl) > w) & - ((ay(2) - y(i)) > h)) { - - maxr = area - rect(1) = tl - rect(2) = y(i) - rect(3) = tr - rect(4) = ay(2) - } - ## check open rectangles below (x(i), y(i)) - ri = ax(2); li = ax(1) - do k=1,n { - if(y(k) < y(i) & x(k) > x(i)) ri = dmin1(ri, x(k)) - if(y(k) < y(i) & x(k) < x(i)) li = dmax1(li, x(k)) - } - area = (ri - li) * (ay(2) - y(i)) - if (area > maxr & ((ri - li) > w) & - ((y(i) - ay(1)) > h)) { - - - maxr = area - rect(1) = li - rect(2) = ay(1) - rect(3) = ri - rect(4) = y(i) - } -} -area = maxr -return -end diff -Nru adagio-0.7.1/src/ratfor/maxsub_rat.for adagio-0.8.4/src/ratfor/maxsub_rat.for --- adagio-0.7.1/src/ratfor/maxsub_rat.for 2012-02-17 15:08:10.000000000 +0000 +++ adagio-0.8.4/src/ratfor/maxsub_rat.for 1970-01-01 00:00:00.000000000 +0000 @@ -1,25 +0,0 @@ -subroutine maxsubf(x, n, s, i1, i2) -integer n, i1, i2, j1, j2 -double precision x(n) -double precision s, s1, s2 - -s1 = 0; s2 = 0 -i1 = 0; i2 = 0 -j1 = 1; j2 = 1 - -do i = 1,n { - if (s2 > -x(i)) { - s2 = s2 + x(i) - j2 = i - if (s2 > s1) { - s1 = s2 - i1 = j1; i2 = j2 - } - } else { - s2 = 0 - j1 = i+1; j2 = i+1 - } -} - -return -end