diff --git a/doc/specs/stdlib_io.md b/doc/specs/stdlib_io.md index 40cb2b426..c33511c5e 100644 --- a/doc/specs/stdlib_io.md +++ b/doc/specs/stdlib_io.md @@ -36,7 +36,7 @@ program demo_loadtxt use stdlib_io, only: loadtxt implicit none real, allocatable :: x(:,:) - call loadtxt('example.dat', x) + call loadtxt('example.dat', x) end program demo_loadtxt ``` @@ -128,6 +128,98 @@ program demo_savetxt use stdlib_io, only: savetxt implicit none real :: x(3,2) = 1 - call savetxt('example.dat', x) + call savetxt('example.dat', x) end program demo_savetxt ``` + + +## `load_npy` + +### Status + +Experimental + +### Description + +Loads an `array` from a npy formatted binary file. + +### Syntax + +`call [[stdlib_io_npy(module):load_npy(interface)]](filename, array[, iostat][, iomsg])` + +### Arguments + +`filename`: Shall be a character expression containing the file name from which to load the `array`. + This argument is `intent(in)`. + +`array`: Shall be an allocatable array of any rank of type `real`, `complex` or `integer`. + This argument is `intent(out)`. + +`iostat`: Default integer, contains status of loading to file, zero in case of success. + It is an optional argument, in case not present the program will halt for non-zero status. + This argument is `intent(out)`. + +`iomsg`: Deferred length character value, contains error message in case `iostat` is non-zero. + It is an optional argument, error message will be dropped if not present. + This argument is `intent(out)`. + +### Return value + +Returns an allocated `array` with the content of `filename` in case of success. + +### Example + +```fortran +program demo_loadnpy + use stdlib_io_npy, only: load_npy + implicit none + real, allocatable :: x(:,:) + call loadtxt('example.npy', x) +end program demo_loadnpy +``` + + +## `save_npy` + +### Status + +Experimental + +### Description + +Saves an `array` into a npy formatted binary file. + +### Syntax + +`call [[stdlib_io_npy(module):save_npy(interface)]](filename, array[, iostat][, iomsg])` + +### Arguments + +`filename`: Shall be a character expression containing the name of the file that will contain the `array`. + This argument is `intent(in)`. + +`array`: Shall be an array of any rank of type `real`, `complex` or `integer`. + This argument is `intent(in)`. + +`iostat`: Default integer, contains status of saving to file, zero in case of success. + It is an optional argument, in case not present the program will halt for non-zero status. + This argument is `intent(out)`. + +`iomsg`: Deferred length character value, contains error message in case `iostat` is non-zero. + It is an optional argument, error message will be dropped if not present. + This argument is `intent(out)`. + +### Output + +Provides a npy file called `filename` that contains the rank-2 `array`. + +### Example + +```fortran +program demo_savenpy + use stdlib_io_npy, only: save_npy + implicit none + real :: x(3,2) = 1 + call save_npy('example.npy', x) +end program demo_savenpy +``` diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index bcb4931b3..d3e177a2f 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -7,6 +7,9 @@ set(fppFiles stdlib_bitsets_64.fypp stdlib_bitsets_large.fypp stdlib_io.fypp + stdlib_io_npy.fypp + stdlib_io_npy_load.fypp + stdlib_io_npy_save.fypp stdlib_kinds.fypp stdlib_linalg.fypp stdlib_linalg_diag.fypp diff --git a/src/Makefile.manual b/src/Makefile.manual index f573fa6cf..78e1508e5 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -4,6 +4,9 @@ SRCFYPP = \ stdlib_bitsets_large.fypp \ stdlib_bitsets.fypp \ stdlib_io.fypp \ + stdlib_io_npy.fypp \ + stdlib_io_npy_load.fypp \ + stdlib_io_npy_save.fypp \ stdlib_kinds.fypp \ stdlib_linalg.fypp \ stdlib_linalg_diag.fypp \ @@ -87,6 +90,16 @@ stdlib_io.o: \ stdlib_optval.o \ stdlib_kinds.o \ stdlib_ascii.o +stdlib_io_npy.o: \ + stdlib_kinds.o +stdlib_io_npy_load.o: \ + stdlib_io_npy.o \ + stdlib_error.o \ + stdlib_strings.o +stdlib_io_npy_save.o: \ + stdlib_io_npy.o \ + stdlib_error.o \ + stdlib_strings.o stdlib_linalg.o: \ stdlib_kinds.o stdlib_linalg_diag.o: \ diff --git a/src/stdlib_io_npy.fypp b/src/stdlib_io_npy.fypp new file mode 100644 index 000000000..bf69a6a0c --- /dev/null +++ b/src/stdlib_io_npy.fypp @@ -0,0 +1,126 @@ +! SPDX-Identifer: MIT + +#:include "common.fypp" +#:set RANKS = range(1, MAXRANK + 1) +#:set KINDS_TYPES = REAL_KINDS_TYPES + INT_KINDS_TYPES + CMPLX_KINDS_TYPES + +!> Description of the npy format taken from +!> https://numpy.org/doc/stable/reference/generated/numpy.lib.format.html +!> +!>## Format Version 1.0 +!> +!> 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 +!> of len(magic string) + 2 + len(length) + HEADER_LEN be evenly divisible by 64 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, the dictionary keys are sorted in alphabetic order. +!> This is for convenience only. A writer SHOULD implement this if possible. A reader MUST +!> NOT depend on this. +!> +!> Following the header comes the array data. If the dtype contains Python objects +!> (i.e. dtype.hasobject is True), then the data is a Python pickle of the array. +!> Otherwise the data is the contiguous (either C- or Fortran-, depending on fortran_order) +!> bytes of the array. Consumers can figure out the number of bytes by multiplying the +!> number of elements given by the shape (noting that shape=() means there is 1 element) +!> by dtype.itemsize. +!> +!>## Format Version 2.0 +!> +!> The version 1.0 format only allowed the array header to have a total size of 65535 bytes. +!> This can be exceeded by structured arrays with a large number of columns. +!> The version 2.0 format extends the header size to 4 GiB. numpy.save will automatically +!> save in 2.0 format if the data requires it, else it will always use the more compatible +!> 1.0 format. +!> +!> The description of the fourth element of the header therefore has become: +!> “The next 4 bytes form a little-endian unsigned int: the length of the header data +!> HEADER_LEN.” +!> +!>## Format Version 3.0 +!> +!> This version replaces the ASCII string (which in practice was latin1) with a +!> utf8-encoded string, so supports structured types with any unicode field names. +module stdlib_io_npy + use stdlib_kinds, only : int8, int16, int32, int64, sp, dp, xdp, qp + implicit none + private + + public :: save_npy, load_npy + + + !> Version: experimental + !> + !> Save multidimensional array in npy format + !> ([Specification](../page/specs/stdlib_io.html#save_npy)) + interface save_npy + #:for k1, t1 in KINDS_TYPES + #:for rank in RANKS + module subroutine save_npy_${t1[0]}$${k1}$_${rank}$(filename, array, iostat, iomsg) + character(len=*), intent(in) :: filename + ${t1}$, intent(in) :: array${ranksuffix(rank)}$ + integer, intent(out), optional :: iostat + character(len=:), allocatable, intent(out), optional :: iomsg + end subroutine save_npy_${t1[0]}$${k1}$_${rank}$ + #:endfor + #:endfor + end interface save_npy + + !> Version: experimental + !> + !> Load multidimensional array in npy format + !> ([Specification](../page/specs/stdlib_io.html#load_npy)) + interface load_npy + #:for k1, t1 in KINDS_TYPES + #:for rank in RANKS + module subroutine load_npy_${t1[0]}$${k1}$_${rank}$(filename, array, iostat, iomsg) + character(len=*), intent(in) :: filename + ${t1}$, allocatable, intent(out) :: array${ranksuffix(rank)}$ + integer, intent(out), optional :: iostat + character(len=:), allocatable, intent(out), optional :: iomsg + end subroutine load_npy_${t1[0]}$${k1}$_${rank}$ + #:endfor + #:endfor + end interface load_npy + + + character(len=*), parameter :: nl = achar(10) + + character(len=*), parameter :: & + type_iint8 = " Implementation of loading npy files into multidimensional arrays +submodule (stdlib_io_npy) stdlib_io_npy_load + use stdlib_error, only : error_stop + use stdlib_strings, only : to_string, starts_with + implicit none + +contains + +#:for k1, t1 in KINDS_TYPES + #:for rank in RANKS + !> Load a ${rank}$-dimensional array from a npy file + module subroutine load_npy_${t1[0]}$${k1}$_${rank}$(filename, array, iostat, iomsg) + !> Name of the npy file to load from + character(len=*), intent(in) :: filename + !> Array to be loaded from the npy file + ${t1}$, allocatable, intent(out) :: array${ranksuffix(rank)}$ + !> Error status of loading, zero on success + integer, intent(out), optional :: iostat + !> Associated error message in case of non-zero status code + character(len=:), allocatable, intent(out), optional :: iomsg + + character(len=*), parameter :: vtype = type_${t1[0]}$${k1}$ + integer, parameter :: rank = ${rank}$ + + integer :: io, stat + character(len=:), allocatable :: msg + + open(newunit=io, file=filename, form="unformatted", access="stream", iostat=stat) + catch: block + character(len=:), allocatable :: this_type + integer, allocatable :: vshape(:) + + call get_descriptor(io, filename, this_type, vshape, stat, msg) + if (stat /= 0) exit catch + + if (this_type /= vtype) then + stat = 1 + msg = "File '"//filename//"' contains data of type '"//this_type//"', "//& + & "but expected '"//vtype//"'" + exit catch + end if + + if (size(vshape) /= rank) then + stat = 1 + msg = "File '"//filename//"' contains data of rank "//& + & to_string(size(vshape))//", but expected "//& + & to_string(rank) + exit catch + end if + + call allocator(array, vshape, stat) + if (stat /= 0) then + msg = "Failed to allocate array of type '"//vtype//"' "//& + & "with total size of "//to_string(product(vshape)) + exit catch + end if + + read(io, iostat=stat) array + end block catch + close(io) + + if (present(iostat)) then + iostat = stat + else if (stat /= 0) then + if (allocated(msg)) then + call error_stop("Failed to read array from file '"//filename//"'"//nl//& + & msg) + else + call error_stop("Failed to read array from file '"//filename//"'") + end if + end if + + if (present(iomsg).and.allocated(msg)) call move_alloc(msg, iomsg) + contains + + !> Wrapped intrinsic allocate to create an allocation from a shape array + subroutine allocator(array, vshape, stat) + !> Instance of the array to be allocated + ${t1}$, allocatable, intent(out) :: array${ranksuffix(rank)}$ + !> Dimensions to allocate for + integer, intent(in) :: vshape(:) + !> Status of allocate + integer, intent(out) :: stat + + allocate(array( & + #:for i in range(rank-1) + & vshape(${i+1}$), & + #:endfor + & vshape(${rank}$)), & + & stat=stat) + + end subroutine allocator + + end subroutine load_npy_${t1[0]}$${k1}$_${rank}$ + #:endfor +#:endfor + + + !> Read the npy header from a binary file and retrieve the descriptor string. + subroutine get_descriptor(io, filename, vtype, vshape, stat, msg) + !> Unformatted, stream accessed unit + integer, intent(in) :: io + !> Filename for error reporting + character(len=*), intent(in) :: filename + !> Type of data saved in npy file + character(len=:), allocatable, intent(out) :: vtype + !> Shape descriptor of the + integer, allocatable, intent(out) :: vshape(:) + !> Status of operation + integer, intent(out) :: stat + !> Associated error message in case of non-zero status + character(len=:), allocatable, intent(out) :: msg + + integer :: major, header_len, i + character(len=:), allocatable :: dict + character(len=8) :: header + character :: buf(4) + logical :: fortran_order + + read(io, iostat=stat) header + if (stat /= 0) return + + call parse_header(header, major, stat, msg) + if (stat /= 0) return + + read(io, iostat=stat) buf(1:merge(4, 2, major > 1)) + if (stat /= 0) return + + if (major > 1) then + header_len = ichar(buf(1)) & + & + ichar(buf(2)) * 2**8 & + & + ichar(buf(3)) * 2**16 & + & + ichar(buf(4)) * 2**32 + else + header_len = ichar(buf(1)) & + & + ichar(buf(2)) * 2**8 + end if + allocate(character(header_len) :: dict, stat=stat) + if (stat /= 0) return + + read(io, iostat=stat) dict + if (stat /= 0) return + + if (dict(header_len:header_len) /= nl) then + stat = 1 + msg = "Descriptor length does not match" + return + end if + + if (scan(dict, achar(0)) > 0) then + stat = 1 + msg = "Nul byte not allowed in descriptor string" + return + end if + + call parse_descriptor(trim(dict(:len(dict)-1)), filename, & + & vtype, fortran_order, vshape, stat, msg) + if (stat /= 0) return + + if (.not.fortran_order) then + vshape = [(vshape(i), i = size(vshape), 1, -1)] + end if + end subroutine get_descriptor + + + !> Parse the first eight bytes of the npy header to verify the data + subroutine parse_header(header, major, stat, msg) + !> Header of the binary file + character(len=*), intent(in) :: header + !> Major version of the npy format + integer, intent(out) :: major + !> Status of operation + integer, intent(out) :: stat + !> Associated error message in case of non-zero status + character(len=:), allocatable, intent(out) :: msg + + integer :: minor + + if (header(1:1) /= magic_number) then + stat = 1 + msg = "Expected z'93' but got z'"//to_string(ichar(header(1:1)))//"' "//& + & "as first byte" + return + end if + + if (header(2:6) /= magic_string) then + stat = 1 + msg = "Expected identifier '"//magic_string//"'" + return + end if + + major = ichar(header(7:7)) + if (.not.any(major == [1, 2, 3])) then + stat = 1 + msg = "Unsupported format major version number '"//to_string(major)//"'" + return + end if + + minor = ichar(header(8:8)) + if (minor /= 0) then + stat = 1 + msg = "Unsupported format version "// & + & "'"//to_string(major)//"."//to_string(minor)//"'" + return + end if + end subroutine parse_header + + !> Parse the descriptor in the npy header. This routine implements a minimal + !> non-recursive parser for serialized Python dictionaries. + subroutine parse_descriptor(input, filename, vtype, fortran_order, vshape, stat, msg) + !> Input string to parse as descriptor + character(len=*), intent(in) :: input + !> Filename for error reporting + character(len=*), intent(in) :: filename + !> Type of the data stored, retrieved from field `descr` + character(len=:), allocatable, intent(out) :: vtype + !> Whether the data is in left layout, retrieved from field `fortran_order` + logical, intent(out) :: fortran_order + !> Shape of the stored data, retrieved from field `shape` + integer, allocatable, intent(out) :: vshape(:) + !> Status of operation + integer, intent(out) :: stat + !> Associated error message in case of non-zero status + character(len=:), allocatable, intent(out) :: msg + + enum, bind(c) + enumerator :: invalid, string, lbrace, rbrace, comma, colon, & + lparen, rparen, bool, literal, space + end enum + + type :: token_type + integer :: first, last, kind + end type token_type + + integer :: pos + character(len=:), allocatable :: key + type(token_type) :: token, last + logical :: has_descr, has_shape, has_fortran_order + + has_descr = .false. + has_shape = .false. + has_fortran_order = .false. + pos = 0 + call next_token(input, pos, token, [lbrace], stat, msg) + if (stat /= 0) return + + last = token_type(pos, pos, comma) + do while (pos < len(input)) + call get_token(input, pos, token) + select case(token%kind) + case(space) + continue + case(comma) + if (token%kind == last%kind) then + stat = 1 + msg = make_message(filename, input, token%first, token%last, & + & "Comma cannot appear at this point") + return + end if + last = token + case(rbrace) + exit + case(string) + if (token%kind == last%kind) then + stat = 1 + msg = make_message(filename, input, token%first, token%last, & + & "String cannot appear at this point") + return + end if + last = token + + key = input(token%first+1:token%last-1) + call next_token(input, pos, token, [colon], stat, msg) + if (stat /= 0) return + + if (key == "descr" .and. has_descr & + & .or. key == "fortran_order" .and. has_fortran_order & + & .or. key == "shape" .and. has_shape) then + stat = 1 + msg = make_message(filename, input, last%first, last%last, & + & "Duplicate entry for '"//key//"' found") + return + end if + + select case(key) + case("descr") + call next_token(input, pos, token, [string], stat, msg) + if (stat /= 0) return + + vtype = input(token%first+1:token%last-1) + has_descr = .true. + + case("fortran_order") + call next_token(input, pos, token, [bool], stat, msg) + if (stat /= 0) return + + fortran_order = input(token%first:token%last) == "True" + has_fortran_order = .true. + + case("shape") + call parse_tuple(input, pos, vshape, stat, msg) + + has_shape = .true. + + case default + stat = 1 + msg = make_message(filename, input, last%first, last%last, & + & "Invalid entry '"//key//"' in dictionary encountered") + return + end select + case default + stat = 1 + msg = make_message(filename, input, token%first, token%last, & + & "Invalid token encountered") + return + end select + end do + + if (.not.has_descr) then + stat = 1 + msg = make_message(filename, input, 1, pos, & + & "Dictionary does not contain required entry 'descr'") + end if + + if (.not.has_shape) then + stat = 1 + msg = make_message(filename, input, 1, pos, & + & "Dictionary does not contain required entry 'shape'") + end if + + if (.not.has_fortran_order) then + stat = 1 + msg = make_message(filename, input, 1, pos, & + & "Dictionary does not contain required entry 'fortran_order'") + end if + + contains + + function make_message(filename, input, first, last, message) result(str) + !> Filename for context + character(len=*), intent(in) :: filename + !> Input string to parse + character(len=*), intent(in) :: input + !> Offset in the input + integer, intent(in) :: first, last + !> Error message + character(len=*), intent(in) :: message + !> Final output message + character(len=:), allocatable :: str + + character(len=*), parameter :: nl = new_line('a') + + str = message // nl // & + & " --> " // filename // ":1:" // to_string(first) // "-" // to_string(last) // nl // & + & " |" // nl // & + & "1 | " // input // nl // & + & " |" // repeat(" ", first) // repeat("^", last - first + 1) // nl // & + & " |" + end function make_message + + !> Parse a tuple of integers into an array of integers + subroutine parse_tuple(input, pos, tuple, stat, msg) + !> Input string to parse + character(len=*), intent(in) :: input + !> Offset in the input, will be advanced after reading + integer, intent(inout) :: pos + !> Array representing tuple of integers + integer, allocatable, intent(out) :: tuple(:) + !> Status of operation + integer, intent(out) :: stat + !> Associated error message in case of non-zero status + character(len=:), allocatable, intent(out) :: msg + + type(token_type) :: token + integer :: last, itmp + + allocate(tuple(0), stat=stat) + if (stat /= 0) return + + call next_token(input, pos, token, [lparen], stat, msg) + if (stat /= 0) return + + last = comma + do while (pos < len(input)) + call get_token(input, pos, token) + select case(token%kind) + case(space) + continue + case(literal) + if (token%kind == last) then + stat = 1 + msg = make_message(filename, input, token%first, token%last, & + & "Invalid token encountered") + return + end if + last = token%kind + read(input(token%first:token%last), *, iostat=stat) itmp + if (stat /= 0) then + return + end if + tuple = [tuple, itmp] + case(comma) + if (token%kind == last) then + stat = 1 + msg = make_message(filename, input, token%first, token%last, & + & "Invalid token encountered") + return + end if + last = token%kind + case(rparen) + exit + case default + stat = 1 + msg = make_message(filename, input, token%first, token%last, & + & "Invalid token encountered") + return + end select + end do + end subroutine parse_tuple + + !> Get the next allowed token + subroutine next_token(input, pos, token, allowed_token, stat, msg) + !> Input string to parse + character(len=*), intent(in) :: input + !> Current offset in the input string + integer, intent(inout) :: pos + !> Last token parsed + type(token_type), intent(out) :: token + !> Tokens allowed in the current context + integer, intent(in) :: allowed_token(:) + !> Status of operation + integer, intent(out) :: stat + !> Associated error message in case of non-zero status + character(len=:), allocatable, intent(out) :: msg + + stat = pos + do while (pos < len(input)) + call get_token(input, pos, token) + if (token%kind == space) then + continue + else if (any(token%kind == allowed_token)) then + stat = 0 + exit + else + stat = 1 + msg = make_message(filename, input, token%first, token%last, & + & "Invalid token encountered") + exit + end if + end do + end subroutine next_token + + !> Tokenize input string + subroutine get_token(input, pos, token) + !> Input strin to tokenize + character(len=*), intent(in) :: input + !> Offset in input string, will be advanced + integer, intent(inout) :: pos + !> Returned token from the next position + type(token_type), intent(out) :: token + + character :: quote + + pos = pos + 1 + select case(input(pos:pos)) + case("""", "'") + quote = input(pos:pos) + token%first = pos + pos = pos + 1 + do while (pos <= len(input)) + if (input(pos:pos) == quote) then + token%last = pos + exit + else + pos = pos + 1 + end if + end do + token%kind = string + case("0", "1", "2", "3", "4", "5", "6", "7", "8", "9") + token%first = pos + do while (pos <= len(input)) + if (.not.any(input(pos:pos) == ["0", "1", "2", "3", "4", "5", "6", "7", "8", "9"])) then + pos = pos - 1 + token%last = pos + exit + else + pos = pos + 1 + end if + end do + token%kind = literal + case("T") + if (starts_with(input(pos:), "True")) then + token = token_type(pos, pos+3, bool) + pos = pos + 3 + else + token = token_type(pos, pos, invalid) + end if + case("F") + if (starts_with(input(pos:), "False")) then + token = token_type(pos, pos+4, bool) + pos = pos + 4 + else + token = token_type(pos, pos, invalid) + end if + case("{") + token = token_type(pos, pos, lbrace) + case("}") + token = token_type(pos, pos, rbrace) + case(",") + token = token_type(pos, pos, comma) + case(":") + token = token_type(pos, pos, colon) + case("(") + token = token_type(pos, pos, lparen) + case(")") + token = token_type(pos, pos, rparen) + case(" ", nl) + token = token_type(pos, pos, space) + case default + token = token_type(pos, pos, invalid) + end select + + end subroutine get_token + + end subroutine parse_descriptor + +end submodule stdlib_io_npy_load diff --git a/src/stdlib_io_npy_save.fypp b/src/stdlib_io_npy_save.fypp new file mode 100644 index 000000000..4dccadeb6 --- /dev/null +++ b/src/stdlib_io_npy_save.fypp @@ -0,0 +1,137 @@ +! SPDX-Identifer: MIT + +#:include "common.fypp" +#:set RANKS = range(1, MAXRANK + 1) +#:set KINDS_TYPES = REAL_KINDS_TYPES + INT_KINDS_TYPES + CMPLX_KINDS_TYPES + +!> Implementation of saving multidimensional arrays to npy files +submodule (stdlib_io_npy) stdlib_io_npy_save + use stdlib_error, only : error_stop + use stdlib_strings, only : to_string + implicit none + +contains + + + !> Generate magic header string for npy format + pure function magic_header(major, minor) result(str) + !> Major version of npy format + integer, intent(in) :: major + !> Minor version of npy format + integer, intent(in) :: minor + !> Magic string for npy format + character(len=8) :: str + + str = magic_number // magic_string // achar(major) // achar(minor) + end function magic_header + + + !> Generate header for npy format + pure function npy_header(vtype, vshape) result(str) + !> Type of variable + character(len=*), intent(in) :: vtype + !> Shape of variable + integer, intent(in) :: vshape(:) + !> Header string for npy format + character(len=:), allocatable :: str + + integer, parameter :: len_v10 = 8 + 2, len_v20 = 8 + 4, block_size = 64 + + str = & + "{'descr': '"//vtype//& + "', 'fortran_order': True, 'shape': "//& + shape_str(vshape)//", }" + + if (len(str) + len_v10 >= 65535) then + str = str // & + & repeat(" ", block_size - mod(len(str) + len_v20 + 1, block_size)) // nl + str = magic_header(2, 0) // to_bytes_i4(int(len(str))) // str + else + str = str // & + & repeat(" ", block_size - mod(len(str) + len_v10 + 1, block_size)) // nl + str = magic_header(1, 0) // to_bytes_i2(int(len(str))) // str + end if + end function npy_header + + !> Write integer as byte string in little endian encoding + pure function to_bytes_i4(val) result(str) + !> Integer value to convert to bytes + integer, intent(in) :: val + !> String of bytes + character(len=4) :: str + + str = achar(mod(val, 2**8)) // & + & achar(mod(val, 2**16) / 2**8) // & + & achar(mod(val, 2**32) / 2**16) // & + & achar(val / 2**32) + end function to_bytes_i4 + + + !> Write integer as byte string in little endian encoding, 2-byte truncated version + pure function to_bytes_i2(val) result(str) + !> Integer value to convert to bytes + integer, intent(in) :: val + !> String of bytes + character(len=2) :: str + + str = achar(mod(val, 2**8)) // & + & achar(mod(val, 2**16) / 2**8) + end function to_bytes_i2 + + + !> Print array shape as tuple of int + pure function shape_str(vshape) result(str) + !> Shape of variable + integer, intent(in) :: vshape(:) + !> Shape string for npy format + character(len=:), allocatable :: str + + integer :: i + + str = "(" + do i = 1, size(vshape) + str = str//to_string(vshape(i))//", " + enddo + str = str//")" + end function shape_str + + +#:for k1, t1 in KINDS_TYPES + #:for rank in RANKS + !> Save ${rank}$-dimensional array in npy format + module subroutine save_npy_${t1[0]}$${k1}$_${rank}$(filename, array, iostat, iomsg) + !> Name of the npy file to load from + character(len=*), intent(in) :: filename + !> Array to be loaded from the npy file + ${t1}$, intent(in) :: array${ranksuffix(rank)}$ + !> Error status of loading, zero on success + integer, intent(out), optional :: iostat + !> Associated error message in case of non-zero status code + character(len=:), allocatable, intent(out), optional :: iomsg + + character(len=*), parameter :: vtype = type_${t1[0]}$${k1}$ + integer :: io, stat + + open(newunit=io, file=filename, form="unformatted", access="stream", iostat=stat) + if (stat == 0) then + write(io, iostat=stat) npy_header(vtype, shape(array)) + end if + if (stat == 0) then + write(io, iostat=stat) array + end if + close(io, iostat=stat) + + if (present(iostat)) then + iostat = stat + else if (stat /= 0) then + call error_stop("Failed to write array to file '"//filename//"'") + end if + + if (present(iomsg)) then + iomsg = "Failed to write array to file '"//filename//"'" + end if + end subroutine save_npy_${t1[0]}$${k1}$_${rank}$ + #:endfor +#:endfor + +end submodule stdlib_io_npy_save diff --git a/src/tests/io/CMakeLists.txt b/src/tests/io/CMakeLists.txt index d9a8bc5fc..bfac1b257 100644 --- a/src/tests/io/CMakeLists.txt +++ b/src/tests/io/CMakeLists.txt @@ -13,5 +13,6 @@ ADDTEST(savetxt_qp) set_tests_properties(loadtxt_qp PROPERTIES LABELS quadruple_precision) set_tests_properties(savetxt_qp PROPERTIES LABELS quadruple_precision) +ADDTEST(npy) ADDTEST(open) ADDTEST(parse_mode) diff --git a/src/tests/io/Makefile.manual b/src/tests/io/Makefile.manual index e35beeccb..b07c0ee47 100644 --- a/src/tests/io/Makefile.manual +++ b/src/tests/io/Makefile.manual @@ -6,6 +6,7 @@ SRCGEN = $(SRCFYPP:.fypp=.f90) PROGS_SRC = test_loadtxt.f90 \ test_savetxt.f90 \ + test_npy.f90 \ test_parse_mode.f90 \ test_open.f90 \ $(SRCGEN) diff --git a/src/tests/io/test_npy.f90 b/src/tests/io/test_npy.f90 new file mode 100644 index 000000000..405cad3b3 --- /dev/null +++ b/src/tests/io/test_npy.f90 @@ -0,0 +1,658 @@ +module test_npy + use stdlib_kinds, only : int8, int16, int32, int64, sp, dp + use stdlib_io_npy, only : save_npy, load_npy + use testdrive, only : new_unittest, unittest_type, error_type, check + implicit none + private + + public :: collect_npy + +contains + + !> Collect all exported unit tests + subroutine collect_npy(testsuite) + !> Collection of tests + type(unittest_type), allocatable, intent(out) :: testsuite(:) + + testsuite = [ & + new_unittest("read-rdp-r2", test_read_rdp_rank2), & + new_unittest("read-rdp-r3", test_read_rdp_rank3), & + new_unittest("read-rsp-r1", test_read_rsp_rank1), & + new_unittest("read-rsp-r2", test_read_rsp_rank2), & + new_unittest("write-rdp-r2", test_write_rdp_rank2), & + new_unittest("write-rsp-r2", test_write_rsp_rank2), & + new_unittest("write-i2-r4", test_write_int16_rank4), & + new_unittest("invalid-magic-number", test_invalid_magic_number, should_fail=.true.), & + new_unittest("invalid-magic-string", test_invalid_magic_string, should_fail=.true.), & + new_unittest("invalid-major-version", test_invalid_major_version, should_fail=.true.), & + new_unittest("invalid-minor-version", test_invalid_minor_version, should_fail=.true.), & + new_unittest("invalid-header-len", test_invalid_header_len, should_fail=.true.), & + new_unittest("invalid-nul-byte", test_invalid_nul_byte, should_fail=.true.), & + new_unittest("invalid-key", test_invalid_key, should_fail=.true.), & + new_unittest("invalid-comma", test_invalid_comma, should_fail=.true.), & + new_unittest("invalid-string", test_invalid_string, should_fail=.true.), & + new_unittest("duplicate-descr", test_duplicate_descr, should_fail=.true.), & + new_unittest("missing-descr", test_missing_descr, should_fail=.true.), & + new_unittest("missing-fortran_order", test_missing_fortran_order, should_fail=.true.), & + new_unittest("missing-shape", test_missing_shape, should_fail=.true.) & + ] + end subroutine collect_npy + + subroutine test_read_rdp_rank2(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + character(len=*), parameter :: dict = & + "{'descr': ' Error handling + type(error_type), allocatable, intent(out) :: error + + character(len=*), parameter :: dict = & + "{'descr': ' Error handling + type(error_type), allocatable, intent(out) :: error + + character(len=*), parameter :: dict = & + "{'descr': ' Error handling + type(error_type), allocatable, intent(out) :: error + + character(len=*), parameter :: dict = & + "{'descr':' Error handling + type(error_type), allocatable, intent(out) :: error + + integer :: stat + character(len=*), parameter :: filename = ".test-rdp-r2-rt.npy" + real(dp), allocatable :: input(:, :), output(:, :) + + allocate(input(10, 4)) + call random_number(input) + call save_npy(filename, input, stat) + + call check(error, stat, "Writing of npy file failed") + if (allocated(error)) return + + call load_npy(filename, output, stat) + call delete_file(filename) + + call check(error, stat, "Reading of npy file failed") + if (allocated(error)) return + + call check(error, size(output), size(input)) + if (allocated(error)) return + + call check(error, any(abs(output - input) <= epsilon(1.0_dp)), & + "Precision loss when rereading array") + end subroutine test_write_rdp_rank2 + + subroutine test_write_rsp_rank2(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + integer :: stat + character(len=*), parameter :: filename = ".test-rsp-r2-rt.npy" + real(sp), allocatable :: input(:, :), output(:, :) + + allocate(input(12, 5)) + call random_number(input) + call save_npy(filename, input, stat) + + call check(error, stat, "Writing of npy file failed") + if (allocated(error)) return + + call load_npy(filename, output, stat) + call delete_file(filename) + + call check(error, stat, "Reading of npy file failed") + if (allocated(error)) return + + call check(error, size(output), size(input)) + if (allocated(error)) return + + call check(error, any(abs(output - input) <= epsilon(1.0_dp)), & + "Precision loss when rereading array") + end subroutine test_write_rsp_rank2 + + subroutine test_write_int16_rank4(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + integer :: stat, i + character(len=*), parameter :: filename = ".test-i2-r4-rt.npy" + integer(int16), allocatable :: input(:, :, :, :), output(:, :, :, :) + + input = reshape([(i*(i+1)/2, i = 1, 40)], [2, 5, 2, 2]) + call save_npy(filename, input, stat) + + call check(error, stat, "Writing of npy file failed") + if (allocated(error)) return + + call load_npy(filename, output, stat) + call delete_file(filename) + + call check(error, stat, "Reading of npy file failed") + if (allocated(error)) return + + call check(error, size(output), size(input)) + if (allocated(error)) return + + call check(error, all(abs(output - input) == 0), & + "Precision loss when rereading array") + end subroutine test_write_int16_rank4 + + subroutine test_invalid_magic_number(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + character(len=*), parameter :: dict = & + "{'descr': ' Error handling + type(error_type), allocatable, intent(out) :: error + + character(len=*), parameter :: dict = & + "{'descr': ' Error handling + type(error_type), allocatable, intent(out) :: error + + character(len=*), parameter :: dict = & + "{'descr': ' Error handling + type(error_type), allocatable, intent(out) :: error + + character(len=*), parameter :: dict = & + "{'descr': ' Error handling + type(error_type), allocatable, intent(out) :: error + + character(len=*), parameter :: dict = & + "{'descr': ' Error handling + type(error_type), allocatable, intent(out) :: error + + character(len=*), parameter :: dict = & + "{'descr': ' Error handling + type(error_type), allocatable, intent(out) :: error + + character(len=*), parameter :: dict = & + "{'fortran_order': True, 'shape': (10, 4, ), 'descr': ' Error handling + type(error_type), allocatable, intent(out) :: error + + character(len=*), parameter :: dict = & + "{'fortran_order': True,, 'shape': (10, 4, ), 'descr': ' Error handling + type(error_type), allocatable, intent(out) :: error + + character(len=*), parameter :: dict = & + "{'fortran_order': True, 'shape': (10, 4, ), 'descr': ' Error handling + type(error_type), allocatable, intent(out) :: error + + character(len=*), parameter :: dict = & + "{'descr': ' Error handling + type(error_type), allocatable, intent(out) :: error + + character(len=*), parameter :: dict = & + "{'fortran_order': True, 'shape': (10, 4, ), } " // & + char(10) + character(len=*), parameter :: header = & + char(int(z"93")) // "NUMPY" // char(1) // char(0) // & + char(len(dict)) // char(0) // dict + + integer :: io, stat + character(len=:), allocatable :: msg + character(len=*), parameter :: filename = ".test-missing-descr.npy" + real(dp), allocatable :: array(:, :) + + open(newunit=io, file=filename, form="unformatted", access="stream") + write(io) header + write(io) spread(0.0_dp, 1, 40) + close(io) + + call load_npy(filename, array, stat, msg) + call delete_file(filename) + + call check(error, stat, msg) + end subroutine test_missing_descr + + subroutine test_missing_fortran_order(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error + + character(len=*), parameter :: dict = & + "{'descr': ' Error handling + type(error_type), allocatable, intent(out) :: error + + character(len=*), parameter :: dict = & + "{'fortran_order': True, 'descr': ' 0) then + write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!" + error stop + end if +end program