#OK to post homework
#Jingze Li(RUID:195000313),04/14/19
#############READ ME###########
#the homework 21 consists of 4 parts, corresponding to 1 to 4 problems of hw21.

# a pdf file is attached for the 4th part ploting.

################
Help:=proc() print(`EvalFCdiffk(L,x,x0,k),CountTerm2(N),CountCoeff2(N),CountTerm3(N),CountCoeff3(N),EvalNR2(L,c1, c2, x0,y0),WeightRand(N,K)`): end:


#######  1 ##########


#EvalFCdiffk(L,x,x0,k): inputs a functional chain L, variable x, and number or symbol x0
#outputs the k th derivative of EvalFC(L,x,y) w.r.t y and y=x0
#

EvalFCdiffk:=proc(L,x,x0,k) local i,Ld,M,L0:

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

Ld:=[seq(diff(L[i],[x$k]),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:


########### 2 ################

CountTerm2 :=proc(N)
local f,g,i:
[seq(nops(convert(diff(f(g(x)),x$i),list)),i=1..N)]:
end:

# the list of the number of terms in df(g(x))dx is not in the OEIS.
# I think it  make little sense since when we take different f and g specificly, the number of terms in df(g(x))dx varies. e.g. f = x*exp(x) give two terms when taking derivative.  


CountCoeff2 :=proc(N)
local f,g,i:
[seq(SumCoeff(diff(f(g(x)),x$i)),i=1..N)]:
end:

# the list of the sum of coefficients in df(g(x))dx is in the OEIS

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


############# 3 ################


CountTerm3 :=proc(N)
local f,g,i:
[seq(nops(convert(diff(f(g(h(x))),x$i),list)),i=1..N)]:
end:

# similar to 2 , no result for the sequence of term number in OEIS

CountCoeff3 :=proc(N)
local f,g,i:
[seq(SumCoeff(diff(f(g(h(x))),x$i)),i=1..N)]:
end:

# it is in OEIS: a way to interpret is number of 3-level labeled rooted trees with n leaves

########### 4 #############

EvalNR2:=proc(L,c1, c2, x0,y0)
local x,y,i,A,B:
x:=x0;
y:=y0:
for i from 1 to nops(L) do 
  A :=L[i][1]:
  B :=L[i][2]:
  x := max(A[1][1]*x+A[1][2]*y+B[1][1],0):
  y := max(A[2][1]*x+A[2][2]*y+B[1][2],0):
od:

max(c1*x+c2*y,0):

end:


# WeightRand(N,K) output a list of random weight matrix 

WeightRand:=proc(N,K)
local i,A,B:
A:= [seq(RandMat(2,2,K),i=1..N)]:
B:= [seq(RandMat(1,2,K),i=1..N)]:
[seq([A[i],B[i]],i=1..N)]:
end:


# I plot some functions with different random weight matrix and different number of layers in the pdf file attached. 


######### from c20.txt and c21.txt #######



#C20.txt,April 8, 2019
Help:=proc(): print(`BT(n) ,NuLeaves(T) , RandMat(m1,m2,K) `): 
print(`RandMatList(P,K) , EvalBT(L,T), PriceM(P,T) `):
end:

with(linalg):

#Lemma: The cost (in multiplication operations) of multiplying an mxn matrix by an nxp matrix
#is mnp
#Cor. The cost of A=[m1xm2], B=[m2xm3] C=[m3xm4] 
#the cost of (AB)C=m1m2m3+ m1m3m4
#The cost of A(BC)=m2m3m4+ m1m2m4
#A1, A2, ..., An, [m1,m2,m3, ..., m[n+1]] 

#input a positive integer n
#returns the set of all binary trees with exactly n leaves 
BT:=proc(n) local S,a,b,i:

	option remember:

	if n<=0 then
		return {}:
	fi:

	if n=1 then 
		return {[]}:
	fi:
	
	S:={}:
	
	#to get a tree on n leaves, you can join any two trees whose numbers of leaves sum to n 
	for i from 1 to n-1 do 
	
		for a in BT(i) do
		
			for b in BT(n-i) do
			
				S:=S union {[a,b]}:
			
			od:
		
		od:
	
		
	
	od:
	
	return S:

end:

#NuLeaves(T): the number of leaves in the binary tree T
NuLeaves:=proc(T)
option remember:
if T=[] then
 RETURN(1):
else
 NuLeaves(T[1])+NuLeaves(T[2]):
fi:

end:



#RandMat(m1,m2,K): a random m1 by m2 matrix with entries from -K to K
#given as a list of lists 
RandMat:=proc(m1,m2,K) local ra,i,j:
ra:=rand(-K..K):
[seq( [seq( ra(), j=1..m2)],i=1..m1)]:
end:

#RandMatList(P,K): given a list of positive integers P of length n+1
#and a pos. integer K outputs a list of n compatible matrices ready to
#be multiplied
RandMatList:=proc(P,K) local n,i:
n:=nops(P)-1:
[seq(RandMat(P[i],P[i+1],K),i=1..n)]:
end:


#EvalBT(inputs a list L of matrices and a binary tree with the number of leaves equal to
#the length of L, outputs the matrix product L[1]...L[n] suggested by the tree
EvalBT:=proc(L,T) local n,P,i,P1,k,A,B:
option remember:
n:=nops(L):
if NuLeaves(T)<>n then
 RETURN(FAIL):
fi:

P:=[seq([nops(L[i]),nops(L[i][1])],i=1..n)]:

for i from 1 to n-1 do
 if P[i][2]<>P[i+1][1] then
   ERROR(L[i], `is not compatible with`, L[i+1] ):
 fi:
od:

#P1:=[seq(P[i][1],i=1..n-1),P[n][1],P[n][2]]:

if n=1 then
  RETURN(L[1]):
fi:

#k:=number of leaves of the left tree of the root
k:=NuLeaves(T[1]):
A:=EvalBT( [op(1..k,L)] ,  T[1]):
B:=EvalBT([op(k+1...n ,L)] ,T[2]):
convert(multiply(A,B),list):
end:

#PriceM(P,T): inputs a profile of lenght n+1 and a binary tree with n leaves
#outputs the number of multiplications
PriceM:=proc(P,T) local n,k:
n:=NuLeaves(T):
if n=1 then
 RETURN(0):
fi:

if n=2 then
 RETURN(P[1]*P[2]*P[3]):
fi:

k:=NuLeaves(T[1]):


PriceM([op(1..k+1,P)],T[1])+PriceM([op(k+1..n+1,P)],T[2])+ P[1]*P[k+1]*P[n+1]:
   
end:




##### c21.txt #########


#C21.txt, April 11, 2019, Functional chains (and networks)
Help:=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: