# Programme generant un tableau de reponses a partir d'un tableau de Burt
# FGC - Janvier 2007

message <- "Voulez-vous une saisie interactive ou la demo (I/D) > "

reponse <- readline(prompt = paste(message))

# Saisie  et initialisation des donnŽes
# Q designe le nombre de variables, appelees "questions"

if (reponse == "D")
    { Q <- 5
    } else
    { Q <- as.integer(readline(prompt=paste("Nombre de variables > "))) }

# Le tableau K contient le nombre de modalites de chaque question
K <- c(rep(0,Q))

if (reponse == "D")
   { K[1] <- 4
     K[2] <- 6
     K[3] <- 4
     K[4] <- 3
     K[5] <- 5
    } else
       { for (q in  1:Q) 
           { message <- paste("Nombre de modalites de la variable",q," > " )
             K[q] <- as.integer(readline (prompt = message))
           }
       }

message <- "Voulez-vous enregister les resultats au format .csv (O/N) > "
enreg   <- readline(prompt = paste(message))

if (enreg == "O")
   { prefix <- readline(prompt = paste("Quel prefixe pour les noms de fichiers ? > "))}


# Calcul du nombre total de modalites
 
# La variable KK contient le nombre total de modalites
KK <- 0
for (q in 1:Q)
    KK <- KK + K[q]
 
# Calcul du nombre NP de patrons de reponses
 
NP <- 1
for (q in 1:Q)
    NP <- NP*K[q]
 
cat(paste("Nombre de modalites de reponses :", KK,"\n"))
 
cat(paste("Nombre de patrons               :", NP,"\n"))

# On genere l'ensemble P des patrons, sous forme d'une matrice NP*KK
 
P <- matrix(c(rep (0,NP*KK)),nrow=NP, ncol=KK)
 
# On remplit par colonne
# On remplit pour la premiere variable
 
# Initialisation des variables auxilliaires
dec <- 0
nbrep <- NP
nbbloc <- 1
lgbloc <- NP
 
for (q in 1:Q)
    { nbrep <- nbrep / K[q]
        for (k in 1:nbbloc)
          { for (i in 1:K[q])
                {for (j in 1:nbrep)
                    P[lgbloc*(k-1)+(i-1)*nbrep + j,i+dec] <- 1
                }
          }
 
      dec <- dec + K[q]
      nbbloc <- nbbloc * K[q]
      lgbloc <- NP / nbbloc
     }
 
# Saisie du tableau de Burt

if (reponse == "D")
   { BB <- matrix(c(185, 0, 0, 0, 34, 90, 40, 13, 3, 5, 61, 60, 59, 5, 70, 	
       101, 14, 32, 61, 71, 21, 0,
       0, 53, 0, 0, 5, 6, 32, 2, 3, 5, 13, 18, 21, 1, 13, 33, 7, 1, 4, 8, 40, 0,
       0, 0, 50, 0, 2, 10, 34, 0, 3, 1, 5, 14, 28, 3, 15, 23, 12, 2, 2, 14, 32, 0,
       0, 0, 0, 95, 0, 2, 10, 5, 67, 11, 1, 3, 2, 89, 5, 34, 56, 0, 1, 2, 57, 35,
       34, 5, 2, 0, 41, 0, 0, 0, 0, 0, 17, 13, 9, 2, 15, 23, 3, 27, 9, 1, 4, 0, 90,
       6, 10, 2, 0, 108, 0, 0, 0, 0, 29, 33, 45, 1, 41, 61, 6, 1, 33, 57, 17, 0, 40,
       32, 34, 10, 0, 0, 116, 0, 0, 0, 23, 35, 47, 11, 37, 62, 17, 1, 10, 29, 74, 2,
       13, 2, 0, 5, 0, 0, 0, 20, 0, 0, 6, 6, 3, 5, 6, 10, 4, 4, 7, 5, 4, 0, 3, 3, 3,
       67, 0, 0, 0, 0, 76, 0, 2, 4, 4, 66, 2, 29, 45, 0, 1, 1, 50, 24, 5, 5, 1, 11,
       0, 0, 0, 0, 0, 22, 3, 4, 2, 13, 2, 6, 14, 2, 8, 2, 1, 9, 61, 13, 5, 1, 17, 29,
       23, 6, 2, 3, 80, 0, 0, 0, 30, 44, 6, 14, 26, 24, 16, 0, 60, 18, 14, 3, 13, 33,
       35, 6, 4, 4, 0, 95, 0, 0, 25, 60, 10, 11, 22, 28, 32, 2, 59, 21, 28, 2, 9, 45,
       47, 3, 4, 2, 0, 0, 110, 0, 43, 53, 14, 10, 14, 41, 45, 0, 5, 1, 3, 89, 2, 1,
       11, 5, 66, 13, 0, 0, 0, 98, 5, 34, 59, 0, 6, 2, 57, 33, 70, 13, 15, 5, 15, 41,
       37, 6, 2, 2, 30, 25, 43, 5, 103, 0, 0, 12, 26, 38, 26, 1, 101, 33, 23, 34, 23,
       61, 62, 10, 29, 6, 44, 60, 53, 34, 0, 191, 0, 20, 35, 47, 82, 7, 14, 7, 12, 56,
       3, 6, 17, 4, 45, 14, 6, 10, 14, 59, 0, 0, 89, 3, 7, 10, 42, 27, 32, 1, 2, 0, 27,
       1, 1, 4, 0, 2, 14, 11, 10, 0, 12, 20, 3, 35, 0, 0, 0, 0, 61, 4, 2, 1, 9, 33, 10,
       7, 1, 8, 26, 22, 14, 6, 26, 35, 7, 0, 68, 0, 0, 0, 71, 8, 14, 2, 1, 57, 29, 5,
       1, 2, 24, 28, 41, 2, 38, 47, 10, 0, 0, 95, 0, 0, 21, 40, 32, 57, 4, 17, 74, 4,
       50, 1, 16, 32, 45, 57, 26, 82, 42, 0, 0, 0, 150, 0, 0, 0, 0, 35, 0, 0, 2, 0, 24,
       9, 0, 2, 0, 33, 1, 7, 27, 0, 0, 0, 0, 35), nrow=22, ncol=22)
    } else
      { BB <- matrix(c(rep (0,KK*KK)),nrow=KK, ncol=KK)
        for (i in 1:KK)
            { for (j in 1:i)
                 { BB[i,j] <- as.integer(readline(prompt=paste("Effectif" , i , j, " >")))
                   BB[j,i] <- BB[i,j]}
            }
      }
 
# Calcul du nombre total NN d'individus statistiques	
 
NN <- 0
for (i in 1:K[1])
     NN <- NN + BB[i,i]
 
cat(paste("Nombre total d'individus        :", NN, "\n"))

##################################################################
# Initialisation des tableaux de solutions
# ################################################################

# Le tableau disjonctif complet Sol
 
Sol <- matrix(c(rep(0,NN*KK)),nrow=NN, ncol=KK)
 
# Le tableau NP*(KK+1) donnant l'effectif de chaque patron
# La derniere colonne est l'effectif du patron indique dans les colonnes precedentes

SolEff <- matrix(c(rep(0,NP*(KK+1))),nrow=NP, ncol=KK+1)

SolEff[,1:KK] <- P

##################################################################
# Fonctions auxilliaires
##################################################################

# Cette definition  determine le minimum non nul d'un tableau KK*KK
minnn <- function(X)
   { if (min(X) < 0 ) { return(min(X))}
       m <- 99999
       for (i in 1:KK)
           {for (j in i:KK)
              { if (X[i,j] > 0) {m <- min(m,X[i,j])}
           }}
        return(m)}
        
 # Fonction permettant de retrouver le numero d'un patron a partir de sa valeur
 
Numpat <- function (Pat)
  { for (ip in 1:NP)
       { flag <- 1
         for (jm in 1:KK)
            { if (Pat[jm] != P[ip,jm]) 
                 { flag <- 0
                   break}
            }
          if (flag == 1) {break}
         }

         return(ip)    
    }


##################################################################
# Quelques initialisations avant d'entrer dans l'algorithme
##################################################################

# La matrice TT des taux de liaison

TT <- matrix(c(rep(0,KK*KK)),nrow=KK, ncol=KK)

# Le compteur du nombre de reponses trouvees

iii <- 0
etape <- 1

# La liste des patrons exclus

exclus <- c(0)
 	

##################################################################
# L'algorithme de recherche des reponses
##################################################################


while (iii < NN) {

#  Calcul du tableau des taux de liaison
   BBs <- apply(BB,1,sum)
   BBss <- sum(BB)

   for (i in 1:KK)
       { for (j in 1:KK)
           {TT[i,j] <- BBss * BB[i,j]/BBs[i]/BBs[j] - 1
            if ((BBs[i] == 0) | (BBs[j] == 0)) { TT[i,j] <- -1 }}}

# Choix du patron de score maximum

   Patchoisi <- 0
   Scourant <- -1

   for (p in 1:NP)
      { BP <- P[p,] %*% t(P[p,])
        for (i in 1:KK) {BP[i,i] <- 0}
        BT <- BP * TT
        Sc <- minnn(BT)
        if (Sc > Scourant ) 
            { Scourant <- Sc
              Patchoisi <- p}
      } 
      
# Retrait du patron - On reste dans l'etape 1
      
   if (Patchoisi != 0)
       { BB <- BB - P[Patchoisi,] %*% t(P[Patchoisi,])
         iii <- iii + 1
         Sol[iii,] <- P[Patchoisi,]
         SolEff[Patchoisi,KK+1] <- SolEff[Patchoisi,KK+1] + 1
         etape <- 1 
         cat(".")   }

# Si on n'a pas pu trouver de patron, on passe dans l'etape 2

    if (Patchoisi == 0) 
        { etape <- 2 }

    if (etape == 2)
       { 
# Recherche du patron realisant le minimum de liens brises
         minbrise <- Q * ( Q - 1 )
         for (p in 1:NP)
             { exclu <- 0
               for (j in exclus) {if (p == j) { exclu <- 1}}
               if (exclu == 0)
                 { BP <- P[p,] %*% t(P[p,])
                   for (i in 1:KK) {BP[i,i] <- 0}
                   BT <- BP * TT
                   brise <- 0
                   for (i in 1:KK)
                      { for (j in 1:KK)
                          { if (BT[i,j] == - 1) { brise <- brise + 1 } }
                      }
                   if (brise < minbrise) 
                    { minbrise <- brise
                      PatCandidat <- p}
                  }
# Si on trouve un patron avec 1 seul lien brise, on arrete la recherche
               if (minbrise == 2 ) {break}      
              }

# Trouver le lien brise (i1,j1)

      BP <- P[PatCandidat,] %*% t(P[PatCandidat,])
      for (i in 1:KK) {BP[i,i] <- 0}
      BT <- BP * TT
      for (i in 1:KK)
          { for (j in i:KK)
             { if (BT[i,j] == - 1) 
                  { i1 <- i 
                    j1 <- j } 
              }
          }
          
# Trouver les patrons Y1 et Y2 a retirer des solutions actuelles

      succes <- 0
 
      for (i in 1:iii)
          {  
            if ((Sol[i,i1]== 1) & (Sol[i,j1] == 1))
              { for (j in 1:iii) 
                 { if ((Sol[j,i1] == 0) & (Sol[j,j1] == 0) & (sum(Sol[i,]*Sol[j,])==(Q-2)))
                    { s1 <- i
                      s2 <- j
                      i2 <- 0
                      j2 <- 0
                      for (ij in 1:KK)
                         { if ((Sol[s2,ij] == 1) & (Sol[s1,ij] == 0))
                             if (i2 == 0 ) 
                                { i2 <- ij} 
                             else
                                { j2 <- ij }
                         }
                      if ((BB[i1,j2] >= 1) & (BB[j1,i2] >=1)) {succes <- 1}
                    }
# On s'arrete des que l'on a trouve un couple (Y1, Y2) convenable
                 if (succes == 1) {break}}}
           if (succes == 1 ) {break}}

# S'il n'existe pas de couple (Y1, Y2) convenable, le patron candidat est exclu

      if (succes == 0) {exclus <- c(exclus,PatCandidat)}

# Si on a trouve (Y1, Y2)
     
     if (succes == 1)
        { BB <- BB + Sol[s1,] %*% t(Sol[s1,])
          BB <- BB + Sol[s2,] %*% t(Sol[s2,])
          i2 <- 0
          j2 <- 0
          for (i in 1:KK)
             { if ((Sol[s2,i] == 1) & (Sol[s1,i] == 0))
                  if (i2 == 0 ) 
                     { i2 <- i} 
                  else
                     {j2 <- i }
             }

# On forme les patrons Z1 et Z2
         P1 <- Sol[s1,] * Sol[s2,]
         P2 <- P1
         P1[i1] <- 1
         P1[j2] <- 1
         P2[i2] <- 1
         P2[j1] <- 1

# On enleve Z1 et 22

         BB <- BB - P1 %*% t(P1)
         BB <- BB - P2 %*% t(P2)

# Mise a jour de SolEff
        pm1 <- Numpat(Sol[s1,])
        pm2 <- Numpat(Sol[s2,])
        ps1 <- Numpat(P1)
        ps2 <- Numpat(P2)
        SolEff[pm1, KK+1] <- SolEff[pm1, KK+1] - 1
        SolEff[pm2, KK+1] <- SolEff[pm2, KK+1] - 1
        SolEff[ps1, KK+1] <- SolEff[ps1, KK+1] + 1
        SolEff[ps2, KK+1] <- SolEff[ps2, KK+1] + 1


# On substitue Z1 et Z2 a Y1 et Y2 dans l'ensemble des solutions
        Sol[s1,] <- P1
        Sol[s2,] <- P2

# Si minbrise = 2, on retire aussi de B le patron candidat 
# et on l'ajoute a l'ensemble des solutions


        if (minbrise == 2 ) 
           { BB <- BB - P[PatCandidat,] %*% t(P[PatCandidat,])
             iii <- iii + 1
             Sol[iii,] <- P[PatCandidat,]
             SolEff[PatCandidat,KK+1] <- SolEff[PatCandidat,KK+1] + 1
             cat(".")   
           }
        }

# On repasse dans l'etape 1

      etape <- 1
      } # Fin de l'etape 2
     

   if (any(BB < 0)) { print("Attention BB n'est plus positif")}

  } # Fin de la boucle while (iii < NN)


# Sortie du programme. Demander a l'utilisateur s'il veut enregistrer ses resultats
# Prevoir aussi une sortie lorsque le nombre d'essais infructueux est trop grand ?? 

cat("\nExecution terminee. Les resultats sont dans les tableaux Sol et SolEff\n")

if (enreg == "O")
   { write.csv2(Sol, paste(prefix,"Sol.csv",sep=""))
     write.csv2(SolEff,paste(prefix,"SolEff.csv",sep=""))
     cat("Les resultats sont enregistres dans les fichiers ",prefix,"Sol.csv\n",sep="")
     cat("et ",prefix,"SolEff.csv\n",sep="")
   }