Algoritmus pomocí semidefinitního programování

algo/semidef_prog.cpp

+#include "../main.h"
+#include <stdlib.h>
+#include <assert.h>
+#include <cstring>
+#include <sdpa_call.h>
+#include <random>
+
+// ****** ARGS ********
+// <int>     number of tested cuts
+
+using std::pair;
+
+const int VERSION = 2;
+const bool ALLOW_MISTAKES = 0;
+
+template<class T>
+auto d2vec(int n, int m, T val=0)
+{
+	return vector<vector<T>>(n, vector<T>(m, val));
+}
+
+void printVector(vector<double> v, FILE * fpout)
+{
+		fprintf(fpout,"[ ");
+		for(auto jt:v)
+			fprintf(fpout,"%3d ",int(jt*100+0.5));
+		fprintf(fpout,"]; \n");
+}
+void printMatrix(vector<vector<double>> v, FILE * fpout)
+{
+	fprintf(fpout,"[\n");
+	for(auto it:v) printVector(it, fpout);
+	fprintf(fpout,"]; \n");
+}
+
+void printVector(vector<int> v, FILE * fpout)
+{
+		fprintf(fpout,"[ ");
+		for(auto jt:v)
+			fprintf(fpout,jt==0?" _ ":"%2d ", jt);
+		fprintf(fpout,"]; \n");
+}
+void printMatrix(vector<vector<int>> v, FILE * fpout)
+{
+	fprintf(fpout,"[\n");
+	for(auto it:v) printVector(it, fpout);
+	fprintf(fpout,"]; \n");
+}
+
+vector<bool> solve(int n, vector<int> input,
+		int argc, char const* const* argv)
+{
+	assert(argc==1);
+	int tested_cuts_count = atoi(argv[0]);
+	assert(tested_cuts_count > 0);
+
+	SDPA::printSDPAVersion(stderr);
+	SDPA sdpa;
+	sdpa.setDisplay(stderr);
+	sdpa.inputConstraintNumber(n);
+	sdpa.inputBlockNumber(1);
+	sdpa.inputBlockSize(1,n);
+	sdpa.inputBlockType(1,SDPA::SDP);
+
+	sdpa.initializeUpperTriangleSpace();
+
+	fo(i,n) sdpa.inputCVec(i+1, 1);
+	fo(i,n) sdpa.inputElement(i+1, 1, i+1, i+1,  1);
+
+	auto other = lib::gen_arr_other_car_of_type(n, input);
+
+	vector<vector<int> > c(n, vector<int>(n));
+	fo(i, 2*n-1)
+		c[input[i]][input[i+1]] += -1+2*((other[i]>i) == (other[i+1]>i+1));
+	fo(i,n) fo(j,i)
+		sdpa.inputElement(0, 1, 1+i, 1+j, c[i][j] + c[j][i]);
+
+
+	sdpa.initializeUpperTriangle();
+	sdpa.initializeSolve();
+
+	sdpa.solve();
+
+	auto m = d2vec<double>(n, n);
+	{
+		auto r = sdpa.getResultYMat(1);
+		fo(i,n) fo(j,n) m[i][j] = r[i*n + j];
+	}
+
+	//printMatrix(c, stderr);
+	//printMatrix(m, stderr);
+
+	auto sqrt_m = d2vec<double>(n, n);
+	fo(c, n)
+	{
+		{
+			double sum=0;
+			for(int i=c-1; i>=0; i--)
+				sum += sqrt_m[c][i] * sqrt_m[c][i];
+			sqrt_m[c][c] = sqrt(m[c][c] - sum);
+		}
+		for(int r=c+1; r<n; r++)
+		{
+			double sum=0;
+			for(int i=c-1; i>=0; i--)
+				sum += sqrt_m[r][i]*sqrt_m[c][i];
+			sqrt_m[r][c] = (m[r][c] - sum) / sqrt_m[c][c];
+		}
+	}
+	//printMatrix(sqrt_m, stderr);
+
+	vector<bool> out(2*n);
+
+	std::default_random_engine generator;
+	std::normal_distribution<double> distribution(0,1);
+
+	vector<double> rand_vector(n);
+
+	int best_score = 2*n; // +inf
+	fo(_, tested_cuts_count)
+	{
+		for (int i=0; i<n; ++i) {
+			rand_vector[i] = distribution(generator);
+		}
+		//printVector(rand_vector, stderr);
+
+		vector<double> multip(n);
+		fo(i,n) fo(j, n) multip[i] += sqrt_m[i][j]*rand_vector[j];
+		//printVector(multip, stderr);
+
+		vector<bool> current_out(2*n);
+		fo(i, 2*n)
+			current_out[i] = (other[i]>i) ^ (multip[input[i]]>=0);
+		int current_score=0;
+		fo(i, 2*n-1) current_score += current_out[i] != current_out[i+1];
+		if (current_score < best_score)
+		{
+			best_score = current_score;
+			out = current_out;
+		}
+	}
+
+	double val = 0;
+	double opt_score = 0;
+	double expected_score = 0;
+	fo(i,n) fo(j,n)
+	{
+		val += c[i][j] * m[i][j];
+	}
+	fo(i,2*n-1)
+	{
+		double scalar_product = m[input[i]][input[i+1]] * (-1+2*((other[i]>i) == (other[i+1]>i+1)));
+		if(scalar_product <-1) scalar_product =-1;
+		if(scalar_product > 1) scalar_product = 1;
+		double angle = acos(scalar_product);
+		double p = angle/M_PI;
+		//double opt_p = 0.5 - scalar_product/2;
+		// fprintf(stderr, "%3d: %3d %4d | %3d %3d\n", i, int(scalar_product*100+0.5), int(angle*100+0.5), int(p*100+0.5), int(opt_p*100+0.5));
+		expected_score += p;
+		//opt_score += opt_p;
+	}
+	lib::log_data("sdpa_primal_score", "%lf", sdpa.getPrimalObj());
+	lib::log_data("sdpa_dual_score", "%lf", sdpa.getDualObj());
+	//lib::log_data("opt_val", "%lf", 2*val);
+	//lib::log_data("opt_score", "%lf", opt_score);
+	lib::log_data("lower_bound", "%lf", (2*n-1)/2.0 - sdpa.getPrimalObj()/2/2);
+	lib::log_data("worst_case_score", "%lf", opt_score + 0.10527 * (2*n-1));
+	lib::log_data("expected_score", "%lf", expected_score);
+	fprintf(stderr, "opt_val: %lf\n", 2*val);
+	fprintf(stderr, "opt_score: %lf   %lf\n", (2*n-1 - val)/2, opt_score);
+	fprintf(stderr, "worst_case: %lf\n", opt_score + 0.10527 * (2*n-1));
+	fprintf(stderr, "expected_score: %lf\n", expected_score);
+
+	return out;
+}