#!/usr/bin/env python3
"``naginterfaces.library.opt.handle_solve_ipopt`` Python Example."

# nAG Copyright 2018-2020.

# pylint: disable=invalid-name,too-many-arguments,too-many-locals,too-many-statements

from math import sqrt

import numpy as np

from naginterfaces.library import opt

def main():
    """
    Example for :func:`naginterfaces.library.opt.handle_solve_ipopt`.

    Interior-point solver for sparse NLP.

    >>> main()
    naginterfaces.library.opt.handle_solve_ipopt Python Example Results.
    Solving a problem based on Hock and Schittkowski Problem 73.
    Solving with a nonlinear objective.
    At the solution the objective function is 2.9894378e+01.
    """

    print(
        'naginterfaces.library.opt.handle_solve_ipopt Python Example Results.'
    )
    print('Solving a problem based on Hock and Schittkowski Problem 73.')
    print('Solving with a nonlinear objective.')

    # The 'non'linear objective:
    cb_objfun = lambda x, inform: (
        24.55*x[0] + 26.75*x[1] + 39.*x[2] + 40.5*x[3], inform,
    )

    def cb_objgrd(x, fdx, inform):
        fdx[:] = [24.55, 26.75, 39., 40.5]
        return inform

    cb_confun = lambda x, _ncnln, inform: (
        [
            12.0*x[0] + 11.9*x[1] + 41.8*x[2] + 52.1*x[3] - 1.645 *
            sqrt(
                0.28*x[0]**2 + 0.19*x[1]**2 + 20.5*x[2]**2 + 0.62*x[3]**2
            )
        ],
        inform,
    )

    def cb_congrd(x, gdx, inform):
        """The Jacobian of the nonlinear constraints."""
        tmp = sqrt(
            0.62*x[3]**2 + 20.5*x[2]**2 + 0.19*x[1]**2 + 0.28*x[0]**2
        )
        gdx[:] = [
            (12.0*tmp-0.4606*x[0])/tmp,
            (11.9*tmp-0.31255*x[1])/tmp,
            (41.8*tmp-33.7225*x[2])/tmp,
            (52.1*tmp-1.0199*x[3])/tmp,
        ]
        return inform

    def cb_hess(x, idf, sigma, lamda, hx, inform):
        """Hessian function."""
        nnzh = hx.size
        if idf == 0:
            # Objective is linear.
            for i in range(nnzh):
                hx[i] = 0.0
            return inform
        tmp = sqrt(
            0.62*x[3]**2 + 20.5*x[2]**2 + 0.19*x[1]**2 + 0.28*x[0]**2
        )
        tmp = tmp*(
            x[3]**2 + 33.064516129032258064*x[2]**2 +
            0.30645161290322580645*x[1]**2 + 0.45161290322580645161*x[0]**2
        )
        hx[0] = (
            -0.4606*x[3]**2 - 15.229516129032258064*x[2]**2 -
            0.14115161290322580645*x[1]*x[1])/tmp
        hx[1] = (0.14115161290322580645*x[0]*x[1])/tmp
        hx[2] = (15.229516129032258064*x[0]*x[2])/tmp
        hx[3] = (0.4606*x[0]*x[3])/tmp
        hx[4] = (
            -0.31255*x[3]**2 - 10.334314516129032258*x[2]**2 -
            0.14115161290322580645*x[0]**2
        )/tmp
        hx[5] = (10.334314516129032258*x[1]*x[2])/tmp
        hx[6] = (0.31255*x[1]*x[3])/tmp
        hx[7] = (
            -33.7225*x[3]**2 - 10.334314516129032258*x[1]**2 -
            15.229516129032258065*x[0]**2
        )/tmp
        hx[8] = (33.7225*x[2]*x[3])/tmp
        hx[9] = (
            -33.7225*x[2]**2 - 0.31255*x[1]**2 - 0.4606*x[0]**2
        )/tmp
        if idf == -1:
            for i in range(nnzh):
                hx[i] = hx[i] * lamda[0]
        else:
            assert idf == 1
        return inform

    # The initial guess:
    x = [1., 1., 1., 1.]
    nvar = len(x)
    nnzu = 2*nvar
    # Create a handle for the problem:
    handle = opt.handle_init(nvar)

    # Define the bounds:
    opt.handle_set_simplebounds(
        handle,
        bl=[0.]*nvar, bu=[1.e20]*nvar,
    )

    # Define the non-linear objective:
    opt.handle_set_nlnobj(
        handle,
        idxfd=[1, 2, 3, 4],
    )

    # Define the linear constraints:
    opt.handle_set_linconstr(
        handle,
        bl=[5., 1.], bu=[1.e20, 1.],
        irowb=[1, 1, 1, 1, 2, 2, 2, 2],
        icolb=[1, 2, 3, 4, 1, 2, 3, 4],
        b=[2.3, 5.6, 11.1, 1.3, 1.0, 1.0, 1.0, 1.0],
    )
    nnzu += 2*2

    # Define the nonlinear constraints:
    opt.handle_set_nlnconstr(
        handle,
        bl=[21.], bu=[1.e20],
        irowgd=[1, 1, 1, 1], icolgd=[1, 2, 3, 4],
    )
    nnzu += 2*1

    # Define the structure of the Hessian, dense upper triangle:
    irowh = np.empty(10, dtype=int)
    icolh = np.empty(10, dtype=int)
    idx = 0
    for i in range(1, nvar + 1):
        for j in range(i, nvar + 1):
            icolh[idx] = j
            irowh[idx] = i
            idx += 1
    # For the objective function:
    opt.handle_set_nlnhess(
        handle,
        idf=0, irowh=irowh, icolh=icolh,
    )
    # For the constraint function:
    opt.handle_set_nlnhess(
        handle,
        idf=1, irowh=irowh, icolh=icolh,
    )

    opt.handle_opt_set(
        handle, 'Stop Tolerance 1 = 2.5e-8',
    )
    opt.handle_opt_set(
        handle, 'Print Level = 0',
    )

    # Solve the problem:
    ip_soln = opt.handle_solve_ipopt(
        handle, x,
        objfun=cb_objfun, objgrd=cb_objgrd, confun=cb_confun, congrd=cb_congrd,
        hess=cb_hess,
    )

    print(
        'At the solution the objective function is {:.7e}.'.format(
            ip_soln.rinfo[0]
        )
    )

    # Destroy the handle:
    opt.handle_free(handle)

if __name__ == '__main__':
    import doctest
    import sys
    sys.exit(
        doctest.testmod(
            None, verbose=True, report=False,
            optionflags=doctest.ELLIPSIS | doctest.REPORT_NDIFF,
        ).failed
    )