A glitch in Jorrie’s Cosmo-Calculator?

[JimJCW]

Here is a patched version of Lightcone7, dated 2002-05-14, which numerically corrects for the offset in the low-z regime.

First, the output z range is still different from the input z range:

I did the following experiment:

If I calculate D_now using Gnedin's and your 2022 calculator and plot 'D_Gnedin - D_Jorrie' vs z, I expect to see points fluctuating around the D_nedin - D_Jorrie = 0 line. But I got the following result:

I wonder what clues we can get out of this.

 

[Ibix]

@JimJCW - You asked: Please let me know the line number and file name concerning the '1.001 + z' problem

In my python version there's a line near the bottom with a comment saying something like "to replicate a bug in Lightcone7 uncomment this line". In the Javascript version, no idea. I noticed it in the output - if you turn on both the $S$ and $z$ columns and bump up the number of dps they display with you can verify it yourself by calculating $S-z$, but I don't know why it happens.

 

[Ibix]

I have used Tamara Davis' "Fundamental Aspects of the Expansion of the The universe and Cosmic Horizons" equations in the Tutorial, but the implementation in Lightcone7 was a bit different and quite possibly inaccurate.

Note that $(D_{par}+D_{hor})S=\int_0^\infty T_Hds=\mathrm{const}$ only holds for comoving distances and not for the proper distances used in Lightcone7, where $D_{par}$ tends to infinity and $D_{hor}$ tends to a finite constant value. This is clear from the Davis Three-panel Cosmological Expansion graphs in my sig.

Right. I'll do a bit more reading.

 

[Jorrie]

First, the output z range is still different from the input z range:

View attachment 301506
View attachment 301507

Hi Jim, yes there is that difference, but that's not the problem, the inputs range must just be a little larger than the output range. As long as the distance now output is correct for the displayed z, that's acceptable for now.

For a trial, I've just subtracted 0.001 (line 274 in the 2022-05-14 version) from the output z in order to cancel the origin's horizontal offset, which is caused by the way the distances are computed in the legacy algorithm. The question is: are those (new) values correct? Perhaps you can check it against Gnedin's.

I know that this is not the right way to handle it, so we should be investigating a purer algorithm, like @Ibix is looking at. My own time is slightly limited atm, but I will keep watching and contributing as it hopefully develops as a community project.

Thanks for contributing!

 

[pbuk]

Still pondering an overhaul of the legacy-ridden calculation.js module.

That might be an idea. Is the duplication of the line result.H_t = result.H_t*Ynow; in ScaleResultsan error?

      case "Normalized":
          result.Tnow = result.Tnow/Tage;
          result.Y = result.Y/Ynow;
          result.Dnow  = result.Dnow/Ynow;
          result.Dthen = result.Dthen/Ynow;
          result.Dhor = result.Dhor/Ynow;
          result.Dpar = result.Dpar/Ynow;
          result.H_t = result.H_t*Ynow;
          result.H_t = result.H_t*Ynow;
          result.TemperatureT = result.TemperatureT / 2.725 ;  
        break;

 

[pbuk]

You could put the code to github? I expect github won't last forever (nothing does), but it has a huge amount of inertia. Your code wouldn't run there, but it would always be available for someone to host should your version disappear.

GitHub provides a (free for open source) service called GitHub Pages which can host this so that it WILL run there (as long as it doesn't go over the free usage limits which should be adequate).

Now GitHub is owned by Microsoft, and they have moved a lot of their code (both open source and proprietory) on to it, it is not going to disappear for a LONG time.

 

[Vanadium 50]

they [NMicrosoft] have moved a lot of their code (both open source and proprietory) on to it, it is not going to disappear for a LONG time.

As much as we might want it to. <sigh>

 

[Jorrie]

That might be an idea. Is the duplication of the line result.H_t = result.H_t*Ynow; in ScaleResultsan error?

      case "Normalized":
          result.Tnow = result.Tnow/Tage;
          result.Y = result.Y/Ynow;
          result.Dnow  = result.Dnow/Ynow;
          result.Dthen = result.Dthen/Ynow;
          result.Dhor = result.Dhor/Ynow;
          result.Dpar = result.Dpar/Ynow;
          result.H_t = result.H_t*Ynow;
          result.H_t = result.H_t*Ynow;
          result.TemperatureT = result.TemperatureT / 2.725 ; 
        break;

Yes, good spot! Case "Normalized" is not often used, so that's maybe why the mistake escaped detection.

 

[JimJCW]

@JimJCW - You asked: Please let me know the line number and file name concerning the '1.001 + z' problem

In my python version there's a line near the bottom with a comment saying something like "to replicate a bug in Lightcone7 uncomment this line". In the Javascript version, no idea. I noticed it in the output - if you turn on both the $S$ and $z$ columns and bump up the number of dps they display with you can verify it yourself by calculating $S-z$, but I don't know why it happens.

Would you please calculate the following D_now values using your calculator so I can do some comparison?

z​	Gnedin​	Jorrie, 2022​	ibix​
0.1	1.41004	1.4084
0.09	1.27214	1.270609
0.08	1.13354	1.132106
0.07	0.99426	0.992931
0.06	0.85427	0.85305
0.05	0.7136	0.712477
0.04	0.57225	0.571225
0.03	0.43023	0.429278
0.02	0.28748	0.286651
0.01	0.14408	0.14336
0	0	0.000621
-0.01	0.14474	0.14526
-0.02	0.29015	0.290567
-0.03	0.43623	0.436536
-0.04	0.58295	0.583161
-0.05	0.73032	0.730433
-0.06	0.87834	0.878358
-0.07	1.02699	1.026926
-0.08	1.17631	1.176124
-0.09	1.32623	1.325968
-0.1	1.47679	1.476416

I am puzzled by the plot shown in Post #36.

 

[Jorrie]

I think the puzzles will go away if we get the distances algorithm sorted out.

 

[Jorrie]

GitHub provides a (free for open source) service called GitHub Pages which can host this so that it WILL run there (as long as it doesn't go over the free usage limits which should be adequate).

I will give GitHub Pages a try, because it seems a good place for a collaborative project.

 

[Ibix]

Would you please calculate the following D_now values using your calculator so I can do some comparison?

Install python, man! If you're remotely interested in any scientific computing it's a useful tool. I believe Anaconda is the recommended Windows distribution, and I believe it comes with the scipy library I used in my code. Then you can play around with my code to your heart's content.

Here's the results. $S$ is set to the correct $1+z$ not the compatibility $1.001+z$ mode.

z​	Dnow Gly​
0.1​	1.41​
0.09​	1.27​
0.08​	1.13​
0.07​	0.994​
0.06​	0.854​
0.05​	0.713​
0.04​	0.572​
0.03​	0.43​
0.02​	0.287​
0.01​	0.144​
0​	0​
-0.01​	0.145​
-0.02​	0.29​
-0.03​	0.436​
-0.04​	0.583​
-0.05​	0.73​
-0.06​	0.878​
-0.07​	1.03​
-0.08​	1.18​
-0.09​	1.33​
-0.1​	1.48​

 

[pbuk]

I think the puzzles will go away if we get the distances algorithm sorted out.

Yes, although this is not going to be as simple as porting @Ibix's code which leverages a scipy module (which in turn leverages a FORTRAN module) to do integration with infinite limits and automatic step size adjustment. This module is not available in the browser so it is necessary to implement a similar algorithm "by hand".

I have some experience of collaborative Javascript projects, @Jorrie do you mind if I set something up as a illustration of how to transition the calculation module to a robust modern approach?

 

[Ibix]

I think the puzzles will go away if we get the distances algorithm sorted out.

I've been having a look at your Calculation.js module. Your $D_{hor}$ is just $1/H(z)$ multiplied by the appropriate unit conversion factor, so is actually the same as your $R$ column - I can do that.

However your $D_{par}$ seems to be $a\int_0^a\frac{da'}{a'H(a')}$ which, noting that $a=1/S$, $a'=1/s$ and $da'=-ds/s^2$, is $\frac 1S\int_S^\infty\frac{ds}{sH(s)}$, which is what I've implemented, so I'm a bit confused. Also full of viruses today, so possible I've screwed something up.

 

[Ibix]

This module is not available in the browser so it is necessary to implement a similar algorithm "by hand".

This has an ODE solver: https://ccc-js.github.io/numeric2/index.html. I think that could be slightly abused to do integration. Or we could just use the approach in the existing module, which seems to work fine with the possible exception of the $D_{par}$ calculation that we need to revisit. I don't think it's a bad approach and seems pretty responsive (although I'd move some stuff around to make it a bit easier to see what's going on, and you can get rid of one pass of the integrator, I think). It's just grown messy and needs a dead code pruning round.

 

[pbuk]

This has an ODE solver: https://ccc-js.github.io/numeric2/index.html. I think that could be slightly abused to do integration.

I don't think that is a very good idea: (i) the project is even more stale than Lightcone7 (no commits for 10 years); (ii) robust methods for quadrature are very different from methods for initial value problems.

Or we could just use the approach in the existing module.

I think that would be better.

I don't think it's a bad approach and seems pretty responsive (although I'd move some stuff around to make it a bit easier to see what's going on, and you can get rid of one pass of the integrator, I think). It's just grown messy and needs a dead code pruning round.

I agree with all of that. An approach I have found successful in the past is to extract the existing code into a stand-alone module and implement continuous integration with some high level functional tests to check you don't break anything before refactoring the existing code so that you can figure out how it works/why it doesn't.

Once you have done this you can prioritize what needs to be done next (quick fix to address bugs? rewrite algorithm for better stability/accuracy/speed? Etc...).

 

[PAllen]

(ii) robust methods for quadrature are very different from methods for initial value problems.
.

While I have no comment on the referenced project, far more development of adaptive and variable step size has occurred for ODEs than quadrature. In particular, for spikey functions poorly approximated by polynomials, a given accuracy is often achieved with much less computation by casting the quadrature as an ODE, and using modern ODE methods.

 

[pbuk]

Here it is https://lightcone7.github.io/LightCone7.html with the calculations split out into a separate module (source here) with the bare minimum of changes to get it to work independently.

Needs some tidying up and documentation: PM me with github user details if you want to use/take over this.

 

[pbuk]

In particular, for spikey functions poorly approximated by polynomials, a given accuracy is often achieved with much less computation by casting the quadrature as an ODE, and using modern ODE methods.

I didn't realize that, although these functions are generally smooth.

Far more development of adaptive and variable step size has occurred for ODEs than quadrature.

Yes this is true in general, but there is not much implemented in Javascript for ODEs either.

 

[PAllen]

I didn't realize that, although these functions are generally smooth.

Agreed.

Yes this is true in general, but there is not much implemented in Javascript for ODEs either.

Probably true, but in other languages, publicly available implementations of most any method should be available. The SciPy quadrature interface @Ibix mentioned looks good. It uses adaptive step size. Scipy also has ODE methods.

[edit: I see this was all discussed before - browser availability is a requirement. I have no idea what is available open source in this context.]

 

[Jorrie]

I have some experience of collaborative Javascript projects, @Jorrie do you mind if I set something up as a illustration of how to transition the calculation module to a robust modern approach?

By all means, do that.

 

[JimJCW]

Install python, man! . . .

That sounds like a good idea. Thanks for calculating the D_now values! A comparison is shown below:

z	Gnedin	Jorrie, 2022	ibix
0.1	1.41004	1.4084	1.41
0.09	1.27214	1.270609	1.27
0.08	1.13354	1.132106	1.13
0.07	0.99426	0.992931	0.994
0.06	0.85427	0.85305	0.854
0.05	0.7136	0.712477	0.713
0.04	0.57225	0.571225	0.572
0.03	0.43023	0.429278	0.43
0.02	0.28748	0.286651	0.28
0.01	0.14408	0.14336	0.144
0	0	0.000621	0
-0.01	0.14474	0.14526	0.145
-0.02	0.29015	0.290567	0.29
-0.03	0.43623	0.436536	0.436
-0.04	0.58295	0.583161	0.583
-0.05	0.73032	0.730433	0.73
-0.06	0.87834	0.878358	0.878
-0.07	1.02699	1.026926	1.03
-0.08	1.17631	1.176124	1.18
-0.09	1.32623	1.325968	1.33
-0.1	1.47679	1.476416	1.48

If we plot the three D_now vs z curves, they are just about indistinguishable. If we add the (D_Gnedin – D_ibix) result to the plot shown in Post #36, we get the following:

In general, (D_Gnedin - D_ibix) fluctuates around the y = 0 line. The greater deviations from the line near the two ends are probably due to round-off effect of the calculated D_ibix. I am still puzzled by the slightly different behaviors of the two data sets.

 

[pbuk]

I am still puzzled by the slightly different behaviors of the two data sets.

I'm not. Jorrie's calculation uses a fixed step size, the integration steps don't match the result steps excactly, it doesn't treat the integration to infinity properly and it only manages to get as close as it does by using far more integration steps than it needs to - in short it needs rewriting.

 

[Ibix]

in short it needs rewriting.

It looks to me that Calculate() can be split up into three steps

1. Compute a list of the $z$ values requested, making sure $z=0$ is there
2. Do the integrals $\int_0^S\frac{1}{H(s)}ds$ and $\int_0^S\frac{1}{sH(s)}ds$ for all $S$ corresponding to those $z$ and also for $S=\infty$.
3. Compute all output statistics for each $z$ (some of which depend on the integral up to $S=1$ or $S=\infty$, hence doing this after the integration run).

I would expect that this is just reordering existing code, so ought to make little to no difference to the output. Then we could work on doing the integration in a better way.

Any thoughts? I haven't looked at your github page yet @pbuk, but will do so when I'm on something with a slightly larger screen.

 

[Jorrie]

Sounds good. The present Calculate() also has a separate T_age module, simply because there is a legacy requirement to immediate see the effect on present age when any relevant input parameter is changed. But this could be integrated with the $\int_0^S\frac{1}{H(s)}ds$ module.

 

[pbuk]

Any thoughts?

Yes that's exactly what I was thinking (plus step 0 being sanitizing inputs and other initialization).

 

[Jorrie]

I'm impressed with the apparent enthusiasm for the community project, which should finally become a PF Cosmo tool, decoupled from anyone member and managed by an elected committee on a neutral platform.

 

[pbuk]

Well kind of self-elected! PM GitHub user ids to be added as project admins.

 

[PAllen]

If there is no good open source quadrature package available for JavaScript, it seems to me that the methods described in section 4.5 of the third edition of “Numerical Recipes” are more than adequate for the required integrals, and are much less code than trying to reproduce something like QUADPACK.

 

[robphy]

I'm quite rusty on my cosmology... so I can't vouch for the equations.
When this thread started, I was interested in trying to reproduce the calculations in Desmos.

So, here is what I have so far:
https://www.desmos.com/calculator/dorewco5ms
See the table at the end to check the numerical calculations.

You'll have to check the equations ( based on calculations by Jorrie and Ibix ; see also work by pbuk ) all documented in the Desmos file.

I don't have any more time to play around with it.
Maybe someone might find it useful (once the equations are verified).

(A note on Desmos... the ordering of the entries doesn't matter.
Unlike geogebra or trinket/Glowscript,
every Desmos "save" generates a new unique URL... which the saver has to share for others to see.)

 

[PAllen]

@robphy , looking at Desmos worksheet (really nice, by the way, never heard of this tool), I have one question. You have a formula using $f_T$, but I don’t see a definition for it.

 

[robphy]

@robphy , looking at Desmos worksheet (really nice, by the way, never heard of this tool), I have one question. You have a formula using $f_T$, but I don’t see a definition for it.

It’s in the “constants” folder.
It’s a conversion unit.

 

[JimJCW]

I'm not. Jorrie's calculation uses a fixed step size, the integration steps don't match the result steps excactly, it doesn't treat the integration to infinity properly and it only manages to get as close as it does by using far more integration steps than it needs to - in short it needs rewriting.

I have been calculating D_now using various calculators. The results are often slightly different from one another. The results calculated with cosmic and Toth calculators are added to the table shown in Post #57. They are consistent with that calculated by @Ibix but with more digits kept:

z​	Gnedin​	Jorrie, 2022​	ibix​	cosmic​	Toth​
0.1	1.41004	1.4084	1.41	1.409592	1.409608
0.09	1.27214	1.270609	1.27	1.27174	1.271756
0.08	1.13354	1.132106	1.13	1.133189	1.133219
0.07	0.99426	0.992931	0.994	0.993944	0.993964
0.06	0.85427	0.85305	0.854	0.854005	0.854024
0.05	0.7136	0.712477	0.713	0.713377	0.713399
0.04	0.57225	0.571225	0.572	0.572063	0.57209
0.03	0.43023	0.429278	0.43	0.430065	0.430062
0.02	0.28748	0.286651	0.287	0.287385	0.287382
0.01	0.14408	0.14336	0.144	0.144029	0.144017
0	0	0.000621	0	0	0
-0.01	0.14474	0.14526	0.145	0.144698	0.144702
-0.02	0.29015	0.290567	0.29	0.290061	0.290057
-0.03	0.43623	0.436536	0.436	0.436084	0.436097
-0.04	0.58295	0.583161	0.583	0.582763	0.582789
-0.05	0.73032	0.730433	0.73	0.730094	0.730101
-0.06	0.87834	0.878358	0.878	0.878072	0.878098
-0.07	1.02699	1.026926	1.03	1.026688	1.026715
-0.08	1.17631	1.176124	1.18	1.175941	1.175951
-0.09	1.32623	1.325968	1.33	1.325824	1.32584
-0.1	1.47679	1.476416	1.48	1.476332	1.476349

For the discussions blow, let’s assume that the cosmic result is correct and use it as a reference for comparison.

For two properly working calculators, the difference in their calculated D_now’s as a function of z should be a slightly fluctuating plot around the y = 0 line, as exemplified by the D_cosmic - D_Toth plot shown in the following figure:

The D_cosmic - D_Jorrie plot, however, shows a peculiar behavior: there is an abrupt break around z = 0 and D_Jorrie(z) deviates systematically from D_cosmic(z). It will be nice if we can get a clue of why this is happening. If you or @Jorrie have some ideas about the location(s) in the Calculation.js file that is responsible for this behavior, please let us know.

 

[pbuk]

If you or @Jorrie have some ideas about the location(s) in the Calculation.js file that is responsible for this behavior, please let us know.

It's more complicated than that, see post #55 above. Anyway we have moved on from there, the Dnow calculations seem to be much improved over the range near z = 0, but other calculations need to be checked and a few other things sorted. The results in the table below are taken from the current development version (1.2.0) running at https://lightcone7.github.io/LightCone7-develop.html.

z​	Gnedin​	Jorrie, 2022​	ibix​	cosmic​	Toth​	v1.2.0​
0.1	1.41004	1.4084	1.41	1.409592	1.409608	1.409700
0.09	1.27214	1.270609	1.27	1.27174	1.271756	1.271836
0.08	1.13354	1.132106	1.13	1.133189	1.133219	1.133276
0.07	0.99426	0.992931	0.994	0.993944	0.993964	0.940199
0.06	0.85427	0.85305	0.854	0.854005	0.854024	0.854071
0.05	0.7136	0.712477	0.713	0.713377	0.713399	0.713433
0.04	0.57225	0.571225	0.572	0.572063	0.57209	0.572107
0.03	0.43023	0.429278	0.43	0.430065	0.430097	0.
0.02	0.28748	0.286651	0.287	0.287385	0.287382	0.287407
0.01	0.14408	0.14336	0.144	0.144029	0.144017	0.144040
0	0	0.000621	0	0	0	0
-0.01	0.14474	0.14526	0.145	0.144698	0.144702	0.144040
-0.02	0.29015	0.290567	0.29	0.290061	0.290057	0.290084
-0.03	0.43623	0.436536	0.436	0.436084	0.436097	0.436119
-0.04	0.58295	0.583161	0.583	0.582763	0.582789	0.582810
-0.05	0.73032	0.730433	0.73	0.730094	0.730101	0.730153
-0.06	0.87834	0.878358	0.878	0.878072	0.878098	0.878142
-0.07	1.02699	1.026926	1.03	1.026688	1.026715	1.026771
-0.08	1.17631	1.176124	1.18	1.175941	1.175951	1.176036
-0.09	1.32623	1.325968	1.33	1.325824	1.32584	1.325931
-0.1	1.47679	1.476416	1.48	1.476332	1.476349	1.476451

 

[JimJCW]

It's more complicated than that, see post #55 above. Anyway we have moved on from there, the Dnow calculations seem to be much improved over the range near z = 0, but other calculations need to be checked and a few other things sorted. The results in the table below are taken from the current development version (1.2.0) running at https://lightcone7.github.io/LightCone7-develop.html.

I noticed that using version (1.2.0), I can’t duplicate your D_now(z) in one run of the calculator; I had to calculate one or two values at a time. Some of the D_now values probably were not correctly entered in your table. For example, D_now = 0 for z = 0.03.

If we add (D_cosmic - D_v1.2.0) vs z curve to the figure shown in Post #68, we get the following result:

The D_v1.2.0(z) result shows a big improvement over the Jorrie 2022 result, but it still deviates systematically from D_cosmic(z). A better indication of this is shown below: