Mathematical Modeling for Software Developers

import numpy as np
from dataclasses import dataclass
from typing import List, Tuple

@dataclass
class LinearRegressionModel:
    coefficients: np.ndarray
    intercept: float
    
    def predict(self, X: np.ndarray) -> np.ndarray:
        return np.dot(X, self.coefficients) + self.intercept
    
    def r_squared(self, X: np.ndarray, y: np.ndarray) -> float:
        y_pred = self.predict(X)
        ss_res = np.sum((y - y_pred) ** 2)
        ss_tot = np.sum((y - np.mean(y)) ** 2)
        return 1 - (ss_res / ss_tot)


def gradient_descent_linear_regression(
    X: np.ndarray, 
    y: np.ndarray,
    learning_rate: float = 0.01,
    iterations: int = 1000
) -> LinearRegressionModel:
    """Train linear regression using gradient descent."""
    m, n = X.shape
    
    # Initialize parameters
    coefficients = np.zeros(n)
    intercept = 0.0
    
    for _ in range(iterations):
        predictions = np.dot(X, coefficients) + intercept
        errors = predictions - y
        
        # Update
        coefficients -= learning_rate * (1/m) * np.dot(X.T, errors)
        intercept -= learning_rate * (1/m) * np.sum(errors)
    
    return LinearRegressionModel(coefficients, intercept)


def normal_equation_linear_regression(
    X: np.ndarray, 
    y: np.ndarray
) -> LinearRegressionModel:
    """Closed-form solution using normal equation."""
    X_with_bias = np.column_stack([np.ones(X.shape[0]), X])
    theta = np.linalg.lstsq(X_with_bias, y, rcond=None)[0]
    
    return LinearRegressionModel(theta[1:], theta[0])


def polynomial_regression(
    X: np.ndarray, 
    y: np.ndarray,
    degree: int = 2
) -> LinearRegressionModel:
    """Polynomial regression for non-linear relationships."""
    # Create polynomial features
    X_poly = np.column_stack([X ** d for d in range(degree + 1)])
    
    return normal_equation_linear_regression(X_poly, y)

def sigmoid(z: np.ndarray) -> np.ndarray:
    """Sigmoid activation function."""
    return 1 / (1 + np.exp(-np.clip(z, -500, 500)))


def logistic_regression(
    X: np.ndarray,
    y: np.ndarray,
    learning_rate: float = 0.1,
    iterations: int = 1000
) -> Tuple[np.ndarray, float]:
    """Binary classification using logistic regression."""
    m, n = X.shape
    weights = np.zeros(n)
    bias = 0.0
    
    for _ in range(iterations):
        linear_output = np.dot(X, weights) + bias
        predictions = sigmoid(linear_output)
        
        # Gradient descent
        dw = (1/m) * np.dot(X.T, (predictions - y))
        db = (1/m) * np.sum(predictions - y)
        
        weights -= learning_rate * dw
        bias -= learning_rate * db
    
    return weights, bias


def predict_logistic(
    X: np.ndarray,
    weights: np.ndarray,
    bias: float,
    threshold: float = 0.5
) -> np.ndarray:
    """Make predictions with logistic regression."""
    probabilities = sigmoid(np.dot(X, weights) + bias)
    return (probabilities >= threshold).astype(int)


class LogisticRegression:
    def __init__(self, learning_rate: float = 0.1, iterations: int = 1000):
        self.learning_rate = learning_rate
        self.iterations = iterations
        self.weights = None
        self.bias = None
        
    def fit(self, X: np.ndarray, y: np.ndarray):
        self.weights, self.bias = logistic_regression(
            X, y, self.learning_rate, self.iterations
        )
        
    def predict_proba(self, X: np.ndarray) -> np.ndarray:
        return sigmoid(np.dot(X, self.weights) + self.bias)
    
    def predict(self, X: np.ndarray, threshold: float = 0.5) -> np.ndarray:
        return predict_logistic(X, self.weights, self.bias, threshold)

class ExponentialSmoothing:
    """Exponential smoothing for time series forecasting."""
    
    def __init__(self, alpha: float = 0.3):
        self.alpha = alpha
        self.level = None
        self.smoothing_level = None
        
    def fit(self, data: np.ndarray):
        self.level = data[0]
        for value in data[1:]:
            self.level = self.alpha * value + (1 - self.alpha) * self.level
        return self
    
    def forecast(self, steps: int) -> np.ndarray:
        return np.full(steps, self.level)


class HoltWinters:
    """Holt-Winters triple exponential smoothing."""
    
    def __init__(self, alpha: float = 0.3, beta: float = 0.1, gamma: float = 0.1, 
                 seasonal_period: int = 12):
        self.alpha = alpha
        self.beta = beta
        self.gamma = gamma
        self.seasonal_period = seasonal_period
        
    def fit(self, data: np.ndarray):
        n = len(data)
        self.level = np.mean(data[:self.seasonal_period])
        self.trend = (np.mean(data[self.seasonal_period:2*self.seasonal_period]) - 
                      np.mean(data[:self.seasonal_period])) / self.seasonal_period
        
        # Initialize seasonal components
        self.seasonal = np.zeros(self.seasonal_period)
        for i in range(self.seasonal_period):
            self.seasonal[i] = data[i] - self.level
            
        # Fit remaining data
        for i in range(self.seasonal_period, n):
            last_level = self.level
            self.level = self.alpha * (data[i] - self.seasonal[i % self.seasonal_period]) + \
                        (1 - self.alpha) * (last_level + self.trend)
            self.trend = self.beta * (self.level - last_level) + \
                        (1 - self.beta) * self.trend
            self.seasonal[i % self.seasonal_period] = \
                self.gamma * (data[i] - self.level) + \
                (1 - self.gamma) * self.seasonal[i % self.seasonal_period]
        
        return self
    
    def forecast(self, steps: int) -> np.ndarray:
        predictions = np.zeros(steps)
        for h in range(steps):
            predictions[h] = (self.level + (h + 1) * self.trend + 
                           self.seasonal[(self.seasonal_period - 1 + h) % self.seasonal_period])
        return predictions


class ARIMA:
    """Simplified ARIMA model for demonstration."""
    
    def __init__(self, p: int = 1, d: int = 0, q: int = 1):
        self.p = p  # AR order
        self.d = d  # Differencing
        self.q = int(q)  # MA order
        self.ar_params = None
        self.ma_params = None
        
    def difference(self, data: np.ndarray, d: int) -> np.ndarray:
        """Apply differencing d times."""
        diff = data.copy()
        for _ in range(d):
            diff = np.diff(diff)
        return diff
    
    def fit(self, data: np.ndarray):
        differenced = self.difference(data, self.d)
        
        # Simplified: use last p values as AR coefficients
        if self.p > 0:
            self.ar_params = np.zeros(self.p)
            for i in range(self.p):
                self.ar_params[i] = 0.5 / self.p
                
        # Simplified: use last q values as MA coefficients
        if self.q > 0:
            self.ma_params = np.zeros(self.q)
            for i in range(self.q):
                self.ma_params[i] = 0.3 / self.q
                
        self.last_values = differenced[-max(self.p, self.q):]
        return self
    
    def forecast(self, steps: int) -> np.ndarray:
        predictions = np.zeros(steps)
        
        for h in range(steps):
            pred = 0
            
            # AR component
            if self.p > 0 and len(self.last_values) >= self.p:
                for i in range(self.p):
                    pred += self.ar_params[i] * self.last_values[-(i + 1)]
            
            # MA component
            if self.q > 0 and len(self.last_values) >= self.q:
                for i in range(self.q):
                    pred += self.ma_params[i] * self.last_values[-(i + 1)]
            
            predictions[h] = pred
            
        return predictions

import random

class MonteCarloSimulation:
    """Monte Carlo simulation framework."""
    
    def __init__(self, num_simulations: int = 10000):
        self.num_simulations = num_simulations
        self.results = []
        
    def run_simulation(self, model_fn):
        """Run simulation with given model function."""
        self.results = [model_fn() for _ in range(self.num_simulations)]
        return self.results
    
    def mean(self) -> float:
        return np.mean(self.results)
    
    def std(self) -> float:
        return np.std(self.results)
    
    def confidence_interval(self, confidence: float = 0.95) -> Tuple[float, float]:
        """Calculate confidence interval."""
        alpha = 1 - confidence
        lower = np.percentile(self.results, 100 * alpha / 2)
        upper = np.percentile(self.results, 100 * (1 - alpha / 2))
        return lower, upper
    
    def percentile(self, p: float) -> float:
        return np.percentile(self.results, p)


def simulate_stock_price(S0: float, mu: float, sigma: float, 
                        T: float, steps: int) -> float:
    """Simulate stock price using geometric Brownian motion."""
    dt = T / steps
    price = S0
    
    for _ in range(steps):
        dW = np.random.normal(0, np.sqrt(dt))
        price *= np.exp((mu - 0.5 * sigma**2) * dt + sigma * dW)
    
    return price


def option_pricing_simulation(
    S0: float,      # Current stock price
    K: float,       # Strike price
    T: float,       # Time to maturity
    r: float,       # Risk-free rate
    sigma: float,   # Volatility
    num_simulations: int = 100000,
    option_type: str = 'call'
) -> float:
    """
    Price options using Monte Carlo simulation.
    """
    payoffs = []
    
    for _ in range(num_simulations):
        # Simulate final stock price
        ST = simulate_stock_price(S0, r, sigma, T, 100)
        
        # Calculate payoff
        if option_type == 'call':
            payoff = max(ST - K, 0)
        else:
            payoff = max(K - ST, 0)
            
        payoffs.append(payoff)
    
    # Discount to present value
    option_price = np.exp(-r * T) * np.mean(payoffs)
    
    return option_price


def simulate_inventory_system(
    daily_demand_mean: float,
    daily_demand_std: float,
    ordering_cost: float,
    holding_cost: float,
    stockout_cost: float,
    initial_stock: int,
    num_days: int
) -> Dict[str, float]:
    """
    Simulate inventory management with uncertainty.
    """
    stock = initial_stock
    total_cost = 0
    
    for day in range(num_days):
        demand = max(0, np.random.normal(daily_demand_mean, daily_demand_std))
        demand = int(round(demand))
        
        if stock >= demand:
            stock -= demand
            holding = holding_cost * stock
            total_cost += holding
        else:
            stock = 0
            stockout = stockout_cost * (demand - stock)
            total_cost += stockout
            
        # Order to restore stock
        order_quantity = initial_stock - stock
        if order_quantity > 0:
            total_cost += ordering_cost * order_quantity
            stock = initial_stock
    
    return {
        'total_cost': total_cost,
        'average_daily_cost': total_cost / num_days,
        'final_stock': stock
    }

from dataclasses import dataclass
from queue import PriorityQueue
from typing import Optional

@dataclass
class Event:
    time: float
    event_type: str
    data: dict = None
    
    def __lt__(self, other):
        return self.time < other.time


class DiscreteEventSimulation:
    """Discrete event simulation framework."""
    
    def __init__(self):
        self.events = PriorityQueue()
        self.current_time = 0
        self.results = {}
        
    def schedule(self, time: float, event_type: str, data: dict = None):
        event = Event(self.current_time + time, event_type, data)
        self.events.put(event)
        
    def run(self, end_time: float = float('inf')):
        while not self.events.empty():
            event = self.events.get()
            
            if event.time > end_time:
                break
                
            self.current_time = event.time
            self.process_event(event)
            
    def process_event(self, event: Event):
        """Override in subclass."""
        pass


class QueueSimulation(DiscreteEventSimulation):
    """Simulate a queueing system (M/M/1)."""
    
    def __init__(self, arrival_rate: float, service_rate: float):
        super().__init__()
        self.arrival_rate = arrival_rate
        self.service_rate = service_rate
        self.customers_served = 0
        self.total_wait_time = 0
        self.queue_length = 0
        self.max_queue_length = 0
        
    def process_event(self, event: Event):
        if event.event_type == 'arrival':
            self.handle_arrival(event)
        elif event.event_type == 'departure':
            self.handle_departure(event)
            
    def handle_arrival(self, event: Event):
        if self.queue_length == 0:
            # Start service immediately
            service_time = np.random.exponential(1 / self.service_rate)
            self.schedule(service_time, 'departure')
        else:
            # Add to queue
            self.queue_length += 1
            self.max_queue_length = max(self.max_queue_length, self.queue_length)
            
        # Schedule next arrival
        inter_arrival = np.random.exponential(1 / self.arrival_rate)
        self.schedule(inter_arrival, 'arrival')
        
    def handle_departure(self, event: Event):
        self.customers_served += 1
        self.total_wait_time += self.queue_length * event.data.get('wait_time', 0)
        
        if self.queue_length > 0:
            self.queue_length -= 1
            service_time = np.random.exponential(1 / self.service_rate)
            self.schedule(service_time, 'departure')
    
    def get_statistics(self):
        return {
            'customers_served': self.customers_served,
            'max_queue_length': self.max_queue_length,
            'avg_wait_time': self.total_wait_time / max(1, self.customers_served),
            'utilization': self.customers_served / (self.current_time * self.service_rate)
        }

from scipy.optimize import linprog

def linear_programming_example():
    """
    Maximize: 3x + 2y
    Subject to:
        x + y <= 4
        x <= 2
        y <= 3
        x, y >= 0
    """
    # Coefficients of objective function (minimize)
    c = [-3, -2]  # Negative because linprog minimizes
    
    # Inequality constraints (A_ub @ x <= b_ub)
    A_ub = [
        [1, 1],  # x + y <= 4
        [1, 0],  # x <= 2
        [0, 1],  # y <= 3
    ]
    b_ub = [4, 2, 3]
    
    # Bounds
    bounds = [(0, None), (0, None)]
    
    result = linprog(c, A_ub, b_ub, bounds=bounds)
    
    print(f"Optimal solution: x={result.x[0]:.2f}, y={result.x[1]:.2f}")
    print(f"Maximum value: {-result.fun:.2f}")
    

def integer_programming_example():
    """Knapsack problem using integer programming."""
    from scipy.optimize import milp, LinearConstraint, Bounds
    
    # Values and weights
    values = [60, 100, 120, 90, 80, 50]
    weights = [10, 20, 15, 8, 25, 12]
    capacity = 50
    
    # Maximize value (minimize negative)
    c = [-v for v in values]
    
    # Constraint: sum(weights * x) <= capacity
    A_ub = [weights]
    b_ub = [capacity]
    
    # Binary variables (0 or 1)
    constraints = LinearConstraint(A_ub, -np.inf, b_ub)
    bounds = Bounds(lb=0, ub=1)
    
    # Integer constraints
    from scipy.optimize import LinearConstraint, Bounds
    integrality = [1] * len(values)  # 1 means integer
    
    result = milp(c, constraints=constraints, bounds=bounds, 
                  integrality=integrality)
    
    selected = [i for i, x in enumerate(result.x) if x > 0.5]
    print(f"Selected items: {selected}")
    print(f"Total value: {-result.fun}")
    print(f"Total weight: {sum(weights[i] for i in selected)}")

class GeneticAlgorithm:
    """Genetic algorithm for optimization."""
    
    def __init__(self, population_size: int = 100, 
                 mutation_rate: float = 0.01,
                 crossover_rate: float = 0.8,
                 generations: int = 100):
        self.population_size = population_size
        self.mutation_rate = mutation_rate
        self.crossover_rate = crossover_rate
        self.generations = generations
        
    def create_individual(self) -> np.ndarray:
        """Override to define individual representation."""
        raise NotImplementedError
        
    def fitness(self, individual: np.ndarray) -> float:
        """Override to define fitness function."""
        raise NotImplementedError
        
    def crossover(self, parent1: np.ndarray, parent2: np.ndarray) -> np.ndarray:
        """Single-point crossover."""
        if np.random.random() > self.crossover_rate:
            return parent1.copy()
        
        point = np.random.randint(1, len(parent1))
        child = np.concatenate([parent1[:point], parent2[point:]])
        return child
    
    def mutate(self, individual: np.ndarray) -> np.ndarray:
        """Random mutation."""
        mutated = individual.copy()
        for i in range(len(mutated)):
            if np.random.random() < self.mutation_rate:
                mutated[i] += np.random.normal(0, 0.1)
        return mutated
    
    def selection(self, population: list, fitness: np.ndarray) -> np.ndarray:
        """Tournament selection."""
        tournament_size = 3
        selected = []
        
        for _ in range(len(population)):
            tournament_idx = np.random.choice(len(population), tournament_size, replace=False)
            tournament_fitness = fitness[tournament_idx]
            winner_idx = tournament_idx[np.argmax(tournament_fitness)]
            selected.append(population[winner_idx])
        
        return np.array(selected)
    
    def run(self) -> tuple:
        """Run genetic algorithm."""
        # Initialize population
        population = np.array([self.create_individual() 
                              for _ in range(self.population_size)])
        
        best_individual = None
        best_fitness = -np.inf
        
        for generation in range(self.generations):
            # Evaluate fitness
            fitness = np.array([self.fitness(ind) for ind in population])
            
            # Track best
            if np.max(fitness) > best_fitness:
                best_fitness = np.max(fitness)
                best_individual = population[np.argmax(fitness)].copy()
            
            # Selection
            selected = self.selection(population, fitness)
            
            # Create next generation
            new_population = []
            for i in range(0, len(selected), 2):
                if i + 1 < len(selected):
                    child = self.crossover(selected[i], selected[i+1])
                    child = self.mutate(child)
                    new_population.append(child)
                    
            population = np.array(new_population[:self.population_size])
        
        return best_individual, best_fitness


class KnapsackGA(GeneticAlgorithm):
    """Genetic algorithm for knapsack problem."""
    
    def __init__(self, values, weights, capacity, **kwargs):
        super().__init__(**kwargs)
        self.values = values
        self.weights = weights
        self.capacity = capacity
        
    def create_individual(self) -> np.ndarray:
        return np.random.randint(0, 2, len(self.values))
    
    def fitness(self, individual: np.ndarray) -> float:
        total_value = np.sum(individual * self.values)
        total_weight = np.sum(individual * self.weights)
        
        if total_weight > self.capacity:
            return 0
        
        return total_value

def k_fold_cross_validation(
    model, X: np.ndarray, y: np.ndarray, k: int = 5
) -> dict:
    """K-fold cross-validation."""
    n = len(y)
    indices = np.random.permutation(n)
    fold_size = n // k
    
    scores = []
    
    for fold in range(k):
        # Split data
        val_start = fold * fold_size
        val_end = val_start + fold_size if fold < k - 1 else n
        
        val_idx = indices[val_start:val_end]
        train_idx = np.concatenate([indices[:val_start], indices[val_end:]])
        
        X_train, X_val = X[train_idx], X[val_idx]
        y_train, y_val = y[train_idx], y[val_idx]
        
        # Train and evaluate
        model.fit(X_train, y_train)
        predictions = model.predict(X_val)
        
        # Calculate score (MSE for regression)
        mse = np.mean((predictions - y_val) ** 2)
        scores.append(mse)
    
    return {
        'mean_mse': np.mean(scores),
        'std_mse': np.std(scores),
        'scores': scores
    }


def bootstrap_validation(
    model, X: np.ndarray, y: np.ndarray, n_bootstrap: int = 1000
) -> dict:
    """Bootstrap validation for confidence intervals."""
    n = len(y)
    bootstrap_scores = []
    
    for _ in range(n_bootstrap):
        # Sample with replacement
        idx = np.random.choice(n, n, replace=True)
        X_boot, y_boot = X[idx], y[idx]
        
        # Train model
        model.fit(X_boot, y_boot)
        predictions = model.predict(X_boot)
        
        # Score
        mse = np.mean((predictions - y_boot) ** 2)
        bootstrap_scores.append(mse)
    
    return {
        'mean_mse': np.mean(bootstrap_scores),
        'ci_lower': np.percentile(bootstrap_scores, 2.5),
        'ci_upper': np.percentile(bootstrap_scores, 97.5)
    }