Particle Swarm Optimization: Collective Intelligence for Global Search

import Foundation
import BusinessMath

// MARK: - Portfolio with Sector Constraints
// Portfolio with sector constraints (creates local minima)
let numAssets = 10
let sectors = [0, 0, 0, 1, 1, 1, 1, 2, 2, 2]  // 3 Tech, 4 Finance, 3 Healthcare
let sectorLimits = [0.40, 0.40, 0.30]  // Max 40% Tech, 40% Finance, 30% Healthcare

// Generate covariance matrix with sector correlation structure
let covariance = generateCovarianceMatrix(size: numAssets, avgCorrelation: 0.25)

// Portfolio objective with constraints
func portfolioObjective(_ rawWeights: VectorN) -> Double {
	// Normalize weights to sum to 1.0 (simplex projection)
	let sum = rawWeights.toArray().reduce(0, +)
	guard sum > 0 else { return 1e10 }  // Avoid division by zero

	let weights = VectorN(rawWeights.toArray().map { $0 / sum })

	// Calculate portfolio variance
	var variance = 0.0
	for i in 0.. sectorLimits[sectorID] {
			sectorPenalty += pow(sectorWeight - sectorLimits[sectorID], 2) * 100.0
		}
	}

	return variance + sectorPenalty
}

// Particle Swarm Optimizer
// Search space: [0, 1] for each asset (will be normalized to sum to 1)
let pso = ParticleSwarmOptimization>(
	config: ParticleSwarmConfig(
		swarmSize: 50,
		maxIterations: 30,
		inertiaWeight: 0.7,
		cognitiveCoefficient: 1.5,
		socialCoefficient: 1.5
	),
	searchSpace: Array(repeating: (0.0, 1.0), count: numAssets)  // 10 dimensions
)

print("Sector-Constrained Portfolio Optimization")
print("═══════════════════════════════════════════════════════════")

// Start from equal weights
let initialGuess = VectorN(Array(repeating: 1.0 / Double(numAssets), count: numAssets))

let result = try pso.minimize(
	portfolioObjective,
	from: initialGuess
)

// Normalize final solution
let finalSum = result.solution.toArray().reduce(0, +)
let finalWeights = VectorN(result.solution.toArray().map { $0 / finalSum })

print("Optimization Results:")
print("  Best Variance: \(result.value.number(6))")
print("  Iterations: \(result.iterations)")
print("  Swarm Size: 50")
print("  Total Evaluations: \(result.iterations * 50)")

// Verify constraints
let totalWeight = finalWeights.toArray().reduce(0, +)
print("\nPortfolio Weights (total: \(totalWeight.percent(1))):")
for (i, weight) in finalWeights.toArray().enumerated() {
	print("  Asset \(i) (Sector \(sectors[i])): \(weight.percent(1))")
}

print("\nSector Allocations:")
for sectorID in 0..<3 {
	let sectorWeight = finalWeights.toArray().enumerated()
		.filter { sectors[$0.offset] == sectorID }
		.map { $1 }
		.reduce(0, +)

	let limit = sectorLimits[sectorID]
	let status = sectorWeight <= limit + 0.01 ? "✓" : "✗"  // Small tolerance

	print("  Sector \(sectorID): \(sectorWeight.percent(1)) (limit: \(limit.percent(0))) \(status)")
}

// MARK: - Hyperparameter Tuning
print("\n\nPattern 2: PSO Configuration Comparison")
print("═══════════════════════════════════════════════════════════")

// Rastrigin function (highly multimodal)
func rastrigin(_ x: VectorN) -> Double {
	let A = 10.0
	let n = Double(5)  // 5 dimensions
	return A * n + (0..<5).reduce(0.0) { sum, i in
		sum + (x[i] * x[i] - A * cos(2 * .pi * x[i]))
	}
}

let searchSpace2 = (0..<5).map { _ in (-5.12, 5.12) }

let configs2: [(name: String, config: ParticleSwarmConfig)] = [
	("Small Swarm", ParticleSwarmConfig(
		swarmSize: 20,
		maxIterations: 100,
		seed: 101
	)),
	("Default", .default),
	("Large Swarm", ParticleSwarmConfig(
		swarmSize: 100,
		maxIterations: 200,
		seed: 101
	))
]

print("\nComparing configurations on 5D Rastrigin function")
print("Known minimum: [0,0,0,0,0] with value 0.0")
print("\nConfig          Swarm  Iters  Final Value  Converged")
print("────────────────────────────────────────────────────────")

for (name, config) in configs2 {
	let pso = ParticleSwarmOptimization>(
		config: config,
		searchSpace: searchSpace2
	)

	let result = pso.optimizeDetailed(
		objective: rastrigin
	)

	print("\(name.padding(toLength: 14, withPad: " ", startingAt: 0))  " +
		  "\(Double(config.swarmSize).number(0).paddingLeft(toLength: 5))   " +
		  "\(Double(config.maxIterations).number(0).paddingLeft(toLength: 4))   " +
		  "\(result.fitness.number(6).paddingLeft(toLength: 10))   " +
		  "\(result.converged ? "✓" : "✗")")
}

print("\n💡 Observation: Larger swarms find better solutions but take longer")

//// Optimize model hyperparameters: [learningRate, regularization, hiddenLayers, batchSize]
//func modelPerformance(_ hyperparameters: VectorN) -> Double {
//	let learningRate = hyperparameters[0]
//	let regularization = hyperparameters[1]
//	let hiddenLayers = Int(hyperparameters[2].rounded())  // Discrete!
//	let batchSize = Int(hyperparameters[3].rounded())     // Discrete!
//
//	// Train model with these hyperparameters (expensive!)
//	let model = trainModel(
//		lr: learningRate,
//		reg: regularization,
//		layers: hiddenLayers,
//		batch: batchSize
//	)
//
//	// Return validation error (minimize)
//	return model.validationError
//}
//
//let hyperparamPSO = ParticleSwarmOptimization>(
//	config: ParticleSwarmConfig(
//	swarmSize: 50,
//	inertiaWeight: 0.8,
//	cognitiveCoefficient: 2.0,
//	socialCoefficient: 2.0
//	),
//	searchSpace: [(-10.0, -10.0), (10.0, 10.0)]
//)
//
//let hyperparamResult = try hyperparamPSO.minimize(
//	modelPerformance,
//	bounds: [
//		(0.0001, 0.1),    // Learning rate
//		(0.0, 0.01),      // Regularization
//		(1.0, 10.0),      // Hidden layers (will round)
//		(16.0, 256.0)     // Batch size (will round)
//	],
//	maxIterations: 50
//)
//
//print("\nHyperparameter Optimization:")
//print("  Learning Rate: \(hyperparamResult.position[0].number(6))")
//print("  Regularization: \(hyperparamResult.position[1].number(6))")
//print("  Hidden Layers: \(Int(hyperparamResult.position[2].rounded()))")
//print("  Batch Size: \(Int(hyperparamResult.position[3].rounded()))")
//print("  Validation Error: \(hyperparamResult.value.number(4))")


//// MARK: - Pattern 1: Multi-Modal Portfolio with Sector Constraints
//
//print("Pattern 1: Portfolio Optimization with Sector Constraints")
//print("═══════════════════════════════════════════════════════════")
//
//let numAssets1 = 30
//
//// Create tiered expected returns
//let returns1 = (0.. Double in
//    if i < 10 { return Double.random(in: 0.06...0.09) }  // Low
//    else if i < 20 { return Double.random(in: 0.09...0.12) }  // Medium
//    else { return Double.random(in: 0.12...0.16) }  // High
//}
//
//// Sector assignments (3 sectors: Tech, Finance, Energy)
//let sectors1 = (0.. Int in
//    i % 3  // Distribute across 3 sectors
//}
//
//// Volatilities by sector
//let sectorVolatility: [Double] = [0.25, 0.18, 0.22]  // Tech, Finance, Energy
//
//func portfolioObjective1(_ weights: VectorN) -> Double {
//    // Calculate portfolio return (negative because we minimize)
//    let portfolioReturn = zip(weights.toArray(), returns1).reduce(0.0) { $0 + $1.0 * $1.1 }
//
//    // Calculate portfolio variance
//    var variance = 0.0
//    for i in 0..>(
//    config: ParticleSwarmConfig(
//        swarmSize: 100,
//        maxIterations: 200,
//        inertiaWeight: 0.7,
//        cognitiveCoefficient: 1.5,
//        socialCoefficient: 1.5,
//        seed: 42
//    ),
//    searchSpace: searchSpace1
//)
//
//let initialWeights1 = VectorN(Array(repeating: 1.0 / Double(numAssets1), count: numAssets1))
//
//// Constraints
//let constraints1: [MultivariateConstraint>] = [
//    // Budget: sum to 1
//    .equality { weights in
//        weights.toArray().reduce(0.0, +) - 1.0
//    },
//
//    // Sector limits: no sector > 40%
//    .inequality { weights in
//        let techWeight = (0.. $1.weight }.prefix(5)
//
//for holding in holdings1 {
//    print("  Asset \(holding.index) (\(holding.sector)): \(holding.weight.percent()) @ \(holding.return.percent())")
//}
//
// MARK: - Pattern 2: Hyperparameter Search



// MARK: - Pattern 3: Hybrid PSO + Local Refinement

print("\n\nPattern 3: Hybrid PSO + L-BFGS Refinement")
print("═══════════════════════════════════════════════════════════")

// Rosenbrock function (smooth but narrow valley)
func rosenbrock2D(_ x: VectorN) -> Double {
    let a = x[0], b = x[1]
    return (1.0 - a) * (1.0 - a) + 100.0 * (b - a * a) * (b - a * a)
}

let searchSpace3 = [(-5.0, 5.0), (-5.0, 5.0)]

print("\nPhase 1: Global search with PSO")
let pso3 = ParticleSwarmOptimization>(
    config: ParticleSwarmConfig(
        swarmSize: 50,
        maxIterations: 100,
        seed: 42
    ),
    searchSpace: searchSpace3
)

let psoResult = pso3.optimizeDetailed(objective: rosenbrock2D)

print("  PSO Solution: [\(psoResult.solution[0].number(4)), \(psoResult.solution[1].number(4))]")
print("  PSO Value: \(psoResult.fitness.number(6))")
print("  Iterations: \(psoResult.iterations)")

print("\nPhase 2: Local refinement with L-BFGS")
let lbfgs = MultivariateLBFGS>()
let refinedResult = try lbfgs.minimizeLBFGS(
    function: rosenbrock2D,
    initialGuess: psoResult.solution
)

print("  Refined Solution: [\(refinedResult.solution[0].number(6)), \(refinedResult.solution[1].number(6))]")
print("  Refined Value: \(refinedResult.value.number(10))")
print("  L-BFGS Iterations: \(refinedResult.iterations)")

let improvement = ((psoResult.fitness - refinedResult.value) / psoResult.fitness)
print("\n  Improvement from refinement: \(improvement.percent(2))")
//
// MARK: - Pattern 4: Real-World Wind Farm Layout

print("\n\nPattern 4: Wind Farm Turbine Layout Optimization")
print("═══════════════════════════════════════════════════════════")

let numTurbines = 10
let farmWidth = 2000.0  // meters
let farmHeight = 1500.0  // meters
let minSpacing = 200.0  // meters (wake effect)

func windFarmPower(_ positions: VectorN) -> Double {
    // positions: [x1, y1, x2, y2, ..., x10, y10]
    // Simplified model: maximize total power considering wake effects

    var totalPower = 0.0

    for i in 0.. x_j {  // Assuming wind from west (left)
                let lateralDistance = abs(y_i - y_j)
                if lateralDistance < 300 {  // In wake zone
                    let wakeEffect = max(0, 1.0 - distance / 1000.0)
                    turbinePower *= (1.0 - 0.3 * wakeEffect)
                }
            }

            // Penalty for being too close
            if distance < minSpacing {
                turbinePower *= 0.5
            }
        }

        totalPower += turbinePower
    }

    return -totalPower  // Negative because minimizing
}

// Search space: x ∈ [0, farmWidth], y ∈ [0, farmHeight]
let searchSpace4 = (0..>(
    config: ParticleSwarmConfig(
        swarmSize: 80,
        maxIterations: 150,
        seed: 42
    ),
    searchSpace: searchSpace4
)

print("\nOptimizing layout for \(numTurbines) turbines")
print("Farm dimensions: \(farmWidth.number(0))m × \(farmHeight.number(0))m")
print("Minimum spacing: \(minSpacing.number(0))m")

let result4 = pso4.optimizeDetailed(
    objective: windFarmPower
)

print("\nResults:")
print("  Total Power: \((-result4.fitness).number(4)) MW (normalized)")
print("  Iterations: \(result4.iterations)")
print("  Converged: \(result4.converged)")

print("\nTurbine Positions:")
for i in 0..