import numpy as np
import scipy.optimize as opt
from scipy.optimize import differential_evolution, basinhopping
import itertools
from typing import List, Tuple, Dict, Any
import time
class FunctionOptimizer:
"""
Optimizer for the function:
f(i1,i2,...,ik) = (sum_j p_j * (1-a_j^i_j)/(1-a_j)) * prod_j q_j^i_j
Where 0 < p_j, q_j, a_j < 1 and i_j are positive integers
"""
def __init__(self, p: List[float], a: List[float], q: List[float]):
"""
Initialize with parameters
Args:
p: List of p_j values
a: List of a_j values
q: List of q_j values
"""
self.p = np.array(p)
self.a = np.array(a)
self.q = np.array(q)
self.k = len(p)
def evaluate(self, indices: List[int]) -> float:
"""
Evaluate the objective function at given indices
Args:
indices: List of positive integers [i_1, i_2, ..., i_k]
Returns:
Function value
"""
indices = np.array(indices, dtype=int)
# Calculate sum term: sum_j p_j * (1-a_j^i_j)/(1-a_j)
geometric_terms = (1 - self.a**indices) / (1 - self.a)
sum_term = np.sum(self.p * geometric_terms)
# Calculate product term: prod_j q_j^i_j
product_term = np.prod(self.q**indices)
return sum_term * product_term
def coordinate_descent(self, start_point: List[int] = None,
max_iter: int = 50, max_value: int = 25) -> Dict[str, Any]:
"""
Coordinate descent optimization
Args:
start_point: Starting point (defaults to all 1s)
max_iter: Maximum iterations
max_value: Maximum value to try for each coordinate
Returns:
Dictionary with solution, value, and history
"""
if start_point is None:
start_point = [1] * self.k
current = np.array(start_point, dtype=int)
current_value = self.evaluate(current)
history = []
for iteration in range(max_iter):
improved = False
for j in range(self.k):
best_local = current[j]
best_local_value = current_value
# Try different values for coordinate j
for test_val in range(1, max_value + 1):
test_solution = current.copy()
test_solution[j] = test_val
test_value = self.evaluate(test_solution)
if test_value > best_local_value:
best_local = test_val
best_local_value = test_value
if best_local_value > current_value:
current[j] = best_local
current_value = best_local_value
improved = True
history.append({
'iteration': iteration + 1,
'solution': current.copy(),
'value': current_value
})
if not improved:
break
return {
'solution': current.tolist(),
'value': current_value,
'history': history,
'converged': not improved
}
def grid_search_around(self, center: List[int], radius: int = 3) -> Dict[str, Any]:
"""
Grid search around a center point
Args:
center: Center point for search
radius: Search radius in each dimension
Returns:
Dictionary with best solution and all results
"""
center = np.array(center)
best_solution = center.copy()
best_value = self.evaluate(center)
all_results = []
# Generate all offset combinations
offsets = list(itertools.product(range(-radius, radius + 1), repeat=self.k))
for offset in offsets:
candidate = center + np.array(offset)
# Ensure all coordinates are positive
candidate = np.maximum(candidate, 1)
value = self.evaluate(candidate)
all_results.append({
'solution': candidate.tolist(),
'value': value,
'offset': list(offset)
})
if value > best_value:
best_solution = candidate.copy()
best_value = value
# Sort results by value (descending)
all_results.sort(key=lambda x: x['value'], reverse=True)
return {
'best_solution': best_solution.tolist(),
'best_value': best_value,
'all_results': all_results
}
def scipy_optimize(self, method='differential_evolution', **kwargs) -> Dict[str, Any]:
"""
Use scipy optimization methods with continuous relaxation
Args:
method: Optimization method ('differential_evolution', 'basinhopping', 'minimize')
**kwargs: Additional arguments for the optimizer
Returns:
Dictionary with results
"""
# Define bounds (continuous relaxation)
bounds = [(1, 25)] * self.k # Reasonable bounds for each variable
# Objective function for minimization (negative of our function)
def objective(x):
# Round to nearest integers
indices = np.round(x).astype(int)
return -self.evaluate(indices)
start_time = time.time()
if method == 'differential_evolution':
# Global optimization
result = differential_evolution(
objective,
bounds,
seed=42,
maxiter=kwargs.get('maxiter', 100),
popsize=kwargs.get('popsize', 15)
)
elif method == 'basinhopping':
# Basin hopping
x0 = np.array([5] * self.k) # Starting point
result = basinhopping(
objective,
x0,
niter=kwargs.get('niter', 100),
T=kwargs.get('T', 1.0),
stepsize=kwargs.get('stepsize', 2.0)
)
else:
# Local optimization from multiple starting points
best_result = None
best_value = float('inf')
start_points = [[1]*self.k, [3]*self.k, [5]*self.k, [10]*self.k]
for start in start_points:
try:
res = opt.minimize(
objective,
start,
bounds=bounds,
method=kwargs.get('optimizer', 'L-BFGS-B')
)
if res.fun < best_value:
best_result = res
best_value = res.fun
except:
continue
result = best_result
if result is not None:
# Convert back to integers and get actual function value
int_solution = np.round(result.x).astype(int)
actual_value = self.evaluate(int_solution)
return {
'solution': int_solution.tolist(),
'value': actual_value,
'continuous_solution': result.x.tolist(),
'scipy_result': result,
'time': time.time() - start_time,
'method': method
}
else:
return {'error': 'Optimization failed'}
def comprehensive_search(self, max_coord_iter: int = 50,
grid_radius: int = 3) -> Dict[str, Any]:
"""
Run comprehensive optimization using multiple methods
Args:
max_coord_iter: Maximum iterations for coordinate descent
grid_radius: Radius for grid search
Returns:
Dictionary with all results
"""
results = {}
start_time = time.time()
print("Running comprehensive optimization...")
# 1. Coordinate descent from multiple starting points
print("1. Coordinate descent...")
start_points = [[1]*self.k, [2]*self.k, [5]*self.k, [10]*self.k]
coord_results = []
for i, start in enumerate(start_points):
print(f" Starting point {i+1}: {start}")
result = self.coordinate_descent(start, max_coord_iter)
coord_results.append({
'start_point': start,
'result': result
})
results['coordinate_descent'] = coord_results
# Find best coordinate descent result
best_coord = max(coord_results, key=lambda x: x['result']['value'])
# 2. Grid search around best coordinate result
print("2. Grid search refinement...")
grid_result = self.grid_search_around(best_coord['result']['solution'], grid_radius)
results['grid_search'] = grid_result
# 3. Scipy optimization methods
print("3. Scipy optimization methods...")
scipy_methods = ['differential_evolution', 'basinhopping']
scipy_results = {}
for method in scipy_methods:
print(f" {method}")
try:
scipy_results[method] = self.scipy_optimize(method)
except Exception as e:
scipy_results[method] = {'error': str(e)}
results['scipy_methods'] = scipy_results
# 4. Find overall best
all_solutions = []
# Add coordinate descent solutions
for coord_res in coord_results:
all_solutions.append({
'method': 'coordinate_descent',
'solution': coord_res['result']['solution'],
'value': coord_res['result']['value']
})
# Add grid search best
all_solutions.append({
'method': 'grid_search',
'solution': grid_result['best_solution'],
'value': grid_result['best_value']
})
# Add scipy results
for method, res in scipy_results.items():
if 'solution' in res:
all_solutions.append({
'method': f'scipy_{method}',
'solution': res['solution'],
'value': res['value']
})
# Find best overall
if all_solutions:
best_overall = max(all_solutions, key=lambda x: x['value'])
results['best_overall'] = best_overall
results['total_time'] = time.time() - start_time
print(f"Optimization complete! Total time: {results['total_time']:.2f}s")
return results
def main():
"""Test the optimizer with the known test case"""
# Test parameters
p = [100, 50, 250,150,100,50]
a = [0.8, 0.75, 0.825,0.9,0.8,99]
q = [0.95, 0.995, 0.9,0.925,0.95,0.975]
known_solution = [1,7,2,2,1,13]
expected_value = 750.440997
print("=" * 60)
print("FUNCTION OPTIMIZER TEST")
print("=" * 60)
print(f"Parameters:")
print(f" p = {p}")
print(f" a = {a}")
print(f" q = {q}")
print(f"Known solution: {known_solution}")
print(f"Expected value: {expected_value}")
print()
# Create optimizer
optimizer = FunctionOptimizer(p, a, q)
# Verify known solution
calculated_value = optimizer.evaluate(known_solution)
print(f"Verification:")
print(f" Calculated value: {calculated_value:.6f}")
print(f" Expected value: {expected_value:.6f}")
print(f" Difference: {abs(calculated_value - expected_value):.6f}")
print()
# Run comprehensive optimization
results = optimizer.comprehensive_search()
# Display results
print("=" * 60)
print("RESULTS SUMMARY")
print("=" * 60)
if 'best_overall' in results:
best = results['best_overall']
print(f"Best solution found: {best['solution']}")
print(f"Best value: {best['value']:.6f}")
print(f"Method: {best['method']}")
if best['value'] > expected_value:
improvement = best['value'] - expected_value
print(f"✨ IMPROVEMENT: +{improvement:.6f}")
else:
print(f"Difference from known: {best['value'] - expected_value:.6f}")
print(f"\nTotal computation time: {results['total_time']:.2f}s")
# Detailed results
print("\n" + "=" * 60)
print("DETAILED RESULTS")
print("=" * 60)
print("\nCoordinate Descent Results:")
for i, coord_res in enumerate(results['coordinate_descent']):
res = coord_res['result']
print(f" Start {coord_res['start_point']}: "
f"{res['solution']} → {res['value']:.6f} "
f"({len(res['history'])} iterations)")
print(f"\nGrid Search Best: {results['grid_search']['best_solution']} "
f"→ {results['grid_search']['best_value']:.6f}")
print(f"\nTop 5 grid search candidates:")
for i, candidate in enumerate(results['grid_search']['all_results'][:5]):
print(f" {i+1}. {candidate['solution']} → {candidate['value']:.6f}")
print("\nScipy Methods:")
for method, res in results['scipy_methods'].items():
if 'solution' in res:
print(f" {method}: {res['solution']} → {res['value']:.6f} "
f"({res['time']:.2f}s)")
else:
print(f" {method}: {res.get('error', 'Failed')}")
main()
To embed this program on your website, copy the following code and paste it into your website's HTML: