kmyk / competitive-programming-library

kimiyuki's library for competitive programming with C++
https://kmyk.github.io/competitive-programming-library
MIT License
66 stars 9 forks source link

one more question about segment trees #3

Open bethandtownes opened 5 years ago

bethandtownes commented 5 years ago

Hi,

I have one more question using your segment tree library. I am trying to use it to solve this problem (https://cn.vjudge.net/problem/HDU-1542). Basically, it is a sweepline + segment tree with coordinate compression in the process to find out the area of union of rectangles.

I chose the full lazy segment tree to attack this problem. However, I could not solve the problem and maintain the integrity of template code at the same time, i.e., I have to specialized the code a little bit in order to solve the problem, mainly due to coordinate compression.

Could you, if possible, solve this problem without changing the template code? I really appreciate it!

Jason

bethandtownes commented 5 years ago

PS: my email address is ds653@cornell.edu:) looking forward to your constructive insights!

I learnt quite a bit from your way of programming.

kmyk commented 5 years ago

Coordinate compression on segment trees can be expressed without modification of the code of segment tree, by adding information of ranges to elements of the monoid. It is like following:

struct extended_monoid {
    typedef struct {
        int l, r;
        base_monoid m;
    } underlying_type;
    underlying_type append(underlying_type a, underlying_type b) {
        underlying_type c;
        c.l = a.l;
        c.r = a.r;
        ...
    }
    ...
};

The given problem can be solved in O(n^2 log n) by simple monoid below, if my understanding is correct. Also note that, it seems we can solve this by O(n^2) using imos's method (いもす法 https://imoz.jp/algorithms/imos_method.html ; I don't know how this technique is called in English, but it's very famous in Japanese community).

struct count_monoid {
    array<double, 101> data;
    ....
};
struct increment_operator_monoid {
    ...
    target_type apply(int a, array<double, 101> b) const {
        if (a == 0) return b;
        array<double, 101> c = {};
        REP (i, N - a) c[i + a] = b[i];
        return c;
    }
};
kmyk commented 5 years ago

By the way, if you know that how the imos's method is called in English, please tell me that. I couldn't find it.

bethandtownes commented 5 years ago

Could you give me a demo on how to use the extended monoid? I could not understand how to update the range information in the range_apply function. I also don’t understand how to incorporate them with Magma Endomorphisms...

kmyk commented 5 years ago

Here is a full code of a demo, applying a linear function over a compressed range and computing the sum.

#include <bits/stdc++.h>
#define REP(i, n) for (int i = 0; (i) < (int)(n); ++ (i))
#define REP3(i, m, n) for (int i = (m); (i) < (int)(n); ++ (i))
#define REP_R(i, n) for (int i = int(n) - 1; (i) >= 0; -- (i))
#define REP3R(i, m, n) for (int i = int(n) - 1; (i) >= (int)(m); -- (i))
#define ALL(x) begin(x), end(x)
using ll = long long;
using namespace std;

template <class Monoid, class OperatorMonoid>
struct lazy_propagation_segment_tree { // on monoids
    static_assert (is_same<typename Monoid::underlying_type, typename OperatorMonoid::target_type>::value, "");
    typedef typename Monoid::underlying_type underlying_type;
    typedef typename OperatorMonoid::underlying_type operator_type;
    const Monoid mon;
    const OperatorMonoid op;
    int n;
    vector<underlying_type> a;
    vector<operator_type> f;
    lazy_propagation_segment_tree() = default;
    lazy_propagation_segment_tree(int a_n, underlying_type initial_value = Monoid().unit(), Monoid const & a_mon = Monoid(), OperatorMonoid const & a_op = OperatorMonoid())
            : mon(a_mon), op(a_op) {
        n = 1; while (n <= a_n) n *= 2;
        a.resize(2 * n - 1, mon.unit());
        fill(a.begin() + (n - 1), a.begin() + ((n - 1) + a_n), initial_value); // set initial values
        REP_R (i, n - 1) a[i] = mon.append(a[2 * i + 1], a[2 * i + 2]); // propagate initial values
        f.resize(max(0, (2 * n - 1) - n), op.identity());
    }
    void point_set(int i, underlying_type z) {
        assert (0 <= i and i < n);
        point_set(0, 0, n, i, z);
    }
    void point_set(int i, int il, int ir, int j, underlying_type z) {
        if (i == n + j - 1) { // 0-based
            a[i] = z;
        } else if (ir <= j or j+1 <= il) {
            // nop
        } else {
            range_apply(2 * i + 1, il, (il + ir) / 2, 0, n, f[i]);
            range_apply(2 * i + 2, (il + ir) / 2, ir, 0, n, f[i]);
            f[i] = op.identity();
            point_set(2 * i + 1, il, (il + ir) / 2, j, z);
            point_set(2 * i + 2, (il + ir) / 2, ir, j, z);
            a[i] = mon.append(a[2 * i + 1], a[2 * i + 2]);
        }
    }
    void range_apply(int l, int r, operator_type z) {
        assert (0 <= l and l <= r and r <= n);
        range_apply(0, 0, n, l, r, z);
    }
    void range_apply(int i, int il, int ir, int l, int r, operator_type z) {
        if (l <= il and ir <= r) { // 0-based
            a[i] = op.apply(z, a[i]);
            if (i < f.size()) f[i] = op.compose(z, f[i]);
        } else if (ir <= l or r <= il) {
            // nop
        } else {
            range_apply(2 * i + 1, il, (il + ir) / 2, 0, n, f[i]);
            range_apply(2 * i + 2, (il + ir) / 2, ir, 0, n, f[i]);
            f[i] = op.identity();
            range_apply(2 * i + 1, il, (il + ir) / 2, l, r, z);
            range_apply(2 * i + 2, (il + ir) / 2, ir, l, r, z);
            a[i] = mon.append(a[2 * i + 1], a[2 * i + 2]);
        }
    }
    underlying_type range_concat(int l, int r) {
        assert (0 <= l and l <= r and r <= n);
        return range_concat(0, 0, n, l, r);
    }
    underlying_type range_concat(int i, int il, int ir, int l, int r) {
        if (l <= il and ir <= r) { // 0-based
            return a[i];
        } else if (ir <= l or r <= il) {
            return mon.unit();
        } else {
            return op.apply(f[i], mon.append(
                    range_concat(2 * i + 1, il, (il + ir) / 2, l, r),
                    range_concat(2 * i + 2, (il + ir) / 2, ir, l, r)));
        }
    }
};

struct plus_monoid_with_range {
    typedef struct {
        int l, r;
        ll value;
    } underlying_type;
    underlying_type unit() const {
        return (underlying_type) { -1, -1, 0ll };
    }
    underlying_type append(underlying_type a, underlying_type b) const {
        underlying_type c;
        c.l = a.l;
        c.r = b.r;
        c.value = a.value + b.value;
        return c;
    }
};
struct linear_operator_monoid_with_range {
    typedef pair<ll, ll> underlying_type;  // \lambda x. ax + b
    typedef plus_monoid_with_range::underlying_type target_type;
    underlying_type identity() const {
        return make_pair(1ll, 0ll);
    }
    target_type apply(underlying_type f, target_type x) const {
        ll a, b; tie(a, b) = f;
        target_type y = x;
        y.value = a * x.value + b * (x.r - x.l);
        return y;
    }
    underlying_type compose(underlying_type f, underlying_type g) const {
        return make_pair(f.first * g.first, f.second + f.first * g.second);
    }
};

int main() {
    vector<ll> ary(100);
    ary[0] = 1;
    REP (i, 100 - 1) {
        ary[i + 1] = ary[i] * 3 % 1000000007;
    }

    lazy_propagation_segment_tree<plus_monoid_with_range, linear_operator_monoid_with_range> segtree(10);
    REP (i, 10) {
        int l = 10 * i;
        int r = 10 * (i + 1);
        ll x = accumulate(ary.begin() + l, ary.begin() + r, 0ll);
        segtree.point_set(i, (plus_monoid_with_range::underlying_type) { l, r, x });
    }

    REP3 (i, 20, 60) {
        ary[i] = 3 * ary[i] + 5;
    }
    segtree.range_apply(2, 6, make_pair(3, 5));

    assert (accumulate(ALL(ary), 0ll) == segtree.range_concat(0, 10).value);
    return 0;
}
bethandtownes commented 5 years ago

thank you so much for taking the time to answer my questions.

I have a few follow up questions: where did you pick these techniques up? I can sense that you are a mathematician since your code totally has that flavor and also that you are an AMAZING functional programmer. I am trying to pick these up, could you recommend me some sources?

Second questions is regarding to the monoids: suppose I want a Sgtree that supports multiple operations(more than two), the way to design it per your template would be make a product monoids I guess? Assuming there is little interaction between each of the coordinate monoids(i.e. they are independent (at least weakly)).

PS. That technique you introduced me seems to have very little english documentation. I could not track it down with my effort. I think that person develops it so much that it goes beyond competitive programming scope and it also has something to do with his master thesis (there is an youtube video of his talking about it and his master research). If I understand it correctly(aka. google translate did a decent job:)), it involves some kind of approximation theory. In the website you showed me, he used it to approximate Gaussian functions. Coming from a mathematical background, I suppose if you would agree that if something can be related to gaussian family, then its implications are most often far reaching mostly because gaussian stuff is an quite important branch of mathematical analysis. I guess it would be interesting to see more theoretical results about this technique and also to have its pure mathematical framework well documented in English.

kmyk commented 5 years ago

About these techniques, I did nothing particular. I found the technique of monoids just as I solve a problem of competitive programming. Also, I am studying mathematics now, but almost all of mathematical tools in this library are found by googling. Therefore, I'm sorry that, I cannot recommend any particular source. Basic skills of mathematics are useful, so I recommend reading a book of mathematical logic and a book of algebras, if you are not familiar to them.

About multiple operations: Yes. You can use product monoids for this purpose. I note that the technique extending by range can be described as a product monoid of a base monoid and the monoid of ranges.

Thank you for researching for imos's method. I am surprised that this is not known as named (it means typical) technique in English. I had forgot that it can be said as inverse of cumulative sums or integral of discrete differential, and I should have explained using such words. I don't have much knowledge for Gaussian functions, but I think Gaussian functions may be just one of the applications. Not only Gaussian functions, this technique can be applied to functions which are suitable to approximation with piecewise polynomials.

bethandtownes commented 5 years ago

Hi I spent two days working on the rectangle area problem I mentioned to you. Unfortunately, I could not adapt your template to this problem and solve it in O(nlog(n)) time. I fathom the reason could be because, in this template, when doing the range apply method, it could not access the left and right child of the node. I tried to work around this by adding the child information in the Monoid structure; however, I don't this it will go in the right recursive order and thus the lchild and rchild info are not updated correctly. This indeed can be worked around by adding an array<pair<double,double>,101> of the sort. But this would make this problem O(n^2logn).

Here is an O(NlogN) solution(https://github.com/morris821028/UVa/blob/master/OnlineJudge/POJ/POJ%201151%20-%20Atlantis.cpp) using segment tree with recursive structure. Note that in his code, Node.clen is updated in sort of like tree DP manner.

Could you take a look at it and tell me your thoughts on this? I really would like to see this problem solved in O(NlogN) using your template...

kmyk commented 5 years ago

This O(N log N) solution is very interesting! Thank you! I didn't know that kind of segtrees, and this cannot be expressed on my templete.

I've formalized lazy-propagation segtrees with monoids and operator monoids. Here, I implicitly require the totality of functions. However, it seems that the solution you mentioned uses a partial operator monoid. Therefore we cannot be express this tree with my templete. This partiality-extension seems useful, but it will make operations not O(log n) if it is too partial.

Here is my code written based on this understanding, with modification of my templete of segtree: https://cn.vjudge.net/solution/16804839 Given code and my code work with almost same speed at N = 10^5, 10^6, so I think they does the same thing correctly.

bethandtownes commented 5 years ago

Thanks for generalizing it. Let me try to see if there is a way to get rid of the try catch block and maintain the structural integrity of the template.

By the looks of it, I think this new template is downward compatible to the original template?

bethandtownes commented 5 years ago

I think the try catch block can be replaced by "options". To my understanding, C++ supports this idea since C++17? (https://en.cppreference.com/w/cpp/utility/optional)

bethandtownes commented 5 years ago

forget about it.. I think given that this OJ can only support C11...The option concept is out of question...

kmyk commented 5 years ago

My new template is completely backward compatible for my old one if I rename the identifier of exception. However, compatibility may not be required, since this new template is useful but too generalized and I won't remove the old one. Also, as you said, we want to use optional and we cannot use it on some OJs.

I'm thinking about the condition that the complexity becomes O(log n). Apparently it may becomes O(n), so this question is important. For example, we can define the all cases undefined except identity one and leaf nodes. I think one of sufficient condition is that: undefined doesn't appears recursively. The problem you mentioned satisfies this condition. Please try to think about better sufficient conditions and/or necessary conditions.

bethandtownes commented 5 years ago

what does it mean by "change the name of the exception"?

kmyk commented 5 years ago

I've used exceptions defining on monoids like below and this breaks the compatibility. I mentioned this. But of course, this is not important.

            try {
                ...
            } catch (typename OperatorMonoid::domain_error e) {
                ...
            }
bethandtownes commented 5 years ago

how would you rename it and maintain the compatibility?

kmyk commented 5 years ago

I'll define it on other namespace or just use std::domain_error .

bethandtownes commented 5 years ago

oh..I see what you are saying. thx for clearing it up. Are you doing grad school now? I see you work late very often..(Isn't it 2am now in Japan?..)

By the way, I am a grad student in statistics. I studied math in college. It indeed would be a pleasure to get to know you a bit:)

kmyk commented 5 years ago

Yes, I am in a graduate school and it is 2 am in Japan (and it becomes 3 am just now).

Also, I’m glad that you’re saying that. Thank you.

bethandtownes commented 5 years ago

I think it won't be too hard to come up with a probabilistic running time of the partiality-extension sgtree by assuming some kind of uniform distribution on query and concats. Probably we need to define a few new things, though. In this case, my conjecture is that the cutoff between nlogn and n^2 would depend on the distribution and some other params. The theorem, ideally, would state that with Prob 1-\exp^{-{function on degree?distribution param? of the sorts}}, the algo runs nlogn.

As for non-probabilistic amortized run time of the range update, I am not an expert in this area.

I am quite busy with my phd research this month. But I will keep on thinking about this as a research problem in my spare time. Perhaps if the result looks promising, would you be interested to write a short paper about it sometime early next year?

Although it is not a significant thing, but it does sound interesting and the result, if could be given, does look cute.

bethandtownes commented 5 years ago

In this specific problem, it can be easily shown that f[i] won't be pushed recursively due to the nature of the problem. But I doubt this can be easily shown or even exists in other problems. So it would be a good idea to relax the condition a little. Maybe we could cook up something like the "distribution of recursive depth"(or more generally, a combinatorial representation (possibly with a probabilistic component, hence the name "distribution") of the recursive behavior) of the undefined behavior and use this as bridge to calculate the cut-off between qlogn and qn. I think the tool needed here is some analytic combinatorics.

kmyk commented 5 years ago

Thank you for your idea. And I'm sorry for my late reply. I've thought about only non-probabilistic one, so your idea is interesting for me. Also, I'm interested at writing a short paper if we can find some nice conditions.

In terms of relaxing, we may be able to use the maximal depth of recursion. Let d be the maximum, time complexities becomes roughly O(n log n * 2^d). Also we may replace this maximum d with some average d'. But this is not so interesting.

I think we need to find problems which we can solve with partiality-ext segtree. The problem of counting values with >= 1 is, as you said, too simple to research about the condition of O(n log n). It seems there are some problems whose intended solution uses this technique implicitly. Or, they may appears as imos-method, since many segtrees can be replaced with imos-method when there is no update-queries after concat-queries.

I cannot use much time for this problem. However, I will keep thinking about this too. For now, I'll try to find and solve advanced tasks of competitive programming using segtrees and imos.

bethandtownes commented 5 years ago

Hi It's been a while and I hope things are going well! Just want to point out that this link for your vjudge code piece is no longer working. I am hoping you can put it back on somewhere so that it can serve as an illustrative example:)

bethandtownes commented 5 years ago

Could you share it here? I could not find my copy of it nor could I reproduce it...Really appreciate it!

kmyk commented 4 years ago

Sorry for late to reply :bow: Here is the code.

I can see my code with logging in. I don't know why the code on vjudge have been invisible.

#include <bits/stdc++.h>
#define REP(i, n) for (int i = 0; (i) < (int)(n); ++ (i))
#define REP_R(i, n) for (int i = int(n) - 1; (i) >= 0; -- (i))
#define ALL(x) begin(x), end(x)
using namespace std;

template <class Monoid, class OperatorMonoid>
struct lazy_propagation_segment_tree { // on monoids
    static_assert (is_same<typename Monoid::underlying_type, typename OperatorMonoid::target_type>::value, "");
    typedef typename Monoid::underlying_type underlying_type;
    typedef typename OperatorMonoid::underlying_type operator_type;
    const Monoid mon;
    const OperatorMonoid op;
    int n;
    vector<underlying_type> a;
    vector<operator_type> f;
    lazy_propagation_segment_tree() = default;
    lazy_propagation_segment_tree(int a_n, underlying_type initial_value = Monoid().unit(), Monoid const & a_mon = Monoid(), OperatorMonoid const & a_op = OperatorMonoid())
            : mon(a_mon), op(a_op) {
        n = 1; while (n <= a_n) n *= 2;
        a.resize(2 * n - 1, mon.unit());
        fill(a.begin() + (n - 1), a.begin() + ((n - 1) + a_n), initial_value); // set initial values
        REP_R (i, n - 1) a[i] = mon.append(a[2 * i + 1], a[2 * i + 2]); // propagate initial values
        f.resize(max(0, (2 * n - 1) - n), op.identity());
    }
    void point_set(int i, underlying_type z) {
        assert (0 <= i and i < n);
        point_set(0, 0, n, i, z);
    }
    void point_set(int i, int il, int ir, int j, underlying_type z) {
        if (i == n + j - 1) { // 0-based
            a[i] = z;
        } else if (ir <= j or j+1 <= il) {
            // nop
        } else {
            flush(i, il, ir, false);
            point_set(2 * i + 1, il, (il + ir) / 2, j, z);
            point_set(2 * i + 2, (il + ir) / 2, ir, j, z);
            a[i] = mon.append(a[2 * i + 1], a[2 * i + 2]);
        }
    }
    void range_apply(int l, int r, operator_type z) {
        assert (0 <= l and l <= r and r <= n);
        if (z == op.identity()) return;
        range_apply(0, 0, n, l, r, z);
    }
    void range_apply(int i, int il, int ir, int l, int r, operator_type z) {
        if (l <= il and ir <= r) { // 0-based
            if (i < f.size()) f[i] = op.compose(z, f[i]);
            try {
                a[i] = op.apply(z, a[i]);
            } catch (typename OperatorMonoid::domain_error e) {
                if (i < f.size()) {
                    flush(i, il, ir, false);
                    a[i] = mon.append(a[2 * i + 1], a[2 * i + 2]);
                } else {
                    assert (false);
                }
            }
        } else if (ir <= l or r <= il) {
            // nop
        } else {
            flush(i, il, ir, false);
            range_apply(2 * i + 1, il, (il + ir) / 2, l, r, z);
            range_apply(2 * i + 2, (il + ir) / 2, ir, l, r, z);
            a[i] = mon.append(a[2 * i + 1], a[2 * i + 2]);
        }
    }
    underlying_type range_concat(int l, int r) {
        assert (0 <= l and l <= r and r <= n);
        return range_concat(0, 0, n, l, r);
    }
    underlying_type range_concat(int i, int il, int ir, int l, int r) {
        if (l <= il and ir <= r) { // 0-based
            return a[i];
        } else if (ir <= l or r <= il) {
            return mon.unit();
        } else {
            try {
                return op.apply(f[i], mon.append(
                    range_concat(2 * i + 1, il, (il + ir) / 2, l, r),
                    range_concat(2 * i + 2, (il + ir) / 2, ir, l, r)));
            } catch (typename OperatorMonoid::domain_error e) {
                flush(i, il, ir, true);
                return mon.append(
                    range_concat(2 * i + 1, il, (il + ir) / 2, l, r),
                    range_concat(2 * i + 2, (il + ir) / 2, ir, l, r));
            }
        }
    }
private:
    void flush(int i, int il, int ir, bool pred) {
        if (f[i] == op.identity()) return;
        range_apply(2 * i + 1, il, (il + ir) / 2, 0, n, f[i]);
        range_apply(2 * i + 2, (il + ir) / 2, ir, 0, n, f[i]);
        f[i] = op.identity();
        if (pred) a[i] = mon.append(a[2 * i + 1], a[2 * i + 2]);
    }
};

template <typename T>
struct count_monoid_with_range {
    typedef struct {
        T l, r;  // of range
        bool is_atomic;
        int multiplicity;
        T covered;
    } underlying_type;
    underlying_type unit() const {
        return (underlying_type) { -1, -1, false, -1, -1 };
    }
    underlying_type append(underlying_type a, underlying_type b) const {
        if (a.multiplicity == -1) return b;
        if (b.multiplicity == -1) return a;
        underlying_type c;
        c.l = a.l;
        c.r = b.r;
        c.is_atomic = false;
        c.multiplicity = min(a.multiplicity, b.multiplicity);
        c.covered = a.covered + b.covered;
        return c;
    }
};

template <typename T>
struct increment_operator_monoid_with_range {
    typedef int underlying_type;
    typedef typename count_monoid_with_range<T>::underlying_type target_type;
    typedef struct {} domain_error;
    underlying_type identity() const {
        return 0;
    }
    target_type apply(underlying_type a, target_type b) const {
        if (a == 0) return b;
        b.multiplicity += a;
        if (b.multiplicity > 0) {
            b.covered = b.r - b.l;
        } else if (b.is_atomic) {
            b.covered = 0;
        } else {
            throw (domain_error) {};
        }
        return b;
    }
    underlying_type compose(underlying_type a, underlying_type b) const {
        return a + b;
    }
};

double solve(int n, vector<double> const & x1, vector<double> const & y1, vector<double> const & x2, vector<double> const & y2) {
    // coordinate compression
    vector<double> xs;
    xs.insert(xs.end(), ALL(x1));
    xs.insert(xs.end(), ALL(x2));
    sort(ALL(xs));
    xs.erase(unique(ALL(xs)), xs.end());
    auto compress = [&](double x) { return lower_bound(ALL(xs), x) - xs.begin(); };

    // prepare segtree
    lazy_propagation_segment_tree<count_monoid_with_range<double>, increment_operator_monoid_with_range<double> > segtree(xs.size() - 1);
    REP (i, xs.size() - 1) {
        typename count_monoid_with_range<double>::underlying_type a;
        a.l = xs[i];
        a.r = xs[i + 1];
        a.is_atomic = true;
        a.multiplicity = 0;
        a.covered = 0;
        segtree.point_set(i, a);
    }

    // prepare events
    vector<tuple<double, int, int, int> > events;
    REP (i, n) {
        events.emplace_back(y1[i], compress(x1[i]), compress(x2[i]), +1);
        events.emplace_back(y2[i], compress(x1[i]), compress(x2[i]), -1);
    }
    sort(ALL(events));

    // use sweepline
    double acc = 0;
    double last_y = 0;
    for (auto event : events) {
        double y; int l, r, delta; tie(y, l, r, delta) = event;
        acc += (y - last_y) * segtree.range_concat(0, xs.size()).covered;
        segtree.range_apply(l, r, delta);
        last_y = y;
    }

    return acc;
}

int main() {
    for (int testcase = 1; ; ++ testcase) {
        int n; cin >> n;
        if (n == 0) break;
        vector<double> x1(n), y1(n), x2(n), y2(n);
        REP (i, n) {
            cin >> x1[i] >> y1[i] >> x2[i] >> y2[i];
        }
        cout << "Test case #" << testcase << endl;
        cout << "Total explored area: " << fixed << setprecision(2) << solve(n, x1, y1, x2, y2) << endl;
        cout << endl;
    }
    return 0;
}
kmyk commented 4 years ago

By the way, I still haven't found any other usage of segment trees similar to this one. This tree may be understood by separately from trees with monoids. e.g. also Li-Chao tree is understood as unrelated to monoids.

kmyk commented 4 years ago

I recently found a way to solve this problem https://vjudge.net/problem/HDU-1542 with a usual lazy-propagation segtree.

This problem is counting values with >= 1 from a sequence of non-negative integers. And this problem is equivalent to counting values = 0. It's possible with segtree using the non-negativity of the sequence, as follows:

  1. find a minimum value x in the sequence and the number k of occurrences of x
  2. the number of zeros is k if x is zero; else the number is 0
bethandtownes commented 4 years ago

Great! Would you like to share the code as a demonstrative example?

kmyk commented 4 years ago

Please read this https://atcoder.jp/contests/jsc2019-final/submissions/8691511 (around L279).

bethandtownes commented 4 years ago

this is beautiful code. I just want to say thank you since I learnt quite a bit from your way of programming. I am a machine learning researcher and currently developing new model package in Julia. It is mostly mathematics in its pure form but written with care to guarantee scalability. Many of the package design idea originated from your library (it is surprising, I know). I hope if you don't mind I can add your name as a special mention to the to be published paper soon.

Japan is a such a nice place. I have been there multiple times. Clean air, great people. I hope one day I can visit there again and perhaps we can get a cup of coffee:)

Jason

kmyk commented 4 years ago

I'm very glad that my code help you. That is worth opening this library to the public :smile: I cannot imagine how your design originates from my library, but it's my pleasure to mention my name.

Also, please tell me when you come Japan. Let's take coffee :coffee: