MineLib in AMPL and amplpy#

Description: A notebook that works with the open-pit mining problems defined in the MineLib library and shows how to read and parse the instances using Python and amplpy

Tags: educational, mining, amplpy, gurobi, minelib

Notebook author: Eduardo Salazar <eduardo@ampl.com>

Model author: Daniel Espinoza (Universidad de Chile), Marcos Goycoolea (Universidad Adolfo Ibañez), Eduardo Moreno (Universidad Adolfo Ibañez), Alexandra Newman (Colorado School of Mines).

References:

Help on Minelib: https://mansci-web.uai.cl/minelib/Help.xhtml
MineLib: A Library of Open Pit Mining Problems. Daniel Espinoza, Marcos Goycoolea, Eduardo Moreno, Alexandra Newman https://mgoycool.github.io/papers/13espinoza_aor.pdf.

# Install dependencies
%pip install -q amplpy requests

# Google Colab & Kaggle integration
from amplpy import AMPL, ampl_notebook

ampl = ampl_notebook(
    modules=["gurobi"],  # modules to install
    license_uuid="default",  # license to use
)  # instantiate AMPL object and register magics

1. Introduction to Open Pit Mining #

Open pit mining is a surface mining technique used to extract minerals from the earth by removing layers of rock and soil. The extracted material is discretized into blocks, each containing estimated amounts of ore (valuable mineral) and waste. The key decisions involve determining which blocks to extract and when, subject to various operational and geometric constraints.

Key Concepts#

Block Model: The mineral deposit is represented as a three-dimensional grid of blocks. Each block has associated attributes such as position $(x, y, z)$, tonnage, ore grade, and estimated economic value.

Precedence Constraints: Due to geotechnical requirements for pit wall stability and equipment access, blocks cannot be extracted in arbitrary order. If block $a$ lies above block $b$, then $a$ must be extracted before $b$. These relationships form a directed acyclic graph where edges represent “must be extracted before” relationships.

Pit Slope: The angle at which the pit walls are maintained to prevent collapse. This angle determines which blocks must be removed to access deeper blocks.

Cutoff Grade: The minimum ore concentration that makes a block economically viable for processing. Blocks below this threshold are typically sent to waste dumps.

Why Optimization?#

Mining operations represent massive capital investments spanning decades. A typical open pit mine may contain millions of blocks, and the extraction sequence significantly impacts profitability due to the time value of money. Early extraction of high-value ore maximizes net present value (NPV). Mathematical optimization helps mine planners make informed decisions about pit design and production scheduling.

2. The MineLib Library#

MineLib was introduced by Espinoza, Goycoolea, Moreno, and Newman (2012) as a standardized library of open pit mining problem instances. Similar to how MIPLIB serves the mixed-integer programming community, MineLib provides benchmark data for evaluating and comparing optimization algorithms.

Available Instances#

Instance	Blocks	Precedences	Periods	Description
Newman1	1,060	3,922	6	Small academic dataset
Zuck_small	9,400	145,640	20	Fictitious mine
D	14,153	219,778	12	North American copper deposit
Zuck_medium	29,277	1,271,207	15	Fictitious mine
P4HD	40,947	738,609	10	North American gold/copper mine
Marvin	53,271	650,631	20	Whittle software test mine
W23	74,260	764,786	12	North American gold mine
Zuck_large	96,821	1,053,105	30	Fictitious mine
SM2	99,014	96,642	30	Based on Brazilian nickel mine
McLaughlin_limit	112,687	3,035,483	15	California gold mine (final pit)
McLaughlin	2,140,342	73,143,770	20	California gold mine (full)

Problem Types#

MineLib supports three problem formulations of increasing complexity:

UPIT (Ultimate Pit Limit): Determines the optimal pit boundary ignoring time.
CPIT (Constrained Pit Limit): Adds time periods and resource constraints with fixed cutoff grades.
PCPSP (Precedence Constrained Production Scheduling): Variable cutoff grades with multiple destinations and general side constraints.

3. Mathematical Formulations#

We now present the mathematical formulations for the three problem types. The following notation is used throughout:

3.1 Notation#

Indices and Sets:

Symbol	Description
$t \in \mathcal{T}$	Set of time periods $t$ in the horizon
$b \in \mathcal{B}$	Set of blocks $b$
$b' \in \mathcal{B}_b$	Set of blocks $b'$ that are predecessor blocks for block $b$
$r \in \mathcal{R}$	Set of operational resource types $r$
$d \in \mathcal{D}$	Set of destinations $d$

Parameters:

Symbol	Description
$p_b$ ($\hat{p}_{bt}$, $\check{p}_{bd}$, $\tilde{p}_{bdt}$)	Profit obtained from extracting (and processing) block $b$ (at time period $t$ and/or sending it to destination $d$) ($)
$\alpha$	Discount rate used in computing the objective function (profit) coefficients
$q_{br}$ ($\hat{q}_{brd}$)	Amount of operational resource $r$ used to extract and, if applicable, process block $b$ (when sent to destination $d$) (tons)
$\underline{R}_{rt}$	Minimum availability of operational resource $r$ in time period $t$ (tons)
$\overline{R}_{rt}$	Maximum availability of operational resource $r$ in time period $t$ (tons)
$A$	Arbitrary constraint coefficients on general side constraints
$\underline{a}, \overline{a}$	Arbitrary lower and upper bounds, respectively, on general side constraints

Variables:

Symbol	Description
$\hat{x}_b$	1 if block $b$ is in the final pit design, 0 otherwise
$x_{bt}$	1 if we extract block $b$ in time period $t$, 0 otherwise
$y_{bdt}$	The amount of block $b$ sent to destination $d$ in time period $t$ (%)

3.2 The Ultimate Pit Problem (UPIT)#

The simplest model we consider is known as the ultimate pit limit problem, or the maximum-weight closure problem. The problem entails determining only the envelope of profitable blocks within the orebody and, hence, there is no temporal dimension and there are no operational resource constraints. The constraint set consists merely of precedences between blocks; the corresponding matrix of left-hand-side coefficients is totally unimodular, rendering this problem a network flow problem.

Formulation:

\[ \text{(UPIT)} \quad \max \sum_{b \in \mathcal{B}} p_b \, \hat{x}_b \]

subject to:

\[ \hat{x}_b \leq \hat{x}_{b'} \quad \forall b \in \mathcal{B}, \; \forall b' \in \mathcal{B}_b \tag{1} \]

\[ \hat{x}_b \in \{0, 1\} \quad \forall b \in \mathcal{B} \tag{2} \]

Notes:

The objective maximizes the undiscounted value of all extracted blocks.
Constraints (1) ensure that each block is extracted only if its predecessor blocks are extracted.
Variables need not be restricted to be binary because of the total unimodularity of the constraint matrix.
The Lerchs-Grossmann algorithm (1965) and successors provide efficient specialized solutions.

3.3 The Constrained Pit Limit Problem (CPIT)#

The constrained pit limit problem generalizes the ultimate pit limit problem by introducing a time dimension, and associated constraints, into the model. The underlying assumption is that a block can be mined in its entirety in a single time period. Not only are precedence constraints considered, but per-period operational resource restrictions are present as well.

Formulation:

\[ \text{(CPIT)} \quad \max \sum_{b \in \mathcal{B}} \sum_{t \in \mathcal{T}} \hat{p}_{bt} \, x_{bt} \]

subject to:

\[ \sum_{\tau \leq t} x_{b\tau} \leq \sum_{\tau \leq t} x_{b'\tau} \quad \forall b \in \mathcal{B}, \; b' \in \mathcal{B}_b, \; t \in \mathcal{T} \tag{3} \]

\[ \sum_{t \in \mathcal{T}} x_{bt} \leq 1 \quad \forall b \in \mathcal{B} \tag{4} \]

\[ \underline{R}_{rt} \leq \sum_{b \in \mathcal{B}} q_{br} \, x_{bt} \leq \overline{R}_{rt} \quad \forall t \in \mathcal{T}, \; r \in \mathcal{R} \tag{5} \]

\[ x_{bt} \in \{0, 1\} \quad \forall b \in \mathcal{B}, \; t \in \mathcal{T} \tag{6} \]

where $\hat{p}_{bt} = \dfrac{p_b}{(1+\alpha)^t}$.

Constraint Interpretation:

Precedence (3): If block $b'$ is an immediate predecessor of block $b$, then $b'$ must be extracted in the same time period as or prior to $b$.
Single extraction (4): Each block can be extracted no more than once.
Resource bounds (5): Minimum and maximum operational resource constraints are satisfied each period.

Complexity: Unlike UPIT, CPIT is strongly NP-hard.

3.4 The Precedence Constrained Production Scheduling Problem (PCPSP)#

A generalization of CPIT determines whether a block, if extracted, is sent to the processing plant or to the waste dump. In addition to variable $x_{bt}$, we employ a second variable $y_{bdt}$, which equals the amount of block $b$ we send to destination $d$ in time period $t$. This enables variable cutoff grades where the destination is optimized rather than predetermined.

Formulation:

\[ \text{(PCPSP)} \quad \max \sum_{b \in \mathcal{B}} \sum_{d \in \mathcal{D}} \sum_{t \in \mathcal{T}} \tilde{p}_{bdt} \, y_{bdt} \]

subject to:

\[ \sum_{\tau \leq t} x_{b\tau} \leq \sum_{\tau \leq t} x_{b'\tau} \quad \forall b \in \mathcal{B}, \; b' \in \mathcal{B}_b, \; t \in \mathcal{T} \tag{7} \]

\[ x_{bt} = \sum_{d \in \mathcal{D}} y_{bdt} \quad \forall b \in \mathcal{B}, \; t \in \mathcal{T} \tag{8} \]

\[ \sum_{t \in \mathcal{T}} x_{bt} \leq 1 \quad \forall b \in \mathcal{B} \tag{9} \]

\[ \underline{R}_{rt} \leq \sum_{b \in \mathcal{B}} \sum_{d \in \mathcal{D}} \hat{q}_{brd} \, y_{bdt} \leq \overline{R}_{rt} \quad \forall r \in \mathcal{R}, \; t \in \mathcal{T} \tag{10} \]

\[ \underline{a} \leq A y \leq \overline{a} \tag{11} \]

\[ y_{bdt} \in [0, 1] \quad \forall b \in \mathcal{B}, \; d \in \mathcal{D}, \; t \in \mathcal{T} \tag{12} \]

\[ x_{bt} \in \{0, 1\} \quad \forall b \in \mathcal{B}, \; t \in \mathcal{T} \tag{13} \]

where $\tilde{p}_{bdt} = \dfrac{\check{p}_{bd}}{(1+\alpha)^t}$.

Constraint Interpretation:

Precedence (7): Enforces precedence requirements for all blocks and time periods.
Linking (8): If a block is not extracted, its contents cannot be sent to any destination; if extracted, the entirety must be sent somewhere.
Single extraction (9): A block can be extracted at most once over the horizon.
Resource bounds (10): No more operational resource than available is used for extraction purposes.
General side constraints (11): Can model blending, grade constraints, stockpiling, etc.

Example - Grade Constraint:

Letting $g_{bm}$ be the amount of mineral $m$ in block $b$, and $\underline{G}_m$, $\overline{G}_m$ be the minimum/maximum acceptable average grades:

\[ \underline{G}_m \sum_{b \in \mathcal{B}} y_{bdt} \leq \sum_{b \in \mathcal{B}} g_{bm} \, y_{bdt} \leq \overline{G}_m \sum_{b \in \mathcal{B}} y_{bdt} \quad \forall d \in \mathcal{D}, \; m \in \mathcal{M}, \; t \in \mathcal{T} \tag{14} \]

4. MineLib File Format#

MineLib instances are distributed in three file types:

4.1 Block-Model Descriptor File (`.blocks`)#

Contains geological data for each block:

<id> <x> <y> <z> [additional fields...]

id: Unique block identifier (0-indexed)
x, y, z: Block coordinates (z=0 is the bottom)
Additional fields may include tonnage, grade, etc.

4.2 Block-Precedence Descriptor File (`.prec`)#

Defines precedence relationships:

<block_id> <num_predecessors> <pred_1> <pred_2> ... <pred_n>

Example:

100 5 10 11 12 20 21

Block 100 has 5 predecessors: blocks 10, 11, 12, 20, and 21 must be extracted first.

4.3 Optimization-Model Descriptor File (`.upit`, `.cpit`, `.pcpsp`)#

Contains optimization parameters:

NAME: instance_name
TYPE: UPIT|CPIT|PCPSP
NBLOCKS: n
NPERIODS: T                    (CPIT/PCPSP only)
NDESTINATIONS: D               (PCPSP only)
NRESOURCE_SIDE_CONSTRAINTS: R  (CPIT/PCPSP only)
DISCOUNT_RATE: alpha           (CPIT/PCPSP only)

OBJECTIVE_FUNCTION:
<block_id> <value_dest_0> [<value_dest_1> ...]
...

RESOURCE_CONSTRAINT_LIMITS:
<resource> <period> <type> <bound1> [<bound2>]
...

RESOURCE_CONSTRAINT_COEFFICIENTS:
<block> [<dest>] <resource> <coefficient>
...
EOF

The constraint type can be:

L: Less than or equal ($\leq$ bound1)
G: Greater than or equal ($\geq$ bound1)
I: Interval (bound1 $\leq \cdot \leq$ bound2)

In the hidden code cells we define the MineLib parser tools for instances and solutions in Python.

For this notebook, we will be solving the smallest instance which is newman1.

from pathlib import Path
from typing import Dict, List, Optional, Tuple, Union, TextIO
from dataclasses import dataclass, field
import math
import requests
from io import StringIO
from contextlib import contextmanager

# @title


@dataclass
class Precedences:
    """Block precedence relationships."""

    num_blocks: int
    predecessors: Dict[int, List[int]] = field(default_factory=dict)

    def get_predecessors(self, block_id: int) -> List[int]:
        return self.predecessors.get(block_id, [])

    def num_precedence_arcs(self) -> int:
        return sum(len(preds) for preds in self.predecessors.values())


@dataclass
class ResourceLimits:
    """Resource constraint bounds by resource and time period."""

    lower_bounds: Dict[int, Dict[int, float]] = field(default_factory=dict)
    upper_bounds: Dict[int, Dict[int, float]] = field(default_factory=dict)


@dataclass
class UPITData:
    """Data for Ultimate Pit Limit Problem."""

    name: str
    num_blocks: int
    objective: Dict[int, float] = field(default_factory=dict)
    precedences: Optional[Precedences] = None


@dataclass
class CPITData:
    """Data for Constrained Pit Limit Problem."""

    name: str
    num_blocks: int
    num_periods: int
    num_resources: int
    discount_rate: float
    objective: Dict[int, float] = field(default_factory=dict)
    resource_limits: ResourceLimits = field(default_factory=ResourceLimits)
    resource_coefficients: Dict[int, Dict[int, float]] = field(default_factory=dict)
    precedences: Optional[Precedences] = None


@dataclass
class PCPSPData:
    """Data for Precedence Constrained Production Scheduling Problem."""

    name: str
    num_blocks: int
    num_periods: int
    num_destinations: int
    num_resources: int
    num_general_constraints: int
    discount_rate: float
    objective: Dict[int, Dict[int, float]] = field(default_factory=dict)
    resource_limits: ResourceLimits = field(default_factory=ResourceLimits)
    resource_coefficients: Dict[int, Dict[int, Dict[int, float]]] = field(
        default_factory=dict
    )
    general_coefficients: Dict[int, Dict[int, Dict[int, Dict[int, float]]]] = field(
        default_factory=dict
    )
    general_limits: Dict[int, Tuple[float, float]] = field(default_factory=dict)
    precedences: Optional[Precedences] = None

# @title


def parse_precedences(file_source: Union[str, Path, TextIO]) -> Precedences:
    """
    Parse a MineLib precedence file.

    Args:
        file_source: Either a file path (str/Path) or a file-like object (StringIO)
    """
    predecessors: Dict[int, List[int]] = {}
    num_blocks = 0

    # Handle both file paths and file-like objects
    if isinstance(file_source, (str, Path)):
        f = open(file_source, "r")
        should_close = True
    else:
        f = file_source
        should_close = False

    try:
        for line in f:
            line = line.strip()
            if not line or line.startswith("%"):
                continue

            parts = line.split()
            block_id = int(parts[0])
            num_preds = int(parts[1])
            preds = [int(parts[i]) for i in range(2, 2 + num_preds)]

            predecessors[block_id] = preds
            num_blocks = max(num_blocks, block_id + 1)
    finally:
        if should_close:
            f.close()

    return Precedences(num_blocks=num_blocks, predecessors=predecessors)


def _parse_float(value: str) -> float:
    """Parse float, handling infinity."""
    value = value.strip().lower()
    if value in ("infinity", "inf", "+infinity", "+inf"):
        return math.inf
    if value in ("-infinity", "-inf"):
        return -math.inf
    return float(value)


@contextmanager
def open_source(source: Union[str, Path, TextIO]):
    """Context manager that handles both file paths and file-like objects."""
    if isinstance(source, (str, Path)):
        f = open(source, "r")
        try:
            yield f
        finally:
            f.close()
    else:
        # Already a file-like object (StringIO), just yield it
        yield source


def parse_upit(
    opt_source: Union[str, Path, TextIO],
    prec_source: Optional[Union[str, Path, TextIO, Precedences]] = None,
) -> UPITData:
    """Parse a UPIT optimization file."""
    data = UPITData(name="", num_blocks=0)

    # Use context manager instead of Path()
    with open_source(opt_source) as f:
        for line in f:
            line = line.strip()
            if not line or line.startswith("%"):
                continue

            if ":" in line:
                key, value = line.split(":", 1)
                key, value = key.strip().upper(), value.strip()

                if key == "NAME":
                    data.name = value
                elif key == "NBLOCKS":
                    data.num_blocks = int(value)
                elif key == "OBJECTIVE_FUNCTION":
                    for _ in range(data.num_blocks):
                        obj_line = f.readline().strip()
                        while not obj_line or obj_line.startswith("%"):
                            obj_line = f.readline().strip()
                        parts = obj_line.split()
                        block_id = int(parts[0])
                        data.objective[block_id] = _parse_float(parts[1])

    # Handle precedences flexibly
    if prec_source is not None:
        if isinstance(prec_source, Precedences):
            data.precedences = prec_source
        else:
            data.precedences = parse_precedences(prec_source)

    return data


def parse_cpit(
    opt_source: Union[str, Path, TextIO],
    prec_source: Optional[Union[str, Path, TextIO, Precedences]] = None,
) -> CPITData:
    """Parse a CPIT optimization file."""
    data = CPITData(
        name="", num_blocks=0, num_periods=0, num_resources=0, discount_rate=0.0
    )

    with open_source(opt_source) as f:
        for line in f:
            line = line.strip()
            if not line or line.startswith("%"):
                continue

            if ":" in line:
                key, value = line.split(":", 1)
                key, value = key.strip().upper().replace(" ", "_"), value.strip()

                if key == "NAME":
                    data.name = value
                elif key == "NBLOCKS":
                    data.num_blocks = int(value)
                elif key == "NPERIODS":
                    data.num_periods = int(value)
                elif key in (
                    "NRESOURCE_SIDE_CONSTRAINTS",
                    "NRESOURCE_SIDE_CONSTRAINTS",
                ):
                    data.num_resources = int(value)
                elif key == "DISCOUNT_RATE":
                    data.discount_rate = _parse_float(value)
                elif key == "OBJECTIVE_FUNCTION":
                    for _ in range(data.num_blocks):
                        obj_line = f.readline().strip()
                        while not obj_line or obj_line.startswith("%"):
                            obj_line = f.readline().strip()
                        parts = obj_line.split()
                        data.objective[int(parts[0])] = _parse_float(parts[1])
                elif key == "RESOURCE_CONSTRAINT_LIMITS":
                    data.resource_limits.lower_bounds = {
                        r: {} for r in range(data.num_resources)
                    }
                    data.resource_limits.upper_bounds = {
                        r: {} for r in range(data.num_resources)
                    }
                    for _ in range(data.num_resources * data.num_periods):
                        lim_line = f.readline().strip()
                        while not lim_line or lim_line.startswith("%"):
                            lim_line = f.readline().strip()
                        parts = lim_line.split()
                        r, t, ctype = int(parts[0]), int(parts[1]), parts[2].upper()
                        if ctype == "L":
                            data.resource_limits.lower_bounds[r][t] = -math.inf
                            data.resource_limits.upper_bounds[r][t] = _parse_float(
                                parts[3]
                            )
                        elif ctype == "G":
                            data.resource_limits.lower_bounds[r][t] = _parse_float(
                                parts[3]
                            )
                            data.resource_limits.upper_bounds[r][t] = math.inf
                        elif ctype == "I":
                            data.resource_limits.lower_bounds[r][t] = _parse_float(
                                parts[3]
                            )
                            data.resource_limits.upper_bounds[r][t] = _parse_float(
                                parts[4]
                            )
                elif key == "RESOURCE_CONSTRAINT_COEFFICIENTS":
                    while True:
                        coef_line = f.readline()
                        if not coef_line:
                            break
                        coef_line = coef_line.strip()
                        if not coef_line or coef_line.startswith("%"):
                            continue
                        if coef_line.upper() == "EOF":
                            break
                        parts = coef_line.split()
                        if len(parts) < 3:
                            break
                        block_id, resource_id = int(parts[0]), int(parts[1])
                        coef_value = _parse_float(parts[2])
                        if block_id not in data.resource_coefficients:
                            data.resource_coefficients[block_id] = {}
                        data.resource_coefficients[block_id][resource_id] = coef_value

    # Handle precedences flexibly
    if prec_source is not None:
        if isinstance(prec_source, Precedences):
            data.precedences = prec_source
        else:
            data.precedences = parse_precedences(prec_source)

    return data

def parse_pcpsp(
    opt_source: Union[str, Path, TextIO],
    prec_source: Optional[Union[str, Path, TextIO, Precedences]] = None,
) -> PCPSPData:
    """Parse a PCPSP optimization file."""
    data = PCPSPData(
        name="",
        num_blocks=0,
        num_periods=0,
        num_destinations=0,
        num_resources=0,
        num_general_constraints=0,
        discount_rate=0.0,
    )

    with open_source(opt_source) as f:
        for line in f:
            line = line.strip()
            if not line or line.startswith("%"):
                continue

            if ":" in line:
                key, value = line.split(":", 1)
                key, value = key.strip().upper().replace(" ", "_"), value.strip()

                if key == "NAME":
                    data.name = value
                elif key == "NBLOCKS":
                    data.num_blocks = int(value)
                elif key == "NPERIODS":
                    data.num_periods = int(value)
                elif key == "NDESTINATIONS":
                    data.num_destinations = int(value)
                elif key in ("NRESOURCE_SIDE_CONSTRAINTS",):
                    data.num_resources = int(value)
                elif key in ("NGENERAL_SIDE_CONSTRAINTS",):
                    data.num_general_constraints = int(value)
                elif key == "DISCOUNT_RATE":
                    data.discount_rate = _parse_float(value)
                elif key == "OBJECTIVE_FUNCTION":
                    for _ in range(data.num_blocks):
                        obj_line = f.readline().strip()
                        while not obj_line or obj_line.startswith("%"):
                            obj_line = f.readline().strip()
                        parts = obj_line.split()
                        block_id = int(parts[0])
                        data.objective[block_id] = {}
                        for d in range(data.num_destinations):
                            data.objective[block_id][d] = _parse_float(parts[1 + d])
                elif key == "RESOURCE_CONSTRAINT_LIMITS":
                    data.resource_limits.lower_bounds = {
                        r: {} for r in range(data.num_resources)
                    }
                    data.resource_limits.upper_bounds = {
                        r: {} for r in range(data.num_resources)
                    }
                    for _ in range(data.num_resources * data.num_periods):
                        lim_line = f.readline().strip()
                        while not lim_line or lim_line.startswith("%"):
                            lim_line = f.readline().strip()
                        parts = lim_line.split()
                        r, t, ctype = int(parts[0]), int(parts[1]), parts[2].upper()
                        if ctype == "L":
                            data.resource_limits.lower_bounds[r][t] = -math.inf
                            data.resource_limits.upper_bounds[r][t] = _parse_float(
                                parts[3]
                            )
                        elif ctype == "G":
                            data.resource_limits.lower_bounds[r][t] = _parse_float(
                                parts[3]
                            )
                            data.resource_limits.upper_bounds[r][t] = math.inf
                        elif ctype == "I":
                            data.resource_limits.lower_bounds[r][t] = _parse_float(
                                parts[3]
                            )
                            data.resource_limits.upper_bounds[r][t] = _parse_float(
                                parts[4]
                            )
                elif key == "RESOURCE_CONSTRAINT_COEFFICIENTS":
                    while True:
                        coef_line = f.readline()
                        if not coef_line:
                            break
                        coef_line = coef_line.strip()
                        if not coef_line or coef_line.startswith("%"):
                            continue
                        if coef_line.upper() == "EOF":
                            break
                        parts = coef_line.split()
                        if len(parts) < 4:
                            break
                        block_id, dest_id = int(parts[0]), int(parts[1])
                        resource_id, coef_value = int(parts[2]), _parse_float(parts[3])
                        if block_id not in data.resource_coefficients:
                            data.resource_coefficients[block_id] = {}
                        if dest_id not in data.resource_coefficients[block_id]:
                            data.resource_coefficients[block_id][dest_id] = {}
                        data.resource_coefficients[block_id][dest_id][
                            resource_id
                        ] = coef_value

    # Handle precedences flexibly
    if prec_source is not None:
        if isinstance(prec_source, Precedences):
            data.precedences = prec_source
        else:
            data.precedences = parse_precedences(prec_source)

    return data


UPIT_MODEL = """
# Ultimate Pit Limit Problem (UPIT)
# Maximum-weight closure problem - polynomial time solvable

param NBLOCKS;
set BLOCKS = 0..NBLOCKS-1;
param PROFIT {BLOCKS};
set PREC {BLOCKS};  # Predecessors for each block

var x {BLOCKS} binary;

maximize TotalProfit:
    sum {b in BLOCKS} PROFIT[b] * x[b];

subject to Precedence {b in BLOCKS, i in PREC[b]}:
    x[b] <= x[i];
"""

CPIT_MODEL = """
# Constrained Pit Limit Problem (CPIT)
# NP-hard - maximizes NPV with resource constraints

param NBLOCKS;
param NPERIODS;
param NRESOURCES;
param DISCOUNT_RATE;

set BLOCKS = 0..NBLOCKS-1;
set PERIODS = 0..NPERIODS-1;
set RESOURCES = 0..NRESOURCES-1;

param PROFIT {BLOCKS};
param RES_LB {RESOURCES, PERIODS};
param RES_UB {RESOURCES, PERIODS};
param RES_COEF {BLOCKS, RESOURCES} default 0;
set PREC {BLOCKS};

var x {BLOCKS, PERIODS} binary;

maximize NPV:
    sum {t in PERIODS, b in BLOCKS}
        (1/(1 + DISCOUNT_RATE))^t * PROFIT[b] * x[b,t];

# Strong formulation for precedence
subject to Precedence {b in BLOCKS, i in PREC[b], t in PERIODS}:
    sum {s in 0..t} x[b,s] <= sum {s in 0..t} x[i,s];

# Each block extracted at most once
subject to ExtractOnce {b in BLOCKS}:
    sum {t in PERIODS} x[b,t] <= 1;

# Resource capacity constraints
subject to ResourceBounds {r in RESOURCES, t in PERIODS}:
    RES_LB[r,t] <= sum {b in BLOCKS} RES_COEF[b,r] * x[b,t] <= RES_UB[r,t];
"""

PCPSP_MODEL = """
# Precedence Constrained Production Scheduling Problem (PCPSP)
# Variable cutoff grade with multiple destinations

param NBLOCKS;
param NPERIODS;
param NDESTINATIONS;
param NRESOURCES;
param DISCOUNT_RATE;

set BLOCKS = 0..NBLOCKS-1;
set PERIODS = 0..NPERIODS-1;
set DESTINATIONS = 0..NDESTINATIONS-1;
set RESOURCES = 0..NRESOURCES-1;

param PROFIT {BLOCKS, DESTINATIONS};
param RES_LB {RESOURCES, PERIODS};
param RES_UB {RESOURCES, PERIODS};
param RES_COEF {BLOCKS, RESOURCES, DESTINATIONS} default 0;
set PREC {BLOCKS};

var x {BLOCKS, PERIODS} binary;
var y {BLOCKS, DESTINATIONS, PERIODS} >= 0, <= 1;

maximize NPV:
    sum {t in PERIODS, d in DESTINATIONS, b in BLOCKS}
        (1/(1 + DISCOUNT_RATE))^t * PROFIT[b,d] * y[b,d,t];

subject to Precedence {b in BLOCKS, i in PREC[b], t in PERIODS}:
    sum {s in 0..t} x[b,s] <= sum {s in 0..t} x[i,s];

subject to ExtractOnce {b in BLOCKS}:
    sum {t in PERIODS} x[b,t] <= 1;

# Link x and y: extraction equals sum of destinations
subject to DestinationLink {b in BLOCKS, t in PERIODS}:
    x[b,t] = sum {d in DESTINATIONS} y[b,d,t];

subject to ResourceBounds {r in RESOURCES, t in PERIODS}:
    RES_LB[r,t] <= sum {b in BLOCKS, d in DESTINATIONS} RES_COEF[b,r,d] * y[b,d,t] <= RES_UB[r,t];
"""

# @title


@dataclass
class SolutionUPIT:
    """Solution for UPIT problem."""

    objective: float
    solve_time: float
    status: str
    extracted: Dict[int, bool]  # block_id -> in pit?

    def num_extracted(self) -> int:
        return sum(1 for v in self.extracted.values() if v)


@dataclass
class SolutionCPIT:
    """Solution for CPIT problem."""

    objective: float
    solve_time: float
    status: str
    extraction_period: Dict[
        int, Optional[int]
    ]  # block_id -> period (None if not extracted)

    def num_extracted(self) -> int:
        return sum(1 for v in self.extraction_period.values() if v is not None)

    def blocks_in_period(self, t: int) -> List[int]:
        return [b for b, period in self.extraction_period.items() if period == t]


@dataclass
class SolutionPCPSP:
    """Solution for PCPSP problem."""

    objective: float
    solve_time: float
    status: str
    extraction_period: Dict[int, Optional[int]]
    destination_allocation: Dict[
        int, Dict[int, float]
    ]  # block_id -> {dest_id: fraction}

    def num_extracted(self) -> int:
        return sum(1 for v in self.extraction_period.values() if v is not None)

# @title


def _format_bound(value: float) -> float:
    """Convert infinity to large number for AMPL."""
    if value == math.inf:
        return 1e20
    elif value == -math.inf:
        return -1e20
    return value


def solve_upit(
    data: UPITData,
    solver: str = "gurobi",
    solver_options: str = "",
    verbose: bool = True,
) -> SolutionUPIT:
    """Solve a UPIT problem."""
    if data.precedences is None:
        raise ValueError("Precedences must be loaded")

    ampl = AMPL()
    ampl.eval(UPIT_MODEL)

    # Load parameters
    ampl.param["NBLOCKS"] = data.num_blocks
    for b in range(data.num_blocks):
        ampl.param["PROFIT"][b] = data.objective.get(b, 0.0)

    # Load precedences as indexed sets
    for b in range(data.num_blocks):
        ampl.set["PREC"][b] = data.precedences.get_predecessors(b)

    # Configure and solve
    ampl.option["solver"] = solver
    if solver_options:
        ampl.option[f"{solver}_options"] = solver_options

    if verbose:
        print(f"Solving UPIT: {data.name}")
        print(
            f"  Blocks: {data.num_blocks}, Arcs: {data.precedences.num_precedence_arcs()}"
        )

    ampl.solve()

    # Extract solution
    status = ampl.get_value("solve_result")
    solve_time = ampl.get_value("_solve_time")
    objective = ampl.obj["TotalProfit"].value()

    extracted = {}
    X = ampl.var["x"]
    for b in range(data.num_blocks):
        extracted[b] = X[b].value() > 0.5

    if verbose:
        print(
            f"  Status: {status}, Objective: ${objective:,.2f}, Time: {solve_time:.2f}s"
        )

    return SolutionUPIT(
        objective=objective, solve_time=solve_time, status=status, extracted=extracted
    )


def solve_cpit(
    data: CPITData,
    solver: str = "gurobi",
    solver_options: str = "",
    verbose: bool = True,
) -> SolutionCPIT:
    """Solve a CPIT problem."""
    if data.precedences is None:
        raise ValueError("Precedences must be loaded")

    ampl = AMPL()
    ampl.eval(CPIT_MODEL)

    # Load scalar parameters
    ampl.param["NBLOCKS"] = data.num_blocks
    ampl.param["NPERIODS"] = data.num_periods
    ampl.param["NRESOURCES"] = data.num_resources
    ampl.param["DISCOUNT_RATE"] = data.discount_rate

    # Load profits
    for b in range(data.num_blocks):
        ampl.param["PROFIT"][b] = data.objective.get(b, 0.0)

    # Load resource limits
    for r in range(data.num_resources):
        for t in range(data.num_periods):
            lb = data.resource_limits.lower_bounds.get(r, {}).get(t, -math.inf)
            ub = data.resource_limits.upper_bounds.get(r, {}).get(t, math.inf)
            ampl.param["RES_LB"][r, t] = _format_bound(lb)
            ampl.param["RES_UB"][r, t] = _format_bound(ub)

    # Load resource coefficients
    for b, res_dict in data.resource_coefficients.items():
        for r, coef in res_dict.items():
            ampl.param["RES_COEF"][b, r] = coef

    # Load precedences
    for b in range(data.num_blocks):
        ampl.set["PREC"][b] = data.precedences.get_predecessors(b)

    # Configure and solve
    ampl.option["solver"] = solver
    if solver_options:
        ampl.option[f"{solver}_options"] = solver_options

    if verbose:
        print(f"Solving CPIT: {data.name}")
        print(
            f"  Blocks: {data.num_blocks}, Periods: {data.num_periods}, Resources: {data.num_resources}"
        )

    ampl.solve()

    # Extract solution
    status = ampl.get_value("solve_result")
    solve_time = ampl.get_value("_solve_time")
    objective = ampl.obj["NPV"].value()

    extraction_period = {}
    X = ampl.var["x"]
    for b in range(data.num_blocks):
        extraction_period[b] = None
        for t in range(data.num_periods):
            if X[b, t].value() > 0.5:
                extraction_period[b] = t
                break

    if verbose:
        print(f"  Status: {status}, NPV: ${objective:,.2f}, Time: {solve_time:.2f}s")

    return SolutionCPIT(
        objective=objective,
        solve_time=solve_time,
        status=status,
        extraction_period=extraction_period,
    )


def solve_pcpsp(
    data: PCPSPData,
    solver: str = "gurobi",
    solver_options: str = "",
    verbose: bool = True,
) -> SolutionPCPSP:
    """Solve a PCPSP problem."""
    if data.precedences is None:
        raise ValueError("Precedences must be loaded")

    ampl = AMPL()
    ampl.eval(PCPSP_MODEL)

    # Load scalar parameters
    ampl.param["NBLOCKS"] = data.num_blocks
    ampl.param["NPERIODS"] = data.num_periods
    ampl.param["NDESTINATIONS"] = data.num_destinations
    ampl.param["NRESOURCES"] = data.num_resources
    ampl.param["DISCOUNT_RATE"] = data.discount_rate

    # Load profits (per block and destination)
    for b in range(data.num_blocks):
        for d in range(data.num_destinations):
            val = data.objective.get(b, {}).get(d, 0.0)
            ampl.param["PROFIT"][b, d] = val

    # Load resource limits
    for r in range(data.num_resources):
        for t in range(data.num_periods):
            lb = data.resource_limits.lower_bounds.get(r, {}).get(t, -math.inf)
            ub = data.resource_limits.upper_bounds.get(r, {}).get(t, math.inf)
            ampl.param["RES_LB"][r, t] = _format_bound(lb)
            ampl.param["RES_UB"][r, t] = _format_bound(ub)

    # Load resource coefficients (per block, resource, destination)
    for b, dest_dict in data.resource_coefficients.items():
        for d, res_dict in dest_dict.items():
            for r, coef in res_dict.items():
                ampl.param["RES_COEF"][b, r, d] = coef

    # Load precedences
    for b in range(data.num_blocks):
        ampl.set["PREC"][b] = data.precedences.get_predecessors(b)

    # Configure and solve
    ampl.option["solver"] = solver
    if solver_options:
        ampl.option[f"{solver}_options"] = solver_options

    if verbose:
        print(f"Solving PCPSP: {data.name}")
        print(
            f"  Blocks: {data.num_blocks}, Periods: {data.num_periods}, Destinations: {data.num_destinations}"
        )

    ampl.solve()

    # Extract solution
    status = ampl.get_value("solve_result")
    solve_time = ampl.get_value("_solve_time")
    objective = ampl.obj["NPV"].value()

    extraction_period = {}
    destination_allocation = {}
    X = ampl.var["x"]
    Y = ampl.var["y"]

    for b in range(data.num_blocks):
        extraction_period[b] = None
        for t in range(data.num_periods):
            if X[b, t].value() > 0.5:
                extraction_period[b] = t
                destination_allocation[b] = {}
                for d in range(data.num_destinations):
                    y_val = Y[b, d, t].value()
                    if y_val > 1e-6:
                        destination_allocation[b][d] = y_val
                break

    if verbose:
        print(f"  Status: {status}, NPV: ${objective:,.2f}, Time: {solve_time:.2f}s")

    return SolutionPCPSP(
        objective=objective,
        solve_time=solve_time,
        status=status,
        extraction_period=extraction_period,
        destination_allocation=destination_allocation,
    )

# Base URL for raw GitHub content
BASE_URL = "https://raw.githubusercontent.com/ampl/colab.ampl.com/master/authors/eduardosalaz/minelib/data"


def fetch_minelib_file(instance: str, extension: str) -> StringIO:
    """
    Fetch a MineLib file from GitHub and return as file-like object.

    Args:
        instance: Instance name (e.g., 'newman1', 'zuck_small')
        extension: File extension (e.g., 'prec', 'upit', 'cpit', 'pcpsp', 'blocks')

    Returns:
        StringIO object that can be used like a file
    """
    url = f"{BASE_URL}/{instance}/{instance}.{extension}"
    response = requests.get(url)
    response.raise_for_status()  # Raise exception if download fails
    return StringIO(response.text)


# Example usage
INSTANCE = "newman1"

# Fetch files as file-like objects
prec_file = fetch_minelib_file(INSTANCE, "prec")
upit_file = fetch_minelib_file(INSTANCE, "upit")
cpit_file = fetch_minelib_file(INSTANCE, "cpit")
pcpsp_file = fetch_minelib_file(INSTANCE, "pcpsp")

precedences = parse_precedences(prec_file)
# Parse each problem type and assign precedences
upit_data = parse_upit(upit_file)
upit_data.precedences = precedences

cpit_data = parse_cpit(cpit_file)
cpit_data.precedences = precedences

pcpsp_data = parse_pcpsp(pcpsp_file)
pcpsp_data.precedences = precedences

print(f"Instance: {upit_data.name}")
print(f"Blocks: {upit_data.num_blocks}")
print(f"Precedence arcs: {upit_data.precedences.num_precedence_arcs()}")

# Count ore vs waste blocks
positive = sum(1 for v in upit_data.objective.values() if v > 0)
negative = sum(1 for v in upit_data.objective.values() if v < 0)
print(f"Ore blocks (positive value): {positive}")
print(f"Waste blocks (negative value): {negative}")

# %%
# Solve UPIT
upit_solution = solve_upit(upit_data, solver="gurobi", verbose=True)

print(
    f"\nBlocks in optimal pit: {upit_solution.num_extracted()} / {upit_data.num_blocks}"
)

Instance: Newman1
Blocks: 1060
Precedence arcs: 3922
Ore blocks (positive value): 545
Waste blocks (negative value): 515
Solving UPIT: Newman1
  Blocks: 1060, Arcs: 3922
Gurobi 13.0.0:Gurobi 13.0.0: optimal solution; objective 26086899.03
0 simplex iterations
  Status: solved, Objective: $26,086,899.03, Time: 0.05s

Blocks in optimal pit: 1059 / 1060

print(f"Instance: {cpit_data.name}")
print(f"Blocks: {cpit_data.num_blocks}, Periods: {cpit_data.num_periods}")
print(f"Resources: {cpit_data.num_resources}, Discount rate: {cpit_data.discount_rate}")

# %%
# Solve CPIT (may take longer due to time-indexed variables)
cpit_solution = solve_cpit(
    cpit_data, solver="gurobi", solver_options="timelimit=300 outlev=1", verbose=True
)

print(f"\nBlocks extracted: {cpit_solution.num_extracted()} / {cpit_data.num_blocks}")

# %%
# Analyze extraction schedule
print("\nExtraction Schedule:")
print(f"{'Period':<10} {'Blocks':<15} {'Cumulative':<15}")
print("-" * 40)

cumulative = 0
for t in range(cpit_data.num_periods):
    count = len(cpit_solution.blocks_in_period(t))
    cumulative += count
    if count > 0:
        print(f"{t:<10} {count:<15} {cumulative:<15}")

Instance: Newman1
Blocks: 1060, Periods: 6
Resources: 2, Discount rate: 0.08
Solving CPIT: Newman1
  Blocks: 1060, Periods: 6, Resources: 2
Gurobi 13.0.0:   lim:time = 300
Set parameter LogToConsole to value 1
  tech:outlev = 1

AMPL MP initial flat model has 6360 variables (0 integer, 6360 binary);
Objectives: 1 linear; 
Constraints:  24604 linear;

AMPL MP final model has 6360 variables (0 integer, 6360 binary);
Objectives: 1 linear; 
Constraints:  24604 linear;


Set parameter InfUnbdInfo to value 1
Gurobi Optimizer version 13.0.0 build v13.0.0rc1 (linux64 - "Ubuntu 22.04.5 LTS")

CPU model: Intel(R) Xeon(R) CPU @ 2.20GHz, instruction set [SSE2|AVX|AVX2]
Thread count: 1 physical cores, 2 logical processors, using up to 2 threads

Non-default parameters:
TimeLimit  300
InfUnbdInfo  1

Optimize a model with 24604 rows, 6360 columns and 180876 nonzeros (Max)
Model fingerprint: 0xe72d45d4
Model has 6360 linear objective coefficients
Variable types: 0 continuous, 6360 integer (0 binary)
Coefficient statistics:
  Matrix range     [1e+00, 6e+03]
  Objective range  [5e+00, 2e+05]
  Bounds range     [1e+00, 1e+00]
  RHS range        [1e+00, 2e+06]
Found heuristic solution: objective 989307.55494
Found heuristic solution: objective 1822786.3759
Presolve removed 605 rows and 147 columns
Presolve time: 0.92s
Presolved: 23999 rows, 6213 columns, 173942 nonzeros
Variable types: 0 continuous, 6213 integer (6213 binary)
Found heuristic solution: objective 2014319.6047

Root relaxation: objective 2.441235e+07, 7977 iterations, 3.20 seconds (2.53 work units)

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0 2.4412e+07    0 2270 2014319.60 2.4412e+07  1112%     -    5s
     0     0 2.4401e+07    0 2252 2014319.60 2.4401e+07  1111%     -    5s
     0     0 2.4396e+07    0 2199 2014319.60 2.4396e+07  1111%     -    6s
H    0     0                    5009329.3557 2.4396e+07   387%     -    6s
H    0     0                    5083504.7013 2.4396e+07   380%     -    6s
H    0     0                    7971178.5077 2.4396e+07   206%     -    6s
H    0     0                    1.108309e+07 2.4396e+07   120%     -    6s
H    0     0                    1.175160e+07 2.4396e+07   108%     -    6s
H    0     0                    1.175792e+07 2.4396e+07   107%     -    6s
H    0     0                    1.175991e+07 2.4396e+07   107%     -    6s
H    0     0                    1.176253e+07 2.4396e+07   107%     -    6s
H    0     0                    1.176367e+07 2.4396e+07   107%     -    6s
H    0     0                    1.176437e+07 2.4396e+07   107%     -    6s
H    0     0                    1.774323e+07 2.4396e+07  37.5%     -    6s
H    0     0                    1.994616e+07 2.4396e+07  22.3%     -    6s
H    0     0                    2.006290e+07 2.4396e+07  21.6%     -    6s
H    0     0                    2.012337e+07 2.4396e+07  21.2%     -    6s
H    0     0                    2.015188e+07 2.4396e+07  21.1%     -    6s
H    0     0                    2.152247e+07 2.4396e+07  13.4%     -    6s
H    0     0                    2.178737e+07 2.4396e+07  12.0%     -    6s
     0     0 2.4374e+07    0 2189 2.1787e+07 2.4374e+07  11.9%     -    9s
H    0     0                    2.187666e+07 2.4374e+07  11.4%     -    9s
H    0     0                    2.229493e+07 2.4374e+07  9.33%     -    9s
H    0     0                    2.320495e+07 2.4374e+07  5.04%     -    9s
H    0     0                    2.321282e+07 2.4374e+07  5.00%     -    9s
     0     0 2.4364e+07    0 2153 2.3213e+07 2.4364e+07  4.96%     -   11s
     0     0 2.4364e+07    0 2153 2.3213e+07 2.4364e+07  4.96%     -   12s
H    0     0                    2.397581e+07 2.4364e+07  1.62%     -   13s
H    0     0                    2.400214e+07 2.4364e+07  1.51%     -   13s
H    0     0                    2.400453e+07 2.4364e+07  1.50%     -   13s
H    0     0                    2.400466e+07 2.4364e+07  1.50%     -   13s
H    0     0                    2.402174e+07 2.4364e+07  1.42%     -   14s
H    0     0                    2.402250e+07 2.4364e+07  1.42%     -   14s
H    0     0                    2.402271e+07 2.4364e+07  1.42%     -   14s
H    0     0                    2.402577e+07 2.4364e+07  1.41%     -   16s
H    0     0                    2.402630e+07 2.4364e+07  1.40%     -   16s
H    0     0                    2.402762e+07 2.4364e+07  1.40%     -   16s
H    0     0                    2.402834e+07 2.4364e+07  1.40%     -   16s
H    0     0                    2.403287e+07 2.4364e+07  1.38%     -   16s
H    0     0                    2.405813e+07 2.4364e+07  1.27%     -   16s
H    0     0                    2.406321e+07 2.4364e+07  1.25%     -   16s
H    0     0                    2.406449e+07 2.4364e+07  1.24%     -   16s
H    0     0                    2.406449e+07 2.4364e+07  1.24%     -   17s
H    0     0                    2.406510e+07 2.4364e+07  1.24%     -   17s
H    0     0                    2.406510e+07 2.4364e+07  1.24%     -   17s
H    0     0                    2.406710e+07 2.4364e+07  1.23%     -   17s
     0     2 2.4364e+07    0 2153 2.4067e+07 2.4364e+07  1.23%     -   18s
H    2     2                    2.406911e+07 2.4364e+07  1.22%   302   18s
H    2     2                    2.408161e+07 2.4364e+07  1.17%   302   18s
H    6     6                    2.408365e+07 2.4364e+07  1.16%   215   19s
H    6     6                    2.408514e+07 2.4364e+07  1.16%   215   19s
H    6     6                    2.408516e+07 2.4364e+07  1.16%   215   19s
H    6     6                    2.408593e+07 2.4364e+07  1.15%   215   19s
H    6     6                    2.408689e+07 2.4364e+07  1.15%   215   19s
H    8     8                    2.410380e+07 2.4364e+07  1.08%   247   20s
H    8     8                    2.416307e+07 2.4364e+07  0.83%   247   20s
H    8     8                    2.416659e+07 2.4364e+07  0.82%   247   20s
H    8     8                    2.417281e+07 2.4364e+07  0.79%   247   20s
H    8     8                    2.417423e+07 2.4364e+07  0.78%   247   20s
H    8     8                    2.417553e+07 2.4364e+07  0.78%   247   20s
H    8     8                    2.417562e+07 2.4364e+07  0.78%   247   20s
H   29    29                    2.417570e+07 2.4364e+07  0.78%  85.1   20s
H   36    31                    2.417573e+07 2.4364e+07  0.78%  74.8   20s
H   36    30                    2.417686e+07 2.4364e+07  0.77%  74.8   20s
    65    35     cutoff    9      2.4177e+07 2.4314e+07  0.57%   118   25s
   116    35 2.4264e+07    5 1420 2.4177e+07 2.4283e+07  0.44%   134   30s

Cutting planes:
  Gomory: 2
  Cover: 2

Explored 172 nodes (33288 simplex iterations) in 34.94 seconds (25.38 work units)
Thread count was 2 (of 2 available processors)

Solution count 10: 2.41769e+07 2.41757e+07 2.41757e+07 ... 2.41038e+07

Optimal solution found (tolerance 1.00e-04)
Best objective 2.417686090020e+07, best bound 2.417867161503e+07, gap 0.0075%
Gurobi 13.0.0: optimal solution; objective 24176860.9
33288 simplex iterations
172 branching nodes
absmipgap=1810.71, relmipgap=7.48945e-05
  Status: solved, NPV: $24,176,860.90, Time: 34.74s

Blocks extracted: 1059 / 1060

Extraction Schedule:
Period     Blocks          Cumulative     
----------------------------------------
0          406             406            
1          366             772            
2          287             1059           

print(f"Instance: {pcpsp_data.name}")
print(f"Blocks: {pcpsp_data.num_blocks}, Destinations: {pcpsp_data.num_destinations}")

# %%
# Solve PCPSP
pcpsp_solution = solve_pcpsp(
    pcpsp_data, solver="gurobi", solver_options="timelimit=300 outlev=1", verbose=True
)

print(f"\nBlocks extracted: {pcpsp_solution.num_extracted()}")

# Destination analysis
dest_names = {0: "Processing", 1: "Waste"}
dest_counts = {d: 0 for d in range(pcpsp_data.num_destinations)}
for b, alloc in pcpsp_solution.destination_allocation.items():
    for d, frac in alloc.items():
        if frac > 0.5:
            dest_counts[d] += 1

print("\nDestination breakdown:")
for d, count in dest_counts.items():
    name = dest_names.get(d, f"Dest {d}")
    print(f"  {name}: {count} blocks")

Instance: Newman1
Blocks: 1060, Destinations: 2
Solving PCPSP: Newman1
  Blocks: 1060, Periods: 6, Destinations: 2
Gurobi 13.0.0:   lim:time = 300
Set parameter LogToConsole to value 1
  tech:outlev = 1

AMPL MP initial flat model has 19080 variables (0 integer, 6360 binary);
Objectives: 1 linear; 
Constraints:  30964 linear;

AMPL MP final model has 19080 variables (0 integer, 6360 binary);
Objectives: 1 linear; 
Constraints:  30964 linear;


Set parameter InfUnbdInfo to value 1
Gurobi Optimizer version 13.0.0 build v13.0.0rc1 (linux64 - "Ubuntu 22.04.5 LTS")

CPU model: Intel(R) Xeon(R) CPU @ 2.20GHz, instruction set [SSE2|AVX|AVX2]
Thread count: 1 physical cores, 2 logical processors, using up to 2 threads

Non-default parameters:
TimeLimit  300
InfUnbdInfo  1

Optimize a model with 30964 rows, 19080 columns and 209244 nonzeros (Max)
Model fingerprint: 0x4b2e1105
Model has 12720 linear objective coefficients
Variable types: 12720 continuous, 6360 integer (0 binary)
Coefficient statistics:
  Matrix range     [1e+00, 6e+03]
  Objective range  [5e+00, 6e+06]
  Bounds range     [1e+00, 1e+00]
  RHS range        [1e+00, 2e+06]
Found heuristic solution: objective -0.0000000
Presolve removed 3653 rows and 6249 columns
Presolve time: 1.36s
Presolved: 27311 rows, 12831 columns, 144650 nonzeros
Variable types: 6618 continuous, 6213 integer (6213 binary)

Root relaxation: objective 2.441860e+07, 12388 iterations, 2.96 seconds (1.87 work units)

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0 2.4419e+07    0 1910   -0.00000 2.4419e+07      -     -    5s
H    0     0                    2.290322e+07 2.4419e+07  6.62%     -    6s
H    0     0                    2.291816e+07 2.4419e+07  6.55%     -    6s
H    0     0                    2.306184e+07 2.4419e+07  5.88%     -    6s
H    0     0                    2.306357e+07 2.4419e+07  5.88%     -    6s
H    0     0                    2.309448e+07 2.4419e+07  5.73%     -    6s
H    0     0                    2.309491e+07 2.4419e+07  5.73%     -    6s
H    0     0                    2.309866e+07 2.4419e+07  5.71%     -    6s
     0     0 2.4401e+07    0 1953 2.3099e+07 2.4401e+07  5.64%     -    8s
     0     0 2.4401e+07    0 1953 2.3099e+07 2.4401e+07  5.64%     -   11s
     0     0 2.4401e+07    0 1953 2.3099e+07 2.4401e+07  5.64%     -   11s
H    0     0                    2.404638e+07 2.4401e+07  1.47%     -   12s
H    0     0                    2.405396e+07 2.4401e+07  1.44%     -   12s
H    0     0                    2.405402e+07 2.4401e+07  1.44%     -   12s
H    0     0                    2.405558e+07 2.4401e+07  1.43%     -   12s
H    0     0                    2.405570e+07 2.4401e+07  1.43%     -   12s
H    0     0                    2.405846e+07 2.4401e+07  1.42%     -   12s
     0     2 2.4401e+07    0 1953 2.4058e+07 2.4401e+07  1.42%     -   13s
     2     2 2.4397e+07    1 1893 2.4058e+07 2.4397e+07  1.41%   558   15s
H    7     5                    2.405850e+07 2.4330e+07  1.13%   379   18s
H    7     5                    2.408752e+07 2.4330e+07  1.01%   379   18s
H    7     5                    2.408975e+07 2.4330e+07  1.00%   379   18s
H    7     5                    2.409005e+07 2.4330e+07  1.00%   379   18s
H    9     7                    2.409712e+07 2.4321e+07  0.93%   404   19s
H    9     7                    2.409746e+07 2.4321e+07  0.93%   404   19s
H    9     7                    2.410061e+07 2.4321e+07  0.91%   404   19s
H   14    12                    2.410101e+07 2.4321e+07  0.91%   332   20s
H   14    12                    2.411557e+07 2.4321e+07  0.85%   332   20s
H   14    12                    2.411642e+07 2.4321e+07  0.85%   332   20s
H   14    12                    2.411839e+07 2.4321e+07  0.84%   332   20s
H   14    12                    2.411882e+07 2.4321e+07  0.84%   332   20s
H   14    12                    2.411952e+07 2.4321e+07  0.84%   332   20s
H   14    12                    2.412359e+07 2.4321e+07  0.82%   332   20s
H   18    16                    2.412369e+07 2.4321e+07  0.82%   320   21s
H   18    16                    2.412392e+07 2.4321e+07  0.82%   320   21s
H   18    16                    2.412431e+07 2.4321e+07  0.82%   320   21s
H   18    16                    2.412453e+07 2.4321e+07  0.81%   320   21s
*  144   122              99    2.414971e+07 2.4321e+07  0.71%  74.5   24s
   147   121 2.4249e+07    3 1220 2.4150e+07 2.4321e+07  0.71%  76.6   25s
H  149    47                    2.416327e+07 2.4321e+07  0.65%  79.2   25s
H  149    18                    2.417519e+07 2.4321e+07  0.60%  79.2   26s
   175    11     cutoff    9      2.4175e+07 2.4264e+07  0.37%   101   30s
H  190     4                    2.417577e+07 2.4196e+07  0.08%   104   32s

Cutting planes:
  Gomory: 1

Explored 235 nodes (34098 simplex iterations) in 33.10 seconds (22.52 work units)
Thread count was 2 (of 2 available processors)

Solution count 10: 2.41758e+07 2.41752e+07 2.41633e+07 ... 2.41195e+07

Optimal solution found (tolerance 1.00e-04)
Best objective 2.417576894988e+07, best bound 2.417813514902e+07, gap 0.0098%
Gurobi 13.0.0: optimal solution; objective 24175768.95
34098 simplex iterations
235 branching nodes
absmipgap=2366.2, relmipgap=9.78748e-05
  Status: solved, NPV: $24,175,768.95, Time: 32.88s

Blocks extracted: 1059

Destination breakdown:
  Processing: 487 blocks
  Waste: 572 blocks

print("=" * 65)
print("Comparison of Formulations")
print("=" * 65)
print(f"{'Model':<10} {'Objective ($)':<20} {'Blocks':<12} {'Time (s)':<10}")
print("-" * 65)
print(
    f"{'UPIT':<10} {upit_solution.objective:>18,.0f} {upit_solution.num_extracted():>10} {upit_solution.solve_time:>9.2f}"
)
print(
    f"{'CPIT':<10} {cpit_solution.objective:>18,.0f} {cpit_solution.num_extracted():>10} {cpit_solution.solve_time:>9.2f}"
)
print(
    f"{'PCPSP':<10} {pcpsp_solution.objective:>18,.0f} {pcpsp_solution.num_extracted():>10} {pcpsp_solution.solve_time:>9.2f}"
)
print("-" * 65)

# NPV reduction due to discounting and constraints
reduction = (
    (upit_solution.objective - cpit_solution.objective) / upit_solution.objective * 100
)
print(f"\nNPV reduction (UPIT -> CPIT): {reduction:.1f}%")

=================================================================
Comparison of Formulations
=================================================================
Model      Objective ($)        Blocks       Time (s)  
-----------------------------------------------------------------
UPIT               26,086,899       1059      0.05
CPIT               24,176,861       1059     34.74
PCPSP              24,175,769       1059     32.88
-----------------------------------------------------------------

NPV reduction (UPIT -> CPIT): 7.3%

7. Conclusions#

This notebook demonstrated solving MineLib open pit mining problems with AMPL:

UPIT: Polynomial-time pit boundary optimization (network flow)
CPIT: NP-hard scheduling with time periods and resource constraints
PCPSP: Variable cutoff grades with destination optimization

Key observations:

Strong precedence formulation improves LP bounds
Discounting and constraints reduce achievable NPV vs. UPIT
Variable cutoff (PCPSP) can recover some value vs. fixed cutoff (CPIT)

Symbol	Description
\(t \in \mathcal{T}\)	Set of time periods \(t\) in the horizon
\(b \in \mathcal{B}\)	Set of blocks \(b\)
\(b' \in \mathcal{B}_b\)	Set of blocks \(b'\) that are predecessor blocks for block \(b\)
\(r \in \mathcal{R}\)	Set of operational resource types \(r\)
\(d \in \mathcal{D}\)	Set of destinations \(d\)

Symbol	Description
\(p_b\) (\(\hat{p}_{bt}\), \(\check{p}_{bd}\), \(\tilde{p}_{bdt}\))	Profit obtained from extracting (and processing) block \(b\) (at time period \(t\) and/or sending it to destination \(d\)) ($)
\(\alpha\)	Discount rate used in computing the objective function (profit) coefficients
\(q_{br}\) (\(\hat{q}_{brd}\))	Amount of operational resource \(r\) used to extract and, if applicable, process block \(b\) (when sent to destination \(d\)) (tons)
\(\underline{R}_{rt}\)	Minimum availability of operational resource \(r\) in time period \(t\) (tons)
\(\overline{R}_{rt}\)	Maximum availability of operational resource \(r\) in time period \(t\) (tons)
\(A\)	Arbitrary constraint coefficients on general side constraints
\(\underline{a}, \overline{a}\)	Arbitrary lower and upper bounds, respectively, on general side constraints

Symbol	Description
\(\hat{x}_b\)	1 if block \(b\) is in the final pit design, 0 otherwise
\(x_{bt}\)	1 if we extract block \(b\) in time period \(t\), 0 otherwise
\(y_{bdt}\)	The amount of block \(b\) sent to destination \(d\) in time period \(t\) (%)

MineLib in AMPL and amplpy#

1. Introduction to Open Pit Mining #

Key Concepts#

Why Optimization?#

2. The MineLib Library#

Available Instances#

Problem Types#

3. Mathematical Formulations#

3.1 Notation#

3.2 The Ultimate Pit Problem (UPIT)#

3.3 The Constrained Pit Limit Problem (CPIT)#

3.4 The Precedence Constrained Production Scheduling Problem (PCPSP)#

4. MineLib File Format#

4.1 Block-Model Descriptor File (.blocks)#

4.2 Block-Precedence Descriptor File (.prec)#

4.3 Optimization-Model Descriptor File (.upit, .cpit, .pcpsp)#

7. Conclusions#

4.1 Block-Model Descriptor File (`.blocks`)#

4.2 Block-Precedence Descriptor File (`.prec`)#

4.3 Optimization-Model Descriptor File (`.upit`, `.cpit`, `.pcpsp`)#