Assertion failure when using larger windows with more points

[hn-88]

I wanted to model ps pumped SCG over many metres of propagation. Thinking that running out of memory would be the bottleneck, I implemented an iterative script, saving the pulse-out to file after short lengths, and reading from it in the next iteration. Attached
SCG-Head-by-parts2.py.txt

But I saw that the temporal pulse goes outside the 20 ps window when propagated for 20 metres, as in SCG-Head-by-partsNL.py.txt Wanting to increase the time window, if I increase the time window by a factor of two and increase the number of points also by a factor of two, I start getting assertion failure errors.

Time taken so far - 0 hours 0 minutes 30 seconds
30.3910000324
Iteration 2 of 2000
Pulse energy before NL-1050-NEG-1.txt : 11.9999464432 nJ
Traceback (most recent call last):
  File "SCG-Head-by-partsNL2.py", line 80, in <module>
    y, AW, AT, pulse_out = evol.propagate(pulse_in=pulse, fiber=fiber1, n_steps=
Steps)
  File "D:\Anaconda\lib\site-packages\pynlo\interactions\FourWaveMixing\SSFM.py"
, line 485, in propagate
    self.setup_fftw(pulse_in, fiber, output_power)
  File "D:\Anaconda\lib\site-packages\pynlo\interactions\FourWaveMixing\SSFM.py"
, line 225, in setup_fftw
    self.A[:]       = self.conditional_fftshift(pulse_in.AT, verify=True)
  File "D:\Anaconda\lib\site-packages\pynlo\interactions\FourWaveMixing\SSFM.py"
, line 672, in conditional_fftshift
    assert abs(chksum - np.sum(abs(x))) <= np.finfo(float).eps
AssertionError 

Can you help me avoid this problem? The file which gave this error is SCG-Head-by-partsNL2.py.txt and it uses the dispersion from NL-1050-NEG-1.txt

[DanHickstein]

@ycasg can probably give us a better idea of what is happening, but in general, this assertion is to check that the pulse energy is maintained after a Fourier transform. I have received a similar error when trying to model ps pulses with added noise, and I just turned the assertion off, which is not ideal.

My guess is that there are small numerical round-off errors that mean that the energy is not precisely maintained after the FFT. The fact that there can be small numerical differences should be taken care of by the fact that we are comparing to np.finfo(float).eps, which is basically the numerical precision of the machine.

However, since we are taking the sum of np.abs(x), then I think that the error is much higher than .eps. On my system, the machine precision (np.finfo(float).eps) is 2.22044604925e-16.

For the first propagation with your script, the error (abs(chksum - np.sum(abs(x)))) is zero. However, on the second propagation, with the "messy" pulse, the error is 2.1827872842550278e-11, which is still very small compared to np.sum(np.abs(x)), which is 126082.16....

If I multiple the machine precision by the length of x, I get np.finfo(float).eps*x.shape[0] = 1.4551915228366852e-11, which is close to 1/2 of the observed error. I'm guessing that the possible error is doubled because x is complex-valued. I'm a little confused why the observed error seems so close to the maximum possible value, but numerical errors are weird, I guess.

In any case, I think that there is a solution, the offending code (line 673 in SSFM.py) should change to acknowledge that the errors can add. So, the previous line

assert abs(chksum - np.sum(abs(x))) <= np.finfo(float).eps

should change to:

assert abs(chksum - np.sum(abs(x))) <= np.finfo(float).eps*x.shape[0]*2

I made a PR for this: https://github.com/pyNLO/PyNLO/pull/31

But, for now, you can just change the line in your own version of SSFM.py.

Also, I'm not sure if you need to run multiple propagation steps. I have never run into an issue with memory.

[hn-88]

Thanks @DanHickstein - will try again with the change to SSFM.py

Regarding memory, my machine has only 2 GB RAM and 4 GB swap, probably that's why I get memory errors and you don't :)

Trying 10 metres propagation - 10,000 steps directly with 2^17 points,

$ python SCG-Head-17.py
Traceback (most recent call last):
  File "SCG-Head-17.py", line 74, in <module>
    y, AW, AT, pulse_out = evol.propagate(pulse_in=pulse, fiber=fiber1, n_steps=Steps)
  File "/home/hari/anaconda2/lib/python2.7/site-packages/pynlo/interactions/FourWaveMixing/SSFM.py", line 477, in propagate
    AW = np.complex64(np.zeros((pulse_in.NPTS, n_steps)))
MemoryError

[hn-88]

Nope, even with the change suggested by @DanHickstein - assertion failure.

Step: 499 Distance remaining: 0.001
Pulse energy after: 11.9843154261 nJ
Pulse energy after NL-1050-NEG-1.txt : 11.9891360822 nJ
Iteration 1 of 20
Pulse energy before NL-1050-NEG-1.txt : 11.9891360822 nJ
Traceback (most recent call last):
  File "SCG-Head-by-partsNL2.py", line 80, in <module>
    y, AW, AT, pulse_out = evol.propagate(pulse_in=pulse, fiber=fiber1, n_steps=Steps)
  File "D:\Anaconda\lib\site-packages\pynlo\interactions\FourWaveMixing\SSFM.py"
, line 485, in propagate
    self.setup_fftw(pulse_in, fiber, output_power)
  File "D:\Anaconda\lib\site-packages\pynlo\interactions\FourWaveMixing\SSFM.py"
, line 225, in setup_fftw
    self.A[:]       = self.conditional_fftshift(pulse_in.AT, verify=True)
  File "D:\Anaconda\lib\site-packages\pynlo\interactions\FourWaveMixing\SSFM.py"
, line 672, in conditional_fftshift
    assert abs(chksum - np.sum(abs(x))) <= np.finfo(float).eps*x.shape[0]*2
AssertionError

Slightly modified script below, SCG-Head-by-partsNL2.py.txt

[DanHickstein]

Hi @hn-88,

Interestingly, your script runs fine on my laptop, even with the master version of PyNLO (without the changes). It's possible that the numerical error is platform dependent.

You can always just disable this check temporarily (comment out that line), but I am curious what's going on.

Can you hack SSFM.py to include a print statement that compares the values of abs(chksum - np.sum(abs(x))) and np.finfo(float).eps*x.shape[0]?

Edit: and yes, I have 16 GB of ram, so I might not notice if the calculation took up several GB ;). The problem is likely the 2D arrays that save all of the frequency/time values at every step. These are basically huge images. With 2^17 points and 10k steps, you have over 1e9 points. Multiply that by the number of bytes per complex128 number (16 bytes) and even I am out of memory.

Why take so many steps? The beauty of this algorithm is that it has an adaptive step size: it will automatically detect when there are changes to the spectrum rapidly taking place and change to a very small step size to correctly model the physics. Then, when there is very little happening (typically in the later stages of propagation, after soliton fission) the integrator will use a larger step size and run faster.

So, do not worry about using a lot of steps to ensure a good result. Simply request the number of steps that you wish to display (~200 usually makes a nice plot) and let the integrator take care of situations where finer steps are needed.

[hn-88]

Thanks @DanHickstein - I did not realize that the step size is only for display - will check with smaller step sizes. My problem with assert fails may be due to precision loss during the save/read - this link talks of precision loss while printing, http://stackoverflow.com/questions/12956333/printing-numpy-float64-with-full-precision

[hn-88]

With print statements before the Assert,

Iteration 1 of 20
Pulse energy before NL-1050-NEG-1.txt : 11.9891360822 nJ
abs(chksum - np.sum(abs(x))) is
2.91038304567e-11
np.finfo(float).eps*x.shape[0]*2 is
1.45519152284e-11
Traceback (most recent call last):
  File "SCG-Head-by-partsNL2.py", line 80, in <module>
    y, AW, AT, pulse_out = evol.propagate(pulse_in=pulse, fiber=fiber1, n_steps=Steps)
  File "/home/hari/anaconda2/lib/python2.7/site-packages/pynlo/interactions/FourWaveMixing/SSFM.py", line 485, in propagate
    self.setup_fftw(pulse_in, fiber, output_power)
  File "/home/hari/anaconda2/lib/python2.7/site-packages/pynlo/interactions/FourWaveMixing/SSFM.py", line 225, in setup_fftw
    self.A[:]       = self.conditional_fftshift(pulse_in.AT, verify=True)
  File "/home/hari/anaconda2/lib/python2.7/site-packages/pynlo/interactions/FourWaveMixing/SSFM.py", line 676, in conditional_fftshift
    assert abs(chksum - np.sum(abs(x))) <= np.finfo(float).eps*x.shape[0]*2
AssertionError

[hn-88]

Just to check, I ran the same thing on Linux - Lubuntu 16.04 64 bit - and Windows XP 32 bit - running on the same machine. The linux python crashed as above, at the Iteration 1 propagate call, after loading the pulse saved by Iteration 0.

Interestingly, on Windows, this crash did not happen. Iteration 1 continued until the point of saving the pulse, when it ran out of Memory - I watched with task manager as Windows increased the memory allocated from 600 to 800 MB, and then gave up :)

Step: 499 Distance remaining: 0.001
Pulse energy after: 11.983763174 nJ
Pulse energy after NL-1050-NEG-1.txt : 11.9863579478 nJ
Traceback (most recent call last):
  File "SCG-Head-by-partsNL2.py", line 87, in <module>
    pickle.dump(AW, fid)
  File "D:\Anaconda\lib\pickle.py", line 1370, in dump
    Pickler(file, protocol).dump(obj)
  File "D:\Anaconda\lib\pickle.py", line 224, in dump
    self.save(obj)
  File "D:\Anaconda\lib\pickle.py", line 331, in save
    self.save_reduce(obj=obj, *rv)
  File "D:\Anaconda\lib\pickle.py", line 419, in save_reduce
    save(state)
  File "D:\Anaconda\lib\pickle.py", line 286, in save
    f(self, obj) # Call unbound method with explicit self
  File "D:\Anaconda\lib\pickle.py", line 562, in save_tuple
    save(element)
  File "D:\Anaconda\lib\pickle.py", line 286, in save
    f(self, obj) # Call unbound method with explicit self
  File "D:\Anaconda\lib\pickle.py", line 488, in save_string
    self.write(STRING + repr(obj) + '\n')
MemoryError

On Windows also the eps value seems to be the same. But during iteration 1, the result on Windows was

abs(chksum - np.sum(abs(x))) is
0

Perhaps on Linux the file is being read before it is completely written? Anyway, now I will try running without the saves etc, with just fewer display steps.

[hn-88]

Just for reference, I also ran on Windows with fewer number of steps and more save/load iterations. On some iterations, the error is exactly equal to the eps value! And on some, zero!

[DanHickstein]

Hmmm, so you found a situation where abs(chksum - np.sum(abs(x))) is almost exactly 4x np.finfo(float).eps*x.shape[0]? I'm not sure that I understand this - is it related to saving and loading the pulses?

Still, it seems like the error remains on same order-of-magnitude as the numerical noise.

I don't want to make the tolerance for errors in the FFT too high, but I don't think that it should be a problem to increase the tolerance by a few times.

Alternatively, I think that the difference becomes always zero if we change np.sum to math.fsum (at least this was the case in my limited testing) https://docs.python.org/2/library/math.html#math.fsum https://github.com/numpy/numpy/issues/2448 I think that the only downside is that it increases the time for the sum from 200 us to 2 ms. I'll make a PR for that change.

[hn-88]

Some more interesting info.

Since the plot with 2^15 points and a 40 fs window was not wide enough - the time-based plot was going out wider than 40 ps in the 10 m of fiber - I tried with 2^18 points and a window of 160 ps. Interestingly, the time taken was almost the same as for the 2^15 plot run - 2 hrs 45 min or so. But it gave a warn mesg at the beginning:

D:\Anaconda\Saved data\SCG-Head-NL4>python NL4.py
Pulse energy before NL-1050-NEG-1.txt : 12.0 nJ
D:\Anaconda\lib\site-packages\pynlo\interactions\FourWaveMixing\SSFM.py:271: Run
timeWarning: overflow encountered in exp
  np.exp(-T/tau2)*np.sin(T/tau1)
D:\Anaconda\lib\site-packages\pynlo\interactions\FourWaveMixing\SSFM.py:271: Run
timeWarning: overflow encountered in multiply
  np.exp(-T/tau2)*np.sin(T/tau1)
abs(chksum - np.sum(abs(x))) =
2.9103830456733704e-11
np.finfo(float).eps*x.shape[0]*2 =
5.8207660913467407e-11
pulse energy in  126079.778504
copied as   126079.778504
Step: 0 Distance remaining: 10.0

The results were: which seem to indicate that a good SC had been generated by 20 cm itself, without the need to propagate 10 metres as in the paper http://www.orc.soton.ac.uk/publications/65xx/6571.pdf Will try that tomorrow.

[DanHickstein]

I think that those warnings might be associated with np.exp(-T/tau2) getting really close to zero with high values of T. Basically, I think numpy is trying to warn you that you requested a really small number and the number that you are getting might not be distinct from zero. Just a guess...

It looks like you are starting with a >5 ps pulse. Once this pulse has spread to several times it's original temporal width, most of the nonlinear processes will be over and it's just broadening temporally due to the dispersion of the fiber. So, in this case, it looks like there is no point in simulating more than a few cm of fiber.

You may be able to use fewer points. Just make sure that the simulation results are converged to your satisfaction with the number of points.

[hn-88]

@DanHickstein , the pulsewidth is 400 fs at the beginning. In the paper http://www.orc.soton.ac.uk/publications/65xx/6571.pdf , results after 10 m of propagation through the NL-1050-NEG-1 fiber are reported. I was trying to re-create that. As noted above, it looks like most of the SCG is done in the first 50 cm of the fiber, results for 50 cm below. With 2^15 points, 40 ps window. These were obtained without the multiple iterations - only saving at the end of the run instead of displaying immediately. Displayed using a separate script. From the time plot, I see that I can reduce the window further, anyway my purpose - recreating the paper - is served :)

[DanHickstein]

Did you obtain satisfactory agreement with the paper? It looks like you predict a broader spectrum with additional features, where the paper just shows the pure SPM in the center. Potentially they have some pedestal of uncompressed light that makes their "short pulse power" lower than what is in the manuscript. Alternatively, they may not have been perfectly compressing their pulse at the entrance to the HNLF fiber. In trying to reproduce the experimental data, I would consider most published pulse durations and powers to be optimistic estimates, and adjust your simulations accordingly.

The fact that they used 10 meters of fiber when 1 meter may have sufficed is probably a result of the HNLF fiber being expensive, and they did not want to cut it into pieces. :)

[hn-88]

Yes @DanHickstein , you're right. Also, the resolution of their spectrometer was only 5 nm. So, the oscillations at the centre of the spectrum have been smoothed out - we can use smoothing functions to get similar results.

I did these as a preliminary test, to simulate what sort of lengths we would need if we pump a commercially available PCF with the lasers we have available - 10 ns Nd:YAG which can give >200 mJ, and a 532nm pumped tunable Ti-sapphire, giving 10ns pulses tunable 700 to 950 nm. In both cases, we would be limited by the damage threshold of the fiber to pumping < 50 kW peak power.

I tried with 50 cm, not much spectral broadening with SC-5.0-1040 fiber. That simulation run took less than half an hour. Now I'm trying 2 metres propagation, looks like the simulation is going to take 50+ hours! I'm simulating with a 3 ps pulse, with 50 kW peak power, assuming 10 ns with the same peak power will have similar effects. 10^18 points, 160 ps window.

[ycasg]

With those pulse parameters, you are in a sort of "quasi CW regime." The fact that the simulation becomes extremely slow suggests to me that you are seeing modulation instability -- the spontaneous and stochastic formation of many solitons. This would be seen most clearly in the time domain, where the electric field would become spikey. Simulating this gets tricky, although it should hopefully be accurate.

Out of curiosity, why are you looking at SCG with long pulses, as opposed to something shorter?

[hn-88]

Those are the lasers which we have currently in our lab :) A femtosecond laser will be coming, but only next year or maybe later.

[DanHickstein]

Yes, simulations that suddenly slow tend to be a sign that modulation instability or other chaotic behavior is occurring. If the light does not need to be coherent (or overly stable) then operating in the modulation instability regime might be okay...