Hi Nathan,

thank you for your great work on UTF8 strings and their integration in Numpy. This is a very important dtype to support, especially with the widespread use of large language models (LLM) nowadays.
However, I would like to comment on the serialization. Hope it's not too late at this point (the last time I looked I think the serialization was less detailed), but the approach with sidecar_size has the disadvantage that Numpy arrays would not be efficiently appendable anymore. Also, while it's comfortable to not need to specify how the sidecar data looks like, it makes the format more proprietary, less open and more difficult to debug. To be precise, I see this part of NEP 1 violated:

Be reverse engineered. Datasets often live longer than the programs that created them. A competent developer should be able to create a solution in his preferred programming language to read most NPY files that he has been given without much documentation.

On top of that, what if the data format for the sidecar data changes? Is it then still possible to read old files with a newer Numpy version?

To overcome this issues, first and foremost I suggest not to introduce a .npy file format version 4.0 with this kind of fundamental change. I suggest to instead of adding the sidecar data to the same .npy file to add it to an extra, standard .npy file. Here an example: if the user wants to save to "mystrings.npy", the files "mystrings.npy" and "mystrings.npy.idx" could be created where latter is also just a regular .npy file and contains indices/offsets into mystrings.npy that would otherwise end up in the sidecar data. When the file "mystrings.npy" is loaded, Numpy checks whether there is a mystrings.npy.idx and tries to load it as well. This is an approach I use pretty regularly with video data: the single frames are all in one big (appendable) array and the index array contains the begin indices/offsets and (redundantly, for fast lookups, sort etc.) the lengths and end offsets. This concept is very generic, can be used for all sorts of ragged arrays including text, satisfies the requirements of NEP 1 and is efficiently appendable (or at least moves the burden from the programmer to the filesystem).

Best, Michael

Also, a dedicated index array can come with a custom shape, which would allow for multidimensional ragged arrays (in the future). Wouldn't that be great?

Thankfully, I haven't started work on the serialization yet. I wasn't aware that there was already a capability to add extra data in the npy format. I will look at this and make sure it's usable for my case and if so, update the NEP to drop the part about updating the npy format. I only added that because I thought there was no other option - I'm glad to hear I was wrong!

Okay, good to hear I'm not too late 😅 It's just an alternative idea where to house the sidecar data and of course for discussion. I'll send it to the Numpy user group, too. Also, how those .idx files would be handled would still need to be defined. Just as a definition of how the sidecar data would look like for the String dtype and all aspects around that. Two options come to my mind:

1. check if some dtype is String and then look for an .npy.idx file. If that fails, raise some error
2. check for .npy.idx files every time a .npy file is loaded (could be variant for ragged arrays)

Ragged arrays are not defined at all in Numpy yet, afaik. Just described a solution I used frequently (which could be generalized to do the trick). Let's see what others think.

Actually, I've spend quite some time thinking about this entire topic and I think I now came up with a solution that could fit very well with how things are done in Numpy.

Above, I've mentioned that I've used index files for Numpy files to create some sort of ragged array. That's production proven and works and users are happy. However, it would probably make more sense to do it the other way round: to have the index file as the main .npy file and declare what used to be the main .npy file as some sort of heap. I've mentioned that idea in a Numpy community meeting once. The advantage would be a stronger hierarchy: one main file with indices/offsets/pointers that point to heaps in other files. Let me give an example what that could mean for strings: let's say the main .npy file "main.npy" has a shape of (100, 100) and that new String dtype. Then, the actual strings would all be lined up one after another in "main.npy.strings" (also a regular Numpy file, e.g. dtype byte or char) and "main.npy" could e.g. contain int64 values that can be interpreted as byte offsets to "main.npy.strings". Or dual int64 values that contain byte offsets of the string begin and string end (or string length). Just as a very small example to get some opinions.

What do you think?

Can we open a new issue for the storage disucssion? It's interesting, but I don't think it is central to this PR.

I'll open an issue in the coming week.

Maybe reviewing the current tests in the numpy-user-dtype repo may actually be useful.

Those tests are coming over here too.

One thing I am curious: Is there a reason to keep static_string.c in the multiarray folder rather than the new stringdtype one?

No particular reason. I wrote static_string.c such that the dtype doesn't need to know anything about the implementation of the string library and thought it might be generally useful elsewhere. Happy to move it back in with the dtype since it does logically belong there.

(just missed if/where we reclaim unused memory).

The memory gets reclaimed by the OS in NpyString_free_allocator, which is called in stringdtype_dealloc.

Can we open a new issue for the storage disucssion? It's interesting, but I don't think it is central to this PR.

Okay, I've done that, see #25374. However, I think serialization is important for the String dtype, because otherwise one could just keep using object arrays if I understand correctly.

It looks like all the tests pass except on pypy. Woohoo! I know at one point we had some concerns about how things would work on a 32 bit python or a big endian system so it's good to see that everything seems to be working with only very minor adjustments (see the last few commits).

The pypy failures look like they're happening inside the stringdtype tests and are crashing the pypy interpreter (I think?), so I'm going to need to figure out how I'm breaking pypy. Maybe @mattip we could pair on this sometime if you're interested in taking a look. Some RPython tracebacks showing up in the build log:

RPython traceback:
  File "pypy_interpreter.c", line 47668, in BuiltinCode_funcrun_obj
  File "pypy_module_cpyext_5.c", line 64539, in generic_cpy_call__StdObjSpaceConst_funcPtr_SomeI_11
  File "pypy_module_cpyext_3.c", line 52360, in CpyTypedescr_str_realize
  File "pypy_module_cpyext_6.c", line 19973, in unicode_realize
  File "pypy_interpreter.c", line 32661, in _unibuf_to_utf8

Next up after that are docs and serialization.

The only remaining TODO comment in the codebase is because I never implemented casts to and from extended precision types like longdouble. Do I need to do that? It's not terribly hard but also I have other things to worry about.

Also I could add a cast to and from bytes, erroring out if the bytes aren't valid UTF-8 (say).

Also I could add a cast to and from bytes, erroring out if the bytes aren't valid UTF-8 (say).

Bytes or void? For bytes, I think it is either latin1 or error out if not ascii? For void, you could allow it but check for valid utf-8. Bytes cast seems like it should exist to me.

Now coercing Python byte objects could and maybe should raise, I think.

I never implemented casts to and from extended precision types like longdouble.

I would consider it a bug to miss them, but I can live with considering it an unfixed bug. Would be nice to have casts for all builtin dtypes though.

@mattip take a look at 743084e. Happy to entertain alternative approaches.

I also tried just now to report a bug on the pypy bug tracker about the incorrect implementation of Py_Is on pypy 3.10, but heptapod thought my bug report was spam and refused to let me post it...

I isolated the windows pypy crash and reported it upstream: https://foss.heptapod.net/pypy/pypy/-/issues/4046

It turns out we can rely on pointer comparison for singletons, just not for NaN. I refactored things to first do a pointer comparison, then do a NaN check, and finally check for == equality. That seems to work fine on CPython and pypy, so no pypy-specific hacks needed. I also added a new pypy-specific hack to work around the crash I saw last week, which should fix CI and avoid crashes on older unpatched pypy versions.

Bytes or void? For bytes, I think it is either latin1 or error out if not ascii? For void, you could allow it but check for valid utf-8. Bytes cast seems like it should exist to me.

Coming back to this - I'm fine with it, we just need to document it.

I could say something like "If you need to work with memory-mapped UTF-8 rows, use a void dtype to load it as a numpy array and then cast it to stringdtype if you need to treat it as strings. If you have Latin-1 strings that don't contain embedded nulls (FITS?), memory-map it as the bytes dtype and then cast it to stringdtype if you want to treat it as variable-width python strings". And then mention that both casts need a copy since stringdtype's storage model is fundamentally different from the other dtypes.

I finished out all the casts, please feel free to look over the last few commits I just pushed. I ended up deciding to make the bytes casts only work with ASCII data, to match the existing unicode to bytes casts. I also decided to allow the string to void cast to truncate a UTF-8 character without warning or erroring.

There should be a github reward for doing careful code review on a 10,000 line PR. Thank you so so so much for your time on this. I'm very impressed with your ability to spot C bugs. I'll finish up the comments on the tests tomorrow and hopefully we'll be able to merge this soon!

Thanks! It has actually been not too bad reviewing, since your code is really quite clear (but indeed doing the lot meant getting sometimes a bit grumpy about boilerplate!), and I felt that after looking over the development of the new dtype and array method machinery, it made sense to try to understand an actual application.

Also a question: the hardest part about reviewing is to see what is missing and I realized that I don't think I saw many tests with views (either implicit as one gets by slicing, or explicit .view()). For an implicit view (say a[:5]) presumably the StringDType is shared (or at least the allocator); this matters a bit for in-place operations like a += a[:] - for that, the iterator will make internal copies, so sin != sout; is your dtype indeed also not shared? But they have the same allocator? What about a[a.startswith("a")] (copies array data but what happens to the arena? copied as well?).

Though really the above came up beause I wondered what happens if I do a.view(StringDType(...)), perhaps to replace the NaN object. Does the storage get properly shared?

Anyway, can obviously look at the code again, but not at this hour...

beause I wondered what happens if I do a.view(StringDType(...))

We can follow up on this, I think. But would be good to have a test that this just fails due to the HASREFS flag being checked on the dtype instance (I am assuming it fails, I didn't test).

It looks like we need more tests with views...

In [1]: a = np.array([['0', '±'*20], [';'*10, 'X']], 'T')

In [2]: a += a[0]
---------------------------------------------------------------------------
MemoryError                               Traceback (most recent call last)
Cell In[2], line 1
----> 1 a += a[0]

MemoryError: Failed to allocate string in string to string cast

In [3]: a
Out[3]: 
array([['00',
        '\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00'],
       [';;;;;;;;;;', 'X']], dtype=StringDType())

Things that work:

a += a[0].copy()
a += a
a += a[:]

Also no problem for short strings only:

 In [1]: a = np.array([['0', '±'*3], [';'*6, 'X']], 'T')

In [2]: a += a.T

In [3]: a
Out[3]: 
array([['00', '±±±;;;;;;'],
       [';;;;;;±±±', 'XX']], dtype=StringDType())

In [4]: a += a.T
---------------------------------------------------------------------------
MemoryError                               Traceback (most recent call last)
Cell In[4], line 1
----> 1 a += a.T

MemoryError: Failed to allocate string in string to string cast

For an implicit view (say a[:5]) presumably the StringDType is shared (or at least the allocator); this matters a bit for in-place operations like a += a[:] - for that, the iterator will make internal copies, so sin != sout; is your dtype indeed also not shared? But they have the same allocator? What about a[a.startswith("a")] (copies array data but what happens to the arena? copied as well?).

Thanks for pointing this out, I wasn't aware of this wrinkle with the iterators.

I'd rather not completely ban views, since I found it convenient with pandas in a few places to change the na_object temporarily without a copy.

It turns out there are two issues I need to fix to get this to work. First, I need a way to detect whether a descriptor actually owns the array data, or if it is a view into another array. It looks like setting array_owned=2 and checking for that in a few places avoids double-frees when cleaning up descriptors and allocators.

The other issue is that there's no way to easily determine in the ufunc loop if two different packed strings happen to share a heap allocation. As you said, the iterator makes a copy, so the packed strings pointers are distinct and the input and output array pointers are distinct, but they still share an underlying heap allocation.

One way to avoid this issue would be to assume if there are any shared allocators between the input and output arrays that there may be an in-place operation happening, and in that case ufuncs should allocate a temporary output buffer, do the operation, and only then pack the string data into the output array.

Of course that means we'll be doing copies unnecessarily inside ufunc loops for any possibly in-place operation, even if it isn't actually in-place.

Another option would be to add a new C API function along the lines of:

int NpyString_share_memory(npy_packed_static_string *s1, npy_string_allocator *alloc1,
                           npy_packed_static_string *s2, npy_string_allocator *alloc2)

which would return a positive value if the two strings share a heap allocation, 0 if they don't, and a negative value to indicate memory corruption or a logic error (e.g. s1 or s2 have flags set indicating a heap allocation, but alloc1 or alloc2 have a null arena pointer).

Do either of you two (or anyone else still following along) have an opinion on this? I think the second option is probably preferable and doable, but I'll need to make sure I audit all the ufunc loops to make sure they have such a check and add various tests that do operations with an array and a slice into the same array and an array and a view into that array.

Should I fix this immediately before this PR is merged or should we wait until it's easier to do code review on the fix?

Maybe disabling the broken views for now could make sense in this PR to undo that later? But if not also fine. I have to think about this more to make a good input, but making a larger iteration here also doesn't seem helpful.

I similarly would need more time... But maybe you were right all along to compare the dtype to see whether things were in-place? It is true that views share the dtype, i.e., a.T.dtype is a.dtype. Though I cannot immediately reason it through, it may be worth reverting to what you had before and check whether my examples "just work".

But maybe you were right all along to compare the dtype to see whether things were in-place? It is true that views share the dtype, i.e., a.T.dtype is a.dtype. Though I cannot immediately reason it through, it may be worth reverting to what you had before and check whether my examples "just work".

Thanks for this hint! This does mostly work, although I also need to fix the issue alluded to above in the string to string cast, since in that case I do need to do the more thorough memory sharing check that looks at the unpacked representation. So for now I've added a private version of NpyString_share_memory and if we need to expose it publicly in the future we can come back to making it public. It turns out the identical descriptor check is sufficient in the ufunc loops, at least in the tests I have so far.

While doing this after applying some of the other test changes you suggested above, I uncovered two different bugs in the arena allocator implementation related to handling mutated strings. Glad I caught these!

Unfortunately that took all day and I didn't have enough time to add the tests you wanted with views and slices, but I did verify they work now and the errors you ran into before are fixed.

Great to hear that things mostly work! I'm now reviewing an actual paper, but can have a last look at the ufunc loops/tests later. It will be very nice to have this in and be able to look at smaller bits & pieces...

@ngoldbaum - I slightly lost track of where we are with this -- let me know when it is ready for final review on the ufuncs and associated tests...

It's ready for final review whenever you're ready.

OK, great! Tomorrow hopefully!

The last push addresses today's comments. Thanks @seberg for the suggestion to refactor the multiply promoters, it's much simpler now with a single promoter.

Once Marten approves I'm happy to squash this into one enormous commit. There might be some history in the PR that would be useful in the future though so I'm also happy to leave it as-is. It should at least build on most systems on every commit.

Sounds good. However you like it, I would normally squash, but some big histories are just neat to see how much effort went in ;).
I thought git bisect is smart enough to prefer splitting on merge commits anyway, so even if there is a bad commit it doesn't destroy bisecting a lot? (But I admit, I don't bisect often in practice).

Yeah, you just git bisect skip if you hit a commit that doesn't build (there are definitely commits in the numpy history and most open source projects on github that don't build...). It also can handle arbitrary history graphs with merging.

OK, in that case... can I hit the merge button? 😁

Yes go ahead! Let's have this wonderful addition!!

OK, squashed everything into one commit. I'll merge it tonight once all the tests are green, just in case.