Quality assesment

[christianbrinch]

I have people emailing me with questions about what is wrong with LIME. I will pass this on to the community since I have had little to do with the recent updates. I have no idea of what is going on, but something is clearly badly broken.

Here is a continuum image of a very simple model, a box of constant density and temperature. This is what I get from LIME 1.41. The size is correct and the intensity is correct too:

The exact same model in version 1.7. Notice the factor of ~100 in peak intensity.

This is rather concerning. Has anyone actually bothered to test the code?

[christianbrinch]

Guys! Stuff needs to be tested before it is released. People expect to be able to do science with LIME and they tend to become angry if they have to write errata to their papers because of buggy codes.

Version 1.6.1 is a different kind of wrong:

Version 1.5 seems to be OK (for this simple test). We should tag 1.6 and 1.7 as buggy and let users know that they shouldn't use them.

[tlunttil]

Can you send the model that you used for this? I'll have a look.

For what it's worth, in my opinion ALL versions of LIME have major issues.

[christianbrinch]

I beg your pardon? You don't think that not being able to calculate the simplest model correctly is a significantly bigger issue? In my book, completely wrong results rank pretty high on the list of issue severity. Let me remind you that there is an analytic solution to the problem above and LIME 1.41 gets it right to within 1%.

The model I ran is pasted below.

`#include "lime.h" void input(inputPars par, image img){ par->radius = 2000AU; par->minScale = 0.1AU; par->dust = "jena_thin_e6.tab"; par->pIntensity = 60000; par->sinkPoints = 4000; par->sampling = 1;

img[0].freq = 267e9; img[0].pxls = 511; img[0].imgres = 0.1; img[0].theta = 0.0; img[0].distance = 100*PC; img[0].unit = 1; img[0].filename = "box1.6.1.fits"; }

void density(double x, double y, double z, double density){ if(x<500AU && x>-500AU && y<500AU && y>-500AU && z<500AU && z>-500AU){ density[0] = 1e41e6; } else { density[0] = 0.; } }

void temperature(double x, double y, double z, double temperature){ if(x<500AU && x>-500AU && y<500AU && y>-500AU && z<500AU && z>-500*AU){ temperature[0] = 60; } else { temperature[0] = 2.7; } }`

[tlunttil]

Take a look at #147. LIME results are all over the place.

[christianbrinch]

Seriously, could we maybe stick to the problem at hand?

And FYI, when I run a level population comparison with RATRAN I get this. Red points are LIME points and black points are RATRAN points. A K-S test will show you that both sets are drawn from the same distribution. I do not see what you call a major issue here. I am not saying that LIME prioer to version 1.6 has no problems, but I am saying that LIME 1.4 is passing the most basic tests. 1.6 and 1.7 is clearly not.

[allegroLeiden]

Pre-xmas is probably not the best time maybe to be delving deep into this but my take is as follows:

• I agree it is very necessary for 1.7 to be thoroughly tested with a variety of models.
• The code user community is best placed to do this, and this is exactly what they seem to be doing.

If users are discontented then I would say this is rather an issue of false expectations. If Lime had been my project from the start I would certainly not have released it when it was; I would be thinking of making a pre-release for beta testing about now. The genie is out of the bottle unfortunately and we all have to cope with that as best we can. Simply sticking with 1.4 is not an option.

By all means label 1.7 and 1.6.1 as buggy if you think that is a good idea but in fairness I think you then ought to label all the releases as buggy. That might dispel the (false) impression that Lime was polished, production software to which some wild fool has recently added a wodge of waffly code which has gone and broken everything.

Of far more concern to me is the apparently broken feedback loop. I can't fix bugs if I never get to hear about them. So far there have been a raft of bugs reported by students of Michiel Hogerheijde, all of which have been fixed; and #162, which I hope is now fixed with the 1.7 release. So if people mail you I would suggest you either get them to raise issues here, raise them yourself, or get users to mail me direct so I can do so.

[christianbrinch]

There are levels of buggyness! Some bugs are minor, some bugs are show stoppers. Take a look at my previous post. That is a 2D version of the model used in van Zadelhoff 2002 tested against the latest version of RATRAN. Please point out the show stopping error in that solution, because I don't see it. In fact, if you raytrace both solutions with either SKY or LIME you get identical spectra. The model at the top of this thread is constructed to be a 10 arcsec square with an intensity of 2e-7 jansky/pixel. That is exactly what LIME 1.4 gives you. Not a circle and not two orders of magnitude off. We may have different opinions, but in my view, having a buggy code that gives the correct result is by far better than a buggy code that gives wrong results. It is well established (ask the users!) that LIME 1.3/1.4 gives (mostly) correct results. Therefore, regression test should be made before any new release. New versions have to produce the same (or better) results than the previous version before they are publicly released; not worse results. LIME is a production code! People who download LIME from the official repository expect to get a working code. It is my name which is associated with the LIME code and even if I have had nothing to do with the latest releases, it is my credibility which is at stake. It has taken me many years to build a solid reputation around LIME and so far, I only (almost) had satisfied users. Last week I heard people refering to the "lame code". I am not amused!

[aschmiedeke]

Hi Christian, I just tried to run the model you posted further up to see if I can reproduce the behaviour. Unfortunately the model setup does not compile out of the box. Mostly due to missing multiplications. (On a sidenote, are you using 1e4 or 1e6 as the density? -->)

So I modified the model.c to read as given further below. Lime compiles, but then misses the velocity field setup. What setup are you using?

Edit: Ok, I see what was going on - the formatting here removes the * unless there are spaces.

include "lime.h"

void input(inputPars par, image img) { par->radius = 2000 AU; par->minScale = 0.1 AU; par->dust = "jena_thin_e6.tab"; par->pIntensity = 60000; par->sinkPoints = 4000; par->sampling = 1;

img[0].freq = 267e9;
img[0].pxls = 511;
img[0].imgres = 0.1;
img[0].theta = 0.0;
img[0].distance = 100*PC;
img[0].unit = 1;
img[0].filename = "box1.7.fits";

}

void density(double x, double y, double z, double density) { if(x < 500 AU && x > -500 AU && y < 500 AU && y > -500 AU && z < 500 AU && z > -500 AU) { density[0] = 1e4 1e6; } else { density[0] = 0.; } }

void temperature(double x, double y, double z, double temperature) { if(x < 500 AU && x > -500 AU && y < 500 AU && y > -500 AU && z < 500 AU && z > -500 * AU) { temperature[0] = 60; } else { temperature[0] = 2.7; } }

[christianbrinch]

Apparently github removed all multiplication signs. The density is 1e41e6, 1e6 being the conversion factor between cm^-3 -> m^-3. The version of LIME which I am using doesn't require the functions which is not used. This is continuum only, so abundance, doppler, and velocity is not used. If your version doesn't support this, just add void abundance(double x, double y, double z, double abundance){} and similarly for doppler and velocity.

[christianbrinch]

Michiel seems to have found the bug, so we can close this issue again when it is fixed.

[MichielHogerheijde]

Indeed, I did. The bug is that the CMB is added at the end of the raytracing, but not subtracted. If you are analysing continuum-subtracted line cubes, the results are ok. Otherwise the continuum has the CMB contribution still in there for part of the image (only where the raytracing was done).

This bug is present in versions 1.6 an 1.7 but not in 1.5. I've included the solution for src/raytrace.c below. I'm sure there is a more elegant way to communicate this via github, but I am not sufficiently versatile with github for that.

Solution: in src/raytrace.c from line 235, change:

/ Add or subtract cmb. / if(par->polarization){ / just add it to Stokes I /

ifdef FASTEXP

ray.intensity[stokesIi]+=FastExp(ray.tau[stokesIi])*local_cmb;

else

ray.intensity[stokesIi]+=exp(   -ray.tau[stokesIi])*local_cmb;

endif

}else{

ifdef FASTEXP

for(ichan=0;ichan<img[im].nchan;ichan++){
  ray.intensity[ichan]+=FastExp(ray.tau[ichan])*local_cmb;
}

else

for(ichan=0;ichan<img[im].nchan;ichan++){
  ray.intensity[ichan]+=exp(-ray.tau[ichan])*local_cmb;
}

endif

}

to:

/ Add or subtract cmb. / if(par->polarization){ / just add it to Stokes I /

ifdef FASTEXP

ray.intensity[stokesIi]+=(FastExp(ray.tau[stokesIi])-1.)*local_cmb;

else

ray.intensity[stokesIi]+=(exp(   -ray.tau[stokesIi])-1.)*local_cmb;

endif

}else{

ifdef FASTEXP

for(ichan=0;ichan<img[im].nchan;ichan++){
  ray.intensity[ichan]+=(FastExp(ray.tau[ichan])-1.)*local_cmb;
}

else

for(ichan=0;ichan<img[im].nchan;ichan++){
  ray.intensity[ichan]+=(exp(-ray.tau[ichan])-1.)*local_cmb;
}

endif

}

[tlunttil]

We actually discussed this with Ian, but we were not sure what the desired behaviour is (and then at least I kinda forgot about it). Should the subtraction/non-subtraction be an option for the user? It seems to me that at least for polarization stuff automatic CMB subtraction may not always be desired.

[tlunttil]

There's of course also a third option: including only the emission from the model, not the background. This doesn't make much sense when the background is CMB, but if at some point a non-isotropic background is implemented, this could be useful.

Technically all this would be trivial, it's just a matter of what people would find useful.

[MichielHogerheijde]

I cannot think of a situation where the flat CMB is relevant for LIME output. Any single-dish telescope will subtract this out, and interferometers can't see it. LIME does not have a model for the CMB polarization or any CMB intensity fluctuations. The only imprint that LIME can see is absorption by the source against the flat CMB. If the absorption is polarized and therefore what is left after CMB subtraction is a polarized signal, this will still be implemented correctly (it will have a negative intensity, just like the real measurements will have). If the intention was to leave in the CMB, it needs to be everywhere in the image, not just where the ray tracing is done. However, this suggests one could use LIME to model the kind of measurements a satelite link Planck does; I would not advice to use LIME for something like that, and think closing that option by hardwiring CMB subtraction is the safe thing to do.

Not adding/subtracting CMB is also not an option: there are real cases where molecular lines will absorb against the CMB leaving a negative intensity. The CMB is a given. Anyone who wants to experiment with it, should simply change the code themselves. We do not need to accomodate that.

Whether to include other varieties of uniform background intensities is a different matter, and should not be mixed up with a CMB background.

[tlunttil]

If the intention was to leave in the CMB, it needs to be everywhere in the image, not just where the ray tracing is done.

That reminds me of another question I had: what should be the default value for image pixels that are 'outside the model'? Currently they are set (hardwired) to value 1e-30, but wouldn't it be better to use NaN so that it would be more obvious that they do not correspond to anything in the model?

[MichielHogerheijde]

NaN is indistinguishable from "could not calculate" or "tried to divide by 0" or something else illegal. Setting it to 0 or 1e-30 is unambiguous. And since the model has an R_max, it is not undefined; it is empty space. I don't think users would appreciate finding NaNs there.

[christianbrinch]

Certainly not NANs. That would crash CASA or miriad. You want to be able to do arithmetic operations on the entire image. 0 is also no good, because then you cant ratio two images. Leave it as "sufficiently small".

[allegroLeiden]

Re the background: FITS has machinery for defining null values (which is not quite the same as NaN), and ideally I would suggest that these are the correct values to use outside the model boundary. Downstream programs ought to be able to cope with that. Obviously, if they cannot, then we should shrug our shoulders and resort to the otherwise very ugly expedient of writing 1e-30s. I think this ought to be checked out however.

[allegroLeiden]

Re code quality/reputation etc, there is a lot of heat being generated here and I don't propose to add to it. I will just state that I cannot agree that the quality of code should be judged by the numerical size of an error it makes on a single test case. I will further reiterate that in my view, whatever the version number, Lime is now essentially in a beta test phase: i.e. that all the heavy lifting of code fixing (involving 11k lines of new or edited code from myself alone) has now pretty much been done, and we can turn to the job of bringing the results in line with expectations. For this to be an efficient process we need users testing Lime with a variety of models and reporting any bugs to the development team.