https://github.com/RadioAstronomySoftwareGroup/pyuvdata
Tip revision: c08b766220b432f1fc1f3ca21a505217c2a20423 authored by Bryna Hazelton on 13 January 2022, 01:02:29 UTC
update the changelog for v2.2.6
update the changelog for v2.2.6
Tip revision: c08b766
miriad_wrap.pyx
# -*- mode: python; coding: utf-8 -*-
# Copyright (c) 2020 Radio Astronomy Software Group
# Licensed under the 2-clause BSD License
# distutils: language = c++
# python imports
import numpy as np
# cython imports
cimport cython
cimport numpy
cimport libcpp.complex
from libc.string cimport strncmp
DEF PREAMBLE_SIZE = 5
DEF MAXVAR = 32768
# This initializes the numpy 1.7 c-api.
# cython 3.0 will do this by default.
# We may be able to just remove this then.
numpy.import_array()
cdef extern from "miriad.h":
cdef int H_BYTE "H_BYTE"
cdef int H_INT "H_INT"
cdef int H_INT2 "H_INT2 "
cdef int H_REAL "H_REAL "
cdef int H_DBLE "H_DBLE "
cdef int H_TXT "H_TXT "
cdef int H_CMPLX "H_CMPLX"
cdef int H_BYTE_SIZE "H_BYTE_SIZE"
cdef int H_INT_SIZE "H_INT_SIZE"
cdef int H_INT2_SIZE "H_INT2_SIZE"
cdef int H_INT8_SIZE "H_INT8_SIZE"
cdef int H_REAL_SIZE "H_REAL_SIZE"
cdef int H_DBLE_SIZE "H_DBLE_SIZE"
cdef int H_TXT_SIZE "H_TXT_SIZE"
cdef int H_CMPLX_SIZE "H_CMPLX_SIZE"
void uvclose_c(int tno)
void uvopen_c(int *tno, const char *name, const char *status)
void uvset_c(int tno, const char *object, const char *type, int n, double p1, double p2, double p3)
void uvrewind_c(int tno)
void uvwrite_c(int tno, double *preamble, float *data, int *flags, int n)
void uvread_c (int tno, double *preamble, float *data, int *flags, int n, int *nread);
void uvcopyvr_c (int tin, int tout)
void uvtrack_c (int tno, const char *name, const char *switches)
void uvprobvr_c (int tno, const char *var, char *type, int *length, int *updated);
void uvgetvr_c (int tno, int type, const char *var, char *data, int n);
void uvputvr_c (int tno, int type, const char *var, const char *data, int n);
void uvselect_c (int tno, const char *object, double p1, double p2, int datasel);
void haccess_c(int tno, int *ihandle, const char *keyword, const char *status, int *iostat);
void hdaccess_c(int ihandle, int *iostat);
void hwriteb_c(int item, char buf[4], int offset, int length, int *iostat)
void hreadb_c(int item, char buf[4], int offset, int length, int *iostat)
void hreadi_c(int item, int *buf, int offset, int length, int *iostat)
void hwritei_c(int item, int *buf, int offset, int length, int *iostat)
void hwritej_c(int item, short *buf, int offset, int length, int *iostat)
void hreadj_c(int item, short *buf, int offset, int length, int *iostat)
void hwritel_c(int item, long *buf, int offset, int length, int *iostat)
void hreadl_c(int item, long *buf, int offset, int length, int *iostat)
void hwriter_c(int item, float *buf, int offset, int length, int *iostat)
void hreadr_c(int item, float *buf, int offset,int length, int *iostat)
void hreadd_c(int item, double *buf, int offset, int length, int *iostat)
void hwrited_c(int item, double *buf, int offset, int length, int *iostat)
void hwritec_c(int item, float buf[2], int offset, int length, int *iostat)
void hreadc_c(int item, float buf[2], int offset, int length, int *iostat)
cdef extern from "io.h":
cdef const int ITEM_HDR_SIZE "ITEM_HDR_SIZE"
cdef char binary_item[ITEM_HDR_SIZE]
cdef char real_item[ITEM_HDR_SIZE]
cdef char int_item[ITEM_HDR_SIZE]
cdef char int2_item[ITEM_HDR_SIZE]
cdef char int8_item[ITEM_HDR_SIZE]
cdef char char_item[ITEM_HDR_SIZE]
cdef char dble_item[ITEM_HDR_SIZE]
cdef char cmplx_item[ITEM_HDR_SIZE]
cdef extern from "hio.h":
int mroundup(int a, int b)
cdef extern from "maxdimc.h":
cpdef int _MAXCHAN "MAXCHAN"
MAXCHAN = _MAXCHAN
ctypedef numpy.int_t DTYPE_t
ctypedef numpy.complex64_t DTYPE_c
ctypedef numpy.float64_t DTYPE_f64
cdef inline int GETI(int bl):
return (bl - 65536) // 2048 - 1 if bl > 65536 else (bl >> 8) - 1
cdef inline int GETJ(int bl):
return (bl - 65536) % 2048 - 1 if bl > 65536 else (bl & 255) - 1
cdef inline float MKBL(int i, int j):
return (i + 1) << 8 | (j + 1) if (i + 1 < 256 and j + 1 < 256) else (i + 1) * 2048 + (j + 1 + 65536)
cdef inline void CHK_IO(int i) except *:
if (i != 0):
raise IOError("IO failed.")
cpdef void hdaccess(int item_hdl) except +:
cdef int iostat
hdaccess_c(item_hdl, &iostat)
return
cdef int INIT(int item_hdl, char type_item[4], int size):
cdef int offset, iostat
hwriteb_c(item_hdl, type_item, 0, ITEM_HDR_SIZE, &iostat)
CHK_IO(iostat)
offset = mroundup(ITEM_HDR_SIZE, size)
return offset
cpdef int hwrite_init(int item_hdl, str type) except +:
cdef int offset
if type[0] == "a":
offset = INIT(item_hdl, char_item, H_BYTE_SIZE)
elif type[0] == "b":
offset = INIT(item_hdl, binary_item, ITEM_HDR_SIZE)
elif type[0] == "i":
offset = INIT(item_hdl, int_item, H_INT_SIZE)
elif type[0] == "j":
offset = INIT(item_hdl, int2_item, H_INT2_SIZE)
elif type[0] == "l":
offset = INIT(item_hdl, int8_item, H_INT8_SIZE)
elif type[0] == "r":
offset = INIT(item_hdl, real_item, H_REAL_SIZE)
elif type[0] == "d":
offset = INIT(item_hdl, dble_item, H_DBLE_SIZE)
elif type[0] == "c":
offset = INIT(item_hdl, cmplx_item, H_CMPLX_SIZE)
else:
raise ValueError(f"Unknown type {type[0]}")
return offset
# hread_init surmises the type of a header item from the first few bytes
cdef int FIRSTINT(char s[4]):
return (<int *>s)[0]
cpdef tuple hread_init(int item_hdl) except +*:
cdef int offset, iostat, code
cdef char s[ITEM_HDR_SIZE]
hreadb_c(item_hdl, s, 0, ITEM_HDR_SIZE, &iostat)
CHK_IO(iostat)
code = FIRSTINT(s)
if (code == FIRSTINT(char_item)):
offset = mroundup(ITEM_HDR_SIZE, H_BYTE_SIZE)
return "a", offset
elif (code == FIRSTINT(binary_item)):
return "b", ITEM_HDR_SIZE
elif (code == FIRSTINT(int_item)):
offset = mroundup(ITEM_HDR_SIZE, H_INT_SIZE)
return "i", offset
elif (code == FIRSTINT(int2_item)):
offset = mroundup(ITEM_HDR_SIZE, H_INT2_SIZE)
return "j", offset
elif (code == FIRSTINT(int8_item)):
offset = mroundup(ITEM_HDR_SIZE, H_INT8_SIZE)
return "l", offset
elif (code == FIRSTINT(real_item)):
offset = mroundup(ITEM_HDR_SIZE, H_REAL_SIZE)
return "r", offset
elif (code == FIRSTINT(dble_item)):
offset = mroundup(ITEM_HDR_SIZE, H_DBLE_SIZE)
return "d", offset
elif (code == FIRSTINT(cmplx_item)):
offset = mroundup(ITEM_HDR_SIZE, H_CMPLX_SIZE)
return "c", offset
else:
raise RuntimeError("unknown item type.")
cpdef int hwrite(int item_hdl, int offset, val, str type) except *:
cdef int iostat
cdef int int_1
cdef long lg
cdef float fl
cdef double db
cdef float cx[2]
cdef char *st
if type[0] == "a":
if not isinstance(val, (str, bytes)):
raise ValueError("expected a string")
if isinstance(val, str):
val = val.encode()
st = <char *>val
int_1 = len(val)
hwriteb_c(item_hdl, st, offset, int_1, &iostat)
CHK_IO(iostat)
offset = H_BYTE_SIZE * int_1
elif type[0] == "i":
if not isinstance(val, (int, np.int_, np.intc)):
raise ValueError("expected an int")
int_1 = <int>val
hwritei_c(item_hdl, &int_1, offset, H_INT_SIZE, &iostat)
CHK_IO(iostat)
offset = H_INT_SIZE
elif type[0] == "j":
if not isinstance(val, (int, np.int_, np.intc)):
raise ValueError("expected an int")
sh = <short>val
hwritej_c(item_hdl, &sh, offset, H_INT2_SIZE, &iostat)
CHK_IO(iostat)
offset = H_INT2_SIZE
elif type[0] == "l":
if not isinstance(val, (int, np.intc, np.int_)):
raise ValueError("expected a long")
lg = <long>val
hwritel_c(item_hdl, &lg, offset, H_INT8_SIZE, &iostat)
CHK_IO(iostat)
offset = H_INT8_SIZE
elif type[0] == "r":
if not isinstance(val, (float, np.float32)):
raise ValueError("expected a float")
fl = <float>val
hwriter_c(item_hdl, &fl, offset, H_REAL_SIZE, &iostat)
CHK_IO(iostat)
offset = H_REAL_SIZE
elif type[0] == "d":
if not isinstance(val, (float, np.float32, np.float64, np.float_)):
raise ValueError("expected a double")
db = <double>val
hwrited_c(item_hdl, &db, offset, H_DBLE_SIZE, &iostat)
CHK_IO(iostat)
offset = H_DBLE_SIZE
elif type[0] == "c":
if not isinstance(val, np.complex64):
raise ValueError("expected a complex")
cx[0] = <float>val.real
cx[1] = <float>val.imag
hwritec_c(item_hdl, cx, offset, H_CMPLX_SIZE, &iostat)
CHK_IO(iostat)
offset = H_CMPLX_SIZE
else:
raise ValueError(f"unknown item type: {type [0]}")
return offset
cpdef tuple hread(int item_hdl, int offset, str type) except +:
cdef int iostat
cdef int int_1
cdef short sh
cdef long lg
cdef float fl
cdef double db
cdef float cx[2]
cdef char st[2]
if type[0] == "a":
hreadb_c(item_hdl, st, offset, H_BYTE_SIZE, &iostat)
CHK_IO(iostat)
return st, H_BYTE_SIZE
elif type[0] == "i":
hreadi_c(item_hdl, &int_1, offset, H_INT_SIZE, &iostat)
CHK_IO(iostat)
return int_1, H_INT_SIZE
elif type[0] == "j":
hreadj_c(item_hdl, &sh, offset, H_INT2_SIZE, &iostat)
CHK_IO(iostat)
return sh, H_INT2_SIZE
elif type[0] == "l":
hreadl_c(item_hdl, &lg, offset, H_INT8_SIZE, &iostat)
CHK_IO(iostat)
return lg, H_INT8_SIZE
elif type[0] == "r":
hreadr_c(item_hdl, &fl, offset, H_REAL_SIZE, &iostat)
CHK_IO(iostat)
return fl, H_REAL_SIZE
elif type[0] == "d":
hreadd_c(item_hdl, &db, offset, H_DBLE_SIZE, &iostat)
CHK_IO(iostat)
return db, H_DBLE_SIZE
elif type[0] == "c":
hreadc_c(item_hdl, cx, offset, H_CMPLX_SIZE, &iostat)
CHK_IO(iostat)
# check this line here
return <float>cx[0] + 1j * <float>cx[1], H_CMPLX_SIZE
else:
raise ValueError(f"unknown item type: {type[0]}")
cdef class UV:
cpdef int tno
cpdef long decimate
cpdef long decphase
cpdef long intcnt
cpdef double curtime
def __init__(self, str filename, str status, str corrmode):
self.tno = -1
self.decimate = 1
self.decphase = 0
self.intcnt = -1
self.curtime = -1
if corrmode[0] not in ["r", "j"]:
raise ValueError(f"UV corrmode must be 'r' or 'j' but received {corrmode}.")
try:
uvopen_c(&self.tno, filename.encode(), status.encode())
uvset_c(self.tno, "preamble", "uvw/time/baseline", 0, 0., 0., 0.)
uvset_c(self.tno, "corr", corrmode.encode(), 0, 0., 0., 0.)
except:
self.tno = -1
raise
cpdef void close(self):
if self.tno != -1:
uvclose_c(self.tno)
self.tno = -1
return
def __dealloc__(self):
self.close()
return
@cython.boundscheck(False)
cdef _get_j_type(self, int htype, char *name, int length):
cdef numpy.ndarray[dtype=int, ndim=1] arr = np.zeros((length,), dtype=np.int16)
uvgetvr_c(self.tno, htype, name, <char *>&arr[0], length)
if length == 1:
return arr.item(0)
return arr
@cython.boundscheck(False)
cdef _get_i_type(self, int htype, char *name, int length):
cdef numpy.ndarray[dtype=numpy.int32_t, ndim=1] arr = np.zeros((length,), dtype=np.int32)
uvgetvr_c(self.tno, htype, name, <char *>&arr[0], length)
if length == 1:
return arr.item(0)
return arr
@cython.boundscheck(False)
cdef _get_d_type(self, int htype, char *name, int length):
cdef numpy.ndarray[dtype=DTYPE_f64, ndim=1] arr = np.zeros((length,), dtype=np.float64)
uvgetvr_c(self.tno, htype, name, <char *>&arr[0], length)
if length == 1:
return arr.item(0)
return arr
@cython.boundscheck(False)
cdef _get_r_type(self, int htype, char *name, int length):
cdef numpy.ndarray[dtype=numpy.float32_t, ndim=1] arr = np.zeros((length,), dtype=np.float32)
uvgetvr_c(self.tno, htype, name, <char *>&arr[0], length)
if length == 1:
return arr.item(0)
return arr
@cython.boundscheck(False)
cdef _get_c_type(self, int htype, char *name, int length):
cdef numpy.ndarray[dtype=DTYPE_c, ndim=1] arr = np.zeros((length,), dtype=np.complex64)
uvgetvr_c(self.tno, htype, name, <char *>&arr[0], length)
if length == 1:
return arr.item(0)
return arr
@cython.boundscheck(False)
cdef void _store_j_type(self, int htype, char *name, numpy.ndarray[dtype=int] value):
uvputvr_c(self.tno, htype, name, <char *>&value[0], value.size)
return
@cython.boundscheck(False)
cdef void _store_i_type(self, int htype, char *name, numpy.ndarray[dtype=numpy.int32_t] value):
uvputvr_c(self.tno, htype, name, <char *>&value[0], value.size)
return
@cython.boundscheck(False)
cdef void _store_d_type(self, int htype, char *name, numpy.ndarray[dtype=DTYPE_f64] value):
uvputvr_c(self.tno, htype, name, <char *>&value[0], value.size)
return
@cython.boundscheck(False)
cdef void _store_r_type(self, int htype, char *name, numpy.ndarray[dtype=numpy.float32_t] value):
uvputvr_c(self.tno, htype, name, <char *>&value[0], value.size)
return
@cython.boundscheck(False)
cdef void _store_c_type(self, int htype, char *name, numpy.ndarray[dtype=DTYPE_c] value):
uvputvr_c(self.tno, htype, name, <char *>&value[0], value.size)
return
cpdef void rewind(self):
uvrewind_c(self.tno)
self.intcnt = -1
self.curtime = -1
return
@cython.boundscheck(False)
cpdef raw_read(self, int n2read) except +:
cdef int nread, i, j
cdef double preamble[PREAMBLE_SIZE]
cdef numpy.ndarray[numpy.complex64_t , ndim=1] data = np.zeros((n2read,), dtype=np.complex64)
cdef numpy.ndarray[int, ndim=1] flags = np.zeros((n2read,), dtype=np.intc)
cdef numpy.ndarray[DTYPE_f64, ndim=1] uvw = np.zeros((3,), dtype=np.float64)
while True:
uvread_c(self.tno, preamble, <float *>&data[0], <int *>&flags[0], n2read, &nread)
if (preamble[3] != self.curtime):
self.intcnt += 1
self.curtime = preamble[3]
if ((self.intcnt - self.decphase) % self.decimate == 0 or nread == 0):
break
uvw[0] = preamble[0]
uvw[1] = preamble[1]
uvw[2] = preamble[2]
i = GETI(<int>preamble[4])
j = GETJ(<int>preamble[4])
return (uvw, preamble[3], (i, j)), data, flags, nread
cpdef void raw_write(self, object input_preamble, numpy.ndarray[dtype=DTYPE_c, ndim=1] data, numpy.ndarray[dtype=int, ndim=1] flags) except +:
cdef int nread
cdef double preamble[PREAMBLE_SIZE]
cdef double t = input_preamble[1]
cdef int i = input_preamble[2][0], j = input_preamble[2][1]
if len(input_preamble[0]) != 3:
raise ValueError(f"uvw must have shape (3,) but got {len(input_preamble[0])}")
preamble[0] = input_preamble[0][0]
preamble[1] = input_preamble[0][1]
preamble[2] = input_preamble[0][2]
preamble[3] = t
preamble[4] = MKBL(i, j)
uvwrite_c(self.tno, preamble, <float *>&data[0], <int *>&flags[0], data.size)
return
cpdef void copyvr(self, UV uv):
uvcopyvr_c(uv.tno, self.tno)
return
cpdef void trackvr(self, str name, str switches):
uvtrack_c(self.tno, name.encode(), switches.encode())
return
cpdef _rdvr(self, str name, str type) except +*:
cdef char value[MAXVAR]
cdef int length, updated, elem_size
uvprobvr_c(self.tno, name.encode(), value, &length, &updated)
if type[0] == "a":
elem_size = 1
elif type[0] == "j":
elem_size = 2
elif type[0] in ["i", "r"]:
elem_size = 4
elif type[0] in ["d", "c"]:
elem_size = 8
else:
raise ValueError(f"unknown type of UV variable {name}: {type[0]}")
if length * elem_size > MAXVAR:
raise ValueError(f"UV variable {name} is too big for pyuvdata's internal buffers")
if type[0] == "a":
uvgetvr_c(self.tno, H_BYTE, name.encode(), value, length + 1)
return value.decode("utf-8")
elif type[0] == "j":
return self._get_j_type(H_INT2, name.encode(), length)
elif type[0] == "i":
return self._get_i_type(H_INT, name.encode(), length)
elif type[0] == "r":
return self._get_r_type(H_REAL, name.encode(), length)
elif type[0] == "d":
return self._get_d_type(H_DBLE, name.encode(), length)
elif type[0] == "c":
return self._get_c_type(H_CMPLX, name.encode(), length)
return
cpdef _wrvr(self, str name, str type, value) except +*:
cdef char c_value[MAXVAR]
cdef char *st
if isinstance(value, np.ndarray):
if value.ndim != 1:
raise ValueError(f"Input variable {name} has rank {value.ndim} but must be rank 1.")
if type[0] == "a":
if not isinstance(value, (str, bytes)):
raise ValueError("Expected a string")
if isinstance(value, str):
value = value.encode()
st = <char *>value
uvputvr_c(self.tno, H_BYTE, name.encode(), st, len(value) + 1)
elif type[0] == "j":
return self._store_j_type(H_INT2, name.encode(), np.atleast_1d(value).astype(np.int16))
elif type[0] == "i":
return self._store_i_type(H_INT, name.encode(), np.atleast_1d(value).astype(np.int32))
elif type[0] == "r":
return self._store_r_type(H_REAL, name.encode(), np.atleast_1d(value).astype(np.float32))
elif type[0] == "d":
return self._store_d_type(H_DBLE, name.encode(), np.atleast_1d(value).astype(np.float64))
elif type[0] == "c":
return self._store_c_type(H_CMPLX, name.encode(), np.atleast_1d(value).astype(np.complex64))
else:
raise ValueError(f"Unkown UV variable type {type[0]}")
return
cpdef void _select(self, str name, numpy.float64_t ind1, numpy.float64_t ind2, int include_flag) except +:
# we used to only call strncmp(name, decimation, 5) so only look at first 5 letters
if strncmp(name.encode(), "decimation", 5) == 0:
self.decimate = <long> ind1
self.decphase = <long> ind2
else:
uvselect_c(self.tno, name.encode(), <double>ind1, <double>ind2, include_flag)
return
cpdef int haccess(self, str name, str mode) except *:
cdef int item_hdl, iostat
haccess_c(self.tno, &item_hdl, name.encode(), mode.encode(), &iostat)
CHK_IO(iostat)
return item_hdl