From c1692f67f2f6c612088447a4d72b3ab62fcb5706 Mon Sep 17 00:00:00 2001 From: Nick Wogan Date: Thu, 15 Aug 2024 17:06:28 -0700 Subject: [PATCH] futils wrapper --- clima/cython/futils.pyx | 130 +++++++++++++++++++++++++++++--- clima/cython/futils_pxd.pxd | 36 ++++++++- clima/fortran/futils.f90 | 125 ++++++++++++++++++++++++++++-- src/dependencies/CMakeLists.txt | 4 +- 4 files changed, 273 insertions(+), 22 deletions(-) diff --git a/clima/cython/futils.pyx b/clima/cython/futils.pyx index 694acbe..fd146e1 100644 --- a/clima/cython/futils.pyx +++ b/clima/cython/futils.pyx @@ -1,5 +1,26 @@ cimport futils_pxd as f_pxd +cdef void rebin_error_message(int ierr): + + cdef int err_length + cdef void *err_ptr + cdef ndarray err_ + + if ierr != 0: + f_pxd.futils_rebin_error_message_wrapper1( + &ierr, + &err_length, + &err_ptr + ) + err_ = np.zeros((),dtype=np.dtype(('S', err_length+1))) + f_pxd.futils_rebin_error_message_wrapper2( + &err_length, + err_ptr, + err_.data + ) + raise Exception(err_.item().decode()) + + cpdef rebin(ndarray[double, ndim=1] old_bins, ndarray[double, ndim=1] old_vals, ndarray[double, ndim=1] new_bins): """Rebins `old_vals` defined on `old_bins` to `new_bins`. An example is rebinning a high resolution spectra of infrared emission of Earth @@ -20,19 +41,106 @@ cpdef rebin(ndarray[double, ndim=1] old_bins, ndarray[double, ndim=1] old_vals, new_vals : ndarray[double,ndim=1] Values defined on new_bins. """ - cdef int n_old = old_bins.shape[0] - 1 - cdef int n_new = new_bins.shape[0] - 1 - cdef ndarray new_vals = np.empty(n_new, np.double) - cdef int ierr - if old_vals.shape[0] != n_old: - raise ValueError('"old_bins" and "old_vals" arguments to "rebin" have incompatible shapes') + cdef int old_bins_len = old_bins.shape[0] + cdef int old_vals_len = old_vals.shape[0] + cdef int new_bins_len = new_bins.shape[0] + cdef int new_vals_len = new_bins_len - 1 + cdef ndarray new_vals = np.empty(new_vals_len, np.double) + cdef int ierr - f_pxd.futils_rebin_wrapper(&n_old, old_bins.data, old_vals.data, - &n_new, new_bins.data, new_vals.data, &ierr) - if ierr < 0: - raise Exception("rebin returned error code: "+str(ierr)) + f_pxd.futils_rebin_wrapper( + &old_bins_len, old_bins.data, + &old_vals_len, old_vals.data, + &new_bins_len, new_bins.data, + &new_vals_len, new_vals.data, + &ierr + ) + + rebin_error_message(ierr) return new_vals - +cpdef rebin_with_errors(ndarray[double, ndim=1] old_bins, ndarray[double, ndim=1] old_vals, ndarray[double, ndim=1] old_errs, ndarray[double, ndim=1] new_bins): + """Rebins `old_vals` and `old_errs` defined on `old_bins` to `new_bins`. This function has + the same behavior as `rebin`, except it also rebins errors. + + Parameters + ---------- + old_bins : ndarray[double,ndim=1] + Edges of bins for which old_vals are defined + old_vals : ndarray[double,ndim=1] + Values defined on old_bins. + old_errs : ndarray[double,ndim=1] + Standard deviations defined on old_bins. + new_bins : ndarray[double,ndim=1] + Edges of target bin that you want to rebin to. + + Returns + ------- + new_vals : ndarray[double,ndim=1] + Values defined on new_bins. + new_errs : ndarray[double,ndim=1] + Standard deviations defined on new_bins. + """ + + cdef int old_bins_len = old_bins.shape[0] + cdef int old_vals_len = old_vals.shape[0] + cdef int old_errs_len = old_errs.shape[0] + cdef int new_bins_len = new_bins.shape[0] + cdef int new_vals_len = new_bins_len - 1 + cdef ndarray new_vals = np.empty(new_vals_len, np.double) + cdef int new_errs_len = new_bins_len - 1 + cdef ndarray new_errs = np.empty(new_errs_len, np.double) + cdef int ierr + + f_pxd.futils_rebin_with_errors_wrapper( + &old_bins_len, old_bins.data, + &old_vals_len, old_vals.data, + &old_errs_len, old_errs.data, + &new_bins_len, new_bins.data, + &new_vals_len, new_vals.data, + &new_errs_len, new_errs.data, + &ierr + ) + + rebin_error_message(ierr) + + return new_vals, new_errs + +cpdef grid_at_resolution(double wv_min, double wv_max, double R): + + cdef int wavl_len + cdef void *wavl_ptr + cdef int ierr + + f_pxd.futils_grid_at_resolution_wrapper1( + &wv_min, &wv_max, &R, + &wavl_len, &wavl_ptr, + &ierr + ) + cdef ndarray wavl = np.empty(wavl_len, np.double) + f_pxd.futils_grid_at_resolution_wrapper2( + &wavl_len, wavl_ptr, wavl.data + ) + + rebin_error_message(ierr) + + return wavl + +cpdef make_bins(ndarray[double, ndim=1] wv): + + cdef int wv_len = wv.shape[0] + cdef int wavl_len = wv_len + 1 + cdef ndarray wavl = np.empty(wavl_len, np.double) + cdef int ierr + + f_pxd.futils_make_bins_wrapper( + &wv_len, wv.data, + &wavl_len, wavl.data, + &ierr + ) + + rebin_error_message(ierr) + + return wavl diff --git a/clima/cython/futils_pxd.pxd b/clima/cython/futils_pxd.pxd index d615992..4721ea8 100644 --- a/clima/cython/futils_pxd.pxd +++ b/clima/cython/futils_pxd.pxd @@ -1,2 +1,34 @@ -cdef extern void futils_rebin_wrapper(int *n_old, double *old_bins, double *old_vals, - int *n_new, double *new_bins, double *new_vals, int *ierr) \ No newline at end of file +cdef extern void futils_rebin_wrapper( + int *old_bins_len, double *old_bins, + int *old_vals_len, double *old_vals, + int *new_bins_len, double *new_bins, + int *new_vals_len, double *new_vals, + int *ierr +) + +cdef extern void futils_rebin_with_errors_wrapper( + int *old_bins_len, double *old_bins, + int *old_vals_len, double *old_vals, + int *old_errs_len, double *old_errs, + int *new_bins_len, double *new_bins, + int *new_vals_len, double *new_vals, + int *new_errs_len, double *new_errs, + int *ierr +) + +cdef extern void futils_grid_at_resolution_wrapper1( + double *wv_min, double *wv_max, double *R, int *wavl_len, + void *wavl_ptr, int *ierr +) +cdef extern void futils_grid_at_resolution_wrapper2(int *wavl_len, void *wavl_ptr, double *wavl) + +cdef extern void futils_make_bins_wrapper( + int *wv_len, double *wv, int *wavl_len, double *wavl, int *ierr +) + +cdef extern void futils_rebin_error_message_wrapper1( + int *ierr, int *err_length, void *err_ptr +) +cdef extern void futils_rebin_error_message_wrapper2( + int *err_length, void *err_ptr, char *err +) \ No newline at end of file diff --git a/clima/fortran/futils.f90 b/clima/fortran/futils.f90 index 4121ad5..6fce7b5 100644 --- a/clima/fortran/futils.f90 +++ b/clima/fortran/futils.f90 @@ -1,13 +1,124 @@ ! futils -subroutine futils_rebin_wrapper(n_old, old_bins, old_vals, n_new, new_bins, new_vals, ierr) bind(c) +subroutine futils_rebin_wrapper(old_bins_len, old_bins, old_vals_len, old_vals, & + new_bins_len, new_bins, new_vals_len, new_vals, ierr) bind(c) use futils, only: rebin - integer(c_int), intent(in) :: n_old - real(c_double), intent(in) :: old_bins(n_old+1) - real(c_double), intent(in) :: old_vals(n_old) - integer(c_int), intent(in) :: n_new - real(c_double), intent(in) :: new_bins(n_new+1) - real(c_double), intent(out) :: new_vals(n_new) + integer(c_int), intent(in) :: old_bins_len + real(c_double), intent(in) :: old_bins(old_bins_len) + integer(c_int), intent(in) :: old_vals_len + real(c_double), intent(in) :: old_vals(old_vals_len) + integer(c_int), intent(in) :: new_bins_len + real(c_double), intent(in) :: new_bins(new_bins_len) + integer(c_int), intent(in) :: new_vals_len + real(c_double), intent(out) :: new_vals(new_vals_len) integer(c_int), intent(out) :: ierr call rebin(old_bins, old_vals, new_bins, new_vals, ierr) +end subroutine + +subroutine futils_rebin_with_errors_wrapper(old_bins_len, old_bins, old_vals_len, old_vals, old_errs_len, old_errs, & + new_bins_len, new_bins, new_vals_len, new_vals, new_errs_len, new_errs, ierr) bind(c) + use futils, only: rebin_with_errors + integer(c_int), intent(in) :: old_bins_len + real(c_double), intent(in) :: old_bins(old_bins_len) + integer(c_int), intent(in) :: old_vals_len + real(c_double), intent(in) :: old_vals(old_vals_len) + integer(c_int), intent(in) :: old_errs_len + real(c_double), intent(in) :: old_errs(old_errs_len) + integer(c_int), intent(in) :: new_bins_len + real(c_double), intent(in) :: new_bins(new_bins_len) + integer(c_int), intent(in) :: new_vals_len + real(c_double), intent(out) :: new_vals(new_vals_len) + integer(c_int), intent(in) :: new_errs_len + real(c_double), intent(out) :: new_errs(new_errs_len) + integer(c_int), intent(out) :: ierr + call rebin_with_errors(old_bins, old_vals, old_errs, new_bins, new_vals, new_errs, ierr) +end subroutine + +subroutine futils_grid_at_resolution_wrapper1(wv_min, wv_max, R, wavl_len, wavl_ptr, ierr) bind(c) + use futils, only: grid_at_resolution, dp + real(c_double), intent(in) :: wv_min + real(c_double), intent(in) :: wv_max + real(c_double), intent(in) :: R + integer(c_int), intent(out) :: wavl_len + type(c_ptr), intent(out) :: wavl_ptr + integer(c_int), intent(out) :: ierr + + real(dp), allocatable :: wavl(:) + real(c_double), pointer :: wavl_p(:) + + call grid_at_resolution(wv_min, wv_max, R, wavl, ierr) + + if (allocated(wavl)) then + wavl_len = size(wavl) + allocate(wavl_p(size(wavl))) + wavl_p = wavl + wavl_ptr = c_loc(wavl_p) + else + wavl_len = 0 + wavl_ptr = c_null_ptr + endif + +end subroutine + +subroutine futils_grid_at_resolution_wrapper2(wavl_len, wavl_ptr, wavl) bind(c) + integer(c_int), intent(in) :: wavl_len + type(c_ptr), value, intent(in) :: wavl_ptr + real(c_double), intent(out) :: wavl(wavl_len) + + real(c_double), pointer :: wavl_p(:) + + if (.not.c_associated(wavl_ptr)) return + + call c_f_pointer(wavl_ptr, wavl_p, [wavl_len]) + wavl = wavl_p + deallocate(wavl_p) + +end subroutine + +subroutine futils_make_bins_wrapper(wv_len, wv, wavl_len, wavl, ierr) bind(c) + use futils, only: make_bins + integer(c_int), intent(in) :: wv_len + real(c_double), intent(in) :: wv(wv_len) + integer(c_int), intent(in) :: wavl_len + real(c_double), intent(out) :: wavl(wavl_len) + integer(c_int), intent(out) :: ierr + call make_bins(wv, wavl, ierr) +end subroutine + +subroutine futils_rebin_error_message_wrapper1(ierr, err_length, err_ptr) bind(c) + use futils, only: rebin_error_message + integer(c_int), intent(in) :: ierr + integer(c_int), intent(out) :: err_length + type(c_ptr), intent(out) :: err_ptr + + character(:), allocatable :: err + character(kind=c_char), pointer :: err_p(:) + + err = rebin_error_message(ierr) + + if (allocated(err)) then + err_length = len(err) + allocate(err_p(err_length+1)) + call copy_string_ftoc(err, err_p) + err_ptr = c_loc(err_p) + else + err_length = 0 + err_ptr = c_null_ptr + endif + +end subroutine + +subroutine futils_rebin_error_message_wrapper2(err_length, err_ptr, err) bind(c) + integer(c_int), intent(in) :: err_length + type(c_ptr), value, intent(in) :: err_ptr + character(c_char), intent(out) :: err(err_length+1) + + character(kind=c_char), pointer :: err_p(:) + + if (.not.c_associated(err_ptr)) return + + call c_f_pointer(err_ptr, err_p, [err_length+1]) + err = err_p + deallocate(err_p) + end subroutine \ No newline at end of file diff --git a/src/dependencies/CMakeLists.txt b/src/dependencies/CMakeLists.txt index 78ccb38..65d39b8 100644 --- a/src/dependencies/CMakeLists.txt +++ b/src/dependencies/CMakeLists.txt @@ -32,9 +32,9 @@ CPMAddPackage( # futils CPMAddPackage( NAME futils - VERSION 0.1.11 + VERSION 0.1.13 GITHUB_REPOSITORY "Nicholaswogan/futils" - GIT_TAG "v0.1.11" + GIT_TAG "v0.1.13" EXCLUDE_FROM_ALL ON )