// Test001.cpp
/*
(C) Datasim Education BV 2022
2022-10-20 DJD #8
Super-strippped down sample code to model all kinds of mathematical functions that are neeeded
in numerical applications. Here we reduce the scope by introducing scalar-valued
functions of a scalar argument.
Rationale/preview:
1. Introducing some C++ features to support functional programming style.
2. (mathematically-defined) Vector spaces of scalar, vector and vector-valued functions.
3. Applications of Vector spaces in numerical computing (optimisation, PDE, ML etc. etc.)
(e.g. build complex functions from simpler functions)
4. Heterogeneous container, functions and algorithms (Boost C++ Fusian and Hana).
5. Integration with Design Patterns, e.g. delegation (stateful/stateless versions).
Ideally, this should really in some new Boost C++ library.
** No futher comments for the moment. Later variadic arguments.
*/
#include <functional>
#include <cmath>
#include <iostream>
#include <string>
#include <boost/math/tools/minima.hpp>
template <typename T>
using FType = std::function<T (T)>;
template <typename T>
FType<T> operator + (const FType<T>& f, const FType<T>& g)
{ // Addition
return [=](T x)
{
return f(x) + g(x);
};
}
template <typename T>
FType<T> operator + (const FType<T>& f, T a)
{ // Scalar addition
return [=](T x)
{
return f(x) + a;
};
}
template <typename T>
FType<T> operator * (T a, const FType<T>& f)
{ // Scalar multiplication
return [=](T x)
{
return a*f(x);
};
}
template <typename T>
FType<T> operator - (const FType<T>& f, FType<T>& g)
{ // Subtraction
return [=](T x)
{
return f(x) - g(x);
};
}
template <typename T>
FType<T> operator - (const FType<T>& f)
{ // Unary negation
return [=](T x)
{
return -f(x);
};
}
template <typename T>
FType<T> operator + (const FType<T>& f)
{ // Unary plus
return [=](T x)
{
return +f(x);
};
}
template <typename T>
FType<T> operator << (FType<T>& f, FType<T>& g)
{ // Composition of functions
return [=](T x)
{
return f(g(x));
};
}
template <typename T>
FType<T> assign(T constant)
{ // Create a function from a scalar constant
return [=](T x)
{
return constant;
};
}
template <typename T>
FType<T> exp(const FType<T>& f)
{ // exp
return [=](T x)
{
return std::exp(f(x));
};
}
template <typename T>
FType<T> Exp(const FType<T>& f)
{ // exp
return [=](T x)
{
return std::exp(f(x));
};
}
template <typename T>
FType<T> log(const FType<T>& f)
{ // log
return [=](T x)
{
return std::log(f(x));
};
}
template <typename T>
FType<T> pow(const FType<T>& f, int n)
{ // powers of a function f(x)*...*f(x) (n times)
return [=](T x)
{
return std::pow(f(x), n);
};
}
//////////////////////////////////////////////////////////////////////////////////////////////////
//
// Variadic arguments
//
template <typename RT, typename ...T>
using FTypeV = std::function<RT(T... args)>;
template <typename RT, typename ...T>
FTypeV<RT, T...> operator + (const FTypeV<RT, T...>& f, const FTypeV<RT, T...>& g)
{ // Addition
return [=](T... x)
{
return f(x...) + g(x...);
};
}
// User functions
double G(double x)
{
return x+1;
}
double H(double x, double y)
{
return x + y;
}
double AckleyGlobalFunction(double x)
{
double a = 20.0; double b = 0.2; double c = 2.0*3.14159;
return -a*std::exp(-b*x) - std::exp(std::cos(c*x)) + a + std::exp(1.0);
}
double AckleyGlobalFunctionII(double x)
{ // Building the functions in a HOF assembly process
double a = 20.0; double b = 0.2; double c = 2.0*3.14159;
// Bit unwieldy; used to show function vector spaces.
FType<double> f11 = [&](double x) {return -b*x; };
FType<double> f1 = -a*exp(f11);
FType<double> f22 = [&](double x) {return std::cos(c*x); };
FType<double> f2 = -exp(f22);
// FType<double> f2A = f2 << f1 << f11;
//FType<double> f3 = assign(a + std::exp(1.0));
FType<double> f2B = [&](double x) {return std::exp(1.0) ; };
FType<double> f3 = f2B + a;
FType<double> Ackley = f1 + f2 + f3;
return Ackley(x);
}
double TestFunction(double x)
{ // To show how function algebra works
// Bit unwieldy; used to show function vector spaces.
FType<double> f1 = [](double x) {return x*x; };
FType<double> f2 = [&](double x) {return f1(x)*f1(x); };
FType<double> f3 = [&](double x) {return f2(x)*f2(x); };
FType<double> sumFunc = f1 + f2 + f3;
return sumFunc(x);
}
double RastriginGlobalFunction(double x)
{
double A = 10.0;
return A + x*x - A*std::cos(2.0*3.14159);
}
double SphereGlobalFunction(double x)
{
return x*x;
}
// Testing a function in terms of components
double SphereGlobalFunction(double x, double y, double z)
{
return x * x + y * y + z * z;
}
double SphereX(double x, double y, double z)
{
return x * x;
}
double SphereY(double x, double y, double z)
{
return y * y;
}
double SphereZ(double x, double y, double z)
{
return z * z;
}
int main()
{
FType<double> e = [](double x) { return 5.0; };
FType<double> f = [](double x) { return G(x); };
FType<double> g = [](double x) { return x*x; };
// Currying, I -> H(double x, double y)
FType<double> bindFun = std::bind(H, 1.0, std::placeholders::_1);
auto func1 = g + bindFun;
std::cout << "13, func1? " << func1(3.0) << '\n';
// Currying, II
double y = 1.0;
FType<double> capturdFun = [&](double x) { return H(x, y); };
auto func2 = g + capturdFun;
std::cout << "13, func2? " << func2(3.0) << '\n';
auto h1 = f + g;
FType<double> h2 = f + exp(e);
FType<double> h3 = g + f;
// Experimenting
FType<double> h4 = exp(exp(f));
FType<double> h5 = exp(exp(f) + g);
FType<double> h6 = log(exp(h5) + h5);
FType<double> h7 = h6 + 1.0;
double x = 5.0;
std::cout << "1 " << h1(x) << "\n";
std::cout << "2 " << h2(x) << "\n";
std::cout << "3 " << h3(x) << "\n";
x = 0.0;
std::cout << "4 " << h4(x) << "\n";
std::cout << "5 " << h5(x) << "\n";
std::cout << "6 " << h6(x) << "\n";
std::cout << "7 " << h7(x) << "\n";
{ // Powers
FType<double> f = [](double x) { return x*x; };
int n = 4;
auto prodFun = pow(f, n);
double x = 2.0;
std::cout << "Power: " << prodFun(x) << "\n";
}
{ // Unary minus and plus
auto f = assign(5.0);
auto f1 = exp(f);
auto f2 = +exp(f);
auto f3 = -exp(f);
auto f4 = -(-exp(f));
auto f5 = +(-exp(f));
auto f6 = -(+exp(f));
double x = 5.0;
std::cout << "exp stuff: " << f1(x) << ", " << f2(x) << ", "
<< f3(x) << ", " << f4(x) << ", " << f5(x) << ", " << f6(x) << '\n';
}
FType<double> f11 = [&](double x) {return -b*x; };
FType<double> f1 = -a*exp(f11);
FType<double> f22 = [&](double x) {return std::cos(c*x); };
FType<double> f2 = -exp(f22);
FType<double> f2AB = f2 << f1;
FType<double> f2A = f2AB << f11;
FType<double> f3 = assign(a + std::exp(1.0));
FType<double> Ackley = f1 + f2 + f3;
FType<double> Ackley2 = [](double x) { return AckleyGlobalFunction(x); };
x = -1.0;
while (x < 1.0)
{ // Sanity clause: compute the value in different ways
std::cout << "Error " << x << " " <<
AckleyGlobalFunction(x) - AckleyGlobalFunctionII(x) << " ";
std::cout << "Value " << AckleyGlobalFunction(x) << ","
<< AckleyGlobalFunctionII(x) << ", " << Ackley2(x) << '\n';
x += 0.1;
}
// C. Rastrigin
double A = 10.0;
FType<double> fr1 = assign(A);
FType<double> fr2 = [&](double x) {return x*x; };
FType<double> fr3 = [&](double x) {return std::cos(2.0*3.14159*x); };
FType<double> fr4 = -A*fr3;
FType<double> func = fr1 + fr2 + fr4;
x = -1.0;
while (x < 1.0)
{
std::cout << x << ", " << RastriginGlobalFunction(x) << ", " << func(x) << '\n';
x += 0.1;
}
{ // RT x y z
FTypeV<double, double, double, double> f1 = SphereX;
FTypeV<double, double, double, double> f2 = SphereY;
FTypeV<double, double, double, double> f3 = SphereZ;
FTypeV<double, double, double, double> fSum = f1 + f2 + f3;
std::cout << "Sphere " << fSum(2, 2, 2) << '\n'; // 12
}
}