#include <oxstd.h>
#import <maxsqp>

func(const vP, const adFunc, const avScore, const amHessian)
{
	decl f = vP[0] * sqr(vP[1]);
	adFunc[0] = -f;
	if (isarray(avScore))
	{
		decl score = sqr(vP[1]) | 2 * vP[0] * vP[1];
		avScore[0] = -score;
	}
	return 1;
}
cfunc_gt0(const avF, const vP)
{
	avF[0] = matrix(2 - sqr(vP[0]) - sqr(vP[1]));
	return 1;
}
cfunc_gt0_jac(const amJacobian, const vP)
{
	amJacobian[0] = -2 * vP[0] ~ -2 * vP[1];
	return 1;
}

main()
{
	MaxControl(-1, 50, TRUE);
	MaxControlEps(1e-6, 1e-3);
	decl vp, f, vlo, vhi;

	vp = <0;0>;
	vlo = <-2;-2>;
	vhi = <-0.1;-0.1>;

	println("Problem:");
	println("    min  x[0]*x[1]^2");
	println("    s.t. x[0]^2 + x[1]^2 <= 2");
	println("In format for MaxSQP/MaxSQPF:");
	println("    max -x[0]*x[1]^2");
	println("    s.t. 2 - x[0]^2 - x[1]^2 >= 0");
	println("\n------ MaxSQPF using analytical derivatives");
	MaxSQPF(func, &vp, &f, 0, FALSE, cfunc_gt0, 0, vlo, vhi, cfunc_gt0_jac);
	println("Maximum found at x=", vp);

	println("\n------ MaxSQPF using numerical derivatives for score and Jacobian)");
	vp = <0;0>;
	MaxSQPF(func, &vp, &f, 0, TRUE, cfunc_gt0, 0, vlo, vhi);
	println("Maximum found at x=", vp);

	println("\n------ MaxSQP using analytical derivatives");
	vp = <0;0>;
	MaxSQP(func, &vp, &f, 0, FALSE, cfunc_gt0, 0, vlo, vhi, cfunc_gt0_jac);
	println("Maximum found at x=", vp);
}