(!*********************************************************************
   Mosel NL examples
   =================
   file polygon8.mos
   `````````````````
   Maximize the area of polygon of N vertices and diameter of 1. 
   
   The position of vertices is indicated as (rho,theta) coordinates
   where rho denotes the distance to the base point (vertex with number N) 
   and theta the angle from the x-axis.
   
   -- Formulation using a simple Mosel user function redirected into Java --
 
   !!! Before running this model, compile Polygon.java into Polygon.class 
 
   (c) 2018 Fair Isaac Corporation
       author: Z.Csizmadia, Sep. 2018
*********************************************************************!)

model "Polygon 8 JavaSimpleFunction"
 uses "mmxnlp"
 uses 'mosjvm'

 parameters
  N=5                               ! Number of vertices
 end-parameters

 declarations
  Vertices = 1..N
  Polar = {"rho","theta"}
  Area: nlctr
  rho: array(Vertices) of mpvar     ! Distance of vertex from the base point
  theta: array(Vertices) of mpvar   ! Angle from x-axis
  D: array(Vertices,Vertices) of nlctr         ! Limit on side length
  FunctionArg: array(Vertices,Polar) of nlctr  ! User function arguments
  AreaFunction: userfunc                       ! User function definition
 end-declarations

! Objective - sum of areas. Definition of a user function
 AreaFunction := userfuncMosel("AreaInJava")

! Create function arguments
 forall(i in Vertices) do
  FunctionArg(i, "rho")   := rho(i)
  FunctionArg(i, "theta") := theta(i)
 end-do

! Use the Mosel user function in a formula for the objective
 Area := F(AreaFunction,FunctionArg)

! Bounds and start values for decision variables
 forall(i in 1..N-1) do
  rho(i) >= 0.1
  rho(i) <= 1
  setinitval(rho(i),4*i*(N + 1 - i)/((N+1)^2))
  setinitval(theta(i),M_PI*i/N)
 end-do

! Third side of all triangles <= 1
 forall (i in 1..N-2, j in i+1..N-1)
  D(i,j) := rho(i)^2 + rho(j)^2 - rho(i)*rho(j)*2*cos(theta(j)-theta(i)) <= 1

! Vertices in increasing order
 forall(i in 2..N-1) theta(i) >= theta(i-1) +.01

! Boundary conditions (last vertex above x-axis)
 theta(N-1) <= M_PI

! Abort model if we encounter a Java exception
 setparam('jvmabortonexception', true)

! Tell Java to look for classes in working directory
 setparam('jvmclasspath', getparam('workdir'))

! Optional parameter settings
 setparam("xnlp_verbose", true)         ! Enable XNLP output log

 setparam("xnlp_solver", 0)             ! Select SLP as the solver

! Uncomment to display user function info
 userfuncinfo(AreaFunction)

! Solve the problem
 maximise(Area)

! Solution output
 writeln("Area = ", getobjval)
 forall(i in 1..N-1)
  writeln("V", i, ": r=", getsol(rho(i)), " theta=", getsol(theta(i)))

! **** Definition of the Mosel user function ****
 public function AreaInJava(I: array(Indices: set of integer, Types: set of string) of real): real
  declarations
  ! Structures for communicating the input variable values to Java
   IndicesToJava=0..N-1
   rhoToJava: array(IndicesToJava) of real
   thetaToJava: array(IndicesToJava) of real
  end-declarations
  ! Project "rho" and "theta" values to structures formatted for Java
  forall(i in Vertices) rhoToJava(i-1) := I(i,"rho")
  forall(i in Vertices) thetaToJava(i-1) := I(i,"theta")

  returned := jvmcallreal( "Polygon.CalculateArea", rhoToJava, thetaToJava )
 end-function

end-model