From f9b899e6f07681f72fae8e9a8c722eb0bfaf356b Mon Sep 17 00:00:00 2001 From: =?utf8?q?Bj=C3=B8rn=20Rustad?= Date: Wed, 6 Dec 2017 22:53:05 +0100 Subject: [PATCH] LATEST --- ng.cpp | 93 +++++---- test.cpp | 587 ++++++++++++++++++++++++++++++++++++++++++++----------- 2 files changed, 532 insertions(+), 148 deletions(-) diff --git a/ng.cpp b/ng.cpp index de9da27..1164fa7 100644 --- a/ng.cpp +++ b/ng.cpp @@ -22,21 +22,20 @@ using namespace SIMDCompressionLib; double target_mult = 1e300; double bestest = 1e300; double total_target = 0.0; -vector> bsets; +vector> bsets; long long dbad = 0; long long dbad_limit = 0; long long without_success = 0; long long success_limit = 0; double min_score = 1e100; -vector > > ncombs; +vector > > ncombs; vector prods; vector primes; int n; +intersectionfunction inter; void test_slib() { int N = 1000; - intersectionfunction inter = - IntersectionFactory::getFromName("simd"); vector mydata1(N); vector mydata2(N); vector mydata3(N); @@ -45,7 +44,7 @@ void test_slib() { mydata2[i] = 2 * i; } clock_t start = clock(); - for (int j = 0; j < 100000; ++j) { + for (int j = 0; j < 1000; ++j) { size_t intersize = inter(mydata1.data(), mydata1.size(), mydata2.data(), mydata2.size(), mydata3.data()); @@ -59,7 +58,7 @@ void test_slib() { set s2(mydata2.begin(), mydata2.end()); set s3; start = clock(); - for (int j = 0; j < 100000; ++j) { + for (int j = 0; j < 1000; ++j) { set_intersection(s1.begin(), s1.end(), s2.begin(), s2.end(), inserter(s3, s3.begin())); } end = clock(); @@ -108,10 +107,10 @@ double get_target_mult(vector& primes, int n) { double upper = 100; double lower = 100; -multimap > combs; +multimap > combs; int siz = 0; -int try_comb(int placed, int start, double prod, bitset<36> b) { +int try_comb(int placed, int start, double prod, bitset<66> b) { //cout << "placed: " << placed << endl; if (placed == n-1) { if (prod > lower && prod < upper) { @@ -119,7 +118,9 @@ int try_comb(int placed, int start, double prod, bitset<36> b) { siz++; if (siz % 100 == 0) cout << "Added 1, now: " << siz << endl; return 0; - } + } else if (prod > upper) { + return 1; + } return 0; } else if (prod > upper) { return 1; @@ -135,7 +136,7 @@ int try_comb(int placed, int start, double prod, bitset<36> b) { } void smart_search() { - bitset<36> b; + bitset<66> b; try_comb(0, 0, 1.0, b); } @@ -146,7 +147,7 @@ void exhaustive_search(vector& primes, int n, double target) { // print integers and permute bitmask do { double prod = 1.0; - bitset<36> b; + bitset<66> b; for (int i = 0; i < primes.size(); ++i) // [0..N-1] integers { if (bitmask[i]) { @@ -186,7 +187,7 @@ double score_arr(vector& arr) { } -void print_bitset(bitset<36> b) { +void print_bitset(bitset<66> b) { double prod = 1.0; cout << "{"; int printed = 0; @@ -200,7 +201,7 @@ void print_bitset(bitset<36> b) { cout << "}"; } -double score_bitset(bitset<36> b) { +double score_bitset(bitset<66> b) { double prod = 1.0; for (int j = 0; j < b.size(); ++j) { if (b[j]) { @@ -234,16 +235,32 @@ void complete_rest( } } +int calls = 0; +clock_t start_call; +int max_mults = 0; + int find_best( double score, vector& arr, - set& allowed, - bitset<36> rest, - bitset<36> current, - bitset<36> not_allowed, + const vector& allowed, + const bitset<66>& rest, + const bitset<66>& current, + const bitset<66>& not_allowed, int maxi, int mults) { + calls++; + if (calls % 1000000000 == 0) { + calls = 0; + clock_t end = clock(); + cout << "1000000000 calls in " << double(end - start_call) / 1000000.0 << " max_mults: " << max_mults << endl; + start_call = clock(); + } + + if (mults > max_mults) { + max_mults = mults; + } + if (mults >= n) { double ss = score_arr(arr); if (ss - total_target < bestest) { @@ -276,17 +293,19 @@ int find_best( } return 0; } else { - for (auto it = allowed.begin(); it != allowed.end(); ++it) { + for (vector::const_iterator it = allowed.begin(); it != allowed.end(); ++it) { if (*it <= maxi) continue; if ((not_allowed & ncombs[*it].second).count() > 0) continue; - auto p = ncombs[*it]; + auto p = ncombs[*it].second; - set new_allowed; + vector new_allowed(min(allowed.size(), bsets[*it].size())); - set_intersection(allowed.begin(), allowed.end(), - bsets[*it].begin(), bsets[*it].end(), - inserter(new_allowed, new_allowed.begin())); + size_t intersize = inter(allowed.data(), allowed.size(), bsets[*it].data(), + bsets[*it].size(), new_allowed.data()); + + new_allowed.resize(intersize); + //new_allowed.shrink_to_fit(); arr.push_back(*it); @@ -294,9 +313,9 @@ int find_best( score + prods[*it], arr, new_allowed, - p.second ^ rest, - current | ncombs[*it].second, - not_allowed | (ncombs[*it].second & current), + p ^ rest, + current | p, + not_allowed | (p & current), *it, mults + 1); @@ -311,12 +330,12 @@ void get_sets() { double sum = 0; int outer = 0; for (auto it = ncombs.begin(); it != ncombs.end(); ++it) { - set b; + vector b; short i = -1; for (auto jt = ncombs.begin(); jt != ncombs.end(); ++jt) { ++i; if ((it->second & jt->second).count() == 1) { - b.insert(i); + b.push_back(i); } } sum += b.size(); @@ -344,18 +363,18 @@ void print_graph(multimap >& combs) { void print_all() { for (int i = 0; i < ncombs.size(); ++i) { for (auto it = bsets[i].begin(); it != bsets[i].end(); ++it) { - if (*it < i) - cout << i+1 << " " << *it+1 << endl; - } + if (*it < i) { + cout << (i+1) << " " << (*it) + 1 << endl; + } + } } } int main(int argc, char* argv[]) { n = atoi(argv[1]); primes = get_primes((n*(n-1))/2); + inter = IntersectionFactory::getFromName("simd"); test_slib(); - char c; - cin >> c; cout << "PRIMES:" << endl; cout << setprecision(90); for (auto it = primes.begin(); it != primes.end(); ++it) { @@ -371,7 +390,7 @@ int main(int argc, char* argv[]) { cout << "Number of primes: " << primes.size() << endl; cout << "Target node product: " << target_mult << endl; - //auto combs = exhaustive_search<36>(primes, n, target_mult); + //auto combs = exhaustive_search<66>(primes, n, target_mult); smart_search(); cout << "GOT COMBS: " << combs.size() << endl; @@ -394,12 +413,14 @@ int main(int argc, char* argv[]) { //print_graph(combs); cout << "GOTEM " << bsets.size() << endl; //print_all(); + start_call = clock(); vector arr; i = 0; - bitset<36> em; + bitset<66> em; + assert(em.size() == n*(n-1)/2); for (auto it = ncombs.begin(); it != ncombs.end(); ++it) { - if (i % 10 == 0) cout << i << endl; + cout << i << endl; arr.push_back(i); int ret = find_best(prods[i], arr, bsets[i], it->second, it->second, em, i, 1); arr.clear(); diff --git a/test.cpp b/test.cpp index a02d5bf..3a400ce 100644 --- a/test.cpp +++ b/test.cpp @@ -1,11 +1,20 @@ #include #include #include +#include #include #include +#include +#include +#include #include #include +using namespace std; + +double gtarget_row; + + class Edge { public: int u; @@ -15,67 +24,75 @@ public: Edge() {} }; -std::ostream& operator<<(std::ostream& o, Edge& e) { +ostream& operator<<(ostream& o, Edge& e) { o << e.u << "-" << e.v; return o; } +void initialize_primes(vector& primes, int n) { + int i = 0; + for (int p = 2; ; ++p) { + bool isprime = true; + for (int j = 0; j < i; ++j) { + if (p % primes[j] == 0) { + isprime = false; + break; + } + } + if (isprime) { + primes.push_back(p); + ++i; + } + if (i >= n*(n-1)/2) break; + } +} + + class Graph { public: - long long *primes; - long long n; - double *row_pd; - mpz_t *row_prod; - Edge *edges; + vector primes; + int n; + vector app_row_prod; + mpz_t *row_prod; + vector edges; + mpz_t total_target; + mpz_t row_target; + mpz_t row_target_adj; Graph(long long n) : n(n) { - row_prod = new mpz_t[n]; - row_pd = new double[n]; + row_prod = new mpz_t[n]; for (int i = 0; i < n; ++i) { - mpz_init(row_prod[i]); - } - int i = 0; - primes = new long long[n*(n-1)/2]; - for (int p = 2; ; ++p) { - bool isprime = true; - for (int j = 0; j < i; ++j) { - if (p % primes[j] == 0) { - isprime = false; - break; - } - } - if (isprime) { - primes[i] = p; - ++i; - } - if (i >= n*(n-1)/2) break; + mpz_init(row_prod[i]); } - edges = new Edge[n*(n-1)/2]; - int idx = 0; + initialize_primes(primes, n); + for (int i = 0; i < n; i++) { for (int j = i + 1; j < n; ++j) { - edges[idx] = Edge(i, j); - idx++; + edges.push_back(Edge(i, j)); } } - calc_row_prod(); + calc_row_prods(); + mpz_init(total_target); + mpz_init(row_target); + mpz_init(row_target_adj); + set_total_target(); } void print() { - std::cout << "THE PRIMES" << std::endl; + cout << "THE PRIMES" << endl; for (int i = 0; i < n*(n-1)/2; ++i) { - std::cout << primes[i] << std::endl; + cout << primes[i] << endl; } - std::cout << "THE EDGES" << std::endl; + cout << "THE EDGES" << endl; for (int i = 0; i < n*(n-1)/2; ++i) { - std::cout << edges[i] << std::endl; + cout << edges[i] << endl; } } - void calc_row_prod() { + void calc_row_prods() { for (int i = 0; i < n; ++i) { mpz_t mult; mpz_init(mult); @@ -89,7 +106,14 @@ public: } } - double score() { + void score(mpz_t *res) { + mpz_set_si(*res, 0); + for (int i = 0; i < n; ++i) { + mpz_add(*res, *res, row_prod[i]); + } + } + + double app_score() { mpz_t score; mpz_init(score); mpz_set_si(score, 0); @@ -102,31 +126,67 @@ public: return d; } + double app_score2() { + mpz_t score; + mpz_init(score); + mpz_set_si(score, 0); + for (int i = 0; i < n; ++i) { + mpz_add(score, score, row_prod[i]); + } + mpz_sub(score, score, total_target); + double d = mpz_get_d(score); + mpz_clear(score); + return d; + } + + vector min_max_indexes() { + double mmin = 1e300; + double mmax = 0; + double tt; + int maxi, mini; + vector ple; + for (int i = 0; i < n; ++i) { + tt = mpz_get_d(row_prod[i]); + if (tt > mmax) { + mmax = tt; + maxi = i; + } + if (tt < mmin) { + mmin = tt; + mini = i; + } + } + for (int j = 0; j < n*(n-1)/2; ++j) { + if (edges[j].u == maxi || edges[j].v == maxi + || edges[j].u == mini || edges[j].v == mini) + ple.push_back(j); + } + return ple; + } + void print_info() { - std::vector v; + vector v; for (int i = 0; i < n; ++i) { double d = mpz_get_d(row_prod[i]); v.push_back(d); - std::cout << d << ", "; + //cout << d << ", "; } - std::cout << std::endl; - double sum = std::accumulate(v.begin(), v.end(), 0.0); + //cout << endl; + double sum = accumulate(v.begin(), v.end(), 0.0); double mean = sum / v.size(); - double sq_sum = std::inner_product(v.begin(), v.end(), v.begin(), 0.0); - double stdev = std::sqrt(sq_sum / v.size() - mean * mean); - std::cout << "AVG: " << mean << std::endl; - std::cout << "STD: " << stdev << std::endl; - std::cout << "MIN: " << *std::min_element(v.begin(), v.end()) << std::endl; - std::cout << "MAX: " << *std::max_element(v.begin(), v.end()) << std::endl; + double sq_sum = inner_product(v.begin(), v.end(), v.begin(), 0.0); + double stdev = sqrt(sq_sum / v.size() - mean * mean); + cout << "AVG: " << mean << endl; + cout << "STD: " << stdev << endl; + cout << "MIN: " << *min_element(v.begin(), v.end()) << endl; + cout << "MAX: " << *max_element(v.begin(), v.end()) << endl; } - void ddiv(int i, long long div) { - //std::cout << "dividing by " << div << std::endl; + void ddiv(int i, int div) { mpz_divexact_ui(row_prod[i], row_prod[i], div); } - void mmul(int i, long long mul) { - //std::cout << "mul by " << mul << std::endl; + void mmul(int i, int mul) { mpz_mul_si(row_prod[i], row_prod[i], mul); } @@ -153,29 +213,263 @@ public: edges[to] = e; } + void jiggle_cycle(vector& cycle) { + Edge e = edges[cycle[0]]; + for (int i = 1; i < cycle.size(); ++i) { + int to = cycle[i]; + int from = cycle[i-1]; + mmul(edges[from].u, primes[to]); + mmul(edges[from].v, primes[to]); + ddiv(edges[from].u, primes[from]); + ddiv(edges[from].v, primes[from]); + edges[from] = edges[to]; + } + int to = cycle[0]; + int from = cycle.back(); + mmul(edges[from].u, primes[to]); + mmul(edges[from].v, primes[to]); + ddiv(edges[from].u, primes[from]); + ddiv(edges[from].v, primes[from]); + edges[from] = e; + } + + double optimize(int nn) { + //for (int i = 0; i < nn; ++i) { + // for (int from = 0; from < n*(n-1) / 2; ++from) { + // for (int to = from + 1; to < n*(n-1) / 2; ++to) { + // if (jiggle_score(from, to) < 0.0) { + // jiggle(from, to); + // } + // } + // } + //} + //for (int i = 0; i < nn; ++i) { + // for (int from = 0; from < n*(n-1) / 2; ++from) { + // for (int to = from + 1; to < n*(n-1) / 2; ++to) { + // for (int to2 = to + 1; to2 < n*(n-1) / 2; ++to2) { + // vector cycle; + // cycle.push_back(from); + // cycle.push_back(to); + // cycle.push_back(to2); + // if (jiggle_cycle_score(cycle) < 0.0) { + // jiggle_cycle(cycle); + // } + // } + // } + // } + //} + mpz_t oldscore, newscore; + mpz_init(oldscore); + mpz_init(newscore); + //for (int f = 0; f < n*(n-1) / 2; ++f) { + for (int i = 0; i < nn; ++i) { + int f = rand() % (n * (n-1)/2); + for (int t = 0; t < n*(n-1) / 2; ++t) { + for (int f2 = t; f2 < n*(n-1) / 2; ++f2) { + for (int t2 = f2; t2 < n*(n-1) / 2; ++t2) { + score(&oldscore); + jiggle(f, t); + jiggle(f2, t2); + score(&newscore); + mpz_sub(newscore, newscore, oldscore); + double diff = mpz_get_d(newscore); + if (diff >= 0.0) { + jiggle(t2, f2); + jiggle(t, f); + } + } + } + } + } + mpz_clear(oldscore); + mpz_clear(newscore); + } + + void set_fixed(vector& fixed) { + int ati = 0; + for (int i = 1; i < n; ++i) { + while (fixed[ati] == false) { + ati++; + } + int from; + int to = ati; + for (int tt = 0; tt < n*(n-1)/2; ++tt) { + if (edges[tt].u == 0 && edges[tt].v == i) { + from = tt; + break; + } + } + jiggle(from, to); + ati++; + } + } + + void print_table() { + for (int i = 0; i < edges.size(); ++i) { + cout << edges[i].u << " - " << primes[i] << " - " << edges[i].v << endl; + } + } + + void print_vector(vector& fixed) { + for (int i = 0; i < fixed.size(); ++i) { + if (fixed[i]) cout << primes[i] << ", "; + } + cout << endl; + } + + // take old product, divide by prime to remove + // then multiply by new prime + // subtract old value + void jiggle_diff_u(mpz_t *res, int index, int prev, int nnew) { + mpz_divexact_ui(*res, row_prod[edges[index].u], prev); + mpz_mul_si(*res, *res, nnew); + mpz_sub(*res, *res, row_prod[edges[index].u]); + } + + void jiggle_diff_v(mpz_t *res, int index, int prev, int nnew) { + mpz_divexact_ui(*res, row_prod[edges[index].v], prev); + mpz_mul_si(*res, *res, nnew); + mpz_sub(*res, *res, row_prod[edges[index].v]); + } + + // a - a * b / c + double jiggle_score(int from, int to) { + mpz_t res; + mpz_t sum; + mpz_init(res); + mpz_init(sum); + mpz_set_si(sum, 0); + jiggle_diff_u(&res, from, primes[from], primes[to]); + mpz_add(sum, sum, res); + jiggle_diff_v(&res, from, primes[from], primes[to]); + mpz_add(sum, sum, res); + jiggle_diff_u(&res, to, primes[to], primes[from]); + mpz_add(sum, sum, res); + jiggle_diff_v(&res, to, primes[to], primes[from]); + mpz_add(sum, sum, res); + double d = mpz_get_d(sum); + mpz_clear(res); + mpz_clear(sum); + return d; + } + + //TODO:DO this accurately maaan + void jiggle_diff_u2(mpz_t *res, int index, int prev, int nnew) { + mpz_t sum; + mpz_init(sum); + + mpz_divexact_ui(*res, row_prod[edges[index].u], prev); + mpz_mul_si(*res, *res, nnew); + + mpz_sub(*res, *res, row_target_adj); + mpz_abs(*res, *res); + + mpz_sub(sum, row_prod[edges[index].u], row_target_adj); + mpz_abs(sum, sum); + + mpz_sub(*res, *res, sum); + mpz_clear(sum); + } + + void jiggle_diff_v2(mpz_t *res, int index, int prev, int nnew) { + mpz_t sum; + mpz_init(sum); + + mpz_divexact_ui(*res, row_prod[edges[index].v], prev); + mpz_mul_si(*res, *res, nnew); + + mpz_sub(*res, *res, row_target_adj); + mpz_abs(*res, *res); + + mpz_sub(sum, row_prod[edges[index].v], row_target_adj); + mpz_abs(sum, sum); + + mpz_sub(*res, *res, sum); + mpz_clear(sum); + } + + // a - a * b / c + double jiggle_score2(int from, int to) { + mpz_t res; + mpz_t sum; + mpz_init(res); + mpz_init(sum); + mpz_set_si(sum, 0); + jiggle_diff_u2(&res, from, primes[from], primes[to]); + mpz_add(sum, sum, res); + jiggle_diff_v2(&res, from, primes[from], primes[to]); + mpz_add(sum, sum, res); + jiggle_diff_u2(&res, to, primes[to], primes[from]); + mpz_add(sum, sum, res); + jiggle_diff_v2(&res, to, primes[to], primes[from]); + mpz_add(sum, sum, res); + mpz_clear(res); + double d = mpz_get_d(sum); + mpz_clear(sum); + return d; + } + + + double jiggle_cycle_score(vector& cycle) { + mpz_t res; + mpz_t sum; + mpz_init(res); + mpz_init(sum); + mpz_set_si(sum, 0); + for (int i = 1; i < cycle.size(); ++i) { + int f = cycle[i]; + int t = cycle[i-1]; + jiggle_diff_u(&res, f, primes[f], primes[t]); + mpz_add(sum, sum, res); + jiggle_diff_v(&res, f, primes[f], primes[t]); + mpz_add(sum, sum, res); + } + int f = cycle[0]; + int t = cycle.back(); + jiggle_diff_u(&res, f, primes[f], primes[t]); + mpz_add(sum, sum, res); + jiggle_diff_v(&res, f, primes[f], primes[t]); + mpz_add(sum, sum, res); + double d = mpz_get_d(sum); + mpz_clear(res); + mpz_clear(sum); + return d; + } + + void set_total_target() { + mpz_set_si(total_target, 1); + for (int i = 0; i < n*(n-1) / 2; ++i) { + mpz_mul_si(total_target, total_target, primes[i]); + } + mpz_pow_ui(total_target, total_target, 2); + mpz_root(total_target, total_target, n); + mpz_set(row_target, total_target); + mpz_mul_si(total_target, total_target, n); + } + }; -std::ostream& operator<<(std::ostream& o, const Graph& g) { +ostream& operator<<(ostream& o, const Graph& g) { for (int i = 0; i < g.n; ++i) { - std::cout << "{"; + cout << "{"; int pp = 0; for (int j = 0; j < (g.n * (g.n-1)/2); ++j) { if (g.edges[j].u == i || g.edges[j].v == i) { - std::cout << g.primes[j]; + cout << g.primes[j]; pp++; if (pp < g.n-1) { - std::cout << ", "; + cout << ", "; } } } - std::cout << "}"; - if (i != g.n-1) std::cout << ", "; + cout << "}"; + if (i != g.n-1) cout << ", "; } return o; } -double get_target_mult(long long *primes, int n) { +double get_target_mult(vector& primes, int n) { mpz_t target_energy; mpz_init(target_energy); mpz_set_si(target_energy, 1); @@ -191,81 +485,150 @@ double get_target_mult(long long *primes, int n) { } - int main(int argc, char* argv[]) { - Graph g(atoi(argv[1])); srand(time(NULL)); - std::cout << g << std::endl; - std::cout << "SCORE: " << g.score() << std::endl; - std::cout << std::setprecision(10); - std::cout << "SHUFFLIGN" << std::endl; + cout << setprecision(10); + + Graph g(atoi(argv[1])); + + double target = get_target_mult(g.primes, atoi(argv[1])); + cout << "TARGET IS: " << target << endl; + gtarget_row = atof(argv[7]) + target / g.n; + cout << atof(argv[7]) << endl; + mpz_set_d(g.row_target_adj, atof(argv[7])); + mpz_add(g.row_target_adj, g.row_target_adj, g.row_target); + gmp_printf("Adjusted: %Zd\n", g.row_target_adj); + + cout << "SCORE: " << g.app_score2() << endl; + cout << "SHUFFLING" << endl; g.shuffle(10000000); - std::cout << "SHUFFLIGN DONE" << std::endl; - double cur = g.score(); - double min = g.score(); + cout << "SHUFFLING DONE" << endl; + + double cur = g.app_score2(); + double min = g.app_score2(); + double T = atof(argv[2]); - double target = get_target_mult(g.primes, atoi(argv[1])); - std::cout << "TARGET IS: " << target << std::endl; + int i = 0; int moves = 1; - g.print(); int fromi = 0; - std::vector take_from; + vector take_from; for (int ii = 0; ii < g.n*(g.n-1)/2; ++ii) { take_from.push_back(ii); } + int shuffle = 0; + int span = atoi(argv[4]); + vector scores; + vector ple = g.min_max_indexes(); + while (true) { - fromi++; - if (fromi >= g.n * (g.n-1) / 2) { - std::random_shuffle(take_from.begin(), take_from.end()); - fromi = 0; - } - int from = take_from[fromi]; - //int from = rand() % (g.n * (g.n-1) / 2); - //int to = rand() % (g.n * (g.n-1) / 2); - //int to2 = rand() % (g.n * (g.n-1) / 2); - int to = from + (rand() % 7); - int to2 = from - (rand() % 7); - if (to >= g.n * (g.n-1)/2) to = g.n * (g.n-1)/2 - 1; - if (to2 >= g.n * (g.n-1)/2) to2 = g.n * (g.n-1)/2 - 1; - if (to < 0) to = 0; - if (to2 < 0) to2 = 0; - if (to == to2) continue; - if (to == from) continue; - if (to2 == from) continue; - g.jiggle(from, to); - g.jiggle(from, to2); - double score = g.score(); - double change = log(score) - log(cur); - double acc_prob = exp(-(change)/T); + //fromi++; + //if (fromi >= g.n * (g.n-1) / 2) { + // if (shuffle % 10 == 0) + // random_shuffle(take_from.begin(), take_from.end()); + // shuffle++; + // fromi = 0; + //} + //int from = take_from[fromi]; + //if (double(rand()) / MAX_RAND > 0.4) { + //if (i % 10 == 0) { + //} + int from; + int to; + if (double(rand()) / RAND_MAX > 0.99) { + ple = g.min_max_indexes(); + from = ple[rand() % ple.size()]; + to = ple[rand() % ple.size()]; + } else { + from = rand() % (g.n*(g.n-1)/2); + to = from + rand() % span - (span-1) / 2; + } + if (0) { + from = rand() % (g.n*(g.n-1)/2); + for (int i = 0; i < g.n*(g.n-1)/2; ++i) { + if (i == from) continue; + //need to somehow pick a two other random edges + //it should preferrable be close + //one should go from.v == x.u => x.v which + //means x.v is higher index than x.u + //then we need the last edge from.u => x.v + //where again x.v is higher than from.u + } + } + + //TODO: SHOULD SUbtract target when we get score!! + //but since we mostly look at change it maybe doesn't matter so much... + + //int to = from + (rand() % span) - (span-1)/2; + if (to >= g.n * (g.n-1)/2) to = g.n * (g.n-1)/2 - 1; + if (to < 0) to = 0; + if (to == from) continue; + + double change = g.jiggle_score2(from, to); + double old_score = log(cur + 100000); + double new_score = log(cur + change + 100000); + double logchange = new_score - old_score; + double acc_prob = exp(-(logchange)/T); if (change < 0.0 || (acc_prob > double(rand()) / RAND_MAX)) { if (fabs(change) > 0.0000000000001) moves += 1; - cur = score; + g.jiggle(from, to); + cur += change; + if (moves % 1000 == 0) { + cur = g.app_score2(); + } + scores.push_back(cur); } else { - g.jiggle(to2, from); - g.jiggle(to, from); } + //if (double(rand()) / RAND_MAX > 0.9999999) { + // cout << "Tryin a little bit of optimize" << endl; + // g.optimize(5); + // cur = g.app_score(); + //} if (cur < min) { + cur = g.app_score2(); + } + //IDEA: heuristic for how many are below and above the target, added to the "loss" + //IDEA: loss that is more about getting each row close to the target + if (cur < min) { + //cout << "NEW MIN: " << cur - target << endl; min = cur; - std::cout << "NEW MIN: " << min - target << std::endl; - std::cout << g << std::endl; - std::cout << "temp = " << T << std::endl; - std::cout << "prob = " << acc_prob << std::endl; + cout << "NEW MIN: " << min << endl; + cout << g << endl; + cout << "temp = " << T << endl; + cout << "prob = " << acc_prob << endl; //g.print(); - g.print_info(); - } - T *= atof(argv[3]); - if (i % 10000000 == 0) { - if (moves <= 0) { - T /= pow(atof(argv[3]), 100000000); - std::cout << "No moves lately " << T << ", score:" << cur - target << std::endl; - std::cout << "Shuffling" << std::endl; - g.shuffle(10); - cur = g.score(); - } - moves = 0; + g.print_info(); + //if (min < 1e24) { + // cout << "OPTIMIZING" << endl; + // g.optimize(10); + // cur = g.app_score(); + // cout << "SCORE NOW: " << cur << endl; + // min = cur; + //} } + if (i % 5000000 == 0) { + cout << "Score: " << cur << ", min: " << min << ", T: " << T << endl; + } + int num_scores = 500000; + if (scores.size() == num_scores) { + double avgall = accumulate(scores.begin(), scores.end(), 0.0) / scores.size(); + double avgend = 0.0; + int nnum = 0; + for (int jj = num_scores * 2 / 3; jj < scores.size(); ++jj) { + avgend += scores[jj]; + nnum++; + } + avgend /= nnum; + //cout << "TEMP CHECK: all: " << avgall << " end: " << avgend << endl; + if (avgend > avgall) { + //cout << "OK reducing tempt from " << T << " to " << T * atof(argv[3]) << endl; + T *= atof(argv[3]); + } else { + //cout << "TEMP CHECK didnt make it: " << avgall << " > " << avgend << endl; + } + scores.clear(); + } i++; } return 0; -- 2.47.3