Production Forecaster

Production Forecaster

Forecast oil and gas production using industry-standard decline curve analysis methods. Supports Arps decline models, type curve generation, and EUR (Estimated Ultimate Recovery) calculations for reserve estimation.

When to Use

• Forecasting future production for oil and gas wells
• Estimating EUR (Estimated Ultimate Recovery)
• Generating type curves for field development planning
• Fitting decline parameters to historical production data
• Comparing well performance across different fields
• Supporting reserve booking and economic evaluation
• Analyzing production trends and anomalies

Core Pattern

Historical Production → Decline Curve Fit → Parameter Estimation → Forecast → EUR

Decline Curve Models

Arps Decline Equations

Exponential Decline (b = 0):

q(t) = qi * exp(-Di * t)

Hyperbolic Decline (0 < b < 1):

q(t) = qi / (1 + b * Di * t)^(1/b)

Harmonic Decline (b = 1):

q(t) = qi / (1 + Di * t)

Where:

• q(t) = Production rate at time t
• qi = Initial production rate
• Di = Initial decline rate
• b = Decline exponent (b-factor)
• t = Time

Modified Hyperbolic Decline

For long-term forecasts, use terminal decline rate to prevent over-optimistic recoveries:

q(t) = qi / (1 + b * Di * t)^(1/b)     when D(t) > D_min
q(t) = q_switch * exp(-D_min * t')     when D(t) <= D_min

Where D_min is the terminal decline rate (typically 5-8%/year) and t' is time from switch point.

Duong Decline Model (Unconventional)

For unconventional wells with transient linear flow:

q(t) = q1 * t^(-m) * exp(-a/(1-m) * (t^(1-m) - 1))

Where:

• q1 = Rate at t=1 day
• a = Intercept (controls EUR)
• m = Slope (controls decline steepness)

Stretched Exponential Decline

Alternative for unconventional with varying flow regimes:

q(t) = qi * exp(-(t/τ)^n)

Where:

• τ = Characteristic time constant
• n = Stretch parameter (0 < n < 1)

Type Curve Analysis Framework

Normalization Methods

Type Curve Bins

Group wells for representative type curves:

Type Curve Confidence Intervals

Implementation

Data Models

from dataclasses import dataclass, field
from datetime import date
from typing import Optional, List, Dict, Tuple
from enum import Enum
import numpy as np
import pandas as pd

class DeclineType(Enum):
    """Decline curve types."""
    EXPONENTIAL = "exponential"       # b = 0
    HYPERBOLIC = "hyperbolic"         # 0 < b < 1
    HARMONIC = "harmonic"             # b = 1
    MODIFIED_HYPERBOLIC = "modified"  # Hyperbolic with terminal decline
    DUONG = "duong"                   # Unconventional transient
    STRETCHED_EXP = "stretched"       # Stretched exponential

class NormalizationMethod(Enum):
    """Type curve normalization methods."""
    PEAK = "peak"
    FIRST_MONTH = "first_month"
    IP_30 = "30_day_ip"
    IP_90 = "90_day_ip"
    MOVING_AVG = "moving_average"
    CUMULATIVE = "cumulative_at_time"

class TypeCurveBinType(Enum):
    """Type curve binning categories."""
    FORMATION = "formation"
    COMPLETION = "completion"
    WATER_DEPTH = "water_depth"
    LATERAL_LENGTH = "lateral_length"
    PROPPANT = "proppant"
    VINTAGE = "vintage"

class ProductionPhase(Enum):
    """Well production phases."""
    BUILDUP = "buildup"
    PLATEAU = "plateau"
    DECLINE = "decline"
    ABANDONMENT = "abandonment"

@dataclass
class DeclineParameters:
    """Decline curve parameters."""
    qi: float  # Initial rate (bbls/day or mcf/day)
    di: float  # Initial decline rate (fraction/year)
    b: float   # Decline exponent
    decline_type: DeclineType = DeclineType.HYPERBOLIC

    # Optional constraints
    min_rate: float = 0.0  # Economic limit
    max_time: float = 50.0  # Years
    d_min: Optional[float] = None  # Terminal decline rate

    # Duong model parameters (unconventional)
    a: Optional[float] = None  # Duong intercept
    m: Optional[float] = None  # Duong slope

    # Stretched exponential parameters
    tau: Optional[float] = None  # Time constant
    n: Optional[float] = None    # Stretch parameter

    def __post_init__(self):
        """Validate parameters and set decline type."""
        if self.a is not None and self.m is not None:
            self.decline_type = DeclineType.DUONG
        elif self.tau is not None and self.n is not None:
            self.decline_type = DeclineType.STRETCHED_EXP
        elif self.d_min is not None and self.b > 0:
            self.decline_type = DeclineType.MODIFIED_HYPERBOLIC
        elif self.b == 0:
            self.decline_type = DeclineType.EXPONENTIAL
        elif self.b == 1:
            self.decline_type = DeclineType.HARMONIC
        else:
            self.decline_type = DeclineType.HYPERBOLIC

@dataclass
class TypeCurveBin:
    """Type curve bin definition for well grouping."""
    bin_type: TypeCurveBinType
    bin_name: str
    criteria: Dict[str, Any] = field(default_factory=dict)
    wells: List[str] = field(default_factory=list)

    def matches(self, well_metadata: Dict[str, Any]) -> bool:
        """Check if well matches bin criteria."""
        for key, value in self.criteria.items():
            if key not in well_metadata:
                return False
            well_value = well_metadata[key]
            if isinstance(value, tuple):
                # Range check (min, max)
                if not (value[0] <= well_value <= value[1]):
                    return False
            elif isinstance(value, list):
                # In list check
                if well_value not in value:
                    return False
            else:
                # Exact match
                if well_value != value:
                    return False
        return True

@dataclass
class TypeCurveResult:
    """Type curve analysis results."""
    months: List[int]
    p10_rates: List[float]
    p50_rates: List[float]
    p90_rates: List[float]
    mean_rates: List[float]
    well_counts: List[int]

    # Fitted parameters for each percentile
    p10_params: Optional[DeclineParameters] = None
    p50_params: Optional[DeclineParameters] = None
    p90_params: Optional[DeclineParameters] = None

    # Validation metrics
    r_squared: Optional[float] = None
    rmse: Optional[float] = None
    mae: Optional[float] = None

    # EUR estimates
    eur_p10: Optional[float] = None
    eur_p50: Optional[float] = None
    eur_p90: Optional[float] = None

    def to_dataframe(self) -> pd.DataFrame:
        """Convert to DataFrame."""
        return pd.DataFrame({
            'month': self.months,
            'P10': self.p10_rates,
            'P50': self.p50_rates,
            'P90': self.p90_rates,
            'mean': self.mean_rates,
            'well_count': self.well_counts
        })

@dataclass
class ProductionRecord:
    """Single production record."""
    date: date
    oil_rate: float = 0.0      # bbls/day
    gas_rate: float = 0.0      # mcf/day
    water_rate: float = 0.0    # bbls/day
    days_on: int = 30          # Days producing

    @property
    def oil_volume(self) -> float:
        """Monthly oil volume in barrels."""
        return self.oil_rate * self.days_on

    @property
    def gas_volume(self) -> float:
        """Monthly gas volume in mcf."""
        return self.gas_rate * self.days_on

    @property
    def gor(self) -> float:
        """Gas-oil ratio (mcf/bbl)."""
        if self.oil_rate > 0:
            return self.gas_rate / self.oil_rate
        return 0.0

@dataclass
class ForecastResult:
    """Production forecast results."""
    dates: List[date]
    oil_rates: List[float]
    gas_rates: List[float]
    cumulative_oil: List[float]
    cumulative_gas: List[float]

    # Decline parameters used
    parameters: DeclineParameters

    # Key metrics
    eur_oil: float = 0.0       # Estimated Ultimate Recovery (bbls)
    eur_gas: float = 0.0       # EUR gas (mcf)
    remaining_oil: float = 0.0  # Remaining reserves
    remaining_gas: float = 0.0

    def to_dataframe(self) -> pd.DataFrame:
        """Convert to DataFrame."""
        return pd.DataFrame({
            'date': self.dates,
            'oil_rate': self.oil_rates,
            'gas_rate': self.gas_rates,
            'cumulative_oil': self.cumulative_oil,
            'cumulative_gas': self.cumulative_gas
        })

Decline Curve Analyzer

import numpy as np
from scipy.optimize import curve_fit, minimize
from typing import Tuple, Optional
import pandas as pd

class DeclineCurveAnalyzer:
    """
    Analyze and fit decline curves to production data.
    """

    def __init__(self, production_data: pd.DataFrame):
        """
        Initialize with production data.

        Args:
            production_data: DataFrame with 'date' and 'rate' columns
        """
        self.data = production_data.copy()
        self._prepare_data()

    def _prepare_data(self):
        """Prepare data for analysis."""
        self.data['date'] = pd.to_datetime(self.data['date'])
        self.data = self.data.sort_values('date')

        # Calculate time from first production
        first_date = self.data['date'].min()
        self.data['time_years'] = (
            (self.data['date'] - first_date).dt.days / 365.25
        )

    @staticmethod
    def exponential_decline(t: np.ndarray, qi: float, di: float) -> np.ndarray:
        """Exponential decline model."""
        return qi * np.exp(-di * t)

    @staticmethod
    def hyperbolic_decline(t: np.ndarray, qi: float, di: float,
                          b: float) -> np.ndarray:
        """Hyperbolic decline model."""
        return qi / np.power(1 + b * di * t, 1/b)

    @staticmethod
    def harmonic_decline(t: np.ndarray, qi: float, di: float) -> np.ndarray:
        """Harmonic decline model."""
        return qi / (1 + di * t)

    def fit_exponential(self) -> DeclineParameters:
        """Fit exponential decline curve."""
        t = self.data['time_years'].values
        q = self.data['rate'].values

        # Initial guesses
        qi_guess = q[0]
        di_guess = 0.3

        try:
            popt, _ = curve_fit(
                self.exponential_decline,
                t, q,
                p0=[qi_guess, di_guess],
                bounds=([0, 0], [qi_guess * 2, 2.0]),
                maxfev=5000
            )
            return DeclineParameters(qi=popt[0], di=popt[1], b=0.0)
        except:
            return DeclineParameters(qi=qi_guess, di=di_guess, b=0.0)

    def fit_hyperbolic(self, b_range: Tuple[float, float] = (0.1, 0.9)
                      ) -> DeclineParameters:
        """Fit hyperbolic decline curve."""
        t = self.data['time_years'].values
        q = self.data['rate'].values

        # Initial guesses
        qi_guess = q[0]
        di_guess = 0.3
        b_guess = 0.5

        try:
            popt, _ = curve_fit(
                self.hyperbolic_decline,
                t, q,
                p0=[qi_guess, di_guess, b_guess],
                bounds=(
                    [0, 0, b_range[0]],
                    [qi_guess * 2, 2.0, b_range[1]]
                ),
                maxfev=5000
            )
            return DeclineParameters(qi=popt[0], di=popt[1], b=popt[2])
        except:
            return DeclineParameters(qi=qi_guess, di=di_guess, b=b_guess)

    @staticmethod
    def duong_decline(t: np.ndarray, q1: float, a: float,
                     m: float) -> np.ndarray:
        """Duong decline model for unconventional wells."""
        # Avoid t=0 issues
        t = np.maximum(t, 1e-6)
        return q1 * np.power(t, -m) * np.exp(-a / (1 - m) * (np.power(t, 1 - m) - 1))

    @staticmethod
    def stretched_exponential(t: np.ndarray, qi: float, tau: float,
                             n: float) -> np.ndarray:
        """Stretched exponential decline model."""
        return qi * np.exp(-np.power(t / tau, n))

    @staticmethod
    def modified_hyperbolic(t: np.ndarray, qi: float, di: float,
                           b: float, d_min: float) -> np.ndarray:
        """Modified hyperbolic with terminal decline rate."""
        result = np.zeros_like(t)
        for i, ti in enumerate(t):
            # Current decline rate
            d_t = di / (1 + b * di * ti)
            if d_t > d_min:
                # Still in hyperbolic phase
                result[i] = qi / np.power(1 + b * di * ti, 1/b)
            else:
                # Switch to exponential at d_min
                t_switch = (di / d_min - 1) / (b * di)
                q_switch = qi / np.power(1 + b * di * t_switch, 1/b)
                result[i] = q_switch * np.exp(-d_min * (ti - t_switch))
        return result

    def fit_duong(self, t_range: Tuple[float, float] = (0.5, 3.0)
                 ) -> DeclineParameters:
        """Fit Duong decline model (unconventional)."""
        t = self.data['time_years'].values
        q = self.data['rate'].values

        # Convert to days for Duong
        t_days = t * 365.25 + 1  # Add 1 to avoid t=0

        # Initial guesses
        q1_guess = q[0]
        a_guess = 1.0
        m_guess = 1.2

        try:
            popt, _ = curve_fit(
                self.duong_decline,
                t_days, q,
                p0=[q1_guess, a_guess, m_guess],
                bounds=(
                    [0, 0.1, 0.5],
                    [q1_guess * 3, 5.0, 2.0]
                ),
                maxfev=10000
            )
            return DeclineParameters(
                qi=popt[0], di=0, b=0,
                a=popt[1], m=popt[2],
                decline_type=DeclineType.DUONG
            )
        except:
            return DeclineParameters(
                qi=q1_guess, di=0, b=0,
                a=a_guess, m=m_guess,
                decline_type=DeclineType.DUONG
            )

    def fit_modified_hyperbolic(self, d_min: float = 0.06
                               ) -> DeclineParameters:
        """Fit modified hyperbolic with terminal decline."""
        t = self.data['time_years'].values
        q = self.data['rate'].values

        # First fit standard hyperbolic
        hyp_params = self.fit_hyperbolic()

        # Use those as starting point
        try:
            def mod_hyp_fixed_dmin(t, qi, di, b):
                return self.modified_hyperbolic(t, qi, di, b, d_min)

            popt, _ = curve_fit(
                mod_hyp_fixed_dmin,
                t, q,
                p0=[hyp_params.qi, hyp_params.di, hyp_params.b],
                bounds=(
                    [0, 0, 0.1],
                    [hyp_params.qi * 2, 2.0, 1.5]
                ),
                maxfev=5000
            )
            return DeclineParameters(
                qi=popt[0], di=popt[1], b=popt[2],
                d_min=d_min,
                decline_type=DeclineType.MODIFIED_HYPERBOLIC
            )
        except:
            return DeclineParameters(
                qi=hyp_params.qi, di=hyp_params.di, b=hyp_params.b,
                d_min=d_min,
                decline_type=DeclineType.MODIFIED_HYPERBOLIC
            )

    def fit_stretched_exponential(self) -> DeclineParameters:
        """Fit stretched exponential decline model."""
        t = self.data['time_years'].values
        q = self.data['rate'].values

        # Initial guesses
        qi_guess = q[0]
        tau_guess = 1.0
        n_guess = 0.7

        try:
            popt, _ = curve_fit(
                self.stretched_exponential,
                t, q,
                p0=[qi_guess, tau_guess, n_guess],
                bounds=(
                    [0, 0.1, 0.1],
                    [qi_guess * 2, 20.0, 1.0]
                ),
                maxfev=5000
            )
            return DeclineParameters(
                qi=popt[0], di=0, b=0,
                tau=popt[1], n=popt[2],
                decline_type=DeclineType.STRETCHED_EXP
            )
        except:
            return DeclineParameters(
                qi=qi_guess, di=0, b=0,
                tau=tau_guess, n=n_guess,
                decline_type=DeclineType.STRETCHED_EXP
            )

    def fit_best_model(self, include_unconventional: bool = False
                      ) -> Tuple[DeclineParameters, str]:
        """
        Fit all models and return best fit.

        Args:
            include_unconventional: Include Duong and stretched exp models

        Returns:
            Tuple of (best parameters, model name)
        """
        models = {
            'exponential': self.fit_exponential(),
            'hyperbolic': self.fit_hyperbolic(),
        }

        if include_unconventional:
            models['duong'] = self.fit_duong()
            models['stretched'] = self.fit_stretched_exponential()

        t = self.data['time_years'].values
        q = self.data['rate'].values

        best_model = None
        best_rmse = float('inf')
        best_params = None

        for name, params in models.items():
            if name == 'exponential':
                predicted = self.exponential_decline(t, params.qi, params.di)
            elif name == 'duong':
                predicted = self.duong_decline(t * 365.25 + 1, params.qi, params.a, params.m)
            elif name == 'stretched':
                predicted = self.stretched_exponential(t, params.qi, params.tau, params.n)
            else:
                predicted = self.hyperbolic_decline(
                    t, params.qi, params.di, params.b
                )

            rmse = np.sqrt(np.mean((q - predicted) ** 2))

            if rmse < best_rmse:
                best_rmse = rmse
                best_model = name
                best_params = params

        return best_params, best_model

    def calculate_validation_metrics(self, params: DeclineParameters
                                    ) -> Dict[str, float]:
        """
        Calculate validation metrics for fitted parameters.

        Returns:
            Dict with R², RMSE, MAE, MAPE
        """
        t = self.data['time_years'].values
        q = self.data['rate'].values

        # Calculate predicted values
        if params.decline_type == DeclineType.EXPONENTIAL:
            predicted = self.exponential_decline(t, params.qi, params.di)
        elif params.decline_type == DeclineType.DUONG:
            predicted = self.duong_decline(t * 365.25 + 1, params.qi, params.a, params.m)
        elif params.decline_type == DeclineType.STRETCHED_EXP:
            predicted = self.stretched_exponential(t, params.qi, params.tau, params.n)
        elif params.decline_type == DeclineType.MODIFIED_HYPERBOLIC:
            predicted = self.modified_hyperbolic(t, params.qi, params.di, params.b, params.d_min)
        else:
            predicted = self.hyperbolic_decline(t, params.qi, params.di, params.b)

        # R-squared
        ss_res = np.sum((q - predicted) ** 2)
        ss_tot = np.sum((q - np.mean(q)) ** 2)
        r_squared = 1 - (ss_res / ss_tot) if ss_tot > 0 else 0

        # RMSE
        rmse = np.sqrt(np.mean((q - predicted) ** 2))

        # MAE
        mae = np.mean(np.abs(q - predicted))

        # MAPE (Mean Absolute Percentage Error)
        mask = q > 0
        mape = np.mean(np.abs((q[mask] - predicted[mask]) / q[mask])) * 100

        return {
            'r_squared': r_squared,
            'rmse': rmse,
            'mae': mae,
            'mape': mape
        }

    def calculate_eur(self, params: DeclineParameters,
                     economic_limit: float = 10.0,
                     max_years: float = 50.0) -> float:
        """
        Calculate Estimated Ultimate Recovery.

        Args:
            params: Decline parameters
            economic_limit: Minimum economic rate
            max_years: Maximum forecast period
        """
        if params.b == 0:
            # Exponential: EUR with economic limit and time cap
            # If already below economic limit, no recoverable reserves
            if params.qi <= economic_limit:
                return 0.0
            # Time to reach economic limit: t = ln(qi/q_limit) / di
            t_econ = np.log(params.qi / economic_limit) / params.di
            t_max = min(max_years, max(0.0, t_econ))
            # Cumulative: Np = (qi/di) * (1 - exp(-di*t))
            return params.qi * 365.25 / params.di * (1 - np.exp(-params.di * t_max))
        elif params.b == 1:
            # Harmonic: EUR is infinite, use time limit
            # Cumulative = qi/di * ln(1 + di*t)
            t_max = min(max_years, (params.qi / economic_limit - 1) / params.di)
            return params.qi * 365.25 / params.di * np.log(1 + params.di * t_max)
        else:
            # Hyperbolic
            # Cumulative = qi^b / (di * (1-b)) * (qi^(1-b) - q_limit^(1-b))
            q_limit = max(economic_limit, 0.1)
            factor = params.qi ** params.b / (params.di * (1 - params.b))
            cum = factor * (params.qi ** (1 - params.b) - q_limit ** (1 - params.b))
            return cum * 365.25

Production Forecaster

from datetime import date, timedelta
from typing import List, Optional
import numpy as np
import pandas as pd

class ProductionForecaster:
    """
    Generate production forecasts using decline curve analysis.
    """

    def __init__(self, parameters: DeclineParameters,
                 start_date: date = None,
                 cumulative_to_date: float = 0.0):
        """
        Initialize forecaster.

        Args:
            parameters: Decline curve parameters
            start_date: Forecast start date
            cumulative_to_date: Production already recovered
        """
        self.params = parameters
        self.start_date = start_date or date.today()
        self.cum_to_date = cumulative_to_date

    def _calculate_rate(self, t: float) -> float:
        """Calculate rate at time t (years)."""
        if self.params.b == 0:
            return self.params.qi * np.exp(-self.params.di * t)
        elif self.params.b == 1:
            return self.params.qi / (1 + self.params.di * t)
        else:
            return self.params.qi / np.power(
                1 + self.params.b * self.params.di * t,
                1 / self.params.b
            )

    def _calculate_cumulative(self, t: float) -> float:
        """Calculate cumulative production at time t (years)."""
        if self.params.b == 0:
            # Exponential
            return self.params.qi / self.params.di * (
                1 - np.exp(-self.params.di * t)
            ) * 365.25
        elif self.params.b == 1:
            # Harmonic
            return self.params.qi / self.params.di * np.log(
                1 + self.params.di * t
            ) * 365.25
        else:
            # Hyperbolic
            q_t = self._calculate_rate(t)
            factor = self.params.qi ** self.params.b / (
                self.params.di * (1 - self.params.b)
            )
            return factor * (
                self.params.qi ** (1 - self.params.b) -
                q_t ** (1 - self.params.b)
            ) * 365.25

    def forecast(self, years: float = 30.0,
                interval_months: int = 1,
                economic_limit: float = 10.0,
                gor: float = 1.0) -> ForecastResult:
        """
        Generate production forecast.

        Args:
            years: Forecast period in years
            interval_months: Output interval in months
            economic_limit: Minimum economic rate (bbls/day or mcf/day)
            gor: Gas-oil ratio for gas calculation

        Returns:
            ForecastResult with forecast data
        """
        dates = []
        oil_rates = []
        gas_rates = []
        cumulative_oil = []
        cumulative_gas = []

        current_date = self.start_date
        total_months = int(years * 12)

        for month in range(0, total_months, interval_months):
            t = month / 12.0  # Time in years

            rate = self._calculate_rate(t)

            # Check economic limit
            if rate < economic_limit:
                break

            cum = self._calculate_cumulative(t) + self.cum_to_date

            dates.append(current_date)
            oil_rates.append(rate)
            gas_rates.append(rate * gor)
            cumulative_oil.append(cum)
            cumulative_gas.append(cum * gor)

            # Advance date by interval_months
            new_month = current_date.month + interval_months
            new_year = current_date.year + (new_month - 1) // 12
            new_month = ((new_month - 1) % 12) + 1
            current_date = date(new_year, new_month, 1)

        # Calculate EUR
        analyzer = DeclineCurveAnalyzer.__new__(DeclineCurveAnalyzer)
        eur_oil = analyzer.calculate_eur(
            self.params, economic_limit, years
        ) + self.cum_to_date

        return ForecastResult(
            dates=dates,
            oil_rates=oil_rates,
            gas_rates=gas_rates,
            cumulative_oil=cumulative_oil,
            cumulative_gas=cumulative_gas,
            parameters=self.params,
            eur_oil=eur_oil,
            eur_gas=eur_oil * gor,
            remaining_oil=eur_oil - self.cum_to_date,
            remaining_gas=(eur_oil - self.cum_to_date) * gor
        )

    def to_dataframe(self, forecast: ForecastResult) -> pd.DataFrame:
        """Convert forecast to DataFrame."""
        return forecast.to_dataframe()

Type Curve Generator

from typing import List, Dict, Optional, Tuple, Callable
import numpy as np
import pandas as pd
from scipy import stats
from dataclasses import dataclass

class TypeCurveGenerator:
    """
    Generate type curves from multiple well production histories.
    Supports multiple normalization methods, binning, and probabilistic analysis.
    """

    # Standard water depth bins for GOM
    WATER_DEPTH_BINS = {
        'shallow': (0, 1000),
        'deep': (1000, 5000),
        'ultra_deep': (5000, 15000)
    }

    def __init__(self, wells: List[pd.DataFrame],
                 metadata: Optional[List[Dict]] = None):
        """
        Initialize with list of well production DataFrames.

        Args:
            wells: List of DataFrames with 'date' and 'rate' columns
            metadata: Optional list of metadata dicts per well
        """
        self.wells = wells
        self.metadata = metadata or [{} for _ in wells]
        self.normalized_wells = []
        self.normalization_factors = []
        self.bins = {}

    def normalize_wells(self, method: NormalizationMethod = NormalizationMethod.PEAK,
                       time_reference: int = 30) -> List[pd.DataFrame]:
        """
        Normalize wells for type curve generation.

        Args:
            method: Normalization method to use
            time_reference: Days for IP calculation or cumulative reference
        """
        normalized = []
        self.normalization_factors = []

        for well in self.wells:
            df = well.copy()
            df['date'] = pd.to_datetime(df['date'])
            df = df.sort_values('date')

            # Calculate days and months from start
            first_date = df['date'].min()
            df['days'] = (df['date'] - first_date).dt.days
            df['month'] = (df['days'] / 30.44).astype(int)

            # Calculate normalization factor based on method
            if method == NormalizationMethod.PEAK:
                norm_factor = df['rate'].max()
            elif method == NormalizationMethod.FIRST_MONTH:
                norm_factor = df[df['month'] == 0]['rate'].mean()
            elif method == NormalizationMethod.IP_30:
                mask = df['days'] <= 30
                norm_factor = df[mask]['rate'].mean() if mask.sum() > 0 else df['rate'].iloc[0]
            elif method == NormalizationMethod.IP_90:
                mask = df['days'] <= 90
                norm_factor = df[mask]['rate'].mean() if mask.sum() > 0 else df['rate'].iloc[0]
            elif method == NormalizationMethod.MOVING_AVG:
                df['rate_smooth'] = df['rate'].rolling(window=30, min_periods=1).mean()
                norm_factor = df['rate_smooth'].max()
            elif method == NormalizationMethod.CUMULATIVE:
                # Normalize so all wells have same cumulative at reference time
                ref_cum = df[df['days'] <= time_reference]['rate'].sum()
                norm_factor = ref_cum / time_reference if ref_cum > 0 else df['rate'].iloc[0]
            else:
                norm_factor = df['rate'].max()

            # Avoid division by zero
            norm_factor = max(norm_factor, 0.1)

            df['normalized_rate'] = df['rate'] / norm_factor
            self.normalization_factors.append(norm_factor)

            normalized.append(df)

        self.normalized_wells = normalized
        return normalized

    def create_bins(self, bin_type: TypeCurveBinType,
                   custom_bins: Optional[Dict[str, Any]] = None) -> Dict[str, TypeCurveBin]:
        """
        Create type curve bins based on well metadata.

        Args:
            bin_type: Type of binning to apply
            custom_bins: Optional custom bin definitions
        """
        bins = {}

        if bin_type == TypeCurveBinType.WATER_DEPTH:
            bin_defs = custom_bins or self.WATER_DEPTH_BINS
            for name, (min_wd, max_wd) in bin_defs.items():
                bins[name] = TypeCurveBin(
                    bin_type=bin_type,
                    bin_name=name,
                    criteria={'water_depth_ft': (min_wd, max_wd)}
                )

        elif bin_type == TypeCurveBinType.FORMATION:
            # Extract unique formations from metadata
            formations = set()
            for meta in self.metadata:
                if 'formation' in meta:
                    formations.add(meta['formation'])
            for formation in formations:
                bins[formation] = TypeCurveBin(
                    bin_type=bin_type,
                    bin_name=formation,
                    criteria={'formation': formation}
                )

        elif bin_type == TypeCurveBinType.VINTAGE:
            # Bin by year of first production
            for i, well in enumerate(self.wells):
                df = well.copy()
                df['date'] = pd.to_datetime(df['date'])
                year = df['date'].min().year
                year_str = str(year)
                if year_str not in bins:
                    bins[year_str] = TypeCurveBin(
                        bin_type=bin_type,
                        bin_name=year_str,
                        criteria={'vintage_year': year}
                    )

        # Assign wells to bins
        for i, meta in enumerate(self.metadata):
            for name, bin_def in bins.items():
                if bin_def.matches(meta):
                    bin_def.wells.append(i)

        self.bins[bin_type] = bins
        return bins

    def generate_type_curve(self,
                           percentiles: List[float] = [10, 50, 90],
                           min_wells: int = 3,
                           well_indices: Optional[List[int]] = None
                           ) -> TypeCurveResult:
        """
        Generate type curve with percentile ranges.

        Args:
            percentiles: Percentiles to calculate (P10, P50, P90)
            min_wells: Minimum wells required per month
            well_indices: Optional subset of well indices to use

        Returns:
            TypeCurveResult with percentile curves
        """
        if not self.normalized_wells:
            self.normalize_wells()

        # Select wells
        if well_indices:
            wells_to_use = [self.normalized_wells[i] for i in well_indices]
        else:
            wells_to_use = self.normalized_wells

        # Find maximum months
        max_months = max(df['month'].max() for df in wells_to_use)

        # Collect rates by month
        monthly_rates = {m: [] for m in range(int(max_months) + 1)}

        for df in wells_to_use:
            for _, row in df.iterrows():
                month = int(row['month'])
                if month in monthly_rates:
                    monthly_rates[month].append(row['normalized_rate'])

        # Calculate percentiles and statistics
        months = []
        p10_rates, p50_rates, p90_rates = [], [], []
        mean_rates, well_counts = [], []

        for month in sorted(monthly_rates.keys()):
            rates = monthly_rates[month]
            if len(rates) >= min_wells:
                months.append(month)
                # Note: P10 = high side, P90 = low side (exceedance)
                p10_rates.append(np.percentile(rates, 90))
                p50_rates.append(np.percentile(rates, 50))
                p90_rates.append(np.percentile(rates, 10))
                mean_rates.append(np.mean(rates))
                well_counts.append(len(rates))

        return TypeCurveResult(
            months=months,
            p10_rates=p10_rates,
            p50_rates=p50_rates,
            p90_rates=p90_rates,
            mean_rates=mean_rates,
            well_counts=well_counts
        )

    def fit_type_curve(self, percentile: float = 50,
                      model: DeclineType = DeclineType.HYPERBOLIC,
                      well_indices: Optional[List[int]] = None
                      ) -> Tuple[DeclineParameters, Dict[str, float]]:
        """
        Fit decline curve to type curve percentile.

        Args:
            percentile: Percentile to fit (10, 50, or 90)
            model: Decline model type to fit
            well_indices: Optional subset of wells

        Returns:
            Tuple of (fitted parameters, validation metrics)
        """
        result = self.generate_type_curve([percentile], well_indices=well_indices)

        # Get rate column
        if percentile == 10:
            rates = result.p10_rates
        elif percentile == 90:
            rates = result.p90_rates
        else:
            rates = result.p50_rates

        # Create analyzer with type curve data
        tc_df = pd.DataFrame({
            'date': pd.date_range(start='2020-01-01',
                                  periods=len(result.months), freq='MS'),
            'rate': rates
        })

        analyzer = DeclineCurveAnalyzer(tc_df)

        if model == DeclineType.EXPONENTIAL:
            params = analyzer.fit_exponential()
        elif model == DeclineType.HYPERBOLIC:
            params = analyzer.fit_hyperbolic()
        elif model == DeclineType.DUONG:
            params = analyzer.fit_duong()
        elif model == DeclineType.MODIFIED_HYPERBOLIC:
            params = analyzer.fit_modified_hyperbolic()
        else:
            params, _ = analyzer.fit_best_model()

        # Calculate validation metrics
        t = np.array(result.months) / 12.0
        if params.decline_type == DeclineType.EXPONENTIAL:
            predicted = analyzer.exponential_decline(t, params.qi, params.di)
        elif params.decline_type == DeclineType.DUONG:
            predicted = analyzer.duong_decline(t * 365.25, params.qi, params.a, params.m)
        else:
            predicted = analyzer.hyperbolic_decline(t, params.qi, params.di, params.b)

        actual = np.array(rates)
        ss_res = np.sum((actual - predicted) ** 2)
        ss_tot = np.sum((actual - np.mean(actual)) ** 2)
        r_squared = 1 - (ss_res / ss_tot) if ss_tot > 0 else 0
        rmse = np.sqrt(np.mean((actual - predicted) ** 2))
        mae = np.mean(np.abs(actual - predicted))

        metrics = {
            'r_squared': r_squared,
            'rmse': rmse,
            'mae': mae
        }

        return params, metrics

    def generate_bin_type_curves(self, bin_type: TypeCurveBinType
                                ) -> Dict[str, TypeCurveResult]:
        """
        Generate type curves for each bin.

        Args:
            bin_type: Type of binning to use

        Returns:
            Dict mapping bin name to TypeCurveResult
        """
        if bin_type not in self.bins:
            self.create_bins(bin_type)

        results = {}
        for name, bin_def in self.bins[bin_type].items():
            if len(bin_def.wells) >= 3:
                results[name] = self.generate_type_curve(
                    well_indices=bin_def.wells
                )

        return results

    def scale_type_curve(self, type_curve: TypeCurveResult,
                        target_ip: float,
                        percentile: str = 'P50') -> pd.DataFrame:
        """
        Scale a normalized type curve to a target initial production rate.

        Args:
            type_curve: TypeCurveResult to scale
            target_ip: Target initial production rate
            percentile: Which percentile to scale (P10, P50, P90)

        Returns:
            DataFrame with scaled production rates
        """
        tc_df = type_curve.to_dataframe()

        # Get the rate column for the percentile
        rate_col = percentile

        # Scale to target IP
        first_rate = tc_df[rate_col].iloc[0]
        scale_factor = target_ip / first_rate if first_rate > 0 else target_ip

        scaled_df = pd.DataFrame({
            'month': tc_df['month'],
            f'{rate_col}_scaled': tc_df[rate_col] * scale_factor,
            'well_count': tc_df['well_count']
        })

        return scaled_df

    def calculate_eur_distribution(self, economic_limit: float = 10.0,
                                   max_years: float = 30.0
                                   ) -> Dict[str, float]:
        """
        Calculate EUR distribution (P10/P50/P90) from type curves.

        Args:
            economic_limit: Minimum economic rate
            max_years: Maximum forecast period

        Returns:
            Dict with P10, P50, P90 EUR estimates
        """
        eur_values = []

        for i, well in enumerate(self.normalized_wells):
            # De-normalize to get actual rates
            norm_factor = self.normalization_factors[i]
            actual_rates = well['normalized_rate'] * norm_factor

            # Fit decline and calculate EUR
            well_df = pd.DataFrame({
                'date': well['date'],
                'rate': actual_rates
            })
            analyzer = DeclineCurveAnalyzer(well_df)
            params, _ = analyzer.fit_best_model()
            eur = analyzer.calculate_eur(params, economic_limit, max_years)
            eur_values.append(eur)

        return {
            'P10': np.percentile(eur_values, 90),  # High side
            'P50': np.percentile(eur_values, 50),
            'P90': np.percentile(eur_values, 10),  # Low side
            'mean': np.mean(eur_values),
            'std': np.std(eur_values)
        }

Report Generator

import plotly.graph_objects as go
from plotly.subplots import make_subplots
from pathlib import Path
impor

…(truncated)