martes, 25 de septiembre de 2018

A Fast Functional Algorithm to Generate Random Latin Squares

Introduction


In this blog post, I will describe briefly an algorithm I have developed in the Scala language, which is useful to generate Random Latin Squares (LSs). In fact, it is the functional, stateless and recursive version of another Java algorithm that I had done as part of my master thesis, back in 2015.
 I called this algorithm "the SeqGen (sequential-generation type) with a replacement chain conflict resolution strategy" or more succinctly "SeqGen with replacement chain", or just "the replacement chain algorithm". Its main advantage is that is faster than (by a factor of around 2) the Jacobson & Matthews' (J&M) algorithm, the most widely known and accepted algorithm in the area. On the other hand, its main disadvantage is that it is not uniformly distributed (as the J&M's), but its distribution is acceptable under certain conditions, as it is approximately uniform.


LS Definition


A Latin Square (LS) of order n, is a matrix of n x n symbols (to simplify just numbers from 1 to n), in which every symbol appears exactly once in each row and column. For example:

1 2 4 3 
2 4 3 1
4 3 1 2
3 1 2 4

is a random LS of order 4.

LSs are used for many applications, for example for games, combinatorics and some cryptographic applications.

A remarkable property is that L(n), the number of different LSs of order n, is only known when n<=11. For larger orders, the number L(n) is unknown up to now, and it is so high that only upper and lower bounds can be given:



The "random"-LS part means that the numbers are "mixed" as much as possible along the matrix as if they were a face of a dice in each position. An algorithm is "uniformly distributed" if every possible LS has the same probability of being generated by the algorithm.


The Scala algorithm


We will explain a generator class, where the algorithm is divided into different methods. 

We start with an empty matrix, a matrix full of 0s or nulls, as you prefer. We begin by row with index 0, extracting a random number one at a time. When we finish the row 0, we can continue with row 1, in order. As the matrix is immutable, we need a function that takes a partially completed LS with n rows filled and a row index and returns a new LS with n+1 rows completed. There is a theorem that says that this can always be done, row by row, up to the last row.

How can we model a computation like this in FP? This was the first "eureka" moment: using the foldlLeft function and the empty matrix as the base element. This is the "main" method to generate an LS:


override def generateLS(): LatinSquare[Int] = {
  val ls0 = VectorLatinSquare.getEmptyMatrix(order)

  val rows = (0 to (order - 1)).toList
  rows.foldLeft (ls0) { (lsi, i) => this.generateRow(i, lsi) }
}


With the 0-th index, beginning with the empty matrix ls0, we generate row 0. Similarly, with the index i-th, and an LS filled up to the row i-1, we generate row i.

The getEmptyMatrix auxiliary method is simply:

def getEmptyMatrix(order: Int) = {
  val vector = Vector.tabulate(order, order)
                       { (i, j) => RandomUtils.getNullElem() }
  new VectorLatinSquare(vector)

}

Where we generate a Vector full of a special element called "null" element. Then we create and return a vector-based implementation of an LS.

We continue with the generateRow method:

def generateRow(i: Int, partialLatinSquare: VectorLatinSquare[Int]): VectorLatinSquare[Int] = {
  val emptyRow = Vector()
  val row = this.generateRowWithPartialRow(i, 0, partialLatinSquare, emptyRow)
  partialLatinSquare.setRow(i, row)

}

This method is just the wrapper and initialization for a recursive method:

protected def generateRowWithPartialRow(i: Int, j: Int, partialLS: ImmutableLatinSquare[Int], partialRow: Vector[Int]): Vector[Int] = {

  if (partialRow.size == partialLS.order) {
    partialRow
  } else {
    val elem = this.generateElem(i, j, partialLS, partialRow)
    if (!partialRow.contains(elem)) {
      this.generateRowWithPartialRow(i, j + 1, partialLS, partialRow :+ elem)
    } else { //conflict: elem is in partialRow
      val emptyPath = Vector()
      val partialRowWithElem = partialRow :+ elem
      val availInRow = partialLS.allNumsVector.diff(partialRow)
      val idxOfFirstElem = partialRow.indexOf(elem) //search for the root of the problem
      val newPartialRow = ControlStructures.retryUntilSucceed(
                            this.makeRoomForElem(elem, idxOfFirstElem, partialLS, partialRowWithElem, availInRow, emptyPath)
                               )
      //now elem fits in position (i,j), continue with column (j + 1)
      this.generateRowWithPartialRow(i, j + 1, partialLS, newPartialRow)
    }//conflict
  }//not the last elem

}//method

Here we have a lot to explain. But let us explain by cases. The first condition checks if the row is completed (the base case), that is, if its length is equal to the order of the LS. If this condition holds, we return the completely generated row.

Then, if the row is incomplete, we try to generate the element in position (i,j) randomly, possibly causing a conflict (repetition) in that position with previously generated elements in the same row. If a conflict is created (there is not another possibility in that position) we trigger the resolution strategy called "replacement chain". Basically, we make room for the element elem, taking into account, that this chain of replacements can fail. So, we repeat the procedure until we can accommodate the element elem without repetitions. In that case, we can continue the generation of the element in position (i, j + 1) recursively.

We highlight here the creation of a new control structure retryUntilSucceed that repeats an operation in case of an exception, extending in a way the Scala language:

import util._

object ControlStructures {

  def retryUntilSucceed[T](fn: => T): T = {
    util.Try { fn } match {
      case Success(x) => x
      case Failure(e) => retryUntilSucceed(fn)
    }
  }
... }

Just for the sake of completeness, we list here the generateElem method, using set operations like diff (for set subtract), intersect (for set intersection). It returns a new random element that can be positioned at the position (i, j), using the java.SecureRandom class, that is secure for cryptographic applications:

private def generateElem(i: Int, j: Int, partialLS: ImmutableLatinSquare[Int], partialRow: Vector[Int]): Int = {

  val availInRow = partialLS.allNumsVector.diff(partialRow)
  val availInCol = partialLS.availInCol(j)
  val availInPos = availInRow.intersect(availInCol)
  if (availInPos.isEmpty) {
    RandomUtils.randomChoice(availInCol) //there is a conflict here, allow a repetition in row but never in column j
  } else {
    RandomUtils.randomChoice(availInPos)
  }

}

Last, here is the algorithm that makes the replacement chain, that is, the sequence of replacements to make room for the selected element:

private def makeRoomForElem(elem: Int, prevIdx: Int, partialLS: ImmutableLatinSquare[Int], currentRow: Vector[Int], availInRow: Vector[Int], path: Vector[Int]): Vector[Int] = {

  //calculate available to choose at idxOld
  val availInCol = partialLS.availInCol(prevIdx)
  val available = (availInCol).diff(path).diff(Vector(elem)) //must be in available in col but not in path, and must not use elem
  
  if (available.isEmpty) {
    // Path no good, begin again
    throw new Exception("No good path... failure")
  }

  //choose a new element for inxOld position
  val newElem = RandomUtils.randomChoice(available);
  val idxNew = currentRow.indexOf(newElem); //index of this newElem before replacement because it will be repeated
  //make the replacement:
  val newRow = currentRow.updated(prevIdx, newElem) //replace and generate a repetition of elem newElem
  //store in path
  val newPath: Vector[Int] = path :+ newElem
  //remove from available in row
  val newAvailInRow = availInRow.diff(Vector(newElem))
  if (idxNew == -1) { //elem is now available in row and there are no repetitions
    return newRow
  } else {
    return this.makeRoomForElem(elem, idxNew, partialLS, newRow, newAvailInRow, newPath)
  }

}

In each recursive call, we take the previously generated and repeated element, and replace it by some other possible element in that position, taking care of not repeating elements in the column of that element. The base case of the recursion is when the updated row does not contain the newElem element, and the generation can proceed from that point.


Conclusion


We have seen a complete algorithm to generate Random LSs. The algorithm is very efficient and behaves better than the most widely known algorithm of J&M's. Although the Java iterative version of the replacement chain algorithm is more efficient, the Scala recursive version is more clear and less error-prone, composable, more robust and easier to test.

Ignacio Gallego Sagastume
Blue Montag Software
@igallegosagas

miércoles, 13 de diciembre de 2017

First complete Scala application

Hello from the world of the Scala language...

I've just reengineered a Java project from my Master thesis into a new project in Scala. The code can be found at:

https://github.com/bluemontag/ScalaLSGP

I've called the new project ScalaLSGP (Scala Latin Square Generation Package)

To test the code, just run the Scala tests in the project to see the Jacobson & Matthews' algorithm in action, and the OpenGL graphics in a new Window.

There are also some more tests to see the other algorithms in action.

Let me know if you have any problems. See you soon!

sábado, 1 de julio de 2017

Sobre Blue Montag Software

Bienvenido al blog de Blue Montag Software!

  En el año 2015, luego de terminar mi Magister en Ing. de Software, nació la idea de este proyecto, con el fin de seguir investigando en temas de ciencias de la computación que considero interesantes. Esta no es mi principal actividad, aunque pretende ser más que un hobby, con el desarrollo de productos útiles para la comunidad científica e industria en general.

  Los temas de mi interés son los siguientes:
  • Algoritmos para resolver problemas de combinatoria y matemática en general
  • Generación de estructuras aleatorias mediante algoritmos
  • Estructuras algebraicas aplicables en problemas de informática
  • Computación gráfica de algoritmos (Open GL 3D)
  • Computación científica (implementación de métodos conocidos en lenguajes de programación)
  • Matemáticas para ciencias de la computación
  • Semántica de lenguajes de programación
  • Cálculos formales (lambda, sigma, etc.)
  • Programación funcional (Haskell, Clean, Scala, etc.)
  • Métodos formales (verificación de programas)
  • Graph Analytics?

  Si sos matemático o informático y tenes conocimientos sobre alguno de mis temas de interés y te interesa aportar trabajo o comentarios, no dejes de enviarme un e-mail a bluemontag@gmail.com. Si tenés algún problema relacionado y querés resolverlo con un algoritmo, también contactate!


  Si te pareció interesante o te sirvió alguno de los proyectos de Blue Montag Software, podés contribuir donando la cantidad que quieras, haciendo click en el botón "Donate" de más abajo. Con tu aporte podremos seguir trabajando y desarrollando mejores productos. Muchas gracias!

miércoles, 28 de diciembre de 2016

Drawing data structures in 3D with Haskell

In this small tutorial I will explain all I can about how to draw a Graph in 3D using OpenGL bindings for Haskell and GLUT. Additionally, I will use the System.Random module to select random 3D points in space for drawing the graph.
  1. Seting up the environment

  2. First, you have to download and install the Haskell platform from:


    For the creation and configuration of a cabal sandbox, I strongly recomend following this tutorial:


  3. Compiling and running the project

  4. Every time you change one or more files from the folder that composes the sandbox, you have to run the following command:

    cabal build

    or 

    cabal run

  5. Defining the data structure to draw

  6. First, we name the module "DataGraph":

    module DataGraph where

    Then, import the required libraries:

    import Graphics.Rendering.OpenGL
    import Graphics.UI.GLUT
    import Data.Graph
    import System.Random


    For convenience, we create the synonym type "GraphPoint", to be a triple of GLfloat (the float type for OpenGL bindings):

    type GraphPoint = (GLfloat,GLfloat,GLfloat)

    Then, we use the built-in function of module Data.Graph "buildG" to generate the graph to be drawn. The function takes an interval (n,m) of ints that will be the graph's vertex set, and a list of directed edges (v1,v2) to build the graph:

    graph1 :: Graph
    graph1 = buildG (0,10) [(1,2),
                            (2,3),
                            (3,1), (3,5), (3,4),
                            (4,1), (4,2),
                            (5,4),
                            (6,2), (6,5),
                            (7,1),
                            (8,2), (8,5),
                            (9,3),
                            (9,4), (9,5),
                            (10,9), (10,8)]



    The function points3D will take every vertex from the graph previously defined, and for each one, it will return a 3D point to be drawn:

    points3D :: Graph -> Int -> Int -> [ GraphPoint ]
    points3D graph r1 r2 = [ getPoint r1 r2 v | v <- (vertices graph) ]


    The next function in the module is getPoint. It takes two int seeds to generate two random lists, and returns a triple, where the x-coordinate is the point divided by 10 (to obtain a float between 0 and 1) and the other two coordinates are random floats also between 0 and 1, extracted as the v-th component of each of the random lists... I think that was not very clear, I will clarify soon:

    getPoint :: Int -> Int -> Data.Graph.Vertex -> GraphPoint
    getPoint r1 r2 v = let
                         randList1 = randomList r1
                         randList2 = randomList r2
                       in ((fromIntegral v)/10, randList1!!v, randList2!!v)


    The function myPoints is just a synonym of points3D, but without the graph parameter in the middle:

    myPoints :: Int -> Int -> [ GraphPoint ]
    myPoints r1 r2 = points3D graph1 r1 r2

    The function myEdges is a list of pairs of GraphPoints, where each component is the point previously generated for the corresponding vertex. The unit cube is just a reference that will be useful when we rotate the graph:

    myEdges :: Int -> Int -> [ (GraphPoint, GraphPoint) ]
    myEdges r1 r2 = [ (getPoint r1 r2 v1, getPoint r1 r2 v2) | (v1,v2) <- (edges graph1) ] ++ unitCube


    The origin is represented by a constant, not to repeat too many times the same triple:

    origin = (0.0, 0.0, 0.0)


    These three functions add one to each component of a point:

    plusX (x, y, z) = (x+1, y, z)
    plusY (x, y, z) = (x, y+1, z)
    plusZ (x, y, z) = (x, y, z+1)

    These functions are used to generate the reference axis:

    unitCube = [ (origin , plusZ origin ), -- +z
                 (plusZ origin , plusX $ plusZ origin), -- +x
                 (plusX $ plusZ origin , plusX origin ), -- -z
                 (plusX origin , origin ), -- -x
                 (origin , plusX origin ), -- +x
                 (plusX origin , plusY $ plusX origin), -- +y
                 (plusY $ plusX origin , plusY origin ), -- -x
                 (plusY origin , origin ), -- -y
                 (origin , plusY origin ), -- +y
                 (plusY origin , plusZ $ plusY origin), -- +z
                 (plusZ $ plusY origin , plusZ origin ), -- -y
                 (plusZ origin , origin ) -- -z
               ]

    Finally, the function randomList takes a seed and returns an infinite list of floats:

    randomList :: Int -> [Float]
    randomList seed = randoms (mkStdGen seed) :: [Float]


More code and graphics comming soon...

jueves, 11 de agosto de 2016

Desafío AlexCipher agosto 2016!!!

DESAFIO agosto de 2016:

  Te dejo un pequeño desafío. El siguiente texto fue encriptado utilizando un cuadrado Latino aleatorio de orden 256. El código del cifrador AlexCipher está en mi proyecto en GitHub.com. El texto original es en español. Si lográs descifrar el texto, dejame un comentario o enviame un mail por favor a bluemontag@gmail.com !!

Enlace al texto cifrado

 El ganador recibirá... bueno, unas buenas felicitaciones!! 

(Excepto que alguien quiera financiar esto :D)

miércoles, 13 de abril de 2016

Sobre Ignacio Gallego Sagastume


En abril de 2004 obtuvo la Licenciatura en Informática en la Facultad de Informática de la Universidad Nacional de La Plata (UNLP), eligiendo un perfil especializado en métodos formales en materias optativas (investigación operativa, teoría de grafos, verificación de programas, teoría de la computación, lenguajes de descripción de hardware, etc.). Realizó su tesis de licenciatura sobre una semántica de lenguajes web (como PHP, ASP, y JSP), bajo la dirección de la Dra. Claudia Pons.

Desde entonces ha trabajado primero en el Ministerio de Economía de la provincia de Buenos Aires y  luego en el centro de cómputo de ARBA (Agencia de Recaudación de la provincia de Buenos Aires), como analista, desarrollador y arquitecto web, en tecnologías Java, J2EE, JSP, Struts, Hibernate, ASP.NET, ASP 3.0, Oracle, SQL Server.

También ha trabajado en diversos proyectos freelance en PHP y MySQL.

En agosto de 2015 obtuvo el título de Magister en Ingeniería de Software, también en la UNLP. Realizó su trabajo de tesis de Magister en temas relacionados con la aplicación y generación de cuadrados Latinos en criptografía, también bajo la dirección de la Dra. Claudia Pons.

Es miembro de la ACM (Association for Computing Machinery) desde 2013.

Ir a Publicaciones

Artículos y Publicaciones

Cuadrados Latinos en criptografía