Julia Community 🟣: Matthijs Cox

How to solve the two language problem?

Matthijs Cox — Mon, 08 May 2023 08:03:46 +0000

My professional obsession is solving the Two Culture Problem. How can scientists optimally join forces with software engineers and their principles, so that we can work on the same problems together? How to accelerate the cycle from idea to product? The Two Culture Problem requires a solution to the related Two Language Problem, which has a technical nature. A solution to the technical problem does not guarantee a solution to the organizational problem, but when it comes to engineering cultures you first need to prove the technical solution before you can even begin to tackle the social implications. I have a strong opinion on the best technical solution, but let's review all our options.

As far as I can tell, we have the following alternatives:

Accept the status quo: use a slow and a fast (usually harder) language
Code generation using a look-a-like framework inside the slow language
Using (LLVM-based) optimization frameworks that look like the slow language
Speed up the slow language itself, working around its limitations
Design a new language that is both easy and fast

There are many tutorials about all of these options. Here I'd like to write a short overview of all of these.

For another similar technical overview, see Martin Maas's blog posts about Julia vs Python vs Numba vs Cython.

The two language problem - Python and C++ as a primary example

Your scientists or domain experts write prototypes in a simple language, let's say Python, where they can rapidly explore, do dynamic data analysis, model desired behavior and gather requirements from users. When they find something valuable, software engineers convert the prototype into high-performant code, let's say C++, and integrate it into production systems to sell as professional services to those users. This is my assumption of the status quo.

(I will continue to refer to the modeling culture as "scientists" whether or not they are actual scientists, domain experts, requirements engineers, data analysts, quants or any other kind of expert whose primary job is modeling the behavior of your product without actually writing and deploying the final source code.)

Depending on the size of your organization and the skill level of your engineers, you may end up with several configurations:

Teams of highly skilled scientific engineers who can do all the work
Teams with a mix of scientists and software engineers
Separate teams of scientists and separate teams of software engineers

You may have any combination of the above. The first option is a team of unicorns, which I have seen the least, but is amazing to work with.

Perhaps you have accepted this status quo. As the organization grows, separate code bases may evolve for the two types of tasks. In my experience, the production code rarely gets re-integrated into the analysis code, because it's not worth the effort in the short term. Long term you may get inconsistencies and other issues, but that's typically for someone else to worry about. Or perhaps people notice the problems, but profit margins are good, so why worry?

If you integrate your fast code (C++) as embedded libraries into the slow language (Python), you typically need some intermediate glue code or language in between. This requires yet more technical expertise from your people. See for example this blog about How to Call C++ from Python.

One stated benefit of keeping the two-language culture intact is that your prototypes and your scientists never mess up your production systems. The production systems are brittle and valuable, so this is a valid concern, but I think there are better ways to teach people to write better code than by blocking them.

For this article, I assume you are looking for alternatives. Maybe the problems have grown too big, or you want to avoid them early on, or you simply cannot hire enough senior software engineers. Thus for one reason or another, you need your scientists to be deeply involved in the software development.

Learning C/C++ is still a good idea to grow your expertise or the competence of your scientists, but it can take a long time to develop. At a minimum, I advise learning what it means to compile and link libraries. And learn a bit about computers by reading great summaries such as What Scientists Should Know About Hardware to Write Fast Code.

Still want to find a technology that's easier to use, yet brings some of the hardcore software benefits? Let's see what's possible!

Code generation - Cython example

Generating low-level code, most likely C, from a high-level language, most likely Python, is typically done to try to avoid some of the disadvantages of the two-language problem. Maybe you want to compile static libraries to embed into devices. Or you have some other reason. Unless your examples are very simple, do not expect big performance boosts though, the generated code still needs to make similar kinds of assumptions as the high-level language. Also, most code generators do not support the complete language semantics, so you will have to make sure your high-level code adheres to the capabilities of the generator.

In Python you can use the Cython 'compiler' to help you generate C code. On the surface it looks a lot like Python, yet with C types and certain decorators. This means you need to rewrite the parts of your codebase that you want to speed up. The process of turning Python into Cython is sometimes called "cythonizing". You get such examples in the Cython quickstart tutorial:

@cython.cfunc
@cython.exceptval(-2, check=True)
def f(x: cython.double) -> cython.double:
    return x ** 2 - x

def integrate_f(a: cython.double, b: cython.double, N: cython.int):
    i: cython.int
    s: cython.double
    dx: cython.double
    s = 0
    dx = (b - a) / N
    for i in range(N):
        s += f(a + i * dx)
    return s * dx

You have to build this cythonized code. As I mentioned, this happens in two stages:

The .py or .pyx file is converted by Cython to a .c file.
The .c file is compiled by a C compiler to a .so file (or .pyd on Windows) which can be import-ed back into Python with setuptools.

The downside of Cython, and any similar C code generator, is that you obtain rather obscure C code. If you want code obfuscation, you can consider that a benefit, but trouble begins when you have to work with that C code. It can be hard to debug once deployed in the field. Make sure to add lots of clear error messages and logging. If you integrate the generated code inside existing C/C++ codebases, your software engineers may dislike writing the necessary glue-code (I learned that from experience). Finally, naive Cython is not very performant and writing optimized Cython can be as difficult as writing regular C code. But the benefit is that you can move gradually up in complexity.

How to make a standalone-ish distribution? Cython generates code that interfaces with the python runtime. You can create a binary executable, but you also need to distribute it with libpython.so which is the python runtime. Moreover, you also need to add all the python dependencies and .so/.dlls that those packages are using. This might be a bit tedious using Cython, but it is certainly possible. Other packages like Nuitka make this process a bit less painless by figuring out all your dependencies.

Fun fact: Code generation is sometimes referred to as "transpiling", since you translate your code to another language that's ready for compiling.

Interlude: LLVM

What if we do not want to write or generate C code? Do we have any other options? Yes, we can generate something else: LLVM code! Before we go into such frameworks, let's do a quick introduction into LLVM itself.

LLVM is a middleman between your source code and the compiled native code. Compilers typically consist of two stages: byte code and native code. The byte code is an intermediate representation that is agnostic of the CPU or GPU architecture. LLVM is an attempt to standardize the byte code definition, which will then be compiled for you to any architecture you want. In some frameworks or languages (like Julia) you can ask to see the LLVM code, and the eventual native code, the processor instructions, which is typically assembly code (that's just before it becomes those zeros and ones you always hear about).

Frameworks that use LLVM may store the compiled native code in memory. In that case, if you want to distribute the compiled code, you need the option to package the native code into a library (that's a .so on Linux or a .dll on Windows) together with all of its dependencies. This packaging option may be important to investigate for your deployment strategy.

Fun fact: the clang compiler from C also compiles via LLVM. So if you write Cython and then compile via clang, you are taking an interesting route.

Optimization Frameworks - Numba example

Numba is an LLVM code generator that integrates directly with Python. What you need to do is add decorators to every Python function you want to optimize. In principle it looks simple:

@jit(int32(int32, int32), nopython=True)
def f(x, y):
    return x + y

JIT stands for Just-In-Time, as it compiles to LLVM at the moment of calling the function, and inferring which types you used, just in time before executing. You can optionally provide the types yourself, as I did above. And there are lots of other settings for the @jit decorator, like the nopython mode to get faster performance.

Only a subset of Python is supported with Numba. It works well with your NumPy code, they made sure of that. It doesn't work with other packages such as Pandas, because those work differently. Even dictionaries are not supported by Numba. Nobody writes a custom file format parser in Numba, it's purely for numerical code.

If you want to compile to GPU: there are other decorators, such as cuda.@jit. This suggests you need to edit your code for GPU.

If you want to compile ahead of time, you will again have to replace all your @jit decorators with the @cc.export decorator and be explicit about your types.

You cannot debug the jitted Numba code itself, you'll have to change the decorator setting to debug mode and use the gdb tool, so be careful there. That's another disadvantage of Numba.

We have never tried to make a standalone distribution of Numba compiled code, but I assume you ship the entire Python environment, with the ahead-of-time compiled code. If someone has experience with the nitty-gritty details of distributing Numba code, then let me know!

Codon

Codon is a recent attempt similar to Numba, except it claims zero-overhead; you do not necessarily have to decorate your code. Well, except if you want to use it inside larger Python codebases (you probably do), then you have the @codon.jit decorator and other decorators depending on your use-case.

Codon has only 9 contributors at the moment, and it has a non-permissive license, so you'll have to pay to use Codon commercially in production. It's interesting but looks more like a startup than a regular open-source project. Similar to Numba it only supports a subset of the Python language, which may get better over time (or worse, if Python evolves, yet the developers do not update Codon).

Jax, TensorFlow, PyTorch

Every scientific computing and machine learning framework in Python implements its own optimized numerical libraries it seems. Some of them, like JAX, have a @jit decorator like Numba. All these frameworks look like Python, but to get performant code you'll have to use their API, not Python itself. Often you write Python in a more complicated directed-acyclic-graph (DAG) structure that can be fed to the underlying libraries for execution. Don't ask me how to debug these things. Please consider whether you are really writing Python or another language.

Also see this section from the Mojo language comparing such Python improvements.

Boost the slow language - PyPy

There is a continuous effort to improve the performance of slow interpreted languages like Python and R. In a blog post called Python 3.14 Will be Faster than C++ the author joked that linear extrapolation of Python improvements will soon surpass C++ performance. Let's see how that graph evolves in the next Python versions.

PyPy

An alternative to waiting for Python to improve is PyPy. This is a replacement for CPython. Note, CPython is not Cython. CPython is essentially Python itself, as the Python interpreter is written in the C language. PyPy is an attempt to make the entire Python language faster with a better interpreter. PyPy can optionally use LLVM as a backend, to use similar tricks as Numba, and also has a JIT decorator.

In general, the Python language design creates limitations on the performance, see this video for example on How Python was Shaped by Leaky Internals. If you don't want to change language, you may hope that Python 4 ever comes around with syntax that can actually be optimized.

An optimized language - Julia as example

If you don't want to wait for Python 4, then there's Julia instead. Julia is a language that is optimized for talking to LLVM, while looking as similar as possible to high-level languages like Python and MATLAB. In short, it's an LLVM whisperer. This makes it fast and easy. See Why Was Julia created? to get an impression of the rationale.

Similar to Cython and Numba, you can optionally add type information to Julia, which can help the compiler, though Julia is good at type interference.

function f(x::Int, y::Int)
    return x + y
end

Is Julia better than Numba and Cython? For an opinionated and long blog post read Why Numba and Cython are not substitutes for Julia. There are also lengthy discussions in this discourse on Why weren't Numpy, Numba, SciPy good enough?. And I also like Martin Maas's blog post about Julia vs Numba and Cython.

I would summarize the benefits as: You don't have to decorate your code, Julia is the JIT decorator. You can write the same code for CPU and GPU. When compiling ahead-of-time, it's again the same code. The compiler can optimize across all Julia code, not just a single package like NumPy that you are currently using. Composability is often praised: Julia packages work easily together.

The downside of Julia, if you are coming from another language like Python, is obviously that you have to learn another language. Though I wonder how much more difficult Julia is compared to learning a complex framework like Numba. And writing optimized Cython can be considered similar to writing another language. Julia has many similarities with Python, check for example this cheat sheet to compare MATLAB to Python to Julia, except Julia bypasses the problems that make Python difficult to compile to LLVM.

Similar to Cython and Numba, naive Julia code is good, but not necessarily as performant as optimized C. Read the performance tips to get the most out of your code.

If you want to move gradually to Julia, you can embed Julia into Python via PyJulia or the more recent two-way package PythonCall. Or you can re-use existing Python code inside Julia, via PyCall.jl, or the aforementioned PythonCall. This way you can use Julia as if it's yet another Python framework, instead of a completely new language. You can even send Pandas dataframes to Julia and back.

Similar to Numba, Julia stores the compiled native code in memory. In the upcoming Julia 1.9 release, this code will also be automatically cached on disk per Julia package, so compilation happens only once, instead of the first call in every new Julia session. Ahead-of-time compilation was always possible with PackageCompiler.jl. The PackageCompiler is not actually compiling (Julia and LLVM do that), it simply gathers the in-memory compiled code and stores it in a .so library (or .dll on windows). This can be used for a standalone library of your compiled code, and will automatically include all dependent libraries. I have written a long tutorial on how to embed such Julia libraries inside C++ on my private website.

Static compilation of Julia, into tiny libraries fully independent of the runtime, is in an early stage with StaticCompiler.jl and StaticTools.jl, but needs more investment. Once you try out static compilation, you will notice that it enforces limitations on your code, because you cannot use all the fancy dynamic language features. I believe this is an unavoidable trade-off in any of the discussed technologies so far, but I'd love to be surprised on this point.

Other attempts to make a fast language easier to use are Zig, Swift and GoLang to a certain extent. Rust is very interesting, but I would not call it easy for scientists. None of them are targeting numerical computing as much as Julia.

Mojo

A new language that was revealed very recently is Mojo. In their article Why Mojo? they rephrase the two language problem as a Two World Problem, or even three world problem (Python, C++ and CUDA) for machine learning. From the code snippets on their website it looks like they want a Python compatible language that has features of Rust. Note that Mojo is not yet released to the public, the mojo github repository is empty at this time of writing, so we don't even know if Mojo will have a permissive license. Ambitious, but very young. We'll keep an eye on this one.

Final Comparison

What are good comparison criteria? I have chosen a few below. The development community size can be used as an estimate of how much effort goes into each project. Other than that I have tried to compare the usage and technology choices. I don't want to give quantitative performance comparisons here, they are heavily dependent on your use case, but from the benchmarks on complex examples that I have seen, Julia typically performs best. However, performance might not be your main criterion. Other unlisted aspects, such as debugging, profiling or cloud deployment may be more relevant for your use case.

	Cython	Numba	Julia
Contributors	430	298	1386
Github stars	7.9k	8.6k	42.2k
Backend tech	C transpiler	LLVM	LLVM
Usage	Decorators and Cython types	Decorators everywhere	Learn another language
Python interoperability	Import cythonized modules	Just-in-Time (JIT) decorators	PythonCall package
Performance	Decent	Good	Best
Distribution	Ship the .so with all dependencies	Ahead-of-Time compilation decorators	Ahead-of-Time compilation via PackageCompiler.jl

While I have a preference for Julia, I tried to stay as unbiased as possible in this blog post. All options are amazing open source projects and are maintained by mostly voluntary developer communities. Investigating them all is a humbling experience in the complexity of software technology.

There are other software engineering requirements that I have not yet included in this post, but might be important for your use case:

Package management. How easy is it to create and install a package, with all of its dependencies, in the chosen technology. I believe Julia's Pkg is superior here.
How does dependency management and distribution of binary artifacts work exactly? These nitty-gritty details can slow down your project. I have not yet tried this extensively for Numba. Cython involves some manual work. Julia has an artifact manager that works together with the package manager inside PackageCompiler.
Complex cases, like multi-threading inside the framework or external threads calling from another language. (Note: Julia doesn't have the global interpreter lock (GIL) like Python. Cython can release the GIL in certain cases.) I haven't gone into such topics yet, but there are many complex use cases that you may want to gradually add to your codebase. How far can you go with each technology before hitting a wall?

Finally, remember that a technical solution does not necessarily result in a cultural improvement. If you hired a lot of scientists or analysts or domain experts, and none of them have the necessary software skillset, it is difficult to improve collaboration with software engineers by forcing a 'better' technology onto them. You will have to empower your scientists to learn the necessary software development tools and processes, such as version control, test-driven development and continuous integration. Vice versa, your software engineers can learn the business domain and the tricks of numerical computing with the help of your scientists, to know exactly what code to write. By bridging these gaps, you can create a more effective team that can leverage the full potential of the technology investments you make.

Thanks to Jorge Vieyra and Jeroen van der Meer for reviewing and suggesting excellent improvements to the article.

These long posts take me quite some time and effort to write. If you like them, please encourage me by leaving a comment with suggestions or subscribe to my email newsletter. With enough support, I intend to write a book about building and deploying professional numerical computing applications. With Julia examples.

Automate your Code Quality in Julia

Matthijs Cox — Fri, 05 May 2023 14:21:25 +0000

Code quality is a topic in Julia that I believe deserves more attention from both users and developers. The Julia language originated in academia and focused heavily on helping scientists write better code, which is going great and deserves much praise! However, to onboard more software engineers and professional organizations we're going to have to invest even further into code quality and automated code quality tools and other methods such as used in the field of "quality assurance". In this article I'll explore the current state in the Julia ecosystem.

At our workplace we have investigated the following tools and practices. I'll start from generic practices and then move on to more advanced tools.

Package structure
Unit testing with Pkg.jl
Automated testing and Continuous Integration (CI)
Code Coverage with Pkg.jl
Documentation testing with Documenter.jl
Style guides and JuliaFormatter.jl
Static Code Analysis with StaticLint.jl
Quality Assurance with Aqua.jl
Type stability with JET.jl

Let's have a look at all of them.

Packages

The Julia community uses a standardized package structure and has a single package manager called Pkg.jl. There's also a single documentation system and a single testing system. This consensus alone helps tremendously with automating any workflows in your projects and organizations.

Please make sure you share professional code with others via packages. It's straightforward to adhere to the package structure. To get you started with creating well-defined (open source) packages, you can look at PkgTemplates.jl.

When setting up my first open source Julia package, I enjoyed the documentation of the JuliaTemplateRepo.jl, which goes through all the basic steps and configurations for a Julia package.

Unit Testing

Nowadays unit testing is a common practice in professional software engineering. Developers in Julia should be no exception. Fortunately, according to Viral Instruction, 89% of all open source Julia packages have tests, including a lot of beginner packages. It's safe to say that the Julia community puts a lot of emphasis on testing, which I think is remarkable for a language that originated in academia. This really sets a good example.

All unit testing use the Test.jl package, which is shipped with the base language. There are some extensions like TestSetExtensions.jl and ReTest.jl, but I believe you can do most of your work with Test.jl.

Getting started with testing is trivial in Julia, just add a test/runtests.jl file to your package and add code like this:

using Test
using MyPackage

@testset "MyPackage tests" begin
    @test 1==2
end

You can then run the tests with Pkg.test("MyPackage") which starts an isolated sandbox environment for the tests.

If you are creating a reproducible package for others, including your future self, then there is absolutely no reason to not write tests. However, writing good descriptive tests that cover all your bases is a more advanced art. Many books are written on this topic. I would advise to just get started and improve your testing strategy as you grow your codebase.

Automated Testing & Continuous Integration

Once you have unit tests defined, this aspect is a low hanging fruit for automation. It's very easy to forget to run the unit tests before committing. Automatically testing the code will save you from simple mistakes.

You can setup Github Actions or use a tool like Jenkins, to automatically run the tests upon every commit and block developers from merging if the unit tests do not pass. PkgTemplates.jl will typically already generate this Github Action for your open source package.

The tools and practice of frequently and automatically checking your code development is called Continuous Integration (CI). Inside software organization this is often combined with Continuous Deployment (CI/CD). All continuous integration (CI) tools that the Julia community uses can be found in JuliaCI · GitHub, I will address a few of those. If you want to configure your own Github actions you be inspired by the examples in Julia Actions · GitHub.

Code Coverage

A straightforward metric for monitoring your code quality is to check the fraction of code covered by your tests. Similar to automating your tests, this is a very low hanging fruit in the Julia community.

The code coverage generation itself is embedded inside the package manager via Pkg.test("MyPackage", coverage=true) . This will generate .jl.cov files with information about how often a line of code is touched by the tests.

Analyzing the code coverage visually line-by-line, for example inside VS Code, can help you identify where you are lacking tests, or help you find out where you can delete un-used code that is never called by your functions. You can automatically send the files to a service, like Coveralls.io or Codecov.io, and analyze in the browser. Here's an example in Flux's functions.jl file, that's has one uncovered line (note that it's very common to forget to test our errors or other corner cases):

You can calculate statistics on these coverage files, for example with Coverage.jl or with one of the services above. That way you can monitor coverage statistics over time. And if you use such a service for your open source package you can add a shiny badge to your readme to show-off your coverage.

Other commercial tools are busy adopting Julia's code coverage, check with your preferred supplier if they already support Julia. If not, please request them to do so.

Documentation Testing

Good documentation is incredible important for the users of your package, both in the open source community as well as inside organizations. Unfortunately, documentation that includes code examples can run out of sync with your code if you forget to test those. But it's very easy to automatically test those code examples with Documenter.jl doctesting.

Whenever you write a docstring or write code snippets in your docs folder, just add jldoctest and the expected output. (I replaced the triple ticks with double ticks, don't know how to properly format that here on Forem.)

``jldoctest
a = 1
b = 2
a + b

# output

3
``

Now just add Documenter.doctest(MyPackage) to your automated tests, and you know immediately when your examples no longer work. Easy!

Style Guides & Code Formatting

One of the challenges when working with many people on a single codebase is to adhere to a consistent coding style that is pleasant and unambiguous for everyone. This is where style guides help, together with formatting tools that make it easy to adhere to such a style guide.

The primary open source Julia style guides are at the moment are:

YAS (yet another style)
Blue style
SciML style

You can automatically format your Julia files in VS Code with the click of a button. If you are a command line user, or want to automate the formatting in your CI, you can use the underlying JuliaFormatter.jl package.

This should get you started with code formatting in no time. Discuss with your colleagues which style guide you prefer. Personally I use the BlueStyle so far, but the SciML style is relatively new, so looking into that one as well.

Static Code Analysis

Static code analysis will look at the code without executing it. One package we found is StaticLint.jl, which is used primarily by the Julia VS Code plugin LanguageServer.jl to report on potential problems in your code. These are normally reported under the "Problems" tab. Here I found a few potential problems inside the DataFrames.jl package:

The VS Code plugin also reports on potential problems when hovering over code, such as a warning about this unused function argument. These are not reported in the "Problems" tab.

StaticLint still misses some documentation for users, but you can use the following script and read my discussion in an issue here. After some fiddling with the code and environments I am able to obtain the same "diagnostics" on my REPL for a given file:

julia> docs[3]
Document: file:///c%3A/Users/matcox/Documents/Julia/static_lint/src/abstractdataframe/selection.jl

julia> docs[3].diagnostics[1]
LanguageServer.Diagnostic(LanguageServer.Range(LanguageServer.Position(223, 15), LanguageServer.Position(223, 36)), 4, "UnusedFunctionArgument", missing, "Julia", "An argument is included in a function signature but 
not used within its body.", [1], missing)

So StaticLint.jl can be used, but it's not yet user friendly for integration into any command line interfaces or automated tooling.

Quality Assurance with Aqua.jl

The package Aqua.jl is developed to automate quality assurance for Julia. The readme is clear on what it checks:

There are no method ambiguities.
There are no undefined exports.
There are no unbound type parameters.
There are no stale dependencies listed in Project.toml.
Check that test target of the root project Project.toml and test project (test/Project.toml) are consistent.
Check that all external packages listed in deps have corresponding compat entry.
Project.toml formatting is compatible with Pkg.jl output.
There are no "obvious" type piracies

Aqua provides a function to use in your testing environment, which will report issues and throws an error so your tests cannot pass unless all Aqua checks pass. Here's a snippet of what we get for DataFrames:

julia> using Aqua, DataFrames

julia> Aqua.test_all(DataFrames)
17 ambiguities found
Ambiguity #1
<=(a::Integer, b::SentinelArrays.ChainedVectorIndex) in SentinelArrays at C:\Users\matcox\.julia\packages\SentinelArrays\BcfVF\src\chainedvector.jl:208
<=(x::BigInt, i::Integer) in Base.GMP at gmp.jl:696

Possible fix, define
  <=(::BigInt, ::SentinelArrays.ChainedVectorIndex)

... 

Test Summary:    | Fail  Total  Time
Method ambiguity |    1      1  8.3s
ERROR: Some tests did not pass: 0 passed, 1 failed, 0 errored, 0 broken.

Unfortunately I only find ambiguities for DataFrames, maybe I should find a package with more problems. You can also run all the underlying checks independently if you read the Aqua documentation. Note that if you only want to check for ambiguities, you can also choose to call Test.detect_ambiguities directly from the standard Julia Test package.

A nice addition to Aqua would be a way to store the found issues in a standardized file format instead of printing them on the REPL. Similar to code coverage reporting, this can help to make overviews in automated systems. Now we would have to capture the printed output and parse that somehow.

Finding Type Instability with JET.jl

There is non-stop activity in the Julia community to analyze our own code for improvements. An advanced package is JET.jl, which uses the Julia compiler itself to detect potential bugs and type instabilities.

What is type instability? This happens when the type of a variable changes unexpectedly. Here's a simple example that outputs either an integer or a floating point variable:

function foo(value)
    if value < 1
        return 1
    else
        return 1.0
    end
end

Type instability is bad for performance because the compiler cannot infer the types and generate optimal native code. It may also point at bugs in your code, if you did not intend to have such instability. Julia does not enforce type stability like certain languages, because it wants to remain an easy language to use. Sometimes you don't care about performance and don't want to worry about types, in which cases it's easy if you are not forced to worry.

If you just want to check whether the output value can be inferred, you can use Test.@inferred in your tests:

julia> using Test

julia> Test.@inferred foo(0.5)
ERROR: return type Int64 does not match inferred return type Union{Float64, Int64}

julia> Test.@inferred foo(1.5)
ERROR: return type Float64 does not match inferred return type Union{Float64, Int64}

However, when you want more certainty about the internals of your code then you can turn to JET. Most of JET is doing specific method analysis with @report_opt and @report_call . JET can do some limited static analysis of your entire package with the report_package function. Unlike @report_call , this means JET doesn't know what types you want to input into your methods, so it has to make some assumptions.

I do warn that the output of JET can be rather intimidating. Here's what you get when executing the example @report_call sum("julia") :

And that's just the example from the Quick Start page of JET.

We're still investigating how to use JET, because it is pretty advanced tooling. If you just started with types and Julia, I wouldn't dive right into this. Take your time to think about what type inference really means, and read the documentation of JET if you want to know more.

Reducing Compile Time with SnoopCompile.jl

Optimizing your code such that the compilation times are reduced is maybe not the first thing that comes to mind when thinking about "code quality", but it can improve the user experience of your package. Nobody likes to wait long to import your package. The SnoopCompile.jl package helps you with analyzing your code for such improvements. It "snoops" on the compiler and reports on it's findings.

There is a lengthy blog post from the SciML community on how they improved their compilation times with SnoopCompile and other tools, called How Julia ODE Solve Compile Time Was Reduced From 30 Seconds to 0.1. Definitely read that one for more information, I will not go into details here, but I do think SnoopCompile is worth a mention.

Conclusion

You have plenty of options to check the code quality of your Julia packages and improve the quality over time. If this feels like a daunting task as a beginning (Julia) developer, don't worry, you can slowly add these tools to your workflow over time. The most important thing is to start with a good package structure and basic testing. The fact that the Julia ecosystem is so focused on making quality easy for beginners is truly praise-worthy and will help us all in the long run.

For senior developers and managers looking into these tools, one thing to remember is that lot's of code quality tooling in Julia is written with the human developer in mind. This currently limits some of the integration in automated CI tools. I believe this topic deserves some more attention in the Julia community and more support from commercial code quality tooling vendors. The good thing is that due to the standardization of the Julia package management it is very easy to get started with a uniform automation system in your organization. As the tools improve for these systems, it will be easy to incrementally add such tools to any open source or internal CI workflows.

Thanks to my colleague Matthijs den Otter for helping with the investigation. If we find better ways to monitor your Julia code quality, I intend to share that on my blog, so don't forget to subscribe.

The Art of Multiple Dispatch

Matthijs Cox — Thu, 20 Apr 2023 13:50:45 +0000

Read the original article and subscribe to my blog on The Scientific Coder.

I love thinking visually by drawing doodles and schematics for my work. It's one of my favorite things to do, next to coding. When working with the Julia language, one visualization I enjoy is seeing the type space of a method that you are dispatching on. Normally I do this in my mind's eye, but let me clarify this by drawing some actual figures.

To start with the basics; Julia has functions and methods. A function is simply the name, like push! or read. Methods are specific definitions of a function, for certain types of arguments. Take for example push!(s::Set, x) or read(io::IO). From an object-oriented perspective you could say that methods are instances of functions.

For any given method you can consider the dispatching as slicing a part of the entire possible type space of that given function. For a given set of arguments of course. If you increase the number of arguments in the function definition, then more dimensions get added to the type space. I don't even know how to find the best written words for this, the visualization above just feels intuitive to me.

Let's take the function f and imagine for a moment that there are only 3 types in the whole Julia type universe: the Float64, Int64and the String. The Float64and the Int64are a subtype of Number, which is obvious I hope. By default in Julia if you specify no type in your function argument, then it will be assumed you mean the Anytype, of which every other type is a subtype.

A method f(::Any, ::Any) thus describes the entire space of all possible types for the function named f. On the other hand, a method like f(::Int64, ::String) is super concrete, it's a singular point in the type space.

You can use abstract types like Number or unions like Union{Float64, Int64} to capture a subset of the discrete type space. This way you can choose which part you want to define for your function, with the chosen set of types you will be dispatching on at runtime. Abstract types in Julia exist only for this dispatching purpose, to dispatch on a set of subtypes, they have no other influence on their subtypes what so ever. They are not dictators like classes in other languages.

I like these visuals. Some junior engineers wonder what "diagonal dispatch" is. I don't have any other way of explaining the concept then by just drawing it. The figure is immediately obvious. Diagonal dispatch happens when the type of all arguments is forced to be equal with f(::T, ::T) where T. This truly represent a diagonal through the type space. You can see it in the example above. You can also limit the diagonal dispatch to a subset with f(::T, ::T) where T<:Number and in higher dimensions you can be fancy like f(::T, ::T, ::S) where {T<:Number, S<:AbstractString} by adding multiple of these parametric types.

Note that when you define a method twice, you have to take care that it is clear which method gets dispatched on. The compiler will prioritize the one that is most concrete, so the one that is most specific about the types. In the figure above, I ordered them from most specific to least specific. You can try for yourself to see if I ordered them correctly.

For example if you define the following:

f(::Any, ::Any) = println("any & any")
f(::Int64, ::Int64) = println("int & int")

Then most function calls will run the broadest method because that one is defined for the whole type space, but when you input two integers you will call the very specific method f(::Int64, ::Int64). Let's give it a go:

julia> f("string", 5)
any & any

julia> f(4, 5)
int & int

From a visual perspective, we have created an overlapping dispatch, where one method is specifically defined for the integer case f(::Int64, ::Int64) and will be called when only integers are used as arguments.

There are some caveats here. If you are not careful the methods can become ambiguous and Julia won't like that. For example if you define the following:

f(::Any, ::String) = println("any & string")
f(::String, ::Any) = println("string & any")

Which one of these methods should be called with f("string", "string")?

julia> f("string", 5)
string & any

julia> f(5, "string")
any & string

julia> f("string", "string")
ERROR: MethodError: f(::String, ::String) is ambiguous.  Candidates:
  f(::Any, ::String) in Main at REPL[8]:1
  f(::String, ::Any) in Main at REPL[9]:1
Possible fix, define
  f(::String, ::String)

Yikes, that's impossible! Fortunately there is a fix proposed, by explicitly defining the ambiguous case. Though perhaps you should reconsider what your actual intentions are in this design. The visual representation below hopefully makes the mistake more clear. At first there is confusion because the two dispatches overlap and neither is more specific than the other, but we can fix it by defining a more concrete method in the conflicting area.

When you define a lot of methods, you are creating a colorful patchwork in the type space of your function. You can come up with the craziest designs in your methods, but be careful. Finding the right balance of a few big broad abstract methods versus multiple tiny concrete methods is a true art in Julia.

People do not often share how they visualize the code design in their mind, while I believe this really shapes the creative process. The closest visual representation in Julia I have seen is the article about Julia's dispatch with Pokemon types. You can read that for more detailed examples with Julia's multiple dispatch.

That concludes this short artsy post, but I hope it helps the visual thinkers in the programming community! Let me know if you use different kinds of visualizations in your coding work. And don't forget to subscribe to my newsletter.

Embedding Julia libraries in C++

Matthijs Cox — Thu, 21 Jul 2022 09:57:03 +0000

Embedding compiled Julia libraries inside a foreign environment with a C-callable interface is an advanced topic on the border of my expertise. It's somewhat underdocumented and non-trivial, so I've made this tutorial by writing down the steps I followed myself.

The fundamentals are explained in the Embedding section in the Julia manual. For the datatypes that can be passed between C and Julia, see calling-c-and-fortran-code.

I will make things hard for myself and compile on Windows. Also note that I am using c++ instead of c, but most steps will be similar for c.

High level the steps involved are:

Setup a c/c++ compiler
Make a Julia package
Write the Julia c-interface functions
Write the c++ code calling those functions
Compile Julia code to a library with PackageCompiler.jl
Compile c++ and link it to the Julia library

I'll explain how I did these steps and especially touch upon the Julia-C interface functions and types.

Compiling C++

I will compile C++ on Windows. It's notoriously difficult on Windows to find the right compiler. I got burned once on some generic MinGW compiler, creating all kinds of wrong string conversions, took us days to find out. So far the MinGW x86_64-8.1.0-posix-seh-rt_v6-rev0 version is working fine on my personal system. You can also use the Microsoft Visual C++ compiler tools, you can download the command line tools separate from the Visual Studio IDE. Make sure the right tool is added to your windows path. Use where g++ or where gcc to find out which one you are using.

Let's start at the real basics. So make a file called example.cpp.

#include <iostream>

int main()
{
    std::cout << "hello world" << std::endl;
    return 0;
}

Now on your command line run: g++ example.cpp -o example. This will compile an example.exe, which you can then run on your command line with example.exe and see the output string "hello world".

Use gcc for c and g++ for c++. These functions come with a bunch of compile and link options. There are so many options... The most common ones are:

-c (Compilation option). Compile only. Produces .o files from source files without doing any linking.
-I[/path/to/headers] (Compilation option). Include a folder with header files, like julia.h.
-o file-name (Link option, usually). Use file-name as the name of the file produced by g++ (usually, this is an executable file).
-l[library-name] (Link option). Link in the specified library, for example -ljulia for libjulia.dll.
-L[/path/to/shared-libraries] -l[library-name] (Link option). Link in the specified library from a given folder.

Create a Julia package

I assume you have basic Julia knowledge, including package development. Julia libraries are created from a package, so go ahead and make one in a folder. Simply run:

julia> cd("my/desired/package/path")
julia> import Pkg; Pkg.generate("MyLibPackage")

Compiling Julia

PackageCompiler.jl is your primary friend and documentation is improving rapidly. Perhaps in the future we will have an incredible StaticCompiler.jl, but that's for another time.

PackageCompiler does take a few minutes to compile any Julia code, even a simple hello world. That's because all of Julia base is included, regardless of whether actually need all base functionality or not.

We will be using the create_library functionality of PackageCompiler. The easiest way is to find some example build scripts from others and add those to your package, so let's go find one.

Existing examples?

A skeleton Julia compilation example is already available here.

For a pure C example implementation see Simon Byrne's libcg. See if you can already run that example.

Steps:

Clone the repository in a folder. You know: git clone https://github.com/simonbyrne/libcg.git.
Run the Makefile. Uh oh...

OK, running the Makefile on Windows isn't trivial either. StackOverflow provides some answers. My c/c++ mingw installation comes with mingw32-make, but that doesn't work with this specific Makefile, see this issue. Advise is to install Cygwin together with make, because the examples repositories use a lot of shell scripting which doesn't work on Windows.

OK, so this example is not so simple on Windows. In the end I decided to write my own Windows Makefile for my own c++ code and run it with mingw32-make.

You can find my personal examples in this repository called embedjuliainc.

Interlude: Makefiles?

Compiling c/c++ projects generally involves a lot of steps to compile multiple files and link libraries together. This can become a complex build process and you don't want to type all commands manually all the time. So people made build tools! Yay! Typically these come with their own kind of scripting language that you need to learn. Hmm, OK, just do it. It's also generally Unix oriented, not Windows. Make is a one common build tool, here's a Makefile tutorial.

As an example, I enjoyed this g++ makefile example. It explains all compilation steps and how to make a Makefile for our simple example c++ program above. Very informative, please read and try it out!

Interlude: Name Mangling

Another thing that is different in all my examples below, is that c++ mangles the names of functions. That means that a function f(int) get's turning into something like __f_i(int). To avoid this we need to use extern C whenever we define function interfaces. This took me a while to figure out, so a lesson learned the hard way!

extern "C"
{
    #include "julia_init.h"
}

Basic Types

I created an example of basic type data transfer between Julia and c++ in this repository. It contains implementations of the following data types:

Booleans
Integers
Doubles
Strings
Structs
Arrays
Enumerations

For a full comparison look at:

The Julia source code file BasicTypes.jl
The C++ file calling into the Julia library in basic.cpp which uses basic.h.

In general Julia's C interface documentation is your friend.

The typical pattern is straightforward, this Julia code:

Base.@ccallable function test_int64(myInt64::Int64)::Int64
    return myInt64
end

You can call with such c++ code:

int64_t test_int64(int64_t myInt64);
int64_t myInt64 = 9006271;
test_int64(myInt64);

Header file

The PackageCompiler.create_library function also likes to receive a header file. So please write all your c/c++ interface types and functions inside this header file. I wonder if this header file could be created automatically from all the Base.@ccallable function signatures? Now you have to manually keep the Julia source code and the header file in sync.

Struct types

As the manual says in the section "Struct Type Correspondences" you can pass structs. Fixed size arrays in c/c++ map onto the NTuple in Julia.

Warning! Be absolutely certain that the Julia struct definition matches the c struct definition!

Nested structs work just fine:

struct ChildStruct
    ChildStructId::Cint
end

struct ParentStruct
    ParentStructId::Cint
    myChildStruct::ChildStruct
    staticArray::NTuple{3, Cint}
end

Base.@ccallable function test_nested_structs(myParentStruct::ParentStruct)::ParentStruct
    return myParentStruct
end

The corresponding C/C++ interface:

struct ChildStruct
{
    int ChildStructId;
};

struct ParentStruct
{
    int ParentStructId;
    ChildStruct myChildStruct;
    int staticArray[3];
};

ParentStruct test_nested_structs(ParentStruct myParentStruct);

Passing by reference

If you want to avoid any copying and additional memory allocation, you will have to pass the data by reference as a pointer. A typical example is to pass an array by reference. In the example code I pass an Array{Cint}. Note that you also need to pass the dimensions of the array, in this case only the length, since we assume it's a vector.

Base.@ccallable function test_array(myArrayPtr::Ptr{Cint}, myArraySize::Cint)::Cvoid
    myArray = unsafe_wrap(Array{Cint}, myArrayPtr, myArraySize, own=false)
    # do stuff, mutating an element will mutate the original C memory, be careful
end

Mutating

Let's pass a struct by reference and mutate it. You can first load the entire struct with unsafe_load.

Base.@ccallable function test_nested_structs_ptr(myParentStructPtr::Ptr{ParentStruct})::Cvoid
    myParentStruct = unsafe_load(myParentStructPtr)
    myParentStruct.myChildStruct = ChildStruct(12)
    return Cvoid()
end

For arrays you can first unsafe_wrap like in the section above. Or you can immediately unsafe_store! on an individual element. As always be very careful with these operations.

Nested variable sized structs

Be careful with placing variable-sized arrays inside structs (this includes strings). You will have to somehow pass along the array size and manually unwrap such complexity. I still have to write a complex example for this. The Julia manual has a very minimal example in the "Calling C and Fortran" section.

For example, I would like to pass such a structure by reference:

struct VarSizedStruct
    len::Cint
    varArray::Array{Cdouble}
    string::Cstring
end
struct NestedVarStruct
    len::Cint
    varArray::Array{VarSizedStruct}
end

I don't know if it is a smart thing to do. You'll have to manually interpret the bytes...

I think the best way to go, is to place pointers inside the structs and to manually unsafe_wrap every array.

struct VarSizedStruct
    lenArray::Cint
    varArrayPtr::Ptr{Cdouble}
    lenString::Cint # in case the string is not NUL-terminated
    string::Ptr{Uint8}
end
struct NestedVarStruct
    len::Cint
    varArrayPtr::Ptr{VarSizedStruct}
end

function Base.@ccallable(nestedPtr::Ptr{NestedVarStruct})::Cvoid
    nested = unsafe_load(nestedPtr)
    len = nested.len
    struct_array = unsafe_wrap(Array{VarSizedStruct}, nested.varArrayPtr, len, own=false)
    last_struct = struct_array[len]
    last_double_array = unsafe_wrap(Array{Cdouble}, last_struct.varArrayPtr, last_struct.lenArray, own=false)
    last_string = unsafe_string(last_struct.string, last_struct.lenString)
    return Cvoid()
end

Basically you are not passing along the struct with data, but a collection of pointers and lengths. You'll have to then manually convert the data to an internal Julia representation of your choosing. It seems error prone and feels like it could be automated.

Garbage Collector

The manual is clear on memory management from within c/c++. You can even disable the garbage collector if you want.

One thing we ran into while testing the c-callable Julia functions from within a Julia script, is that the garbage collector may remove your object even while the function is executing. This can happen when passing pointers instead of objects and leads to horribly unexpected segmentation faults. Please use the GC.@preserve for those cases.

I placed an example in the precompile statements file:

arr = Cint[1,2,3]
arr_pointer = Ptr{Cint}(pointer_from_objref(arr))
len_arr = Base.cconvert(Cint, length(arr))
# please garbage collector, preserve my array during execution
GC.@preserve arr test_array(arr_pointer, len_arr)

Error handling

A typical old C way of error handling is to always return an integer on the C interface. The Julia code is then responsible for catching errors and returning the corresponding error integer. Any other desired output arguments are passed as mutable input parameters on the C interface. If you want the error messages as well, you could also pass along a struct with the error code/type and a Cstring with the error message. I've added an example with an ExceptionHandler.jl package for this case to my repository.

Here's a simple example:

error_code(::Exception)::Cint = 2

Base.@ccallable function something(inputPtr::Ptr{Cint}, outputPtr::Ptr{Cint})::Status
    resultCode::Cint = 0
    try
        # do stuff
    catch e
        resultCode = error_code(e)
    end
    return resultCode
end

But that is not what I am looking for. I want a way to catch Julia exceptions inside c++. The Julia manual embedding section on exceptions is not clear on how to do this for native Base.@ccallable functions. Through experimentation I found out that exceptions cannot be caught by a regular try/catch block inside the c++ code wrapping around the julia library call.

Let's say we have a function that throws errors on the c-interface:

Base.@ccallable function throw_basic_error()::Cint
    throw(ErrorException("this is an error"))
    return 0
end

Missed exceptions

This would be an expected way to catch errors, but it doesn't work:

try
{
    throw_basic_error();
}
catch (const std::exception& e)
{
    std::cout << "\n a standard exception was caught, with message '"
                << e.what() << std::endl;
}
catch (...)
{
    std::cout << "\n unknown exception caught" << std::endl;
}

When running this code, you will see an error like fatal: error thrown and no exception handler available.. Our current assumption is that this is due to Julia initializing in another thread or process than the c++ code itself.

Catch those exceptions

After some digging we found the JL_TRY and JL_CATCH macros in the Julia header file. These can be used to catch Julia exceptions:

JL_TRY {
    throw_basic_error();
}
JL_CATCH {
    jl_value_t *errs = jl_stderr_obj();
    std::cout << "A Julia exception was caught" << std::endl;
    if (errs) {
        jl_value_t *showf = jl_get_function(jl_base_module, "showerror");
        if (showf != NULL) {
            jl_call2(showf, errs, jl_current_exception());
            jl_printf(jl_stderr_stream(), "\n");
        }
    }
    return 1;
}

To print the error type and message, you will have to use functions directly from the Julia runtime. We did not yet find a nice and easy way to convert the Julia exception into a c++ exception.

Outlook

The package CBinding.jl can automatically generate the Julia types from a c header file. There is a lot of knowledge in that package, so I should investigate it better. I also wonder if it's possible to do the opposite: to generate a c header from Julia types and functions.

Other interesting packages are Cxx.jl and CxxWrap.jl. These packages focus on embedding c++ (libraries) inside Julia, but again contain a lot of knowledge and some examples on embedding Julia inside c/c++.

Conclusion

Embedding in c/c++ is non-trivial. I advise to avoid embedding if you can ;) If you cannot avoid embedding, then use the examples from this post, but I do not guarantee that all my examples are safe and robust. I did not yet create examples with complicated nested variable sized structures and arrays, that's for another time. For now I'd advise to keep the c-interface as simple as possible. Good luck!

Acknowledgements

Special acknowledgements to Daan Sperber, Biao Xu and Evangelos Paradas for helping me figure out a lot of the steps involved.

Julia Data Modeling: A Pointy Example

Matthijs Cox — Wed, 20 Jul 2022 12:45:42 +0000

Julia's high performance gives the developer a lot of freedom in their choice of data structures. The downside of this freedom is that you need to spend more time thinking about the optimal data model for your code.

In a low performance language you are often restrained by the interface of the underlying high performance library. For example, in the MATLAB language or the Python Numpy package you are encouraged to make vectors of only numbers.

In Julia you can define your own custom types, which can be just as performant, or even more performant depending on your use case, than the already available types.

Lot's of Points

To illustrate the possibilities, I'll take a simple example that I personally encountered a lot. You've got some data points on an XY plane and the data is labeled as valid or not.

The straightforward way to model this data (coming from a vectorized framework) is to define these as vectors.

To keep them easily together in Julia, you'll then want to put them in a custom type:

struct VectorData
    x::Vector{Float64}
    y::Vector{Float64}
    validx::Vector{Bool}
    validy::Vector{Bool}
end

The reason that vectors (of primitive types) are often performant is due to how your computer stores memory in data. When you apply operations that need to iterate over x for example, the compact memory layout is quite ideal. You can read this excellent article about What scientists must know about hardware to write fast code.

So this can be a very good data structure. For example for multiplying your x values with a matrix.

If you spend a lot of time iterating over the data points however, like filtering the valid data based on some function, things become a bit more messy.

Also some developer might suddenly decide not all vectors in VectorData need to be equally long. Or the user of your data structure accidentally only changes one vector. If you want the vectors to be mutable, yet always of equal length, then you should check the vector lengths in all your functions to be safe :(

Let's try an alternative.

struct PointData
    x::Float64
    y::Float64
    validx::Bool
    validy::Bool
end
struct PointVector
   points::Vector{PointData}
end

Note that I put the vector of PointData into a struct, which just makes function dispatching more pleasant.

Let's write our desired filter function and try them both. Decide for yourself which solution you find most readable.

function within_threshold(x::T, y::T, threshold::T) where T<:Real
    return x^2 + y^2 < threshold^2
end

function filter_points(pts::VectorData, threshold=0.4)
    keep = pts.validx .& pts.validy .& within_threshold.(pts.x, pts.y, Ref(threshold))
    return VectorData(pts.x[keep], pts.y[keep], pts.validx[keep], pts.validy[keep])
end

function keep_point(p::PointData, threshold)
    return p.validx && p.validy && within_threshold(p.x, p.y, threshold)
end

function filter_points(points::AbstractArray{PointData}, threshold=0.4)
    return filter(p -> keep_point(p, threshold), points)
end

function filter_points(pts::PointVector, threshold=0.4)
    return PointVector(filter_points(pts.points))
end

We can quickly try both functions out by constructing some example data:

function Base.convert(::Type{PointVector}, v::VectorData)
    points = PointData.(v.x, v.y, v.validx, v.validy)
    return PointVector(points)
end

using BenchmarkTools
N = 10_000_000 # making the vectors do not fit in my computer's cache ;)
vd = VectorData(rand(N), rand(N), rand(Bool, N), rand(Bool, N))

@benchmark filter_points(vd)

BenchmarkTools.Trial: 175 samples with 1 evaluation.
 Range (min … max):  25.246 ms … 62.595 ms  ┊ GC (min … max): 0.00% … 27.22%
 Time  (median):     26.840 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   28.581 ms ±  4.728 ms  ┊ GC (mean ± σ):  1.55% ±  5.52%

 Memory estimate: 6.60 MiB, allocs estimate: 14.

pv = convert(PointVector, vd)
@benchmark filter_points(pv)

BenchmarkTools.Trial: 170 samples with 1 evaluation.
 Range (min … max):  21.897 ms … 85.568 ms  ┊ GC (min … max):  0.00% … 8.27%
 Time  (median):     28.346 ms              ┊ GC (median):     0.72%
 Time  (mean ± σ):   29.538 ms ±  7.879 ms  ┊ GC (mean ± σ):  10.11% ± 9.53%

 Memory estimate: 228.88 MiB, allocs estimate: 4.

In this case it turns the times are roughly equal. My PointVector solution is slightly faster, but does use way more memory.

Interlude: Filtering better

I am a little annoyed at the memory usage of the filter function. So I am exploring how to improve that here.

The Base.filter() function might be the memory consumer. If you write a list comprehension it actually use Iterators.filter, so these two solutions behave equally. Less memory usage, but slightly slower:

function filter_points(points::AbstractArray{PointData}, threshold=0.4)
    return [p for p in points if keep_point(p, threshold)]
end
function filter_points(points::AbstractArray{PointData}, threshold=0.4)
    return collect(Iterators.filter(p -> keep_point(p, threshold), points))
end

@benchmark filter_points(pv)

BenchmarkTools.Trial: 117 samples with 1 evaluation.
 Range (min … max):  36.909 ms … 101.812 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     39.738 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   42.855 ms ±   7.999 ms  ┊ GC (mean ± σ):  1.13% ± 3.73%

 Memory estimate: 14.26 MiB, allocs estimate: 14.

I can also define the functionality myself. You can see it's pretty much equal to the list comprehension.

function filter_points(points::T, threshold=0.4) where T<:AbstractArray{PointData}
    new_points = T()
    for p in points
        if keep_point(p, threshold)
            push!(new_points, p)
        end
    end
    return new_points
end

@benchmark filter_points(pv)

BenchmarkTools.Trial: 116 samples with 1 evaluation.
 Range (min … max):  41.746 ms … 57.218 ms  ┊ GC (min … max): 0.00% … 22.16%
 Time  (median):     42.415 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   43.246 ms ±  2.438 ms  ┊ GC (mean ± σ):  1.21% ±  4.03%

 Memory estimate: 14.26 MiB, allocs estimate: 14.

It seems I can reduce memory usage at the cost of some runtime performance. In the source code of Base.filter I see that it pre-allocates a new vector and shrinks it afterwards, while in this code we grow the new vector dynamically. I'm keeping this info here for future reference.

More point possibilities

Having to choose between an array of structs, or a struct of arrays, is a typical scenario while coding in Julia. The best choice really depends on what you intend to do most with your data. If performance is not the issue, then I prefer the array of structs for reasons of clarity (that means the PointVector in the above example).

But let's not stop there and explore some more options together.

Why not accept the fact that we have two points intertwined in our data model?

struct Point{T}
    x::T
    y::T
end
struct TwoPoints
    values::Point{Float64}
    valid::Point{Bool}
end
struct VectorTwoPoints
    data::Vector{TwoPoints}
end

Or pair up the value and valid field. Perhaps this is a common pattern in your codebase to flag any value as valid or not?

struct ValidValue{T}
    value::T
    valid::Bool
end
struct VectorPointValues
    data::Vector{Point{ValidValue}}
end

Or accept that we always have a value and a valid per direction. The interesting thing here is that we can extend the directions easily in the future if we want to go 3D for example. To label the direction, we could use an enum, but those are not popular in Julia. One reason is that you cannot easily dispatch on enums. So you might find something like this:

abstract type Direction end
struct X <: Direction end
struct Y <: Direction end

struct DirectionalValidValue
    direction::Direction
    value::Float64
    valid::Bool
end
struct VectorDirectionalValues
    data::Vector{DirectionalValidValue}
end

And now we could go full circle and just make a struct array again:

struct DirectionalVectorData
    direction::Vector{<:Direction}
    values::Vector{Float64}
    valid::Vector{Bool}
end

This last one might be best if we often find ourselves concatenating the x and y value from the initial structure. For example if we want to fit the x and y data together in a 2D model.

I'm not going to test the filtering performance of all of these. I expect roughly equal performance as the PointVector solution, with a little extra overhead if we use structs inside structs.

The ideal structures will depend on your own use case; what you do with the data in the functions you write. I'm still learning what is best for my own use cases. A structure of arrays seems to be the first choice for performance, assuming your data fits into memory as continuous vectors. DataFrames.jl stores each column as a vector for this purpose. But the desire to want both a struct array and an array of structs is a common pattern as well. So common that it has it's own package: StructArrays.jl.

StructArrays

Let's have a brief look at StructArrays. I'm re-using the data generated in the example from the beginning.

julia> using StructArrays

julia> sa = StructArray(pv.points)
10000000-element StructArray(::Vector{Float64}, ::Vector{Float64}, ::Vector{Bool}, ::Vector{Bool}) with eltype PointData:
 PointData(0.8342412444868733, 0.7156389385059589, false, false)
 ⋮
 PointData(0.21446860650049493, 0.23189128521923086, false, false)

# you can get the properties again
julia> sa.x
10000000-element Vector{Float64}:
 0.8342412444868733
 ⋮
 0.21446860650049493

# we can apply our filter function on it
julia> filter_points(sa)
314598-element StructArray(::Vector{Float64}, ::Vector{Float64}, ::Vector{Bool}, ::Vector{Bool}) with eltype PointData:
 PointData(0.09048345510701994, 0.320553989305556, true, true)
 ⋮
 PointData(0.11419461083293259, 0.34645035134668634, true, true)

julia> @benchmark filter_points(sa)
BenchmarkTools.Trial: 73 samples with 1 evaluation.
 Range (min … max):  63.260 ms … 92.607 ms  ┊ GC (min … max): 0.00% … 14.19%
 Time  (median):     66.733 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   68.442 ms ±  5.277 ms  ┊ GC (mean ± σ):  0.48% ±  2.31%

 Memory estimate: 14.26 MiB, allocs estimate: 14.

Unfortunately it's slower than the direct implementations, StructArrays seems to add a bit of overhead. But the convenience might be worth it for you, and maybe you do not spend that much time iterating over the rows/elements of the struct array.

Custom struct array implementation

If you want to build your own custom struct array, you'll have to write your own iterators and other functionality. It's more work, but is it worth it? Let's see how it goes.

Basically I will just extend the VectorData type and link it to the PointData type via iterate and getindex and such. I will copy the structures here for clarity again.

struct PointData
    x::Float64
    y::Float64
    validx::Bool
    validy::Bool
end
struct VectorData
    x::Vector{Float64}
    y::Vector{Float64}
    validx::Vector{Bool}
    validy::Vector{Bool}
end

# remember there is some danger in this assumption
Base.length(v::VectorData) = length(v.x)
Base.eachindex(v::VectorData) = eachindex(v.x)
Base.IteratorSize(::Type{VectorData}) = Base.HasLength()

# I stole this function from AbstractArrays.jl, after reading StructArrays.jl
# we could also choose to make VectorData <: AbstractArray
function Base.iterate(A::VectorData, state=(eachindex(A),))
    y = iterate(state...)
    y === nothing && return nothing
    A[y[1]], (state[1], Base.tail(y)...)
end

function Base.getindex(A::VectorData, idx::Integer)
    return PointData(A.x[idx], A.y[idx], A.validx[idx], A.validy[idx])
end

# defining my own filter function, not sure if there is a better way
function Base.filter(f::Function, A::VectorData)
    keep = Int64[]
    for (index, p) in enumerate(A)
        if f(p)
            push!(keep, index)
        end
    end
    return VectorData(A.x[keep], A.y[keep], A.validx[keep], A.validy[keep])
end

Now we can use iteration to filter the points:

function filter_points(pts::VectorData, threshold=0.4)
    return filter(p -> keep_point(p, threshold), pts)
end

Let's look at the filtering performance.

using BenchmarkTools
N = 10_000_000
vd = VectorData(rand(N), rand(N), rand(Bool, N), rand(Bool, N))

@benchmark filter_points(vd)

BenchmarkTools.Trial: 82 samples with 1 evaluation.
 Range (min … max):  54.039 ms … 83.368 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     60.138 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   61.772 ms ±  5.472 ms  ┊ GC (mean ± σ):  1.02% ± 3.27%

 Memory estimate: 10.61 MiB, allocs estimate: 21.

You see it's hard to do a much better job than StructArrays. And definitely hard to beat an actual struct of arrays.