Fix bug in Droop quota set "method", add Droop region code to multiwi…

…nner.cc

This is still pretty ugly - there should be two multiwinner modes or
programs, one for VSE calculation and one for feasible region calculation.
But that will require more refactoring.

I'm not sure it's correct either - will need more testing against the
results from STV, which currently seem to be outside the feasible region
despite passing Droop proportionality.

CMakeLists.txt

@@ -208,7 +208,8 @@ add_library(qe_multiwinner_methods
 	src/multiwinner/rusty/aux/dsc.cc
 	src/multiwinner/shuntsstv.cc
 	src/multiwinner/stv.cc
-	src/stats/multiwinner/mwstats.cc)
+	src/stats/multiwinner/mwstats.cc
+	src/stats/multiwinner/vse_region_mapper.cc)
 
 add_library(quadelect_lib src/common/ballots.cc
 	src/bandit/lilucb.cc

src/main/multiwinner.cc

@@ -15,7 +15,7 @@
 
 #include "stats/multiwinner/mwstats.h"
 #include "stats/multiwinner/vse.h"
-//#include "stats/multiwinner/vse_region_mapper.h"
+#include "stats/multiwinner/vse_region_mapper.h"
 
 #include "multiwinner/helper/errors.cc"
 
@@ -406,12 +406,14 @@ void get_limits(int num_candidates, int council_size,
 
 // This is currently very hacky.
 
-/*std::vector<VSE_point>*/ void get_droop_council_results(
+std::vector<VSE_point> get_droop_council_results(
 	const election_t & election,
 	size_t num_candidates, size_t council_size,
 	const std::vector<vector<bool> > & population_profiles,
 	const std::vector<double> & whole_pop_profile,
-	const std::vector<double> & utilities) {
+	const std::vector<double> & utilities,
+	double best_disprop, double random_disprop,
+	double best_utility, double random_utility) {
 
 	quota_helper droop_sifter(council_size);
 
@@ -426,6 +428,8 @@ void get_limits(int num_candidates, int council_size,
 	std::vector<std::vector<size_t> > output = optimum.
 		get_optimal_solutions();
 
+	std::vector<VSE_point> droop_councils;
+
 	for (std::vector<size_t> & council: output) {
 		// TODO: Deal with this annoying thing where it only accepts
 		// lists, later. I might just change the signature so that
@@ -442,9 +446,20 @@ void get_limits(int num_candidates, int council_size,
 		std::cout << "Valid council: ";
 		std::copy(converted.begin(), converted.end(),
 			std::ostream_iterator<size_t>(std::cout, " "));
-		std::cout << " disproportionality " << disprop << ", utility "
-			<< utility << "\n";
+
+		VSE disprop_VSE;
+		disprop_VSE.add_result(random_disprop, disprop, best_disprop);
+		VSE utility_VSE;
+		utility_VSE.add_result(random_utility, utility, best_utility);
+
+		std::cout << " disprop. " << disprop << "(VSE " << disprop_VSE.get() <<
+			"),"
+			" utility " << utility << "(VSE " << utility_VSE.get() << ")\n";
+
+		droop_councils.push_back({disprop_VSE, utility_VSE});
 	}
+
+	return droop_councils;
 }
 
 // n * tries instead of plain n(as we would have with a vector of ballots.
@@ -788,6 +803,8 @@ std::vector<multiwinner_stats> get_multiwinner_methods() {
 
 int main(int argc, char * * argv) {
 
+	vse_region_mapper region_mapper(200);
+
 	// Used for inlining.
 	int maxnum = 0;
 	bool run_forever = false;
@@ -891,15 +908,24 @@ int main(int argc, char * * argv) {
 
 		double cur_worst_disprop, cur_mean_disprop, cur_best_disprop;
 
-		get_droop_council_results(ballots,
-			num_candidates, council_size,
-			population_profiles, pop_opinion_profile,
-			utilities);
-
 		get_limits(num_candidates, council_size, population_profiles,
 			pop_opinion_profile, 75000, cur_worst_disprop,
 			cur_mean_disprop, cur_best_disprop);
 
+		// Calculate the feasible region for Droop proportionality
+		// (almost; it's slightly too generous). TODO: Explain why.
+		// And perhaps put it somewhere else once I've refactored
+		// properly.
+
+		region_mapper.add_points(get_droop_council_results(ballots,
+				num_candidates, council_size,
+				population_profiles, pop_opinion_profile,
+				utilities, cur_best_disprop, cur_mean_disprop,
+				cur_best_utility, cur_random_utility));
+		region_mapper.update();
+		std::ofstream coords_out("droop_coords.txt");
+		region_mapper.dump_coordinates(coords_out);
+
 		// The stuff below should use proper multiwinner_stats objects.
 		// TODO!!!
 

src/multiwinner/exhaustive/quota_helper.h

@@ -45,16 +45,27 @@ class quota_helper : public exhaustive_method {
 double quota_helper::evaluate(combo::it & start, combo::it & end) {
 	std::set<size_t> council(start, end);
 
+	// Assumes that we can get the total weight from the first solid
+	// coalition, i.e. that they're in sorted order. TODO, fix so that
+	// if I change this later, it will trigger instead of silently fail.
+	// Or get the total weight (num voters) from somewhere else.
+
+	double total_weight = solid_coalitions.begin()->support;
+
 	if (council_size == 0) {
 		throw std::runtime_error("Quota helper: council size not set!");
 	}
 
 	for (auto cc_pos = solid_coalitions.begin();
 		cc_pos != solid_coalitions.end(); ++cc_pos) {
 
+		// See PSC-CLE.
+		size_t minimum = ceil(cc_pos->support *
+				(council_size + 1) / total_weight - 1);
+
 		size_t elect_constraint = std::min(
-				(double)cc_pos->coalition.size(),
-				cc_pos->support/(council_size + 1.0));
+				cc_pos->coalition.size(),
+				minimum);
 
 		if (elect_constraint == 0) {
 			continue;

src/stats/multiwinner/vse.h

@@ -1,5 +1,9 @@
 #pragma once
 
+#include <vector>
+#include <cstddef>
+#include <stdexcept>
+
 // VSE container.
 
 // This works both for scores and penalties; it doesn't matter if higher values

src/stats/multiwinner/vse_region_mapper.cc

@@ -1,13 +1,19 @@
 #include "vse_region_mapper.h"
+#include "tools/tools.h"
+
+#include <cmath>
+#include <vector>
+
+#include <iostream>
 
 std::vector<double> vse_region_mapper::get_quantized_position(
 	const VSE_point & pt) const {
 
 	std::vector<double> quantized_VSEs;
 
 	for (const VSE & v: pt) {
-		quantized_VSEs.push_back(round(
-				pt.get()*elements_per_dimension)/(double)elements_per_dimension);
+		quantized_VSEs.push_back(round(v.get()*
+				elements_per_dimension)/(double)elements_per_dimension);
 	}
 
 	return quantized_VSEs;
@@ -18,43 +24,67 @@ double vse_region_mapper::get_quantization_error(
 
 	double squared_error = 0;
 
-	for (const VSE & v: pt) {
-		double quantized = round(pt.get()*elements_per_dimension)/
-			(double)elements_per_dimension;
+	for (const VSE & coordinate: pt) {
+		double quantized = round(coordinate.get()*
+				elements_per_dimension)/(double)elements_per_dimension;
 
-		squared_error += square(quantized-pt.get());
+		squared_error += square(quantized-coordinate.get());
 	}
 
 	return sqrt(squared_error);
 }
 
-void vse_region_mapper::add_point(const VSE_point & cur_round_VSE) {
-	for (auto pos = VSE_cloud.begin(); pos != VSE_cloud.end(); ++pos) {
+void vse_region_mapper::add_augmented_point(
+	const VSE_point & augmented_point) {
 
-		// Copy the cloud point, add the new data, and insert.
-		VSE new_cloud_point = pos->second;
-		new_cloud_point.add_last(cur_round_VSE);
+	std::vector<double> quant_pos =
+		get_quantized_position(augmented_point);
+
+	if (next_VSE_cloud.find(quant_pos) == next_VSE_cloud.end()) {
+		next_VSE_cloud[quant_pos] = augmented_point;
+	} else {
+		VSE_point incumbent = next_VSE_cloud.find(quant_pos)->second;
 
-		std::vector<double> quant_pos =
-			get_quantized_position(new_cloud_point);
+		if (get_quantization_error(augmented_point) <
+			get_quantization_error(incumbent)) {
 
-		if (next_VSE_cloud.find(quant_pos) == next_VSE_cloud.end()) {
-			next_VSE_cloud[quant_pos] = new_cloud_point;
-		} else {
-			VSE incumbent = next_VSE_cloud.find(quant_pos)->second;
+			next_VSE_cloud[quant_pos] = augmented_point;
+		}
+	}
+}
 
-			if (get_quantization_error(new_cloud_point) <
-				get_quantization_error(incumbent)) {
+void vse_region_mapper::add_point(const VSE_point & cur_round_VSE) {
 
-				next_VSE_cloud[quant_pos] = new_cloud_point;
-			}
+	// If the cloud is empty, just add the point directly.
+	if (VSE_cloud.empty()) {
+		add_augmented_point(cur_round_VSE);
+		return;
+	}
+
+	// Otherwise, augment every existing cloud point from the previous
+	// round with this point, and add them all to the new cloud.
+	for (auto pos = VSE_cloud.begin(); pos != VSE_cloud.end(); ++pos) {
+
+		// Copy the cloud point, add the new data, and insert.
+		VSE_point new_cloud_point = pos->second;
+		for (size_t i = 0; i < new_cloud_point.size(); ++i) {
+			new_cloud_point[i].add_last(cur_round_VSE[i]);
 		}
+
+		add_augmented_point(new_cloud_point);
+	}
+}
+
+void vse_region_mapper::add_points(const std::vector<VSE_point> & points) {
+	for (const VSE_point & p: points) {
+		add_point(p);
 	}
 }
 
-void vse_region_mapper:update() {
+void vse_region_mapper::update() {
 	VSE_cloud = filter_augmented_points(next_VSE_cloud);
 	next_VSE_cloud.clear();
+	std::cout << "Update: now we have " << VSE_cloud.size() << " points.\n";
 }
 
 std::vector<VSE_point> vse_region_mapper::get() const {
@@ -65,4 +95,19 @@ std::vector<VSE_point> vse_region_mapper::get() const {
 	}
 
 	return out;
+}
+
+void vse_region_mapper::dump_coordinates(std::ostream & where) const {
+
+	for (auto pos = VSE_cloud.begin(); pos != VSE_cloud.end(); ++pos) {
+		bool first = true;
+		for (const VSE & v: pos->second) {
+			if (!first) {
+				where << " ";
+			}
+			where << v.get();
+			first = false;
+		}
+		where << "\n";
+	}
 }

src/stats/multiwinner/vse_region_mapper.h

@@ -14,6 +14,9 @@
 // binned in n-space to produce "voxels"; multiple points within the same
 // hypercube are discarded before the new point cloud is admitted.
 
+// Even so, 2D convolution with lots of points can be pretty expensive.
+// Don't use too high an elements per dimension value.
+
 #include "vse.h"
 
 #include <map>
@@ -43,19 +46,23 @@ class vse_region_mapper {
 		double get_quantization_error(
 			const VSE_point & pt) const;
 
-		// For enforcing additional constraints, e.g. Pareto frontier.
+		// For enforcing additional constraints, e.g. convex hull or
+		// Pareto frontier.
 		virtual std::map<std::vector<double>, VSE_point>
 		filter_augmented_points(const
 			std::map<std::vector<double>, VSE_point> & proposed) const {
 			return proposed;
 		}
 
+		void add_augmented_point(const VSE_point & augmented_point);
+
 	public:
 		// Call this with a vector of VSEs; the mapper will use
 		// their last round results. Each n-vector defines a new
 		// point in n-space, so this would be called once per
 		// new point.
 		void add_point(const VSE_point & cur_round_VSE);
+		void add_points(const std::vector<VSE_point> & points);
 		// Then call this to update the cloud.
 		void update();
 
@@ -64,9 +71,7 @@ class vse_region_mapper {
 
 		vse_region_mapper(int elements_per_dim_in) {
 			elements_per_dimension = elements_per_dim_in;
-
-			// Add a null point so that we don't have to special-
-			// case a bunch of stuff in add_point.
-			VSE_cloud[ {}] = VSE();
 		}
+
+		void dump_coordinates(std::ostream & where) const;
 };