#OK to post homework
#TianHao Zhang,04/14/19
#################################Help#######################################
Help:=proc():print(`EvalFCdiffk(L,x,x0,k) EvalFCdiffk_sym(L,x,k)`) end:
#################################Part1######################################
#EvalFCdiffk(L,x,x0,k): inputs a functional chain L, variable x, and number or symbol x0 and outputs the k-th derivative of EvalFC(L,x,y) w.r.t. y and y=x0.
#By Bruno's formula d^nf(g(x))/dx^n = sum(f^(k)(g(x))*B_{n,k}(g'(x),..,g^(n-k+1)(x)))

EvalFCdiffk:=proc(L,x,x0,k) local i, Ld, M, L1, n, g, j:

n:=nops(L):
if n = 1 then
  return(subs(x=x0,diff(L[1], [x$k]))):
fi:

L1:=[op(2..n,L)]:
g:=EvalFC(L1,x,x0):
add(subs(x=g,diff(L[1], [x$i]))*IncompleteBellB(k, i, seq(EvalFCdiffk(L1,x,x0,j),j=1..k-i+1)),i=1..k):

end:

#################################Part2######################################
#Find out whether the number of terms of diff(f(g(x)),x$i) and the sum of its coefficients are in the OEIS.
#Before dealing with the problem, we transfer the function in Part 1 into a symbolic #function which use 'Diff' instead of 'diff' 
EvalFCdiffk_sym:=proc(L,x,k) local i, Ld, M, L1, n, g, j:

n:=nops(L):
if n = 1 then
  return(Diff(L[1], [x$k])):
fi:

L1:=[op(2..n,L)]:
g:=EvalFC(L1,x,x):
expand(add(Diff(L[1], [x$i])*IncompleteBellB(k, i, seq(EvalFCdiffk_sym(L1,x,j),j=1..k-i+1)),i=1..k)):

end:

#####
#seq(SumCoeff(EvalFCdiffk_sym([f, g], x, k)), k = 2 .. 10);
#          2, 5, 15, 52, 203, 877, 4140, 21147, 115975

#In OEIS:
#A000110		
#Bell or exponential numbers: number of ways to partition a set of n labeled elements. 

#################################Part3######################################
#Find out whether the number of terms of diff(f(g(h(x))),x$i) and the sum of its coefficients are in the OEIS.

#seq(SumCoeff(EvalFCdiffk_sym([f, g, h], x, k)), k = 2 .. 10);
#     3, 12, 60, 358, 2471, 19302, 167894, 1606137, 16733779

#In OEIS:
#A000258		
#E.g.f.: exp(exp(exp(x)-1)-1).  (Formerly M2932 N1178)
###############################From C21#####################################
#C21.txt, April 11, 2019, Functional chains (and networks)
Help21:=proc(): print(`  EvalFC1(L,x,0), EvalFC(L,x,x0) , RP(d,K,x) , RC(d,K,n,x) `): 
print(`EvalFCdiff(L,x,x0), SumCoeff(A) `):
end:

#y=a*x+b #Data-Set [x1,y1], ..., [xn,yn]  L(a,b):-=Sum(a*xi+b-yi)^2,i=1..n);
#diff(L,a)=0, diff(L,b)=0   a*x_{n+1}+b
# [x1,..., xm; y]: Loss(a1, ..., am,a0):=Sum((a1*x1+...+am*xm+a0-y)^2 , data points in the set)
#[m1,m2,m3, ..., mk] m1xm2 +m2xm3+....
#CNN
#functional chain x->f1(x)->f2(f1(x))
#Def: A functional chain is a list of expressions [f1, ..., fn] in x

#EvalFC1(L,x,x0): inputs L=[f1,..., fn] expressions in x and a number (or symbol) x0
#outputs f1(f2(...fn(x0))))) 
EvalFC1:=proc(L,x,x0) local i,L1:

if nops(L)=1 then
   RETURN(subs(x=x0,L[1])):
fi:

L1:=[op(1..nops(L)-1,L)]:
 expand(subs( x=EvalFC1(L1,x,x0),L[nops(L)])):


end:

#RP(d,K,x): a random polynomial of degree d in x with coeff. between -K and K
RP:=proc(d,K,x) local ra,i:
ra:=rand(-K..K):
add(ra()*x^i,i=0..d):
end:

#RC(d,K,n,x): a random functional chain of length n
RC:=proc(d,K,n,x) local i:
[seq(RP(d,K,x),i=1..n)]:
end:

#f(g(x))'= f'(g(x))*g'(x)
#f1(f2(f3(x))'=f1'(f2(f3(x))*f2'(f3(x))*f3'(x)
#(fk(f(k-1)(f(k-2)))'= fk'(f(k-1)....) [(f(k-1)(f(k-2)....) f1'(x0)]

#x0->f1(x0)->f2(f1(x0)->.... fn(f_{n-1}(x)) ...
#EvalFC(L,x,x0): inputs L=[f1,..., fn] expressions in x and a number (or symbol) x0
#outputs [f1(x0),f2(f1(x0)), f3(f2(f1(x0)), ..., fn(....)]: 
EvalFC:=proc(L,x,x0) local i,L1:

if nops(L)=1 then
   RETURN([subs(x=x0,L[1])]):
fi:

L1:=[op(1..nops(L)-1,L)]:

L1:=EvalFC(L1,x,x0):


[op(L1),subs(x=L1[nops(L1)] , L[nops(L)])]:

end:

#EvalFCdiff(L,x,x0): inputs a functional chain L, variable x, and number or symbol x0
#outputs the derivative of EvalFC(L,x,y) w.r.t y and y=x0
EvalFCdiff:=proc(L,x,x0) local i,Ld,M,L0:

L0:=EvalFC(L,x,x0):

Ld:=[seq(diff(L[i],x),i=1..nops(L))]:

M:=[subs(x=x0,Ld[1])]:

for i from 2 to nops(L) do
 M:=[op(M),  subs(x=L0[i-1]  ,Ld[i])  ]:
od:

convert(M,`*`):

end:

#SumCoeff(A): the sum of the coefficient in a + expressin
SumCoeff:=proc(A) local su,i, lauren:
su:=0:
if not type(A,`+`) then
  RETURN(FAIL):
fi:
for i from 1 to nops(A) do
   lauren:=op(1,op(i,A)):
    if type(lauren,integer) then
       su:=su+lauren:
       else
       su:=su+1:
     fi:
od:
su:
end: