// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
// Copyright (C) 2011 - DIGITEO - Michael Baudin
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution.  The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt

mprintf("Finds the largest n-polygon with diameter smaller than 1.\n");

function [ f , c , index ] = largesmallpolygon ( x , index )
    // Bibliography
    //   "Benchmarking Optimization Software with Cops"
    //   Dolan, Moré, 2001
    //   Section 1, "Largest Small Polygon"
    //
    //   "Biggest Little Polygon"
    //   http://mathworld.wolfram.com/BiggestLittlePolygon.html
    //
    //   Audet, C. "Optimisation globale structurée: propriétés, équivalences et résolution."
    //   Thèse de Doctorat. Montréal, Canada: École Polytechnique de Montréal, 1997.
    //   http://www.gerad.ca/Charles.Audet.
    //
    // Known optimal values are
    // A6 = 0.677981 (Wolfram Mathworld)
    // A8 = 0.726869 (Wolfram Mathworld)
    // A25 = 0.77974 (Dolan & Moré - SNOPT)
    // A50 = 0.784016 (Dolan & Moré - SNOPT)
    // A75 = 0.784769 (Dolan & Moré - SNOPT)
    // A100 = 0.785040 (Dolan & Moré - SNOPT)
    //
    nv = size(x,"*")/2
    f = []
    c = []
    // nv: number of vertices
    // x : a (2*nv)-by-1 matrix of doubles, where
    // r is in [0,1]
    // t is in [0,pi]
    r = x(1:nv)
    t = x(nv+1:2*nv)

    //polygon_draw(r ,t , nv);
    //pause

    if ( index == 2 | index == 6 ) then
        f = polygon_area(r ,t , nv)
        f = -f
    end
    if ( index == 5 | index == 6 ) then
        c = zeros(nv^2+nv-1,1)
        // Set the diameters
        d = polygon_diameters(r ,t , nv)
        d = matrix(d,nv^2,1)
        c(1:nv^2) = 1-d
        // Set the angles
        c(nv^2+1:nv^2+nv-1) = t(2:nv)-t(1:nv-1)
        // Transpose for neldermead
        c = c'
    end
endfunction

function polygon_draw(r ,t , nv)
    // Draws a polygon with nv vertices.
    // r : a nv-by-1 matrix of doubles, the radius
    // t : a nv-by-1 matrix of doubles, the angles
    x = r.*cos(t);
    y = r.*sin(t);
    i = 1;
    x($+1) = r(i)*cos(t(i));
    y($+1) = r(i)*sin(t(i));
    plot(x,y,"bo-")
    h = gcf()
    h.children.isoview = "on"
endfunction

function polygon_update(h,r ,t , nv)
    // Draws a polygon with nv vertices.
    // r : a nv-by-1 matrix of doubles, the radius
    // t : a nv-by-1 matrix of doubles, the angles
    x = r.*cos(t);
    y = r.*sin(t);
    i = 1;
    x($+1) = r(i)*cos(t(i));
    y($+1) = r(i)*sin(t(i));
    h.data = [x y]
endfunction

function [r,t] = polygon_regular (nv)
    // Returns a regular polygon with nv vertices
    // and unit radius, centered at origin.
    // r : a nv-by-1 matrix of doubles, the radius
    // t : a nv-by-1 matrix of doubles, the angles
    r = ones(nv,1);
    tnv = 2*%pi/nv
    t = linspace(-%pi+tnv,%pi,nv)'
endfunction

function f = polygon_area(r ,t , nv)
    // Returns the area of a polygon with nv vertices.
    // r : a nv-by-1 matrix of doubles, the radius
    // t : a nv-by-1 matrix of doubles, the angles
    f = r(nv) * r(1) * sin(t(1)-t(nv))
    for i = 1: nv-1
        f = f + r(i+1) * r(i) * sin(t(i+1)-t(i))
    end
    f = 0.5 * f
endfunction



function d = polygon_diameters(r ,t , nv)
    // Returns the diameters of a polygon with nv vertices.
    // r : a nv-by-1 matrix of doubles, the radius
    // t : a nv-by-1 matrix of doubles, the angles
    d = zeros(nv,nv)
    for i = 1 : nv
        for j = 1 : nv
            d(i,j) = r(i)^2 + r(j)^2 - 2*r(i)*r(j)*cos(t(i)-t(j))
            d(i,j) = abs(d(i,j))
            d(i,j) = sqrt(d(i,j))
        end
    end
endfunction




/////////////////////////////////////////////////
//
// Maximize the size of the polynomial
//
//
function stop = myoutputcmd(state, data)
    stop = %f
    iter = data.iteration
    if ( state == "init" ) then
        mprintf ( "=================================\n");
        mprintf ( "Initialization\n");
    elseif ( state == "done" ) then
        mprintf ( "=================================\n");
        mprintf ( "End of Optimization\n");
    end
    fc = data.funccount
    fval = data.fval
    x = data.x
    simplex = data.simplex
    step = data.step
    ssize = optimsimplex_size ( simplex )
    //
    // Plot current solution
    if ( %t & modulo(iter,10) == 0) then
        h = findobj ( "user_data" , "nmpolygon" );
        nv = size(x,"*")/2
        r = x(1:nv)
        t = x(nv+1:$)
        polygon_update(h.children.children,r ,t , nv)
        a = polygon_area(r ,t , nv);
        str = msprintf("Largest Small Polygon - Area=%f",a);
        h.title.text = str
    end
    //
    if ( modulo(iter,10) == 0 ) then
        mprintf ( "Iter. #%3d, Feval #%3d, Fval = %f, S = %.1e\n", ..
        iter, fc, fval, ssize);
    end
endfunction





function [A,r,t] = findlargestpolygon (nv)
    // Finds the largest smallest polygon with nv vertices
    // A : a 1-by-1 matrix of doubles, the area
    // r : a nv-by-1 matrix of doubles, the radius
    // t : a nv-by-1 matrix of doubles, the angle
    radius = 0.45;
    [r,t] = polygon_regular (nv);
    r = radius*r;
    scf()
    polygon_draw(r ,t , nv);
    h = gcf();
    h.children.data_bounds = [
    -0.6 -0.6
    0.6 0.6
    ];
    h.children.user_data = "nmpolygon";
    x0 = [r;t];
    index = 6;
    [ f0 , c0 , index ] = largesmallpolygon ( x0 , index );
    mprintf("Current area = %f\n",-f0);
    mprintf("Constraint satisfaction = %f (expected positive)\n",min(c0));
    //
    // Setup bounds
    rmin = zeros(nv,1);
    rmax = ones(nv,1);
    tmin = -ones(nv,1)*%pi;
    tmax = ones(nv,1)*%pi;
    xmin=[rmin;tmin];
    xmax=[rmax;tmax];
    //
    nm = neldermead_new ();
    nm = neldermead_configure(nm,"-numberofvariables",2*nv);
    nm = neldermead_configure(nm,"-function",largesmallpolygon);
    nm = neldermead_configure(nm,"-x0",x0);
    nm = neldermead_configure(nm,"-maxiter",2000);
    nm = neldermead_configure(nm,"-maxfunevals",2000);
    nm = neldermead_configure(nm,"-method","box");
    nm = neldermead_configure(nm,"-boundsmin",xmin');
    nm = neldermead_configure(nm,"-boundsmax",xmax');
    nm = neldermead_configure(nm,"-simplex0method","randbounds");
    nm = neldermead_configure(nm,"-nbineqconst",nv^2+nv-1);
    nm = neldermead_configure(nm,"-outputcommand",myoutputcmd);
    //
    // Check that the cost function is correctly connected.
    [ nm , result ] = neldermead_function ( nm , x0 );
    //
    // Perform optimization
    nm = neldermead_search(nm, "off");
    fopt = neldermead_get(nm,"-fopt")
    A = -fopt
    xopt = neldermead_get(nm,"-xopt")
    r = xopt(1:nv)
    t = xopt(nv+1:$)
    nm = neldermead_destroy(nm)
endfunction

////////////////////////////////////////////////////////
// Checking area for nv=6
// Use a regular hexagon (see Graham, page 5, Fig. 5).
// A = 0.64952...
nv = 6;
radius = 0.5;
[r,t] = polygon_regular (nv);
r = radius*r;
h = scf();
polygon_draw(r ,t , nv);
f = polygon_area(r ,t , nv);
//
// Check this
// A = 3*sqrt(3)*s^2/2, where s is the side length, i.e. the radius
A = 3*sqrt(3)*radius^2/2;

mprintf("Area =%f (expected = %f)\n",f,A);

d = polygon_diameters(r ,t , nv);
dmax = max(d);
mprintf("Maximum diameter=%f (expected = %f)\n",dmax,2*radius);
dmin = min(d(d<>0));
mprintf("Minimum diameter=%f (expected=%f)\n",dmin,radius);
close(h);

////////////////////////////////////////////////////////
//
// Solve problem
//
nv = 6;
rand("seed" , 0);
[A,r,t] = findlargestpolygon (nv);
mprintf("Maximum Area =%f (expected = %f)\n",A,0.677981);