### Quick Example

Source: https://github.com/robin-vjc/nsopy/blob/master/README.md

A basic example demonstrating how to use nsopy to minimize a function.

```APIDOC
## Usage
### Quick Example

We seek to minimize a function obtained by taking the `max` over a set of affine functions. 
The feasible set considered is the set of non-negative real numbers, i.e., <img src="https://latex.codecogs.com/svg.latex?\mathbb{X}&space;=&space;\mathbb{R}^&plus;" title="\mathbb{X} = \mathbb{R}^+"/>
for which the projection operation is straightforward. 
<p align="center">
  <img src="docs/img/basic_example.png" alt="Example" width="65%" href="#"/>
</p>

It is straightforward to see that the optimum is at `x* = 2.25`; we can solve this optimization problem numerically as follows:

```python
import numpy as np

def oracle(x_k):
    # evaluation of the f_i components at x_k
    fi_x_k = [-2*x_k + 2,  -1.0/3*x_k + 1,  x_k - 2]

    f_x_k = max(fi_x_k)  # function value at x_k

    diff_fi = [-2, -1.0/3.0, 1]  # gradients of the components
    max_i = fi_x_k.index(f_x_k)
    # subgradient at x_k is the gradient of the active function component; cast as (1x1 dimensional) np.array
    diff_f_xk = np.array([diff_fi[max_i], ])  

    return 0, f_x_k, diff_f_xk

def projection_function(x_k):
    if x_k is 0:
        return np.array([0,])
    else:
        return np.maximum(x_k, 0)
```

Instantiation of method and logger, solve and print

```python
from nsopy.methods.subgradient import SubgradientMethod
from nsopy.loggers import GenericMethodLogger

method = SubgradientMethod(oracle, projection_function, stepsize_0=0.1, stepsize_rule='constant', sense='min')
logger = GenericMethodLogger(method)

for iteration in range(200):
    method.step()
```

Result:

```
>>> print(logger.x_k_iterates[-5:])

[2.1999999999999904, 2.216666666666657, 2.2333333333333236, 2.2499999999999902, 2.266666666666657]
```
```

--------------------------------

### Installation

Source: https://github.com/robin-vjc/nsopy/blob/master/README.md

Install the nsopy library using pip.

```APIDOC
## Installation

```bash
pip install nsopy
```
```

--------------------------------

### Instantiate Analytical Example

Source: https://github.com/robin-vjc/nsopy/blob/master/notebooks/02. Application to Duality.ipynb

Creates an instance of the AnalyticalExample class to be used with the dual methods.

```python
inner_problem = AnalyticalExample()
```

--------------------------------

### Install Nsopy using Pip

Source: https://github.com/robin-vjc/nsopy/blob/master/docs/usage.md

Install the nsopy library using pip. Ensure you are in your virtual environment.

```console
(.venv) $ pip install lumache
```

--------------------------------

### CuttingPlanesMethod Setup

Source: https://context7.com/robin-vjc/nsopy/llms.txt

Initializes the classical Cutting Planes Method, which approximates the objective function with a polyhedral structure. This method requires the Gurobi solver and defines an optimality gap tolerance.

```python
import numpy as np
from nsopy.methods.bundle import CuttingPlanesMethod
from nsopy.loggers import GenericDualMethodLogger

def oracle(x_k):
    fi_x_k = [-2*x_k + 2, -1.0/3*x_k + 1, x_k - 2]
    f_x_k = max(fi_x_k)
    diff_fi = [-2, -1.0/3.0, 1]
    max_i = fi_x_k.index(f_x_k)
    diff_f_xk = np.array([diff_fi[max_i]])
    return 0, f_x_k, diff_f_xk

def projection_function(x_k):
    return np.maximum(x_k, 0)

method = CuttingPlanesMethod(
    oracle=oracle,
    projection_function=projection_function,
    dimension=1,
    epsilon=0.01,       # Optimality gap tolerance
    search_box_min=0,   # Lower bound on search region
    search_box_max=10,  # Upper bound on search region
    sense='min'
)

# Optional: configure dual domain constraints
# method.set_dual_domain(type='positive orthant')

logger = GenericDualMethodLogger(method)
```

--------------------------------

### Import Libraries and Set Up Environment

Source: https://github.com/robin-vjc/nsopy/blob/master/notebooks/01. Analytical Example.ipynb

Imports necessary libraries like numpy and matplotlib, and sets the working directory. This is a common setup for numerical and plotting tasks.

```python
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
%matplotlib notebook

%cd ..

```

--------------------------------

### UniversalFGM Implementation

Source: https://context7.com/robin-vjc/nsopy/llms.txt

Demonstrates the setup for the Universal Fast Gradient Method, an accelerated variant designed for optimal convergence rates. It uses an epsilon parameter for line search accuracy and can track the adaptive Lipschitz estimate.

```python
import numpy as np
from nsopy.methods.universal import UniversalFGM
from nsopy.loggers import DualDgmFgmMethodLogger

def oracle(x_k):
    fi_x_k = [-2*x_k + 2, -1.0/3*x_k + 1, x_k - 2]
    f_x_k = max(fi_x_k)
    diff_fi = [-2, -1.0/3.0, 1]
    max_i = fi_x_k.index(f_x_k)
    diff_f_xk = np.array([diff_fi[max_i]])
    return 0, f_x_k, diff_f_xk

def projection_function(x_k):
    return np.maximum(x_k, 0)

method = UniversalFGM(
    oracle=oracle,
    projection_function=projection_function,
    dimension=1,
    epsilon=0.5,
    averaging=False,
    sense='min'
)

logger = DualDgmFgmMethodLogger(method)

for _ in range(30):
    method.step()

print(f"Solution: {logger.lambda_k_iterates[-1]}")
print(f"Adaptive L_k: {logger.L_k_iterates[-1]}")  # Track Lipschitz estimate
```

--------------------------------

### SGMTripleAveraging for Quasi-Monotone Framework

Source: https://context7.com/robin-vjc/nsopy/llms.txt

Implements Subgradient Method with Triple Averaging, offering variants within Nesterov's quasi-monotone framework. This example shows the setup for a problem with a piecewise linear objective function.

```python
import numpy as np
from nsopy.methods.quasi_monotone import SGMTripleAveraging
from nsopy.loggers import GenericDualMethodLogger

def oracle(x_k):
    fi_x_k = [-2*x_k + 2, -1.0/3*x_k + 1, x_k - 2]
    f_x_k = max(fi_x_k)
    diff_fi = [-2, -1.0/3.0, 1]
    max_i = fi_x_k.index(f_x_k)
    diff_f_xk = np.array([diff_fi[max_i]])
    return 0, f_x_k, diff_f_xk

def projection_function(x_k):
    return np.maximum(x_k, 0)
```

--------------------------------

### MILP Lagrangian Dual Decomposition with nsopy

Source: https://context7.com/robin-vjc/nsopy/llms.txt

This example defines a MILP dual problem and compares different nsopy methods (Subgradient, Double Simple Averaging, Universal PGM) for solving it. It requires defining an oracle, projection function, and then iterating through the chosen methods.

```python
import numpy as np
from nsopy.methods.subgradient import SubgradientMethod
from nsopy.methods.quasi_monotone import SGMDoubleSimpleAveraging
from nsopy.methods.universal import UniversalPGM
from nsopy.loggers import EnhancedDualMethodLogger

class MILPDualProblem:
    """
    MILP: min -0.5x1 - x2 + x3
          s.t. 0.5x1 + 0.5x2 + x3 >= 1  (lambda_1)
               x1 + x2 <= 1              (lambda_2)
               0 <= x1,x2,x3 <= 1, x1 binary

    Optimal dual: lambda* = [1, ~1.25], d* = -0.5
    """
    def __init__(self):
        self.dimension = 2
        self.instance_name = 'MILP Example'
        self.instance_type = 'analytical'
        self.instance_subtype = 'lagrangian'

    def oracle(self, lambda_k):
        c = np.array([-0.5 - 0.5*lambda_k[0] + lambda_k[1],
                      -1 - 0.5*lambda_k[0] + lambda_k[1],
                      1 - lambda_k[0]])
        x_k = (c < 0).astype(float)

        diff_d_k = np.array([
            1 - 0.5*x_k[0] - 0.5*x_k[1] - x_k[2],  # constraint 1 slack
            x_k[0] + x_k[1] - 1                      # constraint 2 slack
        ])
        d_k = np.dot(c, x_k) + lambda_k[0] - lambda_k[1]
        return x_k, d_k, diff_d_k

    def projection_function(self, lambda_k):
        return np.maximum(lambda_k, 0)

problem = MILPDualProblem()

# Compare different methods
methods = {
    'Subgradient 1/k': SubgradientMethod(
        problem.oracle, problem.projection_function,
        dimension=2, stepsize_rule='1/k', sense='max'
    ),
    'Double Simple Averaging': SGMDoubleSimpleAveraging(
        problem.oracle, problem.projection_function,
        dimension=2, gamma=1.0, sense='max'
    ),
    'Universal PGM': UniversalPGM(
        problem.oracle, problem.projection_function,
        dimension=2, epsilon=0.1, sense='max'
    )
}

results = {}
for name, method in methods.items():
    logger = EnhancedDualMethodLogger(method)
    for _ in range(100):
        method.step()
    results[name] = {
        'lambda': logger.lambda_k_iterates[-1],
        'd_k': logger.d_k_iterates[-1],
        'oracle_calls': logger.oracle_calls[-1]
    }

for name, res in results.items():
    print(f"{name}:")
    print(f"  lambda* = {res['lambda']}")
    print(f"  d* = {res['d_k']:.4f}")
    print(f"  Oracle calls: {res['oracle_calls']}")

```

--------------------------------

### SGMTripleAveraging Method Variants

Source: https://context7.com/robin-vjc/nsopy/llms.txt

Demonstrates the setup for two variants of the SGMTripleAveraging method, differing in their internal parameter updates (a_t and gamma_t).

```python
method_v1 = SGMTripleAveraging(
    oracle=oracle,
    projection_function=projection_function,
    dimension=1,
    variant=1,  # Options: 1 or 2
    gamma=1.0,
    sense='min'
)

method_v2 = SGMTripleAveraging(
    oracle=oracle,
    projection_function=projection_function,
    dimension=1,
    variant=2,
    gamma=1.0,
    sense='min'
)

logger_v1 = GenericDualMethodLogger(method_v1)
logger_v2 = GenericDualMethodLogger(method_v2)

for _ in range(100):
    method_v1.step()
    method_v2.step()

print(f"Variant 1 result: {logger_v1.lambda_k_iterates[-1]}")
print(f"Variant 2 result: {logger_v2.lambda_k_iterates[-1]}")
```

--------------------------------

### Define Analytical Example Class

Source: https://github.com/robin-vjc/nsopy/blob/master/notebooks/02. Application to Duality.ipynb

Defines an AnalyticalExample class to represent an optimization problem with oracle, projection, and other necessary methods. This class is used to instantiate problems for the dual methods.

```python
class AnalyticalExample(object):
    def __init__(self):
        self.dimension = 2

    def oracle(self, lambda_k):
        assert type(lambda_k) == np.ndarray, 'WARNING: lambda_k should be a numpy array.'
        c = np.zeros(3, dtype=float)
        x_k = np.zeros(3, dtype=float)

        c[0] = - 0.5 - 0.5*lambda_k[0] + lambda_k[1]
        c[1] = - 1 - 0.5*lambda_k[0] + lambda_k[1]
        c[2] = + 1 - lambda_k[0]

        x_k[0] = 1 if c[0] < 0 else 0
        x_k[1] = 1 if c[1] < 0 else 0
        x_k[2] = 1 if c[2] < 0 else 0

        diff_d_k = np.zeros(2)
        diff_d_k[0] = 1 - 0.5*x_k[0] - 0.5*x_k[1] - x_k[2]
        diff_d_k[1] = x_k[0] + x_k[1] -1

        d_k = c[0]*x_k[0] + c[1]*x_k[1] + c[2] *x_k[2] + lambda_k[0] - lambda_k[1]

        return x_k, d_k, diff_d_k

    def projection_function(self, lambda_k):
        # initial point lambda_0
        if isinstance(lambda_k, int) and lambda_k == 0:
            return np.zeros(self.dimension)
        
        # project lambda_k on the positive orthant
        return np.maximum(lambda_k, 0)
```

--------------------------------

### Get Random Ingredients with Nsopy

Source: https://github.com/robin-vjc/nsopy/blob/master/docs/usage.md

Use the `get_random_ingredients` function from the lumache library to retrieve a list of random ingredients. The `kind` parameter is optional and can be 'meat', 'fish', or 'veggies'.

```python
>>> import lumache
>>> lumache.get_random_ingredients()
['shells', 'gorgonzola', 'parsley']
```

--------------------------------

### Initialize OptimalShedding and UniversalPGM Optimizer

Source: https://github.com/robin-vjc/nsopy/blob/master/notebooks/03. Dual Decomposition.ipynb

Sets up the OptimalShedding problem and initializes the UniversalPGM optimizer. The optimizer requires the oracle and projection functions from the problem, along with a specified epsilon value.

```python
problem = OptimalShedding(n_customers=20, n_products=2)
optimizer = UniversalPGM(problem.oracle, problem.projection_function, epsilon=0.01)
# optimizer = UniversalPGM(problem.oracle, problem.projection_function, epsilon=0.1)
# optimizer = SubgradientMethod(problem.oracle, problem.projection_function, stepsize_0=0.1, sense='max')

logger = GenericDualMethodLogger(optimizer)
```

--------------------------------

### Initialize and Run Dual Method

Source: https://github.com/robin-vjc/nsopy/blob/master/notebooks/02. Application to Duality.ipynb

Imports various dual methods from nsopy, instantiates UniversalPGM with the analytical problem, and runs it for a specified number of iterations. The GenericDualMethodLogger is used to track the method's progress.

```python
from nsopy.methods.quasi_monotone import SGMDoubleSimpleAveraging as DSA
from nsopy.methods.quasi_monotone import SGMTripleAveraging as TA
from nsopy.methods.subgradient import SubgradientMethod as SG
from nsopy.methods.universal import UniversalPGM as UPGM
from nsopy.methods.universal import UniversalDGM as UDGM
from nsopy.methods.universal import UniversalFGM as UFGM

from nsopy.loggers import GenericDualMethodLogger

dual_method = UPGM(inner_problem.oracle, inner_problem.projection_function, epsilon=0.01)
method_logger = GenericDualMethodLogger(dual_method)

for iteration in range(60):
    dual_method.dual_step()
```

--------------------------------

### Initialize UniversalPGM optimizer and logger

Source: https://github.com/robin-vjc/nsopy/blob/master/notebooks/03. Dual Decomposition.ipynb

Sets up the UniversalPGM optimizer with the defined problem's oracle and projection functions, and initializes a logger to track optimization iterates.

```python
problem = AnalyticalExample()
optimizer = UniversalPGM(problem.oracle, problem.projection_function, epsilon=0.01)
logger = GenericDualMethodLogger(optimizer)
```

--------------------------------

### Instantiate and Run Subgradient Method

Source: https://github.com/robin-vjc/nsopy/blob/master/README.md

Instantiates the SubgradientMethod with the defined oracle and projection function, sets up a logger, and runs the optimization for 200 iterations.

```python
from nsopy.methods.subgradient import SubgradientMethod
from nsopy.loggers import GenericMethodLogger

method = SubgradientMethod(oracle, projection_function, stepsize_0=0.1, stepsize_rule='constant', sense='min')
logger = GenericMethodLogger(method)

for iteration in range(200):
    method.step()
```

--------------------------------

### Initialize and Run Optimization Method

Source: https://github.com/robin-vjc/nsopy/blob/master/notebooks/01. Analytical Example.ipynb

Initializes a UniversalPGM optimization method and runs it for a specified number of iterations. It also sets a custom initial point.

```python
# method = DSA(oracle, projection_function, dimension=2, gamma=0.5)
# method = TA(oracle, projection_function, dimension=2, variant=2, gamma=0.5)
# method = SG(oracle, projection_function, dimension=2)
method = UPGM(oracle, projection_function, dimension=2, epsilon=10, averaging=True)
# method = UDGM(oracle, projection_function, dimension=2, epsilon=1.0)
# method = UFGM(oracle, projection_function, dimension=2, epsilon=1.0)

method_logger = GenericDualMethodLogger(method)
# start from an different initial point
x_0 = np.array([2.01,2.01])
method.lambda_hat_k = x_0

for iteration in range(100):
    method.dual_step()

```

--------------------------------

### Instantiate Bundle Method

Source: https://github.com/robin-vjc/nsopy/blob/master/README.md

Constructor for the Bundle Method. Requires an oracle and projection function. Requires 'gurobipy'. Optional parameters include dimension, epsilon, mu, and optimization sense.

```python
BundleMethod(oracle, projection_function, dimension=0, epsilon=0.01, mu=0.5, sense='min')
```

--------------------------------

### DualMethodsFactory for Instantiating Optimization Methods

Source: https://context7.com/robin-vjc/nsopy/llms.txt

Utilizes DualMethodsFactory to create optimization method instances with a standardized interface. Demonstrates listing available methods and configuring a specific method ('DSA') with parameters. Requires a problem class adhering to a standard interface.

```python
import numpy as np
from nsopy.methods_factory import DualMethodsFactory, AVAILABLE_METHODS
from nsopy.loggers import EnhancedDualMethodLogger

# Define a problem class with standard interface
class MyProblem:
    def __init__(self):
        self.dimension = 2

    def oracle(self, lambda_k):
        c = np.array([-0.5 - 0.5*lambda_k[0] + lambda_k[1],
                      -1 - 0.5*lambda_k[0] + lambda_k[1],
                      1 - lambda_k[0]])
        x_k = (c < 0).astype(float)
        diff_d_k = np.array([1 - 0.5*x_k[0] - 0.5*x_k[1] - x_k[2],
                             x_k[0] + x_k[1] - 1])
        d_k = np.dot(c, x_k) + lambda_k[0] - lambda_k[1]
        return x_k, d_k, diff_d_k

    def projection_function(self, lambda_k):
        return np.maximum(lambda_k, 0)

problem = MyProblem()

# List available methods
print(f"Available methods: {AVAILABLE_METHODS}")
# ('SG 1/k', 'SG const', 'UPGM', 'UDGM', 'UFGM', 'DSA', 'TA 1', 'TA 2', 'CP', 'bundle')

# Create method using factory
method = DualMethodsFactory(problem, method='DSA', param=1.0)

logger = EnhancedDualMethodLogger(method)

for _ in range(50):
    method.step()

print(f"Result: {logger.lambda_k_iterates[-1]}")
print(f"Value: {logger.d_k_iterates[-1]}")
```

--------------------------------

### Solve Optimization Problem with Subgradient Method

Source: https://github.com/robin-vjc/nsopy/blob/master/docs/index.md

Initialize and run the SubgradientMethod from nsopy. Requires the oracle function, projection function, and optimization parameters like dimension and stepsize rule.

```python
from nsopy.methods.subgradient import SubgradientMethod
from nsopy.loggers import GenericMethodLogger

method = SubgradientMethod(oracle, projection_function, dimension=1, stepsize_0=0.1, stepsize_rule='constant', sense='min')
logger = GenericMethodLogger(method)

for iteration in range(200):
    method.step()
```

--------------------------------

### Instantiate Universal Fast Gradient Method

Source: https://github.com/robin-vjc/nsopy/blob/master/README.md

Constructor for the Universal Fast Gradient Method (UFGM). Requires an oracle and projection function. Optional parameters include dimension, epsilon, averaging flag, and optimization sense.

```python
UniversalFGM(oracle, projection_function, dimension=0, epsilon=1.0, averaging=False, sense='min')
```

--------------------------------

### Import Optimization Methods

Source: https://github.com/robin-vjc/nsopy/blob/master/notebooks/01. Analytical Example.ipynb

Imports various optimization methods from the nsopy library, including quasi-monotone and universal methods.

```python
from nsopy.methods.quasi_monotone import SGMDoubleSimpleAveraging as DSA
from nsopy.methods.quasi_monotone import SGMTripleAveraging as TA
from nsopy.methods.subgradient import SubgradientMethod as SG
from nsopy.methods.universal import UniversalPGM as UPGM
from nsopy.methods.universal import UniversalDGM as UDGM
from nsopy.methods.universal import UniversalFGM as UFGM
from nsopy.loggers import GenericDualMethodLogger

```

--------------------------------

### Instantiate Standard Subgradient Method

Source: https://github.com/robin-vjc/nsopy/blob/master/README.md

Constructor for the Standard Subgradient Method. Requires an oracle and projection function. Optional parameters include dimension, initial stepsize, stepsize rule, and optimization sense.

```python
SubgradientMethod(oracle, projection_function, dimension=0, stepsize_0=1.0, stepsize_rule='1/k', sense='min')
```

--------------------------------

### UniversalPGM Implementation

Source: https://context7.com/robin-vjc/nsopy/llms.txt

Sets up and runs the Universal Primal Gradient Method. Requires defining an oracle function and a projection function. The accuracy is controlled by the epsilon parameter.

```python
import numpy as np
from nsopy.methods.universal import UniversalPGM
from nsopy.loggers import EnhancedDualMethodLogger

def oracle(x_k):
    fi_x_k = [-2*x_k + 2, -1.0/3*x_k + 1, x_k - 2]
    f_x_k = max(fi_x_k)
    diff_fi = [-2, -1.0/3.0, 1]
    max_i = fi_x_k.index(f_x_k)
    diff_f_xk = np.array([diff_fi[max_i]])
    return 0, f_x_k, diff_f_xk

def projection_function(x_k):
    return np.maximum(x_k, 0)

method = UniversalPGM(
    oracle=oracle,
    projection_function=projection_function,
    dimension=1,
    epsilon=0.1,      # Accuracy parameter for line search
    averaging=False,  # Use averaged iterates (theoretical guarantees)
    sense='min'
)

logger = EnhancedDualMethodLogger(method)

for _ in range(50):
    method.step()

print(f"Solution: {logger.lambda_k_iterates[-1]}")
print(f"Value: {logger.d_k_iterates[-1]}")
print(f"Oracle calls: {logger.oracle_calls[-1]}")
print(f"Computation time: {logger.iteration_time[-1]:.4f}s")
```

--------------------------------

### Instantiate Universal Dual Gradient Method

Source: https://github.com/robin-vjc/nsopy/blob/master/README.md

Constructor for the Universal Dual Gradient Method (UDGM). Requires an oracle and projection function. Optional parameters include dimension, epsilon, averaging flag, and optimization sense.

```python
UniversalDGM(oracle, projection_function, dimension=0, epsilon=1.0, averaging=False, sense='min')
```

--------------------------------

### Instantiate Cutting Planes Method

Source: https://github.com/robin-vjc/nsopy/blob/master/README.md

Constructor for the Cutting Planes Method. Requires an oracle and projection function. Requires 'gurobipy'. Optional parameters include dimension, epsilon (stopping criterion), search box min/max, and optimization sense.

```python
CuttingPlanesMethod(oracle, projection_function, dimension=0, epsilon=0.01, search_box_min=-10, search_box_max=10, sense='min')
```

--------------------------------

### Instantiate Universal Primal Gradient Method

Source: https://github.com/robin-vjc/nsopy/blob/master/README.md

Constructor for the Universal Primal Gradient Method (UPGM). Requires an oracle and projection function. Optional parameters include dimension, epsilon, averaging flag, and optimization sense.

```python
UniversalPGM(oracle, projection_function, dimension=0, epsilon=1.0, averaging=False, sense='min')
```

--------------------------------

### Run UniversalPGM optimization steps

Source: https://github.com/robin-vjc/nsopy/blob/master/notebooks/03. Dual Decomposition.ipynb

Executes the dual step of the UniversalPGM optimizer for a fixed number of iterations to find the optimal dual solution.

```python
for iteration in range(60):
    optimizer.dual_step()
```

--------------------------------

### Inspect Solution Iterates

Source: https://github.com/robin-vjc/nsopy/blob/master/docs/index.md

Print the last few iterates of the solution found by the optimization method to inspect convergence.

```python
print(logger.x_k_iterates[-5:])
```

--------------------------------

### Instantiate SGM Double Simple Averaging

Source: https://github.com/robin-vjc/nsopy/blob/master/README.md

Constructor for the SGM Double Simple Averaging method. Requires an oracle and projection function. Optional parameters include dimension, gamma, and optimization sense.

```python
SGMDoubleSimpleAveraging(oracle, projection_function, dimension=0, gamma=1.0, sense='min')
```

--------------------------------

### UniversalDGM Implementation

Source: https://context7.com/robin-vjc/nsopy/llms.txt

Configures and executes the Universal Dual Gradient Method. This method operates in the dual space and can be beneficial for specific problem structures. Averaging can be enabled for theoretical guarantees.

```python
import numpy as np
from nsopy.methods.universal import UniversalDGM
from nsopy.loggers import EnhancedDualMethodLogger

def oracle(x_k):
    fi_x_k = [-2*x_k + 2, -1.0/3*x_k + 1, x_k - 2]
    f_x_k = max(fi_x_k)
    diff_fi = [-2, -1.0/3.0, 1]
    max_i = fi_x_k.index(f_x_k)
    diff_f_xk = np.array([diff_fi[max_i]])
    return 0, f_x_k, diff_f_xk

def projection_function(x_k):
    return np.maximum(x_k, 0)

method = UniversalDGM(
    oracle=oracle,
    projection_function=projection_function,
    dimension=1,
    epsilon=0.5,
    averaging=True,  # Enable averaging for theoretical guarantees
    sense='min'
)

logger = EnhancedDualMethodLogger(method)

for _ in range(30):
    method.step()

print(f"Final solution: {logger.lambda_k_iterates[-1]}")
print(f"Final value: {logger.d_k_iterates[-1]}")
```

--------------------------------

### SubgradientMethod for Non-Smooth Convex Optimization

Source: https://context7.com/robin-vjc/nsopy/llms.txt

Implements the standard subgradient method. Supports 'constant', '1/k', and '1/sqrt(k)' stepsize rules. Requires an oracle and a projection function for the feasible set.

```python
import numpy as np
from nsopy.methods.subgradient import SubgradientMethod
from nsopy.loggers import GenericMethodLogger

# Define the oracle for a piecewise linear function: max(-2x+2, -x/3+1, x-2)
def oracle(x_k):
    fi_x_k = [-2*x_k + 2, -1.0/3*x_k + 1, x_k - 2]
    f_x_k = max(fi_x_k)  # function value
    diff_fi = [-2, -1.0/3.0, 1]  # gradients of components
    max_i = fi_x_k.index(f_x_k)
    diff_f_xk = np.array([diff_fi[max_i]])  # subgradient
    return 0, f_x_k, diff_f_xk

# Projection onto non-negative reals (X = R+)
def projection_function(x_k):
    return np.maximum(x_k, 0)

# Initialize method with constant stepsize
method = SubgradientMethod(
    oracle=oracle,
    projection_function=projection_function,
    dimension=1,
    stepsize_0=0.1,
    stepsize_rule='constant',  # Options: 'constant', '1/k', '1/sqrt(k)'
    sense='min'
)

# Attach logger to track iterates
logger = GenericMethodLogger(method)

# Run optimization
for _ in range(200):
    method.step()

# Access results
print(f"Final iterate: {logger.x_k_iterates[-1]}")  # ~2.25 (optimal x*)
print(f"Final value: {logger.f_k_iterates[-1]}")    # ~-0.25 (optimal f(x*))
```

--------------------------------

### Run Optimization Iterations

Source: https://github.com/robin-vjc/nsopy/blob/master/notebooks/03. Dual Decomposition.ipynb

Executes the optimization process for a fixed number of iterations. The dual step is performed in each iteration, and the output shows the current iteration number.

```python
for iteration in range(10):
    print(f"{iteration=}")
    optimizer.dual_step()

# slight buf in the record of first iteration
logger.x_k_iterates[0] = problem.oracle([0, 0])[0]
```

--------------------------------

### Instantiate SGM Triple Averaging

Source: https://github.com/robin-vjc/nsopy/blob/master/README.md

Constructor for the SGM Triple Averaging method. Requires an oracle and projection function. Optional parameters include dimension, variant, gamma, and optimization sense. Variant can be 1 or 2.

```python
SGMTripleAveraging(oracle, projection_function, dimension=0, variant=1, gamma=1.0, sense='min')
```

--------------------------------

### BundleMethod Optimization

Source: https://context7.com/robin-vjc/nsopy/llms.txt

Implements the BundleMethod for optimization with proximal regularization. Requires a custom oracle and projection function. Suitable for problems where stable convergence is prioritized over speed.

```python
import numpy as np
from nsopy.methods.bundle import BundleMethod
from nsopy.loggers import GenericDualMethodLogger

def oracle(x_k):
    fi_x_k = [-2*x_k + 2, -1.0/3*x_k + 1, x_k - 2]
    f_x_k = max(fi_x_k)
    diff_fi = [-2, -1.0/3.0, 1]
    max_i = fi_x_k.index(f_x_k)
    diff_f_xk = np.array([diff_fi[max_i]])
    return 0, f_x_k, diff_f_xk

def projection_function(x_k):
    return np.maximum(x_k, 0)

method = BundleMethod(
    oracle=oracle,
    projection_function=projection_function,
    dimension=1,
    epsilon=0.01,  # Optimality tolerance
    mu=0.5,        # Proximal weight (regularization)
    sense='min'
)

# Optional: set dual domain constraints
# method.set_dual_domain(type='positive orthant')

logger = GenericDualMethodLogger(method)

for _ in range(20):
    method.step()
    if not method.optimizer_not_yet_found:
        break

print(f"Solution: {method.lambda_k}")
print(f"Optimal value: {method.d_k}")
```

--------------------------------

### GenericMethodLogger for Optimization Progress

Source: https://context7.com/robin-vjc/nsopy/llms.txt

Demonstrates the use of GenericMethodLogger to track optimization iterates (x_k) and function values (f_k). Includes plotting the convergence of function values using matplotlib. Other logger types like EnhancedDualMethodLogger and SlimDualMethodLogger are available for different tracking needs.

```python
import numpy as np
from nsopy.methods.subgradient import SubgradientMethod
from nsopy.loggers import (
    GenericMethodLogger,
    GenericDualMethodLogger,
    EnhancedDualMethodLogger,
    SlimDualMethodLogger
)

def oracle(x_k):
    fi_x_k = [-2*x_k + 2, -1.0/3*x_k + 1, x_k - 2]
    f_x_k = max(fi_x_k)
    diff_fi = [-2, -1.0/3.0, 1]
    max_i = fi_x_k.index(f_x_k)
    diff_f_xk = np.array([diff_fi[max_i]])
    return 0, f_x_k, diff_f_xk

def projection_function(x_k):
    return np.maximum(x_k, 0)

method = SubgradientMethod(oracle, projection_function, dimension=1, sense='min')

# Basic logger - tracks x_k and f_k
logger = GenericMethodLogger(method)

# Enhanced logger - also tracks time and oracle calls
# enhanced_logger = EnhancedDualMethodLogger(method)

# Slim logger - only d_k for memory efficiency on large problems
# slim_logger = SlimDualMethodLogger(method)

for _ in range(100):
    method.step()

# Access logged data
print(f"Iterates: {logger.x_k_iterates[-5:]}")
print(f"Function values: {logger.f_k_iterates[-5:]}")

# Plot convergence
import matplotlib.pyplot as plt
plt.plot(logger.f_k_iterates)
plt.xlabel('Iteration')
plt.ylabel('f(x_k)')
plt.title('Convergence Plot')
plt.show()
```

--------------------------------

### Universal Gradient Methods

Source: https://github.com/robin-vjc/nsopy/blob/master/README.md

Implementation of Nesterov's universal gradient methods (primal, dual, and fast versions).

```APIDOC
* **Universal Gradient Methods**

Implementation of Nesterov's [universal gradient methods](http://link.springer.com/article/10.1007/s10107-014-0790-0), primal, dual and fast versions.

```python
UniversalPGM(oracle, projection_function, dimension=0, epsilon=1.0, averaging=False, sense='min')
UniversalDGM(oracle, projection_function, dimension=0, epsilon=1.0, averaging=False, sense='min'):        
UniversalFGM(oracle, projection_function, dimension=0, epsilon=1.0, averaging=False, sense='min'):
```
```

--------------------------------

### Quasi-Monotone Methods

Source: https://github.com/robin-vjc/nsopy/blob/master/README.md

Implementation of double simple averaging, and triple averaging methods for quasi-monotone optimization.

```APIDOC
* **Quasi-Monotone Methods**

Implementation of double simple averaging, and triple averaging methods from Nesterov's paper on [quasi-monotone methods](http://link.springer.com/article/10.1007/s10957-014-0677-5).

```python
SGMDoubleSimpleAveraging(oracle, projection_function, dimension=0, gamma=1.0, sense='min')
SGMTripleAveraging(oracle, projection_function, dimension=0, variant=1, gamma=1.0, sense='min'):
```

Variants of `SGMTripleAveraging` available: `variant: [1, 2]`
```

--------------------------------

### Inspect Earlier Primal Solutions

Source: https://github.com/robin-vjc/nsopy/blob/master/notebooks/02. Application to Duality.ipynb

View a history of primal solutions (x_k) from earlier iterations stored in the method logger. This helps in understanding solution behavior and potential oscillations.

```python
method_logger.x_k_iterates[-10:]
```

--------------------------------

### Cutting Planes Method

Source: https://github.com/robin-vjc/nsopy/blob/master/README.md

The Cutting Planes Method requires the 'gurobipy' library. It uses an epsilon parameter for suboptimality.

```APIDOC
* **Cutting Planes Method**

*Warning*: this method requires `gurobipy`; if you are an academic, you can get a free license [here](http://www.gurobi.com/academia/for-universities]). 

```python
CuttingPlanesMethod(oracle, projection_function, dimension=0, epsilon=0.01, search_box_min=-10, search_box_max=10, sense='min')
```

The parameter `epsilon` is the absolute required suboptimality level `|f_k - f*|` used as a stopping criterion. Note that a search box needs to be specified.
```

--------------------------------

### SGMDoubleSimpleAveraging for Quasi-Monotone Convergence

Source: https://context7.com/robin-vjc/nsopy/llms.txt

Implements Subgradient Method with Double Simple Averaging for quasi-monotone convergence. Suitable for Lagrangian dual problems. Requires a dual oracle and projection onto the dual feasible set.

```python
import numpy as np
from nsopy.methods.quasi_monotone import SGMDoubleSimpleAveraging
from nsopy.loggers import GenericDualMethodLogger

# Define a 2D analytical problem (Lagrangian dual)
class DualProblem:
    def __init__(self):
        self.dimension = 2

    def oracle(self, lambda_k):
        # Inner problem: min -0.5x1 - x2 + x3 with dualized constraints
        c = np.zeros(3)
        x_k = np.zeros(3)
        c[0] = -0.5 - 0.5*lambda_k[0] + lambda_k[1]
        c[1] = -1 - 0.5*lambda_k[0] + lambda_k[1]
        c[2] = 1 - lambda_k[0]

        x_k[0] = 1 if c[0] < 0 else 0
        x_k[1] = 1 if c[1] < 0 else 0
        x_k[2] = 1 if c[2] < 0 else 0

        diff_d_k = np.array([
            1 - 0.5*x_k[0] - 0.5*x_k[1] - x_k[2],
            x_k[0] + x_k[1] - 1
        ])
        d_k = c[0]*x_k[0] + c[1]*x_k[1] + c[2]*x_k[2] + lambda_k[0] - lambda_k[1]
        return x_k, d_k, diff_d_k

    def projection_function(self, lambda_k):
        return np.maximum(lambda_k, 0)  # Project onto positive orthant

problem = DualProblem()

method = SGMDoubleSimpleAveraging(
    oracle=problem.oracle,
    projection_function=problem.projection_function,
    dimension=problem.dimension,
    gamma=1.0,  # Tuning parameter
    sense='max'  # Maximize dual function
)

logger = GenericDualMethodLogger(method)

for _ in range(50):
    method.step()

print(f"Dual multipliers: {logger.lambda_k_iterates[-1]}")
print(f"Dual value: {logger.d_k_iterates[-1]}")  # Should approach -0.5
```

--------------------------------

### Standard Subgradient Method

Source: https://github.com/robin-vjc/nsopy/blob/master/README.md

Details on the Standard Subgradient Method.

```APIDOC
## Available Methods

* **Standard Subgradient Method**

```python
SubgradientMethod(oracle, projection_function, dimension=0, stepsize_0=1.0, stepsize_rule='1/k', sense='min')
```
Stepsize rules valiable: `stepsize_rule: ['constant', '1/k', '1/sqrt(k)']`
```

--------------------------------

### Bundle Method

Source: https://github.com/robin-vjc/nsopy/blob/master/README.md

A basic variant of the Bundle Method, also requiring 'gurobipy'.

```APIDOC
* **Bundle Method**

*Warning*: this method requires `gurobipy`; if you are an academic, you can get a free license [here](http://www.gurobi.com/academia/for-universities]). 

Implementation of a basic variant of the bundle method. 

```python
BundleMethod(oracle, projection_function, dimension=0, epsilon=0.01, mu=0.5, sense='min'):
```
```

--------------------------------

### Inspect last dual iterates

Source: https://github.com/robin-vjc/nsopy/blob/master/notebooks/03. Dual Decomposition.ipynb

Displays the last five dual iterates found by the optimizer, representing approximate solutions to the dual problem.

```python
logger.lambda_k_iterates[-5:]
```

--------------------------------

### Run Optimization Method Until Convergence

Source: https://context7.com/robin-vjc/nsopy/llms.txt

Iteratively calls the step method of an optimization object until convergence criteria are met or a maximum number of iterations is reached. Prints the final solution and optimal value.

```python
for _ in range(20):
    method.step()
    if not method.optimizer_not_yet_found:
        break

print(f"Solution: {method.lambda_k}")
print(f"Optimal value: {method.d_k}")
```

--------------------------------

### Important Remarks

Source: https://github.com/robin-vjc/nsopy/blob/master/README.md

Key considerations for using the nsopy library, including dimension handling, projection functions, and performance.

```APIDOC
## Important Remarks

* Methods have to either be instantiated with the appropriate dimension argument, or implement a special case for 0. 
The basic usage example above illustrates an oracle implementing such a special case. 
For this example, alternatively one could have instantiated the solution method with `dimension = 1`.  

* The first-order oracle must also provide a projection function; [here is a list of cases](docs/img/simple_projections.png) for which 
the projection operation is computationally inexpensive. 

* Currently, all methods are implemented in Python. Numerical performance is not optimized, but they may
be still useful for quick comparisons or for applications in which the main computational burden is in
evaluating the first order oracle. 
```

--------------------------------

### OptimalShedding class for resource allocation

Source: https://github.com/robin-vjc/nsopy/blob/master/notebooks/03. Dual Decomposition.ipynb

Defines a class for optimal shedding problems, modeling customer demand and product availability. It uses OR-Tools for solving subproblems.

```python
class OptimalShedding:
    def __init__(self, n_customers: int = 10, n_products: int = 2):
        self.n_customers = n_customers
        self.n_products = n_products

        low_demand = 2
        high_demand = 10
        mean_demand = low_demand + (high_demand - low_demand)/2

        low_reward = 8
        high_reward = 12
        
        # model data
        self.R_i = np.random.randint(low_reward, high_reward, n_customers)
        self.P_ij = np.random.randint(low_reward*n_customers, high_reward*n_customers, (n_customers, n_products))

        self.M_ij = np.random.randint(low_demand, low_demand+2, (n_customers, n_products))
        self.D_ij = np.random.randint(low_demand, high_demand, (n_customers, n_products))
        self.C_j = mean_demand*0.5*np.random.randint(int(n_customers*1/3), int(n_customers*2/3), n_products)
    
    def solve_subproblem(self, i: int, lambda_k: np.ndarray):
        solver = pywraplp.Solver.CreateSolver("SAT")

        # add variables
        x_i = None
        y_ij = [None for _ in range(self.n_products)]

        x_i = solver.BoolVar(f"x_{i}")
        for j in range(self.n_products):
            y_ij[j] = solver.NumVar(0.0, float(self.D_ij[i, j]), f"D_({i},{j})")

        # add constraints
        for j in range(self.n_products):
            solver.Add(self.M_ij[i, j]*x_i <= y_ij[j], f"Min_qty_({i}, {j})")
            solver.Add(y_ij[j] <= self.D_ij[i, j]*x_i, f"Demand_({i}, {j})")

        # objective
        solver.Minimize(-self.R_i[i]*x_i + sum([(lambda_k[j] - (self.P_ij[i,j]/(self.D_ij[i,j])))*y_ij[j] for j in range(self.n_products)]))

        # solve
        status = solver.Solve()

        if status != pywraplp.Solver.OPTIMAL:
            raise ValueError(f"Could not solve model: {status}")

        # return values
        x_ik = np.array([x_i.solution_value()] + [y.solution_value() for y in y_ij])
        d_ik = solver.Objective().Value()
```

--------------------------------

### Inspect last primal iterates

Source: https://github.com/robin-vjc/nsopy/blob/master/notebooks/03. Dual Decomposition.ipynb

Displays the last five primal iterates, which are approximate solutions to the original primal optimization problem.

```python
logger.x_k_iterates[-5:]
```

--------------------------------

### Import necessary libraries for NSopy

Source: https://github.com/robin-vjc/nsopy/blob/master/notebooks/03. Dual Decomposition.ipynb

Imports essential libraries for numerical operations, plotting, and NSopy's optimization modules. Sets a random seed for reproducibility and enables interactive matplotlib plots.

```python
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
from ortools.linear_solver import pywraplp

from nsopy.methods.universal import UniversalPGM, UniversalDGM, UniversalFGM
from nsopy.methods.subgradient import SubgradientMethod
from nsopy.loggers import GenericDualMethodLogger

np.random.seed(102)
%matplotlib notebook
```

--------------------------------

### Visualize dual space trajectory

Source: https://github.com/robin-vjc/nsopy/blob/master/notebooks/03. Dual Decomposition.ipynb

Plots the trajectory of dual iterates in 3D space, showing the dual objective function surface and the path taken by the optimizer.

```python
box = np.linspace(0, 3, 30)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

ax.plot_trisurf(np.array([lmd_1 for lmd_1 in box for lmd_2 in box]),
                np.array([lmd_2 for lmd_1 in box for lmd_2 in box]),
                np.array([problem.oracle(np.array([lmd_1, lmd_2]))[1] for lmd_1 in box for lmd_2 in box]),
                alpha=0.7)

plt.plot([lmd[0] for lmd in logger.lambda_k_iterates],
         [lmd[1] for lmd in logger.lambda_k_iterates],
         [d_lmd for d_lmd in logger.d_k_iterates], 'r.-')

ax.set_xlabel('$λ_1$')
ax.set_ylabel('$λ_2$')
ax.set_zlabel('$d(λ)$')
```

--------------------------------

### Analyze Customer Service and Shedding

Source: https://github.com/robin-vjc/nsopy/blob/master/notebooks/03. Dual Decomposition.ipynb

Analyzes the last primal iterate to determine which customers are served and how. It identifies customers that are not served (shed) and prints details about the service units and demands for served customers.

```python
x_sol = logger.x_k_iterates[-1]
not_served = []
for i in range(problem.n_customers):
    if x_sol[i,0] == 0:
        not_served.append(i)

print(f'We have {problem.n_customers} customers.')
print(f'We shed customers {not_served}.')

for i in range(problem.n_customers):
    if i not in not_served:
        print(f'Customer {i} is served {x_sol[i,1]} "flop" units (min: {problem.M_ij[i,0]}, demand: {problem.D_ij[i,0]}) and {x_sol[i,2]} "network" units (min: {problem.M_ij[i,1]}, demand: {problem.D_ij[i,1]}) ')
```