From 236f45435733d6c827a21b460b04e2cbf5c9e2cb Mon Sep 17 00:00:00 2001
From: Jiri Kalvoda <jirikalvoda@kam.mff.cuni.cz>
Date: Wed, 24 May 2023 22:46:13 +0200
Subject: [PATCH] =?UTF-8?q?Algoritmus=20pomoc=C3=AD=20semidefinitn=C3=ADho?=
=?UTF-8?q?=20programov=C3=A1n=C3=AD?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
---
algo/semidef_prog.cpp | 175 ++++++++++++++++++++++++++++++++++++++++++
1 file changed, 175 insertions(+)
create mode 100644 algo/semidef_prog.cpp
diff --git a/algo/semidef_prog.cpp b/algo/semidef_prog.cpp
new file mode 100644
index 000000000..60d4b402b
--- /dev/null
+++ b/algo/semidef_prog.cpp
@@ -0,0 +1,175 @@
+#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;
+}
--
GitLab