Julia fn clearer than C++/Fortran (examples sought)

As part of an explanation of Julia’s import in facilitating & speeding the process of bringing someone on into an ongoing larger project and getting them up to speed, I am looking for examples of the same function written in C++ or Fortran (recent versions of each) and written in Julia.

Helpful examples would be where the C++/Fortran is “scary” or not at all clear without prior implementation familiarity, while the Julia version is as a walk in the park on a lovely Spring day. It is helpful if the longest versions are displayable as half-page (top or bottom half) float or in less space. Excerpts from larger functions are ok, if the excerpted parts are purpose matched and the Julia version is self-evident.

Thanks.

4 Likes

Below that post of mine, Chris Elrod posted another version of the same thing which is more elegant (in both languages), although the Fortran implementation in his case is directly thought to be interfaced with Julia, so some syntax is not natural for Fortran in general.

What I find scary about practical C++ code mostly not the functions, but the class boilerplate.

Which, BTW, is completely absent from Julia :wink:

9 Likes

@JeffreySarnoff Whatever example you pick, please do not use Fortran as a punchbag by presenting some archaic F77 or F90 code, just because so many people are unaware of modern Fortran syntax (I believe the example provided in the comments by @lmiq suffers from a similar syndrome) Use Fortran 2018 syntax, otherwise, such unfair language comparisons are examples of scientific dishonesty, in my opinion. Also, please learn about and consider all Fortran/C++ compiler optimization flags before presenting a benchmark.

7 Likes

One particular strength of Julia compared to Fortran or C is that once you have written your function, you can AD thorough it, or perhaps send Sympy symbols through it. I use this “feature” daily and my productivity would be tremendously impacted by not having it.

A simple thing like switching from Float64 to Float32 or BigFloat is also something incredibly powerful.

3 Likes

Yes, it does. That is why I mentioned the rewrite of Chris Elrod, which is much more aware of the modern Fortran syntax, and he is also the one person to provide the correct compiler flags.

One thing to keep in mind, though, is that a lot of people using Fortran use Fortran77 and Fortran90 syntax yet, and comparing the examples fairly will not only showcase Julia, but also be a lesson on the newer Fortran syntax for these people. One might wander if these people (like me) would prefer the effort to learn a new syntax of Fortran or move to another language.

3 Likes

Could you expand on “switching from Float64 to Float32 or BigFloat” feature?

Something tells me this will take some work. Should we encourage side discussions/open a repo to make sure things are fairly compared and people are working together?

I have created a repository for this purpose (with the understanding that it may be transferred at a later date). Any member of the community who would like commit privileges can have them – just ask me on this thread https://github.com/JeffreySarnoff/ComparisonsWithJulia/issues/1

1 Like

I would not do that. Old fortran code is not applicable here.
These are intended for visual comparison, others with real (modern) Fortran/C++ expertise have worked through and published careful comparative benchmarking results. One focused on communication: http://hunoldscience.net/paper/julia_hpc_paper_sahu_2020.pdf. If you know of others, do share.

2 Likes

An example I posted on Facebook a while ago to show how nice Julia is:

C++

#include <boost/random/poisson_distribution.hpp>
#include <boost/random/variate_generator.hpp>
#include <boost/random/mersenne_twister.hpp>
#include <vector>

using boost;
using std;

void create_poisson_vec(vector<int> & result)
    {
    mt19937 gen; 
    poisson_distribution<int> pdist(0.5);
    variate_generator<mt19937, poisson_distribution<int> > rvt(gen, pdist);

    for (int i=0; i<10; i++)
        result.push_back(rvt());
    }

(this might not be idiomatic, current C++ - it has been a while)

Julia:

using Distributions

function create_poisson_vec(result)
    pdist = Poisson(0.5)
    append!(result, rand(pdist, 10))
end

Edit: corrections, functions, typo

2 Likes

Would someone take a look at the foregoing C++ and let us know if alterations should be made to reflect current use and practice?

I think sampling from the poisson distribution can nowadays be done with the C++ std. Not that this any less convoluted:

https://en.cppreference.com/w/cpp/numeric/random/poisson_distribution

3 Likes

A nice example just popped up here :slight_smile:

1 Like

@Juliet I didn’t put in that thread how it looked in the end

Writing it as the following meant that the solver couldn’t handle some very rapid growth on v.

uih = [v0;β0;ϕ0;p0]
tspan = (0.0,3000.0)
prob = ODEProblem(geometria,big.(uih),tspan)
sol = DifferentialEquations.solve(prob,alghints=[:stiff],reltol=1e-6);

Writing it like this solves all the issues.

uih = [v0;β0;ϕ0;p0]
tspan = (big(0.0),big(3000.0))
prob = ODEProblem(geometria,big.(uih),tspan)
sol = DifferentialEquations.solve(prob,alghints=[:stiff],reltol=1e-6);
2 Likes

The key here being that you could easily switch the type of the parameters in your problem to BigFloat and then the ODE library was able to solve the problem with a higher precision. An ODE library written in fortran would have to provide explicit support for solving problems with BigFloat, whereas in julia this comes for free.

3 Likes

looks fine(not knowing the API’s) to me but is missing an end bracket to close the function.

It’s pretty common for people to make a header and a cpp file - but the example still stands.

About to put a series of a few simple examples up (2 files + 1 small written exposition).

Scroll down a bit, it’s hidden due to limited text box size. In any case my example still uses boost and as @simonschoelly said this stuff is in std nowadays, so I (or someone else) should probably adapt it accordingly.

1 Like

Calculating an estimate for \pi by choosing random points in a square and counting the points in the enclosed circle:

#include <iostream>
#include <random>

std::mt19937_64 rng;
std::uniform_real_distribution<double>distUReal;

using namespace std;

double calc_pi(int M) {
    int count = 0;
    double x;
    double y;
    for (int i = 0; i < M; ++i)
    {
        x = distUReal(rng);
        y = distUReal(rng);
        count += x*x + y*y < 1;
    }
    return 4.0 * count / M;
}

int main()
{
    int M = pow(10,9);
    cout << calc_pi(M) << endl;
}

vs

calc_pi(M) = 4 * count(rand()^2 + rand()^2 < 1 for i=1:M) / M

println(calc_pi(10^9))

On my laptop, Julia is even faster.

5 Likes

Defining the comparison operator in C++ looked scary to me. But after the spaceship operator in C++20 it could be an example that it’s actually easier in C++…

3 Likes