diff --git a/Makefile b/Makefile index db3fb775..a516ec81 100644 --- a/Makefile +++ b/Makefile @@ -35,4 +35,6 @@ gokibin32mp: ${CXX} ${CXXFLAGS} ${MPFLAGS} -D_FLOAT_BITS_=32 -o gokibin32mp tools.cc gokibin64mp: ${CXX} ${CXXFLAGS} ${MPFLAGS} -D_FLOAT_BITS_=64 -o gokibin64mp tools.cc +gokibin128mp: + ${CXX} ${CXXFLAGS} ${MPFLAGS} -D_FLOAT_BITS_=128 -o gokibin128mp tools.cc diff --git a/goki.hh b/goki.hh index 370a8643..77832ffa 100644 --- a/goki.hh +++ b/goki.hh @@ -326,52 +326,6 @@ template SimpleMatrix sharpen(const int& size) { return s; } - - -template SimpleMatrix rotate(const SimpleMatrix& d, const T& theta) { - assert(abs(theta) < atan(T(1))); - const auto c(cos(theta)); - const auto s(sin(theta)); - const auto h0(abs(int(c * T(d.rows()) - s * T(d.cols())))); - const auto h1(h0 + abs(int(s * T(d.cols()))) * 2); - const auto w0(abs(int(s * T(d.rows()) + c * T(d.cols())))); - const auto w1(w0 + abs(int(s * T(d.rows()))) * 2); - SimpleMatrix res(h0 < d.rows() ? h1 : h0, - w0 < d.cols() ? w1 : w0); - const T offy(h0 < d.rows() ? abs(int(s * T(d.cols()))) : 0); - const T offx(w0 < d.cols() ? abs(int(s * T(d.rows()))) : 0); - res.O(); - const auto diag(int(sqrt(res.rows() * res.rows() + - res.cols() * res.cols())) + 1); - for(int j = - diag; j < diag; j ++) - for(int k = - diag; k < diag; k ++) { - const int yy(c * T(j) - s * T(k) + offy); - const int xx(s * T(j) + c * T(k) + offx); - if(0 <= yy && yy < res.rows() && - 0 <= xx && xx < res.cols()) { - const auto dyy(getImgPt(j, d.rows())); - const auto dxx(getImgPt(k, d.cols())); - { - res(yy, xx) = res(min(yy + 1, int(res.rows()) - 1), xx) = - res(yy, min(xx + 1, int(res.cols()) - 1)) = - res(min(yy + 1, int(res.rows()) - 1), - min(xx + 1, int(res.cols()) - 1)) = - d(dyy, dxx); - } - } - } - return res; -} - -template static inline SimpleMatrix center(const SimpleMatrix& dr, const SimpleMatrix& d) { - SimpleMatrix res(d.rows(), d.cols()); - for(int i = 0; i < res.rows(); i ++) - for(int j = 0; j < res.cols(); j ++) - res(i, j) = dr(max(0, min(i + (dr.rows() - d.rows()) / 2, dr.rows() - 1)), - max(0, min(j + (dr.cols() - d.cols()) / 2, dr.cols() - 1))); - return res; -} - template static inline SimpleMatrix flip(const SimpleMatrix& d) { auto res(d); #if defined(_OPENMP) diff --git a/lieonn.hh b/lieonn.hh index 91a473f3..c8a136d1 100644 --- a/lieonn.hh +++ b/lieonn.hh @@ -2502,16 +2502,18 @@ template static inline SimpleMatrix exp01(const SimpleMatrix& template static inline SimpleMatrix exp(const SimpleMatrix& m) { const auto p0(ceil(sqrt(norm2M(m)))); +#if defined(_FLOAT_BITS_) + // XXX: + myuint p(p0.operator myint()); +#else myuint p(p0); - if(T(int(1)) < abs(p0 - T(p))) throw "too large abs num in exp matrix"; - auto mm(exp01(m / T(p))); +#endif + if(T(myint(int(1))) < abs(p0 - T(myint(p)))) throw "too large abs num in exp matrix"; + auto mm(exp01(m / T(myint(p)))); auto res(m); - res.I(); - for( ; p; ) { - if(bool(p & myuint(int(1)))) res *= mm; - if(! (p >>= 1)) break; - mm *= mm; - } + for(res.I(); p; mm *= mm, p >>= 1) + if(bool(p & myuint(myint(int(1))))) + res *= mm; return res; } @@ -3639,6 +3641,44 @@ public: bool addp; }; +template class PpersistentOnce { +public: + inline PpersistentOnce() { ; } + inline ~PpersistentOnce() { ; } + inline const T& progression(const SimpleVector& h, const int& idx, const int& count) { + assert(0 <= idx && 0 <= count); + if(! count) return h[idx]; + if(ph[idx][count]) return eh[idx][count]; + ph[idx][count] = 1; + return (eh[idx][count] = progression(h, idx, count - 1) - progression(h, idx - 1, count - 1)); + } + inline T next(const SimpleVector& in0, const int& istat) { + const auto in(in0.size() & 1 ? in0.subVector(1, in0.size() - 1) : in0); + assert(! (in0.size() & 1) && 0 < istat); + { + vector eh0; + vector ph0; + eh0.resize(in.size() / 2, T(int(0))); + ph0.resize(in.size() / 2, false); + eh.resize(in.size(), eh0); + ph.resize(in.size(), ph0); + } + T res(int(0)); + for(int i = 0; i < in.size() / 2 - istat / 2; i ++) { + idFeeder buf(in.size() / 2 - i + istat / 2); + for(int j = i; j < in.size() / 2 + istat / 2; j ++) + buf.next(progression(in, j, i)); + // XXX align: + res += P(in.size() / 2 - istat / 2 - i).next(buf.res); + for(int j = i - 1, k = 1; 0 <= j; j --, k ++) + res += progression(in, in.size() - k, j); + } + return res /= T(in.size() / 2 - istat); + } + vector > eh; + vector > ph; +}; + template class PAthenB { public: inline PAthenB() { ; } @@ -4219,11 +4259,57 @@ template static inline T getImgPt(const T& y, const T& h) { return yy % h; } -template pair, SimpleVector > predv(const vector >& in) { +template SimpleMatrix rotate(const SimpleMatrix& d, const T& +theta) { + assert(abs(theta) < atan(T(1))); + const auto c(cos(theta)); + const auto s(sin(theta)); + const auto h0(abs(int(c * T(d.rows()) - s * T(d.cols())))); + const auto h1(h0 + abs(int(s * T(d.cols()))) * 2); + const auto w0(abs(int(s * T(d.rows()) + c * T(d.cols())))); + const auto w1(w0 + abs(int(s * T(d.rows()))) * 2); + SimpleMatrix res(h0 < d.rows() ? h1 : h0, + w0 < d.cols() ? w1 : w0); + const T offy(h0 < d.rows() ? abs(int(s * T(d.cols()))) : 0); + const T offx(w0 < d.cols() ? abs(int(s * T(d.rows()))) : 0); + res.O(); + const auto diag(int(sqrt(res.rows() * res.rows() + + res.cols() * res.cols())) + 1); + for(int j = - diag; j < diag; j ++) + for(int k = - diag; k < diag; k ++) { + const int yy(c * T(j) - s * T(k) + offy); + const int xx(s * T(j) + c * T(k) + offx); + if(0 <= yy && yy < res.rows() && + 0 <= xx && xx < res.cols()) { + const auto dyy(getImgPt(j, d.rows())); + const auto dxx(getImgPt(k, d.cols())); + { + res(yy, xx) = res(min(yy + 1, int(res.rows()) - 1), xx) = + res(yy, min(xx + 1, int(res.cols()) - 1)) = + res(min(yy + 1, int(res.rows()) - 1), + min(xx + 1, int(res.cols()) - 1)) = + d(dyy, dxx); + } + } + } + return res; +} + +template static inline SimpleMatrix center(const SimpleMatrix& dr, const SimpleMatrix& d) { + SimpleMatrix res(d.rows(), d.cols()); + for(int i = 0; i < res.rows(); i ++) + for(int j = 0; j < res.cols(); j ++) + res(i, j) = dr(max(int(0), min(i + (dr.rows() - d.rows()) / 2, dr.rows() - 1)), + max(int(0), min(j + (dr.cols() - d.cols()) / 2, dr.cols() - 1))); + return res; +} + +template pair >, vector > > predv(const vector >& in) { // N.B. we need to initialize p0 vector. SimpleVector init(3); for(int i = 0; i < init.size(); i ++) init[i] = T(int(i)); + cerr << (const char*)(persistent ? "persistent" : "sloppy") << endl; cerr << "Coherent: P0: " << P0maxRank0().next(init) << endl; SimpleVector seconds(in.size()); seconds.O(); @@ -4233,47 +4319,52 @@ template pair, SimpleVector > predv(const vector for(int i = 0; i < in.size(); i ++) { seconds[i] = makeProgramInvariant(in[i], - T(int(1)), true).second; } - // N.B. we need Pprogression<..., P01...> with maximum range. + // N.B. we need Ppersistent<..., P01...> with maximum range. const int istat(5 * 5 - 4 + 2); - const int p0((in.size() - istat) / 3); - SimpleVector p; - SimpleVector q; - p.resize(in[0].size()); - q.resize(in[0].size()); - p.O(); - q.O(); + const int p0((in.size() - istat) / 2); + vector > p; + vector > q; + p.resize(persistent ? 1 : p0, SimpleVector(in[0].size()).O()); + q.resize(persistent ? 1 : p0, SimpleVector(in[0].size()).O()); #if defined(_OPENMP) #pragma omp parallel for schedule(static, 1) #endif for(int j = 0; j < in[0].size(); j ++) { cerr << j << " / " << in[0].size() << endl; - Pprogression > pf(p0, istat); - Pprogression > pb(p0, istat); - for(int i = 0; i < in.size(); i ++) { - p[j] = pf.next(makeProgramInvariantPartial(in[i][j], - seconds[i], true)); - q[j] = pb.next(makeProgramInvariantPartial(in[in.size() - i - 1][j], - seconds[seconds.size() - i - 1], true)); + idFeeder buf(in.size()); + for(int i = 0; i < in.size(); i ++) + buf.next(makeProgramInvariantPartial(in[i][j], seconds[i], true)); + if(persistent) { + p[0][j] = PpersistentOnce >().next(buf.res, istat); + q[0][j] = PpersistentOnce >().next(buf.res.reverse(), istat); + } else for(int i = 0; i < p.size(); i ++) { + p[i][j] = P01(i + 1).next(buf.res); + q[i][j] = P01(i + 1).next(buf.res.reverse()); } } - Pprogression > pf(p0, istat); - Pprogression > pb(p0, istat); - T rp; - T rq; - for(int i = 0; i < in.size(); i ++) { - rp = pf.next(seconds[i]); - rq = pb.next(seconds[seconds.size() - 1 - i]); - } - for(int j = 0; j < p.size(); j ++) { - p[j] = p[j] == T(int(0)) || rp == T(int(0)) ? T(int(1)) / T(int(2)) - : revertProgramInvariant(make_pair(p[j], rp), true); - q[j] = q[j] == T(int(0)) || rp == T(int(0)) ? T(int(1)) / T(int(2)) - : revertProgramInvariant(make_pair(q[j], rq), true); + if(persistent) { + const auto rp(PpersistentOnce >().next(seconds, istat)); + const auto rq(PpersistentOnce >().next(seconds.reverse(), istat)); + for(int j = 0; j < p.size(); j ++) { + p[0][j] = p[0][j] == T(int(0)) || rp == T(int(0)) ? T(int(1)) / T(int(2)) + : revertProgramInvariant(make_pair(p[0][j], rp), true); + q[0][j] = q[0][j] == T(int(0)) || rp == T(int(0)) ? T(int(1)) / T(int(2)) + : revertProgramInvariant(make_pair(q[0][j], rq), true); + } + } else for(int i = 0; i < p0; i ++) { + const auto rp(P01(i + 1).next(seconds)); + const auto rq(P01(i + 1).next(seconds.reverse())); + for(int j = 0; j < p.size(); j ++) { + p[i][j] = p[i][j] == T(int(0)) || rp == T(int(0)) ? T(int(1)) / T(int(2)) + : revertProgramInvariant(make_pair(p[i][j], rp), true); + q[i][j] = q[i][j] == T(int(0)) || rp == T(int(0)) ? T(int(1)) / T(int(2)) + : revertProgramInvariant(make_pair(q[i][j], rq), true); + } } return make_pair(move(p), move(q)); } -template pair >, vector > > predVec(const vector > >& in0) { +template pair > >, vector > > > predVec(const vector > >& in0) { assert(in0.size() && in0[0].size() && in0[0][0].size()); vector > in; in.resize(in0.size()); @@ -4286,20 +4377,49 @@ template pair >, vector > > in[i].setVector(j * in0[i][0].size(), in0[i][j]); } } - const auto p(predv(in)); - pair >, vector > > res; - res.first.resize(in0[0].size()); - res.second.resize(in0[0].size()); - for(int j = 0; j < in0[0].size(); j ++) { - res.first[j] = - p.first.subVector(in0[0][0].size() * j, in0[0][0].size()); - res.second[j] = - p.second.subVector(in0[0][0].size() * j, in0[0][0].size()); + const auto p(predv(in)); + pair > >, vector > > > res; + res.first.resize(p.first.size()); + res.second.resize(p.second.size()); + for(int i = 0; i < res.first.size(); i ++) { + res.first[i].resize(in0[0].size()); + res.second[i].resize(in0[0].size()); + for(int j = 0; j < in0[0].size(); j ++) { + res.first[i][j] = + p.first[i].subVector(in0[0][0].size() * j, in0[0][0].size()); + res.second[i][j] = + p.second[i].subVector(in0[0][0].size() * j, in0[0][0].size()); + } } return res; } -template pair >, vector > > predMat(const vector > >& in0) { +template pair > >, vector > > > predMat(const vector > >& in0, const int& rot = 20) { + assert(0 <= rot); + if(0 < rot) { + auto res(predMat(in0, 0)); + if(rot <= 1) return res; + for(int i = 0; i < rot; i ++) { + cerr << "rotate(" << i << ")" << endl << flush; + const auto theta((T(i) - T(rot - 1) / T(2)) * atan(T(1)) / (T(rot) / T(2))); + auto ini(in0); + for(int j = 0; j < ini.size(); j ++) + for(int k = 0; k < ini[j].size(); k ++) + ini[j][k] = rotate(ini[j][k], theta); + auto resi(predMat(ini, 0)); + for(int j = 0; j < resi.first.size(); j ++) + for(int k = 0; k < resi.first[j].size(); k ++) { + res.first[j][k] += center(rotate(resi.first[j][k], - theta), res.first[j][k]); + res.second[j][k] += center(rotate(resi.second[j][k], - theta), res.second[j][k]); + } + } + for(int j = 0; j < res.first.size(); j ++) + for(int k = 0; k < res.first[j].size(); k ++) { + res.first[j][k] /= T(rot); + res.second[j][k] /= T(rot); + } + return res; + } assert(in0.size() && in0[0].size() && in0[0][0].rows() && in0[0][0].cols()); vector > in; in.resize(in0.size()); @@ -4314,28 +4434,32 @@ template pair >, vector > > k * in0[i][0].cols(), in0[i][j].row(k)); } } - const auto p(predv(in)); - pair >, vector > > res; - res.first.resize(in0[0].size()); - res.second.resize(in0[0].size()); - for(int j = 0; j < res.first.size(); j ++) { - res.first[j].resize(in0[0][0].rows(), in0[0][0].cols()); - res.second[j].resize(in0[0][0].rows(), in0[0][0].cols()); - for(int k = 0; k < in0[0][0].rows(); k ++) { - res.first[j].row(k) = - p.first.subVector( - j * in0[0][0].rows() * in0[0][0].cols() + k * in0[0][0].cols(), - in0[0][0].cols()); - res.second[j].row(k) = - p.second.subVector( - j * in0[0][0].rows() * in0[0][0].cols() + k * in0[0][0].cols(), - in0[0][0].cols()); + const auto p(predv(in)); + pair > >, vector > > > res; + res.first.resize(p.first.size()); + res.second.resize(p.second.size()); + for(int i = 0; i < res.first.size(); i ++) { + res.first[i].resize(in0[0].size()); + res.second[i].resize(in0[0].size()); + for(int j = 0; j < res.first[i].size(); j ++) { + res.first[i][j].resize(in0[0][0].rows(), in0[0][0].cols()); + res.second[i][j].resize(in0[0][0].rows(), in0[0][0].cols()); + for(int k = 0; k < in0[0][0].rows(); k ++) { + res.first[i][j].row(k) = + p.first[i].subVector( + j * in0[0][0].rows() * in0[0][0].cols() + k * in0[0][0].cols(), + in0[0][0].cols()); + res.second[i][j].row(k) = + p.second[i].subVector( + j * in0[0][0].rows() * in0[0][0].cols() + k * in0[0][0].cols(), + in0[0][0].cols()); + } } } return res; } -template pair, SimpleSparseTensor > predSTen(const vector >& in0, const vector& idx) { +template pair >, vector > > predSTen(const vector >& in0, const vector& idx) { assert(idx.size() && in0.size()); // N.B. we don't do input scaling. // N.B. the data we target is especially string stream corpus. @@ -4366,18 +4490,21 @@ template pair, SimpleSparseTensor > predST in[i][cnt ++] = (in0[i][idx[j]][idx[k]][idx[m]] + T(int(1))) / T(int(2)); } - auto p(predv(in)); + auto p(predv(in)); in.resize(0); - pair, SimpleSparseTensor > res; - for(int j = 0, cnt = 0; j < idx.size(); j ++) - for(int k = 0; k < idx.size(); k ++) - for(int m = 0; m < idx.size(); m ++) { - if(binary_search(absent.begin(), absent.end(), - make_pair(j, make_pair(k, m)))) - continue; - res.first[idx[j]][idx[k]][idx[m]] = p.first[cnt] * T(int(2)) - T(int(1)); - res.second[idx[j]][idx[k]][idx[m]] = p.second[cnt ++] * T(int(2)) - T(int(1)); - } + pair >, vector > > res; + res.first.resize(p.first.size()); + res.second.resize(p.second.size()); + for(int i = 0; i < res.first.size(); i ++) + for(int j = 0, cnt = 0; j < idx.size(); j ++) + for(int k = 0; k < idx.size(); k ++) + for(int m = 0; m < idx.size(); m ++) { + if(binary_search(absent.begin(), absent.end(), + make_pair(j, make_pair(k, m)))) + continue; + res.first[i][idx[j]][idx[k]][idx[m]] = p.first[i][cnt] * T(int(2)) - T(int(1)); + res.second[i][idx[j]][idx[k]][idx[m]] = p.second[i][cnt ++] * T(int(2)) - T(int(1)); + } return res; } diff --git a/test.py b/test.py index 76c6a839..5ea18ee4 100644 --- a/test.py +++ b/test.py @@ -142,9 +142,7 @@ if(argv[2] == "represent" or argv[2] == "collect" or argv[2] == "flarge" or argv[2] == "blink" or argv[2] == "enlarge" or argv[2] == "shrink" or argv[2] == "sharpen" or argv[2] == "limit" or argv[2] == "bit" or argv[2] == "rgb2xyz" or argv[2] == "xyz2rgb"): subprocess.call([argv[1], argv[2], root + ".ppm", root + "-" + argv[2] + ".ppm", str(pixels), str(rot)]) elif(argv[2] == "bump"): - subprocess.call([argv[1], argv[2], root + ".ppm", root + "-" + argv[2] + "0.ppm", str(pixels), str(rot)]) - subprocess.call(["python3", argv[0], argv[1], "cleanq", root + "-bump0.ppm"]) - subprocess.call(["convert", root + "-bump0-cleanq.png", "-format", "ppm", "-compress", "none", root + "-bump.ppm"]) + subprocess.call([argv[1], argv[2], root + ".ppm", root + "-" + argv[2] + ".ppm", "1", str(pixels)]) elif(argv[2] == "obj"): subprocess.call([argv[1], argv[2], root + "-bump.ppm", root + ".obj"]) elif(argv[2] == "jps"): @@ -176,7 +174,7 @@ subprocess.call([argv[1], "reshape", str(pixels), root + ".ppm", root + "-bump.ppm", root + "-illust.ppm", str(pow(float(pixels), - .5))]) elif(argv[2] == "gray"): subprocess.call(["convert", line, "-colorspace", "Gray", "-separate", "-average", root + "-gray.png"]) - elif(argv[2] == "cleanl" or argv[2] == "cleanlc" or argv[2] == "cleanL" or argv[2] == "cleanLc" or argv[2] == "cleanq"): + elif(argv[2] == "cleanl" or argv[2] == "cleanlc" or argv[2] == "cleanL" or argv[2] == "cleanLc"): w = h = 0 with open(root + ".ppm") as f: f.readline() @@ -186,10 +184,10 @@ # N.B. # With ddpmopt/README.md, we have 20kbit upper limit on function entropy. # In practical, 2 bit from MSB expecting input is enough. - upix = pow(19683 / 2., .5) + # However, working with bump command, we use 8 bit for each. + upix = pow(19683 / 8., .5) upixr = upix / pow(w * h, .5) if(argv[2][- 1] == "c"): upixr /= pow(3., .5) - sz = str(int(upixr * w) + 1) + "x" + str(int(upixr * h) + 1) # N.B. # refering en.wikipedia.org/wiki/Condorcet's-jury-theorem # n ~ 11 to get .95 from 2/3 probability. @@ -201,8 +199,6 @@ # From googling, via https://qiita.com/yoya/items/b1590de289b623f18639 . if(argv[2] == "cleanL" or argv[2] == "cleanLc"): subprocess.call(["convert", line, "-filter", "LanczosRadius", "-distort", "Resize", str(int(upixr * w * cj) + 1) + "x" + str(int(upixr * h * cj) + 1), root + "-" + argv[2] + ".png"]) - elif(argv[2] == "cleanl" or argv[2] == "cleanlc"): - subprocess.call(["convert", line, "-filter", "LanczosRadius", "-distort", "Resize", sz, root + "-" + argv[2] + ".png"]) else: - subprocess.call(["convert", line, "-despeckle", "-filter", "LanczosRadius", "-distort", "Resize", sz, "-despeckle", "+sigmoidal-contrast", "7.5", "-filter", "LanczosRadius", "-distort", "Resize", str(w) + "x" + str(h) + "!", "-sigmoidal-contrast", "7.5", root + "-cleanq.png"]) + subprocess.call(["convert", line, "-filter", "LanczosRadius", "-distort", "Resize", str(int(upixr * w) + 1) + "x" + str(int(upixr * h) + 1), root + "-" + argv[2] + ".png"])