Nonlinear Solvers

#!/usr/bin/env python3
"""Analyze residual history to classify convergence type."""
import argparse
import json
import math
import sys
from typing import Any, Dict, List, Optional


def analyze_convergence(
    residuals: List[float],
    tolerance: float = 1e-10,
) -> Dict[str, Any]:
    """Analyze residual history to classify convergence behavior.

    Args:
        residuals: List of residual norms from solver iterations
        tolerance: Convergence tolerance

    Returns:
        Dictionary with convergence analysis results
    """
    if not residuals:
        raise ValueError("residuals must not be empty")

    if any(r < 0 for r in residuals):
        raise ValueError("residuals must be non-negative")

    if tolerance <= 0:
        raise ValueError("tolerance must be positive")

    n = len(residuals)
    final_residual = residuals[-1]
    converged = final_residual <= tolerance

    # Handle single residual case
    if n == 1:
        return {
            "converged": converged,
            "iterations": n,
            "final_residual": final_residual,
            "convergence_type": "unknown",
            "estimated_rate": None,
            "diagnosis": "Insufficient data for convergence analysis.",
            "recommended_action": "Continue iterations to gather more data.",
        }

    # Check for divergence
    if n >= 2 and residuals[-1] > residuals[0] * 1.1:
        return {
            "converged": False,
            "iterations": n,
            "final_residual": final_residual,
            "convergence_type": "diverged",
            "estimated_rate": None,
            "diagnosis": "Residuals are increasing; solver is diverging.",
            "recommended_action": "Reduce step size, add damping, or check Jacobian accuracy.",
        }

    # Check for stagnation
    if n >= 3:
        recent = residuals[-3:]
        if len(recent) == 3:
            rel_change = abs(recent[-1] - recent[0]) / (abs(recent[0]) + 1e-30)
            if rel_change < 0.01 and not converged:
                return {
                    "converged": False,
                    "iterations": n,
                    "final_residual": final_residual,
                    "convergence_type": "stagnated",
                    "estimated_rate": None,
                    "diagnosis": "Residual has stagnated without reaching tolerance.",
                    "recommended_action": "Improve preconditioner, check for near-singularity, or use different solver.",
                }

    # Estimate convergence rate from log-residuals
    rates = []
    for i in range(1, n):
        if residuals[i - 1] > 1e-30 and residuals[i] > 1e-30:
            ratio = residuals[i] / residuals[i - 1]
            if ratio > 0 and ratio < 1:
                rates.append(ratio)

    if not rates:
        convergence_type = "unknown"
        estimated_rate = None
    else:
        avg_rate = sum(rates) / len(rates)
        estimated_rate = avg_rate

        # Try to detect quadratic convergence by looking at log-log behavior
        # For quadratic: log(r_{k+1}) ≈ 2 * log(r_k)
        if n >= 4 and all(r > 1e-30 for r in residuals):
            log_residuals = [math.log10(r) for r in residuals if r > 0]
            if len(log_residuals) >= 4:
                # Check if log(r_{k+1}) / log(r_k) ≈ 2 for quadratic
                log_ratios = []
                for i in range(1, len(log_residuals)):
                    if abs(log_residuals[i - 1]) > 0.1:
                        log_ratios.append(log_residuals[i] / log_residuals[i - 1])

                if log_ratios and all(1.5 < r < 2.5 for r in log_ratios[-3:]):
                    convergence_type = "quadratic"
                elif avg_rate < 0.3:
                    convergence_type = "superlinear"
                elif avg_rate < 0.9:
                    convergence_type = "linear"
                else:
                    convergence_type = "sublinear"
            else:
                if avg_rate < 0.3:
                    convergence_type = "superlinear"
                elif avg_rate < 0.9:
                    convergence_type = "linear"
                else:
                    convergence_type = "sublinear"
        else:
            if avg_rate < 0.3:
                convergence_type = "superlinear"
            elif avg_rate < 0.9:
                convergence_type = "linear"
            else:
                convergence_type = "sublinear"

    # Generate diagnosis and recommendation
    diagnosis_map = {
        "quadratic": "Quadratic convergence indicates optimal Newton behavior.",
        "superlinear": "Superlinear convergence; quasi-Newton methods working well.",
        "linear": "Linear convergence; rate is acceptable but could be improved.",
        "sublinear": "Sublinear convergence; solver is making slow progress.",
        "unknown": "Could not determine convergence type.",
    }

    action_map = {
        "quadratic": "Continue with current solver; convergence is optimal.",
        "superlinear": "Current setup is effective; monitor for stagnation.",
        "linear": "Consider stronger preconditioner or switch to Newton if Jacobian available.",
        "sublinear": "Switch to Newton method, improve globalization, or check problem formulation.",
        "unknown": "Gather more iterations for analysis.",
    }

    diagnosis = diagnosis_map.get(convergence_type, "Unknown convergence behavior.")
    recommended_action = action_map.get(convergence_type, "Review solver configuration.")

    if converged:
        recommended_action = "Solver converged successfully."

    return {
        "converged": converged,
        "iterations": n,
        "final_residual": final_residual,
        "convergence_type": convergence_type,
        "estimated_rate": estimated_rate,
        "diagnosis": diagnosis,
        "recommended_action": recommended_action,
    }


def parse_args() -> argparse.Namespace:
    parser = argparse.ArgumentParser(
        description="Analyze convergence from residual history.",
        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
    )
    parser.add_argument(
        "--residuals",
        type=str,
        required=True,
        help="Comma-separated residual values",
    )
    parser.add_argument(
        "--tolerance",
        type=float,
        default=1e-10,
        help="Convergence tolerance",
    )
    parser.add_argument("--json", action="store_true", help="Emit JSON output")
    return parser.parse_args()


def main() -> None:
    args = parse_args()

    try:
        residuals = [float(x.strip()) for x in args.residuals.split(",")]
    except ValueError:
        print("Error: residuals must be comma-separated numbers", file=sys.stderr)
        sys.exit(2)

    try:
        result = analyze_convergence(
            residuals=residuals,
            tolerance=args.tolerance,
        )
    except ValueError as exc:
        print(str(exc), file=sys.stderr)
        sys.exit(2)

    payload: Dict[str, Any] = {
        "inputs": {
            "residuals": residuals,
            "tolerance": args.tolerance,
        },
        "results": result,
    }

    if args.json:
        print(json.dumps(payload, indent=2, sort_keys=True))
        return

    print("Convergence analysis")
    print(f"  converged: {result['converged']}")
    print(f"  iterations: {result['iterations']}")
    print(f"  final_residual: {result['final_residual']:.2e}")
    print(f"  convergence_type: {result['convergence_type']}")
    if result["estimated_rate"] is not None:
        print(f"  estimated_rate: {result['estimated_rate']:.4f}")
    print(f"  diagnosis: {result['diagnosis']}")
    print(f"  recommended_action: {result['recommended_action']}")


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""Recommend line search vs trust region globalization strategy."""
import argparse
import json
import sys
from typing import Any, Dict, List


def advise_globalization(
    problem_type: str,
    jacobian_quality: str,
    previous_failures: int,
    oscillating_residual: bool,
    step_rejection_rate: float,
) -> Dict[str, Any]:
    """Recommend globalization strategy for nonlinear solvers.

    Args:
        problem_type: Type of problem (root-finding, optimization, least-squares)
        jacobian_quality: Quality of Jacobian (good, ill-conditioned, near-singular)
        previous_failures: Number of previous solver failures
        oscillating_residual: Whether residual history shows oscillation
        step_rejection_rate: Fraction of rejected steps (0.0 to 1.0)

    Returns:
        Dictionary with globalization strategy recommendations
    """
    valid_problem_types = {"root-finding", "optimization", "least-squares"}
    if problem_type not in valid_problem_types:
        raise ValueError(f"problem_type must be one of {valid_problem_types}")

    valid_jacobian_quality = {"good", "ill-conditioned", "near-singular"}
    if jacobian_quality not in valid_jacobian_quality:
        raise ValueError(f"jacobian_quality must be one of {valid_jacobian_quality}")

    if previous_failures < 0:
        raise ValueError("previous_failures must be non-negative")

    if not 0.0 <= step_rejection_rate <= 1.0:
        raise ValueError("step_rejection_rate must be between 0.0 and 1.0")

    notes: List[str] = []

    # Decision logic for strategy selection
    use_trust_region = False

    # Trust region preferred for ill-conditioned or near-singular Jacobians
    if jacobian_quality in {"ill-conditioned", "near-singular"}:
        use_trust_region = True
        notes.append("Trust region provides better stability for poor Jacobian quality.")

    # Trust region preferred for high failure rate
    if previous_failures >= 2:
        use_trust_region = True
        notes.append("Multiple failures suggest trust region for more robust steps.")

    # Trust region preferred for oscillating residuals
    if oscillating_residual:
        use_trust_region = True
        notes.append("Oscillating residuals indicate step size issues; trust region helps.")

    # Trust region preferred for high step rejection
    if step_rejection_rate > 0.3:
        use_trust_region = True
        notes.append("High step rejection rate favors trust region approach.")

    # Problem-type specific adjustments
    if problem_type == "least-squares":
        use_trust_region = True
        notes.append("Trust region is standard for least-squares (Levenberg-Marquardt).")

    if use_trust_region:
        strategy = "trust-region"

        # Select trust region type
        if jacobian_quality == "near-singular":
            trust_region_type = "Levenberg-Marquardt"
            notes.append("LM regularization handles near-singular Jacobian.")
        elif problem_type == "optimization":
            trust_region_type = "Steihaug-CG"
            notes.append("Steihaug-CG efficient for large-scale optimization.")
        else:
            trust_region_type = "dogleg"
            notes.append("Dogleg combines Cauchy and Newton steps effectively.")

        # Initial damping parameter
        if jacobian_quality == "near-singular":
            initial_damping = 1.0
        elif jacobian_quality == "ill-conditioned":
            initial_damping = 0.1
        else:
            initial_damping = 0.01

        parameters = {
            "initial_radius": 1.0,
            "max_radius": 100.0,
            "eta1": 0.25,
            "eta2": 0.75,
            "gamma1": 0.25,
            "gamma2": 2.0,
        }

        return {
            "strategy": strategy,
            "line_search_type": None,
            "trust_region_type": trust_region_type,
            "initial_damping": initial_damping,
            "parameters": parameters,
            "notes": notes,
        }

    else:
        strategy = "line-search"

        # Select line search type
        if problem_type == "optimization":
            line_search_type = "Wolfe"
            notes.append("Wolfe conditions ensure sufficient decrease and curvature.")
        elif previous_failures > 0:
            line_search_type = "backtracking"
            notes.append("Backtracking is simple and robust after failures.")
        else:
            line_search_type = "Armijo"
            notes.append("Armijo sufficient decrease condition for root-finding.")

        # Parameters for line search
        if jacobian_quality == "ill-conditioned":
            initial_damping = 0.5
            notes.append("Starting with reduced step for ill-conditioned problem.")
        else:
            initial_damping = 1.0

        parameters = {
            "c1": 1e-4,  # Armijo/Wolfe sufficient decrease
            "c2": 0.9,  # Wolfe curvature condition
            "alpha_init": 1.0,
            "rho": 0.5,  # Backtracking factor
            "max_backtracks": 20,
        }

        return {
            "strategy": strategy,
            "line_search_type": line_search_type,
            "trust_region_type": None,
            "initial_damping": initial_damping,
            "parameters": parameters,
            "notes": notes,
        }


def parse_args() -> argparse.Namespace:
    parser = argparse.ArgumentParser(
        description="Recommend globalization strategy for nonlinear solvers.",
        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
    )
    parser.add_argument(
        "--problem-type",
        type=str,
        required=True,
        choices=["root-finding", "optimization", "least-squares"],
        help="Type of problem",
    )
    parser.add_argument(
        "--jacobian-quality",
        type=str,
        default="good",
        choices=["good", "ill-conditioned", "near-singular"],
        help="Quality of Jacobian",
    )
    parser.add_argument(
        "--previous-failures",
        type=int,
        default=0,
        help="Number of previous solver failures",
    )
    parser.add_argument(
        "--oscillating",
        action="store_true",
        help="Residual history shows oscillation",
    )
    parser.add_argument(
        "--step-rejection-rate",
        type=float,
        default=0.0,
        help="Fraction of rejected steps (0.0 to 1.0)",
    )
    parser.add_argument("--json", action="store_true", help="Emit JSON output")
    return parser.parse_args()


def main() -> None:
    args = parse_args()

    try:
        result = advise_globalization(
            problem_type=args.problem_type,
            jacobian_quality=args.jacobian_quality,
            previous_failures=args.previous_failures,
            oscillating_residual=args.oscillating,
            step_rejection_rate=args.step_rejection_rate,
        )
    except ValueError as exc:
        print(str(exc), file=sys.stderr)
        sys.exit(2)

    payload: Dict[str, Any] = {
        "inputs": {
            "problem_type": args.problem_type,
            "jacobian_quality": args.jacobian_quality,
            "previous_failures": args.previous_failures,
            "oscillating_residual": args.oscillating,
            "step_rejection_rate": args.step_rejection_rate,
        },
        "results": result,
    }

    if args.json:
        print(json.dumps(payload, indent=2, sort_keys=True))
        return

    print("Globalization strategy")
    print(f"  strategy: {result['strategy']}")
    if result["line_search_type"]:
        print(f"  line_search_type: {result['line_search_type']}")
    if result["trust_region_type"]:
        print(f"  trust_region_type: {result['trust_region_type']}")
    print(f"  initial_damping: {result['initial_damping']}")
    for note in result["notes"]:
        print(f"  note: {note}")


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""Analyze Jacobian matrix quality for nonlinear solvers."""
import argparse
import json
import sys
from typing import Any, Dict, List, Optional

import numpy as np


def diagnose_jacobian(
    matrix: np.ndarray,
    finite_diff_matrix: Optional[np.ndarray] = None,
    tolerance: float = 1e-6,
) -> Dict[str, Any]:
    """Analyze Jacobian matrix quality.

    Args:
        matrix: The Jacobian matrix to analyze
        finite_diff_matrix: Optional finite-difference approximation for comparison
        tolerance: Tolerance for rank deficiency detection

    Returns:
        Dictionary with Jacobian diagnostics
    """
    if matrix.ndim != 2:
        raise ValueError("matrix must be 2-dimensional")

    if matrix.size == 0:
        raise ValueError("matrix must not be empty")

    if tolerance <= 0:
        raise ValueError("tolerance must be positive")

    m, n = matrix.shape
    notes: List[str] = []

    # Compute SVD for condition number and rank analysis
    try:
        singular_values = np.linalg.svd(matrix, compute_uv=False)
    except np.linalg.LinAlgError:
        return {
            "shape": [m, n],
            "condition_number": float("inf"),
            "rank_deficient": True,
            "estimated_rank": 0,
            "singular_value_min": 0.0,
            "singular_value_max": 0.0,
            "jacobian_quality": "singular",
            "finite_diff_error": None,
            "notes": ["SVD computation failed; matrix may be ill-formed."],
        }

    sv_max = float(singular_values[0])
    sv_min = float(singular_values[-1])

    # Condition number
    if sv_min > 1e-30:
        condition_number = sv_max / sv_min
    else:
        condition_number = float("inf")

    # Estimate numerical rank
    rank_tol = max(m, n) * sv_max * np.finfo(float).eps
    estimated_rank = int(np.sum(singular_values > rank_tol))
    rank_deficient = estimated_rank < min(m, n)

    # Classify Jacobian quality
    if condition_number == float("inf") or sv_min < 1e-14:
        jacobian_quality = "near-singular"
        notes.append("Near-singular Jacobian; regularization may be needed.")
    elif condition_number > 1e10:
        jacobian_quality = "ill-conditioned"
        notes.append("Highly ill-conditioned; use iterative refinement or scaling.")
    elif condition_number > 1e6:
        jacobian_quality = "moderately-conditioned"
        notes.append("Moderate conditioning; standard methods should work.")
    else:
        jacobian_quality = "good"
        notes.append("Well-conditioned Jacobian.")

    if rank_deficient:
        notes.append(f"Rank deficient: estimated rank {estimated_rank} < min({m}, {n}).")

    # Compare with finite difference approximation if provided
    finite_diff_error = None
    if finite_diff_matrix is not None:
        if finite_diff_matrix.shape != matrix.shape:
            notes.append("Finite-diff matrix shape mismatch; skipping comparison.")
        else:
            diff = matrix - finite_diff_matrix
            relative_error = np.linalg.norm(diff) / (np.linalg.norm(matrix) + 1e-30)
            finite_diff_error = float(relative_error)

            if relative_error > 0.1:
                notes.append(f"Large discrepancy with finite-diff ({relative_error:.2e}); check analytic Jacobian.")
            elif relative_error > 0.01:
                notes.append(f"Moderate discrepancy with finite-diff ({relative_error:.2e}).")
            else:
                notes.append("Jacobian matches finite-difference approximation well.")

    return {
        "shape": [m, n],
        "condition_number": condition_number,
        "rank_deficient": rank_deficient,
        "estimated_rank": estimated_rank,
        "singular_value_min": sv_min,
        "singular_value_max": sv_max,
        "jacobian_quality": jacobian_quality,
        "finite_diff_error": finite_diff_error,
        "notes": notes,
    }


def load_matrix(path: str) -> np.ndarray:
    """Load matrix from text file."""
    try:
        return np.loadtxt(path)
    except Exception as e:
        raise ValueError(f"Failed to load matrix from {path}: {e}")


def parse_args() -> argparse.Namespace:
    parser = argparse.ArgumentParser(
        description="Analyze Jacobian matrix quality.",
        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
    )
    parser.add_argument(
        "--matrix",
        type=str,
        required=True,
        help="Path to Jacobian matrix file (text format)",
    )
    parser.add_argument(
        "--finite-diff-matrix",
        type=str,
        default=None,
        help="Path to finite-difference Jacobian for comparison",
    )
    parser.add_argument(
        "--tolerance",
        type=float,
        default=1e-6,
        help="Tolerance for rank deficiency detection",
    )
    parser.add_argument("--json", action="store_true", help="Emit JSON output")
    return parser.parse_args()


def main() -> None:
    args = parse_args()

    try:
        matrix = load_matrix(args.matrix)
        if matrix.ndim == 1:
            matrix = matrix.reshape(1, -1)
    except ValueError as exc:
        print(str(exc), file=sys.stderr)
        sys.exit(2)

    finite_diff_matrix = None
    if args.finite_diff_matrix:
        try:
            finite_diff_matrix = load_matrix(args.finite_diff_matrix)
            if finite_diff_matrix.ndim == 1:
                finite_diff_matrix = finite_diff_matrix.reshape(1, -1)
        except ValueError as exc:
            print(str(exc), file=sys.stderr)
            sys.exit(2)

    try:
        result = diagnose_jacobian(
            matrix=matrix,
            finite_diff_matrix=finite_diff_matrix,
            tolerance=args.tolerance,
        )
    except ValueError as exc:
        print(str(exc), file=sys.stderr)
        sys.exit(2)

    payload: Dict[str, Any] = {
        "inputs": {
            "matrix": args.matrix,
            "finite_diff_matrix": args.finite_diff_matrix,
            "tolerance": args.tolerance,
        },
        "results": result,
    }

    if args.json:
        print(json.dumps(payload, indent=2, sort_keys=True))
        return

    print("Jacobian diagnostics")
    print(f"  shape: {result['shape']}")
    print(f"  condition_number: {result['condition_number']:.2e}")
    print(f"  rank_deficient: {result['rank_deficient']}")
    print(f"  estimated_rank: {result['estimated_rank']}")
    print(f"  singular_value_min: {result['singular_value_min']:.2e}")
    print(f"  singular_value_max: {result['singular_value_max']:.2e}")
    print(f"  jacobian_quality: {result['jacobian_quality']}")
    if result["finite_diff_error"] is not None:
        print(f"  finite_diff_error: {result['finite_diff_error']:.2e}")
    for note in result["notes"]:
        print(f"  note: {note}")


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""Monitor residual patterns and detect failure modes."""
import argparse
import json
import sys
from typing import Any, Dict, List, Optional


def monitor_residuals(
    residuals: List[float],
    function_evals: Optional[List[int]] = None,
    step_sizes: Optional[List[float]] = None,
    target_tolerance: float = 1e-10,
) -> Dict[str, Any]:
    """Monitor residual patterns and detect failure modes.

    Args:
        residuals: List of residual norms
        function_evals: Optional list of cumulative function evaluations
        step_sizes: Optional list of step sizes taken
        target_tolerance: Target convergence tolerance

    Returns:
        Dictionary with monitoring results and recommendations
    """
    if not residuals:
        raise ValueError("residuals must not be empty")

    if any(r < 0 for r in residuals):
        raise ValueError("residuals must be non-negative")

    if target_tolerance <= 0:
        raise ValueError("target_tolerance must be positive")

    n = len(residuals)
    patterns_detected: List[str] = []
    alerts: List[str] = []
    recommendations: List[str] = []

    # Basic stats
    initial_residual = residuals[0]
    final_residual = residuals[-1]

    if initial_residual > 1e-30:
        residual_reduction = final_residual / initial_residual
    else:
        residual_reduction = 1.0

    # Check for convergence
    converged = final_residual <= target_tolerance

    # Function evaluation efficiency
    total_function_evals = None
    efficiency = None
    if function_evals and len(function_evals) == n:
        total_function_evals = function_evals[-1]
        if total_function_evals > 0:
            # Orders of magnitude reduced per function eval
            if residual_reduction > 0 and residual_reduction < 1:
                import math
                log_reduction = -math.log10(residual_reduction)
                efficiency = log_reduction / total_function_evals
            else:
                efficiency = 0.0

    # Step size analysis
    step_size_trend = None
    if step_sizes and len(step_sizes) >= 2:
        avg_early = sum(step_sizes[:len(step_sizes)//2 + 1]) / (len(step_sizes)//2 + 1)
        avg_late = sum(step_sizes[len(step_sizes)//2:]) / (len(step_sizes) - len(step_sizes)//2)

        if avg_early > 1e-30:
            ratio = avg_late / avg_early
            if ratio < 0.5:
                step_size_trend = "decreasing"
                patterns_detected.append("step_size_decreasing")
            elif ratio > 2.0:
                step_size_trend = "increasing"
                patterns_detected.append("step_size_increasing")
            else:
                step_size_trend = "stable"

    # Pattern detection
    if n >= 3:
        # Check for oscillation
        oscillations = 0
        for i in range(1, n - 1):
            if (residuals[i] > residuals[i-1] and residuals[i] > residuals[i+1]) or \
               (residuals[i] < residuals[i-1] and residuals[i] < residuals[i+1]):
                oscillations += 1

        if oscillations >= n // 3:
            patterns_detected.append("oscillating")
            alerts.append("Residual oscillation detected; may indicate step size issues.")
            recommendations.append("Consider trust region or reduced initial step size.")

        # Check for plateau
        recent = residuals[-min(5, n):]
        if len(recent) >= 3:
            rel_change = max(recent) / (min(recent) + 1e-30)
            if rel_change < 1.1 and not converged:
                patterns_detected.append("plateau")
                alerts.append("Residual plateau detected; solver may be stagnating.")
                recommendations.append("Try stronger preconditioner or different solver.")

        # Check for divergence
        if n >= 2 and residuals[-1] > residuals[0] * 1.5:
            patterns_detected.append("diverging")
            alerts.append("Residuals are increasing; solver is diverging!")
            recommendations.append("Reduce step size, check Jacobian, or add regularization.")

        # Check for very slow progress
        if n >= 10:
            halfway = n // 2
            early_reduction = residuals[halfway] / (residuals[0] + 1e-30)
            late_reduction = residuals[-1] / (residuals[halfway] + 1e-30)

            if early_reduction > 0.9 and late_reduction > 0.9:
                patterns_detected.append("slow_convergence")
                alerts.append("Very slow convergence detected.")
                recommendations.append("Consider Newton method with good Jacobian.")

    # Check for initial spike (common with bad initial guess)
    if n >= 2 and residuals[1] > residuals[0] * 2:
        patterns_detected.append("initial_spike")
        recommendations.append("Consider better initial guess or damped first step.")

    # Final recommendations
    if not patterns_detected:
        patterns_detected.append("normal")

    if converged:
        recommendations.insert(0, "Solver converged successfully.")
    elif not recommendations:
        recommendations.append("Continue monitoring; no immediate issues detected.")

    result: Dict[str, Any] = {
        "residual_reduction": residual_reduction,
        "iterations": n,
        "converged": converged,
        "patterns_detected": patterns_detected,
        "alerts": alerts,
        "recommendations": recommendations,
    }

    if total_function_evals is not None:
        result["function_evals"] = total_function_evals
    if efficiency is not None:
        result["efficiency"] = efficiency
    if step_size_trend is not None:
        result["step_size_trend"] = step_size_trend

    return result


def parse_args() -> argparse.Namespace:
    parser = argparse.ArgumentParser(
        description="Monitor residual patterns and detect failure modes.",
        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
    )
    parser.add_argument(
        "--residuals",
        type=str,
        required=True,
        help="Comma-separated residual values",
    )
    parser.add_argument(
        "--function-evals",
        type=str,
        default=None,
        help="Comma-separated cumulative function evaluation counts",
    )
    parser.add_argument(
        "--step-sizes",
        type=str,
        default=None,
        help="Comma-separated step sizes",
    )
    parser.add_argument(
        "--target-tolerance",
        type=float,
        default=1e-10,
        help="Target convergence tolerance",
    )
    parser.add_argument("--json", action="store_true", help="Emit JSON output")
    return parser.parse_args()


def main() -> None:
    args = parse_args()

    try:
        residuals = [float(x.strip()) for x in args.residuals.split(",")]
    except ValueError:
        print("Error: residuals must be comma-separated numbers", file=sys.stderr)
        sys.exit(2)

    function_evals = None
    if args.function_evals:
        try:
            function_evals = [int(x.strip()) for x in args.function_evals.split(",")]
        except ValueError:
            print("Error: function-evals must be comma-separated integers", file=sys.stderr)
            sys.exit(2)

    step_sizes = None
    if args.step_sizes:
        try:
            step_sizes = [float(x.strip()) for x in args.step_sizes.split(",")]
        except ValueError:
            print("Error: step-sizes must be comma-separated numbers", file=sys.stderr)
            sys.exit(2)

    try:
        result = monitor_residuals(
            residuals=residuals,
            function_evals=function_evals,
            step_sizes=step_sizes,
            target_tolerance=args.target_tolerance,
        )
    except ValueError as exc:
        print(str(exc), file=sys.stderr)
        sys.exit(2)

    payload: Dict[str, Any] = {
        "inputs": {
            "residuals": residuals,
            "function_evals": function_evals,
            "step_sizes": step_sizes,
            "target_tolerance": args.target_tolerance,
        },
        "results": result,
    }

    if args.json:
        print(json.dumps(payload, indent=2, sort_keys=True))
        return

    print("Residual monitoring")
    print(f"  residual_reduction: {result['residual_reduction']:.2e}")
    print(f"  iterations: {result['iterations']}")
    print(f"  converged: {result['converged']}")
    print(f"  patterns: {', '.join(result['patterns_detected'])}")
    for alert in result["alerts"]:
        print(f"  ALERT: {alert}")
    for rec in result["recommendations"]:
        print(f"  recommendation: {rec}")


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""Select nonlinear solver based on problem characteristics."""
import argparse
import json
import sys
from typing import Any, Dict, List


def select_solver(
    jacobian_available: bool,
    jacobian_expensive: bool,
    problem_size: int,
    spd_hessian: bool,
    smooth_objective: bool,
    constraint_type: str,
    memory_limited: bool,
    high_accuracy: bool,
) -> Dict[str, List[str]]:
    """Select nonlinear solver based on problem characteristics.

    Args:
        jacobian_available: Whether analytic Jacobian is available
        jacobian_expensive: Whether Jacobian computation is expensive
        problem_size: Number of unknowns
        spd_hessian: Whether Hessian is symmetric positive definite
        smooth_objective: Whether objective/residual is smooth
        constraint_type: Type of constraints (none, bound, equality, inequality)
        memory_limited: Whether memory is constrained
        high_accuracy: Whether high accuracy is required

    Returns:
        Dictionary with recommended solvers, alternatives, and notes
    """
    if problem_size <= 0:
        raise ValueError("problem_size must be positive")

    valid_constraints = {"none", "bound", "equality", "inequality"}
    if constraint_type not in valid_constraints:
        raise ValueError(f"constraint_type must be one of {valid_constraints}")

    recommended: List[str] = []
    alternatives: List[str] = []
    notes: List[str] = []

    large_problem = problem_size >= 10_000

    # Handle constrained optimization
    if constraint_type == "inequality":
        recommended.append("SQP (Sequential Quadratic Programming)")
        alternatives.append("Interior Point")
        notes.append("Inequality constraints require specialized solvers.")
        if not jacobian_available:
            notes.append("Consider finite-difference Jacobian or quasi-Newton Hessian.")
        return {"recommended": recommended, "alternatives": alternatives, "notes": notes}

    if constraint_type == "equality":
        recommended.append("SQP")
        alternatives.append("Augmented Lagrangian")
        notes.append("Equality constraints: use Lagrange multipliers or penalty methods.")
        return {"recommended": recommended, "alternatives": alternatives, "notes": notes}

    if constraint_type == "bound":
        if spd_hessian:
            recommended.append("L-BFGS-B")
            alternatives.append("Trust-Region Reflective")
        else:
            recommended.append("Trust-Region Reflective")
            alternatives.append("L-BFGS-B")
        notes.append("Bound constraints handled via projected methods.")
        return {"recommended": recommended, "alternatives": alternatives, "notes": notes}

    # Unconstrained case
    if jacobian_available and not jacobian_expensive:
        if high_accuracy:
            recommended.append("Newton (full)")
            alternatives.append("Modified Newton")
            notes.append("Full Newton provides quadratic convergence near solution.")
        elif large_problem:
            recommended.append("Newton-Krylov (GMRES)")
            alternatives.append("Newton-Krylov (BiCGSTAB)")
            notes.append("Newton-Krylov avoids forming full Jacobian for large problems.")
        else:
            recommended.append("Newton (full)")
            alternatives.append("Modified Newton")
    elif jacobian_available and jacobian_expensive:
        if memory_limited:
            recommended.append("L-BFGS")
            alternatives.append("Broyden")
            notes.append("L-BFGS uses limited memory quasi-Newton updates.")
        else:
            recommended.append("Modified Newton")
            alternatives.append("Broyden")
            notes.append("Modified Newton reuses Jacobian for multiple iterations.")
    else:
        # No Jacobian available
        if smooth_objective:
            if memory_limited or large_problem:
                recommended.append("L-BFGS")
                alternatives.append("Broyden (good)")
                notes.append("Quasi-Newton methods build approximate Jacobian.")
            else:
                recommended.append("BFGS")
                alternatives.append("SR1")
                notes.append("BFGS provides superlinear convergence for smooth problems.")
        else:
            recommended.append("Anderson Acceleration")
            alternatives.append("Picard (fixed-point)")
            notes.append("Non-smooth problems may benefit from fixed-point methods.")

    # Additional recommendations based on problem characteristics
    if spd_hessian and not jacobian_available:
        notes.append("SPD Hessian: BFGS maintains positive definiteness.")

    if large_problem and "Newton (full)" in recommended:
        notes.append("For very large problems, consider Newton-Krylov instead.")

    if not smooth_objective:
        notes.append("Non-smooth: consider subgradient or bundle methods.")

    return {"recommended": recommended, "alternatives": alternatives, "notes": notes}


def parse_args() -> argparse.Namespace:
    parser = argparse.ArgumentParser(
        description="Select a nonlinear solver based on problem characteristics.",
        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
    )
    parser.add_argument(
        "--jacobian-available",
        action="store_true",
        help="Analytic Jacobian is available",
    )
    parser.add_argument(
        "--jacobian-expensive",
        action="store_true",
        help="Jacobian computation is expensive",
    )
    parser.add_argument(
        "--size",
        type=int,
        required=True,
        help="Problem size (number of unknowns)",
    )
    parser.add_argument(
        "--spd-hessian",
        action="store_true",
        help="Hessian is symmetric positive definite",
    )
    parser.add_argument(
        "--smooth",
        action="store_true",
        help="Objective/residual is smooth",
    )
    parser.add_argument(
        "--constraints",
        type=str,
        default="none",
        choices=["none", "bound", "equality", "inequality"],
        help="Type of constraints",
    )
    parser.add_argument(
        "--memory-limited",
        action="store_true",
        help="Memory is constrained",
    )
    parser.add_argument(
        "--high-accuracy",
        action="store_true",
        help="High accuracy is required",
    )
    parser.add_argument("--json", action="store_true", help="Emit JSON output")
    return parser.parse_args()


def main() -> None:
    args = parse_args()
    try:
        result = select_solver(
            jacobian_available=args.jacobian_available,
            jacobian_expensive=args.jacobian_expensive,
            problem_size=args.size,
            spd_hessian=args.spd_hessian,
            smooth_objective=args.smooth,
            constraint_type=args.constraints,
            memory_limited=args.memory_limited,
            high_accuracy=args.high_accuracy,
        )
    except ValueError as exc:
        print(str(exc), file=sys.stderr)
        sys.exit(2)

    payload: Dict[str, Any] = {
        "inputs": {
            "jacobian_available": args.jacobian_available,
            "jacobian_expensive": args.jacobian_expensive,
            "problem_size": args.size,
            "spd_hessian": args.spd_hessian,
            "smooth_objective": args.smooth,
            "constraint_type": args.constraints,
            "memory_limited": args.memory_limited,
            "high_accuracy": args.high_accuracy,
        },
        "results": result,
    }

    if args.json:
        print(json.dumps(payload, indent=2, sort_keys=True))
        return

    print("Solver selection")
    print(f"  recommended: {', '.join(result['recommended'])}")
    if result["alternatives"]:
        print(f"  alternatives: {', '.join(result['alternatives'])}")
    for note in result["notes"]:
        print(f"  note: {note}")


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""Evaluate Newton/quasi-Newton step quality for trust region decisions."""
import argparse
import json
import sys
from typing import Any, Dict, List, Optional


def evaluate_step(
    predicted_reduction: float,
    actual_reduction: float,
    step_norm: float,
    gradient_norm: float,
    trust_radius: Optional[float] = None,
) -> Dict[str, Any]:
    """Evaluate step quality for trust region management.

    Args:
        predicted_reduction: Model-predicted decrease in objective
        actual_reduction: Actual decrease in objective
        step_norm: Norm of the step taken
        gradient_norm: Norm of gradient at current point
        trust_radius: Current trust radius (if using trust region)

    Returns:
        Dictionary with step quality assessment and recommendations
    """
    if predicted_reduction < 0:
        raise ValueError("predicted_reduction must be non-negative (model should predict decrease)")

    if step_norm < 0:
        raise ValueError("step_norm must be non-negative")

    if gradient_norm < 0:
        raise ValueError("gradient_norm must be non-negative")

    if trust_radius is not None and trust_radius <= 0:
        raise ValueError("trust_radius must be positive")

    notes: List[str] = []

    # Compute the ratio (actual/predicted reduction)
    if predicted_reduction > 1e-30:
        ratio = actual_reduction / predicted_reduction
    else:
        # Near zero predicted reduction - step is near stationary point
        if abs(actual_reduction) < 1e-30:
            ratio = 1.0
            notes.append("Near stationary point; both reductions negligible.")
        else:
            ratio = float("inf") if actual_reduction > 0 else float("-inf")
            notes.append("Predicted reduction near zero but actual reduction nonzero.")

    # Classify step quality based on ratio
    # Standard thresholds: eta1 = 0.25, eta2 = 0.75
    if ratio < 0:
        step_quality = "very_poor"
        accept_step = False
        notes.append("Negative ratio: objective increased instead of decreased.")
    elif ratio < 0.1:
        step_quality = "poor"
        accept_step = False
        notes.append("Very poor agreement between model and actual behavior.")
    elif ratio < 0.25:
        step_quality = "marginal"
        accept_step = True
        notes.append("Marginal step quality; accepting but reducing trust radius.")
    elif ratio < 0.75:
        step_quality = "good"
        accept_step = True
        notes.append("Good agreement between model and actual behavior.")
    else:
        step_quality = "excellent"
        accept_step = True
        notes.append("Excellent step quality; model is accurate.")

    # Trust radius adjustment recommendation
    trust_radius_action = None
    suggested_trust_radius = None

    if trust_radius is not None:
        if ratio < 0.1:
            trust_radius_action = "shrink_aggressive"
            suggested_trust_radius = trust_radius * 0.25
            notes.append("Aggressively shrinking trust radius due to poor step.")
        elif ratio < 0.25:
            trust_radius_action = "shrink"
            suggested_trust_radius = trust_radius * 0.5
            notes.append("Shrinking trust radius.")
        elif ratio > 0.75 and step_norm >= 0.9 * trust_radius:
            trust_radius_action = "expand"
            suggested_trust_radius = min(trust_radius * 2.0, 100.0)
            notes.append("Expanding trust radius as step hit boundary.")
        else:
            trust_radius_action = "maintain"
            suggested_trust_radius = trust_radius
            notes.append("Maintaining current trust radius.")
    else:
        # Line search mode
        if ratio < 0.25:
            notes.append("Consider reducing initial step size for line search.")
        elif ratio > 0.9:
            notes.append("Full Newton step successful; line search may be unnecessary.")

    # Check for potential issues
    if step_norm < 1e-14:
        notes.append("Warning: step norm is nearly zero; may be at local minimum or saddle.")

    if gradient_norm < 1e-10 and not accept_step:
        notes.append("Near-zero gradient with rejected step; check for saddle point.")

    # Cauchy decrease check: should get at least some fraction of Cauchy decrease
    cauchy_decrease_expected = gradient_norm * min(step_norm, gradient_norm)
    if accept_step and actual_reduction < 0.1 * cauchy_decrease_expected and gradient_norm > 1e-10:
        notes.append("Step achieved less than expected Cauchy decrease.")

    return {
        "ratio": ratio if not isinstance(ratio, float) or ratio != float("inf") else None,
        "step_quality": step_quality,
        "accept_step": accept_step,
        "trust_radius_action": trust_radius_action,
        "suggested_trust_radius": suggested_trust_radius,
        "notes": notes,
    }


def parse_args() -> argparse.Namespace:
    parser = argparse.ArgumentParser(
        description="Evaluate step quality for trust region decisions.",
        formatter_class=argparse.ArgumentDefaultsHelpFormatter,
    )
    parser.add_argument(
        "--predicted-reduction",
        type=float,
        required=True,
        help="Model-predicted decrease in objective",
    )
    parser.add_argument(
        "--actual-reduction",
        type=float,
        required=True,
        help="Actual decrease in objective",
    )
    parser.add_argument(
        "--step-norm",
        type=float,
        required=True,
        help="Norm of the step taken",
    )
    parser.add_argument(
        "--gradient-norm",
        type=float,
        required=True,
        help="Norm of gradient at current point",
    )
    parser.add_argument(
        "--trust-radius",
        type=float,
        default=None,
        help="Current trust radius (optional)",
    )
    parser.add_argument("--json", action="store_true", help="Emit JSON output")
    return parser.parse_args()


def main() -> None:
    args = parse_args()

    try:
        result = evaluate_step(
            predicted_reduction=args.predicted_reduction,
            actual_reduction=args.actual_reduction,
            step_norm=args.step_norm,
            gradient_norm=args.gradient_norm,
            trust_radius=args.trust_radius,
        )
    except ValueError as exc:
        print(str(exc), file=sys.stderr)
        sys.exit(2)

    payload: Dict[str, Any] = {
        "inputs": {
            "predicted_reduction": args.predicted_reduction,
            "actual_reduction": args.actual_reduction,
            "step_norm": args.step_norm,
            "gradient_norm": args.gradient_norm,
            "trust_radius": args.trust_radius,
        },
        "results": result,
    }

    if args.json:
        print(json.dumps(payload, indent=2, sort_keys=True))
        return

    print("Step quality evaluation")
    if result["ratio"] is not None:
        print(f"  ratio: {result['ratio']:.4f}")
    print(f"  step_quality: {result['step_quality']}")
    print(f"  accept_step: {result['accept_step']}")
    if result["trust_radius_action"]:
        print(f"  trust_radius_action: {result['trust_radius_action']}")
    if result["suggested_trust_radius"] is not None:
        print(f"  suggested_trust_radius: {result['suggested_trust_radius']:.4f}")
    for note in result["notes"]:
        print(f"  note: {note}")


if __name__ == "__main__":
    main()