restart;

#############################################
# PROGRAMME POUR CALCULER LES CHAMBRES
#
# version Feb 08 2003
#############################################

with(linalg):with(combinat):with(simplex):
interface(prettyprint=3):printlevel:=0:
interface(quiet=true):
  
#############################################
# Soit A une matrice N times n.  
# Dans la procedure cocircojesus, l'input est la matrice A. 
# L'output nous donne, pour chaque mur h de A ( c'est a dire un
# hyperplan engendre par des vecteurs colonnes voir explications.tex)
# la liste des ensemble pos(h), zeros(h), neg(h).  On peut le faire
# avec chirotope pour des grandes matrices plutot qu'avec des 
# sousdeterminants.
##############################################
 
cocirco_jesus:=proc(a) 
 local n,N,T,ch,s,zeros,pos,neg,j,aux: 
 option remember:
 N:=coldim(a);n:=rowdim(a) :  
 T:=choose(N,n-1); ch[0]:={};
 for s from 1 to nops(T) do 
  zeros[s]:={};
  pos[s]:={};
  neg[s]:={}; 
  for j from 1 to N do                                                             
    if not(member(j,op(s,T))) then 
      aux[j]:=det(submatrix(a,[$1..n],[op(op(s,T)),j])) ; 
    else
      aux[j]:=0 :
    fi: 
    if aux[j]=0 then 
       zeros[s]:=zeros[s] union {j}  
    elif sign(aux[j])=1 then
           pos[s]:=pos[s] union {j}:
    elif sign(aux[j])=-1 then 
      neg[s]:=neg[s] union {j}:
         fi:
  od:
  if (pos[s]={} and neg[s]={}) then  
    ch[s]:=[op(ch[s-1])]; 
  else 
    ch[s]:=[op(ch[s-1]),[pos[s],zeros[s],neg[s]]]; 
  fi; 
 od;
 RETURN(ch[nops(T)]); 
end:
 
#################################################
#On va tout de suite calculer, grace a la matrice des cocircuits, 
#la chambre lexicographique initiale: c'est fait a l'aide de minimalize, 
#produit, polarized...
#
#
# La procedure minimalize permet de trouver les ensembles minimaux
# d'un ensemble d'ensembles. L'input est un ensemble de sous ensembles
# d'un ensemble Delta (par exemple Delta=(1,2,...,N)). L'ouput est un
# ensemble de sous-ensembles de Delta. C'est l'ensemble des
# sous-ensembles minimaux.
##################################################

sizz := proc(a,b) evalb( nops(a) > nops(b)) end:

minimalize := proc(NF) 
local erg,dada,i,j: 
erg := NF: 
dada := sort(convert(NF,list),sizz): 
for i from 1 to nops(dada)-1 do 
  for j from i+1 to nops(dada) do 
   if ((dada[j] minus dada[i]) = {}) then 
    erg := erg minus {dada[i]}: 
   fi: 
  od:
od: 
RETURN(erg);
end:
 
#################################################### 
#  La procedure produit construit a partir d'une famille d'ensembles
# E_1, E_2,...,E_N l'ensemble constitue par les ensembles ou on prend
# un point dans chaque E_i. L'input est un ensemble d'
# ensembles. L'ouput est un ensemble d'ensembles.
####################################################

produit := proc( fam ) 
local erg,i,j,head,tail:
erg := {}: 
if (fam = {}) then 
  erg := { {} }: 
fi: 
if (fam <> {}) then
  head := fam[1]:
  tail := produit( fam minus {fam[1]}):
  for i from 1 to nops(head) do
    for j from 1 to nops(tail) do
      erg := erg union {{head[i]} union tail[j]}:
    od: 
  od:
fi:
RETURN(erg):
end:
 

#########################################################
# La procedure get_a_chamber a comme input un ensemble de sous-ensembles
# de Delta. Comme output, c'est un ensemble de sous-ensembles de
# Delta. Lorque l'input est appliqué a un ensemble P obtenu en
# "polarisant" les cocircuits d'une matrice A (N fois n) de rang n par
# rapport a un élément tres regulier, l'ensemble obtenu est un sous
# ensemble de sous ensembles de (1,2,...,N) de cardinal n.
########################################################## 

get_a_chamber:= proc(polarized) 
local P,g;
P:=polarized; 
g:=minimalize(produit(minimalize(P)));
RETURN(g);
end:
 
############################### 
# La procedure polarized_lex a comme input un ensemble de cocircuits.
# Un cocircuit est une partition de l'ensemble {1,2,...,N} en 3
# ensembles disjoints c[1],c[2],c[3]. l'output est un ensemble de
# sous-ensembles de (1,2,...,N).  La procedure polarized-lex permet a
# partir d'un ensemble de cocircuits { {c[1],c[2],c[3]}} de choisir soit
# c[1] soit c[3]. On choisit celui des deux ensembles c[1] ou c[3] qui
# contient le plus petit element.  par exemple, la procedure circojesus
# appliquée à une matrice donne comme output une famille de cocircuits.
# Si on applique la procedure polarized_lex a cet ensemble de cocircuits
# cocircojesus(A), on obtient l'ensemble Polarized(A,xi) ou xi est le
# vecteur vecteur_lex=\alpha_1+\epsilon \alpha_2+\epsilon^2\alpha_3+....
############################### 
 
polarized_lex := proc(cocircuits)
local erg,mm,c:
erg := {}:
for c in cocircuits do
  mm := min(convert(c[1] union c[3],list)[]):
  if (member(mm,c[1])) then erg := erg union {c[1]}:
       else erg := erg union {c[3]}: fi:
 od:
 RETURN(erg);
 end:

#############################################
# La procedure lex_chamber produit la chambre lexicographique. Donc
# l'input est une matrice A et l'output est l'ensemble des "basic
# subsets" de Delta contenant l'élément régulier le plus près de 
# \alpha_1+\epsilon \alpha_2+\epsilon^2\alpha_3+....
#############################################

lex_chamber:=proc(A) 
 local CA,PA,c ; 
 CA:=cocirco_jesus(A); 
 PA:=polarized_lex(CA);
 c:=produit(minimalize(PA));
 RETURN(c):
end:
 
################################################## 
# Soit A une matrice N times n. Dans la procedure basis_int_cocirco,
# l'input est la matrice A. L'output nous donne
# la liste des elements de W(A) (c'est a dire des ensembles de (n-1)
# vecteurs colonnes lineairement independants voir explications.tex) et
# pour chaque nu dans W(A) l'ensemble zeros(nu), pos(nu), neg(nu).  On
# aura besoin de ceci que pour les murs interieurs.
#
# Dans l'output le premier element est une base d'un mur interieur,
# puis le mur, puis les elements positifs puis les elements negatifs.
################################################################

basis_int_cocirco:=proc(a) 
local n,N,T,ch,s,zeros,pos,neg,j,aux: 
option remember;
N:=coldim(a);n:=rowdim(a) :  
T:=choose(N,n-1); ch[0]:={};
for s from 1 to nops(T) do 
  zeros[s]:={};pos[s]:={};neg[s]:={};           
  for j from 1 to N do                                                                                  if not(member(j,op(s,T))) then 
       aux[j]:=det(submatrix(a,[$1..n],[op(op(s,T)),j])) ; 
     else
       aux[j]:=0 :
     fi: 
     if aux[j]=0 then 
       zeros[s]:=zeros[s] union {j}  
     elif sign(aux[j])=1 then
           pos[s]:=pos[s] union {j}:
     elif sign(aux[j])=-1 then 
           neg[s]:=neg[s] union {j}:
     fi:
  od:
  if (pos[s]={} or neg[s]={}) then  
     ch[s]:=[op(ch[s-1])]; 
  else ch[s]:=[op(ch[s-1]),
    [op(s,T),zeros[s],pos[s],neg[s]]]; 
  fi;
od;
RETURN(ch[nops(T)]); 
end:
 
####################################################################
# Dans flipper, l'input est une matrice A. L'output est pour chaque
# mur interieur, une base nu de ce mur et les ensembles delta^+(nu) et
# delta^-(nu). Le premier element est nu, puis zeros(nu) puis
# delta^+(nu) qui est un ensemble de simplexes, et de meme
# delta^-(nu) qui est un ensemble de simplexes.
####################################################################
 
int_flipper:=proc(a) 
local flips,cocircuitos, j,pos,zeros,neg,basishyp,deltaplus,deltaminus,i,simplexo,aux: 
option remember: 
flips:=[]:
cocircuitos:= basis_int_cocirco(a); 
for j from 1 to nops(cocircuitos) do
  pos:=op(3,op(j,cocircuitos)): 
  zeros:=op(2,op(j,cocircuitos)):
  neg:=op(4,op(j,cocircuitos)):
  pos := convert(pos,list):
  neg := convert(neg,list): 
  basishyp:=convert(op(1,op(j,cocircuitos)),set):
  deltaplus:={}:
  deltaminus:={}:
  for i from 1 to nops(pos) do
    simplexo:=basishyp union {op(i,pos)}:
    deltaplus:=deltaplus union {simplexo}:
  od:
  for i from 1 to nops(neg)  do
    simplexo:=basishyp union {op(i,neg)}:
    deltaminus := deltaminus union {simplexo}:
  od:
  aux:=[zeros,basishyp,deltaplus,deltaminus];
  flips:=[op(flips),aux]:
od:
end:
 
hyp_flipper:=proc(a,h)
local f,F,j: 
option remember; 
f:=[]: 
F:=int_flipper(a): 
for j from 1 to nops(F) do 
  if op(1,op(j,F))=h then 
    f:=[op(f),op(j,F)]: 
  fi:
od:
RETURN(f):
end: 
 
###############################
# REFLEXION a comme input une matrice A, N fois n, une chambre
# combinatoire c, c'est a dire, un ensemble de sous-ensembles de
# cardinal n de [1,2,...,N] et h qui est un sous ensemble de [1,2,...,N]
# l'output est une chambre combinatoire.
###############################

reflexion:=proc(A,c,h) local f,cc,j: 
option remember; 
cc:=c:
f:=hyp_flipper(A,h); 
for j from 1 to nops(f) do 
   if op(3,op(j,f)) intersect c= op(3,op(j,f)) then 
     cc:=(cc minus op(3,op(j,f))) union op(4,op(j,f)) 
   else if op(4,op(j,f)) intersect c=op(4,op(j,f)) then 
     cc:=(cc minus op(4,op(j,f))) union op(3,op(j,f)) 
   fi: 
  fi: 
 od: 
 cc:
end:

##############################################
# La procedure suivante recopiee de puntos decrit le cone 
# simplicial C(sigma) donne par un
# sous ensemble sigma appele simplexe. L'input est A: une matrice N
# times n et un sous ensemble sigma de n elements de (1,2,..,N).
# l'output est l'ensemble des inegalitees du cone simplicial engendre par
# les vecteurs colonnes correspondants a sigma. En fait on pousse un peu
# le cone C(sigma) vers son interieur, car on voudra utiliser la
# procedure de maple: feasible, qui peut decider si le system d'inegalitees
# a une solution ou non.
##############################################
 
equ_simplicial_cone:=proc(a,simp,LL)
local ineqs,elemento,temp,i,aux,u,facet,ecuacion,  directionineq:
global ineq_vs_point_incidence_table:

 ineqs:={}:
 for elemento in simp do
   temp:=simp minus {elemento}:
   temp:=convert(temp,list):
   aux:=[]:
   for i from 1 to rowdim(a) do
    aux:=[op(aux),x[i]]:
   od:
   u:=convert(aux,vector):
   facet:=stackmatrix(u,transpose(submatrix(a,[$1..rowdim(a)],temp)));
   ecuacion:=det(facet);
   directionineq:=signum(det(stackmatrix(transpose(submatrix(a,[$1..rowdim(a)],[elemento])),transpose(submatrix(a,[$1..rowdim(a)],temp)))));
   if directionineq=-1 then
      ineqs:=ineqs union {ecuacion+1<=0}:
      ineq_vs_point_incidence_table:=ineq_vs_point_incidence_table union {[standardize({ecuacion+1<=0}),convert(temp,set)]}:
   else
      ineqs:=ineqs union {-ecuacion+1<=0}:
      ineq_vs_point_incidence_table:=ineq_vs_point_incidence_table union {[standardize({-ecuacion+1<=0}),convert(temp,set)]}:
   fi:
 od:
 ineqs;
 end:                             


######################################################
# Pour une systeme d'inegalitees X cette procedure trouve
# l'ensemble plus petit possible de inegalitees qui donne le meme
# cone defini par X. Ces sont exactement les murs essentials.
###################################################

findnonredundantequations:=proc(L)
local nonredundant,possible,value,i,new,testset,inequa,maximo;
option remember:
nonredundant:={}:
possible:=L:
for i from 1 to nops(possible) do
  new:={op(1,op(i,possible)) <= op(2,op(i,possible))+1}:
  testset:=possible minus {op(i,possible)} union new:
  inequa:=op(1,(op(i,possible))):
  maximo:=maximize(inequa,testset):
  value:=subs(op(maximo),inequa);
  if value>op(2,op(i,possible)) then
    nonredundant:=nonredundant union {op(i,possible)}:
  fi:
od:
RETURN(nonredundant):
end:

########################################################
# Cette procedure trouve exactement quelle sont les 
# murs essentials. L'input est une matrice A et une
# chambre. L'output est un ensemble d'elements de W(A) que
# engendrent des exactement les murs pour une chambre c.
# L'idee cest que le murs essentiales peut pas etre efface
# sans changer la region.
########################################################## 
 
here_are_the_essential_walls:=proc(A,LL)
local cell,ecuaciones,essential_walls_eqs,combinatorial_walls,aux,c,i,still_not_found,K,L,j:
 global ineq_vs_point_incidence_table:
 option remember:
 
 ecuaciones:={}:
 ineq_vs_point_incidence_table:={}: 

# This set keeps track of which points satisfy with a equality an inequality.
# Next we find the inequalities:

 for cell in LL do
   ecuaciones:=ecuaciones union equ_simplicial_cone(A,cell,LL);
 od:
 ecuaciones:=standardize(ecuaciones):

# Find which are the irredundant equations:

 essential_walls_eqs:=findnonredundantequations(ecuaciones); 

# From inequalities we want to get the zero sets of the wall:

 combinatorial_walls:={}:
 for c in essential_walls_eqs do
  i:=1:
  still_not_found:=true:
  while still_not_found do 
    aux:={c} minus op(1,op(i,ineq_vs_point_incidence_table)):
    if aux ={ } then
      combinatorial_walls:=combinatorial_walls union {op(2,op(i,ineq_vs_point_incidence_table))}:
    still_not_found:=false:
    else
      i:=i+1:
    fi:
  od:
 od: 
RETURN(from_wallbase_to_wallzeros(A,combinatorial_walls));
end: 

from_wallbase_to_wallzeros:=proc(A,combinatorial_walls)
local f,F,h,j;
option remember:
  f:={}:
  F:=int_flipper(A): 
  for h in combinatorial_walls do
   for j from 1 to nops(F) do 
    if op(2,op(j,F))=h then 
     f:=f union {op(1,op(j,F))}: 
    fi:
   od:
 od:
RETURN(f):
end:

 
###############################
# Ici, on calcule les voisins d'une chambre.
###############################
 
 voisins:=proc(A,c) 
  local v,s,r,l,j; 
  option remember: 
 
  v:={};s:=1;r[0]:=c;  
 st:=time():   
  l:=here_are_the_essential_walls(A,r[0]); 
  for j from 1 to nops(l) do  
    r[j]:=reflexion(A,r[0],l[j]); 
     v:=v union {r[j]};
  od;
 RETURN(v);
 end:

  
 get_more_chambers:=proc(A,v) 
 local W,j,V,st: 
 option remember:

 W[0]:=v: 
 for j from 1 to nops(v) do   
   V[j]:=voisins(A,op(j,v)): 
   W[j]:=W[j-1] union V[j]: 
 od:
 RETURN(W[nops(v)]):
 end:
 
 
total_chambers:=proc(A,c) 
  local old, new,tot,i; 
  option remember: 
  old:={}: new:={c}:
  print(c);
  tot:={c}: i:=1:   
  while i>0 do 
    old:=tot: 
    new:=get_more_chambers(A,new) minus old: 
    tot:=old union new: 
    i:=nops(new):  
  od: 
  RETURN(tot):
end:

#############################
# le program principal. On commence avec le matrice
# A et on calcul les chambres du cone (A).
############################# 
 
Chambers:=proc(A) 
local c,C:
 option remember:
 print(`on commence le calcul, voici la premier chambre`); 
 c:=lex_chamber(A):
 C:=total_chambers(A,c):
 print(`fini`);
 RETURN(C):
 end:
 

get_chamber_containing_vector:=proc(a,b)
local n,N,T,bdryfacets_cone,s,simplices_containing_vec,aux: 
N:=coldim(a);n:=rowdim(a) :  
T:=choose(N,n);
aux:=[]:
simplices_containing_vec:={}: 
 for s from 1 to nops(T) do
  if det(submatrix(a,[$1..n],[op(op(s,T))]))<> 0 then
    bdryfacets_cone:=equ_simplicial_cone(a,{op(op(s,T))});
    for i from 1 to rowdim(a) do
     aux:=[op(aux),x[i]=op(i,b)]:
    od:
    if map(z->evalb(z), subs(aux,bdryfacets_cone))={true} then 
      simplices_containing_vec:=simplices_containing_vec union {op(s,T)}:
    fi:
  fi:
 od;
 RETURN(simplices_containing_vec)
end:


#------------EXAMPLES---------------------------- 
#
#moco:=get_chamber_containing_vector(A9,[2,3,7,17]);
 
A9:= matrix([[1, 0, 0, 0, 1, 0, 1,0, 0],
             [0, 1, 0, 0, 1, 1, 1,0, 1], 
             [0, 0, 1, 0, 0, 1, 1,1, 1],
             [0, 0, 0, 1, 0, 0, 0,1, 1]]);
 

nops(Chambers(A9));

#U:=transpose(matrix([[1, -1, 0], [1, 0, -1], [1, 0, 0], [0, 1, -1], [0, 1, 0], [0, 0, 1]]));

#nops(Chamber(U));

#V:=
#matrix([[1, -1, 0, 0], [1, 0, -1, 0], [1, 0, 0, -1], [1,
#0, 0, 0], [0, 1, -1, 0], [0, 1, 0, -1], [0, 1, 0, 0], [0,
#0, 1, -1], [0, 0, 1, 0], [0, 0, 0, 1]]);
#nops(Chambers(transpose(V)));