2019-12-04, 16:17 | #89 | ||
Jan 2017
4F_{16} Posts |
Quote:
By such considerations, it's pretty obvious that ideally you'd want all row vectors to have about the same length (using bigger numbers for already long vectors has diminishing payoff). And dot products between rows to be minimal. I verified that the known best matrix for size 8 from OEIS does fulfill those conditions - near same vector lengths (but not exact - square sum from 202 to 208), and dot products between rows from 155 to 158. It didn't seem like Oleg's solution used much math beyond that. Quote:
Below is the optimization code implemented in Cython: Code:
cdef void add18(double a[18], double b[18], double c): cdef int i for i in range(18): a[i] += c * b[i] cdef int inverse(int mat[81], double r[9][18]): cdef int i, j, k cdef double det = 1 k = 0 for i in range(9): for j in range(9): r[i][j] = mat[k] k += 1 for i in range(9): for j in range(9, 18): r[i][j] = 0 for i in range(9): r[i][i+9] = 1 cdef double a for i in range(9): if abs(r[i][i]) < 1e-6: for j in range(i+1, 9): if abs(r[j][i]) > 2e-6: add18(r[i], r[j], 1) break else: return 0 a = 1/r[i][i] det *= r[i][i] for j in range(18): r[i][j] *= a r[i][i] = 1 # hopefully better numerical stability for j in range(9): if i == j: continue add18(r[j], r[i], -r[j][i]) return int(abs(det)+0.5) cdef void adjust_inverse(double r[9][18], int r1, int c1, int r2, int c2, int v): cdef int i for i in range(9): r[i][c1] += v * r[i][r1+9] for i in range(9): r[i][c2] -= v * r[i][r2+9] cdef double a if c1 == c2: a = 1/r[c1][c1] for i in range(c1, 18): r[c1][i] *= a for i in range(9): if i == c1: continue add18(r[i], r[c1], -r[i][c1]) return if abs(r[c1][c1]) < 1e-6: add18(r[c1], r[c2], 1) a = 1/r[c1][c1] for i in range(18): r[c1][i] *= a r[c1][c1] = 1 # more stable? for i in range(9): if i == c1: continue add18(r[i], r[c1], -r[i][c1]) a = 1/r[c2][c2] for i in range(18): r[c2][i] *= a r[c2][c2] = 1 for i in range(9): if i == c2: continue add18(r[i], r[c2], -r[i][c2]) def optim(arr, pairs_in): cdef int i, j, k, v cdef int m[81] cdef int div[81] cdef short pairs[3240] for i in range(9): for j in range(9): div[i*9+j] = i for i in range(81): m[i] = arr[i] for i in range(3240): pairs[i] = pairs_in[i] cdef double inv[9][18] cdef int det, d2 cdef double r det = inverse(m, inv) if det == 0: return 0 cdef int time, prev, p1, p2 cdef int r1, c1, r2, c2 prev = 3239 # If this is only improvement on first pass, not detected i = -1 while True: i += 1 if i == 3240: i = 0 if i == prev: break p1 = pairs[i] >> 8 p2 = pairs[i] & 255 if m[p1] == m[p2]: continue r1 = div[p1] c1 = p1 - 9*r1 r2 = div[p2] c2 = p2 - 9*r2 v = m[p2] - m[p1] r = (1+inv[c1][r1+9]*v)*(1+inv[c2][r2+9]*-v)+inv[c2][r1+9]*inv[c1][r2+9]*v*v d2 = int(det * abs(r) + .5) if d2 > det: m[p1] += v m[p2] -= v prev = i det = d2 adjust_inverse(inv, r1, c1, r2, c2, v); for i in range(81): arr[i] = m[i] return det Code:
#!/usr/bin/python3 import pyximport pyximport.install(language_level=3) from det_optim import optim import random B = [] for i in range(1, 10): B.extend([i]*9) pairs = [0]*3240 k = 0 for i in range(81): for j in range(i+1, 81): pairs[k] = (i << 8) + j k += 1 pl = [] for _ in range(100): random.shuffle(pairs) pl.append(tuple(pairs)) def unit(B, bestd): for _ in range(100): random.shuffle(B) x = B[:] for pairs in pl: B = x[:] d = optim(B, pairs) if d >= bestd: bestd = d print(d, B) assert d < 935000000 return bestd bestd = 0 #for _ in range(1): while True: bestd = unit(B, bestd) The optimization routine takes the order to try swaps in as an argument. The reason for that was that the optimization of a matrix is fast enough that using Python's random.shuffle() to generate a new random matrix was a significant portion of time used, and trying swaps in a different order seemed to work about as well as always generating a completely new random matrix. If you comment out the 'while True:' in the second file and uncomment the line above, you can benchmark it to see how long it takes to optimize 10000 matrices (note that the Cython code will be compiled to C, so the first run will use time for that, but it should be cached for subsequent runs). |
||
2019-12-04, 19:18 | #90 |
Romulan Interpreter
Jun 2011
Thailand
11^{2}·73 Posts |
Very Nice! Congrats to all.
I am a bit ashamed to say my solution this time hihi |
Thread Tools | |
Similar Threads | ||||
Thread | Thread Starter | Forum | Replies | Last Post |
November 2018 | Batalov | Puzzles | 5 | 2018-12-03 13:31 |
November 2017 | Batalov | Puzzles | 3 | 2017-12-08 14:55 |
November 2016 | Xyzzy | Puzzles | 1 | 2016-12-06 16:41 |
November 2015 | R. Gerbicz | Puzzles | 3 | 2015-12-01 17:48 |
November 2014 | Xyzzy | Puzzles | 1 | 2014-12-02 17:40 |