kantorovich with scala - v00

kantorovich_v00.scala

import breeze.optimize.linear._

object Main { 
  def main(args: Array[String]) {

	  println("\nWelcome.\nOpen the console and type e.g. new kantorovich(Array(1.0/7, 2.0/7, 4.0/7), Array(1.0/4, 1.0/4, 1.0/2)).run()\n")

  }
}


class kantorovich(Mu: Array[Double], Nu: Array[Double]) {
  
  def run() {

	val n = Mu.length

	// ** objective function ** //

	// defines the distance matrix
	val D =  Array.fill(n,n)(0);
	for( i <- 0 to n-1 ) 
		for( j <- 0 to n-1 )
			if( !(i==j) ) D(i)(j) = 1;
	// defines the new Linear Program object
	val lp = new breeze.optimize.linear.LinearProgram();
	// matrix of parameters p_{ij} 
	val P =  Array.fill(n,n)(lp.Real());
	// objective function 
	val Objective = ( 
	for( i <- 0 to 2 ) 
			yield ( 
				for( (x, a) <- P(i) zip D(i) ) 
					yield (x * a) 
				).reduce(_ + _)
	).reduce(_ + _) 
	 
	 
	// ** constraints ** //

	// there are 2n equality constraints and n² positivity constraints
	var Constraints = new Array[lp.Constraint](2*n+n*n)

	// equality constraints 
	val tP=P.transpose
	for (i <- 0 to n-1) {
	Constraints(i) = P(i).reduce[lp.Expression](_ + _)  =:=  Mu(i)  
	Constraints(i+n) = tP(i).reduce[lp.Expression](_ + _)  =:=  Nu(i)
	}
	// positivity constraints 
	for (i <- 0 to n-1) {
		var k = 2*n + n*i
		for (j <- 0 to n-1) Constraints(k+j) = P(i)(j) >=0 
	}
	val lpp = ( 
	-Objective
	subjectTo( Constraints:_* )
	)
	
	val result = lp.maximize(lpp)
	
	println("\n --- \n")
	println(result)
	println("\n --- \n Kantorovich distance: " + -result.value)
  	
  }
}