hfunction[N_,alpha_,beta_,gamma_,x1_,whichc1_]:=Module[{i,b,kc,Ic,c,ci,rorl,a,intx,intervalx,intervalxright,l,j,tao,taoright,x,taoc,out,graph,cumbeta,taon,k,s,isci,positionofx,positionofx2,whichc,whichc2,p,q,cont,S,S2,outtao,outbeta,v,eigen,D,is1eigen,ii,jj,cumbetaci,qci,chiofbeta,outchi,chi,betachi,kk,taonchi,sumend,dsubzero,h},
(*(Subscript[alpha, i]) represents the length of the ith branch (projection).
(Subscript[beta, i])  represents the slope of the ith branch. 
(Subscript[gamma, i]) represents the heigth of the lower endpoint of the ith branch. If (Subscript[gamma, i])=0 that branch is called greedy and if (Subscript[gamma, i])+(Subscript[alpha, i])=1 is called lazy.
N number of branches. 
x1 is the variable x of the function h(x) and whichc1 is used in case one wants to find out taon or cumbeta 
of a specific point.*)

whichc=whichc1;

If[Length[alpha]==Length[beta]&&Length[beta]==Length[gamma]&&Length[gamma]==N,
b=Table[0, {N+1}];

(*The endpoints of the DOMAIN of branches ((Subscript[b, i])). It always starts with 0*)
For[i=1,i≤N,i++,
b[[i+1]]=b[[i]]+alpha[[i]]/Abs[beta[[i]]]
];

b[[N+1]]=1;(*and ends in 1*)

(*Number of (Subscript[c, i]) points Kc*)
kc=0;
For[i=1,i≤N,i++,
If[alpha[[i]]<1,kc++];
If[gamma[[i]]>0 && (alpha[[i]]+gamma[[i]]<1),kc++];
];

(*rorl means right or left, this variable contains 1's (if (Subscript[c, i]) is in (Subscript[U, l])) or 0's (if (Subscript[c, i]) is in (Subscript[U, r]))*)
rorl=Table[,{kc}];
ci=1; (*This is the index of the vector C containing the (Subscript[c, i]). This variable is used when determining the characteristics of (Subscript[c, i])*)
c=Table[0,{kc}];
Ic=Table[0,{kc}]; (*tells us which interval is associated with (Subscript[c, i])*)

(*Finding the values of (Subscript[c, i]) and its characteristics*)
For[i=1,i≤N,i++, 
If[beta[[i]] > 0,(*Increasing branches*)
  If[alpha[[i]] < 1 && gamma[[i]] + alpha[[i]] == 1, (*lazy branches*)  c[[ci]] = b[[i]];
rorl[[ci]] = 1;
Ic[[ci]] = i;
ci++; (*there is only one (Subscript[c, i]) in a lazy branch*)
];
  If[gamma[[i]] > 0 && gamma[[i]] + alpha[[i]] < 1, (*Hanging increasing branches*)   c[[ci]] = b[[i]];
rorl[[ci]] = 1;
Ic[[ci]] = i;
ci++;   (*There are 2 (Subscript[c, i])in this interval so...*)
c[[ci]] = b[[i + 1]];
rorl[[ci]] = 0;
Ic[[ci]] = i;
ci++;
 ];
 If[alpha[[i]] < 1 && gamma[[i]] == 0,(*greedy branches*)   c[[ci]] = b[[i + 1]];
rorl[[ci]] = 0;
Ic[[ci]] = i;
ci++;(*Only one (Subscript[c, i]) in a greedy branch*)
 ];
 ];
If[beta[[i]] < 0,(*Decreasing branches*)  
If[alpha[[i]] < 1 && gamma[[i]] + alpha[[i]] == 1,(*lazy branches*)   c[[ci]] = b[[i+1]];
rorl[[ci]] = 0;
Ic[[ci]] = i;
ci++;
 ];
If[gamma[[i]] > 0 && gamma[[i]] + alpha[[i]] < 1,(*Hanging branches*)   c[[ci]] = b[[i]];
rorl[[ci]] = 1;
Ic[[ci]] = i;
ci++;
c[[ci]] = b[[i + 1]];
Ic[[ci]] = i;
rorl[[ci]] = 0;
ci++;
];
If[alpha[[i]] < 1 && gamma[[i]] == 0, (*greedy branches*)
c[[ci]] = b[[i]];
rorl[[ci]] = 1;
Ic[[ci]] = i;
ci++;
];
];
];


(*Finding (Subscript[a, i]). The "digits" to compute tao(x)*)
a=Table[,{N}];
For[i=1,i≤N,i++,
If[beta[[i]]>0,
a[[i]]=beta[[i]]*b[[i]]-gamma[[i]],
(*else*)
a[[i]]=beta[[i]]*b[[i]]-(gamma[[i]]+alpha[[i]]);
];
];

(*We need to know the interval in which x is in order to graph tao(x)*)
intervalx[x_]:=(
intx=1;
For[i=2,i≤N,i++,
If[x>b[[i]],(*In this case we consider the left interval for the (Subscript[c, i]) *)
intx=i;
];
];
intx
);
intervalxright[x_]:=(
intx=1;
For[i=2,i≤N,i++,
If[x≥b[[i]],(*In this case we consider the right interval for the (Subscript[c, i])*)
intx=i;
];
];
intx
);

(*This is the definition of tao. We have to make two definitions of tao because the definition itself depends on the interval in which the (Subscript[c, i]) is. This means one for the right (Subscript[c, i]) one for the left (Subscript[c, i]). This does not affect the graph of the map*)
tao=((beta[[intervalx[#]]]*#)-a[[intervalx[#]]])&;
taoright=((beta[[intervalxright[#]]]*#)-a[[intervalxright[#]]])&;
(*We create a vector for all the firts images of tao*)
taoc=Table[,{kc}];
(*And now we compute the images of the (Subscript[c, i])'s *)
For[j=1,j≤kc,j++,
If[rorl[[j]]==0,
taoc[[j]]=tao[c[[j]]],
(*else*)
taoc[[j]]=taoright[c[[j]]];
];
];


(*We define a function positionofx to know in which interval the given x is. This will be used in the cumulative beta and in the tao^n. The variable
w is used to know if the (Subscript[c, i]) is the one in the left or in the right.*)
positionofx[x2_,w_]:=(
ii=1;
isci=0; (*This variable tells us if the x is a ci or not, so that we can take care of the tao^n later on. Sometimes it could be 
confused with greedy or lazy branches, but since this variable is just to assure the tale in the function tao^n, then there is no problem*)
positionofx2=intervalx[x2];
If[c[[kc]]==x2&&c[[kc-1]]≠c[[kc]],
If[rorl[[kc]]==1,
If[w==1,
positionofx2=Ic[[kc]],
positionofx2=Ic[[kc]]-1;
];
];
If[rorl[[kc]]==0,
If[w==0,
positionofx2=Ic[[kc]],
positionofx2=Ic[[kc]]+1;
];
];
If[positionofx2==N+1,
positionofx2=positionofx2-1;
];
isci=1;
];

While[ii≤kc-1,
If[c[[ii]]==x2,
isci=1;
If[c[[ii]]≠c[[ii+1]],
If[rorl[[ii]]==1,
If[w==1,
positionofx2=Ic[[ii]],
positionofx2=Ic[[ii]]-1;
];
];
If[rorl[[ii]]==0,
If[w==0,
positionofx2=Ic[[ii]],
positionofx2=Ic[[ii]]+1;
];
];
ii=ii+kc
,
(*If[c[[ii]]==c[[ii+1]],*)
If[w==0,
positionofx2=Ic[[ii]];
ii=ii+kc,
positionofx2=Ic[[ii+1]];
ii=ii+kc;
];
];
];
ii++;
];
positionofx2
);

taon[x_,n_,whichctao_]:=(
outtao=Table[,{2},{n}];(*The first row contains all the taos, starting from tao^1 to tao^n. The second row contains the information on the tales 0 if one is taking the left hand side (Subscript[c, i]) and 1 for the right hand side. The first input of the second row represents the tale of tao^1 and the last input represents the tale of tao^n*)

p=positionofx[x,whichctao];
If[isci==1,
If[whichctao==0,
outtao[[1]][[1]]=tao[x];
If[beta[[p]]>0,
outtao[[2]][[1]]=0,
outtao[[2]][[1]]=1;
],
outtao[[1]][[1]]=taoright[x];
If[beta[[p]]>0,
outtao[[2]][[1]]=1,
outtao[[2]][[1]]=0;
];
]
,
(*else--if the x is not a (Subscript[c, i])*)
outtao[[1]][[1]]=tao[x];
outtao[[2]][[1]]=whichctao;
];
(*for n greater than 1. We need to change the value of the image in the first row and see if we need to change the value of the tale in the second row*)
For[j=2,j≤n,j++,
If[outtao[[2]][[j-1]]==0,(*left hand side (Subscript[c, i])*)
outtao[[1]][[j]]=tao[outtao[[1]][[j-1]]];
p=positionofx[outtao[[1]][[j]],0];
If[isci==0,
outtao[[2]][[j]]=outtao[[2]][[j-1]];
];
If[isci==1,
If[beta[[p]]>0,
outtao[[2]][[j]]=0,
outtao[[2]][[j]]=1;
];
],
(*right hand side (Subscript[c, i])*)
outtao[[1]][[j]]=taoright[outtao[[1]][[j-1]]];
p=positionofx[outtao[[1]][[j]],1];
If[isci==0,
outtao[[2]][[j]]=outtao[[2]][[j-1]];
];
If[isci==1,
If[beta[[p]]>0,
outtao[[2]][[j]]=1,
outtao[[2]][[j]]=0;
];
];
];
];
outtao
);

(*We need the cumulative beta*)
cumbeta[x_,n_,whichcbeta_]:=(
jj=1;
outbeta=Table[,{n}];

q=taon[x,n,whichcbeta];(*we need all the taos to compute the cumulative slope*)
outbeta[[1]]=beta[[positionofx[x,whichcbeta]]];
For[jj=2,jj≤n,jj++,
outbeta[[jj]]=outbeta[[jj-1]]*beta[[positionofx[q[[1]][[jj-1]],q[[2]][[jj-1]]]]];
];
outbeta
);



(*Defining the S matrix and D vector*)
sumend=50; (*There are two infinite sums in the program, sumend is the number of terms we will have on each sum*)

S=Table[,{kc},{kc}];
For[k=1,k≤kc,k++,
For[l=1,l≤kc,l++,
If[rorl[[k]]==0,
cumbetaci=cumbeta[c[[k]],sumend,0];
qci=q; (*these are the tao^n from the function cumbeta*)
S[[k]][[l]]=Sum[(1/Abs[cumbetaci[[s]]])*(Boole[cumbetaci[[s]]>0]*Boole[qci[[1]][[s]]-(10^(-100))>c[[l]]]+Boole[cumbetaci[[s]]<0]*Boole[qci[[1]][[s]]+(10^(-100))<c[[l]]]),{s,1,sumend}]
,
cumbetaci=cumbeta[c[[k]],sumend,1];
qci=q; 
S[[k]][[l]]=Sum[(1/Abs[cumbetaci[[s]]])*(Boole[cumbetaci[[s]]<0]*Boole[qci[[1]][[s]]-(10^(-100))>c[[l]]]+Boole[cumbetaci[[s]]>0]*Boole[qci[[1]][[s]]+(10^(-100))<c[[l]]]),{s,1,sumend}];
];
];
];

(*The transpose plus the identity matrix and define v, a vector of ones*)
S2=Table[,{kc},{kc}];
v=Table[,{kc}];
For[i=1,i≤kc,i++,
For[j=1,j≤kc,j++,
S2[[i]][[j]]=-S[[j]][[i]];
];
S2[[i]][[i]]=S2[[i]][[i]]+1;
v[[i]]=1;
];

(*Now the vector D and the constant (Subscript[D, 0])*)
dsubzero=1;
eigen=Eigenvalues[S];
is1eigen=0;
For[k=1,k≤Length[eigen],k++,
If[ eigen[[k]]==1,
is1eigen=1;
dsubzero=0;
];
];

D=Table[,{kc}];
D=LinearSolve[S2,dsubzero*v];
(*
If[is1eigen==1,
D=LinearSolve[S2,0*v],
D=LinearSolve[S2,v];
];*)


(*the chi function*)
chi=Boole[#1<#2&&#1>#3]&;
chiofbeta[chib_,chix_,inx_]:=(
If[chib>0,
outchi=chi[inx,chix,0],
outchi=chi[inx,1,chix];
];
outchi
);

(*and finally...the h function*)
h=dsubzero;
For[kk=1,kk<=kc,kk++,
If[rorl[[kk]]==0,
betachi=cumbeta[c[[kk]],sumend,0];
taonchi=q;
h=h+(D[[kk]]*(Sum[chiofbeta[betachi[[s]],taonchi[[1]][[s]],x1]/Abs[betachi[[s]]],{s,1,sumend}])),

betachi=cumbeta[c[[kk]],sumend,1];
taonchi=q;
h=h+(D[[kk]]*(Sum[chiofbeta[-betachi[[s]],taonchi[[1]][[s]],x1]/Abs[betachi[[s]]],{s,1,sumend}]));
];
];

(*this is the test for the expansion of x*)
prueba=a[[positionofx[#1,whichc1]]]/beta[[positionofx[#1,whichc1]]]+Sum[a[[positionofx[taon[#1,s-1,whichc1][[1]][[s-1]],whichc1]]]/cumbeta[#1,s,whichc1][[s]],{s,2,sumend}]&;


(*The graph of the map taon[c[[2]],3,whichc]&&cumbeta[c[[2]],3,whichc]&&*)
graph=Plot[{tao[x],x,1},{x,0,1},PlotStyle->{Thick,Green,Black}, Exclusions->b,ImageSize->Medium,AspectRatio->1];
(*graphprueba=Plot[{prueba[x],x,1},{x,0,1},PlotStyle->{Thick,Red,Green},Mesh->All,MaxRecursion->0,PlotPoints->50,Exclusions->b];*)

(*This is the output, we can see h(x1),kc,Ic, c, rorl,taoc,D and the graph, etc*)
out={h,b,kc,Ic,c,rorl,taoc,D,graph,prueba[x1],a,S,taon[x1,5,whichc1](*,graphpruebataon[x1,4,whichc1]*)}

(*out={h,b,kc,Ic,c,rorl,taoc,D,graph,a,S,eigen,taon[x1,1,whichc1],cumbeta[x1,5,whichc1],graphpruebataon[x1,4,whichc1]}*)


];
out
];

total[NT_,alphaT_,betaT_,gammaT_,x1T_,whichc1T_]:=(
(*Computing the normalizing constant*)
normconst=N[Integrate[hfunction[NT,alphaT,betaT,gammaT,x,whichc1T][[1]],{x,0,1}]];

(*Checking the error to see if h is an invariant density*)
numb=10;(*number of tries to check the h funtcion*)
h1=hfunction[NT,alphaT,betaT,gammaT,x1T,whichc1T];
b=h1[[2]];
a=h1[[11]];
y=Table[,{numb}];
err=Table[,{numb}];
pre=Table[,{NT}];
For[j6=1,j6≤numb,j6++,
y[[j6]]=(j6-1)/numb+0.1*RandomReal[];
];
For[j6=1,j6≤numb,j6++,
For[i3=1,i3≤NT,i3++,
pre[[i3]]=(y[[j6]]+a[[i3]])/betaT[[i3]];
];
su=0;
For[i3=1,i3≤NT,i3++,
If[(pre[[i3]]≥b[[i3]])&&(pre[[i3]]≤b[[i3+1]]),
su=su+(hfunction[NT,alphaT,betaT,gammaT,pre[[i3]],whichc1T][[1]])/Abs[betaT[[i3]]];
];
];
err[[j6]]=hfunction[NT,alphaT,betaT,gammaT,y[[j6]],whichc1T][[1]]-su;
];
valuesofh=Table[,{50}];
For[i4=1,i4≤49,i4++,
valuesofh[[i4]]=N[hfunction[NT,alphaT,betaT,gammaT,(i4/50),whichc1T][[1]]]/normconst;
];
valuesofh[[50]]=N[hfunction[NT,alphaT,betaT,gammaT,1-10^(-7),whichc1T][[1]]]/normconst;
(*The output*)

Print[Column[{"alpha=",alphaT,"beta=",betaT,"gamma=",gammaT,"partition points=",b,"graph=",h1[[9]],"(Subscript[c, i]) points=",N[h1[[5]]],"rightorleft endpoint=",N[h1[[6]]],"S=",MatrixForm[N[h1[[12]]]],"eigenvalues of S=",N[h1[[13]]],"D=",N[h1[[8]]],"normalizing constant=",normconst,"errors=",N[err],"values of h(x)=",valuesofh}]];
);

(*First example*)
total[5,{1,35/100,8/10,1,3/10},{3,3,-4,-5,-2},{0,2/10,1/10,0,7/10},2/3,1]
Plot[{hfunction[5,{1,35/100,8/10,1,3/10},{3,3,-4,-5,-2},{0,2/10,1/10,0,7/10},x,0][[1]]/normconst,0},{x,0,1},PlotPoints->25,AspectRatio->1]

(*Second Example*)
total[4,{7/10,2/10,1,45/100},{2,3,4,135/100},{0,2/10,0,55/100},35/100,0]
Plot[{hfunction[4,{7/10,2/10,1,45/100},{2,3,4,135/100},{0,2/10,0,55/100},x,1][[1]]/normconst,0},{x,0,1},PlotPoints->24,AspectRatio->1]

(*Example 3*)
total[4,{1,5/10,5/10,7/10},{4,3,2,21/10},{0,0,0,3/10},7/10,1]
Plot[{hfunction[4,{1,5/10,5/10,7/10},{4,3,2,21/10},{0,0,0,3/10},x,0][[1]]/normconst,0},{x,0,1},PlotPoints->25,AspectRatio->1]

(*Example EM*)
epst=5/1000
total[2,{1/2+epst,1},{1+2*epst,-2},{1/2-epst,0},0,0]
Plot[{hfunction[2,{1/2+epst,1},{1+2*epst,-2},{1/2-epst,0},x,0][[1]]/normconst,0},{x,0,1},PlotPoints->25,AspectRatio->1]

Clear["Global`*"]