# 
# The following toy implementation is intended to illustrate 
# Example 5.6 of ``Modeling and Analysis of Timed Petri Nets using Heaps of Pieces'', S. Gaubert, J. Mairesse. Submitted to IEEE-TAC. The eprint can
# be found on this serveur http://amadeus.inria.fr/gaubert/PAPERS/GM97.ps
# 
# The following code has been tested with Maple version V.3.
# 
# Recall that this is only a fastly designed illustrative implementation
# (the goal is to understand, not to be efficient on real cases).
# It is possible to design much more efficient algorithms,
# based on some of the ideas here. This will be dealt
# with in detail in the DEA report of Eric Hautecloque
# (Ecole Nationale des Techniques Avancees & DEA 
# ``Mod\'elisation et M\'ethodes
# Math\'ematiques en \'Economie'', Universit\'e de Paris I),
# to appear in July 1997.
##################################################################
#
# A word is represented by a list. E.g. w=[a,b,c]
# A finite language is represented by a set of lists
# e.g. {[a,b],[c,d]}
#
#
#
# Compute the language: (shuffle product of u,v)*w
# arguments: u,v,w: words 
shuffle_words:=proc(u,v,w)
if nops(u)=0 then RETURN({[v[],w[]]})
elif nops(v)=0 then RETURN({[u[],w[]]})
fi;
shuffle_words([u[1..(nops(u)-1)]],v,[u[nops(u)],w[]]) union shuffle_words(u,[v[1..(nops(v)-1)]],[v[nops(v)],w[]]);
end:
#
# Compute the shuffle product of u and v
# arguments: u,v: languages
shuffle:=proc(u,v)
map(proc(x,y) map(proc(Y,X) shuffle_words(X,Y,[])[] end, y,x)[] end, u,v);
end:

# L: a language
# INDEP: an independence relation (set of couples of letters that commute).
SIMPLIFY:=proc(L,INDEP)
map(normalform,L,INDEP);
end:

# returns the lexicographic normal form of the word u
# in the trace monoid of independence relation INDEP
# (i.e. the least element of the equivalence class
# of x, for the alphabetic order).
#
# One could do much better than the following coarse (but correct)
# algorithm. More refined algorithms are detailed in the first
# chapter of 
# @BOOK{diekert,
# author={V. Diekert},
# title={Combinatorics on traces},
# year={1990},
# publisher={Springer},
# series={LNCS},
# number={454},
# }
#
normalform:=proc(u,INDEP)
local n,i,j,flag;
n:=nops(u);
if n=1 then RETURN(u)
else
  for i from n to 2 by -1 do
     j:=i-1;
     flag:=member({u[i],u[j]},INDEP);
     while flag do # we look backward to see if there is a letter
		   # that can commute with u[i], and that is 
		   # strictly less for the lexicographic order
	if (lexorder(u[i],u[j])) then # such a letter is found.
	RETURN(normalform([u[1..(j-1)],u[i],u[(j)..(i-1)],u[(i+1)..n]],INDEP));
		   # we exchange u[i] and u[j].
	fi;
	j:=j-1;
        if j=0 then
           flag:=false; # the begining of the word is reached,
			# u[i] commute with any u[j], with j<i,
			# but all such u[j] are less than u[i].
			# ---> continue with u[i-1]
        else
	   flag:=member({u[i],u[j]},INDEP);
	   # continue only if u[j] commutes with u[i].
        fi;
     od;
  od;
 RETURN(u);
fi;
end:

# an approximate but fast simplification, used in a first pass.
coarsenormalform:=proc(u,INDEP)
local n,i;
n:=nops(u);
if n=1 then RETURN(u)
else
  for i to n-1 do
    if (member({u[i],u[i+1]},INDEP) and lexorder(u[i+1],u[i]) ) then
          RETURN(coarsenormalform([u[1..(i-1)],u[i+1],u[i],u[(i+2)..n]],INDEP));
    fi;
  od;
 RETURN(u);
fi;
end:

normalsimpl:=proc(L,INDEP)
map(normalform,map(coarsenormalform,L,INDEP),INDEP);
end:

# returns true if the list x is less than the list y
less2:=proc(x,y) lexorder(cat(x[]),cat(y[])) end:
#

# two words x and y are conjugate if x=uv and y =vu
# minconjugate returns the minimal element in a conjugacy class.
# (this coarse algorithm runs in quadratic time).
minconjugate:=proc(x)
local n,y,z,i,j;
n:=nops(x);
y:=x;
for i to n-1 do 
   z:=[x[(i+1)..n],x[1..i]];
   if less2(z, y) then y:=z; 
   fi;
od;
y;
end:


# 
# 
makeind:=proc(a,b)
local A,B;
A:=convert(a[2],set);
B:=convert(b[2],set);
if (A intersect B)={} then RETURN({{a,b}}) else
RETURN({}); fi;
end:

#
# a: list of tasks
# generates the set of commutation relations. 
MAKEIND:=proc(a)
local i,j,S;
S:={};
for i to nops(a) do for j to nops(a) do
# adds the set { [a[i],b[j]],[b[j],a[i]]} to S if a[i] and b[j] 
# have no common resource.
S:=S union makeind(a[i],a[j]);
od;
od;
end:

# evaluates the performance of all the schedules within the list
# L. N gives the number of resources.
#
# nota bene: karp returns the maximal eigenvalue 
# of a strongly connected component accessible
# from node 1. 
#
# If all the "schedule matrices" are irreducible,
# we thus obtain the throughput.
# Otherwise, some "accessibility" analysis
# should be necessary (not implemented).
PERF:=proc(L,N)
map(proc(x,N) local k; \
k:=Karp(perfs(x,N));\
print(x,'k'=k);\
RETURN([x,k]);\
end, L,N);
end:


# task matrix == array([n,R, T]);
# where: 
# n= number of resources required, 
# R= list of resources
# T= list of corresponding release times.
# Note that the total number of resources N
# is not given in the task matrix.
# (we must have R[i]<= N in the sequel).
#
# 
# Sparse_TR(matrix, task Matrix) --> matrix
# [size of matrix and numbers of resources must be compatible.]
# computes the product of matrix q and task matrix A in a
# sparse way.
#
Sparse_TR:=proc(q,A)
local qq,n,R,i,j,k;
qq:=copy(q);
n:=A[1];
for i to linalg[rowdim](q) do 
	R:=q[i,A[2][1]];
	for k from 2 to n do 
		R:=max(R, q[i,A[2][k]]);
	od;
	for j to n do 
		qq[i,A[2][j]]:= simplify(R + A[3][j]);
	od;
od;
op(qq);
end:

# identity in max algebra
id:=proc(n)
    local c,i,j;
    if type(n,integer) then
    c:=array(1..n,1..n);
	for i from 1 to n do
	    for j from 1 to n do
		c[i,j]:=-infinity;
	    od;
	od;
	for i from 1 to n do c[i,i]:=0;
	od;	
	RETURN(op(c));
    else ERROR(`argument must be an integer`);
    fi;
end:
# perfs(schedule, number of resources) --> matrix
# x: must be a vector representing the schedule
# x[i]=name of the ith task matrix to be
# executed.
perfs:=proc(x,n)
local q,i;
q:=id(n);
for i to linalg[vectdim](x) do
q:=Sparse_TR(q, x[i]);
od;
end:


#
# Computing the eigenvalue of an irreducible matrix a in the
# max+ semiring.
# 
# (More generally, Karp returns the maximal eigenvalue
# of a strongly connected component accessible
# from index 1).
#
Karp:=proc(a)
local f,n,k,i,j,t;
if type(a,matrix) then
        if not(type(a[1,1],realcons)) then ERROR(`entry 11 of A is not a real constant`); fi;
	n:=linalg[rowdim](a);
	f:=array(1..n+2,1..n);
	f[1,1]:=0;
	for i from 2 to n do f[1,i]:=-infinity;
	od;
	# for each value k and node i,
        # we compute the distance from 1 to i in k steps.
	for k from 2 to n+1 do
		for i from 1 to n do
		f[k,i]:=-infinity;
			for j from 1 to n do
			f[k,i]:=max(f[k-1,j]+ a[j,i], f[k,i]);
			od;
		od;
	od;
	# then, the eigenvalue can be derived from the f[k,i] as follows.
	for i from 1 to n do
		f[n+2,i]:=infinity;
		for k from 1 to n do
		f[n+2,i]:=min(f[n+2,i],(1/(n-k+1) *(f[n+1,i] -  f[k,i] ) ));
		od;
	od;
	t:=f[n+2,1];
	for i from 2 to n do
	t:=max(t , f[n+2,i]);
	od;
	RETURN(t);
elif type(a,scalar) then RETURN(a);
else ERROR(`Invalid argument in karp: should be a scalar or a matrix`);
fi;
end: