~/Google Drive/projects/bbnj/software/gurobi_optimizer/gurobi8.1.0_mac64.pkg
By default, the installer will place the Gurobi 8.1.0 files in /Library/gurobi810/mac64
grbgetkey `cat ~/private/gurobi_key`
/Users/bbest/gurobi.lic
/Library/gurobi810/mac64/examples
and described in Example Code and Optimization Model Overview - Gurobigurobi_cl /Library/gurobi810/mac64/examples/data/coins.lp
Install R interface, per Solving the model using the Gurobi command-line interface: Gurobi Optimizer Example Tour: Gurobi Optimizer Quick Start Guide - Mac OS: R Interface: Installing the R Package:
tgz <- list.files('/Library/gurobi810/mac64/R', '.*tgz$', full.names = T)
install.packages(tgz, repos=NULL)
install.packages("slam", repos = "https://cloud.r-project.org")
Try out gurobi in R:
suppressPackageStartupMessages({
library('gurobi')
})
model <- list()
model$A <- matrix(c(1,2,3,1,1,0), nrow=2, ncol=3, byrow=T)
model$obj <- c(1,1,2)
model$modelsense <- 'max'
model$rhs <- c(4,1)
model$sense <- c('<', '>')
model$vtype <- 'B'
params <- list(OutputFlag=0)
#result <- gurobi(model, params)
result <- gurobi(model, list())
## Optimize a model with 2 rows, 3 columns and 5 nonzeros
## Variable types: 0 continuous, 3 integer (3 binary)
## Coefficient statistics:
## Matrix range [1e+00, 3e+00]
## Objective range [1e+00, 2e+00]
## Bounds range [0e+00, 0e+00]
## RHS range [1e+00, 4e+00]
## Found heuristic solution: objective 2.0000000
## Presolve removed 2 rows and 3 columns
## Presolve time: 0.00s
## Presolve: All rows and columns removed
##
## Explored 0 nodes (0 simplex iterations) in 0.00 seconds
## Thread count was 1 (of 12 available processors)
##
## Solution count 2: 3 2
##
## Optimal solution found (tolerance 1.00e-04)
## Best objective 3.000000000000e+00, best bound 3.000000000000e+00, gap 0.0000%
print('Solution:')
## [1] "Solution:"
print(result$objval)
## [1] 3
print(result$x)
## [1] 1 0 1
Solving a prioritzr problem with Gurobi
suppressPackageStartupMessages({
library(prioritizr) # install.packages("prioritizr")
})
# formulate the problem
p <- problem(sim_pu_raster, sim_features) %>%
add_min_set_objective() %>%
add_relative_targets(0.1) %>%
add_gurobi_solver()
# solve the problem
s <- solve(p)
## Optimize a model with 5 rows, 90 columns and 450 nonzeros
## Variable types: 0 continuous, 90 integer (90 binary)
## Coefficient statistics:
## Matrix range [2e-01, 9e-01]
## Objective range [2e+02, 2e+02]
## Bounds range [1e+00, 1e+00]
## RHS range [3e+00, 8e+00]
## Found heuristic solution: objective 2337.9617505
## Presolve time: 0.00s
## Presolved: 5 rows, 90 columns, 450 nonzeros
## Variable types: 0 continuous, 90 integer (90 binary)
## Presolved: 5 rows, 90 columns, 450 nonzeros
##
##
## Root relaxation: objective 1.931582e+03, 12 iterations, 0.00 seconds
##
## Nodes | Current Node | Objective Bounds | Work
## Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
##
## 0 0 1931.58191 0 4 2337.96175 1931.58191 17.4% - 0s
## H 0 0 1985.6818841 1931.58191 2.72% - 0s
##
## Explored 1 nodes (12 simplex iterations) in 0.00 seconds
## Thread count was 1 (of 12 available processors)
##
## Solution count 2: 1985.68 2337.96
##
## Optimal solution found (tolerance 1.00e-01)
## Best objective 1.985681884076e+03, best bound 1.931581908865e+03, gap 2.7245%
# plot solution
plot(s, col = c("grey90", "darkgreen"), main = "Solution",
xlim = c(-0.1, 1.1), ylim = c(-0.1, 1.1))
suppressPackageStartupMessages({
library(tidyverse)
library(glue)
library(raster)
library(sf)
library(leaflet)
library(RColorBrewer)
library(rasterVis)
library(prioritizr) # install.packages("prioritizr")
#library(bbnj)
devtools::load_all()
})
dir_data <- "~/Gdrive Ecoquants/projects/bbnj/data/derived"
cell_res <- 0.5 # n= 113,401
p_boundary_shp <- file.path(
dir_data, glue("boundary/high_seas.shp"))
r_pu_tif <- file.path(
dir_data, glue("boundary/high_seas_cellid_{cell_res}dd.tif"))
# get planning unit raster
r_pu <- raster(r_pu_tif)
pal <- colorRampPalette(brewer.pal(11, "Spectral"))
cols <- rev(pal(255))
plot(r_pu, col = cols, main="pu")
# get biodiversity features
bio_tifs <- list.files(
file.path(dir_data, "biodiversity"),
"spp_.*_be_0\\.5dd\\.tif", full.names = T)
bio_vars <- tibble(
tif = bio_tifs,
var = str_replace(basename(tif), "spp_(.*)_be_0\\.5dd\\.tif", "\\1")) %>%
mutate(
var = recode(var, `tunas.&.billfishes`="tunas.billfishes")) %>%
pull(var)
s_bio <- stack(bio_tifs)
names(s_bio) <- bio_vars
#plot(s_bio, col = cols, main="spp")
levelplot(s_bio, main="spp")
# formulate the problem
# TODO: check that s_bio being continuous is ok, and not binary
p <- problem(r_pu, s_bio) %>%
add_min_set_objective() %>%
add_relative_targets(0.1) %>%
add_gurobi_solver()
# solve the problem
r_sol <- solve(p)
# plot solution
plot(r_sol, col = c("grey90", "darkgreen"), main = "Solution")
# devtools::load_all()
qmap_r(r_sol, "Solution", na0=T)
TODO: show pretty stack like
prioritizr: Systematic Conservation Prioritization in R • prioritizr
# load prioritizr package
library(prioritizr)
# load raster planning unit data
data(sim_pu_raster)
# print description of the data
print(sim_pu_raster)
## class : RasterLayer
## dimensions : 10, 10, 100 (nrow, ncol, ncell)
## resolution : 0.1, 0.1 (x, y)
## extent : 0, 1, 0, 1 (xmin, xmax, ymin, ymax)
## coord. ref. : NA
## data source : in memory
## names : layer
## values : 190.1328, 215.8638 (min, max)
# plot the data
plot(sim_pu_raster)
# load polygon planning unit data
data(sim_pu_polygons)
# print first six rows of attribute table
head(sim_pu_polygons@data)
## cost locked_in locked_out
## 1 215.8638 FALSE FALSE
## 2 212.7823 FALSE FALSE
## 3 207.4962 FALSE FALSE
## 4 208.9322 FALSE TRUE
## 5 214.0419 FALSE FALSE
## 6 213.7636 FALSE FALSE
# plot the planning units
spplot(sim_pu_polygons, zcol = "cost")
# specify file path for planning unit data
pu_path <- system.file("extdata/input/pu.dat", package = "prioritizr")
# load in the tabular planning unit data
# note that we use the data.table::fread function, as opposed to the read.csv
# function, because it is much faster
pu_dat <- data.table::fread(pu_path, data.table = FALSE)
# preview first six rows of the tabular planning unit data
# note that it has some extra columns other than id and cost as per the
# Marxan format
head(pu_dat)
## id cost status xloc yloc
## 1 3 0.0 0 1116623 -4493479
## 2 30 752727.5 3 1110623 -4496943
## 3 56 3734907.5 0 1092623 -4500408
## 4 58 1695902.1 0 1116623 -4500408
## 5 84 3422025.6 0 1098623 -4503872
## 6 85 17890758.4 0 1110623 -4503872
# load feature data
data(sim_features)
# plot the distribution of suitable habitat for each feature
plot(sim_features, main = paste("Feature", seq_len(nlayers(sim_features))),
nr = 2, box = FALSE, axes = FALSE)
# create problem
p1 <- problem(sim_pu_raster, sim_features)
# print problem
print(p1)
## Conservation Problem
## planning units: RasterLayer (90 units)
## cost: min: 190.13276, max: 215.86384
## features: layer.1, layer.2, layer.3, ... (5 features)
## objective: none
## targets: none
## decisions: default
## constraints: <none>
## penalties: <none>
## portfolio: default
## solver: default
# print number of planning units
number_of_planning_units(p1)
## [1] 90
# print number of features
number_of_features(p1)
## [1] 5
# create problem with spatial vector data
# note that we have to specify which column in the attribute table contains
# the cost data
p2 <- problem(sim_pu_polygons, sim_features, cost_column = "cost")
# print problem
print(p2)
## Conservation Problem
## planning units: SpatialPolygonsDataFrame (90 units)
## cost: min: 190.13276, max: 215.86384
## features: layer.1, layer.2, layer.3, ... (5 features)
## objective: none
## targets: none
## decisions: default
## constraints: <none>
## penalties: <none>
## portfolio: default
## solver: default
# set file path for feature data
spec_path <- system.file("extdata/input/spec.dat", package = "prioritizr")
# load in feature data
spec_dat <- data.table::fread(spec_path, data.table = FALSE)
# print first six rows of the data
# note that it contains extra columns
head(spec_dat)
## id prop spf name
## 1 10 0.3 1 bird1
## 2 11 0.3 1 nvis2
## 3 12 0.3 1 nvis8
## 4 13 0.3 1 nvis9
## 5 14 0.3 1 nvis14
## 6 15 0.3 1 nvis20
# set file path for planning unit vs. feature data
puvspr_path <- system.file("extdata/input/puvspr.dat", package = "prioritizr")
# load in planning unit vs feature data
puvspr_dat <- data.table::fread(puvspr_path, data.table = FALSE)
# print first six rows of the data
head(puvspr_dat)
## species pu amount
## 1 26 56 1203448.84
## 2 26 58 451670.10
## 3 26 84 680473.75
## 4 26 85 97356.24
## 5 26 86 78034.76
## 6 26 111 4783274.17
# create problem
p3 <- problem(pu_dat, spec_dat, cost_column = "cost", rij = puvspr_dat)
# print problem
print(p3)
## Conservation Problem
## planning units: data.frame (1751 units)
## cost: min: 0, max: 41569219.38232
## features: bird1, nvis2, nvis8, ... (17 features)
## objective: none
## targets: none
## decisions: default
## constraints: <none>
## penalties: <none>
## portfolio: default
## solver: default
# create a new problem that has the minimum set objective
p3 <- problem(sim_pu_raster, sim_features) %>%
add_min_set_objective()
# print the problem
print(p3)
## Conservation Problem
## planning units: RasterLayer (90 units)
## cost: min: 190.13276, max: 215.86384
## features: layer.1, layer.2, layer.3, ... (5 features)
## objective: Minimum set objective
## targets: none
## decisions: default
## constraints: <none>
## penalties: <none>
## portfolio: default
## solver: default
# create a new problem that has the maximum coverage objective and a budget
# of 5000
p4 <- problem(sim_pu_raster, sim_features) %>%
add_max_cover_objective(5000)
# print the problem
print(p4)
## Conservation Problem
## planning units: RasterLayer (90 units)
## cost: min: 190.13276, max: 215.86384
## features: layer.1, layer.2, layer.3, ... (5 features)
## objective: Maximum coverage objective [budget (5000)]
## targets: none
## decisions: default
## constraints: <none>
## penalties: <none>
## portfolio: default
## solver: default
# create a new problem that has the maximum features objective and a budget
# of 5000
p5 <- problem(sim_pu_raster, sim_features) %>%
add_max_features_objective(budget = 5000)
# print the problem
print(p5)
## Conservation Problem
## planning units: RasterLayer (90 units)
## cost: min: 190.13276, max: 215.86384
## features: layer.1, layer.2, layer.3, ... (5 features)
## objective: Maximum representation objective [budget (5000)]
## targets: none
## decisions: default
## constraints: <none>
## penalties: <none>
## portfolio: default
## solver: default
library(ape) # ?plot.phylo
##
## Attaching package: 'ape'
## The following objects are masked from 'package:raster':
##
## rotate, zoom
# load simulated phylogeny data
data(sim_phylogeny)
plot(sim_phylogeny)
# create a new problem that has the maximum phylogenetic representation
# objective and a budget of 5000
p6 <- problem(sim_pu_raster, sim_features) %>%
add_max_phylo_objective(budget = 5000, tree = sim_phylogeny)
# print the problem
print(p6)
## Conservation Problem
## planning units: RasterLayer (90 units)
## cost: min: 190.13276, max: 215.86384
## features: layer.1, layer.2, layer.3, ... (5 features)
## objective: Phylogenetic representation objective [budget (5000)]
## targets: none
## decisions: default
## constraints: <none>
## penalties: <none>
## portfolio: default
## solver: default
Could this be considered coarse to fine objectives? Habitats vs Species and within Species actual phylogenetic relationships?
# create a new problem that has the maximum utility objective and a budget
# of 5000
p7 <- problem(sim_pu_raster, sim_features) %>%
add_max_utility_objective(budget = 5000)
# print the problem
print(p7)
## Conservation Problem
## planning units: RasterLayer (90 units)
## cost: min: 190.13276, max: 215.86384
## features: layer.1, layer.2, layer.3, ... (5 features)
## objective: Maximum utility objective [budget (5000)]
## targets: none
## decisions: default
## constraints: <none>
## penalties: <none>
## portfolio: default
## solver: default