MOGAwithRandomWalk.cpp

#include <cmath>
// Get gsl.hpp from http://code.google.com/p/igraphhpp/source/browse/trunk/gsl.
#include <gsl/gsl.hpp>
#include <ctime>
#include <algorithm>
#include <cstdio>
#include <tr1/unordered_map>

static const int N = 50;
static const int T = 150;
static const int ELITE_SIZE = N/5;
static const int MAX_WALK_STEPS = 3;

struct QuantumCoin {
	QuantumCoin() {
		std::memset(amplitude, 0, sizeof(amplitude));
	}
	
	QuantumCoin(int i, double val) {
		std::memset(amplitude, 0, sizeof(amplitude));
		amplitude[i] = val;
	}
		
	void randomize(gsl::Random& rng) {
		static gsl::SphericalVectorDistribution svd(32);
		
		double* vector = svd.get(rng);
		std::memcpy(amplitude, vector, sizeof(amplitude));
		delete[] vector;
	}
	
	void grover() {
		double sum = 0;
		
		// TODO: Use or find Kahan's summation algorithm.
		for (int i = 0; i < 32; ++ i)
			sum += amplitude[i];
			
		sum /= 16;
		
		for (int i = 0; i < 32; ++ i)
			amplitude[i] = sum - amplitude[i];
	}
	
	QuantumCoin& operator+=(const QuantumCoin& other) {
		for (int i = 0; i < 32; ++ i)
			amplitude[i] += other.amplitude[i];
		return *this;
	}
	
	void translate(std::tr1::unordered_map<unsigned long, QuantumCoin>& newStates, unsigned long oldPosition) const {
		for (int i = 0; i < 32; ++ i) {
			unsigned long newPosition = oldPosition ^ (1UL << i);
			newStates[newPosition] += QuantumCoin(i, amplitude[i]);
		}
	}
	
	double probability() const {
		double prob = 0;
		for (int i = 0; i < 32; ++ i)
			prob += amplitude[i]*amplitude[i];
		return prob;
	}
	
private:
	double amplitude[32];
};

struct QuantumWalker {
	QuantumWalker(unsigned long loc, gsl::Random& rng) {
		QuantumCoin coin;
		coin.randomize(rng);
		states.insert(std::make_pair(loc, coin));
	}
	
	void walk() {
		std::tr1::unordered_map<unsigned long, QuantumCoin> newStates;
		for (std::tr1::unordered_map<unsigned long, QuantumCoin>::iterator cit = states.begin(); cit != states.end(); ++ cit) {
			cit->second.grover();
		}
		for (std::tr1::unordered_map<unsigned long, QuantumCoin>::const_iterator cit = states.begin(); cit != states.end(); ++ cit) {
			cit->second.translate(newStates, cit->first);
		}
		states.swap(newStates);
	}
	
	unsigned long observe(gsl::Random& rng) const {
		double p = rng.uniform();
		double cur = 0;
		for (std::tr1::unordered_map<unsigned long, QuantumCoin>::const_iterator cit = states.begin(); cit != states.end(); ++ cit) {
			cur += cit->second.probability();
			if (cur >= p)
				return cit->first;
		}
		printf("****************** Oops! Can't observe anything. (%g, %g)\n", cur, p);
		return 0;
	}
	
private:
	std::tr1::unordered_map<unsigned long, QuantumCoin> states;
};

struct Entry {
	unsigned long x;
	double fx;
	
	void randomize(gsl::Random& rng) { x = rng.get(); }
	bool operator< (const Entry& other) const { return fx < other.fx; }
	void flipSpin(gsl::Random& rng) {
		int spin = rng.uniform_int(32);
		x ^= 1UL << spin;
	}
	double value() const {
		double xp = (double)x;
		xp /= 4294967296.0;
		xp *= 3;
		return xp - 1;
	}
	void evaluate() {
		double xp = value();
		fx = -1 - xp * sin(10*M_PI * xp);
	}
};

int main () {
	gsl::Random rng = gsl::Random::mt19937();
	rng.seed();
	
	// Initialize the random population.
	Entry population[N];
	for (int i = 0; i < N; ++ i)
		population[i].randomize(rng);
	
	for (int t = 0; t <= T; ++ t) {
		// Evaluate the population.
		for (int i = 0; i < N; ++ i)
			population[i].evaluate();
		// Sort the population to prepare for selection.
		std::sort(population, population+N);
	
		// Print the best solution.
		std::printf("%.16Lg\n", log10l((long double)population[0].fx - (-2.8502737667680982421L)));
		
		// Discard the last N/5 entries.
		std::memmove(population+ELITE_SIZE*2, population+ELITE_SIZE, (N-2*ELITE_SIZE)*sizeof(*population));
		// Copy the first N/5 entries (for simplicity let's not include order them).
		std::memcpy(population+ELITE_SIZE, population, ELITE_SIZE*sizeof(*population));
		
		// random walk.
		for (int i = 0; i < N; ++ i) {
			int walk_steps = MAX_WALK_STEPS * i / (N/8);
			if (walk_steps > MAX_WALK_STEPS)
				walk_steps = MAX_WALK_STEPS;
#if 0
			for (int n = 0; n < walk_steps; ++ n)
				population[i].flipSpin(rng);
#else
			QuantumWalker walker(population[i].x, rng);
			for (int n = 0; n < walk_steps; ++ n)
				walker.walk();
			population[i].x = walker.observe(rng);
#endif
		}
	}
	
	return 0;
}