I dramatically clean up the code and rewrite the example.py script. Now the code uses a more concise argument passer to pass arguments for simulations, which was using repeated arguments before. Most importantly, Source Extractor now can run in parallel by passing ncpu to the method, which dramatically speeds up the code.
Caveat: I am still using Gaussian and Exponential as the defaul light profiles for simulating galaxies. This is because that I could not find the cause why the SE-derived MAG_AUTO has an offset to input mag_sersic (i.e., I could not use SE to recover the input sersic magnitude). It could be due to background subtraction? Let me know if you have some insights.