Skip to content

Commit fe0a193

Browse files
author
Lennart Noordermeer
committed
Edited buckHpr
1 parent e5c90ef commit fe0a193

File tree

1 file changed

+62
-286
lines changed

1 file changed

+62
-286
lines changed

‎R/buckHpr.R‎

Lines changed: 62 additions & 286 deletions
Original file line numberDiff line numberDiff line change
@@ -1,303 +1,79 @@
11
#' buckHpr
22
#'
3-
#' Calculate optimal bucking for all stems in a hpr file
3+
#' Optimal bucking for all stems in an HPR file using \code{buckStem()}.
44
#'
5-
#' @param XMLNode ouput from getXMLNode()
6-
#' @param PriceMatrices list of price matrices for all ProductKeys (getPriceMatrices())
7-
#' @param ProductData Matrix containing product data (getProductData())
8-
#' @param StemProfile Stem profiles for all stems in hprfile (getStemProfile())
9-
#' @param PermittedGrades list with the same lenght of assortments, each element containing the stemgrades allowed in each assortment (getPermittedGrades())
10-
#' @param SpeciesGroupDefinition getSpeciesGroupDefinition()
11-
#' @return result structure with optimum bucking solution for the stems in the input .hpr file
12-
#' @seealso getPermittedGrades, getPriceMatrices, getProductData
5+
#' @param XMLNode Output from \code{getXMLNode()}.
6+
#' @param PriceMatrices List of price matrices for all ProductKeys (from \code{getPriceMatrices()}).
7+
#' @param ProductData Data frame or data.table with product definitions (from \code{getProductData()}).
8+
#' @param StemProfile Stem profiles for all stems (from \code{getStemprofile()} / \code{getStemProfile()}).
9+
#' @param PermittedGrades Named list where names are ProductKeys and each element defines permitted stem grades
10+
#' (from \code{getPermittedGrades()}).
11+
#' @param SpeciesGroupDefinition Species group definitions (from \code{getSpeciesGroupDefinition()}).
12+
#'
13+
#' @return A \code{data.table} with the bucking solution for all processed stems, row-bound across stems.
14+
#'
15+
#' @seealso \code{\link{buckStem}}, \code{\link{getPermittedGrades}}, \code{\link{getPriceMatrices}}, \code{\link{getProductData}}
1316
#' @author Lennart Noordermeer \email{lennart.noordermeer@nmbu.no}
1417
#' @references Skogforsk 2011. Introduction to StanForD 2010. URL: Skogforsk. https://www.skogforsk.se/contentassets/1a68cdce4af1462ead048b7a5ef1cc06/stanford-2010-introduction-150826.pdf
1518
#' @export
19+
buckHpr <- function(XMLNode,
20+
PriceMatrices,
21+
ProductData,
22+
StemProfile,
23+
PermittedGrades,
24+
SpeciesGroupDefinition) {
1625

17-
buckHpr=function (XMLNode, PriceMatrices, ProductData, StemProfile, PermittedGrades,
18-
SpeciesGroupDefinition)
19-
{
20-
require(XML)
21-
require(plyr)
22-
require(dplyr)
23-
require(data.table)
26+
if (!requireNamespace("XML", quietly = TRUE)) stop("Package 'XML' is required.")
27+
if (!requireNamespace("data.table", quietly = TRUE)) stop("Package 'data.table' is required.")
2428

25-
grdFinder = function(x) {
26-
unique(StemGrade[idxstart:x])
27-
}
28-
asoFinder = function(x) {
29-
names(SGKG)[which(colSums(matrix(sapply(SGKG,
30-
FUN = function(X) all(grd[[x]] %in% X)), ncol = length(SGKG))) >
31-
0)]
32-
}
33-
DiameterValueFinder = function(x) {
34-
DiameterValue[vec[[x]]]
35-
}
36-
Rounder = function(x) {
37-
res = round_any(DV[idx][[x]], 10, floor)
38-
if (sum(idx) > 1) {
39-
res
40-
}
41-
else {
42-
list(res)
43-
}
44-
}
45-
BarkFinder = function(x) {
46-
BarkFunction(DV[idx][[x]], SpeciesGroupKey, SpeciesGroupDefinition,
47-
Top_ob = tab[idx, ][x]$Top_ob, DBH = DBH, LogLength = tab[idx,
48-
][x]$LogLength)
49-
}
50-
rowFinder = function(x) sum(commercial$LogLength[x] >=
51-
rownames[[x]] %>% as.numeric())
52-
colFinder = function(x) sum(commercial$topdiam[x] >=
53-
colnames[[x]] %>% as.numeric())
54-
priceFinder = function(x) lis[[x]][row[x], col[x]]
55-
seqVectozied = Vectorize(seq.default, vectorize.args = c("from",
56-
"to"))
57-
trackTrace = function(res, tt) {
58-
low = min(tt[, "StartPos"])
59-
while (low > 0) {
60-
id_previous = tt$CumulativeValue[order(tt$StartPos)[1]] -
61-
tt$Value[order(tt$StartPos)[1]]
62-
sub = res[res$StopPos == low, ]
63-
prev = sub[which(near(sub$CumulativeValue,
64-
id_previous)), ]
65-
if (!is.vector(prev)) {
66-
prev = prev[1, ]
67-
}
68-
tt = rbind(tt, prev)
69-
low = min(tt$StartPos)
70-
}
71-
tt = tt[nrow(tt):1, ]
72-
if (is.vector(tt)) {
73-
tt[5] = ifelse(tt[5] == 0, 999999, tt[5])
74-
}
75-
return(tt)
29+
ProductData <- ProductData[!is.na(ProductData$ProductName), ]
30+
31+
# Stem order from XML (if available), otherwise from StemProfile
32+
stems_xml <- XMLNode[["Machine"]][names(XML::xmlSApply(XMLNode[["Machine"]], XML::xmlAttrs)) == "Stem"]
33+
if (length(stems_xml) > 0) {
34+
stem_keys <- vapply(
35+
stems_xml,
36+
function(s) as.integer(XML::xmlValue(s[["StemKey"]])),
37+
integer(1)
38+
)
39+
} else {
40+
stem_keys <- unique(as.integer(StemProfile$StemKey))
7641
}
7742

78-
stems <- XMLNode[["Machine"]][names(xmlSApply(XMLNode[["Machine"]],
79-
xmlAttrs)) == "Stem"]
80-
pb <- txtProgressBar(min = 0, max = length(stems), style = 3,
81-
width = 50, char = "=")
82-
ProductData <- ProductData[!is.na(ProductData$ProductName),]
83-
result <- list()
84-
for (i in 1:length(stems)) {
85-
StemKey <- SK <- as.integer(xmlValue(stems[[i]][["StemKey"]]))
86-
stem <- StemProfile[StemProfile$StemKey == SK, ]
87-
if (nrow(stem) > 0) {
88-
diameterPosition <- as.numeric(stem$diameterPosition)
89-
DiameterValue <- as.numeric(stem$DiameterValue)
90-
StemGrade <- as.numeric(stem$StemGrade)
91-
SpeciesGroupKey <- unique(stem$SpeciesGroupKey)
92-
ProductKeys <- ProductData$ProductKey
93-
LengthClassLowerLimit <- as.numeric(ProductData$LengthClassLowerLimit)
94-
LengthClassMAX <- as.numeric(ProductData$LengthClassMAX)
95-
DiameterClassLowerLimit <- as.numeric(ProductData$DiameterClassLowerLimit)
96-
DiameterClassMAX <- as.numeric(ProductData$DiameterClassMAX)
97-
VolumeDiameterCategory <- ProductData$VolumeDiameterCategory
98-
DiameterTopPositions <- ProductData$DiameterTopPosition
99-
DBH <- xmlValue(stems[[i]][["SingleTreeProcessedStem"]][["DBH"]]) %>%
100-
as.numeric()
43+
StemProfileDT <- data.table::as.data.table(StemProfile)
44+
stems_list <- split(StemProfileDT, StemProfileDT$StemKey)
10145

102-
SeqStart = round_any(min(LengthClassLowerLimit[LengthClassLowerLimit >
103-
0]), 10)
104-
SeqStop = ifelse(max(LengthClassMAX) < max(diameterPosition),
105-
max(LengthClassMAX), max(diameterPosition))
106-
DiameterTopPositions = ProductData$DiameterTopPositions
107-
bult = seq(10, 100, 10)
108-
res = data.table(StartPos = -1, StopPos = 0, Top_ub = NA,
109-
LogLength = NA, ProductKey = NA, Volume = 0,
110-
Value = 0, CumulativeValue = 0)
111-
if (SeqStart < SeqStop) {
112-
SeqAsp = seq(SeqStart, SeqStop, 10)
113-
StartPos = 0
114-
while (StartPos < max(diameterPosition) - min(LengthClassLowerLimit[LengthClassLowerLimit >
115-
0])) {
116-
StartPos = sort(res$StopPos[!res$StopPos %in%
117-
res$StartPos])[1]
46+
pb <- utils::txtProgressBar(min = 0, max = length(stem_keys), style = 3, width = 50, char = "=")
47+
on.exit(close(pb), add = TRUE)
11848

119-
if (StartPos == 0) {
120-
StopPos = StartPos + c(bult, SeqAsp)
121-
} else {
122-
StopPos = StartPos + SeqAsp
123-
}
124-
StopPos = StopPos[StopPos <= max(diameterPosition) &
125-
StopPos > 0]
126-
if (length(StopPos) < 1) {
127-
break
128-
}
129-
LogLength = StopPos - StartPos
130-
rotdiam = DiameterValue[which(near(diameterPosition,
131-
StartPos))]
132-
idxstart = as.numeric(which(near(diameterPosition,
133-
StartPos)))
134-
idxstop = as.numeric(match(as.character(StopPos),
135-
as.character(diameterPosition)))
136-
grd = lapply(idxstop, grdFinder)
137-
SGPK = ProductData$ProductKey[ProductData$SpeciesGroupKey ==
138-
SpeciesGroupKey[1]]
139-
m = data.table(idxstart, idxstop, StartPos,
140-
StopPos, LogLength, rotdiam)
141-
m = m[m$StopPos <= max(diameterPosition), ]
142-
SGKG = PermittedGrades[as.character(SGPK)]
143-
asos = lapply(1:length(grd), asoFinder)
144-
lapply(1:length(grd), function(x) {
145-
names(SGKG)[which(colSums(matrix(sapply(SGKG,
146-
FUN = function(X) all(grd[[x]] %in% X)),
147-
ncol = length(SGKG))) > 0)]
148-
})
149-
m$Price = 0
150-
r = rep(idxstop, len = sum(lengths(asos)))
151-
r = r[order(r)]
152-
tab = data.table(idxstop = r, ProductKey = unlist(asos))
153-
idx = match(as.character(tab$ProductKey), as.character(ProductData$ProductKey))
154-
tab = cbind(tab, data.table(DiameterUnderBark = ProductData$DiameterUnderBark[idx],
155-
LengthClassLowerLimit = ProductData$LengthClassLowerLimit[idx],
156-
LengthClassMAX = ProductData$LengthClassMAX[idx],
157-
DiameterClassLowerLimit = ProductData$DiameterClassLowerLimit[idx],
158-
DiameterClassMAX = ProductData$DiameterClassMAX[idx],
159-
VolumeDiameterAdjustment = ProductData$VolumeDiameterAdjustment[idx],
160-
VolumeDiameterCategory = ProductData$VolumeDiameterCategory[idx],
161-
VolumeLengthCategory = ProductData$VolumeLengthCategory[idx],
162-
DiameterTopPosition = as.numeric(ProductData$DiameterTopPositions[idx])))
163-
tab = merge(m, tab, "idxstop", allow.cartesian = TRUE)
164-
tab$StopPosAdj = round((tab$StopPos - tab$DiameterTopPosition)/10) *
165-
10
166-
tab$Top_ob = DiameterValue[match(as.character(tab$StopPosAdj),
167-
as.character(diameterPosition))]
168-
tab$Top_ub = BarkFunction(tab$Top_ob, SpeciesGroupKey,
169-
SpeciesGroupDefinition, Top_ob = Top_ob,
170-
DBH = DBH, LogLength = LogLength)
171-
tab$topdiam = ifelse(tab$DiameterUnderBark,
172-
tab$Top_ub, tab$Top_ob)
173-
tab = tab[tab$LogLength >= tab$LengthClassLowerLimit,]
174-
tab = tab[tab$LogLength <= tab$LengthClassMAX,]
175-
tab = tab[tab$topdiam > tab$DiameterClassLowerLimit,]
176-
tab = tab[tab$rotdiam < tab$DiameterClassMAX,]
177-
if (nrow(tab) > 0){
178-
commercial = tab[tab$ProductKey != "999999",]
179-
if (nrow(commercial) > 0) {
180-
lis = PriceMatrices[commercial$ProductKey]
181-
rownames = lapply(lis, rownames)
182-
colnames = lapply(lis, colnames)
183-
row = sapply(1:length(commercial$LogLength),
184-
rowFinder)
185-
col = sapply(1:length(commercial$topdiam),
186-
colFinder)
187-
tab$Price[tab$ProductKey != "999999"] = sapply(1:length(lis),
188-
priceFinder)
189-
}
190-
tab$idxstop[tab$VolumeLengthCategory == "Rounded downwards to nearest dm-module"] = match(as.character(round_any((tab$StopPos[tab$VolumeLengthCategory ==
191-
"Rounded downwards to nearest dm-module"]),
192-
10, f = floor)), as.character(diameterPosition))
193-
WithLengthClass = tab[tab$VolumeLengthCategory ==
194-
"Length as defined in LengthClasses" &
195-
tab$ProductKey != "999999", ]
196-
lis = LengthClasses[WithLengthClass$ProductKey]
197-
if (nrow(WithLengthClass) > 0) {
198-
l = 1
199-
for (l in 1:nrow(WithLengthClass)) {
200-
LengthClass = LengthClasses[[WithLengthClass$ProductKey[l]]]
201-
WithLengthClass$LogLength[l] = round_any(LengthClass[max(which(WithLengthClass$LogLength[l] >=
202-
LengthClass))], 10, f = ceiling)
203-
WithLengthClass$StopPos[l] = WithLengthClass$StartPos[l] +
204-
WithLengthClass$LogLength[l]
205-
WithLengthClass$idxstop[l] = which(diameterPosition ==
206-
paste(round_any(WithLengthClass$StopPos[l],
207-
10, f = ceiling)))
208-
}
209-
tab$LogLength[tab$VolumeLengthCategory ==
210-
"Length as defined in LengthClasses"] = WithLengthClass$LogLength
211-
tab$StopPos[tab$VolumeLengthCategory ==
212-
"Length as defined in LengthClasses"] = WithLengthClass$StopPos
213-
tab$idxstop[tab$VolumeLengthCategory ==
214-
"Length as defined in LengthClasses"] = WithLengthClass$idxstop
215-
}
216-
vec = seqVectozied(from = tab$idxstart, to = tab$idxstop,
217-
by = 1)
218-
DV = sapply(1:length(vec), DiameterValueFinder)
219-
idx = tab$VolumeDiameterAdjustment == "Measured diameter rounded downwards to cm"
220-
if (sum(idx) > 0) {
221-
if (sum(idx) > 1) {
222-
DV[which(idx > 0)] = lapply(1:sum(idx),
223-
function(x) {
224-
as.vector(sapply(1:length(DV[which(idx >
225-
0)[x]]), Rounder))
226-
})
227-
}
228-
else {
229-
DV[idx] = sapply(1:length(DV[idx]), Rounder)
230-
}
231-
}
232-
idx = tab$DiameterUnderBark == T
233-
if (sum(idx) > 0) {
234-
DV[idx] = sapply(1:length(DV[idx]), BarkFinder)
235-
}
236-
RV = relist(unlist(DV)/2, skeleton = as.relistable(DV))
237-
if (!is.list(RV)) {
238-
RV = list(RV)
239-
}
240-
tab$Volume = -1
241-
idx = tab$VolumeDiameterCategory == "All diameters (solid volume)"
242-
x = 2
243-
tab$Volume[idx] = sapply(1:length(RV), function(x) sum(pi *
244-
(unlist(RV[x])^2) * 10)/1e+08)[idx]
245-
idx = tab$VolumeDiameterCategory == "Calculated Norwegian mid"
246-
Dmid = tab$Top_ub + (tab$LogLength/2 * 0.1) +
247-
0.5
248-
tab$Volume[idx] = ((((Dmid/100) * (Dmid/100)) *
249-
pi/4 * (tab$LogLength/10)) * 0.001)[idx]
250-
idx = tab$VolumeDiameterCategory == "Top" &
251-
tab$DiameterUnderBark == T
252-
r1 = tab$Top_ub/2
253-
r2 = (tab$Top_ub + tab$LogLength * 0.01)/2
254-
tab$Volume[idx] = (((1/3) * pi * (r1^2 +
255-
r2^2 + (r1 * r2)) * tab$LogLength)/1e+08)[idx]
256-
idx = tab$VolumeDiameterCategory == "Top" &
257-
tab$DiameterUnderBark == F
258-
r1 = tab$Top_ob/2
259-
r2 = (tab$Top_ob + tab$LogLength * 0.01)/2
260-
tab$Volume[idx] = (((1/3) * pi * (r1^2 +
261-
r2^2 + (r1 * r2)) * tab$LogLength)/1e+08)[idx]
262-
tab$Value = tab$Volume * tab$Price
263-
head(tab)
264-
m = tab[, c("StartPos", "StopPos", "Top_ub",
265-
"LogLength", "ProductKey", "Volume", "Value")]
266-
idx=paste(res$StopPos)==paste(StartPos)
267-
sub = res[idx,]
49+
result <- vector("list", length(stem_keys))
26850

269-
CumulativeValue = ifelse(nrow(sub)>0&any(!is.na(sub$CumulativeValue)),
270-
max(sub$CumulativeValue,na.rm = T),
271-
0)
272-
m$CumulativeValue = m$Value + CumulativeValue
273-
}else {
274-
m = data.table(StartPos = StartPos, StopPos = StopPos,
275-
Top_ub = NA, LogLength = NA, ProductKey = NA,
276-
Volume = NA, Value = NA, CumulativeValue = NA)
277-
}
278-
res = rbindlist(list(res, m))
51+
for (i in seq_along(stem_keys)) {
52+
SK <- stem_keys[i]
53+
stem <- stems_list[[as.character(SK)]]
54+
55+
if (!is.null(stem) && nrow(stem) > 0) {
56+
sol <- tryCatch(
57+
buckStem(
58+
stem = stem,
59+
ProductData = ProductData,
60+
PriceMatrices = PriceMatrices,
61+
PermittedGrades = PermittedGrades,
62+
StemKey = as.integer(SK)
63+
),
64+
error = function(e) {
65+
warning(sprintf("buckStem failed for StemKey %s: %s", SK, conditionMessage(e)), call. = FALSE)
66+
NULL
27967
}
280-
}
281-
res = res[!is.na(res$LogLength),]
282-
tt = res[which.max(res$CumulativeValue),]
283-
if (nrow(tt) == 1){
284-
res = trackTrace(res, tt)
285-
}else{
286-
res = data.table(StartPos = NA, StopPos = NA,
287-
Top_ub = NA, LogLength = 1, ProductKey = NA,
288-
Volume = NA, Value = 0, CumulativeValue = 0)
289-
}
290-
res = cbind(1:nrow(res), res)
291-
colnames(res)[1] = "LogKey"
292-
out <- cbind(rep(StemKey, nrow(res)), res)
293-
colnames(out)[1] <- c("StemKey")
294-
result[[i]] <- out
68+
)
69+
result[[i]] <- sol
70+
} else {
71+
result[[i]] <- NULL
29572
}
296-
setTxtProgressBar(pb, i)
297-
print(i)
73+
74+
utils::setTxtProgressBar(pb, i)
29875
}
299-
result <- do.call(rbind.data.frame, result)
300-
close(pb)
301-
return(result)
302-
}
30376

77+
out <- data.table::rbindlist(result, fill = TRUE)
78+
return(out)
79+
}

0 commit comments

Comments
 (0)