######################################################################
##Skyscrapers.txt: Save this file as  Skyscrapers.txt                #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read Skyscrapers.txt<Enter>                        #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: June 13, 2016

with(combinat): 
print(`Created:  June 13, 2016`):
print(` This is Skyscrapers.txt `):
print(`A Maple package to solve (and create) Skyscrapers puzzles`):
print(` Written by Doron Zeilberger`):
print(`available from Zeilberger's website at the url`):
print(`http://www.math.rutgers.edu/~zeilberg/tokhniot/Skyscrapes.txt`):
print(``):
print(`This type of puzzles was invented by Wei-Hwa Huang and these puzzle appear (starting Spring 2016) `):
print(`in the New York Times Sunday magazine`):

print(``):
print(`Please report bugs to zeilberg at math dot rutgers dot edu`):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/  .`):

print(`---------------------------------------`):
 print(`For a list of the Supporting procedures type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

with(combinat):

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: Banim, IsLat, IsLegal, Kids, LegalC, LegalCs, LIS, Matim, Ptor, PTOR, Ptor1,RaLS,  Rev, ReLS,  Tav, Yel, YelC, YelR,  `):
 print(``):

else
ezra(args):
fi:

end:



ezra:=proc()

if args=NULL then
 print(`The main procedures are:  DrawP, DrawPS, MakePuzzle,  MakePuzzleB, SOLVE  `):
 print(` `):



elif nops([args])=1 and op(1,[args])=Banim then
print(`Banim(M,a): given a Latin rectangle (of size k by n say) , and an integer a, returns all the legal extensions with one extra row`):
print(`such that it is still a Latin rectangle, and if a=0 then the new row can be anything legal, but`):
print(`if a is between 1 and n then the new (k+1)-th row has length of the longest increasing subsequence equal  to a.`):
print(` Try: `):
print(` Banim([[1,2,3]],2); `):

elif nops([args])=1 and op(1,[args])=DrawP then
print(`DrawP(P): draws a puzzle P`):
print(`For the puzzle of June 19, 2016, type;`):
print(`DrawP([[0,3,0,0],[0,0,2,0],[0,3,2,0],[0,2,3,0]]);`):
print(``):
print(`for the puzzle of June 26, 2016, type`):
print(`DrawP([[2,0,2,0],[0,0,0,0],[0,0,0,0],[2,2,0,2]]);`):
print(``):
print(`for the puzzle of July 3, 2016, type`):
print(`DrawP([[0$4],[0,3,3,0],[0$4],[0,3,3,0]]);`);
print(``):
print(``):
print(`for the puzzle of July 10, 2016, type`):
print(`DrawP([[0,1,0,0],[0,0,1,0],[0$4],[2,1,2,0]]);`);
print(``):

elif nops([args])=1 and op(1,[args])=DrawPS then
print(`DrawPS(P,S): draws a puzzle P and its solution S`):
print(`For the puzzle of June 19, 2016, type;`):
print(`P1:=[[0,3,0,0],[0,0,2,0],[0,3,2,0],[0,2,3,0]]: S1:=PTOR1(P1): DrawPS(P1,S1);`):
print(``):
print(`for the puzzle of June 26, 2016, type`):
print(`P1:=[[2,0,2,0],[0,0,0,0],[0,0,0,0],[2,2,0,2]]:   S1:=PTOR1(P1): DrawPS(P1,S1);  `):
print(``):
print(`for the puzzle of July 3, 2016, type`):
print(`P1:=[[0$4],[0,3,3,0],[0$4],[0,3,3,0]]:    S1:=PTOR1(P1): DrawPS(P1,S1); `);
print(``):
print(``):
print(`for the puzzle of July 10, 2016, type`):
print(`P1:=[[0,1,0,0],[0,0,1,0],[0$4],[2,1,2,0]]:    S1:=PTOR1(P1): DrawPS(P1,S1); `);
print(``):


elif nops([args])=1 and op(1,[args])=IsLat then
print(`IsLat(M): inputs  a list of lists and outputs true if it is a Latin rectangle.`):
print(`Try: `):
print(` IsLat([[1,2,3],[2,3,1]]); `):


elif nops([args])=1 and op(1,[args])=IsLegal then
print(`IsLegal(M): inputs a list of lists M,representing a matrix with non-negative integers checks that every`):
print(`row and every column contains distinct positive integers. Try`):
print(`IsLegal([[0,1,2],[1,0,3],[2,1,0]]);`):

elif nops([args])=1 and op(1,[args])=Kids then
print(`Kids(M,a,b): given a Latin rectangle (of size k by n say) , and  integers a, and b returns all the legal extensions with one extra row`):
print(`such that it is still a Latin rectangle, and if a=0 and b=0 then the new row can be anything legal, but`):
print(`if a is between 1 and n then the new (k+1)-th row has length of the longest increasing subsequence equal  to a.`):
print(`and if b is between 1 and n the revese of the new (k+1)-th row has length of the longest increasing subsequence equal  to b`):
print(`Try: `):
print(`Kids([[1,2,3,4]],2,3); `):


elif nops([args])=1 and op(1,[args])=LegalC then
print(`LegalC(M): Given a Latin rectangle M, finds all the legal continuations. Try:`):
print(`LegalC([[1,2,3,4]]);`):

elif nops([args])=1 and op(1,[args])=LegalCs then
print(`LegalCs(S): Given a set of Latin rectangles S, finds all the legal continuations. Try:`):
print(`LegalCs({[[1,2,3,4]]});`):

elif nops([args])=1 and op(1,[args])=LIS then
print(`LIS(pi): inputs a permutation pi and outputs the size of the largest increasing subsequence (alias the number of`):
print(`sky-scrapers visible from the left)`):

elif nops([args])=1 and op(1,[args])=MakePuzzle then
print(`MakePuzzle(n,K): makes a random puzzle, by trying K times. If it fails, it returns FAIL. Try:`):
print(` MakePuzzle(4,20);`):

elif nops([args])=1 and op(1,[args])=MakePuzzleB then
print(`MakePuzzleB(n,K1,K2): inputs a positive integer n, and positive integers K1, K2, `):
print(`outputs a set of size K1 of pairs [P,S] of [puzzle,solution], by trying MakePuzzle(n,K2) (q.v.) as`):
print(` many times as needed. Try: `):
print(` MakePuzzleB(4,10,20); `):

elif nops([args])=1 and op(1,[args])=Matim then
print(`Matim(Le,L): inputs two lists Le,L, of the same size that outputs true if every non-zero entry of Le equals the corresponding`):
print(`entry of L. Try:`):
print(`Matim([0,3,0,2],[1,3,4,2]);`):

elif nops([args])=1 and op(1,[args])=PTOR then
print(`PTOR(P): Like SOLVE(P), using a different approach `):
print(`inputs a list of length 4, P=[U,D,L,R], each of them a list length n, say, whose entries are integers from 0 to n and`):
print(`prints out the set of Latin n by n squares, `):
print(`such that if `):
print(`U[i] is not 0, then the lenght of the longest increasing sequence of the i-th column is a[i]. `):
print(`D[i] is not 0, then the lenght of the longest increasing sequence of the reverse of the i-th column is a[i]. `):
print(`L[i] is not 0, then the lenght of the longest increasing sequence of  the i-th row is a[i]. `):
print(`R[i] is not 0, then the lenght of the longest increasing sequence of  the reverse of the i-th row is a[i]. `):
print(`For example, to solve the Skyscrapers puzzle in the New York Times Magazine `):
print(`of June 19, 2016, type;`):
print(`PTOR([[0,3,0,0],[0,0,2,0],[0,3,2,0],[0,2,3,0]]);`):
print(``):
print(`for the puzzle of June 26, 2016, type`):
print(`PTOR([[2,0,2,0],[0,0,0,0],[0,0,0,0],[2,2,0,2]]);`):
print(`for the puzzle of July 3, 2016, type`):
print(`PTOR([[0$4],[0,3,3,0],[0$4],[0,3,3,0]]);`);
print(``):





elif nops([args])=1 and op(1,[args])=Ptor1 then
print(`Ptor1(a): inputs a list of length n, say, whose entries are integers from 0 to n outputs the`):
print(`set of Latin n by n squares, such that if a[i] is not 0, then the lenght of the longest increasing`):
print(`sequence of the i-th row is a[i]. For example, to solve the Skyscrapers puzzle in the New York Times`):
print(` magazine of June 12, 2016, type; `):
print(`Ptor1([1,4,3,0]);`):

elif nops([args])=1 and op(1,[args])=RaLS then
print(`RaLS(n): A random n by n Latin Square. Try:`):
print(`RaLS(5);`):


elif nops([args])=1 and op(1,[args])=ReLS then
print(`ReLS(n): the set of all n by n Latin Squares whose first row is 1...n. Try:`):
print(`ReLS(4);`):

elif nops([args])=1 and op(1,[args])=Rev then
print(`Rev(pi): the reverse of the permutation pi`):


elif nops([args])=1 and op(1,[args])=SOLVE then
print(`SOLVE(P): inputs a list of length 4, P=[U,D,L,R], each of them a list length n, say, whose entries are integers from 0 to n outputs the`):
print(`set of Latin n by n squares, `):
print(`such that if `):
print(`U[i] is not 0, then the lenght of the longest increasing sequence of the i-th column is a[i]. `):
print(`D[i] is not 0, then the lenght of the longest increasing sequence of the reverse of the i-th column is a[i]. `):
print(`U[i] is not 0, then the lenght of the longest increasing sequence of the i-th column is a[i]. `):
print(`L[i] is not 0, then the lenght of the longest increasing sequence of  the i-th row is a[i]. `):
print(`R[i] is not 0, then the lenght of the longest increasing sequence of  the reverse of the i-th row is a[i]. `):
print(`For example, to solve the Skyscrapers puzzle in the New York Times Magazine `):
print(`of June 19, 2016, type;`):
print(`SOLVE([[0,3,0,0],[0,0,2,0],[0,3,2,0],[0,2,3,0]]);`):
print(``):
print(`for the puzzle of June 26, 2016, type`):
print(`SOLVE([[2,0,2,0],[0,0,0,0],[0,0,0,0],[2,2,0,2]]);`):
print(``):
print(`for the puzzle of July 3, 2016, type`):
print(`SOLVE([[0$4],[0,3,3,0],[0$4],[0,3,3,0]]);`);
print(``):
print(``):
print(`for the puzzle of July 10, 2016, type`):
print(`SOLVE([[0,1,0,0],[0,0,1,0],[0$4],[2,1,2,0]]);`);
print(``):

print(``):
print(`for the puzzle of July 17, 2016, type`):
print(`SOLVE([[0,5,0,0,0],[0$5],[3,5,0,3,0],[0,0,0,3,3]]);`):
print(``):


print(``):
print(`for the puzzle of July 24, 2016, type`):
print(`SOLVE([[0,0,0,4,0],[3,0,0,0,0],[3,0,4,0,4],[0,3,0,0,0]]);`):
print(``):

print(``):
print(`for the puzzle of July 31, 2016, type`):
print(`SOLVE([[0,4,0,0,0],[0,0,0,4,0],[0,0,0,3,0],[0,4,0,2,0]]);`):
print(``):


print(``):
print(`for the puzzle of Aug. 7, 2016, type`):
print(`SOLVE([[5,5,0,0,0,0],[0,0,0,5,0,0],[0,6,0,0,0,0],[0,0,4,0,0,0]]);`):
print(``):

print(``):
print(`for the puzzle of Aug. 14, 2016, type`):
print(`SOLVE([[0,5,0,0,0,0],[0,0,4,0,0,0],[3,0,4,0,0,6],[0,0,0,0,5,0]]);`):
print(``):


print(``):
print(`for the puzzle of Aug. 21, 2016, type`):
print(`SOLVE([[3,0,5,0,5,0],[0,3,0,0,0,0],[3,0,0,5,0,3],[0,0,5,0,0,0]]);`):
print(``):




elif nops([args])=1 and op(1,[args])=Tav then
print(`Tav(n): inputs a positive integer n, and outputs a list of sets of length n, such that the i-th`):
print(`entry is the set of n-permutations whose length of largest  increasing sequence is i.`):
print(` Try: `):
print(`Tav(4); `):

elif nops([args])=1 and op(1,[args])=Tavla then
print(`Tavla(n): inputs a positive integer n, and outputs a list of (length n) of lists of sets of length n, such that the i-th`):
print(`j-th`):
print(`entry is the set of n-permutations whose length of largest  increasing sequence is i and the length of the reverse is j`):
print(`Tavla(4);`):

elif nops([args])=1 and op(1,[args])=Yel then
print(`Yel(P,M): inputs a puzzle P=[U,D,L,R], and a partial solution M, and a column j from 1 to nops(P[1]) outputs`):
print(`the children from the best row and column. Try`):
print(`Yel([[0,3,0,0],[0,0,2,0],[0,3,2,0],[0,2,3,0]],[[0$4]$4]);`):


elif nops([args])=1 and op(1,[args])=YelC then
print(`YelC(P,M,i): inputs a puzzle P=[U,D,L,R], and a partial solution M, and a column j from 1 to nops(P[1]) outputs`):
print(`all its legal children. Try:`):
print(`YelC([[0,3,0,0],[0,0,2,0],[0,3,2,0],[0,2,3,0]],[[0$4]$4],2);`):

elif nops([args])=1 and op(1,[args])=YelR then
print(`YelR(P,M,i): inputs a puzzle P=[U,D,L,R], and a partial solution M, and a row i from 1 to nops(P) outputs`):
print(`all its legal children. Try:`):
print(`YelR([[0,3,0,0],[0,0,2,0],[0,3,2,0],[0,2,3,0]],[[0$4]$4],2);`):

print(``):


else
print(`There is no ezra for`,args):
fi:
 
end:


#Rev(pi): the reverse of the permutation pi
Rev:=proc(pi) local i,n:
n:=nops(pi):
[seq(pi[n+1-i],i=1..n)]:
end:

#LIS(pi): inputs a permutation pi and outputs the size of the largest increasing subsequence (alias the number of
#sky-scrapers visible from the left)
LIS:=proc(pi) local i,lu:
lu:=[pi[1]]:

for i from 2 to nops(pi) do
 if lu[nops(lu)]<pi[i] then
  lu:=[op(lu),pi[i]]:
 fi:
od:

nops(lu):

end:

with(plots):

#DrawP(P): draws a puzzle P
DrawP:=proc(P) local U1,D1,L1,R1,n,pic,i:
U1:=P[1]: D1:=P[2]: L1:=P[3]: R1:=P[4]:
n:=nops(U1):

pic:=plot([[0,0],[n,0]],axes=none,scaling=CONSTRAINED,color=blue):
#pic:=plot([[0,0],[n,0]],axes=none,color=blue):

for i from 1 to n do
pic:=pic, plot([[0,i],[n,i]],axes=none,scaling=CONSTRAINED,color=blue):
od:

for i from 0 to n do
pic:=pic, plot([[i,0],[i,n]],axes=none,scaling=CONSTRAINED,color=blue):
od:

for i from 1 to n do

  if D1[i]<>0 then
     pic:=pic,textplot([i-1/2,-1/4,D1[i]]):
  fi:
od:

for i from 1 to n do

  if U1[i]<>0 then
     pic:=pic,textplot([i-1/2,n+1/4,U1[i]]):
  fi:
od:

for i from 1 to n do

  if L1[i]<>0 then
     pic:=pic,textplot([-1/2,n-i+1/2,L1[i]]):
  fi:
od:

for i from 1 to n do

  if R1[i]<>0 then
     pic:=pic,textplot([n+1/2,n-i+1/2,R1[i]]):
  fi:
od:
     
display(pic):
end:

#Tav(n): inputs a positive integer n, and outputs a list of sets of length n, such that the i-th
#entry is the set of n-permutations whose length of largest  increasing sequence is i.
#Try:
#Tav(4);

Tav:=proc(n) local mu,i,T,mu1:
option remember:
mu:=permute(n):

for i from 1 to n do
 T[i]:={}:
od:

for mu1 in mu do
 T[LIS(mu1)]:=T[LIS(mu1)] union {mu1}:
od:

[seq(T[i],i=1..n)]:

end:


#Tavla(n): inputs a positive integer n, and outputs a list of (length n) of lists of sets of length n, such that the i-th
#j-th
#entry is the set of n-permutations whose length of largest  increasing sequence is i and the length of the reverse is j
#Tavla(4);

Tavla:=proc(n) local mu,i,j,T,mu1:
option remember:
mu:=permute(n):

for i from 1 to n do
for j from 1 to n do
 T[i,j]:={}:
od:
od:

for mu1 in mu do
 T[LIS(mu1),LIS(Rev(mu1))]:=  T[LIS(mu1),LIS(Rev(mu1))] union {mu1}:
od:

[seq([seq(T[i,j],j=1..n)],i=1..n)]:

end:

#IsLat(M): inputs  a list of lists and outputs true if it is a Latin rectangle.
#Try:
#IsLat([[1,2,3],[2,3,1]]);
IsLat:=proc(M) local i,j,n,k:
k:=nops(M):
n:=nops(M[1]):

for j from 1 to n do
 if nops({seq(M[i][j],i=1..k)})<k then
  RETURN(false):
 fi:
od:
true:
end:



#Banim(M,a): given a Latin rectangle (of size k by n say) , and an integer a, returns all the legal extensions with one extra row
#such that it is still a Latin rectangle, and if a=0 then the new row can be anything legal, but
#if a is between 1 and n then the new (k+1)-th row has length of the longest increasing subsequence equal  to a.
#Try:
#Banim([[1,2,3]],2);
Banim:=proc(M,a) local n,gu,mu,mu1,M1:

n:=nops(M[1]):

gu:={}:
if a=0 then
 mu:=permute(n):
else
 mu:=Tav(n)[a]:
fi:

for mu1 in mu do
 M1:=[op(M),mu1]:
 if IsLat(M1) then
  gu:=gu union {M1}:
 fi:
od:

gu:

end:


#Ptor1(a): inputs a list of length n, say, whose entries are integers from 0 to n outputs the
#set of Latin n by n squares, such that if a[i] is not 0, then the lenght of the longest increasing
#sequence of the i-th row is a[i]. For example, to solve the Skyscrapers puzzle in the New York Times
#magazine of June 12, 2016, type;
#Ptor1([1,4,3,0]);
Ptor1:=proc(a) local n,gu,gu1,i:
n:=nops(a):

if a[1]=0 then
gu:=permute(n):
gu:={seq([gu1],gu1 in gu)}:
else
gu:=Tav(n)[a[1]]:
gu:={seq([gu1],gu1 in gu)}:
fi:

for i from 2 to n do
 gu:={seq(op(Banim(gu1,a[i])),gu1 in gu)}:
od:
gu:
end:





#Kids(M,a,b): given a Latin rectangle (of size k by n say) , and  integers a, and b returns all the legal extensions with one extra row
#such that it is still a Latin rectangle, and if a=0 and b=0 then the new row can be anything legal, but
#if a is between 1 and n then the new (k+1)-th row has length of the longest increasing subsequence equal  to a.
#and if b is between 1 and n the revese of the new (k+1)-th row has length of the longest increasing subsequence equal  to b
#Try:
#Kids([[1,2,3,4]],2,3);
Kids:=proc(M,a,b) local n,gu,mu,mu1,M1:

n:=nops(M[1]):

gu:={}:

if a=0 and b=0 then
 mu:=permute(n):
elif a>0 and b=0 then
 mu:=Tav(n)[a]:
elif a=0 and b>0 then
  mu:=Tav(n)[b]:
  mu:={seq(Rev(mu1),mu1 in mu)}:
else
 mu:=Tavla(n)[a][b]
fi:

for mu1 in mu do
 M1:=[op(M),mu1]:
 if IsLat(M1) then
  gu:=gu union {M1}:
 fi:
od:

gu:

end:

#IsGood(M,U1,D1): Is matrix M compatible with U1 and D1?
IsGood:=proc(M,U1,D1) local i,j,pi,n:
n:=nops(M[1]):

for j from 1 to n do
 pi:=[seq(M[i][j],i=1..n)]:

 if U1[j]<>0 and LIS(pi)<>U1[j] then
   RETURN(false):
 fi:

 if D1[j]<>0 and LIS(Rev(pi))<>D1[j] then
   RETURN(false):
 fi:
od:
true:
end:


#Ptor(P): inputs a list of length 4, P=[U,D,L,R], each of them a list length n, say, whose entries are integers from 0 to n outputs the
#set of Latin n by n squares, 
#such that if 
#U[i] is not 0, then the lenght of the longest increasing sequence of the i-th column is a[i]. 
#D[i] is not 0, then the lenght of the longest increasing sequence of the reverse of the i-th column is a[i]. 
#U[i] is not 0, then the lenght of the longest increasing sequence of the i-th column is a[i]. 
#L[i] is not 0, then the lenght of the longest increasing sequence of  the i-th row is a[i]. 
#R[i] is not 0, then the lenght of the longest increasing sequence of  the reverse of the i-th row is a[i]. 
#
#For example, to solve the Skyscrapers puzzle in the New York Times
#magazine of June 19, 2016, type;
#Ptor([[0,3,0,0],[0,0,2,0],[0,3,2,0],[0,2,3,0]]);
Ptor:=proc(P) local U1,D1,L1,R1,gu1,mu,gu,i,n:
if nops({seq(nops(P[i]),i=1..nops(P))})<>1 then
 RETURN(FAIL):
fi:
U1:=P[1]:
D1:=P[2]:
L1:=P[3]:
R1:=P[4]:


n:=nops(U1):

if L1[1]=0 and R1[1]=0 then
gu:=permute(n):
elif L1[1]=0 and R1[1]>0 then
gu:=Tav(n)[R1[1]]:
  gu:={seq(Rev(gu1),gu1 in gu)}:
elif R1[1]=0 and L1[1]>0 then
gu:=Tav(n)[L1[1]]:

 else
gu:=Tavla(n)[L1[1]][R1[1]]:
fi:

gu:={seq([gu1],gu1 in gu)}:

for i from 2 to n do
 gu:={seq(op(Kids(gu1,L1[i],R1[i])),gu1 in gu)}:
od:

#gu:={seq([gu1],gu1 in gu)}:
mu:={}:

for gu1 in gu do

 if IsGood(gu1,U1,D1) then
  mu:=mu union {gu1}:
 fi:

od:

mu:


end:




#PTOR(P): inputs a list of length 4, P=[U,D,L,R], each of them a list length n, say, whose entries are integers from 0 to n outputs the
#unique Latin n by n squares, 
#such that if 
#U[i] is not 0, then the lenght of the longest increasing sequence of the i-th column is a[i]. 
#D[i] is not 0, then the lenght of the longest increasing sequence of the reverse of the i-th column is a[i]. 
#U[i] is not 0, then the lenght of the longest increasing sequence of the i-th column is a[i]. 
#L[i] is not 0, then the lenght of the longest increasing sequence of  the i-th row is a[i]. 
#R[i] is not 0, then the lenght of the longest increasing sequence of  the reverse of the i-th row is a[i]. 
#
#For example, to solve the Skyscrapers puzzle in the New York Times
#If there are more than one solution, or none, it returns FAIL
#magazine of June 19, 2016, type;
#PTOR([[0,3,0,0],[0,0,2,0],[0,3,2,0],[0,2,3,0]]);

PTOR:=proc(P) local gu,i:
gu:=Ptor(P):

if nops(gu)=0 then
 print(`There are  no solutions`):
RETURN(FAIL):
elif nops(gu)>1 then
 print(`There are  more than one solution, here there are`):
 for i from 1 to nops(gu) do
  print(matrix(gu[i])):
 od:
  RETURN(FAIL):
else
 matrix(gu[1]):
fi:
end:

#PTOR1(P): inputs a list of length 4, P=[U,D,L,R], each of them a list length n, say, whose entries are integers from 0 to n outputs the
#unique Latin n by n squares, 
#such that if 
#U[i] is not 0, then the lenght of the longest increasing sequence of the i-th column is a[i]. 
#D[i] is not 0, then the lenght of the longest increasing sequence of the reverse of the i-th column is a[i]. 
#U[i] is not 0, then the lenght of the longest increasing sequence of the i-th column is a[i]. 
#L[i] is not 0, then the lenght of the longest increasing sequence of  the i-th row is a[i]. 
#R[i] is not 0, then the lenght of the longest increasing sequence of  the reverse of the i-th row is a[i]. 
#
#For example, to solve the Skyscrapers puzzle in the New York Times
#If there are more than one solution, or none, it returns FAIL
#magazine of June 19, 2016, type;
#PTOR1([[0,3,0,0],[0,0,2,0],[0,3,2,0],[0,2,3,0]]);

PTOR1:=proc(P) local gu,i:
gu:=Ptor(P):

if nops(gu)=0 then
 print(`There are  no solutions`):
RETURN(FAIL):
elif nops(gu)>1 then
 print(`There are  more than one solution, here there are`):
 for i from 1 to nops(gu) do
  print(matrix(gu[i])):
 od:
  RETURN(FAIL):
else
 gu[1]:
fi:
end:

 




#LegalC(M): Given a Latin rectangle M, finds all the legal continuations. Try:
#LegalC([[1,2,3,4]]);
LegalC:=proc(M) local n,m,S,lu,lu1,i1,j,gu,gu1:
m:=nops(M):
n:=nops(M[1]):

for j from 1 to n do
 S[j]:={seq(i1,i1=1..n)} minus {seq(M[i1][j],i1=1..m)}:
od:

lu:=permute(n):

gu:={}:

for lu1 in lu do
 if {seq(member(lu1[i1],S[i1]),i1=1..n)}={true} then
  gu:=gu union {lu1}:
 fi:
od:

{seq([op(M),gu1],gu1 in gu)}:

end:


#LegalCs(S): Given a set of Latin rectangles S, finds all the legal continuations. Try:
#LegalCs({[[1,2,3,4]]});
LegalCs:=proc(S) local M:
{seq(op(LegalC(M)), M in S)}:
end:

#ReLS(n): the set of all n by n Latin Squares whose first row is 1...n. Try:
#ReLS(4);
ReLS:=proc(n) local gu,i:
gu:={[[seq(i,i=1..n)]]}:

for i from 2 to n do
 gu:=LegalCs(gu):
od:

gu:

end:

#RaLS(n): A random n by n Latin Square. Try:
#RalLS(5);
RaLS:=proc(n) local gu,L,i:
gu:=permute(n):

L:=[gu[rand(1..nops(gu))()]]:

for i from 2 to n do
 gu:=LegalC(L):
 L:=gu[rand(1..nops(gu))()]:
od:

L:

end:

 
#MakePuzzle1Old(n): makes a random puzzle. Try:
#MakePuzzle1Old(4);
MakePuzzle1Old:=proc(n) local P,U,D,L,R,M,i,j,U1,D1,L1,R1:
M:=RaLS(n):

L:=[seq(LIS(M[i]),i=1..n)]:
R:=[seq(LIS(Rev(M[i])),i=1..n)]:
U:= [seq(LIS([seq(M[i][j],i=1..n)]),j=1..n)]:
D:= [seq(LIS(Rev([seq(M[i][j],i=1..n)])),j=1..n)]:

P:=[U,D,L,R]:

if nops(Ptor(P))<>1 then
  RETURN(FAIL):
fi:

for i from n to 1 by -1  do
U1:=[op(1..i-1,U),0,op(i+1..n,U)]:

if nops(Ptor([U1,D,L,R]))=1 then
 P:=[U1,D,L,R]:
 U:=U1:
fi:
od:

for i from n to 1 by -1  do
L1:=[op(1..i-1,L),0,op(i+1..n,L)]:

if nops(Ptor([U,D,L1,R]))=1 then

 P:=[U,D,L1,R]:
 L:=L1:
fi:
od:


for i from n to 1 by -1 do
D1:=[op(1..i-1,D),0,op(i+1..n,D)]:

if nops(Ptor([U,D1,L,R]))=1 then
 P:=[U,D1,L,R]:
 D:=D1:
fi:
od:


for i from n to 1 by -1 do
R1:=[op(1..i-1,R),0,op(i+1..n,R)]:

if nops(Ptor([U,D,L,R1]))=1 then
 P:=[U,D,L,R1]:
 R:=R1:
fi:
od:

P:

end:


#MakePuzzleS(n): makes a random puzzle. Try:
#MakePuzzleS(4);
MakePuzzleS:=proc(n) local U,D,L,R,M,i,j:
M:=RaLS(n):

L:=[seq(LIS(M[i]),i=1..n)]:
R:=[seq(LIS(Rev(M[i])),i=1..n)]:
U:= [seq(LIS([seq(M[i][j],i=1..n)]),j=1..n)]:
D:= [seq(LIS(Rev([seq(M[i][j],i=1..n)])),j=1..n)]:

[U,D,L,R]:

end:


#DrawPS(P,S): draws a puzzle P and its solution
DrawPS:=proc(P,S) local U1,D1,L1,R1,n,pic,i,j:
U1:=P[1]: D1:=P[2]: L1:=P[3]: R1:=P[4]:
n:=nops(U1):

pic:=plot([[0,0],[n,0]],axes=none,scaling=CONSTRAINED,color=blue):
#pic:=plot([[0,0],[n,0]],axes=none,color=blue):

for i from 1 to n do
pic:=pic, plot([[0,i],[n,i]],axes=none,scaling=CONSTRAINED,color=blue):
od:

for i from 0 to n do
pic:=pic, plot([[i,0],[i,n]],axes=none,scaling=CONSTRAINED,color=blue):
od:

for i from 1 to n do

  if D1[i]<>0 then
     pic:=pic,textplot([i-1/2,-1/4,D1[i]]):
  fi:
od:

for i from 1 to n do

  if U1[i]<>0 then
     pic:=pic,textplot([i-1/2,n+1/4,U1[i]]):
  fi:
od:

for i from 1 to n do

  if L1[i]<>0 then
     pic:=pic,textplot([-1/2,n-i+1/2,L1[i]]):
  fi:
od:

for i from 1 to n do

  if R1[i]<>0 then
     pic:=pic,textplot([n+1/2,n-i+1/2,R1[i]]):
  fi:
od:

for i  from 1 to n do
 for j from 1 to n do
          pic:=pic,textplot([n+1-j-1/2,n+1-i-1/2,S[i][n+1-j]]):
 od:
od:

display(pic):
end:

#MakePuzzleOld(n,K): makes a random puzzle, by trying K times
#MakePuzzleOld(4,20);
MakePuzzleOld:=proc(n,K) local gu,i:

for i from 1 to K do
 gu:=MakePuzzle1Old(n):
 if gu<>FAIL then
   RETURN(gu):
 fi:
od:

FAIL:

end:


KickZ:=proc(L) local i,gu:
gu:=[]:
for i from 1 to nops(L) do
 if L[i]<>0 then
  gu:=[op(gu),L[i]]:
 fi:
od:
gu:
end:

#IsLegal(M): inputs a list of lists M,representing a matrix with non-negative integers checks that every
#row and every column contains distinct positive integers. Try
#IsLegal([[0,1,2],[1,0,3],[2,1,0]]);
IsLegal:=proc(M) local m,n,i,j,gu,i1:
m:=nops(M):
n:=nops(M[1]):

 for i from 1 to m do
  gu:=KickZ(M[i]):
  if nops(convert(gu,set))<>nops(gu) then
   RETURN(false):
  fi:
 od:

 for j from 1 to n do
    gu:=[seq(M[i1][j],i1=1..m)]:
  gu:=KickZ(gu):
  if nops(convert(gu,set))<>nops(gu) then
   RETURN(false):
  fi:
 od:

true:

end:


#Matim(Le,L): inputs two lists Le,L, of the same size that outputs true if every non-zero entry of Le equals the corresponding
#entry of L. Try:
#Matim([0,3,0,2],[1,3,4,2]);
Matim:=proc(Le,L) local i:
if nops(L)<>nops(Le) then RETURN(FAIL): fi:

for i from 1 to nops(Le) do
  if Le[i]<>0 and Le[i]<>L[i] then
   RETURN(false):
  fi:
od:
true:
end:

#IsFilled(L): is the list L zero-free?
IsFilled:=proc(L) local i,x:
if coeff(add(x[L[i]],i=1..nops(L)),x[0])=0 then
 true:
else
 false:
fi:

end:



#YelR(P,M,i): inputs a puzzle P=[U,D,L,R], and a partial solution M, and a row i from 1 to nops(P) outputs
#all its legal children. Try
#YelR([[0,3,0,0],[0,0,2,0],[0,3,2,0],[0,2,3,0]],[[0$4]$4],2);`):
YelR:=proc(P,M,i) local a,b,n,mu,mu1,gu,resh,lu,lu1,gu1:

n:=nops(M):

resh:=M[i]:

if IsFilled(resh) then
  RETURN(FAIL):
fi:


a:=P[3][i]: b:=P[4][i]:



if a=0 and b=0 then
 mu:=permute(n):
elif a>0 and b=0 then
 mu:=Tav(n)[a]:
elif a=0 and b>0 then
  mu:=Tav(n)[b]:
  mu:={seq(Rev(mu1),mu1 in mu)}:
else
 mu:=Tavla(n)[a][b]
fi:

lu:={}:

for mu1 in mu do
 if Matim(resh,mu1) then
  lu:=lu union {mu1}:
 fi:
od:

gu:={}:

for lu1 in lu do
 gu1:=[op(1..i-1,M),lu1,op(i+1..nops(M),M)]:
 if IsLegal(gu1) then
    gu:=gu union {gu1}:
 fi:
od:

gu:

end:


#YelC(P,M,j): inputs a puzzle P=[U,D,L,R], and a partial solution M, and a column j from 1 to nops(P[1]) outputs
#all its legal children. Try
#YelC([[0,3,0,0],[0,0,2,0],[0,3,2,0],[0,2,3,0]],[[0$4]$4],2);`):
YelC:=proc(P,M,j) local i1,a,b,n,mu,mu1,gu,resh,lu,lu1,gu1,m,i:

n:=nops(M):

m:=nops(M[1]):

resh:=[seq(M[i1][j],i1=1..m)]:

if IsFilled(resh) then
  RETURN(FAIL):
fi:

a:=P[1][j]: b:=P[2][j]:



if a=0 and b=0 then
 mu:=permute(n):
elif a>0 and b=0 then
 mu:=Tav(n)[a]:
elif a=0 and b>0 then
  mu:=Tav(n)[b]:
  mu:={seq(Rev(mu1),mu1 in mu)}:
else
 mu:=Tavla(n)[a][b]
fi:

lu:={}:

for mu1 in mu do
 if Matim(resh,mu1) then
  lu:=lu union {mu1}:
 fi:
od:

gu:={}:

for lu1 in lu do
  gu1:=[seq([op(1..j-1,M[i]),lu1[i],op(j+1..nops(M[i]),M[i])], i=1..nops(M))]:

 if IsLegal(gu1) then
    gu:=gu union {gu1}:
 fi:
od:

gu:

end:



#Yel(P,M): inputs a puzzle P=[U,D,L,R], and a partial solution M, and a column j from 1 to nops(P[1]) outputs
#the children from the best row and column. Try
#Yel([[0,3,0,0],[0,0,2,0],[0,3,2,0],[0,2,3,0]],[[0$4]$4]);`):
Yel:=proc(P,M) local i,j,mua,aluf,i1:

for i from 1 to nops(M) while IsFilled(M[i]) do od:

 if i=nops(M)+1 then
  RETURN(FAIL):
fi:

aluf:=YelR(P,M,i):

for i1 from i+1 to nops(M) do
 mua:=YelR(P,M,i1):
  if mua<>FAIL and nops(mua)<nops(aluf) then
     aluf:=mua:
  fi:
od:


for j from 1 to nops(M[1]) do
 if not IsFilled([seq(M[i1][j],i1=1..nops(M))]) then
  mua:=YelC(P,M,j):
  if mua<>FAIL and nops(mua)<nops(aluf) then
     aluf:=mua:
  fi:
 fi:
od:

aluf:

end:


#YelS(P,SetM): inputs a set of partialy-filled puzzles, outputs all their chidren of those not yet fully-filled
YelS:=proc(P,SetM) local gu,M,i:
gu:={}:

for M in SetM do
 if {seq(IsFilled(M[i]),i=1..nops(M))}={true} then
   gu:=gu union {M}:
 else
  gu:=gu union  Yel(P,M):
 fi:
od:

end:



#Sol(P): inputs a list of length 4, P=[U,D,L,R], each of them a list length n, say, whose entries are integers from 0 to n outputs the
#set of Latin n by n squares, 
#such that if 
#U[i] is not 0, then the lenght of the longest increasing sequence of the i-th column is a[i]. 
#D[i] is not 0, then the lenght of the longest increasing sequence of the reverse of the i-th column is a[i]. 
#U[i] is not 0, then the lenght of the longest increasing sequence of the i-th column is a[i]. 
#L[i] is not 0, then the lenght of the longest increasing sequence of  the i-th row is a[i]. 
#R[i] is not 0, then the lenght of the longest increasing sequence of  the reverse of the i-th row is a[i]. 
#
#For example, to solve the Skyscrapers puzzle in the New York Times
#magazine of June 19, 2016, type;
#Sol([[0,3,0,0],[0,0,2,0],[0,3,2,0],[0,2,3,0]]);
Sol:=proc(P) local n,gu,gu1,gu2:
n:=nops(P[1]):

gu:={[[0$n]$n]}:

gu1:=YelS(P,gu):

while gu1<>gu do

gu2:=YelS(P,gu1):

gu:=gu1:
gu1:=gu2:
od:

gu:

end:






#MakePuzzle1(n): makes a random puzzle. Try:
#MakePuzzle1(4);
MakePuzzle1:=proc(n) local P,U,D,L,R,M,i,j,U1,D1,L1,R1:
M:=RaLS(n):

L:=[seq(LIS(M[i]),i=1..n)]:
R:=[seq(LIS(Rev(M[i])),i=1..n)]:
U:= [seq(LIS([seq(M[i][j],i=1..n)]),j=1..n)]:
D:= [seq(LIS(Rev([seq(M[i][j],i=1..n)])),j=1..n)]:

P:=[U,D,L,R]:

if nops(Sol(P))<>1 then
  RETURN(FAIL):
fi:

for i from n to 1 by -1  do
U1:=[op(1..i-1,U),0,op(i+1..n,U)]:

if nops(Sol([U1,D,L,R]))=1 then
 P:=[U1,D,L,R]:
 U:=U1:
fi:
od:

for i from n to 1 by -1  do
L1:=[op(1..i-1,L),0,op(i+1..n,L)]:

if nops(Sol([U,D,L1,R]))=1 then

 P:=[U,D,L1,R]:
 L:=L1:
fi:
od:


for i from n to 1 by -1 do
D1:=[op(1..i-1,D),0,op(i+1..n,D)]:

if nops(Sol([U,D1,L,R]))=1 then
 P:=[U,D1,L,R]:
 D:=D1:
fi:
od:


for i from n to 1 by -1 do
R1:=[op(1..i-1,R),0,op(i+1..n,R)]:

if nops(Sol([U,D,L,R1]))=1 then
 P:=[U,D,L,R1]:
 R:=R1:
fi:
od:

P:

end:

#MakePuzzle(n,K): makes a random puzzle, by trying K times
#MakePuzzle(4,20);
MakePuzzle:=proc(n,K) local gu,i:

for i from 1 to K do
 gu:=MakePuzzle1(n):
 if gu<>FAIL then
   RETURN(gu):
 fi:
od:

FAIL:

end:


#SOLVE(P): inputs a list of length 4, P=[U,D,L,R], each of them a list length n, say, whose entries are integers from 0 to n outputs the
#unique Latin n by n squares, 
#such that if 
#U[i] is not 0, then the lenght of the longest increasing sequence of the i-th column is a[i]. 
#D[i] is not 0, then the lenght of the longest increasing sequence of the reverse of the i-th column is a[i]. 
#U[i] is not 0, then the lenght of the longest increasing sequence of the i-th column is a[i]. 
#L[i] is not 0, then the lenght of the longest increasing sequence of  the i-th row is a[i]. 
#R[i] is not 0, then the lenght of the longest increasing sequence of  the reverse of the i-th row is a[i]. 
#
#For example, to solve the Skyscrapers puzzle in the New York Times
#If there are more than one solution, or none, it returns FAIL
#magazine of June 19, 2016, type;
#SOLVE([[0,3,0,0],[0,0,2,0],[0,3,2,0],[0,2,3,0]]);

SOLVE:=proc(P) local gu,i:
gu:=Sol(P):

if nops(gu)=0 then
 print(`There are  no solutions`):
RETURN(FAIL):
elif nops(gu)>1 then
 print(`There are  more than one solution, here there are`):
 for i from 1 to nops(gu) do
  print(matrix(gu[i])):
 od:
  RETURN(FAIL):
else
 matrix(gu[1]):
fi:
end:

#SOLVE1(P): inputs a list of length 4, P=[U,D,L,R], each of them a list length n, say, whose entries are integers from 0 to n outputs the
#unique Latin n by n squares, 
#such that if 
#U[i] is not 0, then the lenght of the longest increasing sequence of the i-th column is a[i]. 
#D[i] is not 0, then the lenght of the longest increasing sequence of the reverse of the i-th column is a[i]. 
#U[i] is not 0, then the lenght of the longest increasing sequence of the i-th column is a[i]. 
#L[i] is not 0, then the lenght of the longest increasing sequence of  the i-th row is a[i]. 
#R[i] is not 0, then the lenght of the longest increasing sequence of  the reverse of the i-th row is a[i]. 
#
#For example, to solve the Skyscrapers puzzle in the New York Times
#If there are more than one solution, or none, it returns FAIL
#magazine of June 19, 2016, type;
#SOLVE1([[0,3,0,0],[0,0,2,0],[0,3,2,0],[0,2,3,0]]);

SOLVE1:=proc(P) local gu,i:
gu:=Sol(P):

if nops(gu)=0 then
 print(`There are  no solutions`):
RETURN(FAIL):
elif nops(gu)>1 then
 print(`There are  more than one solution, here there are`):
 for i from 1 to nops(gu) do
  print(matrix(gu[i])):
 od:
  RETURN(FAIL):
else
 gu[1]:
fi:
end:



#MakePuzzleB(n,K1,K2): inputs a positive integer n, and positive integers K1, K2, 
#outputs a set of size K1 of pairs [P,S] of [puzzle,solution], by trying MakePuzzle(n,K2) (q.v.) as
#many times as needed. Try:
#MakePuzzleB(4,10,20);
MakePuzzleB:=proc(n,K1,K2) local gu,P,S:

gu:={}:

while nops(gu)<K1 do
P:=MakePuzzle(n,K2):
if gu<>FAIL then
 S:=SOLVE1(P):
 gu:=gu union {[P,S]}:
fi:
od:
gu:
end: