Exactly that is the solution JDC gives. I guess this is the formula should be used always in computational calcs to trust the output.

Exactly that is the solution JDC gives. I guess this is the formula should be used always in computational calcs to trust the output.

- Cuchulainn
**Posts:**64109**Joined:****Location:**Drosophila melanogaster-
**Contact:**

Well, there may be even other solutionsExactly that is the solution JDC gives. I guess this is the formula should be used always in computational calcs to trust the output.

1. iterative (Newton Raphson, bisection)

2. Multi-precision arithmetic in Boost C++

3. Yet another clever closed solution

Sometimes we get lucky, but 99.99% of problems don' have a _computable_ closed solution.

- Cuchulainn
**Posts:**64109**Joined:****Location:**Drosophila melanogaster-
**Contact:**

Maybe you can modify the algorithm to find roots by avoiding (catastrophic) cancellation, etc. After all, the formula was invented in an era when round-off error was in the future

https://en.wikipedia.org/wiki/Completing_the_square

e.g. root x1 of OK then find nasty x2 by

ax^2 + bx + c = (x - x1)(x - x2)

https://en.wikipedia.org/wiki/Completing_the_square

e.g. root x1 of OK then find nasty x2 by

ax^2 + bx + c = (x - x1)(x - x2)

- Cuchulainn
**Posts:**64109**Joined:****Location:**Drosophila melanogaster-
**Contact:**

True.It is much more "complex" to use it with a computer code. Suppose you have a C++ program which output some "cost" function and you want compute the derivative of this cost function to do some optimisation problem with respect to the cost function.Computing derivatives without the pain of 1) differentiation, 2) catastrophic round-off errors. Sample code.

The scalar (almost a 1-liner). I leave the vector case as an exercise.

Code: Select all`// TestComplexStep.cpp // // Complex-step method to compute approximate derivatives. // Example is scalar-valued function of a scalar argument. // // https://pdfs.semanticscholar.org/3de7/e8ae217a4214507b9abdac66503f057aaae9.pdf // // http://mdolab.engin.umich.edu/sites/default/files/Martins2003CSD.pdf // // (C) Datasim Education BV 2018 // #include <functional> #include <complex> #include <iostream> #include <iomanip> #include <cmath> // Notation and function spaces using value_type = double; template <typename T> using FunctionType = std::function < T(const T& c)>; using CFunctionType = FunctionType<std::complex<value_type>>; // Test case from Squire&Trapp 1998 template <typename T> T func(const T& t) { T n1 = std::exp(t); T d1 = std::sin(t); T d2 = std::cos(t); return n1 / (d1*d1*d1 + d2*d2*d2); } template <typename T> T func2(const T& t) { // Derivative of e^t, sanity check return std::exp(std::pow(t,1)); // return std::exp(std::pow(t, 5)); } value_type Derivative(const CFunctionType& f, value_type x, value_type h) { // df/dx at x using tbe Complex step method std::complex<value_type> z(x, h); // x + ih, i = sqrt(-1) return std::imag(f(z)) / h; } int main() { // Squire Trapp double x = 1.5; double h = 0.1; do { std::cout << std::setprecision(12) << Derivative(func<std::complex<value_type>>, x, h) << '\n'; h *= 0.1; } while (h > 1.0e-300); // Exponential function (101 sanity check) x = 5.0; h = 1.0e-10; std::cout << "Exponential 1: " << std::setprecision(12) << Derivative(func2<std::complex<value_type>>, x, h) << '\n'; return 0; }`

You have no choice but complexify your code, e.g define a new class (complex_type) and all the functions that are required for the complex-step method. Replace all "double" with "complex_type" and some cosmetics dealing with IO.

One could write templated code and then instantiate it for compex rtype. This is a potential Pandora's box. for example,I can compute exact BS greeks but this necessitates using Faddeeva's function for N(z) with z complex.

For general functions of several complex variables I compute the gradient and Jacobian without too much effort. The only issue is the fact that I permute the value of the complex arguments in a loop. Maybe there is a more clever way?..

At some stage I would like to test it with (S)GD,

- Cuchulainn
**Posts:**64109**Joined:****Location:**Drosophila melanogaster-
**Contact:**

I think I can live with that.

It is much more "complex" to use it with a computer code. Suppose you have a C++ program which output some "cost" function and you want compute the derivative of this cost function to do some optimisation problem with respect to the cost function.

You have no choice but complexify your code, e.g define a new class (complex_type) and all the functions that are required for the complex-step method. Replace all "double" with "complex_type" and some cosmetics dealing with IO.

Question: It's just during this complexification process the usual inner product in real Euclidean space must become the Hermitian inner product, yes?

- FaridMoussaoui
**Posts:**507**Joined:****Location:**Genève, Genf, Ginevra, Geneva

per definition of the inner product for complex numbers.

- Cuchulainn
**Posts:**64109**Joined:****Location:**Drosophila melanogaster-
**Contact:**

- FaridMoussaoui
**Posts:**507**Joined:****Location:**Genève, Genf, Ginevra, Geneva

if you want to test C++ new features, switch to GNU C++. Concepts are already there.

Use any C++ Compiler with Visual Studio

Use any C++ Compiler with Visual Studio

- Cuchulainn
**Posts:**64109**Joined:****Location:**Drosophila melanogaster-
**Contact:**

Small remark on new(?) C++ behaviour.

Seems that the C++ compiler reacts differently these days to*template protected data t*; if you want to directly access t in derived class you must use this->t (BTW for concrete types you can get away with just using t).

// Of course, protected data is bad practice but is useful for testing with many derived classes

Here's the template code

Seems that the C++ compiler reacts differently these days to

// Of course, protected data is bad practice but is useful for testing with many derived classes

Here's the template code

Code: Select all

```
template <typename T>
class B
{
protected:
T t;
public:
B() {}
B(T n) : t(n) {}
};
template <typename T>
class D1: public B<T>
{
private:
public:
D1(T n) : B<T>(n)
{
this->t = 10;
}
void whatever()
{
std::cout << t << '\n'; //NOT WORK, identifier t not found
// std::cout << this->t << '\n'; // WORKS
}
};
```

Last edited by Cuchulainn on December 22nd, 2019, 9:42 pm, edited 2 times in total.

- Cuchulainn
**Posts:**64109**Joined:****Location:**Drosophila melanogaster-
**Contact:**

No problems if protected data is concrete

Code: Select all

```
class B2
{
protected:
int t;
public:
B2() {}
B2(int n) : t(n) {}
};
class D2: public B2
{
private:
public:
D2(int n) : B2(n)
{
t = 10; // Works
}
virtual void whatever()
{
std::cout << t << '\n'; // WORKS
}
};
```

Last edited by Cuchulainn on December 22nd, 2019, 9:42 pm, edited 1 time in total.

- Cuchulainn
**Posts:**64109**Joined:****Location:**Drosophila melanogaster-
**Contact:**

A nice variant (not) is when 'int t' is both base and derived classes.

- Cuchulainn
**Posts:**64109**Joined:****Location:**Drosophila melanogaster-
**Contact:**

I am solving a 2-factor PDE(x,y,t) using boost::odeint (MOL). The default matrix is boost::ublas but it is incredibly slow.

The question is how to get it working with 1) home-grown minimalist matrix, 2) Eigen.

It is not clear from the documentation how to proceed.

The question is how to get it working with 1) home-grown minimalist matrix, 2) Eigen.

It is not clear from the documentation how to proceed.

- Cuchulainn
**Posts:**64109**Joined:****Location:**Drosophila melanogaster-
**Contact:**

What are your first impressions of a piece of code I am reading? are the anti-patterns here harmful? It's a kind of santy check question:

Also, in Sewep.cpp there is a function foobar() 1) 218 lines of code 2) 4 -deep nested for lops 3) *many* nested if-else conditionals 4) OMP loop parallelisation 5) unreadable 6) it was written 1 week before going into production.

I've got my own views, but what about the reliability, efficiency, readability and maintainability of this code. On a scale [1, 10] I would give

reliability, 4

efficiency, 5?

readability 1

maintainability 2

I value feedback from the members here.

And what kind of mental makeup do you need to have to write this kind of code.

Also, in Sewep.cpp there is a function foobar() 1) 218 lines of code 2) 4 -deep nested for lops 3) *many* nested if-else conditionals 4) OMP loop parallelisation 5) unreadable 6) it was written 1 week before going into production.

I've got my own views, but what about the reliability, efficiency, readability and maintainability of this code. On a scale [1, 10] I would give

reliability, 4

efficiency, 5?

readability 1

maintainability 2

I value feedback from the members here.

And what kind of mental makeup do you need to have to write this kind of code.

Did you mean to include a link there, or are you calling for mind reading?What are your first impressions of a piece of code I am reading? are the anti-patterns here harmful? It's a kind of santy check question:

Also, in Sewep.cpp there is a function foobar() 1) 218 lines of code 2) 4 -deep nested for lops 3) *many* nested if-else conditionals 4) OMP loop parallelisation 5) unreadable 6) it was written 1 week before going into production.

I've got my own views, but what about the reliability, efficiency, readability and maintainability of this code. On a scale [1, 10] I would give

reliability, 4

efficiency, 5?

readability 1

maintainability 2

I value feedback from the members here.

And what kind of mental makeup do you need to have to write this kind of code.

- Cuchulainn
**Posts:**64109**Joined:****Location:**Drosophila melanogaster-
**Contact:**

e.g.

https://github.com/mrc-ide/covid-sim

e.g. Sweep.cpp it changes regularly .. I have already seen 3 versions..

Function DigitalContactTracing() is a ferberg (aka ball of mud).

*With respect ... the suggestions you've been making are all good ones, but I suspect that they're not that useful to the development team - who as has been mentioned before on various threads here are up to their eyeballs in urgent pandemic-related tasks: for better or worse, they simply don't have time themselves for significant refactoring or enhancements of the model. For that reason, they probably won't be acted on.*

Heard this a million times; the question is when this prototype software will be abandoned..

https://github.com/mrc-ide/covid-sim

e.g. Sweep.cpp it changes regularly .. I have already seen 3 versions..

Function DigitalContactTracing() is a ferberg (aka ball of mud).

Heard this a million times; the question is when this prototype software will be abandoned..