search

SyneRBI / SIRF

Main repository for the CCP SynerBI software
http://www.ccpsynerbi.ac.uk
Other
58 stars 29 forks source link

Memory leak in PET_MCIR.py #593

Open rijobro opened 4 years ago

rijobro commented 4 years ago

There is/are error leak(s) in examples/Python/PETMR/PET_MCIR.py.

I ran it through valgrind like this:

valgrind --tool=memcheck --leak-check=full --show-leak-kinds=all \
    python PET_MCIR.py \
        --trans "mfs_i_nii/mf_i_g1.nii" \
        --sinos "EM_g1.hs" \
        --outp mcir_output/recon 
        --trans_type-t disp \
        --iter 1 \
        --save_interval 1 \
        --verbosity 10
        --nifti

The valgrind output is attached, I can't make head or tail of it.

I also tried memory profiler and the output is more concise. Their example shows memory increasing during allocation, and then decreasing when objects go out of scope/destructors are called. My output doesn't show any decrease in memory. I wonder if it means that everything is leaking memory (please don't be this!) or that the profiler isn't smart enough to notice objects are getting deleted on the C-level. Output also attached.

Any ideas? I suppose one thing I need to do is break the code into small chunks and verify each bit, but this will take a long time...

Log files

valgrind_output.txt memprofiler_output.txt

rijobro commented 4 years ago

Also, how should I interpret the end of the valgrind output, with 20GB still reachable?

==10926== LEAK SUMMARY:
==10926==    definitely lost: 1,608 bytes in 12 blocks
==10926==    indirectly lost: 16,460 bytes in 41 blocks
==10926==      possibly lost: 895,780 bytes in 1,615 blocks
==10926==    still reachable: 20,810,150 bytes in 24,290 blocks
==10926==                       of which reachable via heuristic:
==10926==                         stdstring          : 5,374 bytes in 107 blocks
==10926==                         multipleinheritance: 4,448 bytes in 5 blocks
==10926==         suppressed: 0 bytes in 0 blocks
rijobro commented 4 years ago

Lastly, it's also possible that this is a CIL issue as PET_MCIR.py uses both SIRF and CIL. @paskino

rijobro commented 4 years ago

End of the output from valgrind on osem_reconstruction.py (trying to test a smaller component) looks similar, but still undecipherable:

==12421== LEAK SUMMARY:
==12421==    definitely lost: 976 bytes in 5 blocks
==12421==    indirectly lost: 20,172 bytes in 46 blocks
==12421==      possibly lost: 324,318 bytes in 223 blocks
==12421==    still reachable: 13,087,070 bytes in 6,847 blocks
==12421==         suppressed: 0 bytes in 0 blocks
==12421== 
==12421== For counts of detected and suppressed errors, rerun with: -v
==12421== Use --track-origins=yes to see where uninitialised values come from
==12421== ERROR SUMMARY: 4634 errors from 186 contexts (suppressed: 0 from 0)
rijobro commented 4 years ago

@paskino disabled numba on your recommendation, valgrind output attached. Less data "possibly lost". Also no mentions of llvm this time, but I wonder if these were just warnings as opposed to areas of memory leaking? Honestly I've no idea...

valgrind_memcheck_new.txt

KrisThielemans commented 4 years ago

Most of these are Python. Maybe there'd be less of them with a more recent Python (but maybe not). Looking for sirf I found

These re-occur a few times of course. They don't seem to be dramatic though, so not sure if this will help, but I did give up about half-way though the log file. Fixing these small leaks will be good, and also clean-up the log of course.

rijobro commented 4 years ago

I tried simplifying as much as possible to find the leak, so I removed the motion component to simply do at PET reconstruction with PDHG (no corrections with attenuation, randoms, norm, etc.).

Caused my Mac to run out of memory so the leak is still present in this code, will run it through valgrind (on linux as not available on OSX).

The bulk of the code is below, the main objects that could be at fault: AcquisitionData, sirf.STIR.ImageData, AcquisitionModelUsingRayTracingMatrix, KullbackLeibler, BlockFunction, BlockOperator, PDHG, NiftiImageData.

Question for CCPi guys (@paskino, @gfardell, @epapoutsellis): Can I simplify the code further? Do I need BlockFunction and BlockOperator if I only have one acquisition model?

sino = pet.AcquisitionData(sino_file)
sino = make_sino_positive(sino)
image = sino.create_uniform_image(1.0, nxny)

print("Setting up acquisition model...")

acq_model = pet.AcquisitionModelUsingRayTracingMatrix()
acq_model.set_up(sino, image)

print("Setting up reconstructor...")

# Configure the PDHG algorithm
kl = KullbackLeibler(b=sino, eta=(sino * 0 + 1e-5))
f = BlockFunction(kl)
K = BlockOperator(acq_model)
normK = K.norm(iterations=10)

# normK = LinearOperator.PowerMethod(K, iterations=10)[0]
# default values
sigma = 1/normK
tau = 1/normK 
sigma = 0.001
tau = 1/(sigma*normK**2)
print("Norm of the BlockOperator ", normK)

# No regularisation, only positivity constraints
G = IndicatorBox(lower=0)

print("Creating up reconstructor...")

pdhg = PDHG(f=f, g=G, operator=K, sigma=sigma, tau=tau,
            max_iteration=1000,
            update_objective_interval=1)

for i in range(1,num_iters+1):
    print("Running iteration " + str(i) + "...")
    pdhg.run(1, verbose=True)
    reg.NiftiImageData(pdhg.get_output()).write(outp_prefix + "_iters" + str(i))