UAPUtilities.cpp

/*
 * UAP Utilities
 * Universal Accelerator Parser
 * Copyright (C) 2005, 2006 Andy Wolski, Daniel Bates, David Sagan
 * 
 * This library is free software; you can redistribute it and/or modify it
 * under the terms of the GNU Lesser General Public License as published by
 * the Free Software Foundation; either version 2.1 of the License, or (at
 * your option) any later version. 
 * 
 * This library 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 Lesser General Public License
 * for more details. 
 * 
 * You should have received a copy of the GNU Lesser General Public License
 * along with this library; if not, write to the Free Software Foundation, Inc., 
 * 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 * 
 * Direct questions, comments, etc to:
 * Andy Wolski (a.wolski@dl.ac.uk), Daniel Bates (dbates@lbl.gov), David Sagan (dcs16@cornell.edu)
 */

#include "UAPUtilities.hpp"

#include "AML/AMLUtilities.hpp"
#include "antlr/ANTLRException.hpp"
#include <algorithm>
#include "ExpressionLexer.hpp"
#include "ExpressionParser.hpp"
#include "ExpressionTreeWalker.hpp"
#include <iostream>
#include "StringInputBuffer.hpp"
#include "UAPNode.hpp"

using namespace std;

#if defined(_MSC_VER)
#define snprintf _snprintf
#endif

UAPUtilities::UAPUtilities(UAPNode* _UAPRootNode) : UAPRootNode(_UAPRootNode), seed(1L) 
{
    constants.insert(std::make_pair("pi", 3.14159265358979323844));
    constants.insert(std::make_pair("e", 2.71828182845904523536));
    constants.insert(std::make_pair("twopi", 6.28318530717958647688));
    constants.insert(std::make_pair("fourpi", 12.56637061435917295376));
    constants.insert(std::make_pair("sqrt_2", 1.41421356237309504880));
    constants.insert(std::make_pair("degrees", 0.01745329251994329576));
    constants.insert(std::make_pair("Hz", 6.28318530717958647688));
    constants.insert(std::make_pair("m_electron", 0.51099906e6));
    constants.insert(std::make_pair("m_proton", 0.938271998e9));
    constants.insert(std::make_pair("c_light", 2.99792458e8));
    constants.insert(std::make_pair("r_e", 2.8179380e-15));
    constants.insert(std::make_pair("r_p", 1.5346980e-18));
    constants.insert(std::make_pair("e_charge", 1.6021892e-19));
    constants.insert(std::make_pair("h_planck", 6.626196e-34));
    constants.insert(std::make_pair("h_bar_planck", 1.054591e-34));
}

UAPUtilities::~UAPUtilities() { }

void UAPUtilities::setUAPRootNode(UAPNode* _UAPRootNode)
{
  UAPRootNode = _UAPRootNode;
}

UAPNode* UAPUtilities::getUAPRootNode()
{
  return UAPRootNode;
}

int UAPUtilities::getParameterValue(UAPNode* node, string& param_name, string& param_value)
{
    if (!node)
        return 0;
    
  for(NodeVecIter it=node->getChildren().begin(); it!=node->getChildren().end(); ++it)
    if(getParameterValue(*it, param_name, param_value))
      return -1;

  if((node->getName() == "parameter" || node->getName()=="controller") && node->getAttributeString("name") == param_name) {
    param_value = node->getAttributeString("value");
    return -1;
  } else if (node->getName() == "constant" && node->getAttributeString("name") == param_name) {
        param_value = node->getAttributeString("design");
        return -1;
    }

  return 0;
}

int UAPUtilities::evaluate()
{
  return evaluate(UAPRootNode);
}

int UAPUtilities::evaluate(UAPNode* node, UAPNode* context)
{
  list<UAPException*> exceptions;
  
  if (!AMLUtilities::inAMLNameSpace(node))
    return -1;
  
  if (node->getName() == "element")
    context = node;
    
  char new_xpr[20];
  UAPAttribute* attrib;
  try{
    attrib = node->getAttribute("design");
    if (attrib) {
      snprintf(new_xpr,sizeof(new_xpr),"%.8g",evaluate(attrib->getValue(),context));
      attrib->setValue(new_xpr);
    }
    attrib = node->getAttribute("err");
    if (attrib) {
      snprintf(new_xpr,sizeof(new_xpr),"%.8g",evaluate(attrib->getValue(),context));
      attrib->setValue(new_xpr);
    }
    attrib = node->getAttribute("value");
    if (attrib) {
      snprintf(new_xpr,sizeof(new_xpr),"%.8g",evaluate(attrib->getValue(),context));
      attrib->setValue(new_xpr);
    }
    attrib = node->getAttribute("coef");
    if (attrib) {
      snprintf(new_xpr,sizeof(new_xpr),"%.8g",evaluate(attrib->getValue(),context));
      attrib->setValue(new_xpr);
    }
  } catch (UAPParserException& e) {
        exceptions.push_back(new UAPException(e));
  }

  for(NodeVecIter it=node->getChildren().begin(); it!=node->getChildren().end(); ++it)
    evaluate(*it, context);
  
  context = NULL;
    
    for (list<UAPException*>::const_iterator itt = exceptions.begin(); itt != exceptions.end(); ++itt){
        (*itt)->printStackTrace();
        delete *itt;
    }
    
  return -1;
}

double UAPUtilities::evaluate(UAPAttribute& attrib, UAPNode* context) throw(UAPParserException)
{
  return evaluate(attrib.getValue(), context);
}

double UAPUtilities::evaluate(const string attrib, UAPNode* context) throw(UAPParserException)
{
    double result = 0;
    UAPException* missingParams = NULL;
    string end_marker = ";";
    string xpr = attrib;
    try {
        int nsymbols = -1;
        for(int pass=0; pass<100; ++pass){
            string istr = xpr+end_marker;
            StringInputBuffer ib(istr);
            ExpressionLexer xprLexer(ib);
            ExpressionParser xprParser(xprLexer);
            xprParser.setFilename(attrib);
            nsymbols = xprParser.expr();
            RefAST xprt = RefAST(xprParser.getAST());
            ExpressionTreeWalker xprWalker;
            xprWalker.setUtilities(this);
            xprWalker.setMathConstants(constants);
            xprWalker.setContext(context);
            if(nsymbols > 0){
                string new_xpr = xprWalker.subs_params(xprt);
                int fail1 = new_xpr.find("#") + 1;
                if(fail1 > 0){
                    while(fail1 > 0){
                        int n = new_xpr.find("@",fail1) - fail1;
                        if (missingParams) {
                            missingParams = new UAPParameterNotFoundException("error: '"+new_xpr.substr(fail1,n)+"' undeclared", *missingParams);
                        } else 
                            missingParams = new UAPParameterNotFoundException("error: '"+new_xpr.substr(fail1,n)+"' undeclared");
                        fail1 = new_xpr.find("#",fail1) + 1;
                    }
                }
                if (missingParams)
                    throw UAPException(*missingParams);
                xpr = new_xpr;
            } else {
                result = xprWalker.evaluate(xprt); 
                break;
            }
            if (xprParser.getIdentifiers().empty())
                // There are no constant identifiers. To avoid an infinite loop, accelerate the for-loop.
                pass += 10;
        }
        if (nsymbols > 0)
            throw UAPInvalidExpressionException("error: evalution of the constant value for '"+attrib+"' involves a circular definition.");
    } catch (UAPException& e) {
        if (missingParams)
            delete missingParams;
        throw UAPParserException("error: failed to evaluate: "+attrib, e);
    }
    return result;
}

void UAPUtilities::addConstant(const string& name, double value, bool allow_overwrite)
{
    map<string, double>::iterator it = constants.find(name);
    if (it != constants.end() && allow_overwrite)
        it->second = value;
    else
        constants.insert(std::make_pair(name, value));
        
}

int UAPUtilities::removeConstant(const string& name)
{
    map<string, double>::iterator it = constants.find(name);
    if (it == constants.end())
        return 0;
    constants.erase(it);
    return -1;
}

list<string> UAPUtilities::parseIdentifiers(const string& expr) throw(UAPParserException)
{
    string istr = expr+";";
    try {
        StringInputBuffer ib(istr);
        ExpressionLexer xprLexer(ib);
        ExpressionParser xprParser(xprLexer);
        xprParser.setFilename(expr);
        int nsymbols = xprParser.expr();
        return xprParser.getIdentifiers();
    } catch (UAPException& e) {
        throw UAPParserException("error: failed to parse: "+expr, e);
    }
}

double UAPUtilities::nextGaussian()
{
  return static_cast<double>(gasdev(&seed));
}

double UAPUtilities::nextDouble()
{
  return static_cast<double>(ran1(&seed));
}

float UAPUtilities::gasdev(long* idum)
{
  static int iset = 0;
  static float gset;
  float fac, rsq, v1, v2;
  
  if (*idum < 0) iset=0;
  if (iset == 0){
    do{
      v1 = 2.0*ran1(idum)-1.0;
      v2 = 2.0*ran1(idum)-1.0;
      rsq = v1*v1+v2*v2;
    } while (rsq >= 1.0 || rsq == 0.0);
    fac = sqrt(-2.0*log(rsq)/rsq);
    gset = v1*fac;
    iset = 1;
    return v2*fac;
  }else{
    iset = 0;
    return gset;  
  }
}

float UAPUtilities::ran1(long* idum)
{
  int j;
  long k;
  static long iy = 0;
  static long iv[NTAB];
  float temp;
  
  if (*idum <= 0 || !iy){
    if (-(*idum) < 1) *idum = 1;
    else *idum = -(*idum);
    for (j=NTAB+7;j>=0;j--){
      k = (*idum)/IQ;
      *idum = IA*(*idum-k*IQ)-IR*k;
      if (*idum < 0) *idum += IM;
      if (j < NTAB) iv[j] = *idum;
    }
    iy = iv[0];
  }
  k = (*idum)/IQ;
  *idum = IA*(*idum-k*IQ)-IQ*k;
  if (*idum < 0) *idum += IM;
  j = iy/NDIV;
  iy = iv[j];
  iv[j] = *idum;
  if ((temp = AM*iy) > RNMX) return static_cast<float>(RNMX);
  else return temp;
}

---