Latex Equation Compiler eqc

/*******************************************************
    unit.cpp  -  Class for extending GiNaC to handle physical units
                             -------------------
    begin                : Mon Oct 22 2001
    copyright            : (C) 2001 by Jan Rheinlaender
    email                : jrheinlaender@users.sourceforge.net
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   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 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/

#include "unit.h"
#include "msgdriver.h" // *** added in 1.3.1
#include "optstack.h" // *** added in 1.1
#include "printing.h" // *** added in 1.4.1

// *** added print_func in 0.8
GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(Unit, basic, print_func<print_context>(&Unit::do_print).print_func<imathprint>(&Unit::do_print));

// Define static members
  bool Unit::initialized = false;
/*  Unit* Unit::metre; *** removed in 1.4.1
  Unit* Unit::kilogram;
  Unit* Unit::second;
  Unit* Unit::ampere;
  Unit* Unit::kelvin;
  Unit* Unit::celsius;
  Unit* Unit::mole;
  Unit* Unit::candela;*/
  std::map <std::string, Unit*> Unit::factory; // *** added in 1.4.1
  std::map <std::string, expression> Unit::conversions; // *** changed ex to expression in 0.4
  // *** moved initialization of the sets to Unit::init in 1.3.1
//  std::set<const char*, ltstr> Unit::base_units; *** removed in 1.4.1
  //unitvec *Unit::preferred_unit_conversions = new unitvec; *** removed in 1.1
  std::set<const char*, ltstr> Unit::unitnames;
  unit_init unit_init::unit_initializer; //  *** added in 0.5
  bool unit_init::called = false; //  *** added in 0.5

unit_init::unit_init() { // *** added in 0.5
  if (called) {
    // Removed in 1.3.1, msg_error might not be initialized yet 
    //    msg::error(0) << "Attempt to call Unit::init() more than one time!" << endline;
    return;
  } else
    called = true;
  Unit::init();
}

unit_init::~unit_init() { // *** added in 0.5
  Unit::clear();
}

// Required constructors and destructors and other GiNaC-specific methods
  Unit::Unit() {}; // *** changed in 1.2.1, changed in 1.3.1
// TODO: Make this private, outside callers should use Unit::get()
  Unit::Unit(const std::string &n) { // *** changed in 1.3.1
    name = n;
  }
  
  Unit Unit::get(const std::string &n) { // *** added in 1.4.1
    if (factory[n] == NULL) 
      (factory[n] = new Unit(n))->setflag(status_flags::dynallocated); // *** added dynallocated in 1.4.1
    return *factory[n];
  } // Unit::create()
  
  Unit* Unit::getp(const std::string &n) { // *** added in 1.4.1
    if (factory[n] == NULL)
      (factory[n] = new Unit(n))->setflag(status_flags::dynallocated); // *** added dynallocated in 1.4.1
    return factory[n];
  } // Unit::create()
                  
 
  void Unit::archive(archive_node &n) const
    {
      inherited::archive(n);
      n.add_string("Unit", name);
      n.add_string("Unit-printname", printname); // *** added in 1.4.1
    }

/* *** removed in 1.3.1
  Unit::Unit(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst)
    {
      n.find_string("Unit", name);
    }
*/

  // *** changed to read_archive format in 1.3.1
  void Unit::read_archive(const archive_node &n, lst &sym_lst)
    {
      inherited::read_archive(n, sym_lst);
      n.find_string("Unit", name);
      n.find_string("Unit-printname", printname);
//    return (new Unit(n, sym_lst))->setflag(status_flags::dynallocated);//
    }
GINAC_BIND_UNARCHIVER(Unit);  

  int Unit::compare_same_type(const basic &other) const {
    const Unit &o = static_cast<const Unit &>(other);
    if (name == o.name)
      return 0;
    else if (name < o.name)
      return -1;
    else
      return 1;
  }
  
  // Makes GiNaC more efficient
  bool Unit::is_equal_same_type(const basic & other) const { // *** added in 1.4.6
    const Unit & o = static_cast<const Unit &>(other);
          
    return (name == o.name);
  } // Unit::is_equal_same_type()

  void Unit::do_print(const print_context &c, unsigned level) const {
  // *** changed name to do_print in 0.8
    std::string thename = (printname == "" ? name : printname); // added in 1.4.1
  
    if (is_a<print_tree>(c)) {
      c.s << thename;
    } else if (is_a<print_csrc>(c)) {
      c.s << thename;
    } else if (is_a<print_latex>(c)) {
      c.s << ((printname == "") ? (std::string("\\") + thename) : thename);
    } else if (is_a<imathprint>(c)) { // added in 1.4.1 for iMath
      if (thename[0] == '\\')  thename.erase(0,1); // Handle broken names, should not occur at all
      
      if (thename[0] == '{') { // smath unit definition
        c.s << thename;
      } else if (thename[0] == '%') {
        if (thename.length() == 1)
          c.s << "\"%\""; // Percent as a unit
        else
          c.s << thename; // Beginning of a unit name
      } else { // protect the name with quotes
        c.s << "\"" + thename + "\"";
      }
    } else {
      c.s << "[" << thename << "]";
    }
  } // Unit::do_print()

  ex Unit::eval(int level) const
    // It is possible to define simplification or canonicalization rules here, like turning
    // mm to m*10^(-3) etc.
    {
      return this->hold();
    }

  void Unit::set_printname(const std::string& pn) {
    printname = pn;
    // *** added in 1.4.3
    if (pn != name) {
      add_conversion(pn, get(name));
      add_unitname(pn);
    }
  } // set_printname()

  void Unit::add_conversion(const std::string &uname, const expression& other_units) {
    if (conversions[uname] != numeric(0))
      MSG_WARN(0, "Warning: Overwriting unit definition of " << uname << endline);

    conversions[uname] = Unit::canonicalize(other_units);
    MSG_INFO(2, "Adding conversion [" << uname << "] = " << Unit::canonicalize(other_units) << endline);
  }

  void Unit::add_unitname(const std::string &unitname) {
    MSG_INFO(2, "Adding unitname: " << unitname << endline);
    unitnames = set_add(unitnames, unitname);
  } // Unit::add_unitname()

  bool Unit::is_unitname(const std::string &uname) {
    MSG_INFO(4, "Checking whether " << uname << " is a unit." << endline);
    MSG_INFO(5, unitnames << endline);
    return (set_has(unitnames, uname));
  } // Unit::is_unitname()
  
  bool Unit::is_unitnamep(const std::string &uname) {
    MSG_INFO(4, "Checking whether %" << uname << " is a unit." << endline);
    MSG_INFO(5, unitnames << endline);
    return (set_has(unitnames, std::string("%") + uname));
  } // Unit::is_unitnamep()

  unitvec *Unit::create_conversions(const exprvector &units, const bool always) {
  // *** added in 0.3, changed to use vector in 0.9, changed argument 2 to vector in 1.0
    unitvec *result = new unitvec;
    unitvec convs; // *** added in 0.9, changed to vector in 1.0

    for (exprvector::const_iterator i = units.begin(); i < units.end(); i++) {
      MSG_INFO(2, "Creating conversion for " << *i << endline);
      // *** changed to const_iterator in 0.7, to std::vector in 1.0
      expression pref_unit = *i;
      expression conversion = conversions[ex_to<Unit>(pref_unit).get_name()];
      if (conversion.is_zero()) { // *** Added in 1.4.1 because subst_units will crash with this conversion
        MSG_WARN(0, "Warning: No conversion exists for " << *i << ". Ignoring" << endline);
        continue;
      }
      // *** changed in 0.6 to avoid unnecessary substitutions, added algebraic in 0.8
      if (!result->empty()) conversion = conversion.csubs(convs, subs_options::algebraic);

      // Any numerics on the left hand side of the relation need to be moved to the
      // right hand side because 3\mm will not match "1000\mm == \m"
      if (is_a<mul>(conversion)) {
        // Other cases are either impossible or need not be handled
        mul conv = ex_to<mul>(conversion);
        conversion = (conversion / expression(conv.op(conv.nops())));
        pref_unit = pref_unit / expression(conv.op(conv.nops()));
        //TODO: if the numeric is a float, the above code will probably leave a coefficient of
        //0.99999999999995 or 1.00000000000005 on the left side!
      }

      if (always || (conversion != pref_unit)) { // *** added in 0.5
        // If someone says \preferred_units{\metre} no conversion is necessary!
        // But maybe he wants to override a global \preferred_units{\m} ?
        MSG_INFO(2, "Adding " << conversion << "==" << pref_unit << " to list of preferred units." << endline);
        result->push_back(conversion == pref_unit);
        convs.push_back(conversion == pref_unit); // *** changed to push_back() in 1.0
      }
    }

    return result;
  } // create_conversions()

  expression Unit::subst_units(const expression &e, const unitvec& unitlist) { // *** changed in 1.1, added unitlist in 1.4.1
    MSG_INFO(1, "Substituting units in " << e << endline);
    expression result = e;

    // Note that we can NOT substitute all units at once because there might be inter-dependencies
    // between the preferred units (e.g., \mm and N/mm^2)
    // Ginac documentation says all substitutions are done in parallel, but it doesn't work
    for (unitvec::const_iterator i = unitlist.begin(); i < unitlist.end(); i++) {
      MSG_INFO(2, "Substituting unit " << *i << endline);
      
      if (is_a<constant>(i->lhs()) && e.has(ex_to<constant>(i->lhs()))) { 
        // *** handles Pi = 180\degree, added in 1.4.5
        result = result.subs(*i, subs_options::algebraic);
      } else if (is_a<numeric>(expression(i->lhs()).evalf())) {  // *** added in 0.9, handles evalf()'ed cases for Pi = 180\degree
        if (is_a<equation>(e)) { // *** added in 1.4.1, handle special cases. Note equation is a subclass of relational
          // TODO: Add more intelligence to not do the substitution except on expressions containing Pi or a numeric
          // TODO: What about relationals? But relational does not have get_op()
          equation eq = ex_to<equation>(e);
          expression l = eq.lhs();
          expression r = eq.rhs();
          if (!is_a<symbol>(l))  {
            l = l / expression(expression(i->lhs()).evalf() / expression(i->rhs()));
          }
          if (!is_a<symbol>(r)) {
            r = r / expression(expression(i->lhs()).evalf() / expression(i->rhs()));
          }
          result = relational(l, r, eq.getop());
        } else {
          // TODO: This of course does not work if the number given in radians is hidden somewhere
          // inside an expression
          result = result / expression(expression(i->lhs()).evalf() / expression(i->rhs())); // *** bug removed in 1.3.1, was lhs * rhs !!!
        }
      } else {
        result = result.subs(*i, subs_options::algebraic);
      }
      MSG_INFO(2, "Current result: " << result << endline);
    }

    return result;
  } // subst_units()

  expression Unit::canonicalize_from(const std::string &unit) { // *** changed ex to expression in 0.4
    if (conversions[unit] == numeric(0)) {
      MSG_ERROR(0, "Unknown unit " << unit << endline);
      return expression(Unit::get("unknown")); // *** changed in 1.4.1
    } else {
      MSG_INFO(3, "Unit conversion table: " << unit << " = " << conversions[unit] << endline);
      return conversions[unit];
    }
  }

  expression Unit::canonicalize(const expression &units) { // *** changed ex to expression in 0.4
  // This method only works on expressions which are units or a multiplication of units
    MSG_INFO(3, "Canonicalizing " << units << endline);
    expression u = units.expand(); // added in 0.5 to expand functions in unit definitions

    if (is_a<Unit>(u)) {
      return conversions[ex_to<Unit>(u).name];
    } else if (is_a<mul>(u)) {
      expression result(numeric(1));
      //mul m = ex_to_mul(u);
      for (const_iterator i = u.begin(); i != u.end(); i++) {
        // *** changed to const_iterator in 0.7
        result =  result * Unit::canonicalize(*i);
      }
      return result;
    } else if (is_a<power>(u)) {
      power p = ex_to<power>(u);
      return pow(Unit::canonicalize(get_basis(p)), get_exp(p)); // *** changed to use get_... in 1.2
    } else if (is_a<symbol>(u)) {
      // This exception exists to allow for symbols that are constants, i.e. \pi. Of course, it is
      // impossible to avoid allowing other symbols that should not occur in a unit definition!
      return u;
    } else if (!is_a<numeric>(u) && !is_a<constant>(u)) { // *** added constant in 0.9
        MSG_WARN(0, "Warning: Cannot canonicalize " << u << "!" << endline);
    }
    return u;
  } // Unit::canonicalize

  void Unit::init() {
    if (initialized) return;
    
    MSG_INFO(2, "Initializing units" << endline);
    // *** Moved here in 1.3.1
//    const char *bname[8] = {"m", "kg", "s", "A", "K", "XC", "mol", "cd"}; *** removed in 1.4.1
//    base_units = std::set<const char*, ltstr>(bname,bname+8);
    const char *unames[16] = {"\\metre", "\\kilogram", "\\second", "\\ampere",
                             "\\kelvin", "\\celsius", "\\mole", "\\candela",
                             "%metre", "%kilogram", "%second", "%ampere",  // *** added in 1.3.1 for iMath
                            "%kelvin", "%celsius", "%mole", "%candela"};
    unitnames = std::set<const char*, ltstr>(unames,unames+16);

#ifdef CPPU_ENV
    conversions = std::map <std::string, expression>();
#endif

    // Initialize the map of conversions
    // Called during static initialization of class unit_initializer *** changed in 0.5                                       
    MSG_INFO(3, "Initializing unit conversion table..." << endline);
    // base SI units
    // Note that °C must be a base unit because otherwise \eq{\theta = 20\celsius} will print as \theta = 20\K !
    // *** change Unit() to Unit::get() in 1.4.1
/*    metre = Unit::getp("metre"); *** removed in 1.4.1
    kilogram = Unit::getp("kilogram");
    second = Unit::getp("second");
    ampere = Unit::getp("ampere");
    kelvin = Unit::getp("kelvin");
    celsius = Unit::getp("celsius");
    mole = Unit::getp("mole");
    candela = Unit::getp("candela"); */
    
    // added in 1.3.1 for iMath
    // *** changed to Unit::getp in 1.4.1
    Unit::getp("metre")->printname = "m"; // set_printname() is not what we want because "m" is added as a unit name and then we can't use a variable "m" any more
    Unit::getp("kilogram")->printname = "kg"; // two-letter variable names can't occur in latex so this is not a problem
    Unit::getp("second")->printname = "s";
    Unit::getp("ampere")->printname = "A";
    Unit::getp("kelvin")->printname = "K";
    Unit::getp("celsius")->printname = "\\text{\\textdegree{C}}"; // TODO: This is not UTF-8 - safe ?!
    Unit::getp("mole")->printname = "mol";
    Unit::getp("candela")->printname = "cd";
    
    // *** changed to Unit::get in 1.4.1
    conversions["metre"] = expression(Unit::get("metre"));
    conversions["kilogram"] = expression(Unit::get("kilogram"));
    conversions["second"]  = expression(Unit::get("second"));
    conversions["ampere"]  = expression(Unit::get("ampere"));
    conversions["kelvin"]  = expression(Unit::get("kelvin"));
    conversions["celsius"] = expression(Unit::get("celsius")); // we can't have a Latex macro \°C !
    conversions["mole"]= expression(Unit::get("mole"));
    conversions["candela"] = expression(Unit::get("candela"));
    
    MSG_INFO(2, "Finished initializing units" << endline);
    initialized = true;
  } // Unit::init()

  void Unit::print(std::ostream &os) {
    // *** changed ex to expression in 0.4
    std::map<std::string, expression>::iterator conv = conversions.begin();

    os << "List of unit conversions:" << std::endl;
    while (conv != conversions.end()) {
      os << conv->first << ": " << conv->second << std::endl;
      conv++;
    }
    os << std::endl;
  }
 /* 
  void Unit::print(message &ms) {
    std::map<std::string, expression>::iterator conv = conversions.begin();

    ms << "List of unit conversions:" << endline;
    while (conv != conversions.end()) {
      ms << conv->first << ": " << conv->second << endline;
      conv++;
    }
    ms << endline;
  }
*/
  void Unit::clear() {
    conversions.clear();
    // *** changed to Unit::get in 1.4.1 
    conversions["metre"] = expression(Unit::get("metre"));
    conversions["kilogram"] = expression(Unit::get("kilogram"));
    conversions["second"]  = expression(Unit::get("second"));
    conversions["ampere"]  = expression(Unit::get("ampere"));
    conversions["kelvin"]  = expression(Unit::get("kelvin"));
    conversions["celsius"] = expression(Unit::get("celsius")); // we can't have a Latex macro \°C !
    conversions["mole"]= expression(Unit::get("mole"));
    conversions["candela"] = expression(Unit::get("candela"));

/*    conversions["metre"]  = *metre;
    conversions["kilogram"] = *kilogram;
    conversions["second"]  = *second;
    conversions["ampere"]  = *ampere;
    conversions["kelvin"]  = *kelvin;
    conversions["celsius"] = *celsius;
    conversions["mole"]= *mole;
    conversions["candela"] = *candela; */
    unitnames.clear();
    const char *unames[16] = {"\\metre", "\\kilogram", "\\second", "\\ampere", // *** changed in 1.3.1
                             "\\kelvin", "\\celsius", "\\mole", "\\candela",
                             "%metre", "%kilogram", "%second", "%ampere", // *** added in 1.3.1
                             "%kelvin", "%celsius", "%mole", "%candela"};
    unitnames = std::set<const char*, ltstr>(unames,unames+16);
  }