[ot][spam][random][crazy][random][crazy]
Undescribed Horrific Abuse, One Victim & Survivor of Many
gmkarl at gmail.com
Fri Nov 18 04:49:39 PST 2022
battery is at 45% !
0541
i fixed a mistake and two more assertions are passing.
it's paused on the failing one:
269 rfreqs15t = fftfreq(repetition_samples=15, complex=False)
270 rfreqs15f = fftfreq(15, complex=False)
271 irft_from_15 = create_freq2time(freqs=rfreqs15f)
272 rft_from_15 = create_time2freq(15, freqs=rfreqs15t)
273 irft_to_15 = create_freq2time(15, freqs=rfreqs15t)
274 randvec2rtime15 =
randvec[:15].astype(np.complex128).view(float) @ irft_from_15
275 randvec2rfreq15 = randvec[:15] @ rft_from_15
276 randvec2irfft = np.fft.irfft(randvec[:15])
277 randvec2rfft = np.fft.rfft(randvec[:15])
278 randvec2rfreq152randvec = randvec2rfreq15.view(float) @ irft_to_15
279 randvec2rfreq152irfft = np.fft.irfft(randvec2rfreq15, 15)
280 assert np.allclose(rfreqs15t, np.fft.rfftfreq(15))
281 assert np.allclose(rfreqs15f, np.fft.rfftfreq(28))
282 assert np.allclose(randvec2rtime15, randvec2irfft)
283 assert np.allclose(randvec2rfreq15, randvec2rfft)
284 -> assert np.allclose(randvec2rfreq152randvec, randvec2rfreq152irfft)
285 assert np.alloclose(randvec2rfreq152randvec, randvec[:15])
This is likely the time-domain result of complex frequencies. my plan
was to break in the matrix generation, and run test data through the
steps to see when it differed.
134 complex_freqs = np.concatenate([
135 freqs[:tail_idx],
136 -freqs[tail_idx:],
137 -freqs[head_idx:tail_idx][::-1]
138 ])
139 complex_mat = complex_wavelet(wavelet,
np.outer(complex_freqs, offsets))
162 interim_mat = complex_mat.copy()
163 interim_mat[head_idx:tail_idx] += interim_mat[len(freqs):][::-1]
164 mat = np.stack([
165 interim_mat[:len(freqs)].real.T,
166 -interim_mat[:len(freqs)].imag.T
167 ], axis=-1).reshape(len(complex_freqs), len(freqs)*2).T
185 return mat / freq_count # scaled to match numpy convention
note to self: the help said you can change when that / freq_count
normalization happens in numpy, i believe.
(Pdb) break 162
Breakpoint 1 at /shared/src/scratch/fourier.py:162
-> interim_mat = complex_mat.copy()
(Pdb) time_data = np.random.random(complex_mat.shape[-1])
(Pdb) complex_freq_data = time_data @ np.linalg.inv(complex_mat)
(Pdb) real_freq_data = np.concatenate([complex_freq_data[:tail_idx],
complex_freq_data[tail_idx:len(freqs)].conj()])
...
(Pdb) p complex_freq_data @ complex_mat
array([0.0202184 +1.38777878e-16j, 0.83261985-2.77555756e-17j,
0.77815675+2.77555756e-16j, 0.87001215-1.38777878e-16j,
0.97861834-1.24900090e-16j, 0.79915856-1.94289029e-16j,
0.46147936-2.22044605e-16j, 0.78052918-1.80411242e-16j,
0.11827443-1.04083409e-16j, 0.63992102+3.60822483e-16j,
0.14335329-1.04083409e-16j, 0.94466892-1.87350135e-16j,
0.52184832+2.74086309e-16j, 0.41466194+0.00000000e+00j,
0.26455561-1.73472348e-17j, 0.77423369+5.55111512e-17j,
0.45615033-5.13478149e-16j, 0.56843395+2.77555756e-16j,
0.0187898 -4.16333634e-17j, 0.6176355 -1.66533454e-16j,
0.61209572+1.66533454e-16j, 0.616934 -3.33066907e-16j,
0.94374808+6.24500451e-17j, 0.6818203 +2.84494650e-16j,
0.3595079 +2.63677968e-16j, 0.43703195+3.60822483e-16j,
0.6976312 -2.70616862e-16j, 0.06022547+3.60822483e-16j])
(Pdb) p real_freq_data.view(float) @ mat
array([0.0202184 , 0.44642266, 0.73789397, 0.65352205, 0.66906312,
0.74048943, 0.70261372, 0.69873159, 0.36518507, 0.62877826,
0.08107154, 0.75655143, 0.48899933, 0.59444781, 0.26455561,
0.59444781, 0.48899933, 0.75655143, 0.08107154, 0.62877826,
0.36518507, 0.69873159, 0.70261372, 0.74048943, 0.66906312,
0.65352205, 0.73789397, 0.44642266])
i'll pick sample index 6 to debug again i think.
(Pdb) p (complex_freq_data @ complex_mat)[6]
(0.4614793622529322-2.220446049250313e-16j)
(Pdb) p (real_freq_data.view(float) @ mat)[6]
0.7026137203837779
(Pdb) p complex_freq_data @ complex_mat[:,6]
(0.4614793622529323-2.6020852139652106e-16j)
(Pdb) p real_freq_data.view(float) @ mat[:,6]
0.7026137203837779
I tried the nyquist frequency, and it looks like it's working right:
(Pdb) p complex_freq_data[tail_idx] * complex_mat[tail_idx,6]
(-0.09512353353755763-2.275865808140467e-17j)
(Pdb) p real_freq_data.view(float)[tail_idx*2:] @ mat[tail_idx*2:,6]
-0.09512353353755763
so I guess the problem is in combining the conjugates in the interim matrix.
(Pdb) p complex_mat[:,6][head_idx:tail_idx].round(3)
array([ 0.223+0.975j, -0.901+0.434j, -0.623-0.782j, 0.623-0.782j,
0.901+0.434j, -0.223+0.975j, -1. +0.j , -0.223-0.975j,
0.901-0.434j, 0.623+0.782j, -0.623+0.782j, -0.901-0.434j,
0.223-0.975j])
(Pdb) p complex_mat[:,6][len(freqs):][::-1].round(3)
array([ 0.223-0.975j, -0.901-0.434j, -0.623+0.782j, 0.623+0.782j,
0.901-0.434j, -0.223-0.975j, -1. +0.j , -0.223+0.975j,
0.901+0.434j, 0.623-0.782j, -0.623-0.782j, -0.901+0.434j,
0.223+0.975j])
ok, yeah, as conjugates, i lose information if I just sum them in
complex space. i need to negate the imaginary parts of the conjugates.
0556
interim_mat = complex_mat.copy()
interim_mat[head_idx:tail_idx].real +=
interim_mat[len(freqs):][::-1].real
interim_mat[head_idx:tail_idx].imag -=
interim_mat[len(freqs):][::-1].imag
mat = np.stack([
interim_mat[:len(freqs)].real.T,
-interim_mat[:len(freqs)].imag.T
], axis=-1).reshape(len(complex_freqs), len(freqs)*2).T
0600
interim_mat = complex_mat.copy()
interim_mat[head_idx:tail_idx] += interim_mat[len(freqs):][::-1].conj()
mat = np.stack([
interim_mat[:len(freqs)].real.T,
-interim_mat[:len(freqs)].imag.T
], axis=-1).reshape(len(complex_freqs), len(freqs)*2).T
0603
the same assertion is still failing. i figured i had just done the
wrong fix or poorly, but sample 6 is matching now:
(Pdb) p complex_freq_data[:] @ complex_mat[:,6]
(0.4614793622529323-2.6020852139652106e-16j)
(Pdb) p real_freq_data.view(float)[:] @ mat[:,6]
0.46147936225293196
looks like they all are now, which is good:
(Pdb) p (complex_freq_data @ complex_mat).round(3)
array([0.02 +0.j, 0.833-0.j, 0.778+0.j, 0.87 -0.j, 0.979-0.j, 0.799-0.j,
0.461-0.j, 0.781-0.j, 0.118-0.j, 0.64 +0.j, 0.143-0.j, 0.945-0.j,
0.522+0.j, 0.415+0.j, 0.265-0.j, 0.774+0.j, 0.456-0.j, 0.568+0.j,
0.019-0.j, 0.618-0.j, 0.612+0.j, 0.617-0.j, 0.944+0.j, 0.682+0.j,
0.36 +0.j, 0.437+0.j, 0.698-0.j, 0.06 +0.j])
(Pdb) p (real_freq_data.view(float) @ mat).round(3)
array([0.02 , 0.833, 0.778, 0.87 , 0.979, 0.799, 0.461, 0.781, 0.118,
0.64 , 0.143, 0.945, 0.522, 0.415, 0.265, 0.774, 0.456, 0.568,
0.019, 0.618, 0.612, 0.617, 0.944, 0.682, 0.36 , 0.437, 0.698,
0.06 ])
i tested both generations in their generation function that they made
the same output as a complex matrix, and they appeared to to me.
these two values are supposed to be the same but are differing:
278 randvec2rfreq152randvec = randvec2rfreq15.view(float) @ irft_to_15
279 randvec2rfreq152irfft = np.fft.irfft(randvec2rfreq15, 15)
i've bumped into a couple different likely mistakes that could cause
that difference.
0608
0621
summary of failure:
vector.view(float) @ create_freq2time(15,
freqs=fftfreq(repetition_samples=15,complex=False))
is producing different output than
np.fft.irfft(vector, 15)
- there is no nyquist frequency
- the real frequencies used are rfreqs15t
- the test vector is randvec2rfreq15
- the failing matrix is constructed as irft_to_15 on line 273
so, inside the call on line 273, i should expect the final matrix to
have the same effect as np.fft.irfft(_, 15)
(Pdb) list
134 complex_freqs = np.concatenate([
135 freqs[:tail_idx],
136 -freqs[tail_idx:],
137 -freqs[head_idx:tail_idx][::-1]
138 ])
139 -> complex_mat = complex_wavelet(wavelet,
np.outer(complex_freqs, offsets))
140
141 # assuming the input's real and complex components are
broken out into separate elements
142
143 # real_in, imag_in
144 # multiplied by conjugate pair and summed
(Pdb) p is_complex, has_dc_offset, has_nyquist
(False, True, False)
(Pdb) p complex_freqs
array([ 0. , 0.06666667, 0.13333333, 0.2 , 0.26666667,
0.33333333, 0.4 , 0.46666667, -0.46666667, -0.4 ,
-0.33333333, -0.26666667, -0.2 , -0.13333333, -0.06666667])
Something is different before the conversion to real:
(Pdb) time_data = np.random.random(15)
(Pdb) complex_freq_data = np.fft.fft(time_data)
(Pdb) p np.fft.ifft(complex_freq_data).round(3)
array([0.02 -0.j, 0.833-0.j, 0.778-0.j, 0.87 +0.j, 0.979+0.j, 0.799+0.j,
0.461-0.j, 0.781+0.j, 0.118-0.j, 0.64 -0.j, 0.143+0.j, 0.945-0.j,
0.522+0.j, 0.415-0.j, 0.265+0.j])
(Pdb) p (complex_freq_data @ complex_mat / freq_count).round(3)
array([0.022-0.j, 0.892-0.j, 0.834-0.j, 0.932-0.j, 1.049-0.j, 0.856+0.j,
0.494-0.j, 0.836-0.j, 0.127-0.j, 0.686-0.j, 0.154+0.j, 1.012-0.j,
0.559-0.j, 0.444-0.j, 0.283+0.j])
It's strange to encounter this. I wonder what obscure error I have
made for myself.
The frequencies used to make the matrix are the same:
(Pdb) p complex_freqs
array([ 0. , 0.06666667, 0.13333333, 0.2 , 0.26666667,
0.33333333, 0.4 , 0.46666667, -0.46666667, -0.4 ,
-0.33333333, -0.26666667, -0.2 , -0.13333333, -0.06666667])
(Pdb) p np.fft.fftfreq(len(complex_freq_data))
array([ 0. , 0.06666667, 0.13333333, 0.2 , 0.26666667,
0.33333333, 0.4 , 0.46666667, -0.46666667, -0.4 ,
-0.33333333, -0.26666667, -0.2 , -0.13333333, -0.06666667])
Here's the line again it's made with:
139 complex_mat = complex_wavelet(wavelet,
np.outer(complex_freqs, offsets))
I guess it's back down to comparing a single sample or whatnot. I'll
compare sample 6 again.
0634
(Pdb) p (complex_freq_data @ complex_mat[:,6]) /
np.fft.ifft(complex_freq_data)[6]
(15.000000000000009-1.8056631630774395e-15j)
(Pdb) p freq_count
14
Looks like I'm dividing by the wrong denominator.
0636
here's how freq_count is originally made:
117 freq_count = len(freqs)
118 if not has_dc_offset:
119 freq_count += 1
120 if not is_complex:
121 freq_count = (freq_count - 1) * 2
(Pdb) p len(freqs)
8
(Pdb) p has_dc_offset
True
(Pdb) p is_complex
False
(Pdb) p (8 - 1) * 2
14
I'm kind of running into the same situation numpy had, but I've
already established a workaround of checking the frequency content to
handle it.
Where did these freqs come from?
269 rfreqs15t = fftfreq(repetition_samples=15, complex=False)
They were generated as real-domain frequencies for 15 samples.
(Pdb) p np.fft.rfftfreq(15)
array([0. , 0.06666667, 0.13333333, 0.2 , 0.26666667,
0.33333333, 0.4 , 0.46666667])
(Pdb) p len(np.fft.rfftfreq(15))
8
Looking at their values, one can tell that ther is a DC offset, but
there is no nyquist frequency. So the total number would be 7*2
(conjugates) + 1 (dc offset) = 15. If there were a nyquist, it would
be + another 1.
Here:
118 if not has_dc_offset:
119 freq_count += 1
120 if not is_complex:
121 freq_count = (freq_count - 1) * 2
It appears to me to be incorrect to increment freq_count before
doubling it, as it is unaddressed as to whether there is a nyquist
frequency or not. The doubling in general is assuming both a nyquist
and a DC, but in odd frequency lengths there is never a nyquist, even
in the complex domain.
I tihnk I just copied that (n - 1) * 2 expression from the numpy docs
while handling different kinds of confusion.
0642
0643
no, my analysis above is wrong.
118 if not has_dc_offset:
119 freq_count += 1
here, it has a DC offset, so this code isn't hit. otherwise it would
be adjusting for the below to work
120 if not is_complex:
121 freq_count = (freq_count - 1) * 2
here, it ignores the DC offset, and doubles the remainder.
ok umm let me look at it again
8 - 1 = 7. 7 * 2 = 14.
then it would need to add 1 to bring back the DC offset, which it isn't doing.
whta if there's a nyquist?
8(nyquist+dc+6 conjugates) - 1 = 7. 7 * 2 = 14
(6*2+2 = 14)
it's correct if there's a nyquist.
if there is not a nyquist, it should be one more.
I did this, but I'm not confident that it's the right solution, I
don't fully understand it. stlil, it seemed like the immediate
bandaid:
if not is_complex:
has_nyquist = bool(freqs[-1] == 0.5)
freq_count = (freq_count - 1) * 2 + has_nyquist
0648
ok the code is now crashing with data that has a nyquist
ok i had inverted my condition. it took me some exploring in pdb to
identify this. 0651.
if not is_complex:
has_nyquist = bool(freqs[-1] == 0.5)
freq_count = (freq_count - 1) * 2 + (not has_nyquist)
That assertion is now passing. 0652.
0654 the new code is not making the right matrix for the longvec
situation where there are way more samples than frequencies. i feel
much more confident engaging this after the previous work, but i
wanted to note the frequencies.
(Pdb) p longspace_freqs
array([0. , 0.00270655, 0.0054131 , 0.00811965, 0.0108262 ,
0.01353275, 0.0162393 , 0.01894585, 0.0216524 ])
There is a DC offset, no immediately visible nyquist, and the
frequencies are real domain.
It's notable here, that since the upper frequency has changed, the
presence of 0.5 in the middle bin no longer indicates whether there
was an odd or even number of original frequencies :/ . Checking 0.5
was a quick hack to keep using the same code layout, but it was a
quick hack I was hoping worked.
It's notable that in this situation, the repetition_samples (output
sample count of the frequency matrix) is actually a floating point
value. So I am likely using freq_count in a way that isn't appropriate
here. I probably mean to use time_count .
0658
ValueError: cannot reshape array of size 28800 into shape (17,18)
Uncaught exception. Entering post mortem debugging
Running 'cont' or 'step' will restart the program
> /shared/src/scratch/fourier.py(167)create_freq2time()
-> ], axis=-1).reshape(len(complex_freqs), len(freqs)*2).T
(Pdb) list
162 interim_mat = complex_mat.copy()
163 interim_mat[head_idx:tail_idx] +=
interim_mat[len(freqs):][::-1].conj()
164 mat = np.stack([
165 interim_mat[:len(freqs)].real.T,
166 -interim_mat[:len(freqs)].imag.T
167 -> ], axis=-1).reshape(len(complex_freqs), len(freqs)*2).T
The output array has frequencies as columns (first dimension) and
samples as rows (second dimension). It looks like I'm using
len(complex_freqs) where I should be using time_count .
note:
- check for other (n-1)*2 uses and fix them
- check for uses of is_nyquist and consider them for fractional
repetition_samples or manual time_count
looks like that was the only (n-1)*2
i changed the matrix shape to use time_count
regarding is_nyquist, this will mean freq_count is wrong for subsignal
matrices with a nyquist. do they happen? when do they happen? i think
nyquists happen when there is an even number of complex frequency
bins?
if freq_count % 2 == 0:
if complex:
neg_freqs = np.linspace(-max_freq, -min_freq,
num=freq_count // 2, endpoint=True)
pos_freqs = -neg_freqs[:0:-1]
else:
pos_freqs = np.linspace(min_freq, max_freq, num=freq_count
// 2, endpoint=True)
neg_freqs = pos_freqs[:0]
else:
pos_freqs = np.linspace(min_freq, max_freq, num=freq_count //
2, endpoint=True)
neg_freqs = -pos_freqs[::-1] if complex else pos_freqs[:0]
Yes. When an even number of frequency bins was originally specified,
then a nyquist is included. Its value is max_freq, and it has no
complex conjugate.
Compare to when an odd number of frequency bins is originally
specified; now there is no nyquist and max_freq has a complex
conjugate.
When only real-domain frequencies are output, they are all linearly
spaced aside from the DC offset. So it looks like the information on
whether there is a nyquist or not is not within the list of
real-domain frequencies. [a fudgy way to include this data would be to
negate the nyquist, like it is in the complex frequencies]
a way to back negating the nyquist would be that a user generally
doesn't want less data from the same storage. given it is quite
reasonable to take the complex conjugate of the nyquist when it is not
actually the nyquist frequency but simply a middle frequency of
something else, the user may want to assume the larger amount of
frequency bins underlying the data, to extract more information on it
in more detail. does this make a difference ?? ... i think it means
that 1 more sample is needed to process a single repetition. ... i'm
not certain at this time what it means beyond that.
i guess for now i'll issue a warning or raise an exception in that case.
if not complex:
raise NotImplementedError("the time2freq functions do not
presently retain the information that an odd freq_count is associated
with real-domain data")
i raised the exception in this case because even freq_counts are more normative.
so, given the freq_count is even, this means there is indeed a
frequency with no complex conjugate, been calling that nyquist ... so
i can just set that value to True
0734 it looks like my code here:
if not complex:
freq_count = int((freq_count - 1) * 2)
actually usually prevents there from being an odd freq_count in the real domain
0738
trying to comprehend that there is a case where an odd freq_count happens.
when repetition_samples=15 and complex=False are passed to fftfreq,
it basically ends up calculating freq_count = repetition_samples .
this makes for an odd freq_count and real-domain data
it's kind of an obscurity of the test failing. it passes 15 to the
numpy functions too.
uhhh i'll just bring back the odd_freq_count bool parameter for now
0740
0744
ok, now normal assertions are failing again. also i added a fun
warning for when odd+noncomplex happens.
299 -> assert np.allclose(randvec2rfreq152randvec, randvec2rfreq152irfft)
293 randvec2rfreq152randvec = randvec2rfreq15.view(float) @ irft_to_15
294 randvec2rfreq152irfft = np.fft.irfft(randvec2rfreq15, 15)
287 rft_from_15 = create_time2freq(15, freqs=rfreqs15t,
odd_freq_count=True)
288 irft_to_15 = create_freq2time(15, freqs=rfreqs15t)
oops, gotta pass that parameter i made for the second use too.
0748 all the basic asserts are passing again, and additionally the
longvec matrices are created without crashing (the associated asserts
are still failing, and i've changed things they use so that's roughly
expected)
guess i'll attach my mess and poke at the real world a little
-------------- next part --------------
A non-text attachment was scrubbed...
Name: fourier.py
Type: text/x-python
Size: 17114 bytes
Desc: not available
URL: <https://lists.cpunks.org/pipermail/cypherpunks/attachments/20221118/6dacca5f/attachment-0001.py>
More information about the cypherpunks
mailing list