From 965d34b61b768a2374108985a90321ec60f9577c Mon Sep 17 00:00:00 2001
From: Bas van der Tol <tol@astron.nl>
Date: Tue, 7 Jul 2020 13:49:18 +0200
Subject: [PATCH] Run oskar-comparison test in CI/CD
---
.gitignore | 1 +
.gitlab-ci.yml | 12 +-
cpp/oskar/CMakeLists.txt | 6 +
.../oskar_evaluate_spherical_wave_sum.cc | 8 +-
demo/CMakeLists.txt | 2 +
demo/comparison-oskar/CMakeLists.txt | 27 +
demo/comparison-oskar/README.md | 43 ++
.../generate_basefunction_plots.py | 118 ++++
.../generate_basefunction_plots.sh | 5 +
demo/comparison-oskar/generate_oskar_csv.py | 44 ++
demo/comparison-oskar/main.cpp | 38 ++
demo/comparison-oskar/npy.hpp | 516 ++++++++++++++++++
demo/comparison-oskar/read_oskar_beams.py | 85 +++
demo/comparison-oskar/run_oskar.py | 170 ++++++
demo/comparison-oskar/telescope.tm/layout.txt | 2 +
.../telescope.tm/position.txt | 1 +
.../telescope.tm/station000/layout.txt | 2 +
scripts/coeff_scripts/oskar_csv_to_hdf5.py | 93 ++++
18 files changed, 1166 insertions(+), 7 deletions(-)
create mode 100644 demo/comparison-oskar/CMakeLists.txt
create mode 100644 demo/comparison-oskar/README.md
create mode 100644 demo/comparison-oskar/generate_basefunction_plots.py
create mode 100755 demo/comparison-oskar/generate_basefunction_plots.sh
create mode 100755 demo/comparison-oskar/generate_oskar_csv.py
create mode 100644 demo/comparison-oskar/main.cpp
create mode 100644 demo/comparison-oskar/npy.hpp
create mode 100644 demo/comparison-oskar/read_oskar_beams.py
create mode 100644 demo/comparison-oskar/run_oskar.py
create mode 100644 demo/comparison-oskar/telescope.tm/layout.txt
create mode 100644 demo/comparison-oskar/telescope.tm/position.txt
create mode 100644 demo/comparison-oskar/telescope.tm/station000/layout.txt
create mode 100755 scripts/coeff_scripts/oskar_csv_to_hdf5.py
diff --git a/.gitignore b/.gitignore
index 8c766b7b..cd726b63 100644
--- a/.gitignore
+++ b/.gitignore
@@ -4,3 +4,4 @@
.vscode/
build/
test_data/
+__pycache__/
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index ea054eb4..b28bb90c 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -74,23 +74,29 @@ build-oskar:
image: $CI_REGISTRY_IMAGE:$CI_COMMIT_SHORT_SHA
only:
- schedules
+ - master
dependencies:
- build-everybeam
before_script:
- apt-get -y install python3-pip
- - pip3 install numpy
+ # Install python requirements for the OSKAR "integration" test
+ - pip3 install numpy==1.19.0 scipy h5py astropy tqdm matplotlib
- mkdir -p /opt/oskar/build
- cd /opt/oskar && git clone https://github.com/OxfordSKA/OSKAR.git
- cd /opt/oskar/build
script:
- # CPP install
+ # OSKAR cpp install
- cmake -DCMAKE_INSTALL_PREFIX=.. -DCMAKE_BUILD_TYPE=Debug -DCMAKE_EXE_LINKER_FLAGS="-coverage" ../OSKAR/
- make -j8
- make install
# Python install
- export OSKAR_INC_DIR=/opt/oskar/include && export OSKAR_LIB_DIR=/opt/oskar/lib
- cd ./../OSKAR/python && python3 setup.py install && cd $HOME
- - python3 -c "import oskar"
+ - export PATH=/opt/oskar/bin:$PATH
+ # Run OSKAR comparison, set some env variables for this session
+ - export NPIXELS=8 && export APPLY_TRANSPOSE=OFF && MAX_ORDER=3 && TOLERANCE=1e-12
+ - cd /opt/everybeam/build
+ - make VERBOSE=1 comparison-oskar
build-doc:
stage: build-doc
diff --git a/cpp/oskar/CMakeLists.txt b/cpp/oskar/CMakeLists.txt
index 4a166f99..3bfc5678 100644
--- a/cpp/oskar/CMakeLists.txt
+++ b/cpp/oskar/CMakeLists.txt
@@ -14,6 +14,12 @@ add_library(oskar SHARED
string(TOLOWER ${CMAKE_PROJECT_NAME} projectname )
set_target_properties(oskar PROPERTIES LIBRARY_OUTPUT_NAME "${projectname}-oskar")
+# Make sure that when other targets within this project link against the oskar target,
+# they can find the include files.
+target_include_directories(oskar PUBLIC
+ $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>
+)
+
#------------------------------------------------------------------------------
# Link against HDF5 and OpenMP
target_link_libraries(oskar ${HDF5_LIBRARIES} ${HDF5_CXX_LIBRARIES} ${OpenMP_CXX_FLAGS})
diff --git a/cpp/oskar/oskar_evaluate_spherical_wave_sum.cc b/cpp/oskar/oskar_evaluate_spherical_wave_sum.cc
index 26f40781..b7f508ae 100644
--- a/cpp/oskar/oskar_evaluate_spherical_wave_sum.cc
+++ b/cpp/oskar/oskar_evaluate_spherical_wave_sum.cc
@@ -62,7 +62,7 @@ void oskar_evaluate_spherical_wave_sum(int num_points, const FP* theta,
oskar_legendre2(l, abs_m, cos_t, sin_t, p, pds, dpms);
if (abs_m == 0) {
sin_p = (FP)0;
- cos_p = sqrt(f_);
+ cos_p = -sqrt(f_);
const FP4c alpha_ = alpha[ind0];
oskar_sph_wave(pds, dpms, sin_p, cos_p, 0, alpha_.a, alpha_.b, Xt,
Xp);
@@ -74,9 +74,9 @@ void oskar_evaluate_spherical_wave_sum(int num_points, const FP* theta,
d_fact = std::tgamma(d_ + 1);
s_fact = std::tgamma(s_ + 1);
const FP ff = f_ * d_fact / s_fact;
- const FP nf = sqrt(ff);
- const FP4c alpha_m = alpha[ind0 - abs_m];
- const FP4c alpha_p = alpha[ind0 + abs_m];
+ const FP nf = sqrt(ff) * (2 * (abs_m & 1) - 1);
+ const FP4c alpha_m = alpha[ind0 + abs_m];
+ const FP4c alpha_p = alpha[ind0 - abs_m];
p = -abs_m * phi_x_;
oskar_sincos(p, &sin_p, &cos_p);
sin_p *= nf;
diff --git a/demo/CMakeLists.txt b/demo/CMakeLists.txt
index afcb9d06..f5dd3d49 100644
--- a/demo/CMakeLists.txt
+++ b/demo/CMakeLists.txt
@@ -1,3 +1,5 @@
+add_subdirectory(comparison-oskar)
+
set(TEST_MEASUREMENTSET CACHE STRING "measurement set used for testing")
if (TEST_MEASUREMENTSET)
diff --git a/demo/comparison-oskar/CMakeLists.txt b/demo/comparison-oskar/CMakeLists.txt
new file mode 100644
index 00000000..3ac1b9f1
--- /dev/null
+++ b/demo/comparison-oskar/CMakeLists.txt
@@ -0,0 +1,27 @@
+#------------------------------------------------------------------------------
+# CMake file for compiling a comparison between OSKAR and EveryBeam
+add_executable(comparison-oskar-generate-beampattern main.cpp)
+
+target_link_libraries(comparison-oskar-generate-beampattern oskar)
+
+file(COPY "${CMAKE_CURRENT_SOURCE_DIR}/telescope.tm/layout.txt"
+ DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/telescope.tm")
+
+file(COPY "${CMAKE_CURRENT_SOURCE_DIR}/telescope.tm/position.txt"
+ DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/telescope.tm")
+
+file(COPY "${CMAKE_CURRENT_SOURCE_DIR}/telescope.tm/station000/layout.txt"
+ DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/station000/telescope.tm")
+
+#------------------------------------------------------------------------------
+# comparison-oskar knits together the cpp code and the python scripts
+add_custom_target(comparison-oskar
+ COMMAND ${CMAKE_COMMAND} -E env EXTRA_PATH="${CMAKE_CURRENT_BINARY_DIR}:${CMAKE_SOURCE_DIR}/scripts/coeff_scripts"
+ "${CMAKE_CURRENT_SOURCE_DIR}/generate_basefunction_plots.sh"
+ DEPENDS comparison-oskar-generate-beampattern
+)
+
+add_test(NAME comparison-oskar-test
+ CONFIGURATIONS integration
+ COMMAND make comparison-oskar)
+
diff --git a/demo/comparison-oskar/README.md b/demo/comparison-oskar/README.md
new file mode 100644
index 00000000..4efbb8ff
--- /dev/null
+++ b/demo/comparison-oskar/README.md
@@ -0,0 +1,43 @@
+# Introduction
+
+The scripts in this directory compare the implementations of the spherical wave model for the element response in OSKAR and EveryBeam. The comparison is made per basefunction. One basefunction has two complex components.
+
+# Usage
+Comparison plots can be generated by giving the following command in the build directory
+```
+make comparison
+```
+This will generate a series of image files named basefunction*.png in the \<builddir\>/demo/comparison-oskar directory. The number of basefunctions depends on $L$, the maximum order, and is given by $N = 2(L+1)^2-2$.
+
+The maximum order can be set by the environment variable MAX_ORDER, for example
+```
+MAX_ORDER=3 make comparison-oskar
+```
+computes the basefunctions up to order 3, generating 30 image files.
+
+The application of the 'transpose' can be activated by the APPLY_TRANSPOSE environment variable, for example
+```
+APPLY_TRANSPOSE=ON MAX_ORDER=3 make comparison-oskar
+```
+This option was needed to make the results of OSKAR and EveryBeam the same.
+With the fix of the 'transpose' issue the default value was changed from 'ON' to 'OFF.
+
+The environment variable TOLERANCE makes the demo a test that fails if the difference between OSKAR and EveryBeam results is larger its value. For example
+```
+TOLERANCE=1e-16 make comparison-oskar
+```
+fails, while
+```
+TOLERANCE=1e-12 make comparison-oskar
+```
+succeeds.
+
+This test is part of the 'integration' configuration of CTests. It will be run when the integration configuration is run with
+```
+ctest -C integration
+```
+
+
+# Notes
+The read_oskar_beam.py script tranforms the polarization vector from a X,Y coordinate system to a theta,phi coordinate system. It still needs to be verified that this transformation is the same as in OSKAR.
+
diff --git a/demo/comparison-oskar/generate_basefunction_plots.py b/demo/comparison-oskar/generate_basefunction_plots.py
new file mode 100644
index 00000000..b9d048b1
--- /dev/null
+++ b/demo/comparison-oskar/generate_basefunction_plots.py
@@ -0,0 +1,118 @@
+import os
+import sys
+import numpy as np
+import matplotlib.pyplot as plt
+from generate_oskar_csv import generate_oskar_csv
+import run_oskar
+from read_oskar_beams import read_oskar_beams
+import subprocess
+
+# Check and set parameters
+l_max = int(os.environ["MAX_ORDER"]) if "MAX_ORDER" in os.environ else 2
+apply_transpose = (
+ os.environ["APPLY_TRANSPOSE"].upper() in ("1", "ON", "TRUE")
+ if "APPLY_TRANSPOSE" in os.environ
+ else False
+)
+tolerance = float(os.environ["TOLERANCE"]) if "TOLERANCE" in os.environ else 0.0
+npixels = int(os.environ["NPIXELS"]) if "NPIXELS" in os.environ else 256
+
+plt.figure(figsize=(10, 6))
+for em_idx in range(2):
+ for basefunction_idx in range(l_max * (l_max + 2)):
+ plt.clf()
+
+ generate_oskar_csv(basefunction_idx, em_idx)
+ run_oskar.main(npixels)
+ A = read_oskar_beams()
+
+ plt.subplot(2, 4, 1)
+ plt.imshow(np.abs(A[:, :, 0, 0]).T, clim=(0, 0.25), origin="lower")
+ plt.colorbar()
+ plt.title("abs(Etheta)")
+ plt.ylabel("OSKAR")
+
+ plt.subplot(2, 4, 2)
+ plt.imshow(np.abs(A[:, :, 0, 1]).T, clim=(0, 0.25), origin="lower")
+ plt.colorbar()
+ plt.title("abs(Ephi)")
+
+ plt.subplot(2, 4, 3)
+ plt.imshow(
+ np.angle(A[:, :, 0, 0]).T,
+ clim=(-np.pi, np.pi),
+ cmap="twilight",
+ origin="lower",
+ )
+ plt.colorbar()
+ plt.title("angle(Etheta)")
+
+ plt.subplot(2, 4, 4)
+ plt.imshow(
+ np.angle(A[:, :, 0, 1]).T,
+ clim=(-np.pi, np.pi),
+ cmap="twilight",
+ origin="lower",
+ )
+ plt.colorbar()
+ plt.title("angle(Ephi)")
+
+ l = int(np.sqrt(basefunction_idx + 1))
+ m = basefunction_idx - l * l + 1 - l
+ s = em_idx
+
+ if not apply_transpose:
+ generate_oskar_csv(basefunction_idx, em_idx)
+ else:
+ # flip the sign of m
+ generate_oskar_csv(l * l - 1 + l - m, em_idx)
+
+ subprocess.check_call(["oskar_csv_to_hdf5.py", "telescope.tm", "oskar.h5"])
+ subprocess.check_call(["comparison-oskar-generate-beampattern", str(npixels)])
+
+ B = np.load("response.npy")
+
+ if apply_transpose:
+ B *= (-1) ** (m + 1)
+
+ plt.subplot(2, 4, 5)
+ plt.imshow(np.abs(B[:, :, 0, 0]).T, clim=(0, 0.25), origin="lower")
+ plt.colorbar()
+ plt.title("abs(Etheta)")
+ plt.ylabel("EveryBeam")
+
+ plt.subplot(2, 4, 6)
+ plt.imshow(np.abs(B[:, :, 0, 1]).T, clim=(0, 0.25), origin="lower")
+ plt.colorbar()
+ plt.title("abs(Ephi)")
+
+ plt.subplot(2, 4, 7)
+ plt.imshow(
+ np.angle(B[:, :, 0, 0]).T,
+ clim=(-np.pi, np.pi),
+ cmap="twilight",
+ origin="lower",
+ )
+ plt.colorbar()
+ plt.title("angle(Etheta)")
+
+ plt.subplot(2, 4, 8)
+ plt.imshow(
+ np.angle(B[:, :, 0, 1]).T,
+ clim=(-np.pi, np.pi),
+ cmap="twilight",
+ origin="lower",
+ )
+ plt.colorbar()
+ plt.title("angle(Ephi)")
+ plt.gcf().suptitle("l = {}, m = {}, s = {}".format(l, m, s))
+ plt.savefig("basefunction{}-{}".format(basefunction_idx, em_idx))
+
+ if tolerance:
+ difference = np.nanmax(np.abs(A - B))
+ if difference > tolerance:
+ sys.exit(
+ "Difference between OSKAR and EveryBeam spherical wave model is {}, which is larger than the tolerance {}".format(
+ difference, tolerance
+ )
+ )
diff --git a/demo/comparison-oskar/generate_basefunction_plots.sh b/demo/comparison-oskar/generate_basefunction_plots.sh
new file mode 100755
index 00000000..9be17d4a
--- /dev/null
+++ b/demo/comparison-oskar/generate_basefunction_plots.sh
@@ -0,0 +1,5 @@
+#!/bin/sh
+
+export PATH=$EXTRA_PATH:$PATH
+python3 -B `dirname "${0}"`/generate_basefunction_plots.py
+
diff --git a/demo/comparison-oskar/generate_oskar_csv.py b/demo/comparison-oskar/generate_oskar_csv.py
new file mode 100755
index 00000000..91d2d421
--- /dev/null
+++ b/demo/comparison-oskar/generate_oskar_csv.py
@@ -0,0 +1,44 @@
+#!/usr/bin/env python3
+
+import argparse
+import os
+import scipy.io as sio
+import os.path
+import numpy as np
+import h5py
+from tqdm import tqdm
+from glob import glob
+
+output_dir = 'telescope.tm'
+element_id = 0
+
+
+freqs = [50]
+
+def generate_oskar_csv(basefunction_idx, em_idx=0):
+
+ l_max = int(np.sqrt(basefunction_idx+1))
+
+ # Parse all freqs
+ for freq in freqs:
+
+ A = np.zeros( (l_max,2*l_max+1, 2,2), dtype=np.complex128)
+ data = np.zeros((l_max * (l_max+2), 4), dtype=np.complex128)
+ data[basefunction_idx,em_idx] = 1.0
+
+ i = 0
+ for l in range(l_max):
+ n = (l+1) * 2 + 1
+ A[l,:n,:] = data[i:i+n,:].reshape(n,2,2)
+ i += n
+
+ for i, pol in enumerate(('x', 'y')):
+ for j, em in enumerate(('e', 'm')):
+ for reim, s in zip(('re', 'im'), (1.0, 1.0j)) :
+ filename = 'element_pattern_spherical_wave_{}_t{}_{}_{}_{}.txt'.format(pol, em, reim, element_id, freq)
+ alpha = np.real(A[:,:,i,j]/s)
+
+ np.savetxt(os.path.join(output_dir,filename), alpha, delimiter=', ')
+
+if __name__ == "__main__":
+ generate_oskar_csv(0)
diff --git a/demo/comparison-oskar/main.cpp b/demo/comparison-oskar/main.cpp
new file mode 100644
index 00000000..d81e8361
--- /dev/null
+++ b/demo/comparison-oskar/main.cpp
@@ -0,0 +1,38 @@
+#include <cmath>
+#include <iostream>
+#include <cstdlib>
+
+#include <OSKARElementResponse.h>
+
+#include "npy.hpp" // to save arrays in numpy format
+
+int main(int argc, char** argv){
+// int main() {
+ everybeam::OSKARElementResponseSphericalWave element_response("oskar.h5");
+ double freq = 50e6;
+
+ int N;
+ if (argc == 1){
+ N = 256;
+ }
+ else{
+ N = atoi(argv[1]);
+ }
+
+ std::vector<std::complex<double>> result(N*N*2*2);
+ typedef std::complex<double>result_arr_t[N][N][2][2];
+ result_arr_t &result_arr = * (result_arr_t*) result.data();
+
+ for(int i=0; i<N; ++i) {
+ double x = (2.0*i)/(N-1) - 1.0;
+ for(int j=0; j<N; ++j) {
+ double y = (2.0*j)/(N-1) - 1.0;
+ double theta = asin(sqrt(x*x + y*y));
+ double phi = atan2(y,x);
+ element_response.response(0, freq, theta, phi, result_arr[i][j]);
+ }
+ }
+
+ const long unsigned leshape [] = {(long unsigned int) N, N, 2, 2};
+ npy::SaveArrayAsNumpy("response.npy", false, 4, leshape, result);
+}
diff --git a/demo/comparison-oskar/npy.hpp b/demo/comparison-oskar/npy.hpp
new file mode 100644
index 00000000..f8dbf6ca
--- /dev/null
+++ b/demo/comparison-oskar/npy.hpp
@@ -0,0 +1,516 @@
+/*
+ Copyright 2017 Leon Merten Lohse
+
+ Permission is hereby granted, free of charge, to any person obtaining a copy
+ of this software and associated documentation files (the "Software"), to deal
+ in the Software without restriction, including without limitation the rights
+ to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ copies of the Software, and to permit persons to whom the Software is
+ furnished to do so, subject to the following conditions:
+
+ The above copyright notice and this permission notice shall be included in
+ all copies or substantial portions of the Software.
+
+ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+ SOFTWARE.
+*/
+
+#ifndef NPY_H
+#define NPY_H
+
+#include <complex>
+#include <fstream>
+#include <string>
+#include <iostream>
+#include <sstream>
+#include <cstdint>
+#include <cstring>
+#include <vector>
+#include <stdexcept>
+#include <algorithm>
+#include <regex>
+#include <unordered_map>
+
+
+namespace npy {
+
+/* Compile-time test for byte order.
+ If your compiler does not define these per default, you may want to define
+ one of these constants manually.
+ Defaults to little endian order. */
+#if defined(__BYTE_ORDER) && __BYTE_ORDER == __BIG_ENDIAN || \
+ defined(__BIG_ENDIAN__) || \
+ defined(__ARMEB__) || \
+ defined(__THUMBEB__) || \
+ defined(__AARCH64EB__) || \
+ defined(_MIBSEB) || defined(__MIBSEB) || defined(__MIBSEB__)
+const bool big_endian = true;
+#else
+const bool big_endian = false;
+#endif
+
+
+const char magic_string[] = "\x93NUMPY";
+const size_t magic_string_length = 6;
+
+const char little_endian_char = '<';
+const char big_endian_char = '>';
+const char no_endian_char = '|';
+
+constexpr char host_endian_char = ( big_endian ?
+ big_endian_char :
+ little_endian_char );
+
+/* npy array length */
+typedef unsigned long int ndarray_len_t;
+
+inline void write_magic(std::ostream& ostream, unsigned char v_major=1, unsigned char v_minor=0) {
+ ostream.write(magic_string, magic_string_length);
+ ostream.put(v_major);
+ ostream.put(v_minor);
+}
+
+inline void read_magic(std::istream& istream, unsigned char& v_major, unsigned char& v_minor) {
+ char buf[magic_string_length+2];
+ istream.read(buf, magic_string_length+2);
+
+ if(!istream) {
+ throw std::runtime_error("io error: failed reading file");
+ }
+
+ if (0 != std::memcmp(buf, magic_string, magic_string_length))
+ throw std::runtime_error("this file does not have a valid npy format.");
+
+ v_major = buf[magic_string_length];
+ v_minor = buf[magic_string_length+1];
+}
+
+// typestring magic
+struct Typestring {
+ private:
+ char c_endian;
+ char c_type;
+ int len;
+
+ public:
+ inline std::string str() {
+ const size_t max_buflen = 16;
+ char buf[max_buflen];
+ std::sprintf(buf, "%c%c%u", c_endian, c_type, len);
+ return std::string(buf);
+ }
+
+ Typestring(const std::vector<float>& v)
+ :c_endian {host_endian_char}, c_type {'f'}, len {sizeof(float)} {}
+ Typestring(const std::vector<double>& v)
+ :c_endian {host_endian_char}, c_type {'f'}, len {sizeof(double)} {}
+ Typestring(const std::vector<long double>& v)
+ :c_endian {host_endian_char}, c_type {'f'}, len {sizeof(long double)} {}
+
+ Typestring(const std::vector<char>& v)
+ :c_endian {no_endian_char}, c_type {'i'}, len {sizeof(char)} {}
+ Typestring(const std::vector<short>& v)
+ :c_endian {host_endian_char}, c_type {'i'}, len {sizeof(short)} {}
+ Typestring(const std::vector<int>& v)
+ :c_endian {host_endian_char}, c_type {'i'}, len {sizeof(int)} {}
+ Typestring(const std::vector<long>& v)
+ :c_endian {host_endian_char}, c_type {'i'}, len {sizeof(long)} {}
+ Typestring(const std::vector<long long>& v) :c_endian {host_endian_char}, c_type {'i'}, len {sizeof(long long)} {}
+
+ Typestring(const std::vector<unsigned char>& v)
+ :c_endian {no_endian_char}, c_type {'u'}, len {sizeof(unsigned char)} {}
+ Typestring(const std::vector<unsigned short>& v)
+ :c_endian {host_endian_char}, c_type {'u'}, len {sizeof(unsigned short)} {}
+ Typestring(const std::vector<unsigned int>& v)
+ :c_endian {host_endian_char}, c_type {'u'}, len {sizeof(unsigned int)} {}
+ Typestring(const std::vector<unsigned long>& v)
+ :c_endian {host_endian_char}, c_type {'u'}, len {sizeof(unsigned long)} {}
+ Typestring(const std::vector<unsigned long long>& v)
+ :c_endian {host_endian_char}, c_type {'u'}, len {sizeof(unsigned long long)} {}
+
+ Typestring(const std::vector<std::complex<float>>& v)
+ :c_endian {host_endian_char}, c_type {'c'}, len {sizeof(std::complex<float>)} {}
+ Typestring(const std::vector<std::complex<double>>& v)
+ :c_endian {host_endian_char}, c_type {'c'}, len {sizeof(std::complex<double>)} {}
+ Typestring(const std::vector<std::complex<long double>>& v)
+ :c_endian {host_endian_char}, c_type {'c'}, len {sizeof(std::complex<long double>)} {}
+};
+
+inline void parse_typestring( std::string typestring){
+ std::regex re ("'([<>|])([ifuc])(\\d+)'");
+ std::smatch sm;
+
+ std::regex_match(typestring, sm, re );
+
+ if ( sm.size() != 4 ) {
+ throw std::runtime_error("invalid typestring");
+ }
+}
+
+namespace pyparse {
+
+/**
+ Removes leading and trailing whitespaces
+ */
+inline std::string trim(const std::string& str) {
+ const std::string whitespace = " \t";
+ auto begin = str.find_first_not_of(whitespace);
+
+ if (begin == std::string::npos)
+ return "";
+
+ auto end = str.find_last_not_of(whitespace);
+
+ return str.substr(begin, end-begin+1);
+}
+
+
+inline std::string get_value_from_map(const std::string& mapstr) {
+ size_t sep_pos = mapstr.find_first_of(":");
+ if (sep_pos == std::string::npos)
+ return "";
+
+ std::string tmp = mapstr.substr(sep_pos+1);
+ return trim(tmp);
+}
+
+/**
+ Parses the string representation of a Python dict
+
+ The keys need to be known and may not appear anywhere else in the data.
+ */
+inline std::unordered_map<std::string, std::string> parse_dict(std::string in, std::vector<std::string>& keys) {
+
+ std::unordered_map<std::string, std::string> map;
+
+ if (keys.size() == 0)
+ return map;
+
+ in = trim(in);
+
+ // unwrap dictionary
+ if ((in.front() == '{') && (in.back() == '}'))
+ in = in.substr(1, in.length()-2);
+ else
+ throw std::runtime_error("Not a Python dictionary.");
+
+ std::vector<std::pair<size_t, std::string>> positions;
+
+ for (auto const& value : keys) {
+ size_t pos = in.find( "'" + value + "'" );
+
+ if (pos == std::string::npos)
+ throw std::runtime_error("Missing '"+value+"' key.");
+
+ std::pair<size_t, std::string> position_pair { pos, value };
+ positions.push_back(position_pair);
+ }
+
+ // sort by position in dict
+ std::sort(positions.begin(), positions.end() );
+
+ for(size_t i = 0; i < positions.size(); ++i) {
+ std::string raw_value;
+ size_t begin { positions[i].first };
+ size_t end { std::string::npos };
+
+ std::string key = positions[i].second;
+
+ if ( i+1 < positions.size() )
+ end = positions[i+1].first;
+
+ raw_value = in.substr(begin, end-begin);
+
+ raw_value = trim(raw_value);
+
+ if (raw_value.back() == ',')
+ raw_value.pop_back();
+
+ map[key] = get_value_from_map(raw_value);
+ }
+
+ return map;
+}
+
+/**
+ Parses the string representation of a Python boolean
+ */
+inline bool parse_bool(const std::string& in) {
+ if (in == "True")
+ return true;
+ if (in == "False")
+ return false;
+
+ throw std::runtime_error("Invalid python boolan.");
+}
+
+/**
+ Parses the string representation of a Python str
+ */
+inline std::string parse_str(const std::string& in) {
+ if ((in.front() == '\'') && (in.back() == '\''))
+ return in.substr(1, in.length()-2);
+
+ throw std::runtime_error("Invalid python string.");
+}
+
+/**
+ Parses the string represenatation of a Python tuple into a vector of its items
+ */
+inline std::vector<std::string> parse_tuple(std::string in) {
+ std::vector<std::string> v;
+ const char seperator = ',';
+
+ in = trim(in);
+
+ if ((in.front() == '(') && (in.back() == ')'))
+ in = in.substr(1, in.length()-2);
+ else
+ throw std::runtime_error("Invalid Python tuple.");
+
+ std::istringstream iss(in);
+
+ for (std::string token; std::getline(iss, token, seperator);) {
+ v.push_back(token);
+ }
+
+ return v;
+}
+
+template <typename T>
+inline std::string write_tuple(const std::vector<T>& v) {
+ if (v.size() == 0)
+ return "";
+
+ std::ostringstream ss;
+
+ if (v.size() == 1) {
+ ss << "(" << v.front() << ",)";
+ } else {
+ const std::string delimiter = ", ";
+ // v.size() > 1
+ ss << "(";
+ std::copy(v.begin(), v.end()-1, std::ostream_iterator<T>(ss, delimiter.c_str()));
+ ss << v.back();
+ ss << ")";
+ }
+
+ return ss.str();
+}
+
+inline std::string write_boolean(bool b) {
+ if(b)
+ return "True";
+ else
+ return "False";
+}
+
+} // namespace pyparse
+
+
+inline void parse_header(std::string header, std::string& descr, bool& fortran_order, std::vector<ndarray_len_t>& shape) {
+ /*
+ The first 6 bytes are a magic string: exactly "x93NUMPY".
+ The next 1 byte is an unsigned byte: the major version number of the file format, e.g. x01.
+ The next 1 byte is an unsigned byte: the minor version number of the file format, e.g. x00. Note: the version of the file format is not tied to the version of the numpy package.
+ The next 2 bytes form a little-endian unsigned short int: the length of the header data HEADER_LEN.
+ The next HEADER_LEN bytes form the header data describing the array's format. It is an ASCII string which contains a Python literal expression of a dictionary. It is terminated by a newline ('n') and padded with spaces ('x20') to make the total length of the magic string + 4 + HEADER_LEN be evenly divisible by 16 for alignment purposes.
+ The dictionary contains three keys:
+
+ "descr" : dtype.descr
+ An object that can be passed as an argument to the numpy.dtype() constructor to create the array's dtype.
+ "fortran_order" : bool
+ Whether the array data is Fortran-contiguous or not. Since Fortran-contiguous arrays are a common form of non-C-contiguity, we allow them to be written directly to disk for efficiency.
+ "shape" : tuple of int
+ The shape of the array.
+ For repeatability and readability, this dictionary is formatted using pprint.pformat() so the keys are in alphabetic order.
+ */
+
+ // remove trailing newline
+ if (header.back() != '\n')
+ throw std::runtime_error("invalid header");
+ header.pop_back();
+
+ // parse the dictionary
+ std::vector<std::string> keys { "descr", "fortran_order", "shape" };
+ auto dict_map = npy::pyparse::parse_dict(header, keys);
+
+ if (dict_map.size() == 0)
+ throw std::runtime_error("invalid dictionary in header");
+
+ std::string descr_s = dict_map["descr"];
+ std::string fortran_s = dict_map["fortran_order"];
+ std::string shape_s = dict_map["shape"];
+
+ // TODO: extract info from typestring
+ parse_typestring(descr_s);
+ // remove
+ descr = npy::pyparse::parse_str(descr_s);
+
+ // convert literal Python bool to C++ bool
+ fortran_order = npy::pyparse::parse_bool(fortran_s);
+
+ // parse the shape tuple
+ auto shape_v = npy::pyparse::parse_tuple(shape_s);
+ if (shape_v.size() == 0)
+ throw std::runtime_error("invalid shape tuple in header");
+
+ for ( auto item : shape_v ) {
+ ndarray_len_t dim = static_cast<ndarray_len_t>(std::stoul(item));
+ shape.push_back(dim);
+ }
+}
+
+
+inline std::string write_header_dict(const std::string& descr, bool fortran_order, const std::vector<ndarray_len_t>& shape) {
+ std::string s_fortran_order = npy::pyparse::write_boolean(fortran_order);
+ std::string shape_s = npy::pyparse::write_tuple(shape);
+
+ return "{'descr': '" + descr + "', 'fortran_order': " + s_fortran_order + ", 'shape': " + shape_s + ", }";
+}
+
+inline void write_header(std::ostream& out, const std::string& descr, bool fortran_order, const std::vector<ndarray_len_t>& shape_v)
+{
+ std::string header_dict = write_header_dict(descr, fortran_order, shape_v);
+
+ size_t length = magic_string_length + 2 + 2 + header_dict.length() + 1;
+
+ unsigned char version[2] = {1, 0};
+ if (length >= 255*255) {
+ length = magic_string_length + 2 + 4 + header_dict.length() + 1;
+ version[0] = 2;
+ version[1] = 0;
+ }
+ size_t padding_len = 16 - length % 16;
+ std::string padding (padding_len, ' ');
+
+ // write magic
+ write_magic(out, version[0], version[1]);
+
+ // write header length
+ if (version[0] == 1 && version[1] == 0) {
+ char header_len_le16[2];
+ uint16_t header_len = header_dict.length() + padding.length() + 1;
+
+ header_len_le16[0] = (header_len >> 0) & 0xff;
+ header_len_le16[1] = (header_len >> 8) & 0xff;
+ out.write(reinterpret_cast<char *>(header_len_le16), 2);
+ }else{
+ char header_len_le32[4];
+ uint32_t header_len = header_dict.length() + padding.length() + 1;
+
+ header_len_le32[0] = (header_len >> 0) & 0xff;
+ header_len_le32[1] = (header_len >> 8) & 0xff;
+ header_len_le32[2] = (header_len >> 16) & 0xff;
+ header_len_le32[3] = (header_len >> 24) & 0xff;
+ out.write(reinterpret_cast<char *>(header_len_le32), 4);
+ }
+
+ out << header_dict << padding << '\n';
+}
+
+inline std::string read_header(std::istream& istream) {
+ // check magic bytes an version number
+ unsigned char v_major, v_minor;
+ read_magic(istream, v_major, v_minor);
+
+ uint32_t header_length;
+ if(v_major == 1 && v_minor == 0){
+
+ char header_len_le16[2];
+ istream.read(header_len_le16, 2);
+ header_length = (header_len_le16[0] << 0) | (header_len_le16[1] << 8);
+
+ if((magic_string_length + 2 + 2 + header_length) % 16 != 0) {
+ // TODO: display warning
+ }
+ }else if(v_major == 2 && v_minor == 0) {
+ char header_len_le32[4];
+ istream.read(header_len_le32, 4);
+
+ header_length = (header_len_le32[0] << 0) | (header_len_le32[1] << 8)
+ | (header_len_le32[2] << 16) | (header_len_le32[3] << 24);
+
+ if((magic_string_length + 2 + 4 + header_length) % 16 != 0) {
+ // TODO: display warning
+ }
+ }else{
+ throw std::runtime_error("unsupported file format version");
+ }
+
+ auto buf_v = std::vector<char>();
+ buf_v.reserve(header_length);
+ istream.read(buf_v.data(), header_length);
+ std::string header(buf_v.data(), header_length);
+
+ return header;
+}
+
+inline ndarray_len_t comp_size(const std::vector<ndarray_len_t>& shape) {
+ ndarray_len_t size = 1;
+ for (ndarray_len_t i : shape )
+ size *= i;
+
+ return size;
+}
+
+template<typename Scalar>
+inline void SaveArrayAsNumpy( const std::string& filename, bool fortran_order, unsigned int n_dims, const unsigned long shape[], const std::vector<Scalar>& data)
+{
+ Typestring typestring_o(data);
+ std::string typestring = typestring_o.str();
+
+ std::ofstream stream( filename, std::ofstream::binary);
+ if(!stream) {
+ throw std::runtime_error("io error: failed to open a file.");
+ }
+
+ std::vector<ndarray_len_t> shape_v(shape, shape+n_dims);
+ write_header(stream, typestring, fortran_order, shape_v);
+
+ auto size = static_cast<size_t>(comp_size(shape_v));
+
+ stream.write(reinterpret_cast<const char*>(data.data()), sizeof(Scalar) * size);
+}
+
+
+template<typename Scalar>
+inline void LoadArrayFromNumpy(const std::string& filename, std::vector<unsigned long>& shape, std::vector<Scalar>& data)
+{
+ std::ifstream stream(filename, std::ifstream::binary);
+ if(!stream) {
+ throw std::runtime_error("io error: failed to open a file.");
+ }
+
+ std::string header = read_header(stream);
+
+ // parse header
+ bool fortran_order;
+ std::string typestr;
+
+ parse_header(header, typestr, fortran_order, shape);
+
+ // check if the typestring matches the given one
+ Typestring typestring_o {data};
+ std::string expect_typestr = typestring_o.str();
+ if (typestr != expect_typestr) {
+ throw std::runtime_error("formatting error: typestrings not matching");
+ }
+
+
+ // compute the data size based on the shape
+ auto size = static_cast<size_t>(comp_size(shape));
+ data.resize(size);
+
+ // read the data
+ stream.read(reinterpret_cast<char*>(data.data()), sizeof(Scalar)*size);
+}
+
+} // namespace npy
+
+#endif // NPY_H
diff --git a/demo/comparison-oskar/read_oskar_beams.py b/demo/comparison-oskar/read_oskar_beams.py
new file mode 100644
index 00000000..2c6096cf
--- /dev/null
+++ b/demo/comparison-oskar/read_oskar_beams.py
@@ -0,0 +1,85 @@
+from matplotlib import pyplot as plt
+from astropy.io import fits
+import numpy as np
+
+
+def read_oskar_beams():
+ """
+ Read oskar beam and store into
+
+ Returns
+ -------
+ np.array
+ (N, N, 2, 2) array, where N number of image pixels
+ """
+ A = None
+ for i, pol1 in enumerate(["X", "Y"]):
+ for j, pol2 in enumerate(["X", "Y"]):
+ filename_amp = "basefunctions_050_MHz_S0000_TIME_SEP_CHAN_SEP_AMP_{}{}.fits".format(
+ pol1, pol2
+ )
+ filename_phase = "basefunctions_050_MHz_S0000_TIME_SEP_CHAN_SEP_PHASE_{}{}.fits".format(
+ pol1, pol2
+ )
+ d = (
+ fits.getdata(filename_amp) * np.exp(1j * fits.getdata(filename_phase))
+ )[0, 0, :, :].T
+ if A is None:
+ N = d.shape[-1]
+ A = np.zeros((N, N, 2, 2), dtype=np.complex128)
+ A[:, :, i, j] = d
+ N = A.shape[0]
+
+ for i in range(N):
+ x = 2.0 * i / (N - 1) - 1.0
+ for j in range(N):
+ y = 2.0 * j / (N - 1) - 1.0
+ theta = np.arcsin(np.sqrt(x * x + y * y))
+ phi = np.arctan2(y, x)
+
+ # Need to swap the sign of phi to get the transformation right
+ # Still need to figure out why
+ phi = -phi
+
+ e_theta = np.array([[np.cos(phi)], [np.sin(phi)]])
+ e_phi = np.array([[-np.sin(phi)], [np.cos(phi)]])
+
+ # This matrix applies the transformation defined in
+ # https://github.com/OxfordSKA/OSKAR/blob/master/oskar/convert/define_convert_ludwig3_to_theta_phi_components.h
+ T = np.concatenate((e_theta, e_phi), axis=1)
+
+ # Transformation is applied to the right side of A
+ # so the rows a of A are tranformed
+ A[i, j, :, :] = np.dot(A[i, j, :, :], T)
+ return A
+
+
+if __name__ == "__main__":
+
+ A = read_oskar_beams()
+
+ plt.subplot(2, 2, 1)
+ plt.imshow(np.abs(A[:, :, 0, 0]).T, clim=(0, 0.25), origin="lower")
+ plt.colorbar()
+ plt.title("abs(Etheta)")
+
+ plt.subplot(2, 2, 2)
+ plt.imshow(np.abs(A[:, :, 0, 1]).T, clim=(0, 0.25), origin="lower")
+ plt.colorbar()
+ plt.title("abs(Ephi)")
+
+ plt.subplot(2, 2, 3)
+ plt.imshow(
+ np.angle(A[:, :, 0, 0]).T, clim=(-np.pi, np.pi), cmap="twilight", origin="lower"
+ )
+ plt.colorbar()
+ plt.title("angle(Etheta)")
+
+ plt.subplot(2, 2, 4)
+ plt.imshow(
+ np.angle(A[:, :, 0, 1]).T, clim=(-np.pi, np.pi), cmap="twilight", origin="lower"
+ )
+ plt.colorbar()
+ plt.title("angle(Ephi)")
+
+ plt.show()
diff --git a/demo/comparison-oskar/run_oskar.py b/demo/comparison-oskar/run_oskar.py
new file mode 100644
index 00000000..c1740178
--- /dev/null
+++ b/demo/comparison-oskar/run_oskar.py
@@ -0,0 +1,170 @@
+#!/usr/bin/env python3
+"""
+Generates all-sky zenith-centred beam patterns for SKALA-4 and EDA-2 antennas.
+"""
+
+import copy
+import glob
+import os.path
+import re
+import subprocess
+
+from astropy.io import fits
+from astropy.time import Time, TimeDelta
+import matplotlib
+matplotlib.use('Agg')
+# pylint: disable=wrong-import-position
+import matplotlib.pyplot as plt
+import numpy
+import oskar
+
+
+def get_start_time(ra0_deg, length_sec):
+ """Returns optimal start time for field RA and observation length."""
+ t = Time('2000-01-01 00:00:00', scale='utc', location=('116.764d', '0d'))
+ dt_hours = (24.0 - t.sidereal_time('apparent').hour) / 1.0027379
+ dt_hours += (ra0_deg / 15.0)
+ start = t + TimeDelta(dt_hours * 3600.0 - length_sec / 2.0, format='sec')
+ return start.value
+
+
+def plot_panel(ax, image, title, cmap):
+ """Plots a single panel."""
+ im = ax.imshow(numpy.squeeze(image), cmap=cmap)
+ plt.colorbar(im, format='%.2e')
+ plt.tick_params(labelcolor='none', top='off', bottom='off', left='off', right='off')
+ ax.set_xticks([])
+ ax.set_yticks([])
+ ax.invert_yaxis()
+ ax.set_title(title)
+ ax.axis('equal')
+
+
+def make_plots(title, glob_pattern, out_basename):
+ """Generates a plot with four panels."""
+ # Load FITS images matching the glob pattern.
+ files = glob.glob(glob_pattern)
+ files.sort() # Must be sorted!
+ images = []
+ for file in files:
+ images.append(fits.getdata(file))
+
+ # Set titles and colour maps to use.
+ cmap = ''
+ titles = []
+ if '_AUTO_POWER' in glob_pattern:
+ # cmap = 'plasma'
+ cmap = 'Blues_r'
+ titles = ['Stokes I', 'Stokes Q', 'Stokes U', 'Stokes V']
+ elif '_AMP' in glob_pattern:
+ # cmap = 'jet'
+ cmap = 'CMRmap'
+ titles = ['XX', 'XY', 'YX', 'YY']
+
+ # Sub-plots, one for each Stokes parameter.
+ fig = plt.figure(figsize=(8, 6))
+ ax = fig.add_subplot(221, frameon=False)
+ plot_panel(ax, images[0], titles[0], cmap)
+ ax = fig.add_subplot(222, frameon=False)
+ plot_panel(ax, images[1], titles[1], cmap)
+ ax = fig.add_subplot(223, frameon=False)
+ plot_panel(ax, images[2], titles[2], cmap)
+ ax = fig.add_subplot(224, frameon=False)
+ plot_panel(ax, images[3], titles[3], cmap)
+
+ # Add main title.
+ title = title.replace('_', ' ') # Replace underscores with spaces.
+ fig.suptitle(title)
+ fig.tight_layout()
+ fig.subplots_adjust(top=0.88)
+
+ # Save and close.
+ plt.savefig('%s.png' % (out_basename))
+ plt.close('all')
+
+
+def main(npixels):
+ """Main function."""
+ # Name of the application to run, and a settings file for it.
+ app = 'oskar_sim_beam_pattern'
+ settings_path = '_temp_settings.ini'
+
+ # Define some basic observation parameters.
+ ra0_deg = 0.0
+ dec0_deg = -27.0
+ length_sec = 1.0
+
+ # Define base settings dictionary.
+ common_settings = {
+ 'observation': {
+ 'phase_centre_ra_deg': ra0_deg,
+ 'phase_centre_dec_deg': dec0_deg,
+ #'pointing_file': 'station_pointing.txt',
+ 'start_time_utc': get_start_time(ra0_deg, length_sec),
+ 'length': length_sec
+ },
+ 'telescope': {
+ 'normalise_beams_at_phase_centre': False,
+ 'aperture_array/element_pattern/normalise': False
+ },
+ 'beam_pattern': {
+ 'coordinate_frame': 'Horizon',
+ 'beam_image/size': npixels,
+ 'station_outputs/fits_image/amp': True,
+ 'station_outputs/fits_image/phase': True,
+ 'station_outputs/fits_image/auto_power_real': False
+ }
+ }
+
+ # Define frequencies of interest (in MHz).
+ freqs = [50]
+ #freqs = [70]
+
+ # Define telescope models to use, and associated overrides for them.
+ telescopes = {
+ 'basefunctions': {
+ 'telescope/input_directory': 'telescope.tm',
+ 'telescope/aperture_array/element_pattern/enable_numerical': True,
+ 'telescope/aperture_array/element_pattern/swap_xy': True,
+ 'telescope/aperture_array/array_pattern/enable': False
+ },
+ }
+
+ # Loop over telescope models.
+ for tel, tel_params in telescopes.items():
+
+ # Copy the base settings dictionary.
+ current_settings = copy.deepcopy(common_settings)
+
+ # Update current settings with telescope model parameters.
+ current_settings.update(tel_params)
+
+ # Loop over frequencies.
+ for freq in freqs:
+
+ # Create the settings file.
+ open(settings_path, 'w').close()
+ settings = oskar.SettingsTree(app, settings_path)
+ settings.from_dict(current_settings)
+
+ # Update output root path and frequency.
+ tel_root = re.sub(r'[^\w]', '', tel) # Strip symbols from tel.
+ root_path = tel_root + ('_%03d_MHz' % freq)
+ settings['beam_pattern/root_path'] = root_path
+ settings['observation/start_frequency_hz'] = 1e6 * freq
+
+ # Run the app with the settings file.
+ subprocess.call([app, settings_path])
+
+ # Make plots.
+ #title = tel + ' @ ' + str(freq) + ' MHz'
+ #make_plots(title=title,
+ #glob_pattern=root_path+'*_AMP*',
+ #out_basename=root_path+'_amp')
+ #make_plots(title=title,
+ #glob_pattern=root_path+'*_AUTO_POWER*',
+ #out_basename=root_path+'_stokes')
+
+
+if __name__ == '__main__':
+ main()
diff --git a/demo/comparison-oskar/telescope.tm/layout.txt b/demo/comparison-oskar/telescope.tm/layout.txt
new file mode 100644
index 00000000..3d7fd173
--- /dev/null
+++ b/demo/comparison-oskar/telescope.tm/layout.txt
@@ -0,0 +1,2 @@
+0,0
+
diff --git a/demo/comparison-oskar/telescope.tm/position.txt b/demo/comparison-oskar/telescope.tm/position.txt
new file mode 100644
index 00000000..7711436a
--- /dev/null
+++ b/demo/comparison-oskar/telescope.tm/position.txt
@@ -0,0 +1 @@
+0 -50
diff --git a/demo/comparison-oskar/telescope.tm/station000/layout.txt b/demo/comparison-oskar/telescope.tm/station000/layout.txt
new file mode 100644
index 00000000..3d7fd173
--- /dev/null
+++ b/demo/comparison-oskar/telescope.tm/station000/layout.txt
@@ -0,0 +1,2 @@
+0,0
+
diff --git a/scripts/coeff_scripts/oskar_csv_to_hdf5.py b/scripts/coeff_scripts/oskar_csv_to_hdf5.py
new file mode 100755
index 00000000..becae4df
--- /dev/null
+++ b/scripts/coeff_scripts/oskar_csv_to_hdf5.py
@@ -0,0 +1,93 @@
+#!/usr/bin/env python3
+
+import argparse
+import os
+import scipy.io as sio
+import os.path
+import numpy as np
+import h5py
+from tqdm import tqdm
+from glob import glob
+import csv
+
+# Parse command-line arguments
+parser = argparse.ArgumentParser()
+parser.add_argument("input_dir", type=str, help="directory containing coefficients")
+parser.add_argument("output_file", type=str, help="name of the HDF5 output file")
+parser.add_argument("-debug", dest="debug", action="store_true")
+parser.set_defaults(input_dir=".")
+parser.set_defaults(output_file="oskar.h5")
+parser.set_defaults(debug=False)
+args = parser.parse_args()
+
+# Read the input directory
+# This directory should contain subdirectories,
+# with names corresponding to frequency.
+input_dir = args.input_dir
+print("Reading coefficients from: " , input_dir)
+
+# Only read element_id 0, because EveryBeam's implementation of the OSKAR model
+# currently supports only a single model for all elements
+element_id = 0
+nr_elements = 1
+
+files = glob(os.path.join(input_dir, 'element_pattern_spherical_wave_?_t?_??_{}_*.txt').format(element_id))
+
+files = [os.path.basename(f) for f in files]
+
+freqs = sorted(set([int(f[:-4].split('_')[-1]) for f in files]))
+
+
+print("Parsing frequencies %d-%d MHz" % (int(freqs[0]), int(freqs[-1])))
+
+# Create HDF5 file
+output_file = args.output_file
+print("Creating HDF5 file: ", output_file)
+hf = h5py.File(output_file, 'w')
+
+# Debugging
+debug = args.debug
+if (debug):
+ freqs = [freqs[0]]
+ tqdm = lambda x: x
+
+# Parse all freqs
+for freq in tqdm(freqs):
+
+ A = None
+
+ for i, pol in enumerate(('x', 'y')):
+ for j, em in enumerate(('e', 'm')):
+ for reim, s in zip(('re', 'im'), (1.0, 1.0j)) :
+ filename = 'element_pattern_spherical_wave_{}_t{}_{}_{}_{}.txt'.format(pol, em, reim, element_id, freq)
+ alpha = np.loadtxt(os.path.join(input_dir,filename), delimiter=',', ndmin=2)
+ l_max = alpha.shape[0]
+ assert(alpha.shape[1] == 2*l_max + 1)
+ if A is None:
+ A = np.zeros( (l_max,2*l_max+1, 2,2), dtype=np.complex128)
+ else:
+ assert(A.shape == (l_max,2*l_max+1, 2,2))
+ A[:,:,i,j] += alpha*s
+
+ # Create the output matrix
+ data = np.zeros((nr_elements, l_max * (l_max+2), 4), dtype=np.complex128)
+
+ # Fill the output matrix
+ # The inner dimension are set as follows:
+ # (x_te_re, x_te_im), (x_tm_re, x_tm_im),
+ # (y_te_re, y_te_im), (y_tm_re, y_tm_im)
+
+ i = 0
+ for l in range(l_max):
+ n = (l+1) * 2 + 1
+ data[0,i:i+n,:] = A[l,:n,:].reshape(n,4)
+ i += n
+
+ # Add group for current frequency
+ hf.create_dataset(str(freq), data.shape, data=data)
+
+# Close HDF5 file
+hf.close()
+
+# Done
+print("Done!")
--
GitLab