#include "gradient_descent.h"

////////////////////////////////////////////////////////////////////////////////
//////////////////////////////Functii matematice////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
real f2d1(real x, real y)	// functia patratica de ord 2; pt [-2;2]
{
	return x*x+y*y;
}

real f2d2(real x, real y)	// functia Rosenbrock; pt [-2;2]
{
	return sin(0.5*x*x-0.25*y*y+3)*cos(2*x+1-exp(y));
}
real f2d3(real x, real y)   // functia banana
{
	return (1-x)*(1-x) + 100*(y-x*x)*(y-x*x);
}
real f2d4(real x, real y)
{
	return (0.5*x*x+0.2*y*y);
}
real f2d5(real x, real y)
{
	return (0.5*(x-1)*(x-1)+1.5*(y-0.5)*(y-0.5));
}
real f2d6(real x, real y)
{
	return (0.5*(x-1)*(x-1)+0.5*x*y+1.5*(y-0.5)*(y-0.5));
}
////////////////////////////////////////////////////////////////////////////////
//////////////////////////////Optimizare////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
bool check_params(real x1, real x2, real y1, real y2)
{
	if (x1>=x2)	return false;
	if (y1>=y2) return false;
	return true;
}
bool mesh2d(real(*pf)(real,real), real x1, real x2, real y1, real y2, real pasx, real pasy)
{
	FILE *fp, *fp2;
	fp  = fopen("mesh.in","w");
	fp2 = fopen("param.in","w");
	dword 	par_i = 0,
			par_j = 0;
	bool checked = false;
	if (!fp)
	{
		cout<<"mesh: EROARE! Nu am putut deschide fisierul mesh.in!"<<endl;
		fclose(fp);
		return false;
	}
	if (!fp2)
	{
		cout<<"mesh: EROARE! Nu am putut deschide fisierul param.in!"<<endl;
		fclose(fp2);
		fclose(fp);
		return false;
	}
	if (!check_params(x1,x2,y1,y2))
	{
		fclose(fp2);
		fclose(fp);
		return false;
	}
	fprintf(fp2, "%llf %llf %llf %llf %llf %llf\n", x1, x2, y1, y2, pasx, pasy);
	for (real i=x1; i<=x2; i+=pasx)
	{
		par_i++;
		for (real j=y1; j<=y2; j+=pasy)
			{
				if (!checked) par_j++;
				fprintf(fp, "%llf %llf %llf\n",i,j,pf(i,j));
			}
		checked = true;
	}
	fprintf(fp2, "%li %li\n", par_i, par_j);
	fclose(fp2);
	fclose(fp);
	return true;
}
inline real norm2d(real x, real y)
{
	return sqrt(x*x+y*y); // este sqrt(VT*V)
}
inline real dist2d(real x1, real y1, real x2, real y2)
{
	return sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
}
preal diff2d2(real (*pf)(real,real), real x, real y)
{
	preal temp = new real[2];
	temp[0] = (pf(x+h,y)-pf(x-h,y))/h;
	temp[1] = (pf(x,y+h)-pf(x,y-h))/h;
	return temp;
}
preal diff2d5(real (*pf)(real,real), real x, real y) // ret Jacobianul
{
	preal temp = new real[2];

	temp[0] = (-pf(x+2*h,y)+8*pf(x+h,y)-8*pf(x-h,y)+pf(x-2*h,y)) / (12*h);
	temp[1] = (-pf(x,y+2*h)+8*pf(x,y+h)-8*pf(x,y-h)+pf(x,y-2*h)) / (12*h);

	return temp;
}
inline preal genfibvec(dword n)
{
	preal fib = new real[n+1];
	dword i=0;

	fib[0] = 0;
	fib[1] = 1;
	fib[2] = 1;
	for (i=2; i<n+1; i++)
		fib[i] = fib[i-2] + fib[i-1];

	return fib;
}
real optmf1d(real (*pf)(real,real), preal fib, preal i, preal pc, preal sd)
{
	dword 	n = 47,
			k = 0;
	real 	ak = i[0],
			bk = i[1],
			ck = 0.0,
			dk = 0.0;
	bool dbg_cnd1 = fabs(pf(pc[0]+bk*sd[0],pc[1]+bk*sd[1])-pf(pc[0]+ak*sd[0],pc[1]+ak*sd[1]))<dltf,
		 dbg_cnd2 = fabs(bk-ak)<dltx;
	while ((fabs(bk-ak)>dltx) &&
				(fabs(pf(pc[0]+bk*sd[0],pc[1]+bk*sd[1])-pf(pc[0]+ak*sd[0],pc[1]+ak*sd[1]))>dltf) &&
				(k<n-2))
	{
		ck = ak + (1-fib[n-1-k]/fib[n-k]) *(bk-ak);
		dk = ak +  fib[n-1-k]/fib[n-k]    *(bk-ak);
        if (pf(pc[0]+ck*sd[0],pc[1]+ck*sd[1])<=pf(pc[0]+dk*sd[0],pc[1]+dk*sd[1]))
		{
			bk = dk;
		}
		else
		{
			ak = ck;
		}
		k++;
	}
	return fabs(ak+bk)/2;
}
bool optgd2d(real (*pf)(real,real), real x0, real y0, real a, real b)
{
	FILE *fp;
	dword n = 47;
	preal 	grad	 	= NULL,	// gradient intr-un punct
			sdir 		= NULL,	// directie de cautare
			interval 	= NULL, // interval de cautare pt 1D
			pcurent		= NULL, // punct curent
			pnou 		= NULL, // noul punct
			fib			= genfibvec(n);
	real 	alfa = 0,
			alfa_vechi = 0;
	dword 	contor = 0;
	bool 	cndstop = true;

	fp = fopen("grd_dsc.out","w");
	if (!pf)
	{
		cout<<"optgd2d: LIPSA! Lipseste functia ce trebuie optimizata!\n";
		return false;
	}
	if (!fp)
	{
		cout<<"optgd2d: EROARE! Nu am putut deschide fisierul grd_dsc.out!\n";
		fclose(fp);
		return false;
	}

	fprintf(fp,"%llf %llf %llf\n",x0,y0,pf(x0,y0));
	grad = diff2d5(pf, x0, y0);
	if (dist2d(grad[0],grad[1],0,0)>delta)
	{
		sdir = new real[2];
		sdir[0] = -grad[0]/norm2d(-grad[0],-grad[1]);
		sdir[1] = -grad[1]/norm2d(-grad[0],-grad[1]);
		interval = new real[2]; interval[0] = a; interval[1] = b;
		pcurent = new real[2];
		pcurent[0] = x0;
		pcurent[1] = y0;
		alfa = optmf1d(pf, fib, interval, pcurent, sdir);
		pnou = new real[2];
		pnou[0] = pcurent[0] + alfa * sdir[0];
		pnou[1] = pcurent[1] + alfa * sdir[1];
		fprintf(fp,"%llf %llf %llf\n",pnou[0],pnou[1],pf(pnou[0],pnou[1]));
		while ((dist2d(pcurent[0],pcurent[1],pnou[0],pnou[1])>delta) &&
				(contor<N) && cndstop && (alfa-alfa_vechi!=0))
		{
			delete []grad;
			delete []sdir;
			delete []pcurent;
			alfa_vechi = alfa;
			grad = diff2d5(pf,pnou[0],pnou[1]);
			//if (!(grad[0]==0 && grad[1]==0))
			if (dist2d(grad[0],grad[1],0,0)>delta)
			{
				sdir = new real[2];
				sdir[0] = -grad[0]/norm2d(-grad[0],-grad[1]);
				sdir[1] = -grad[1]/norm2d(-grad[0],-grad[1]);
				pcurent = pnou;
				alfa = optmf1d(pf, fib, interval, pcurent, sdir);
				pnou = new real[2];
				pnou[0] = pcurent[0] + alfa * sdir[0];
				pnou[1] = pcurent[1] + alfa * sdir[1];

				fprintf(fp,"%llf %llf %llf\n",pnou[0],pnou[1],pf(pnou[0],pnou[1]));
				contor++;
			} else cndstop = false;
		} //while
		delete []pnou;
		delete []pcurent;
		delete []interval;
		delete []sdir;
		delete []grad;
	} //if

	delete []fib;
	delete []grad;
	fclose(fp);
	return true;
}