-
Notifications
You must be signed in to change notification settings - Fork 0
/
export_R.qmd
172 lines (129 loc) · 6.16 KB
/
export_R.qmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
---
title: "Export to '.tmat.h5' format with R"
author: "baptiste"
date: today
engine: knitr
---
This document showcases a basic R script to export T-matrices in the `.tmat.h5` HDF5 format. For illustration, we start by producing a dummy dataset. This minimal reproducible example file is available for download as a standalone script: [export_R.R](export_R.R).
## Mockup input data
Consistent with the other examples, we generate a `50x30x30` array of dummy data, which repeats 50 times (50 wavelengths) a matrix of 900 entries ranging from $1+1i$ (first element, top left) to $900 + 900i$ (last element, bottom right). Note the expectation of **row-major** ordering in HDF5. The `3x3` top-left block is
```
1.0000 + 1.0000i 2.0000 + 2.0000i 3.0000 + 3.0000i ...
31.0000 +31.0000i 32.0000 +32.0000i 33.0000 +33.0000i ...
61.0000 +61.0000i 62.0000 +62.0000i 63.0000 +63.0000i ...
... ... ... ...
```
The `rhdf5` package has support for a variety of R objects, including lists which are automatically written as grouped objects in HDF5.
```{r, eval=FALSE}
library(rhdf5) # note: dev version to support complex
library(uuid) # uuid
library(glue) # string interpolation
library(purrr) # mapping functions
## dummy data
# possibly multiple wavelengths
wavelength <- seq(400, 800, by=50)
Nl <- length(wavelength)
Lmax <- 3
qmax <- 2*(Lmax*(Lmax+1)+Lmax) # T-matrix size
# dummy 30x30 matrix values for each wavelength
# note the byrow due to HDF5 expecting
# row-major ordering vs R's default column-major
tdata <- matrix(1:qmax^2 + 1i*(1:qmax^2), qmax, qmax, byrow=TRUE)
tmatrix <- array(NA_complex_, c(Nl,qmax,qmax))
for(i in seq_len(Nl))
tmatrix[i,,] <- tdata
tmatrix[1,1:3,1:3]
modes <- list(l = rep(NA_integer_, qmax),
m = rep(NA_integer_, qmax),
polarization = rep(NA_character_, qmax))
i <- 1
for (li in 1:Lmax){
for (mi in -li:li){
for (si in c("electric", "magnetic")){
modes$l[i] <- li
modes$m[i] <- mi
modes$polarization[i] <- si
i <- i+1
}
}
}
embedding <- list('relative_permeability' = 1.0,
'relative_permittivity' = 1.33^2)
scatterer <- list(material = list('relative_permeability' = 1.0,
'relative_permittivity' =
rep(-11.4+1.181i,Nl)),
geometry = list('radiusxy' = 20.0,
'radiusz' = 40.0))
computation <- list(method_parameters = list('Lmax' = Lmax,
'Ntheta' = 100),
files = list(script = paste(readLines('export_R.R'), collapse = "\n")))
```
## Saving to HDF5
The `rhdf5` package provides support for reading/writing HDF5 files into/from R objects. Until my [recent request](https://support.bioconductor.org/p/9156305/#9156408), complex arrays were not supported, but this is now implemented in the dev version of the package.
```{r, eval=FALSE}
library(rhdf5) # note: requires dev version to support complex arrays
# cf https://support.bioconductor.org/p/9156305/#9156408
# install.packages("remotes")
# remotes::install_github("grimbough/rhdf5@devel")
# may need 'crypto' library from openSSL to compile
# e.g. via brew install openssl on macos
```
We can then write the different objects defined above using `hwrite`. Note the important `native = TRUE` parameter for the T-matrix data: R stores arrays column-wise, while HDF5 uses row-major ordering. To avoid confusion between different programming languages, we suggest sticking with the native HDF5 convention (`native = TRUE` ensures that the array is written transposed).
Attributes don't seem to have a convenient high-level interface, and unfortunately they need to be written slightly differently for the root level, datasets, and groups.
```{r, eval=FALSE}
f <- 'ar.tmat.h5'
software = sprintf("SMARTIES=1.1, R=%s, rhdf5=%s", paste0(R.version$major, R.version$minor), packageVersion("rhdf5"))
unlink(f) # delete previous file if it exists
h5createFile(f)
h5closeAll() # in case connections open
h5write(wavelength, file=f, name="/vacuum_wavelength")
h5write(tmatrix, file=f, name="/tmatrix", native=TRUE) # store rowwise
h5write(modes, file=f, name='/modes')
h5write(embedding, file=f, name="/embedding")
h5write(scatterer, file=f, name="/scatterer")
h5write(computation, file=f, name='/computation')
## write attributes
fid <- H5Fopen(f)
# root level
h5writeAttribute("Au prolate spheroid in water", fid, "name")
h5writeAttribute("Computation using SMARTIES, a numerically robust EBCM implementation for spheroids", fid, "description")
h5writeAttribute("gold, spheroid, ebcm, passive, reciprocal, czinfinity, mirrorxyz", fid, "keywords")
h5writeAttribute("v0.01", fid, "storage_format_version")
# wavelength
did <- H5Dopen(fid, "vacuum_wavelength")
h5writeAttribute("nm", did, "unit")
H5Dclose(did)
# embedding
gid <- H5Gopen(fid, "embedding")
h5writeAttribute("H2O, Water", gid, "name")
h5writeAttribute("non-dispersive", gid, "keywords")
H5Gclose(gid)
# material
gid <- H5Gopen(fid, "scatterer/material")
h5writeAttribute("Au, Gold", gid, "name")
h5writeAttribute("dispersive, plasmonic", gid, "keywords")
h5writeAttribute("Au from Raschke et al 10.1103/PhysRevB.86.235147", gid, "reference")
H5Gclose(gid)
# geometry
gid <- H5Gopen(fid, "scatterer/geometry")
h5writeAttribute("nm", gid, "unit")
h5writeAttribute("spheroid", gid, "shape")
h5writeAttribute("homogeneous spheroid with symmetry axis z", gid, "name")
H5Gclose(gid)
# computation
gid <- H5Gopen(fid, "computation")
h5writeAttribute("EBCM, Extended Boundary Condition Method", gid, "method")
h5writeAttribute("Computation using SMARTIES, a numerically robust EBCM implementation for spheroids", gid, "description")
h5writeAttribute(software, gid, "software")
h5writeAttribute("SMARTIES", gid, "name")
H5Gclose(gid)
H5Fclose(fid)
```
<!-- export the code into standalone file -->
`r invisible(knitr::purl(xfun::with_ext(knitr::current_input(), "qmd"),output=xfun::with_ext(knitr::current_input(), "R")))`
`r invisible(source("_fun.R"))`
`r capture.output(lobstr::tree(check_h5("ar.tmat.h5")), file="tmp.out")`
_Summary of the file:_
```
`r paste(readLines("tmp.out"), collapse="\n")`
```