r"""
The associated paper is given at

    M. Festa, A. Sommariva,
    Computing Almost Minimal Formulas on the Square,
    Journal of Computational and Applied Mathematics, Vol. 236, pp. 4296--4302, 2012,
    https://doi.org/10.1016/j.cam.2012.05.021

The quadrature rules extracted here can be found
    https://www.math.unipd.it/~alvise/sets.html

These rules are licensed under the GPL and are not included by default.
"""

import re

import numpy as np
from pytools import download_from_web_if_not_present


_URL = "https://web.archive.org/web/20201029182232/http://www.math.unipd.it/~alvise/POINTSETS/set_amr_square.m"  # noqa

_PYTHON_TEMPLATE = """# GENERATED by modepy/contrib/weights-from-festa-sommariva.py
# DO NOT EDIT

import numpy as np
from modepy.quadrature import Quadrature, QuadratureRuleUnavailable


class FestaSommarivaQuadrature(Quadrature):
    def __init__(self, order):
        try:
            rule = square_data[order]
        except KeyError:
            raise QuadratureRuleUnavailable

        super().__init__(rule["points"], rule["weights"])
        self.exact_to = rule["quad_degree"]


def process_rules(rules):
    result = {}
    for order, rule in rules.items():
        result[order] = {
            k: np.array(v) if isinstance(v, list) else v
            for k, v in rule.items()
            }

    return result


square_data = process_rules(%s)"""


def generate_festa_sommariva_quadrature_rules(outfile):
    filename = f"set_amr_square.m"
    download_from_web_if_not_present(
            url=f"{_URL}",
            local_name=filename)

    with open(filename, "r") as fd:
        mfile = fd.read()

    # rules are hardcoded as part of a MATLAB file
    # the are extracted by looking for the pattern
    #
    #   xw=[
    #       <x> <y> <weight>
    #       ...
    #   ];
    #
    # The degree of each quadrature rule is written in a line above the
    # nodes and weights in the pattern
    #
    #   DEGREE:  <degree>

    _re_degree = re.compile(r"DEGREE:\s+(\d+)")
    _re_xw = re.compile(r"xw\s?=\s?\[(.+?)\];", re.DOTALL)

    # NOTE: some degrees have multiple quadrature rules with a repeated
    # header, so we just take the unique ones here
    degrees = np.unique(
            np.fromiter(_re_degree.findall(mfile), dtype=np.int)
            )

    rules = {}
    for imatch, match in enumerate(_re_xw.findall(mfile)):
        d = degrees[imatch]
        assert d == imatch + 1, (d, imatch + 1)

        xw = np.fromstring(match.strip(), dtype=np.float64, sep=" ").reshape(-1, 3)
        rules[d] = {
                "quad_degree": d,
                "points": xw[:, :-1].T.tolist(),
                "weights": xw[:, -1].tolist(),
                }

        assert abs(np.sum(xw[:, -1]) - 4.0) < 1.0e-12
        print("degree {} points {} / {}".format(
            d, xw.shape[0], (d + 1) * (d + 2) // 2
            ))

    from pprint import pformat
    txt = (_PYTHON_TEMPLATE % pformat(rules, width=80)).replace('"', "")

    if outfile:
        with open(outfile, "w") as fd:
            fd.write(txt)
    else:
        print(txt)


if __name__ == "__main__":
    import argparse

    parser = argparse.ArgumentParser()
    parser.add_argument("outfile", nargs="?", default="")
    args = parser.parse_args()

    generate_festa_sommariva_quadrature_rules(outfile=args.outfile)
