« first day (3887 days earlier)      last day (1061 days later) » 
00:00 - 23:0023:00 - 00:00

11:00 PM
alright folks, I'm calling it a day. Have a good one, everybody! See y'all tomorrow
 
take care
 
Ok, Andras-enhanced prototype is producing good output (as if there was any doubt), now to shave the various yaks that are preventing me from plugging my real data in
 
there's always doubt with untested code of mine :P
 
When I tested it, the waveform collapsed and it became non-untested, which means that it must have had been working the whole time. Ipso facto.
 
Let me share a barely-interesting factoid: it might be surprising at first glance that the eigenvector matrix (array) contains the eigenvectors as columns, because numpy's default memory layout is such that arr[i, :] is a contiguous block of memory rather than arr[:, i]. I suspect the reason for the choice to have vectors[:, i] be a single eigenvector is so that M @ vectors acts on each eigenvector according to their definition.
(i.e. the result has the eigenvectors as columns, but each one is multiplied by the corresponding eigenvalue)
Still I have to look up the documentation or do a manual check every single time to ensure I'm using the right convention. This is extra fun when you're doing diagonalization and you need to do something like vectors.T @ M @ vectors... or the other way around.
 
11:10 PM
Makes sense. Maybe they chose the layout that optimizes and/or simplifies the most common use case of eig.
If multiplying each eigenvector by its respective eigenvalue is in high demand among users, there you go.
I definitely got the row/column ordering wrong in my first prototype. I was relieved that switching my indices around gave me sensible output, because I was nearly out of ideas
 
Huh, now this is interesting
In [62]: arr = np.random.rand(3, 3)
    ...: evals, evects = np.linalg.eig(arr)

In [63]: evals
Out[63]: array([ 1.53688253, -0.01145095,  0.36737775])

In [64]: evals.base
Out[64]: array([ 1.53688253+0.j, -0.01145095+0.j,  0.36737775+0.j])
numpy seems to detect that every eigenvalue is real for this real matrix, and return only the real part of eigenvectors and eigenvalues (which are implemented as views in numpy)
This explains a surprising feature of the eigenvector matrix:
In [66]: evects.flags['C_CONTIGUOUS']
Out[66]: False

In [67]: evects.flags['F_CONTIGUOUS']
Out[67]: False
It's not contiguous in either direction.
(Because complexes are implemented by interleaving real and imaginary parts in the array, if memory serves)
 
Can you access the non-view complex array of the eigenvectors? Is that likely to be contiguous?
 
@Kevin yes, evects.base is that original array, and it's C contiguous
In [75]: evects.base.flags
Out[75]:
  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  WRITEBACKIFCOPY : False
  UPDATEIFCOPY : False
This is the same thing:
In [74]: arr = np.arange(3) + 1j*np.arange(3)
    ...: np.may_share_memory(arr.real, arr.imag), np.shares_memory(arr.real, arr.imag)
Out[74]: (True, False)
 
Hmm neat
 
may_share_memory() only checks that the memory wrapped by the arrays have an overlap in domains. shares_memory() actually solves the Diophantine equations telling us whether there's any coincident memory addresses.
since we have an alternation of real and imaginary double values in the complex128 array, the former is True and the latter is False
 
11:22 PM
I'm surprised that it needs a Diophantine equation to determine the result. And numpy.org/doc/stable/reference/generated/… says the function can be exponentially slow... I guess arrays can have fancier layouts than I imagine
 
@MisterMiyagi do you do numpy? If so you might be interested in mail.python.org/pipermail/numpy-discussion/2021-June/…
@Kevin yeah, when you create a view with a slice that has a step other than 1 you start skipping over elements. That's how you get arrays tjhat are overall non-contiguous.
you have to determine whether a1 + n*b1 can equal a2 + m*b2 for any n, m that belongs to the respective arrays
 
Makes sense... And I guess this specific case is pretty cheap to solve since the steps are equal
 
This is pretty much what's happening with real/imaginary parts:
In [82]: arr = np.arange(3) + 1j*np.arange(4, 7)
    ...: arr.dtype
Out[82]: dtype('complex128')

In [83]: arr.view('float64')[::2]
Out[83]: array([0., 1., 2.])

In [84]: arr.view('float64')[1::2]
Out[84]: array([4., 5., 6.])
@Kevin yeah, although I don't know how brute force the solver is (probably a lot, eventually).
I went source diving a while back exactly to learn how may_share_memory vs shares_memory worked for a question on SO
 
Let us hope they at least have a special case for if step_a == 1 and step_b == 1:
 
11:30 PM
Or maybe the guy in charge of shares_memory is saving that one for a rainy day. Performance review coming up? Check in a change that reduces an O(N) case to O(1), where N is likely to be in the millions in practice
 
"strides" are the term for the step size (in bytes) in memory for a change of 1 in an array index
Oh yeah, and of course with multidimensional arrays you get more variables to find overlaps for.
I suspect that's when issues start to happen
 
I don't know what half these variables mean, but github.com/numpy/numpy/blob/maintenance/1.11.x/numpy/core/src/… emits a vibe of "solves stride-1 problems in O(1) time"
 
Yeah, I guess that would make sense. With trivial strides you just have to check the extents, i.e. what may_share_memory does by default.
 
I'm sort of assuming without proof that diophantine_simplify can reduce nterms to 1 if the arrays' layouts share certain properties
Listener: "elaborate upon which properties you mean"
Me: "No."
 
Hmm, I wouldn't vouch for the meaning of n there either way. I don't know what that corresponds to.
it could mean "trivial strides for a 1d array" but it could also mean "1d strides vs 0d array"
in other words, I don't have intuition for translating the memory overlap problem to the Diophantine equations
 
11:39 PM
Yes, unfortunately my vibe detector can't give me information at that resolution
The UI is just a happy green light and an angry red light
 
efinancialcareers.com/news/2021/06/… apparently python could be 5x faster
Crap sent early
 
I'm pretty satisfied with Python's speed as it is, so any improvement will be a nice bonus
Finally, I can compute all the primes below five million, not just one million!
 
@12944qwerty slowpoke
 
Before trying to speed up the program, first seek to speed up the programmer
 
Hi, Anyone with html expertise, I cannot find a chat room, specialised in html
 
11:49 PM
@YatShan please don't come here asking explicitly python-unrelated questions
 
Ok, sorry for coming
 
Thanks a lot @AndrasDeak
 
I'm trying to figure out how bad the Diophantine equation can get, but I can't slow it down too much
In [145]: base = np.arange(4**12).reshape((4,)*12)
     ...: even_slicer = (slice(None, None, 2),) * 12
     ...: odd_slicer = (slice(1, None, 2),) * 12
     ...: arr1, arr2 = base[even_slicer], base[odd_slicer]  # 12d arrays
     ...: arr3, arr4 = base.ravel()[::4096], base.ravel()[4096//2::4096]  # 1d arrays
     ...: %timeit np.shares_memory(arr1, arr2)
     ...: %timeit np.shares_memory(arr3, arr4)
1.9 µs ± 9.58 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
699 ns ± 9.33 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
the comparison is not entirely fair because the 1d arrays have a more regular stride pattern
 
00:00 - 23:0023:00 - 00:00

« first day (3887 days earlier)      last day (1061 days later) »