"""
Based on:

cyipopt: Python wrapper for the Ipopt optimization package, written in Cython.

Copyright (C) 2012-2015 Amit Aides
Copyright (C) 2015-2017 Matthias Kümmerer
Copyright (C) 2017-2024 cyipopt developers

License: EPL 2.0
"""
# Based on matlab code by Peter Carbonetto.
# ruff: noqa: N801, D101, D102, D103, N806, ARG002, D400, F841, PLR2004, INP001

import logging

import cyipopt

import numpy as np



class hs071:

    def objective(self, x):
        if semaphore[2] == 1 and semaphore[0] >= semaphore[1]:
            value = 1 / 0
        return x[0] * x[3] * np.sum(x[0:3]) + x[2]

    def gradient(self, x):
        if semaphore[2] == 2 and semaphore[0] >= semaphore[1]:
            value = 1 / 0
        return np.array(
            [
                x[0] * x[3] + x[3] * np.sum(x[0:3]),
                x[0] * x[3],
                x[0] * x[3] + 1.0,
                x[0] * np.sum(x[0:3]),
            ]
        )

    def constraints(self, x):
        if semaphore[2] == 3 and semaphore[0] >= semaphore[1]:
            value = 1 / 0
        return np.array((np.prod(x), np.dot(x, x)))

    def jacobian(self, x):
        if semaphore[2] == 4 and semaphore[0] >= semaphore[1]:
            value = 1 / 0
        return np.concatenate((np.prod(x) / x, 2 * x))

    def hessianstructure(self):
        if semaphore[2] == 5:
            value = 1 / 0
        return np.nonzero(np.tril(np.ones((4, 4))))

    def hessian(self, x, lagrange, obj_factor):
        if (semaphore[2] == 6 and semaphore[0] >= semaphore[1]) or semaphore[2] == 7:
            value = 1 / 0
        H = obj_factor * np.array(
            (
                (2 * x[3], 0, 0, 0),
                (x[3], 0, 0, 0),
                (x[3], 0, 0, 0),
                (2 * x[0] + x[1] + x[2], x[0], x[0], 0),
            )
        )
        H += lagrange[0] * np.array(
            (
                (0, 0, 0, 0),
                (x[2] * x[3], 0, 0, 0),
                (x[1] * x[3], x[0] * x[3], 0, 0),
                (x[1] * x[2], x[0] * x[2], x[0] * x[1], 0),
            )
        )
        H += lagrange[1] * 2 * np.eye(4)
        row, col = self.hessianstructure()
        return H[row, col]

    def intermediate(*args):
        semaphore[0] = args[2]
        if semaphore[2] == 8 and semaphore[0] >= semaphore[1]:
            value = 1 / 0
        print(f"Objective value at iteration #{args[2]} is {args[3]}")


def main():
    x0 = [1.0, 5.0, 5.0, 1.0]
    lb = [1.0, 1.0, 1.0, 1.0]
    ub = [5.0, 5.0, 5.0, 5.0]
    cl = [25.0, 40.0]
    cu = [2.0e19, 40.0]
    nlp = cyipopt.Problem(n=len(x0), m=len(cl), problem_obj=hs071(), lb=lb, ub=ub, cl=cl, cu=cu)
    nlp.add_option("max_iter", 15)
    cyipopt.set_logging_level(logging.DEBUG)

    semaphore[:] = [0, 3, 0]
    failure_case = input(
        "Enter number of failure case to consider:\n"
        "0: No failure case\n"
        "1: Fail in objective after iteration 3\n"
        "2: Fail in gradient after iteration 3\n"
        "3: Fail in constraint after iteration 3\n"
        "4: Fail in jacobian after iteration 3\n"
        "5: Fail in hessianstructure (status -1, max iterations exceeded)\n"
        "6: Fail in hessian after iteration 3 (status -1, max iterations exceeded)\n"
        "7: Fail in hessian immediately (status -1, max iterations exceeded)\n"
        "8: Fail in intermediate after iteration 3\n"
    )
    semaphore[2] = int(failure_case)

    # Solve the problem
    x, info = nlp.solve(x0)
    print(f"Solution of the primal variables: {x=}\n")
    print(f"Solution of the dual variables: lambda={info['mult_g']}\n")
    print(f"Objective={info['obj_val']}\n")
    print(f"Status code={info['status']}\n")


semaphore = [0, 0, 0]  # current iteration, iteration to fail at, failure case number


if __name__ == "__main__":
    main()