299 lines
9.8 KiB
C++
299 lines
9.8 KiB
C++
/***************************************************************************
|
|
|
|
PFunction.cpp - Test version (simplified, no PMetaData)
|
|
|
|
Author: Andreas Suter
|
|
e-mail: andreas.suter@psi.ch
|
|
|
|
Modernized test implementation using Boost.Spirit X3.
|
|
|
|
***************************************************************************/
|
|
|
|
/***************************************************************************
|
|
* Copyright (C) 2007-2026 by Andreas Suter *
|
|
* andreas.suter@psi.ch *
|
|
* *
|
|
* 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. *
|
|
* *
|
|
* 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. *
|
|
***************************************************************************/
|
|
|
|
#include <cmath>
|
|
#include <iostream>
|
|
|
|
#include <boost/spirit/home/x3.hpp>
|
|
#include <boost/variant/apply_visitor.hpp>
|
|
|
|
#include "PFunction.h"
|
|
#include "PFunctionGrammar.h"
|
|
|
|
namespace x3 = boost::spirit::x3;
|
|
|
|
//--------------------------------------------------------------------------
|
|
// Constructor
|
|
//--------------------------------------------------------------------------
|
|
PFunction::PFunction(const std::string& input, std::vector<double> param,
|
|
std::vector<int> map, bool debug)
|
|
: fParam(param), fMap(map), fDebug(debug)
|
|
{
|
|
fValid = false;
|
|
fFuncNo = -1;
|
|
|
|
// Parse the input
|
|
musrfit::ast::assignment assignment;
|
|
std::string::const_iterator iter = input.begin();
|
|
std::string::const_iterator end = input.end();
|
|
|
|
// Get the X3 grammar
|
|
auto const& grammar = musrfit::function_grammar();
|
|
|
|
// Try to parse with X3
|
|
try {
|
|
bool success = x3::phrase_parse(iter, end, grammar, x3::space, assignment);
|
|
|
|
if (success && iter == end) {
|
|
fFuncNo = assignment.func_number;
|
|
fAst = assignment.rhs;
|
|
fValid = true;
|
|
|
|
if (fDebug) {
|
|
std::cout << "Parse successful: FUN" << fFuncNo << std::endl;
|
|
}
|
|
} else {
|
|
std::cerr << "**ERROR** Failed to parse: " << input << std::endl;
|
|
if (iter != end) {
|
|
std::cerr << "**ERROR** Stopped at: " << std::string(iter, end) << std::endl;
|
|
}
|
|
fValid = false;
|
|
}
|
|
} catch (x3::expectation_failure<std::string::const_iterator> const& e) {
|
|
std::cerr << "**ERROR** Parsing failed. Expected: " << e.which() << std::endl;
|
|
std::cerr << "**ERROR** At position: " << std::distance(input.begin(), e.where()) << std::endl;
|
|
fValid = false;
|
|
}
|
|
}
|
|
|
|
//--------------------------------------------------------------------------
|
|
// Destructor
|
|
//--------------------------------------------------------------------------
|
|
PFunction::~PFunction()
|
|
{
|
|
fParam.clear();
|
|
fMap.clear();
|
|
}
|
|
|
|
//--------------------------------------------------------------------------
|
|
// Eval
|
|
//--------------------------------------------------------------------------
|
|
double PFunction::Eval()
|
|
{
|
|
if (!fValid) return 0.0;
|
|
|
|
EvalVisitor evaluator(fMap, fParam);
|
|
|
|
// Evaluate first operand
|
|
double result = boost::apply_visitor(evaluator, fAst.first);
|
|
|
|
// Apply operations
|
|
for (const auto& op : fAst.rest) {
|
|
double rhs = boost::apply_visitor(evaluator, op.operand_);
|
|
|
|
switch (op.operator_) {
|
|
case musrfit::ast::op_plus:
|
|
result += rhs;
|
|
break;
|
|
case musrfit::ast::op_minus:
|
|
result -= rhs;
|
|
break;
|
|
case musrfit::ast::op_times:
|
|
result *= rhs;
|
|
break;
|
|
case musrfit::ast::op_divide:
|
|
if (std::abs(rhs) > 1.0e-20) result /= rhs;
|
|
else result = 0.0;
|
|
break;
|
|
}
|
|
}
|
|
|
|
return result;
|
|
}
|
|
|
|
//==========================================================================
|
|
// EvalVisitor Implementation
|
|
//==========================================================================
|
|
|
|
double PFunction::EvalVisitor::operator()(const musrfit::ast::nil&) const
|
|
{
|
|
return 0.0;
|
|
}
|
|
|
|
double PFunction::EvalVisitor::operator()(double val) const
|
|
{
|
|
return val;
|
|
}
|
|
|
|
double PFunction::EvalVisitor::operator()(const musrfit::ast::constant& c) const
|
|
{
|
|
double val = 0.0;
|
|
|
|
switch (c.const_type) {
|
|
case musrfit::ast::constant::pi:
|
|
val = 3.14159265358979323846;
|
|
break;
|
|
|
|
case musrfit::ast::constant::gamma_mu:
|
|
val = 0.01355342; // GAMMA_BAR_MUON equivalent for test
|
|
break;
|
|
|
|
// Test version doesn't support metadata, so these return 0
|
|
case musrfit::ast::constant::field:
|
|
case musrfit::ast::constant::energy:
|
|
case musrfit::ast::constant::temp:
|
|
std::cerr << "**WARNING** Metadata not supported in test version" << std::endl;
|
|
val = 0.0;
|
|
break;
|
|
|
|
default:
|
|
val = 0.0;
|
|
}
|
|
|
|
return c.sign ? -val : val;
|
|
}
|
|
|
|
double PFunction::EvalVisitor::operator()(const musrfit::ast::parameter& p) const
|
|
{
|
|
int idx = p.number - 1;
|
|
|
|
if (idx < 0 || idx >= static_cast<int>(fParam.size())) {
|
|
std::cerr << "**ERROR** Parameter out of range: PAR" << p.number << std::endl;
|
|
return 0.0;
|
|
}
|
|
|
|
double val = fParam[idx];
|
|
return p.sign ? -val : val;
|
|
}
|
|
|
|
double PFunction::EvalVisitor::operator()(const musrfit::ast::map_ref& m) const
|
|
{
|
|
int mapIdx = m.number - 1;
|
|
|
|
if (mapIdx < 0 || mapIdx >= static_cast<int>(fMap.size())) {
|
|
std::cerr << "**ERROR** Map out of range: MAP" << m.number << std::endl;
|
|
return 0.0;
|
|
}
|
|
|
|
int paramIdx = fMap[mapIdx] - 1;
|
|
|
|
if (paramIdx < 0 || paramIdx >= static_cast<int>(fParam.size())) {
|
|
std::cerr << "**ERROR** Mapped parameter out of range" << std::endl;
|
|
return 0.0;
|
|
}
|
|
|
|
double val = fParam[paramIdx];
|
|
return m.sign ? -val : val;
|
|
}
|
|
|
|
double PFunction::EvalVisitor::operator()(const musrfit::ast::function_call& f) const
|
|
{
|
|
// Evaluate argument
|
|
double arg = boost::apply_visitor(*this, f.arg.first);
|
|
for (const auto& op : f.arg.rest) {
|
|
double rhs = boost::apply_visitor(*this, op.operand_);
|
|
switch (op.operator_) {
|
|
case musrfit::ast::op_plus: arg += rhs; break;
|
|
case musrfit::ast::op_minus: arg -= rhs; break;
|
|
case musrfit::ast::op_times: arg *= rhs; break;
|
|
case musrfit::ast::op_divide:
|
|
if (std::abs(rhs) > 1.0e-20) arg /= rhs;
|
|
break;
|
|
}
|
|
}
|
|
|
|
// Apply function
|
|
double result = 0.0;
|
|
switch (f.func_id) {
|
|
case musrfit::ast::fun_cos: result = std::cos(arg); break;
|
|
case musrfit::ast::fun_sin: result = std::sin(arg); break;
|
|
case musrfit::ast::fun_tan: result = std::tan(arg); break;
|
|
case musrfit::ast::fun_cosh: result = std::cosh(arg); break;
|
|
case musrfit::ast::fun_sinh: result = std::sinh(arg); break;
|
|
case musrfit::ast::fun_tanh: result = std::tanh(arg); break;
|
|
case musrfit::ast::fun_acos: result = std::acos(arg); break;
|
|
case musrfit::ast::fun_asin: result = std::asin(arg); break;
|
|
case musrfit::ast::fun_atan: result = std::atan(arg); break;
|
|
case musrfit::ast::fun_acosh: result = std::acosh(arg); break;
|
|
case musrfit::ast::fun_asinh: result = std::asinh(arg); break;
|
|
case musrfit::ast::fun_atanh: result = std::atanh(arg); break;
|
|
case musrfit::ast::fun_log: result = std::log10(std::abs(arg)); break;
|
|
case musrfit::ast::fun_ln: result = std::log(std::abs(arg)); break;
|
|
case musrfit::ast::fun_exp: result = std::exp(arg); break;
|
|
case musrfit::ast::fun_sqrt: result = std::sqrt(std::abs(arg)); break;
|
|
default: result = 0.0;
|
|
}
|
|
|
|
return result;
|
|
}
|
|
|
|
double PFunction::EvalVisitor::operator()(const musrfit::ast::power_call& p) const
|
|
{
|
|
// Evaluate base
|
|
double base = boost::apply_visitor(*this, p.base.first);
|
|
for (const auto& op : p.base.rest) {
|
|
double rhs = boost::apply_visitor(*this, op.operand_);
|
|
switch (op.operator_) {
|
|
case musrfit::ast::op_plus: base += rhs; break;
|
|
case musrfit::ast::op_minus: base -= rhs; break;
|
|
case musrfit::ast::op_times: base *= rhs; break;
|
|
case musrfit::ast::op_divide:
|
|
if (std::abs(rhs) > 1.0e-20) base /= rhs;
|
|
break;
|
|
}
|
|
}
|
|
|
|
// Evaluate exponent
|
|
double exponent = boost::apply_visitor(*this, p.pow.first);
|
|
for (const auto& op : p.pow.rest) {
|
|
double rhs = boost::apply_visitor(*this, op.operand_);
|
|
switch (op.operator_) {
|
|
case musrfit::ast::op_plus: exponent += rhs; break;
|
|
case musrfit::ast::op_minus: exponent -= rhs; break;
|
|
case musrfit::ast::op_times: exponent *= rhs; break;
|
|
case musrfit::ast::op_divide:
|
|
if (std::abs(rhs) > 1.0e-20) exponent /= rhs;
|
|
break;
|
|
}
|
|
}
|
|
|
|
return std::pow(std::abs(base), exponent);
|
|
}
|
|
|
|
double PFunction::EvalVisitor::operator()(const musrfit::ast::expression& e) const
|
|
{
|
|
double result = boost::apply_visitor(*this, e.first);
|
|
|
|
for (const auto& op : e.rest) {
|
|
double rhs = boost::apply_visitor(*this, op.operand_);
|
|
switch (op.operator_) {
|
|
case musrfit::ast::op_plus: result += rhs; break;
|
|
case musrfit::ast::op_minus: result -= rhs; break;
|
|
case musrfit::ast::op_times: result *= rhs; break;
|
|
case musrfit::ast::op_divide:
|
|
if (std::abs(rhs) > 1.0e-20) result /= rhs;
|
|
break;
|
|
}
|
|
}
|
|
|
|
return result;
|
|
}
|