Calculus for C++, a framework for manipulating real functions.

[SourceFORGE Project Main Page][downloads][CVS]

Welcome to the home of calculus-cpp, a library for intuitive functional manipulations in C++. This is NOT a LISP-like C++ engine, but a replacement for <math.h> which supports mathematical associativity and all sorts of manipulations. To some extent, this is a function-pointer wrapping project, but is much, much more than that.

Motivation
Additional goals and non-goals
Our implementation
Sample of usage
BIO of the current implementation
Acknowledgements
Legal notices and such

1. Motivation

It is fascinating to watch engineering students new to programming make their first attempts at writing software. Most have little problem with basic constructs – branching, looping, arithmetic – since, after all, these relate to everyday life. One of their first disappointments with compiled languages is discovering that they are not in any way associative, though they do not necessarily articulate it this way. They expect variables to be bound to early expressions throughout when, in fact, such expressions are only evaluated once and variables are assigned volatile values.

There are indeed great advantages in having languages behave in the way they do which shall not be the topic of discussion here. However, one should argue that this does not preclude mathematical associativity totally. The essential problem in this case is that functions, mathematical or otherwise, are rarely first-class objects – or not objects at all. Furthermore, enabling first-class object manipulations of function for any language is no small feat since most are not designed to be extended in this way. For most of these languages, this would mean developing either an interpreting or a compiling layer. The trade-offs are considerable, either way.

2. Additional goals and non-goals

Primarily, the goal here is to produce functionality which is usable by the target audience (engineering and scientific programmers and students):

It must be easy to use as well as simple and robust.
It must play nice with functionality currently available and cherished by all.

The target audience is vast and everyone involved may have a different interest in this project:

Some will choose to use it for its manipulation capabilities.
Some will choose to use it for generating new functions at run-time.
Some will choose to use it to write generic functional methods.
Some will choose to use it for scripting scientific user interfaces.

Because the scope of our mission could explode into an impossible project, there are things we must not do:

Repeat functionality available elsewhere.
Attempt to provide a formal framework to derive mathematical forms proper.
Attempt to support all platforms, both software and hardware.
Attempt to do everything, at the expense of performance and quality.

3. Our implementation

We chose to implement the functional layer somewhere between the hardware processor and the function pointer abstraction in C++. This forced us to develop some amount of run-time compilation. Either we simply emit binary stubs to object code or produce entire encapsulated binaries for the manipulated functions. Note that most similar projects do not attempt this, where most provide class hierarchies with "operator()" overloads that evaluate all functions by an interpretation scheme.

This implementation allows these operations:

text string -> functional tree
functional tree -> text string
functional tree -> binary executable
functional tree -> functional tree

As well as the construction and manipulation of functional trees by direct programmatic means. The uses and applications are virtually countless.

4. Samples of usage

The only calls that *must* be made in order to use this package is a single call to initialize_calculus before using any functional trees and one call to uninitialize_calculus after using them. The former requires a single argument, though its value is currently ignored, and the latter requires no arguments but returns an error code whose value is current fixed at 0 (successs). You may want to encapsulate these calls into a small object declared at global scope or simply put them in your "main()" function.

The following samples guide you through the use of simple calls to leverage the framework. What you do with it is entirely up to you and, because these samples may seem too simple to be of any use, remember that you may combine them in countless ways to support complex programs.

4.1 Creating a function from programmatic calls

Variable x = "x", y = "y";
Function f = sin(x)*sin(y);

4.2 Creating a function from a string

//From a string literal
Function f = "sin(x)*sin(y)";
//From a run-of-the-mill string
char buf[1024];
scanf("%[^\n]",buf);
Function g = buf;

4.3 Differentiating a function

Funtion dfdx = f->get_partial_derivative(x);

4.4 Building a gradient vector

int num_vars = f->get_number_of_variables();
calculus::variable ** variables = f->get_variables();
Function * grad_f = new Function[num_vars];
for(int i = 0;i<num_vars;i++)
    grad_f[i] = f->get_partial_derivative(variables[i]);

4.5 Obtaining a string from a function

char * pBuffer = new char[f->to_string(NULL)+1]
f->to_string(pBuffer);
printf("f = %s\n",pBuffer);
delete [] pBuffer;

4.6 Obtaining a binary executable function

double(*pExec)(double*) = f;
double point[] = { 1.0, 1.5 };
double val = pExec(point);

4.7 Implementing a robust 2nd-order minimization method

void gradmin(
     Function f,
     double * p,
     double tolerance,
     double & value,
     unsigned int & iterations
     )
{
     //Get the number of Variables in the Function given
     int num_vars = f->get_number_of_variables();
     //Obtain all Variables into a convenient array
     calculus::variable ** variables = f->get_variables();
     //Build the gradient vector for this Function
     Function * grad_f = new Function[num_vars];
     Function * grad_div_f = new Function[num_vars];
     for(int i = 0;i<num_vars;i++) {
         grad_f[i] = f->get_partial_derivative(variables[i]);
         grad_div_f[i] = grad_f[i]->get_partial_derivative(variables[i]);
     }

     //Gradient components may be Functions of fewer (but no more) Variables
     //than the original Function so we need to rebuild local Variable/value
     //sets for each gradient vector component.
     int *gradient_num_vars = new int[num_vars],*gradient_div_num_vars = new int[num_vars];
     int **gradient_var_mask = new int*[num_vars],**gradient_div_var_mask = new int*[num_vars];
     for(int i = 0;i < num_vars;i++) {
         gradient_num_vars[i] = grad_f[i]->get_number_of_variables();
         gradient_div_num_vars[i] = grad_div_f[i]->get_number_of_variables();
         calculus::variable **ppGradientVars = grad_f[i]->get_variables(),**ppGradientDivVars = grad_div_f[i]->get_variables();
         gradient_var_mask[i] = new int[gradient_num_vars[i]];
         gradient_div_var_mask[i] = new int[gradient_div_num_vars[i]];
         for(int j = 0,k = 0;(j < num_vars) && (k < gradient_num_vars[i]);j++) {
             if (ppGradientVars[k] == variables[j])
                 gradient_var_mask[i][k++] = j;
         }
         for(int j = 0,k = 0;(j < num_vars) && (k < gradient_div_num_vars[i]);j++) {
             if (ppGradientDivVars[k] == variables[j])
                 gradient_div_var_mask[i][k++] = j;
         }
     }

     //More variable allocations which I will need in the iteration
     double delta_value = 0,
     *local_p = new double[num_vars],
     *grad_value = new double[num_vars],
     *delta = new double[num_vars],
     *grad_div_value = new double[num_vars];

     for(int i = 0;i < num_vars;i++)
         local_p[i] = grad_div_value[i] = grad_value[i] = delta[i] = 0;
     //Get the local Function value
     value = f(p);
     iterations = 0;
     //Iterate this section to get progressively closer to the minimum
     set_up_report_info(num_vars,f,grad_f);
     do
         {
         report_info(num_vars,iterations,value,p,grad_value,grad_div_value,delta,delta_value);
         //Get the local gradient value
         for(int i = 0;i < num_vars;i++) {
             for(int j = 0;j < gradient_num_vars[i];j++)
                 local_p[j] = p[gradient_var_mask[i][j]];}
             grad_value[i] = grad_f[i](local_p);
             for(int j = 0;j < gradient_div_num_vars[i];j++)
                 local_p[j] = p[gradient_div_var_mask[i][j]];
             grad_div_value[i] = grad_div_f[i](local_p);
         }
         //Get the various directional deltas
         for(int i = 0;i < num_vars;i++)
             delta[i] = -grad_value[i]/fabs(grad_div_value[i]);}
         //Get the next point value guess and determine how much motion resulted
         delta_value = 0;
         for(int i = 0;i < num_vars;i++) {
             p[i] += delta[i];
             delta_value += delta[i]*delta[i];
         }
         delta_value = sqrt(delta_value);
         //Get the local Function value
         value = f(p);
         //Increment iteration count
         if (++iterations > 1000)
             break;
     }
     while(fabs(delta_value) > tolerance);

     //Free all memory allocated
     for(int i = 0;i < num_vars;i++) {
         delete [] gradient_div_var_mask[i];
         delete [] gradient_var_mask[i];
     }
     delete [] local_p;
     delete [] gradient_div_var_mask;
     delete [] gradient_var_mask;
     delete [] gradient_div_num_vars;
     delete [] gradient_num_vars;
     delete [] delta;
     delete [] grad_div_value;
     delete [] grad_value;
     delete [] grad_div_f;
     delete [] grad_f;
}

4.8 Implementing the "subst" functionality from scratch

void subst(
     Function & f,
     Function & ffind,
     Function & freplace
     )
{
     if (is_equal(f,ffind)) {
          f = freplace;
     }
     else
          subst_inner(f,ffind,freplace);
}

void subst_inner(
     calculus::algebraic_operator * pao_f,
     calculus::algebraic_operator * pao_ffind,
     calculus::algebraic_operator * pao_freplace
     )
{
     calculus::unary_operators::unary_operator* puo_f = dynamic_cast<calculus::unary_operators::unary_operator*>(pao_f);
     if (puo_f) {
          if (is_equal(puo_f->get_operand(),pao_ffind)) {
               puo_f->set_operand(pao_freplace);
          }
         else
              subst_inner(puo_f->get_operand(),pao_ffind,pao_freplace);
     return;
     }

     calculus::binary_operators::binary_operator* pbo_f = dynamic_cast<calculus::binary_operators::binary_operator*>(pao_f);
     if (pbo_f) {
          if (is_equal(pbo_f->get_left_operand(),pao_ffind)) {
               pbo_f->set_left_operand(pao_freplace);
          }
          else
               subst_inner(pbo_f->get_left_operand(),pao_ffind,pao_freplace);

          if (is_equal(pbo_f->get_right_operand(),pao_ffind)) {
               pbo_f->set_right_operand(pao_freplace);
          }
          else
               subst_inner(pbo_f->get_right_operand(),pao_ffind,pao_freplace);
     }
}

5. BIO of the current implementation

This software was developed at McGill University (Montreal, 2002) by Olivier Giroux in the course of his studies in Mechanical Engineering. It was presented, along with an accompanying paper, for credit in the fall of 2002. Calculus-cpp was not designed to prove a point or to serve as a formal framework within which exact solutions can be derived. Instead it was created to fill the need for run-time Functional constructions and to accomplish very real and tangible goals. It remains your responsibility to use it properly - as much more sophisticated , which allows Functions to be treated as first-class objects. You are welcome to make any additions you feel are necessary.

6. Acknowledgements

SourceFORGE.net:

This project relies on the generous support of SourceFORGE.net.

7. Legal notices and such

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without any restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, provided that the copyright notice(s) and this permission notice appear in all copies of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT OF THIRD PARTY RIGHTS. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR HOLDERS INCLUDED IN THIS NOTICE BE LIABLE FOR ANY CLAIM, OR ANY SPECIAL INDIRECT OR CONSEQUENTIAL DAMAGES, OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.

Except as contained in this notice, the name of a copyright holder shall not be used in advertising or otherwise to promote the sale, use or other dealings in this Software without prior written authorization of the copyright holder.

THIS SOFTWARE INCLUDES THE NIST'S TNT PACKAGE AND THE LEMON PARSER GENERATOR.

THE FOLLOWING NOTICE APPLIES SOLELY TO THE TNT-->
Template Numerical Toolkit (TNT): Linear Algebra Module
Mathematical and Computational Sciences Division
National Institute of Technology,
Gaithersburg, MD USA

This software was developed at the National Institute of Standards and Technology (NIST) by employees of the Federal Government in the course of their official duties. Pursuant to title 17 Section 105 of the United States Code, this software is not subject to copyright protection and is in the public domain. NIST assumes no responsibility whatsoever for its use by other parties, and makes no guarantees, expressed or implied, about its quality, reliability, or any other characteristic.
<--END NOTICE

THE FOLLOWING NOTICE APPLIES SOLELY TO LEMON-->
Copyright (c) 1991, 1994, 1997, 1998 D. Richard Hipp

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 library; if not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. Author contact information: drh@acm.org http://www.hwaci.com/drh/
<--END NOTICE