pso.go

package eaopt

import (
	"errors"
	"math"
	"math/rand"
	"sync"
)

// A Particle is an element of a Swarm. It tracks it's current position,
// the best position it has encountered, and a velocity vector. It also has a
// pointer to the SPSO which generated it so that it can access the function to
// minimize and the global best position.
type Particle struct {
	CurrentX []float64
	CurrentY float64
	BestX    []float64
	BestY    float64
	Velocity []float64
	SPSO     *SPSO
}

// Evaluate the Particle by computing the value of the function at the current
// position. If the position is better than the best position encountered by
// the Particle then it replaces it. Likewhise, the global best position is
// replaced if the current position is better.
func (p *Particle) Evaluate() (float64, error) {
	p.CurrentY = p.SPSO.F(p.CurrentX)
	// Update the Particle's best position
	if p.CurrentY < p.BestY {
		p.BestX = copyFloat64s(p.CurrentX)
		p.BestY = p.CurrentY
	}
	// Update the global best position. In case of parallelism a lock has to be
	// used to handle concurrent access.
	if !p.SPSO.GA.ParallelEval {
		if p.CurrentY < p.SPSO.BestY {
			p.SPSO.BestX = copyFloat64s(p.CurrentX)
			p.SPSO.BestY = p.CurrentY
		}
	} else {
		p.SPSO.mutex.Lock()
		if p.CurrentY < p.SPSO.BestY {
			p.SPSO.BestX = copyFloat64s(p.CurrentX)
			p.SPSO.BestY = p.CurrentY
		}
		p.SPSO.mutex.Unlock()
	}
	return p.CurrentY, nil
}

// Mutate the Particle by modifying it's velocity and it's current position.
func (p *Particle) Mutate(rng *rand.Rand) {
	var (
		rX = make([]float64, len(p.CurrentX))
		ss float64
	)
	for i, xi := range p.CurrentX {
		G := xi + 1.193*(p.BestX[i]+p.SPSO.BestX[i]-2*xi)
		min, max := xi-G, G-xi
		rX[i] = min + rng.Float64()*(max-min)
		ss += rX[i] * rX[i]
	}
	ss = math.Sqrt(ss)
	for i, xi := range p.CurrentX {
		p.Velocity[i] = p.SPSO.W*p.Velocity[i] + rX[i]/ss - xi
		p.CurrentX[i] += p.Velocity[i]
	}
}

// Crossover doesn't do anything.
func (p *Particle) Crossover(q Genome, rng *rand.Rand) {}

// Clone returns a deep copy of the Particle.
func (p Particle) Clone() Genome {
	return &Particle{
		CurrentX: copyFloat64s(p.CurrentX),
		CurrentY: p.CurrentY,
		BestX:    copyFloat64s(p.BestX),
		BestY:    p.BestY,
		Velocity: copyFloat64s(p.Velocity),
		SPSO:     p.SPSO,
	}
}

// SPSO implements the 2011 version of Standard Particle Swarm Optimization. It
// can optimize single-output real-valued functions.
// Reference: http://clerc.maurice.free.fr/pso/SPSO_descriptions.pdf
type SPSO struct {
	Min, Max float64 // Boundaries for initial values
	W        float64
	NDims    uint
	BestX    []float64
	BestY    float64
	F        func([]float64) float64
	GA       *GA
	mutex    *sync.Mutex
}

// NewSPSO instantiates and returns a SPSO instance after having checked for
// input errors.
func NewSPSO(nParticles, nSteps uint, min, max, w float64, parallel bool, rng *rand.Rand) (*SPSO, error) {
	// Check inputs
	if min >= max {
		return nil, errors.New("min should be stricly inferior to max")
	}
	if rng == nil {
		rng = newRand()
	}
	// Instantiate a GA
	var ga, err = GAConfig{
		NPops:        1,
		PopSize:      nParticles,
		NGenerations: nSteps,
		HofSize:      1,
		Model: ModMutationOnly{
			Strict: false,
		},
		ParallelEval: parallel,
		RNG:          rand.New(rand.NewSource(rng.Int63())),
	}.NewGA()
	if err != nil {
		return nil, err
	}
	return &SPSO{
		Min:   min,
		Max:   max,
		W:     w,
		BestY: math.Inf(1),
		GA:    ga,
		mutex: &sync.Mutex{},
	}, nil
}

// NewDefaultSPSO calls NewSPSO with default values.
func NewDefaultSPSO() (*SPSO, error) {
	return NewSPSO(40, 30, -5, 5, 0.5, false, nil)
}

// newParticle returns a new Particle that has a pointer to the SPSO.
func (pso *SPSO) newParticle(rng *rand.Rand) Genome {
	var (
		x        = InitUnifFloat64(pso.NDims, pso.Min, pso.Max, rng)
		velocity = make([]float64, len(x))
	)
	for i, xi := range x {
		min, max := pso.Min-xi, pso.Max-xi
		velocity[i] = min + rng.Float64()*(max-min)
	}
	return &Particle{
		CurrentX: x,
		BestX:    copyFloat64s(x),
		BestY:    math.Inf(1),
		SPSO:     pso,
		Velocity: velocity,
	}
}

// Minimize finds the minimum of a given real-valued function.
func (pso *SPSO) Minimize(f func([]float64) float64, nDims uint) ([]float64, float64, error) {
	// Set the function to minimize so that the particles can access it
	pso.F = f
	pso.NDims = nDims
	// Run the genetic algorithm
	var err = pso.GA.Minimize(pso.newParticle)
	// Return the best obtained vector along with the associated function value
	return pso.BestX, pso.BestY, err
}