/*
 * Created by SharpDevelop.
 * User: Davide
 * Date: 08/09/2012
 * Time: 20:28
 * 
 * To change this template use Tools | Options | Coding | Edit Standard Headers.
 */
using System;
using System.Collections.Generic;
using SMath.Manager;
using SMath.Math;

namespace NonlinearSolvers.Primitives
{
	/// <summary>
	/// Bisected Direct Quadratic Regula Falsi root-finding method.
	/// </summary>
	public class BDQRF
	{
		static public string SMathPrimitiveName = "BDQRF";
		
		static public TermInfo GetTermInfos_3()
		{
			return
				new TermInfo
				(
					SMathPrimitiveName,
					TermType.Function,
					string.Format(Translator.GetString("RFM_BDQRF_dialog_3"), GlobalParams.DecimalPlaces),
					FunctionSections.Unknown,
					true,
					new ArgumentInfo(ArgumentSections.Function),
					new ArgumentInfo(ArgumentSections.Condition),
					new ArgumentInfo(ArgumentSections.Condition)
				);
		}
		
		static public TermInfo GetTermInfos_4()
		{
			return
				new TermInfo
				(
					SMathPrimitiveName,
					TermType.Function,
					Translator.GetString("RFM_BDQRF_dialog_4"),
					FunctionSections.Unknown,
					true,
					new ArgumentInfo(ArgumentSections.Function),
					new ArgumentInfo(ArgumentSections.Condition),
					new ArgumentInfo(ArgumentSections.Condition),
					new ArgumentInfo(ArgumentSections.Condition)
				);
		}
		
		static public TermInfo GetTermInfos_5()
		{
			return
				new TermInfo
				(
					SMathPrimitiveName,
					TermType.Function,
					Translator.GetString("RFM_BDQRF_dialog_5"),
					FunctionSections.Unknown,
					true,
					new ArgumentInfo(ArgumentSections.Function),
					new ArgumentInfo(ArgumentSections.Condition),
					new ArgumentInfo(ArgumentSections.Condition),
					new ArgumentInfo(ArgumentSections.Condition),
					new ArgumentInfo(ArgumentSections.Condition)
				);
		}
		
		static public TermInfo GetTermInfos_9()
		{
			return
				new TermInfo
				(
					SMathPrimitiveName,
					TermType.Function,
					Translator.GetString("RFM_BDQRF_dialog_9"),
					FunctionSections.Unknown,
					true,
					new ArgumentInfo(ArgumentSections.Function),
					new ArgumentInfo(ArgumentSections.Condition),
					new ArgumentInfo(ArgumentSections.Condition),
					new ArgumentInfo(ArgumentSections.Condition),
					new ArgumentInfo(ArgumentSections.Condition),
					new ArgumentInfo(ArgumentSections.RealNumber),
					new ArgumentInfo(ArgumentSections.Variable),
					new ArgumentInfo(ArgumentSections.Variable),
					new ArgumentInfo(ArgumentSections.Variable)
				);
		}
		
		static public bool IsSupport(Term root)
		{
			if (root.Text == SMathPrimitiveName && ((root.ArgsCount >= 3 && root.ArgsCount <= 5) || root.ArgsCount == 9)) return true;
			
			return false;
		}
		
		static public bool LowLevelEvaluation( Term root, Term[][] args, ref Store context, ref Term[] result )
		{
			Term[]
				arg1 = Computation.Preprocessing(args[0], ref context),
				arg2 = Computation.Preprocessing(args[1], ref context),
				arg3 = Computation.Preprocessing(args[2], ref context),
				// optional arguments
				arg4 = (root.ArgsCount >= 4) ? Computation.Preprocessing(args[3], ref context) : TermsConverter.ToTerms("0"),
				arg5 = (root.ArgsCount >= 5) ? Computation.Preprocessing(args[4], ref context) : TermsConverter.ToTerms("0");
			
			// preprocessing the output pattern, clear LHS defined variables.
			Term[] outputPattern = SharedFunctions.pattern2std(args[1], ref context);
			
			// initialize an empty store
			Store store = context.Scope.CreateStore();
			store.FileName = context.FileName;
			store.Optimization = context.Optimization;
			store.FractionsTypeNumeric = context.FractionsTypeNumeric;
			
			// preprocessing functions and starting values
			arg1 = SharedFunctions.equations2std(context.Scope, arg1, arg2);
			arg2 = SharedFunctions.variables2std(context.Scope, arg2);
			arg3 = SharedFunctions.variables2std(context.Scope, arg3);
			
			// get the variable units
			List<Term> myBuffer = new List<Term>();
			myBuffer.AddRange(arg2);
			myBuffer.Add(new Term(Functions.Abs, TermType.Function, 1));
			myBuffer.AddRange(arg3);
			myBuffer.Add(new Term(Functions.Abs, TermType.Function, 1));
			myBuffer.Add(new Term(Operators.Composition, TermType.Operator, 2));
			Term[] units = UoM.getUnits(context.Scope, myBuffer.ToArray());
			
			// make anything dimensionless
			arg1 = UoM.remUnits(context.Scope, arg1);
			arg2 = UoM.remUnits(context.Scope, arg2);
			arg3 = UoM.remUnits(context.Scope, arg3);
			arg4 = UoM.remUnits(context.Scope, arg4);
			arg5 = UoM.remUnits(context.Scope, arg5);
			
			// set the target precision of the function
			if (root.ArgsCount == 3)
			{
				myBuffer.Clear();
				myBuffer.Add(new Term("10", TermType.Operand, 0));
				myBuffer.Add(new Term(GlobalParams.DecimalPlaces.ToString(), TermType.Operand, 0));
				myBuffer.Add(new Term(Operators.Subtraction, TermType.Operator, 1));
				myBuffer.Add(new Term(Operators.Exponentiation, TermType.Operator, 2));
				arg4 = Computation.Preprocessing(myBuffer.ToArray(), ref store);
			}
			
			// number of iterations (default)
			int arg6 = 100;
			
			// show the number of iterations (default)
			bool arg7 = false;
			
			// show the step-by-step summary (default)
			bool arg8 = false;
			
			// save the step-by-step summary (default)
			bool arg9 = false;
			
			if (root.ArgsCount == 9)
			{
				int a6 = Int32.Parse(TermsConverter.ToString(Computation.SymbolicCalculation(args[5], ref context)));
				arg6 = (a6 > 0) ? a6 : arg6;
				arg7 = (TermsConverter.ToString(Computation.Preprocessing(args[6], ref context)) != "0") ? true : arg7;
				arg8 = (TermsConverter.ToString(Computation.Preprocessing(args[7], ref context)) != "0") ? true : arg8;
				arg9 = (TermsConverter.ToString(Computation.Preprocessing(args[8], ref context)) != "0") ? true : arg9;
			}
			
			// main alghoritm
			Term[][] solution = trySolve(arg1, arg2, arg3, arg4, arg5, arg6, arg8, arg9, ref store);
			
			List<Term> myResult = new List<Term>();
			
			// add the result
			myResult.AddRange(SharedFunctions.PostProcessing(context.Scope, solution[0], units, outputPattern, ref context));
			
			int i = 1;
			
			// add the number of iterations
			if (arg7)
			{
				myResult.AddRange(solution[1]);
				myResult.Add(new Term("1", TermType.Operand, 0));
				myResult.Add(new Term("1", TermType.Operand, 0));
				myResult.Add(new Term(Functions.Mat, TermType.Function, 3));
				i++;
			}
			// add the summary
			if (arg8)
			{
				myResult.AddRange(solution[2]);
				i++;
			}
			// output vector
			if (i > 1)
			{
				myResult.Add(new Term(i.ToString(), TermType.Operand, 0));
				myResult.Add(new Term("1", TermType.Operand, 0));
				myResult.Add(new Term(Functions.Mat, TermType.Function, 2 + i));
			}
			
			// output the result
			result = myResult.ToArray();
			
			return true;
		}
		
		static public Term[][] trySolve(Term[] function, Term[] bound_P1, Term[] bound_P2, Term[] epsilon_fx, Term[] epsilon_Δx, int maxIterations, bool showSummary, bool saveSummary, ref Store context)
		{
			List<Term>
				myBuffer = new List<Term>(),	// Buffer to use for chunks of code
				myResult = new List<Term>();	// Result Term List
			
			// initialize the step-by-step summary
			Summary summary = new Summary("BDQRF", context.FileName);
			
			if (showSummary || saveSummary)
			{
				summary.AddLine(new string[]{"i","Xdn","Xup","X","f(Xdn)","f(Xup)","f(X)"}, saveSummary);
			}
			
			int i = 0;
			
			Term[]
				num_0 = TermsConverter.ToTerms("0"),
				num_1 = TermsConverter.ToTerms("1"),
				num_2 = TermsConverter.ToTerms("2"),
				num_4 = TermsConverter.ToTerms("4"),
				var_X = bound_P1,
				var_Xdn = null,
				var_Xup = null,
				var_Xm = null,
				var_Y = SharedFunctions.setValues(context.Scope, function, var_X),
				var_Ydn = null,
				var_Yup = null,
				var_Ym = null,
				var_Xlast = TermsConverter.ToTerms("10^300"),
				var_D = null,
				var_a = null,
				var_b = null;
			
			// if/else #1, condition: [y<0]
			myBuffer.Clear();
			myBuffer.AddRange(var_Y);
			myBuffer.AddRange(num_0);
			myBuffer.Add(new Term(Operators.BooleanLess, TermType.Operator, 2));
			myBuffer.Add(new Term(Functions.Eval, TermType.Function, 1));
			string myFlag = Computation.Preprocessing(myBuffer.ToArray(), ref context)[0].Text;
			
			if (myFlag == "1")
			{
				var_Xdn = bound_P1;
				var_Xup = bound_P2;					
			}
			else
			{
				var_Xdn = bound_P2;
				var_Xup = bound_P1;	
			}
			
			// [y.dn:f(x.dn)]
			var_Ydn = SharedFunctions.setValues(context.Scope, function, var_Xdn);
			// [y.up:f(x.up)]
			var_Yup = SharedFunctions.setValues(context.Scope, function, var_Xup);

			// check that [sign(f(x.dn))*sign(f(x.up))<=0]
			myBuffer.Clear();
			myBuffer.AddRange(var_Ydn);
			myBuffer.Add(new Term(Functions.Sign, TermType.Function, 1));
			myBuffer.AddRange(var_Yup);
			myBuffer.Add(new Term(Functions.Sign, TermType.Function, 1));
			myBuffer.Add(new Term(Operators.Multiplication, TermType.Operator, 2));
			myBuffer.AddRange(num_0);
			myBuffer.Add(new Term(Operators.BooleanMore, TermType.Operator, 2));
			myBuffer.Add(new Term(Functions.Eval, TermType.Function, 1));
			myFlag = Computation.Preprocessing(myBuffer.ToArray(), ref context)[0].Text;
			
			// wrong delimiters
			if (myFlag == "1")
			{
				throw new Exception(Translator.GetString("error_delimiters"));
			}
			
//			// HACK: set an unit of measurement for X.last
//			// UoM(abs(X.dn)+abs(X.up)) to have a non zero value (zero values in SMath have no units)
//			myBuffer.Clear();
//			myBuffer.AddRange(var_Xdn);
//			myBuffer.Add(new Term(Functions.Abs, TermType.Function, 1));
//			myBuffer.AddRange(var_Xup);
//			myBuffer.Add(new Term(Functions.Abs, TermType.Function, 1));
//			myBuffer.Add(new Term(Operators.Composition, TermType.Operator, 2));
//			Term[] var_UoM = Computation.Preprocessing(myBuffer.ToArray(), ref context);
//			// d*UoM(abs(a)+abs(b))
//			myBuffer.Clear();
//			myBuffer.AddRange(var_Xlast);
//			myBuffer.AddRange(UoM.getUnits(var_UoM, context));
//			myBuffer.Add(new Term(Operators.Multiplication, TermType.Operator, 2));
//			var_Xlast = Computation.Preprocessing(myBuffer.ToArray(), ref context);
			
			bool solved = false;
			
			while (i <= maxIterations)
			{
				// [D:{x.up-x.dn}/2]
				myBuffer.Clear();
				myBuffer.AddRange(var_Xup);
				myBuffer.AddRange(var_Xdn);
				myBuffer.Add(new Term(Operators.Subtraction, TermType.Operator, 2));
				myBuffer.AddRange(num_2);
				myBuffer.Add(new Term(Operators.Division, TermType.Operator, 2));
				var_D = Computation.Preprocessing(myBuffer.ToArray(), ref context);
				
				// [x.m:{x.up+x.dn}/2]
				myBuffer.Clear();
				myBuffer.AddRange(var_Xup);
				myBuffer.AddRange(var_Xdn);
				myBuffer.Add(new Term(Operators.Composition, TermType.Operator, 2));
				myBuffer.AddRange(num_2);
				myBuffer.Add(new Term(Operators.Division, TermType.Operator, 2));
				var_Xm = Computation.Preprocessing(myBuffer.ToArray(), ref context);
				// [y.m:f(x.m)]
				var_Ym = SharedFunctions.setValues(context.Scope, function, var_Xm);
				
				// [a:{y.up+y.dn-2*y.m}/{2*D^2}]
				myBuffer.Clear();
				myBuffer.AddRange(var_Yup);
				myBuffer.AddRange(var_Ydn);
				myBuffer.Add(new Term(Operators.Composition, TermType.Operator, 2));
				myBuffer.AddRange(num_2);
				myBuffer.AddRange(var_Ym);
				myBuffer.Add(new Term(Operators.Multiplication, TermType.Operator, 2));
				myBuffer.Add(new Term(Operators.Subtraction, TermType.Operator, 2));
				myBuffer.AddRange(num_2);
				myBuffer.AddRange(var_D);
				myBuffer.AddRange(num_2);
				myBuffer.Add(new Term(Operators.Exponentiation, TermType.Operator, 2));
				myBuffer.Add(new Term(Operators.Multiplication, TermType.Operator, 2));
				myBuffer.Add(new Term(Operators.Division, TermType.Operator, 2));
				var_a = Computation.Preprocessing(myBuffer.ToArray(), ref context);
				
				// [b:{y.up-y.dn}/{2*D}]
				myBuffer.Clear();
				myBuffer.AddRange(var_Yup);
				myBuffer.AddRange(var_Ydn);
				myBuffer.Add(new Term(Operators.Subtraction, TermType.Operator, 2));
				myBuffer.AddRange(num_2);
				myBuffer.AddRange(var_D);
				myBuffer.Add(new Term(Operators.Multiplication, TermType.Operator, 2));
				myBuffer.Add(new Term(Operators.Division, TermType.Operator, 2));
				var_b = Computation.Preprocessing(myBuffer.ToArray(), ref context);
				
				// [x:x.m-{2*y.m}/{b*(1+sqrt(1-{4*a*y.m}/{b^2}))}]
				myBuffer.Clear();
				myBuffer.AddRange(var_Xm);				
				myBuffer.AddRange(num_2);
				myBuffer.AddRange(var_Ym);
				myBuffer.Add(new Term(Operators.Multiplication, TermType.Operator, 2));				
				myBuffer.AddRange(var_b);				
				myBuffer.AddRange(num_1);				
				myBuffer.AddRange(num_1);				
				myBuffer.AddRange(num_4);
				myBuffer.AddRange(var_a);
				myBuffer.AddRange(var_Ym);
				myBuffer.Add(new Term(Operators.Multiplication, TermType.Operator, 2));
				myBuffer.Add(new Term(Operators.Multiplication, TermType.Operator, 2));				
				myBuffer.AddRange(var_b);
				myBuffer.AddRange(num_2);
				myBuffer.Add(new Term(Operators.Exponentiation, TermType.Operator, 2));				
				myBuffer.Add(new Term(Operators.Division, TermType.Operator, 2));				
				myBuffer.Add(new Term(Operators.Subtraction, TermType.Operator, 2));				
				myBuffer.Add(new Term(Functions.Sqrt, TermType.Function, 1));				
				myBuffer.Add(new Term(Operators.Composition, TermType.Operator, 2));				
				myBuffer.Add(new Term(Operators.Multiplication, TermType.Operator, 2));				
				myBuffer.Add(new Term(Operators.Division, TermType.Operator, 2));				
				myBuffer.Add(new Term(Operators.Subtraction, TermType.Operator, 2));
				myBuffer.Add(new Term(Functions.Eval, TermType.Function, 1));	
				var_X = Computation.Preprocessing(myBuffer.ToArray(), ref context);
				
				// [y:f(x)]
				var_Y = SharedFunctions.setValues(context.Scope, function, var_X);
				
				// thresholds check: [[abs(y)<epsilon|abs(Δx)<epsilon]
				solved = SharedFunctions.isTargetBelowThreshold(context.Scope, var_Y, epsilon_fx);
				
				if (i > 0 && !solved)
				{
					myBuffer.Clear();
					myBuffer.AddRange(var_X);
					myBuffer.AddRange(var_Xlast);
					myBuffer.Add(new Term(Operators.Subtraction, TermType.Operator, 2));
					solved = SharedFunctions.isTargetBelowThreshold(context.Scope, myBuffer.ToArray(), epsilon_Δx);
				}
				
				// record this step
				if (showSummary || saveSummary)
				{
					summary.AddLine(context.Scope, new Term[][]{
					                	TermsConverter.ToTerms(i.ToString()),
					                	var_Xdn,
					                	var_Xup,
					                	var_X,
					                	var_Ydn,
					                	var_Yup,
					                	var_Y
					                }, saveSummary);
				}
				
				// required precision reached
				if (solved)
				{
					break;
				}
				
				// [x.last:x]
				var_Xlast = var_X;

				// if/else #4, condition: [y>0]
				myBuffer.Clear();
				myBuffer.AddRange(var_Y);
				myBuffer.AddRange(num_0);
				myBuffer.Add(new Term(Operators.BooleanMore, TermType.Operator, 2));
				myBuffer.Add(new Term(Functions.Eval, TermType.Function, 1));
				myFlag = Computation.Preprocessing(myBuffer.ToArray(), ref context)[0].Text;

				if (myFlag == "1")
				{
					var_Xup = var_X;
					var_Yup = var_Y;
					
					// if/else #5, condition: [y.m<0]
					myBuffer.Clear();
					myBuffer.AddRange(var_Ym);
					myBuffer.AddRange(num_0);
					myBuffer.Add(new Term(Operators.BooleanLess, TermType.Operator, 2));
					myBuffer.Add(new Term(Functions.Eval, TermType.Function, 1));
					myFlag = Computation.Preprocessing(myBuffer.ToArray(), ref context)[0].Text;
					
					if (myFlag == "1")
					{
						var_Xdn = var_Xm;
						var_Ydn = var_Ym;
					}
				}
				else
				{
					var_Xdn = var_X;
					var_Ydn = var_Y;
					
					// if/else #6, condition: [y.m>0]
					myBuffer.Clear();
					myBuffer.AddRange(var_Ym);
					myBuffer.AddRange(num_0);
					myBuffer.Add(new Term(Operators.BooleanMore, TermType.Operator, 2));
					myBuffer.Add(new Term(Functions.Eval, TermType.Function, 1));
					myFlag = Computation.Preprocessing(myBuffer.ToArray(), ref context)[0].Text;
					
					if (myFlag == "1")
					{
						var_Xup = var_Xm;
						var_Yup = var_Ym;
					}
				}
				
				i++;
			}
			
			if (solved)
			{
				// Output results and iterations
				return new Term[][]
				{
					var_X,
					TermsConverter.ToTerms(i.ToString()),
					summary.GetMatrix()
				};
			}
			else
			{
				// No solutions using the assigned iteration limit
				throw new Exception(string.Format(Translator.GetString("error_iterations"), --i));
			}
			
		}
	}
}