//  This file is part of ff3d - http://www.freefem.org/ff3d
//  Copyright (C) 2001, 2002, 2003 Stphane Del Pino

//  This program is free software; you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation; either version 2, or (at your option)
//  any later version.

//  This program is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.

//  You should have received a copy of the GNU General Public License
//  along with this program; if not, write to the Free Software Foundation,
//  Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  

//  $Id: ScalarFunctionBuilder.cpp,v 1.11 2007/05/20 23:22:26 delpinux Exp $

#include <ScalarFunctionBuilder.hpp>

#include <ScalarFunctionConstant.hpp>

#include <ScalarFunctionCFunction.hpp>

#include <ScalarFunctionUnaryMinus.hpp>
#include <ScalarFunctionNot.hpp>

#include <ScalarFunctionModulo.hpp>

#include <ScalarFunctionSum.hpp>
#include <ScalarFunctionDifference.hpp>
#include <ScalarFunctionProduct.hpp>
#include <ScalarFunctionDivision.hpp>
#include <ScalarFunctionPower.hpp>

#include <ScalarFunctionMin.hpp>
#include <ScalarFunctionMax.hpp>

#include <ScalarFunctionGreaterThan.hpp>
#include <ScalarFunctionGreaterEqual.hpp>
#include <ScalarFunctionLowerThan.hpp>
#include <ScalarFunctionLowerEqual.hpp>

#include <ScalarFunctionNotEqual.hpp>
#include <ScalarFunctionEqual.hpp>

#include <ScalarFunctionAnd.hpp>
#include <ScalarFunctionOr.hpp>
#include <ScalarFunctionXor.hpp>

#include <FEMFunctionBase.hpp>

#include <FEMFunctionBuilder.hpp>

class ScalarFunctionBuilder::Simplifier
{
  struct BinaryOperatorModulo
  {
    /** 
     * Static evaluation of modulo of values
     * 
     * @param x first operand
     * @param y second operand
     * 
     * @return @f$ x%y @f$
     */
    static real_t eval(const real_t& x,const real_t& y) { return static_cast<int>(x)%static_cast<int>(y); }
  };

  struct BinaryOperatorSum
  {
    /** 
     * Static evaluation of sum of values
     * 
     * @param x first operand
     * @param y second operand
     * 
     * @return @f$ x+y @f$
     */
    static real_t eval(const real_t& x,const real_t& y) { return x+y; }
  };

  struct BinaryOperatorDifference
  {
    /** 
     * Static evaluation of difference of values
     * 
     * @param x first operand
     * @param y second operand
     * 
     * @return @f$ x-y @f$
     */
    static real_t eval(const real_t& x,const real_t& y) { return x-y; }
  };

  struct BinaryOperatorProduct
  {
    /** 
     * Static evaluation of product of values
     * 
     * @param x first operand
     * @param y second operand
     * 
     * @return @f$ x\cdot y @f$
     */
    static real_t eval(const real_t& x,const real_t& y) { return x*y; }
  };

  struct BinaryOperatorDivision
  {
    /** 
     * Static evaluation of division of values
     * 
     * @param x first operand
     * @param y second operand
     * 
     * @return @f$ x/y @f$
     */
    static real_t eval(const real_t& x,const real_t& y) { return x/y; }
  };

  struct BinaryOperatorPower
  {
    /** 
     * Static evaluation of power of values
     * 
     * @param x first operand
     * @param y second operand
     * 
     * @return @f$ x^y @f$
     */
    static real_t eval(const real_t& x,const real_t& y) { return std::pow(x,y); }
  };

  struct BinaryOperatorMin
  {
    /** 
     * Static evaluation of min of values
     * 
     * @param x first operand
     * @param y second operand
     * 
     * @return @f$ \min(x,y) @f$
     */
    static real_t eval(const real_t& x,const real_t& y) { return std::min(x,y); }
  };

  struct BinaryOperatorMax
  {
    /** 
     * Static evaluation of max of values
     * 
     * @param x first operand
     * @param y second operand
     * 
     * @return @f$ max(x,y) @f$
     */
    static real_t eval(const real_t& x,const real_t& y) { return std::max(x,y); }
  };

  /** 
   * Computes constant function resulting of two constant functions
   * 
   * @param f1 first constant function
   * @param f2 second constant function
   * 
   * @return @f$ f1~\mbox{op}~f2 @f$
   */
  template <typename BinaryOperatorType>
  ConstReferenceCounting<ScalarFunctionBase>
  __getOperatorF1F2SimplifiedFunction(const ScalarFunctionConstant& f1,
				      const ScalarFunctionConstant& f2) const
  {
    return new ScalarFunctionConstant(BinaryOperatorType::eval(f1(0),
							       f2(0)));
  }

  /** 
   * Computes finite element function resulting of a constant function
   * and a finite element function binary operation
   * 
   * @param f1 constant function
   * @param f2 finite element function
   * 
   * @return @f$ f1~\mbox{op}~f2 @f$
   */
  template <typename BinaryOperatorType>
  ConstReferenceCounting<ScalarFunctionBase>
  __getOperatorF1F2SimplifiedFunction(const ScalarFunctionConstant& f1,
				      const FEMFunctionBase& f2) const
  {
    Vector<real_t> femValues(f2.values().size());

    const real_t constantValue = f1(0);

    const real_t outsideValue = BinaryOperatorType::eval(constantValue,
							 f2.outsideValue());

    for (size_t i=0; i<femValues.size(); ++i) {
      femValues[i] = BinaryOperatorType::eval(constantValue,f2[i]);
    }

    FEMFunctionBuilder builder;
    builder.build(f2.discretizationType(),
		  f2.baseMesh(),
		  femValues,
		  outsideValue);
    return builder.getBuiltScalarFunction();
  }

  /** 
   * Computes finite element function resulting of a finite element
   * function and a constant function binary operation
   * 
   * @param f1 finite element function
   * @param f2 constant function
   * 
   * @return @f$ f1~\mbox{op}~f2 @f$
   */
  template <typename BinaryOperatorType>
  ConstReferenceCounting<ScalarFunctionBase>
  __getOperatorF1F2SimplifiedFunction(const FEMFunctionBase& f1,
				      const ScalarFunctionConstant& f2) const
  {
    Vector<real_t> femValues(f1.values().size());

    const real_t constantValue = f2(0);

    for (size_t i=0; i<femValues.size(); ++i) {
      femValues[i] = BinaryOperatorType::eval(f1[i],constantValue);
    }

    const real_t outsideValue = BinaryOperatorType::eval(f1.outsideValue(),
							 constantValue);

    FEMFunctionBuilder builder;
    builder.build(f1.discretizationType(),
		  f1.baseMesh(),
		  femValues,
		  outsideValue);
    return builder.getBuiltScalarFunction();
  }

  /** 
   * Computes finite element function resulting of two finite element
   * function binary operation
   * 
   * @param f1 first finite element function
   * @param f2 second finite element function
   * 
   * @return @f$ f1~\mbox{op}~f2 @f$
   */
  template <typename BinaryOperatorType>
  ConstReferenceCounting<ScalarFunctionBase>
  __getOperatorF1F2SimplifiedFunction(const FEMFunctionBase& f1,
				      const FEMFunctionBase& f2) const
  {
    Vector<real_t> femValues(f1.values().size());

    ASSERT(f1.discretizationType() == f2.discretizationType());
    ASSERT(f1.baseMesh() == f2.baseMesh());

    for (size_t i=0; i<femValues.size(); ++i) {
      femValues[i] = BinaryOperatorType::eval(f1[i],f2[i]);
    }

    const real_t outsideValue = BinaryOperatorType::eval(f1.outsideValue(),
							 f2.outsideValue());

    FEMFunctionBuilder builder;
    builder.build(f1.discretizationType(),
		  f1.baseMesh(),
		  femValues,
		  outsideValue);
    return builder.getBuiltScalarFunction();
  }


  /** 
   * Specialize the first function type of the binary operation
   * 
   * @param f specialized function
   * @param f2 still base function
   * 
   * @return simplified (or not) binary operation
   */
  template <typename BinaryOperatorType,
	    typename FunctionType1>
  ConstReferenceCounting<ScalarFunctionBase>
  __getOperatorF1SimplifiedFunction(const FunctionType1& f,
				    const ScalarFunctionBase& f2) const
  {
    switch (f2.type()) {
    case ScalarFunctionBase::constant: {
      const ScalarFunctionConstant& g = static_cast<const ScalarFunctionConstant&>(f2);
      return this->__getOperatorF1F2SimplifiedFunction<BinaryOperatorType>(f,g);
    }
    case ScalarFunctionBase::femfunction: {
      const FEMFunctionBase& g = static_cast<const FEMFunctionBase&>(f2);
      return this->__getOperatorF1F2SimplifiedFunction<BinaryOperatorType>(f,g);
    }
    default: {
      throw ErrorHandler(__FILE__,__LINE__,
			 "unexpected function type",
			 ErrorHandler::unexpected);
    }
    }
    return 0;
  }

  /** 
   * Specialize the binary operation type
   * 
   * @param f1 first operand
   * @param f2 second operand
   * 
   * @return simplified (or not) binary operation
   */  
  template <typename BinaryOperatorType>
  ConstReferenceCounting<ScalarFunctionBase>
  __getOperatorSimplifiedFunction(const ScalarFunctionBase& f1,
				  const ScalarFunctionBase& f2) const
  {
    switch (f1.type()) {
    case ScalarFunctionBase::constant: {
      const ScalarFunctionConstant& f = static_cast<const ScalarFunctionConstant&>(f1);
      return this->__getOperatorF1SimplifiedFunction<BinaryOperatorType>(f,f2);
    }
    case ScalarFunctionBase::femfunction: {
      const FEMFunctionBase& f = static_cast<const FEMFunctionBase&>(f1);
      if (f2.type() == ScalarFunctionBase::femfunction) {
	const FEMFunctionBase& g = static_cast<const FEMFunctionBase&>(f2);
	if (not(f.hasSameType(g))) {
	  return 0;
	}
      }
      return this->__getOperatorF1SimplifiedFunction<BinaryOperatorType>(f,f2);
    }
    default: {
      throw ErrorHandler(__FILE__,__LINE__,
			 "unexpected function type",
			 ErrorHandler::unexpected);
    }
    }
    return 0;
  }

public:
  /** 
   * Simplifies a binary operation
   * 
   * @param binaryOperation given binary operation
   * @param f first function
   * @param g second function
   * 
   * @return simplified (if possible) binary operation
   */
  ConstReferenceCounting<ScalarFunctionBase>
  simplify(const BinaryOperation& binaryOperation,
	   const ScalarFunctionBase& f,
	   const ScalarFunctionBase& g) const
  {
    if (not(f.canBeSimplified() and g.canBeSimplified())) {
      return 0;
    }

    switch (binaryOperation.type()) {
    case BinaryOperation::modulo: {
      return __getOperatorSimplifiedFunction<BinaryOperatorModulo>(f,g);
    }
    case BinaryOperation::sum: {
      return __getOperatorSimplifiedFunction<BinaryOperatorSum>(f,g);
    }
    case BinaryOperation::difference: {
      return __getOperatorSimplifiedFunction<BinaryOperatorDifference>(f,g);
    }
    case BinaryOperation::product: {
      return __getOperatorSimplifiedFunction<BinaryOperatorProduct>(f,g);
    }
    case BinaryOperation::division: {
      return __getOperatorSimplifiedFunction<BinaryOperatorDivision>(f,g);
    }
    case BinaryOperation::power: {
      return __getOperatorSimplifiedFunction<BinaryOperatorPower>(f,g);
    }
    case BinaryOperation::min: {
      return __getOperatorSimplifiedFunction<BinaryOperatorMin>(f,g);
    }
    case BinaryOperation::max: {
      return __getOperatorSimplifiedFunction<BinaryOperatorMax>(f,g);
    }
      // these functions are not treated
    case BinaryOperation::gt:
    case BinaryOperation::lt:
    case BinaryOperation::ge:
    case BinaryOperation::le:
    case BinaryOperation::eq:
    case BinaryOperation::ne:
    case BinaryOperation::or_:
    case BinaryOperation::xor_:
    case BinaryOperation::and_: {
      return 0;
    }
    case BinaryOperation::undefined:
    default: {
      throw ErrorHandler(__FILE__,__LINE__,
			 "not implemented",
			 ErrorHandler::unexpected);
      return 0;
    }
    }
    return 0;
  }


  /** 
   * Constructor
   * 
   */
  Simplifier()
  {
    ;
  }

  /** 
   * Destructor
   * 
   */
  ~Simplifier()
  {
    ;
  }
};

void ScalarFunctionBuilder::
setFunction(ConstReferenceCounting<ScalarFunctionBase> f)
{
  ASSERT(__builtFunction == 0);
  __builtFunction = f;
}

void ScalarFunctionBuilder::
setBinaryOperation(const BinaryOperation binaryOperation,
		   ConstReferenceCounting<ScalarFunctionBase> f)
{
  ASSERT(__builtFunction != 0);

  Simplifier simplifier;

  ConstReferenceCounting<ScalarFunctionBase>
    simplifiedFunction = simplifier.simplify(binaryOperation,
					     *__builtFunction,
					     *f);
  if (simplifiedFunction) {
    ffout(4) << " * algebraic simplification of "
	     << *__builtFunction << ' ' << binaryOperation << ' ' << *f << '\n';
    __builtFunction = simplifiedFunction;
  } else {
    switch (binaryOperation.type()) {
    case BinaryOperation::modulo: {
      __builtFunction = new ScalarFunctionModulo(__builtFunction, f);
      break;
    }
    case BinaryOperation::sum: {
      __builtFunction = new ScalarFunctionSum(__builtFunction, f);
      break;
    }
    case BinaryOperation::difference: {
      __builtFunction = new ScalarFunctionDifference(__builtFunction, f);
      break;
    }
    case BinaryOperation::product: {
      __builtFunction = new ScalarFunctionProduct(__builtFunction, f);
      break;
    }
    case BinaryOperation::division: {
      __builtFunction = new ScalarFunctionDivision(__builtFunction, f);
      break;
    }
    case BinaryOperation::power: {
      __builtFunction = new ScalarFunctionPower(__builtFunction, f);
      break;
    }

    case BinaryOperation::min: {
      __builtFunction = new ScalarFunctionMin(__builtFunction, f);
      break;
    }
    case BinaryOperation::max: {
      __builtFunction = new ScalarFunctionMax(__builtFunction, f);
      break;
    }

      // Boolean operations
    case BinaryOperation::gt: {
      __builtFunction = new ScalarFunctionGreaterThan(__builtFunction, f);
      break;
    }
    case BinaryOperation::ge: {
      __builtFunction = new ScalarFunctionGreaterEqual(__builtFunction, f);
      break;
    }
    case BinaryOperation::lt: {
      __builtFunction = new ScalarFunctionLowerThan(__builtFunction, f);
      break;
    }
    case BinaryOperation::le: {
      __builtFunction = new ScalarFunctionLowerEqual(__builtFunction, f);
      break;
    }
    case BinaryOperation::eq: {
      __builtFunction = new ScalarFunctionEqual(__builtFunction, f);
      break;
    }
    case BinaryOperation::ne: {
      __builtFunction = new ScalarFunctionNotEqual(__builtFunction, f);
      break;
    }

    case BinaryOperation::and_: {
      __builtFunction = new ScalarFunctionAnd(__builtFunction, f);
      break;
    }
    case BinaryOperation::or_: {
      __builtFunction = new ScalarFunctionOr(__builtFunction, f);
      break;
    }
    case BinaryOperation::xor_: {
      __builtFunction = new ScalarFunctionXor(__builtFunction, f);
      break;
    }

    default: {
      throw ErrorHandler(__FILE__,__LINE__,
			 "Not implemented binary operation",
			 ErrorHandler::unexpected);
    }
    }
  }
}

void ScalarFunctionBuilder::
setUnaryMinus()
{
  ASSERT(__builtFunction != 0);
  switch (__builtFunction->type()) {
  case ScalarFunctionBase::constant: {
    const ScalarFunctionConstant& f
      = static_cast<const ScalarFunctionConstant&>(*__builtFunction);
    const real_t value = -f(0);

    ffout(4) << " * algebraic simplification of -" << f << '\n';

    __builtFunction = new ScalarFunctionConstant(value);
    break;
  }
  case ScalarFunctionBase::femfunction: {
    const FEMFunctionBase& f
      = static_cast<const FEMFunctionBase&>(*__builtFunction);
    Vector<real_t> femValues(f.values().size());

    for (size_t i=0; i<femValues.size(); ++i) {
      femValues[i] = -f[i];
    }

    const real_t outsideValue = -f.outsideValue();

    ffout(4) << " * algebraic simplification of -" << f << '\n';

    FEMFunctionBuilder builder;
    builder.build(f.discretizationType(),
		  f.baseMesh(),
		  femValues,
		  outsideValue);
    __builtFunction = builder.getBuiltScalarFunction();
    break;
  }
  default: {
    __builtFunction = new ScalarFunctionUnaryMinus(__builtFunction);
  }
  }
}

void ScalarFunctionBuilder::
setNot()
{
  ASSERT(__builtFunction != 0);
  __builtFunction = new ScalarFunctionNot(__builtFunction);
}


void ScalarFunctionBuilder::
setCFunction(const std::string& cfunction)
{
  ASSERT(__builtFunction != 0);
  __builtFunction = new ScalarFunctionCFunction(cfunction, __builtFunction);
}