Line absorption rate estimator #595

Tobychev · 2016-06-21T12:25:33Z

Summary

A start on the adding of calculating spectra using exact source function integration: an line absorption rate estimator.

Background

In Lucy 1999, A&A 345, a method of calculating emerging spectra based on the formal integral of the spectra is suggested, and part of that method requires adding an estimator of the line absorption rate that can be simplified to:

where epsilon is the packet energy, tau is the optical depth of each transition, and the sum is over all packets that pass trough the l -> u transition in each cell. Only the sum explicitly needs to be inside cmontecarlo.c, the other factors can be done in python

To do

Add the sum estimator to cmonecarlo
Fix up the plumbing so the sum estimator can be reached from python
Make unit tests for the sum estimator
Normalize the sums
Multiply with the extinction factor that depends on the tau for each transition

Added a counting estimator to cmontecarlo closely following the example of increment_j_blue_estimator. Updated the various structures and headerfiles to accomodate this and exposing the result the same way as j_blue is exposed.

…folding

unoebauer · 2016-06-21T12:30:22Z

tardis/montecarlo/src/cmontecarlo.h

@@ -13,6 +13,10 @@
 #define LOG_VPACKETS 0
 #endif

+#define RESONANCE 0
+#define DOWNBRANCH 1
+#define MACROATOM 2


I'd rather work with enums and not use preprocessor statements.

unoebauer · 2016-06-21T12:32:13Z

@Tobychev - looks already very good. Maybe let's change the title of the PR because I think this one should be only about the line absorption rate (or at most about the level absorption rate). Also, please make a list of milestones which you want to achieve with this PR

unoebauer · 2016-06-21T12:37:41Z

tardis/montecarlo/src/cmontecarlo.c

+#ifdef WITHOPENMP
+#pragma omp atomic
+#endif
+    storage->line_lists_Edotlu[line_idx] += rpacket_get_energy (packet);


This is not yet the final estimator, right?

It's all that we said is needed in the montecarlo code. There is a factor $1 - e^{\tau_{lu}}$ and some normalization needed after the run to get the complete estimator.

no - the factor 1 - exp(-tau_s) is not global. Instead it is different for each line interaction. As the subscribed lu suggests, it is a quantity which depends on the line transition. In fact, this is the Sobolev line optical depth. It will be different in each shell and for each line.

To clarify: the factor (1 - exp(-tau_s)) describes the interaction probability, i.e. the probability of a packet which comes into resonance with the line to interact with it.

Yes it depends on the transition, I just meant to say that the tau for each line is added to the estimator(s, the code above is also for each transition and each shell) once all monte carlo steps have finished.

@Tobychev - hmm, I guess you could do that. But I think it should be easier to add the factor already during the increment. This way you avoid looping over all entries at the end again.

unoebauer · 2016-06-21T15:00:56Z

@Tobychev - I think one milestone should be the addition of a unittest for your estimator increment routine

unoebauer · 2016-06-23T08:18:59Z

tardis/model.py

+        exptau = 1 - np.exp(- 
+                            self.runner.line_lists_tau_sobolevs.reshape(-1,
+                                self.tardis_config.structure["velocity"]["num"]) ) 
+        self.Edotlu = self.Edotlu_norm_factor*exptau*self.Edotlu_estimators


I wouldn't do this post-processing here. As I've said, the 1- exp(-tau) factor should already be included in the estimator increments in the C part.

unoebauer · 2016-06-23T08:21:17Z

@Tobychev - I wouldn't do any "estimator post-processing" in the python part. We want to have the line absorption rate estimators ready and available on the C or at least on the cython level. This way, we can do the intensive number crunching (calculating the level absorption rates and ultimately performing the integral) on the C or cython level

…dicate a last run of the monte carlo iteration, updating things as required

unoebauer · 2016-06-27T16:56:59Z

@wkerzendorf, @yeganer - it is unclear to us why the unit tests for the estimator doesn't work. Apparently line_lists_Edotlu is not properly allocated. Any ideas?

I'd like to merge this PR as soon as possible, even if the unit tests are broken (we can mark them skip, right?)

kdexd · 2016-06-28T06:23:17Z

tardis/montecarlo/tests/test_cmontecarlo.py

@@ -119,6 +119,9 @@ def model():
        line_lists_j_blues=(c_double * 2)(*([1.e-10] * 2)),
        line_lists_j_blues_nd=0,

+        line_lists_Edotlu=(c_double * 3)(1e-10,1e-10,1.0), # Init to an explicit array 


@Tobychev Please replace this initialization with

line_lists_Edotlu=(c_couble * 3)(*[1.e-10, 1.e-10, 1.0])

Not sure if you can use tuples too.
I think this might be causing the problem. @yeganer what do you think about it ?

yeganer · 2016-06-28T06:37:46Z

@unoebauer I added some comments including where I mention a possible cause for the bug.
Concerning the merge, I have to disagree with you. I'm not familiar with the planned implementation but I'm there are some thing with the current implementation I don't like, which could be improved. As we have the policy don't touch a running system unless it's needed I'd fix them right away.

yeganer · 2016-06-28T06:41:23Z

tardis/montecarlo/base.py

@@ -106,7 +107,7 @@ def legacy_update_spectrum(self, no_of_virtual_packets):
            self.spectrum_virtual.update_luminosity(
                self.montecarlo_virtual_luminosity)

-    def run(self, model, no_of_packets, no_of_virtual_packets=0, nthreads=1):
+    def run(self, model, no_of_packets, no_of_virtual_packets=0, nthreads=1,last_run=False):


I think we could get away without introducing this new argument.
What we should do (if possible), is to move the postprocess to MontecarloRunner.calculate_Edotlu(self, model).
This method could then be called anytime. It would make sense for Simulation to call this function after It's done with it's last iteration.

Hmm - @Tobychev had this idea initially. I was against this approach, since we need the post-processed estimators on the C or cython level. That's why I favoured the current approach.

Postprocessing the results can still be done on the C/Cython level even if we did this using calculate_Edotlu later. I think having a postprocessing class is a sensible step (this also goes for spectral creation out of the packets).

But we can definitely do this in another PR.

unoebauer · 2016-06-28T11:20:13Z

Looks good to me now - I'd like to merge. @tardis-sn/tardis-core, objections?

yeganer · 2016-06-28T11:28:26Z

tardis/montecarlo/tests/test_cmontecarlo.py

@@ -119,8 +119,7 @@ def model():
        line_lists_j_blues=(c_double * 2)(*([1.e-10] * 2)),
        line_lists_j_blues_nd=0,

-        line_lists_Edotlu=(c_double * 3)(1e-10,1e-10,1.0), # Init to an explicit array 
-        line_lists_Edotlu_nd=0,
+        line_lists_Edotlu=(c_double * 3)(*[0.0,0.0,1.0]), # Init to an explicit array 


Why are we initializing this estimator to something nonzero?

not quite sure what you meant.

To see that in properly increments in that case too, might as well test that right?

Tomas Bylund added 2 commits June 15, 2016 14:22

Added simple counting estimator

6f95494

Added a counting estimator to cmontecarlo closely following the example of increment_j_blue_estimator. Updated the various structures and headerfiles to accomodate this and exposing the result the same way as j_blue is exposed.

Added core of Edotlu estimator to cmontecarlo and the supporting scaf…

6ff58aa

…folding

unoebauer reviewed Jun 21, 2016
View reviewed changes

Tobychev changed the title ~~Exact source function integration~~ Line absorption rate estimator Jun 21, 2016

Tomas Bylund and others added 2 commits June 21, 2016 18:03

Defined normalization

467c23a

Added postprocessing of Edotlu estimator

733b694

unoebauer reviewed Jun 23, 2016
View reviewed changes

Tobychev added 4 commits June 23, 2016 14:53

Moved the postprocessing of the estimator to and added argument to in…

e0c5a1d

…dicate a last run of the monte carlo iteration, updating things as required

Merge branch 'tardis-upstream' into Lucy345

b39cae9

Started adding unit tests

598abc0

Removed suggested cleanup

c863f78

kdexd reviewed Jun 28, 2016
View reviewed changes

yeganer reviewed Jun 28, 2016
View reviewed changes

Fixed problem with unit tests

3457c7b

yeganer reviewed Jun 28, 2016
View reviewed changes

wkerzendorf merged commit 82c7088 into tardis-sn:master Jun 28, 2016

ftsamis mentioned this pull request Dec 16, 2016

[TEP006] Aggregate TARDIS restructure and configuration system cleanup #652

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Line absorption rate estimator #595

Line absorption rate estimator #595

Tobychev commented Jun 21, 2016 •

edited

Loading

unoebauer Jun 21, 2016

unoebauer commented Jun 21, 2016

unoebauer Jun 21, 2016 •

edited

Loading

Tobychev Jun 21, 2016 •

edited

Loading

unoebauer Jun 21, 2016

unoebauer Jun 21, 2016

Tobychev Jun 21, 2016

unoebauer Jun 21, 2016

unoebauer commented Jun 21, 2016

unoebauer Jun 23, 2016

unoebauer commented Jun 23, 2016 •

edited

Loading

unoebauer commented Jun 27, 2016

kdexd Jun 28, 2016 •

edited

Loading

yeganer commented Jun 28, 2016

yeganer Jun 28, 2016

unoebauer Jun 28, 2016 •

edited

Loading

wkerzendorf Jun 28, 2016

wkerzendorf Jun 28, 2016

unoebauer commented Jun 28, 2016

yeganer Jun 28, 2016

wkerzendorf Jun 28, 2016

Tobychev Jun 28, 2016

Line absorption rate estimator #595

Line absorption rate estimator #595

Conversation

Tobychev commented Jun 21, 2016 • edited Loading

Summary

Background

To do

Choose a reason for hiding this comment

unoebauer commented Jun 21, 2016

unoebauer Jun 21, 2016 • edited Loading

Choose a reason for hiding this comment

Tobychev Jun 21, 2016 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

unoebauer commented Jun 21, 2016

Choose a reason for hiding this comment

unoebauer commented Jun 23, 2016 • edited Loading

unoebauer commented Jun 27, 2016

kdexd Jun 28, 2016 • edited Loading

Choose a reason for hiding this comment

yeganer commented Jun 28, 2016

Choose a reason for hiding this comment

unoebauer Jun 28, 2016 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

unoebauer commented Jun 28, 2016

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Tobychev commented Jun 21, 2016 •

edited

Loading

unoebauer Jun 21, 2016 •

edited

Loading

Tobychev Jun 21, 2016 •

edited

Loading

unoebauer commented Jun 23, 2016 •

edited

Loading

kdexd Jun 28, 2016 •

edited

Loading

unoebauer Jun 28, 2016 •

edited

Loading