#ifdef _DEBUG
#define _CRTDBG_MAP_ALLOC
#include <stdlib.h>
#include <crtdbg.h>
#endif
 
#include "HLBFGS.h"
 
#include "Lite_Sparse_Matrix.h"
 
#include <iostream>

Lite_Sparse_Matrix* m_sparse_matrix = 0;
 
//////////////////////////////////////////////////////////////////////////
void evalfunc(int N, double* x, double *prev_x, double* f, double* g)
{
	*f = 0;
	for (int i = 0; i < N; i+=2)
	{
		double T1 = 1 - x[i];
		double T2 = 10*(x[i+1]-x[i]*x[i]);
		*f += T1*T1+T2*T2;
		g[i+1]   = 20*T2;
		g[i] = -2*(x[i]*g[i+1]+T1);
	}
}
//////////////////////////////////////////////////////////////////////////
void newiteration(int iter, int call_iter, double *x, double* f, double *g,  double* gnorm)
{
	std::cout << iter <<": " << call_iter <<" " << *f <<" " << *gnorm  << std::endl;
}
//////////////////////////////////////////////////////////////////////////
void evalfunc_h(int N, double *x, double *prev_x, double *f, double *g, HESSIAN_MATRIX& hessian)
{
	//the following code is not optimal if the pattern of hessian matrix is fixed.
	if (m_sparse_matrix)
	{
		delete m_sparse_matrix;
	}
 
	m_sparse_matrix = new Lite_Sparse_Matrix(N, N, SYM_LOWER, CCS, FORTRAN_TYPE, true);
 
	m_sparse_matrix->begin_fill_entry();
 
	static bool first = true;
	double *diag = m_sparse_matrix->get_diag();
 
	if (first)
	{
		// you need to update f and g
		*f = 0;
		double tmp;
		for (int i = 0; i < N; i+=2)
		{
			tmp = x[i]*x[i];
			double T1 = 1 - x[i];
			double T2 = 10*(x[i+1]-tmp);
			*f += T1*T1+T2*T2;
			g[i+1]   = 20*T2;
			g[i] = -2*(x[i]*g[i+1]+T1);
			diag[i] = 2+1200*tmp-400*x[i+1];
			diag[i+1] = 200;
			m_sparse_matrix->fill_entry(i, i+1, -400*x[i]);
		}
	}
	else
	{
		for (int i = 0; i < N; i+=2)
		{
			diag[i] = 2+1200*x[i]*x[i]-400*x[i+1];
			diag[i+1] = 200;
			m_sparse_matrix->fill_entry(i, i+1, -400*x[i]);
		}
	}
 
	m_sparse_matrix->end_fill_entry();
 
	hessian.set_diag(m_sparse_matrix->get_diag());
	hessian.set_values(m_sparse_matrix->get_values());
	hessian.set_rowind(m_sparse_matrix->get_rowind());
	hessian.set_colptr(m_sparse_matrix->get_colptr());
	hessian.set_nonzeros(m_sparse_matrix->get_nonzero());
	first = false;
}
//////////////////////////////////////////////////////////////////////////
void Optimize_by_HLBFGS(int N, double *init_x, int num_iter, int M, int T, bool with_hessian)
{
	double parameter[20];
	int info[20];
	//initialize
	INIT_HLBFGS(parameter, info);
	info[4] = num_iter;
	info[6] = T;
	info[7] = with_hessian?1:0;
	info[10] = 0;
	info[11] = 1;
 
	if (with_hessian)
	{
		HLBFGS(N, M, init_x, evalfunc, evalfunc_h, HLBFGS_UPDATE_Hessian, newiteration, parameter, info);
	}
	else
	{
		HLBFGS(N, M, init_x, evalfunc, 0, HLBFGS_UPDATE_Hessian, newiteration, parameter, info);
	}
 
}
//////////////////////////////////////////////////////////////////////////
void main()
{
#ifdef _DEBUG
	_CrtSetDbgFlag(_CRTDBG_ALLOC_MEM_DF | _CRTDBG_LEAK_CHECK_DF);
#endif
	std::cout.precision(16);
	std::cout << std::scientific;
 
	int N = 1000;
	std::vector<double> x(N);
 
	for (int i = 0; i < N/2; i++)
	{
		x[2*i]   = -1.2;
		x[2*i+1] =  1.0;
	}
 
	int M = 7;
	int T = 0;
 
	//use Hessian
	// if M = 0, T = 0, it is Newton
	//Optimize_by_HLBFGS(N, &x[0], 1000, M, T, true);
 
	//without Hessian
	Optimize_by_HLBFGS(N, &x[0], 1000, M, T, false);  // it is LBFGS(M) actually, T is not used

	if (m_sparse_matrix)
		delete m_sparse_matrix;
}