####THIS PROGRAM RUNS IN SAGEMATH 8.1
####IT WORKS FOR 2 BRANCHES WITHOUT OVERLAP
####IT REQUIRES NEW VERSION FOR 3 OR MORE BRANCHES
####IT REQUIRES NEW VERSION FOR OVERLAPS

################################ definitions of all functions

import cython
import itertools
from mpmath import *

############ Defintion of partitions n_1+n_2+...+n_l=n

def partitions(k, N):
    # return all k-partitions of N
    tab = []
    if k == 1:
        return [N]
    for i in range(1,N):
        resto = partitions(k-1, N-i)
        for r in resto:
            L=[i]
            L.append(r)
            tab.append(L)
    for i in range(len(tab)):
        tab[i] = copy(list(flatten(tab[i]))) 
    return tab
     
def Tabla(N):
    R=[]
    for k in range(2,N+1):
        R.append([partitions(k,i) for i in range(k,N+1)])
    return R
     
############# Definition of indentity function

def identidad(x):
    return x

############# Definition of Dictionary

def Dictionary(m,N):
	aux=str(0)
	for i in range(1,m):
		aux=aux+str(i)
	K=range(N)
	for i in range(N):
		K[i]=["".join(seq) for seq in itertools.product(aux,repeat=i+1)]
	return K

############# Definition of f(T,word) : f(T,'000001')

def f(T,word):
	#word='000'
	m=len(word)
	g=identidad
	for i in range(m):
		g=compose(T[int(word[m-1-i])],g)
	return g

############# Definition of fixed points Tablefixpt(T,K)

def Tablefixpt(T,K):
	N=len(K)
	Lfix=[[] for i in range(N)]
	cont=0
	for words in K:
		for i in range(len(words)):
			def F(x):
				return f(T,words[i])(x)
			def FF(x):
				return F(x)-x
			x=findroot(FF,0)
			Lfix[cont].append(x)
		cont=cont+1
	return Lfix

############# Find the fixpt associated to a word  Findfixpt(Tfixpt,'000')

def Findfixpt(Tfixpt,word):
    index=len(word)
    coordinate=sum([int(word[index-1-i])*2^(i) for i in range(index)])
    return Tfixpt[index-1][coordinate]

############# Definition of shift cycling

def Shift_sec(word,n):
	aux=word[n:len(word)]+word[0:n]
	return aux

############# Definition of p_1*p_2*...*p_N P([1/3,2/3],'0101')

def P(p,word):
	return product([p[int(word[i])] for i in range(len(word))])

############# Definition of shift delete of word Shift_del('01110',2)

def Shift_del(word,n):
	aux=word[n:len(word)]
	return aux

############# Definition of S_ng(z_{i_1},z_{i_2},...,z_{i_n})

def S(Tfixpt,word,g):
	return fsum([g(Findfixpt(Tfixpt,Shift_sec(word,_))) for _ in range(len(word))])

############# Definition of t_m

def t_fun(K,Tfixpt,p,m):
	return fsum([P(p,word)/(1-DT(Tfixpt,word)) for word in K[m]])

############# Definition of tau_m

def tau_fun(K,Tfixpt,p,m,g):
	return fsum([P(p,word)*S(Tfixpt,word,g)/(1-DT(Tfixpt,word)) for word in K[m]])

############# Definition of Integration approximation

def Integration(T,p,g,N):
    R       =   Tabla(N)
    K       =   Dictionary(len(T),N)
    Tfixpt  =   Tablefixpt(T,K)
    t       =   [t_fun(K,Tfixpt,p,m) for m in range(N)]
    tau     =   [tau_fun(K,Tfixpt,p,m,g) for m in range(N)]
    
    alpha   =   [0 for i in range(N)]
    alpha[0]=   (-1)*tau[0]
    a       =   [0 for i in range(N)]
    a[0]    =   (-1)*t[0]
        
    for n in range(2,N+1):
        sum_alpha   =   (-1)*tau[n-1]/n
        sum_a       =   (-1)*t[n-1]/n
        for l in range(2,n+1):
            sum_alpha_l =0
            sum_a_l     =0
            for n_vec in R[l-2][n-l]:
                product =  fprod([t[n_vec[i]-1]/n_vec[i] for i in range(l)])
                for j in range(l):
                    sum_alpha_l = sum_alpha_l+tau[n_vec[j]-1]/(t[n_vec[j]-1])*product
                sum_a_l=sum_a_l+product
            sum_alpha=sum_alpha+(-1)^l*sum_alpha_l/(l*1).factorial()
            sum_a=sum_a+(-1)^l*sum_a_l/(l*1).factorial()
        alpha[n-1]=copy(sum_alpha)
        a[n-1]=copy(sum_a)
    an=[(i+1)*a[i] for i in range(N)]
    res=sum(alpha)/sum(an)
    return res


################# main part of program

################ EXAMPLE IN PAPER WITH NATALIA

################# Regulate accuracy

cm         = sage.structure.element.get_coercion_model
mp.prec    = 1000
mp.dps     = 300 
mp.pretty  = True

################# definition of branches

def f0(x):
	return x/3

def f1(x):
	return x/3+2/3

T= [f0,f1]

################# definition of derivatives DT_{i_1,i_2,...,i_n}(z_{i_1},z_{i_2},...,z_{i_n}) 

def DT(Tfixpt,word):
    
    n=len(word)
    prod=1
    
    for i in range(1,n+1):
        evalx=Findfixpt(Tfixpt,Shift_sec(word,i))
    
        if word[i-1]=='0':
            def F0(x):
                return x/3
            prod=fprod([prod,diff(F0, evalx, addprec=300)])
        
        if word[i-1]=='1':
            def F1(x):
                return x/3+2/3
            prod=fprod([prod,diff(F1, evalx, addprec=300)])
            
    return prod
    
################# definition of probability vector

p=[1/3,2/3]

################# Running times and results for different approximations


def g0(x):
    return 1    

def g1(x):
    return x    

def g2(x):
    return x^2
    
def g3(x):
    return x^3
    
def g4(x):
    return x^4
    
def g5(x):
    return x^5
    
def g6(x):
    return x^6
    
def g7(x):
    return x^7
    
def g8(x):
    return x^8
    
def g9(x):
    return x^9
    
def g10(x):
    return x^10
                                                              
%time Integration(T,p,g0,14)
%time Integration(T,p,g1,14)
%time Integration(T,p,g2,14)
%time Integration(T,p,g3,14)
%time Integration(T,p,g4,14)
%time Integration(T,p,g5,14)
%time Integration(T,p,g6,14)
%time Integration(T,p,g7,14)
%time Integration(T,p,g8,14)
%time Integration(T,p,g9,14)
%time Integration(T,p,g10,14)

%time Integration(T,p,g4,1)
%time Integration(T,p,g4,2)
%time Integration(T,p,g4,3)
%time Integration(T,p,g4,4)
%time Integration(T,p,g4,5)
%time Integration(T,p,g4,6)
%time Integration(T,p,g4,7)
%time Integration(T,p,g4,8)
%time Integration(T,p,g4,9)
%time Integration(T,p,g4,10)
%time Integration(T,p,g4,11)
%time Integration(T,p,g4,12)
%time Integration(T,p,g4,13)
%time Integration(T,p,g4,14)
%time Integration(T,p,g4,15)



################# EXAMPLE IN PAPER WITH NATALIA

################# Regulate accuracy

cm         = sage.structure.element.get_coercion_model()
mp.dps     = 1000; 
mp.pretty  = True

################# definition of branches

def f0(x):
	return 1/(x+2)

def f1(x):
	return 1/(x+4)

T= [f0,f1]

################# definition of derivatives DT_{i_1,i_2,...,i_n}(z_{i_1},z_{i_2},...,z_{i_n}) 

def DT(Tfixpt,word):
    
    n=len(word)
    prod=1
    
    for i in range(1,n+1):
        evalx=Findfixpt(Tfixpt,Shift_sec(word,i))
    
        if word[i-1]=='0':
            def F0(x):
                return 1/(x+2)
            prod=fprod([prod,diff(F0, evalx, addprec=300)])
        
        if word[i-1]=='1':    
            def F1(x):
                return 1/(x+4)
            prod=fprod([prod,diff(F1, evalx, addprec=300)])
            
    return prod
    
################# definition of probability vector

p=[1/2,1/2]


%time Integration(T,p,g0,13)
%time Integration(T,p,g0,14)

%time Integration(T,p,g1,13)
%time Integration(T,p,g1,14)

%time Integration(T,p,g2,13)
%time Integration(T,p,g2,14)

%time Integration(T,p,g3,13)
%time Integration(T,p,g3,14)

%time Integration(T,p,g4,13)
%time Integration(T,p,g4,14)

%time Integration(T,p,g5,13)
%time Integration(T,p,g5,14)

%time Integration(T,p,g6,13)
%time Integration(T,p,g6,14)

%time Integration(T,p,g7,13)
%time Integration(T,p,g7,14)

%time Integration(T,p,g8,13)
%time Integration(T,p,g8,14)

%time Integration(T,p,g9,13)
%time Integration(T,p,g9,14)

%time Integration(T,p,g10,13)
%time Integration(T,p,g10,14)

################# RESULTS OBTAINED

CPU times: user 1min 14s, sys: 713 ms, total: 1min 14s
Wall time: 1min 16s
0.9999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999983
CPU times: user 2min 42s, sys: 1.34 s, total: 2min 43s
Wall time: 2min 46s
0.9999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999998
CPU times: user 1min 13s, sys: 440 ms, total: 1min 14s
Wall time: 1min 15s
0.3304697175264855340801384795184068289815344294101275920335337740232421085586168880132686762410032041
CPU times: user 2min 39s, sys: 877 ms, total: 2min 40s
Wall time: 2min 42s
0.3304697175264855340801384795184068289815344294101275920335337740232421086739865571343423066960105978
CPU times: user 1min 18s, sys: 885 ms, total: 1min 19s
Wall time: 1min 22s
0.1192803776960544798961200249581823359145663180489550186633549589883397537243274811414567108729214671
CPU times: user 2min 49s, sys: 1.81 s, total: 2min 50s
Wall time: 2min 55s
0.1192803776960544798961200249581823359145663180489550186633549589883397537981023761468207871249330553
CPU times: user 1min 13s, sys: 582 ms, total: 1min 14s
Wall time: 1min 16s
0.04612084018579153108812762740892741033130217767374947447203645639401348913411599338286105062750772329
CPU times: user 2min 58s, sys: 1.63 s, total: 3min
Wall time: 3min 10s
0.04612084018579153108812762740892741033130217767374947447203645639401348916954618626747796849682878365
CPU times: user 1min 7s, sys: 309 ms, total: 1min 8s
Wall time: 1min 8s
0.01869566793194042885495853137233781484714063870527561442396421977030979100260303114189385694100625803
CPU times: user 2min 33s, sys: 820 ms, total: 2min 33s
Wall time: 2min 34s
0.01869566793194042885495853137233781484714063870527561442396421977030979101793833918499302783599089063
CPU times: user 1min 12s, sys: 350 ms, total: 1min 13s
Wall time: 1min 13s
0.007807847977063560924555378015913035121359128312247536049064442681158589353249126190734479739667474989
CPU times: user 2min 37s, sys: 752 ms, total: 2min 38s
Wall time: 2min 39s
0.007807847977063560924555378015913035121359128312247536049064442681158589359566459581046167677426733819
CPU times: user 1min 5s, sys: 188 ms, total: 1min 6s
Wall time: 1min 6s
0.003320110573231903768685936564676714748591723927498902758538761678623163985979156848132013950560249321
CPU times: user 2min 21s, sys: 81 ms, total: 2min 21s
Wall time: 2min 22s
0.003320110573231903768685936564676714748591723927498902758538761678623163988507361170715147522992638201
CPU times: user 1min 4s, sys: 26.3 ms, total: 1min 4s
Wall time: 1min 4s
0.001427182118378503652418281093362766634002528335724754419216502382209445217552982761624051598355717609
CPU times: user 2min 22s, sys: 90.1 ms, total: 2min 22s
Wall time: 2min 22s
0.001427182118378503652418281093362766634002528335724754419216502382209445218544509065321157973124890233
CPU times: user 1min 4s, sys: 35.2 ms, total: 1min 4s
Wall time: 1min 5s
0.0006175983077855314124072271759192922823328720962838360688107425294482557463499976211616025745420542472
CPU times: user 2min 26s, sys: 258 ms, total: 2min 26s
Wall time: 7min 18s
0.000617598307785531412407227175919292282332872096283836068810742529448255746732606764133290216605824148
CPU times: user 1min 10s, sys: 159 ms, total: 1min 10s
Wall time: 1min 10s
0.0002684218626959226510757272950170889068103411627636333206312907565016840913860178120157051611828420354
CPU times: user 2min 36s, sys: 426 ms, total: 2min 36s
Wall time: 42min 10s
0.0002684218626959226510757272950170889068103411627636333206312907565016840915315187509076259810635326489
CPU times: user 1min 18s, sys: 395 ms, total: 1min 18s
Wall time: 1min 19s
0.0001170171983606956288423161480535694711087681239194925952930086491322776677776781347154798576946949146
CPU times: user 2min 38s, sys: 611 ms, total: 2min 39s
Wall time: 14min 42s
0.0001170171983606956288423161480535694711087681239194925952930086491322776678321930948023414732521038212





CPU times: user 1min 14s, sys: 713 ms, total: 1min 14s
Wall time: 1min 16s
0.9999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999983
0.9999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999998
CPU times: user 1min 13s, sys: 440 ms, total: 1min 14s
Wall time: 1min 15s
0.330469717526485534080138479518406828981534429410127592033533774023242108 5586168880132686762410032041
0.330469717526485534080138479518406828981534429410127592033533774023242108 6739865571343423066960105978
CPU times: user 1min 18s, sys: 885 ms, total: 1min 19s
Wall time: 1min 22s
0.1192803776960544798961200249581823359145663180489550186633549589883397537 243274811414567108729214671
0.1192803776960544798961200249581823359145663180489550186633549589883397537 981023761468207871249330553
CPU times: user 1min 13s, sys: 582 ms, total: 1min 14s
Wall time: 1min 16s
0.0461208401857915310881276274089274103313021776737494744720364563940134891 3411599338286105062750772329
0.0461208401857915310881276274089274103313021776737494744720364563940134891 6954618626747796849682878365
CPU times: user 1min 7s, sys: 309 ms, total: 1min 8s
Wall time: 1min 8s
0.0186956679319404288549585313723378148471406387052756144239642197703097910 0260303114189385694100625803
0.0186956679319404288549585313723378148471406387052756144239642197703097910 1793833918499302783599089063
CPU times: user 1min 12s, sys: 350 ms, total: 1min 13s
Wall time: 1min 13s
0.0078078479770635609245553780159130351213591283122475360490644426811585893 53249126190734479739667474989
0.0078078479770635609245553780159130351213591283122475360490644426811585893 59566459581046167677426733819
CPU times: user 1min 5s, sys: 188 ms, total: 1min 6s
Wall time: 1min 6s
0.0033201105732319037686859365646767147485917239274989027585387616786231639 85979156848132013950560249321
0.0033201105732319037686859365646767147485917239274989027585387616786231639 88507361170715147522992638201
CPU times: user 1min 4s, sys: 26.3 ms, total: 1min 4s
Wall time: 1min 4s
0.00142718211837850365241828109336276663400252833572475441921650238220944521 7552982761624051598355717609
0.00142718211837850365241828109336276663400252833572475441921650238220944521 8544509065321157973124890233
CPU times: user 1min 4s, sys: 35.2 ms, total: 1min 4s
Wall time: 1min 5s
0.000617598307785531412407227175919292282332872096283836068810742529448255746 3499976211616025745420542472
0.000617598307785531412407227175919292282332872096283836068810742529448255746 732606764133290216605824148
CPU times: user 1min 10s, sys: 159 ms, total: 1min 10s
Wall time: 1min 10s
0.000268421862695922651075727295017088906810341162763633320631290756501684091 3860178120157051611828420354
0.000268421862695922651075727295017088906810341162763633320631290756501684091 5315187509076259810635326489
CPU times: user 1min 18s, sys: 395 ms, total: 1min 18s
Wall time: 1min 19s
0.000117017198360695628842316148053569471108768123919492595293008649132277667 7776781347154798576946949146
0.000117017198360695628842316148053569471108768123919492595293008649132277667 8321930948023414732521038212