# load the package badic (you can install it with "sage -pip install badic [--user]")
from badic import *

# define a function that computes the interior of a beta-adic set,
# assuming that the beta-adic set with
# any word allowed on the same alphabet has zero in its interior
# !! the topology used here is a special one !!
def interior(self, verb=False):
    r"""
    We assume that self.b is a Pisot number.
    Compute a BetaAdicSet describing the interior for the topology for which open sets are
    sets of points of that are in the BetaAdicSet with same beta but with a language that recognize every words over the alphabet of self.a and that projects to an open set of P, for the natural projection on the contracting space (which is a product of copies of R, C and p-adic fields).
    """
    def S(a, verb=False):
        F = []
        for e in a.states:
            ok = True
            for l in range(len(a.alphabet)):
                if a.succ(e, l) != e:
                    ok = False
                    break
                if ok:
                    F.append(e)
                    a2 = a.copy()
                    a2.set_final_states(F)
                    if verb:
                        print("F =",F)
                    return a2
    arel = self.relations_automaton(t=0,couples=True)
    if verb:
        print("arel =",arel)
        #arel.plot()
    ap = dag.AnyWord(self.a.alphabet).product(self.a.concat_zero_star())
    ai = ap.intersection(arel)
    if verb:
        print("ai =",ai)
    aip = ai.proji(0)
    if verb:
        print("aip =", aip)
    aip.zero_complete_op()
    af = S(aip.minimize())
    af.zero_complete_op()
    if verb:
        print("af =",af)
    af2 = af.minimize().intersection(self.a)
    af2 = af2.prune()
    if verb:
        print("af2 =",af2)
        print(af2.equals_langages(self.a))
    return BetaAdicSet(self.b, af2) 

# take the beta-adic set of the Tribonnacci Rauzy fractal
s = WordMorphism('a->ab,b->ac,c->a')
m = DumontThomas(s).mirror()
m.plot()

# compute the interior
mi = interior(m)
mi.is_equal_to(m) # the interior is the whole beta-adic set

True

# take the substitution associated to the smallest Pisot number
s = WordMorphism('a->b,b->c,c->ab')
m = DumontThomas(s).mirror()
m.plot()

# compute the interior
mi = interior(m)
mi.plot()

# allow the user to zoom in this fractal in real-time
l = mi.draw_zoom()

# plot where the user has zoomed
mi.plot(ajust=False)

# take the boundary
mc = mi.complementary(m)
mc.plot()

# compute its Hausdorff dimension
mc.critical_exponent_free()

log(y)/log(1.150963925257758?) where y is the max root of x^13 - x^12 - x^10 + x^9 - 2*x^4 + x^3 - 1, and 1.150963925257758? is root of x^6 - x^2 - 1.

1.9464346032652802

# flipped Tribonnacci
s = WordMorphism('a->ab,b->ca,c->a')
m = DumontThomas(s).mirror()
m.plot()

m.a.plot()

mi = interior(m)
print(mi)
mi.plot()

b-adic set with b root of x^3 - x^2 - x - 1, and an automaton of 79 states and 3 letters

mi.a.plot()

mi.a

DetAutomaton with 79 states and an alphabet of 3 letters

# boundary
mi.complementary(m).plot()

# example with overlaps
s = WordMorphism({1:[1,2], 2:[2,3], 3:[4], 4:[5], 5:[6], 6:[7], 7:[8], 8:[9], 9:[10], 10:[1]})
# choose the direction of projection
pi = x^3-x-1
b = pi.roots(ring=QQbar)[1][0]
print(b)
K.<b> = NumberField(pi, embedding=b)
m = DumontThomas(s, b=b, verb=True).mirror()
m.plot()

-0.6623589786223730? - 0.5622795120623013?*I
(1, b - 1, b^2 - 2*b + 1, -2*b^2 + 2*b + 1, 2*b^2 - b - 2, -b^2 + 2, b - 1, b^2 - b, -b^2 + b + 1, b^2 - 1)
initial state 1
[(1, 1, 0), (1, 2, 1), (2, 2, 0), (2, 3, b - 1), (3, 4, 0), (4, 5, 0), (5, 6, 0), (6, 7, 0), (7, 8, 0), (8, 9, 0), (9, 10, 0), (10, 1, 0)]
DetAutomaton with 10 states and an alphabet of 3 letters

mi = interior(m)
mi

b-adic set with b root of x^3 - x - 1, and an automaton of 505 states and 3 letters

mi.plot()

s.rauzy_fractal_plot()