Skip to content
/ hx Public

A dalliance with C++17 TMP for NMR processing

License

Notifications You must be signed in to change notification settings

geekysuavo/hx

Folders and files

NameName
Last commit message
Last commit date

Latest commit

 

History

40 Commits
 
 
 
 
 
 
 
 
 
 
 
 
 
 

Repository files navigation

hx

A C++17 header-only library supplying razor-thin data structures, in the spirit of 'zero-cost abstractions', for managing multidimensional arrays of multicomplex numbers.

Introduction

The goal is to be able to write simple code for operations on arrays of known configuration. For example:

#include "hx/core.hh"

int main () {
  // heap-allocate a 1024x64x64 array of 3-complex scalars.
  using x_type = hx::array<hx::scalar<3>, 1024, 64, 64>;
  auto xp = std::make_unique<x_type>();

  // declare an index for array traversal.
  x_type::index_type idx;

  // traverse the array.
  do {
    // do something, e.g. conjugate each element:
    (*xp)[idx] = ~(*xp)[idx];
  }
  while (idx++);
}

In the next example, we allocate a 1200x320x125 array of 3-complex numbers on the heap (3 GB!) and apply a 3D mixed-radix multicomplex FFT to it:

#include "hx/core.hh"

int main () {
  using x_type = hx::array<hx::scalar<3>, 1200, 320, 125>;
  auto xp = std::make_unique<x_type>();

  xp->foreach_dim([&xp] (auto dim) {
    auto f = hx::fft::forward<x_type::base_type,
                              x_type::shape<dim.value>,
                              dim.value + 1>{};

    xp->foreach_vector<dim.value>(f);
  });
}

Basic benchmarks show the code generated here to be 3-5x slower than FFTW, but considering the amount of work required to implement multicomplex FFTs by hand using FFTW, this reduction in speed can be accepted. :)

Compiler-optimized processing graphs

The broader goal of hx is to support the construction of processing graphs like this:

// allocate an array to process.
using X = hx::array<hx::scalar<2>, 3, 3>;
auto x = std::make_unique<X>();

// init the array contents (e.g. load from file)...

// build and execute the processing graph.
auto p = hx::proc::node(x)
            .zerofill().zerofill<1>()
            .fft().fft<1>()
            .real();

auto y = p(x);

An interpreted language like Python would then take a script like this:

import hx

@hx.jit
def process(x: hx.array) -> hx.array:
    return (
        x
        .zerofill(dim=0)
        .zerofill(dim=1)
        .fft(dim=0)
        .fft(dim=1)
        .real()
    )

def main():
    x = hx.load("path/to/x.dat")
    y = process(x)

if __name__ == "__main__":
    main()

and write/compile the C++ code with appropriate array sizes.

Licensing

The hx library is released under the MIT license. See the LICENSE.md file for the complete license terms.

Releases

No releases published

Packages

No packages published