#include <oxstd.h>
#include <maximize.h>
#import <solvenle>

BroydenMod(const avF, const vX)
{
	avF[0] = (3 - 2 * vX) .* vX - lag0(vX, 1) - 2 * lag0(vX, -1) + 1;
	// rank deficient modification
	decl n = sizer(vX);
	avF[0][n - 1] = sqr(avF[0][n - 1]);
	return 1;
}
Broyden(const avF, const vX)
{
	avF[0] = (3 - 2 * vX) .* vX - lag0(vX, 1) - 2 * lag0(vX, -1) + 1;
	return 1;
}
BroydenJac(const amJ, const vX)
{
	decl n = sizer(vX);
	amJ[0] = diag(3 - 4 * vX) - lag0(unit(n), 1) - 2 * lag0(unit(n), -1);
	return 1;
}
static decl s_dLambda = -1;
Bratu(const avF, const vX)
{
	decl n = int(round(sqrt(sizer(vX)))), h = 1 / (n + 1);
	decl u = shape(vX, n, n), m = - 4 * u;
	m += (0 | u[ : n - 2][]);
	m += (u[1 : ][] | 0);
	m += (0 ~ u[][ : n - 2]);
	m += (u[][1 : ] ~ 0);
	m /= sqr(h);
	m += s_dLambda * exp(u);
	avF[0] = vec(m);
	return 1;
}

ModifiedRosenbrock(const avF, const vX)
{
	decl n2 = idiv(sizer(vX), 2);
	decl even = 2 * range(0, n2 - 1), odd = even + 1;
	decl f = zeros(2 * n2, 1);
	f[even] = 10 * (vX[even] - sqr(vX[odd]));
	f[odd] = 1 / (1 + exp(-vX[odd])) - 0.73;
	avF[0] = f;
	return 1;
}
AugmentedRosenbrock(const avF, const vX)
{
	decl n4 = idiv(sizer(vX), 4);
	decl i0 = 4 * range(0, n4 - 1), i1 = i0 + 1, i2 = i0 + 2, i3 = i0 + 3;
	decl f = zeros(4 * n4, 1);
	f[i0] = vX[i0];
	f[i1] = 10 * (vX[i2] - sqr(vX[i1]));
	f[i2] = 1 - vX[i0];
	f[i3] = 1.25 * vX[i3] - 0.25 * vX[i3] .^ 3;
	avF[0] = f;
	return 1;
}

main()
{
	decl x, n, time;
	MaxControl(-1, 10, -1);

	x = -ones(100, 1) * 100;
	println("=== Broyden system, n=", sizer(x));
	time = timer();
	SolveNLE(Broyden, &x, 0, BroydenJac);
	println("Broyden system, time=", timespan(time));

	x = -ones(100, 1) * 100;
	println("\n=== Modified Broyden system, n=", sizer(x));
	time = timer();
	SolveNLE(BroydenMod, &x);
	println("time=", timespan(time));
	
	x = vec(zeros(32,32));
//	s_dLambda = -1;
	s_dLambda = 6.80673;
	println("\n=== Bratu system, n=", sizer(x));
	time = timer();
	SolveNLE(Bratu, &x);
	println("time=", timespan(time));

//	x = vec(zeros(320,320));
////	s_dLambda = -1;
//	s_dLambda = 6.80673;
//	println("\n=== large Bratu system, n=", sizer(x));
//	time = timer();
//	SolveNLE(Bratu, &x);
//	println("time=", timespan(time));

//	n = 25;
//	x = vecr(constant(-1.8,n,1) ~ constant(-1,n,1));
//	SolveNLE(ModifiedRosenbrock, &x);
//
//	n = 13;
//	x = vecr(constant(-1.2,n,1) ~ constant(1,n,1) ~ constant(-1,n,1) ~ constant(20,n,1));
//	SolveNLE(AugmentedRosenbrock, &x);
}