Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Bokeh plotter #714

Merged
merged 16 commits into from
Nov 9, 2018
Merged

Bokeh plotter #714

merged 16 commits into from
Nov 9, 2018

Conversation

watsonjj
Copy link
Contributor

@watsonjj watsonjj commented Apr 20, 2018

Created an interactive event plotter that utilises the bokeh package (https://github.com/bokeh/bokeh)

This plotter shows both the camera and the waveform. It allows the selection of events, pixels, and time slices, and also allows the viewing and configuration of various calibration stages.

All EventSources are compatible with it, and all geometries are supported. However, all pixels are displayed as squares at the moment.

This tool requires a bokeh server to run:
bokeh serve ctapipe/tools/bokeh/file_viewer

To pass arguments to the tool:

bokeh serve ctapipe/tools/bokeh/file_viewer --args -f /Volumes/gct-jason/data/170330/onsky-mrk501/Run05477_r1.tio

Also includes
- ctapipe.utils.rgbtohex
- ctapipe.visualization.bokeh
- ctapipe.plotting.bokeh_event_viewer
* master:
  Copy changes to a Factory's trait default value to its subclasses (cta-observatory#713)
@watsonjj
Copy link
Contributor Author

eventdisplayer

An animation of the file_viewer (some updates have been made since this gif's creation)

@maxnoe
Copy link
Member

maxnoe commented Apr 20, 2018

Looks good. Will checkout the code later. 👍

exe.run()


if 'bk_script' in __name__:
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this is very unusual. What is this for?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

To run a bokeh server you have to execute bokeh serve on the script instead of python. When you do this, __name__ equals something like "bk_script_dc88af7ef51b4e0d85854896dd662ed0", where the third term seems to be random.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

However, there may be a way to avoid having to use bokeh serve if I can understand how to utilise this: http://bokeh.pydata.org/en/latest/docs/user_guide/server.html#embedding-bokeh-server-as-a-library

@codecov
Copy link

codecov bot commented Apr 20, 2018

Codecov Report

Merging #714 into master will increase coverage by 1.39%.
The diff coverage is 88.72%.

Impacted file tree graph

@@            Coverage Diff             @@
##           master     #714      +/-   ##
==========================================
+ Coverage   71.93%   73.33%   +1.39%     
==========================================
  Files         204      211       +7     
  Lines       11061    12054     +993     
==========================================
+ Hits         7957     8840     +883     
- Misses       3104     3214     +110
Impacted Files Coverage Δ
ctapipe/utils/tests/test_rgbtohex.py 100% <100%> (ø)
ctapipe/plotting/tests/test_bokeh_event_viewer.py 100% <100%> (ø)
ctapipe/utils/rgbtohex.py 100% <100%> (ø)
ctapipe/visualization/tests/test_bokeh.py 100% <100%> (ø)
ctapipe/tools/tests/test_tools.py 100% <100%> (ø) ⬆️
ctapipe/tools/bokeh/file_viewer.py 72.33% <72.33%> (ø)
ctapipe/plotting/bokeh_event_viewer.py 90.74% <90.74%> (ø)
ctapipe/visualization/bokeh.py 92.27% <92.27%> (ø)
... and 7 more

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 0558f39...b91a275. Read the comment docs.

Copy link
Member

@dneise dneise left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we should not use C-extensions if possible.

#include <stdint.h>
#include <stdio.h>

extern "C" void rgbtohex(const uint8_t* rgb, size_t n_pix, char* hex)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

TL;DR: I believe in this particular case the C-extension is not necessary.

I believe we should avoid C-extensions as long as possible and I assume you believe so too (as opposed to writing C-extension just for fun). So I assume you had no chance to avoid this, since with pure Python this was performing horrible. In such a case I believe C-extensions are fine.

Still, I personally consider writing a C-extension so exceptional that a little discussion of it might not be too annoying. I would have liked the author to present maybe some numbers in a comment or in the PR-discussion to backup the decision to go for a C-extension.

Since these numbers were missing, I tried to understand myself what was so slow here. For this I created a Python implementation of what I understood the C-implementation to do and profiled it by hand (using simply %timeit in my ipython session).

The result was a second implementation of intensity_to_hex which is simply named intensity_to_hex2 which is pure python and happens to be even faster than the version using the C-extension.


Here is the comparison:

In [1]: from rgbtohex import *

In [2]: array = np.random.rand(10000)

In [3]: assert (intensity_to_hex(array) == intensity_to_hex2(array)).all()

In [4]: %timeit intensity_to_hex(array)
2.99 ms ± 13.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

In [5]: %timeit intensity_to_hex2(array)
1.75 ms ± 10.6 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

This is the original intensity_to_hex (rewritten a bit and added the runtimes per line):

def intensity_to_hex(array, minval=None, maxval=None):
    hex_ = np.empty((array.size, 8), dtype='S1')   #   <1 us
    rgb = intensity_to_rgb(array, minval, maxval)  #  200 us
    rgbtohex(rgb, array.size, hex_)                # 1300 us
    return hex_.view('S8').astype('U8')[:, 0]      # 1300 us

This is the other variant:

def intensity_to_hex2(array, minval=None, maxval=None):
    import codecs
    hex_ = np.zeros((array.size, 8), dtype='B')   #    3 us
    rgb = intensity_to_rgb(array, minval, maxval) #  200 us

    hex_encoded = codecs.encode(rgb, 'hex')       #   39 us
    bytes_ = np.frombuffer(hex_encoded, 'B')      #   <1 us
    bytes_2d = bytes_.reshape(-1, 8)              #   <1 us
    hex_[:, 0] = ord('#')                         #    4 us
    hex_[:, 1:7] = bytes_2d[:, 0:6]               #   48 us

    return hex_.view('S8').astype('U8')[:, 0]     # 1300 us

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is great!

Yeah, I couldn't understand how to do this efficiently in Python (I had to use a for loop), so I resorted to the C extension. I love your solution. Thanks!

@dneise
Copy link
Member

dneise commented Apr 23, 2018

I think this looks very nice.

It seems there are quite some viewers popping up here and there these days, which I believe might be necessary, but I also believe viewers are often quite heavy in terms of LOCs and therefore they might be considered expensive in terms of maintenance. I believe ctapipe has 15k lines of code and with 1.7k lines, this PR adds a whopping 10% on top.

I do not think that adding a lot of LOCs alone is a bad thing by itself, but let's assume for a moment that in 6 month from now we realize that of the 3 different existing viewers this particular one is not used by too many people... Then I personally would consider removing it maybe... or not if it turns out that is simply works without any maintenance.

That being said, I have no real clue of Bokeh, but I always wanted to learn it :-D

elif v == 'peakpos':
peakpos = e.dl1.tel[t].peakpos
if peakpos is None:
self.image = np.zeros(self.image.shape)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For the other cases you say:

if samples is None:
    self.image = None

but for peakpos, there is an exceptional handling of this case. Could you explain the difference?

@dneise
Copy link
Member

dneise commented Apr 23, 2018

I dream of a viewer, which does not need to be modified, when someone wants to plot a field from a Container, which was just added. I think, if there is a language which might be able to do this, then it is Python.

I'm not sure how to do this, but as a first step I tried to refactor the case handling for the view_options out of the BokehEventViewerCamera._set_image method.

So _set_image could maybe look like this ..

    def _set_image(self):
        e = self.event
        v = self.view
        t = self.telid
        c = self.channel
        time = self.time
        if not e:
            self.event_viewer.log.warning("No event has been provided")
            return

        tels = list(e.r0.tels_with_data)
        if t is None:
            t = tels[0]
        if t not in tels:
            raise KeyError("Telescope {} has no data".format(t))

        try:
            self.image = self._view_options[v](e, t, c, time)
            self.fig.title.text = '{0} (T = {1})'.format(v, time)
        except:
            self.image = None

If one would consider to combine the name and the action for each of the available view-options in a table ... maybe like this:

class BokehEventViewerCamera(CameraDisplay):
    def __init__(self, event_viewer, fig=None):
        # [...]

        self._view_options = {
            'r0': lambda e, t, c, time: e.r0.tel[t].waveform[c, :, time],
            'r1': lambda e, t, c, time: e.r1.tel[t].waveform[c, :, time],
            'dl0': lambda e, t, c, time: e.dl0.tel[t].waveform[c, :, time],
            'dl1': lambda e, t, c, time: e.dl1.tel[t].image[c, :],
            'peakpos': lambda e, t, c, time: e.dl1.tel[t].peakpos[c, :],
            'cleaned': lambda e, t, c, time: e.dl1.tel[t].cleaned[c, :, time],
        }

If we would go for this way of writing it down, I believe what we see is, that the places in code, we need to modify to allow for plotting of an additional field from the container was reduced.

But we also lost at least one feature:

non-specific title text: self.fig.title.text

Now it says: self.fig.title.text = '{0} (T = {1})'.format(v, time) where before it was very nicely handcrafted for each view_option. If we would be fine with dropping the feature of dedicated title texts, we would gain some flexibility.

The next step would be to dynamically create this "view option lookup table":

       self._view_options = {
            'r0': lambda e, t, c, time: e.r0.tel[t].waveform[c, :, time],
            'r1': lambda e, t, c, time: e.r1.tel[t].waveform[c, :, time],
# ...

from the Container itself. I believe one could write a function, which looks at each content of a container, by traversing it recursively, looking for contents that look like "plotable arrays". As a parameter the function could get the container and the number_of_pixels and it returns the "path" to each array which looks "plotable". Plotable could mean: it has one axis which is equal to the number of pixels...

I am pretty sure that users playing around with new algorithms would love if they could plot "almost anything" as long as it has the right dimensions.

This all being said: I am sure we can merge this PR without this feature I am dreaming about and whoever has time can try to work on that later. But I wanted to tell you guys about this dream and maybe learn if we share this dream or not.

@kbruegge
Copy link
Member

This looks really nice. Good Job. I'm just wondering whether this might be better located in a standalone package outside of ctapipe?

@kosack
Copy link
Contributor

kosack commented Apr 23, 2018

re @dneise - Such a viewer could be even more generic, since Containers can be looped over and accessed by field name string. One could provide a dropdown with the base container name, telescope name, and field name which would automatically be "r0, r1, dl0, dl1...", "1,2,3,4...", "image", "waveform", or even simply search recursively for all fields in any subcontainer that have the same geometry as the image, or something like that. Something to think about for a future feature.

re where it should go: I'm also not sure about whether this should be in ctapipe, or wether there should have a separate "ctapipe-viewers" repo. It would take some time to set another repo up, but might make it more clean and easier to maintain. Its a similar question for the "flow" module, which is a bit out of date and could be easily separated. Perhaps for now, it's ok to have it in the main repo and we can separate things later?

Another comment: I would suggest at the very least adding a top-level Tool for this in the ctapipe/tools/ that can be used by running ctapipe-event-viewer --input <filename> or something like that, so that the code is trivially usable by end users. Otherwise, it will just be forgotten. That's a tool that would be quite nice and which we are missing so far.

@watsonjj
Copy link
Contributor Author

@dneise I like that "dream" too. I have incorporated the updates you have suggested, it should help make the utilisation of an automatically generated "view_options" easier a future update.

@dneise
Copy link
Member

dneise commented Apr 23, 2018

it should help make the utilisation of an automatically generated "view_options" easier a future update.

that's great ... I hope so too

Used bokeh.server.server.Server to run file_viewer through `python`
instead of `bokeh serve`
Created entry_point ctapipe-event-viewer
Rearranged directory structure of tools/bokeh
@watsonjj
Copy link
Contributor Author

I have included all the comments so far. The file viewer can now be ran with either:
python ctapipe/tools/bokeh/file_viewer.py
OR
ctapipe-event-viewer

Codacy is complaining about jinja2.Environment, but I tried the change it suggests and it breaks bokeh entirely.

If everyone is happy, I am ready for this PR to be merged.

@dneise
Copy link
Member

dneise commented Apr 23, 2018

At the moment it seems to break when I feed a digicampipe fits file into it:

For testing I use the example_10_evts.000.fits.fz one can fine here:
https://github.com/cta-sst-1m/protozfitsreader/tree/master/protozfits/tests/resources

(cta-dev) dneise@lair:~/cta/ctapipe$ ctapipe-event-viewer -f /path/to/example_10_evts.000.fits.fz
INFO [BokehFileViewer] (tool/initialize): ctapipe version 0.5.3.post117+git31c39a2
INFO [BokehFileViewer] (tool/run): Starting: BokehFileViewer
INFO: Obtaining SST1MEventSource from EventSourceFactory [BokehFileViewer.EventSourceFactory._product]
INFO: INPUT PATH = /home/dneise/sst_new/digicampipe/digicampipe/tests/resources/example_10_evts.000.fits.fz [BokehFileViewer.SST1MEventSource.__init__]
INFO: Obtaining NeighbourPeakIntegrator from ChargeExtractorFactory [BokehFileViewer.ChargeExtractorFactory._product]
INFO: Obtaining NullWaveformCleaner from WaveformCleanerFactory [BokehFileViewer.WaveformCleanerFactory._product]
INFO: Obtaining NullR1Calibrator from CameraR1CalibratorFactory [BokehFileViewer.CameraR1CalibratorFactory._product]
INFO: Using NullR1Calibrator, if event source is at the R0 level, then r1 samples will equal r0 samples [BokehFileViewer.NullR1Calibrator.__init__]
INFO: Applying no data volume reduction in the conversion from R1 to DL0 [BokehFileViewer.CameraDL0Reducer.__init__]
WARNING: Seeking to event by looping through events... (potentially long process) [BokehFileViewer.EventSeeker.__getitem__]
Traceback (most recent call last):
  File "/home/dneise/anaconda3-5.1.0/envs/cta-dev/bin/ctapipe-event-viewer", line 11, in <module>
    load_entry_point('ctapipe', 'console_scripts', 'ctapipe-event-viewer')()
  File "/home/dneise/cta/ctapipe/ctapipe/tools/bokeh/file_viewer.py", line 414, in main
    exe.run()
  File "/home/dneise/cta/ctapipe/ctapipe/core/tool.py", line 167, in run
    self.start()
  File "/home/dneise/cta/ctapipe/ctapipe/tools/bokeh/file_viewer.py", line 134, in start
    self.event_index = 0
  File "/home/dneise/cta/ctapipe/ctapipe/tools/bokeh/file_viewer.py", line 165, in event_index
    self.event = self.seeker[val]
  File "/home/dneise/cta/ctapipe/ctapipe/tools/bokeh/file_viewer.py", line 218, in event
    self.viewer.event = val
  File "/home/dneise/cta/ctapipe/ctapipe/plotting/bokeh_event_viewer.py", line 423, in event
    sub.change_event(val, self.telid)
  File "/home/dneise/cta/ctapipe/ctapipe/plotting/bokeh_event_viewer.py", line 96, in change_event
    self.event = event
  File "/home/dneise/cta/ctapipe/ctapipe/plotting/bokeh_event_viewer.py", line 89, in event
    self._update_geometry()
  File "/home/dneise/cta/ctapipe/ctapipe/plotting/bokeh_event_viewer.py", line 74, in _update_geometry
    self.geom = e.inst.subarray.tel[t].camera
KeyError: 1

Update:
It breaks similarly when I feed the NectarCam file into it, just with a different KeyError:

(cta-dev) dneise@lair:~/cta/ctapipe$ ctapipe-event-viewer -f /home/dneise/sst_new/protozfitsreader/protozfits/tests/resources/example_9evts_NectarCAM.fits.fz 
INFO [BokehFileViewer] (tool/initialize): ctapipe version 0.5.3.post117+git31c39a2
INFO [BokehFileViewer] (tool/run): Starting: BokehFileViewer
INFO: Obtaining NectarCAMEventSource from EventSourceFactory [BokehFileViewer.EventSourceFactory._product]
INFO: INPUT PATH = /home/dneise/sst_new/protozfitsreader/protozfits/tests/resources/example_9evts_NectarCAM.fits.fz [BokehFileViewer.NectarCAMEventSource.__init__]
INFO: Obtaining NeighbourPeakIntegrator from ChargeExtractorFactory [BokehFileViewer.ChargeExtractorFactory._product]
INFO: Obtaining NullWaveformCleaner from WaveformCleanerFactory [BokehFileViewer.WaveformCleanerFactory._product]
INFO: Obtaining NullR1Calibrator from CameraR1CalibratorFactory [BokehFileViewer.CameraR1CalibratorFactory._product]
INFO: Using NullR1Calibrator, if event source is at the R0 level, then r1 samples will equal r0 samples [BokehFileViewer.NullR1Calibrator.__init__]
INFO: Applying no data volume reduction in the conversion from R1 to DL0 [BokehFileViewer.CameraDL0Reducer.__init__]
WARNING: Seeking to event by looping through events... (potentially long process) [BokehFileViewer.EventSeeker.__getitem__]
Traceback (most recent call last):
  File "/home/dneise/anaconda3-5.1.0/envs/cta-dev/bin/ctapipe-event-viewer", line 11, in <module>
    load_entry_point('ctapipe', 'console_scripts', 'ctapipe-event-viewer')()
  File "/home/dneise/cta/ctapipe/ctapipe/tools/bokeh/file_viewer.py", line 414, in main
    exe.run()
  File "/home/dneise/cta/ctapipe/ctapipe/core/tool.py", line 167, in run
    self.start()
  File "/home/dneise/cta/ctapipe/ctapipe/tools/bokeh/file_viewer.py", line 134, in start
    self.event_index = 0
  File "/home/dneise/cta/ctapipe/ctapipe/tools/bokeh/file_viewer.py", line 165, in event_index
    self.event = self.seeker[val]
  File "/home/dneise/cta/ctapipe/ctapipe/tools/bokeh/file_viewer.py", line 218, in event
    self.viewer.event = val
  File "/home/dneise/cta/ctapipe/ctapipe/plotting/bokeh_event_viewer.py", line 423, in event
    sub.change_event(val, self.telid)
  File "/home/dneise/cta/ctapipe/ctapipe/plotting/bokeh_event_viewer.py", line 96, in change_event
    self.event = event
  File "/home/dneise/cta/ctapipe/ctapipe/plotting/bokeh_event_viewer.py", line 89, in event
    self._update_geometry()
  File "/home/dneise/cta/ctapipe/ctapipe/plotting/bokeh_event_viewer.py", line 74, in _update_geometry
    self.geom = e.inst.subarray.tel[t].camera
KeyError: 0

@watsonjj
Copy link
Contributor Author

At the moment it seems to break when I feed a digicampipe fits file into it

It appears you aren't filling the subarray container in sst1meventsource.py. So not a fault of the file viewer :)

@dneise
Copy link
Member

dneise commented Apr 24, 2018

So not a fault of the file viewer :)

This answer makes me feel ridiculed and not taken seriously. I believe it was meant to be funny, but let me just tell you, it was not funny for me.

Apart from that I would like to mention these points

1. Thank you

You found a problem in SST1MEventSource. I will try to find out how correctly fill this subarray by reviewing other EventSources. I was not aware the inst Field of DataContainer was supposed to be used.

2. Personification of code

It appears you aren't filling the subarray container in sst1meventsource.py.

I was calling: ctapipe-event-viewer -f <some_path> and not filling anything in a container.

I assume by "you aren't filling" you refer to the SST1MEventSource I have been the author of? Here I would like to say, that in my understanding we all become equal "owners" of code as soon as it is merged into the master. As a consequence I believe phrases like "I do this and that in order to ..." or "you forgot to do this and that..." are simply wrong as soon as the code was merged into the master. From that point on it should in my little opinion only be "I think, we forgot to fix this and that" ... self accusations like "I think, I messed up here and there" are not good either, since there was a code review and others are just as responsible as the person trying to take all the blame.

2. Depending on deprecated features

I believe new code should try to not depend on deprecated features of the current code base.

count = Field(0, "number of events processed")
inst = Field(InstrumentContainer(), "instrumental information (deprecated")
pointing = Field(Map(TelescopePointingContainer),
'Telescope pointing positions')

If this viewer you are proposing to merge into the master does make use of deprecated features, then I believe it was good that we spotted it.

3. We should care about our users

So ... do I understand correctly that you find it acceptable to merge this viewer into the master, knowing that any user who wants to try it on a digicam or nectarcam file will see the same exception as I reported? How do you suppose our response should be, if users write issues asking why they cannot use the nice bokeh viewer on their files? Answering "not a fault of the viewer" is in my understand not very helpful.

@watsonjj
Copy link
Contributor Author

I apologise, I wrote that comment in a rush before I caught a bus so that you weren't left hanging without a response (after downloading the file and the protozfits software to double check the problem).

I was not aware this container is deprecated. To my understanding it is the only method available to obtain the pixel positions from the event container.

Yes, I should refer to us all as owners of the code. However, I (and the majority of contributors) aren't familiar with the organisation of SST1M, and were unqualified to make the sst1meventsource, and are subsequently unqualified to figure out the best method to supply the pixel positions that accompany the data. In the TargetIOEventSource, we utilise the TargetCalib library which has knowledge of the pixel positions of our camera.

Again, sorry for the careless response.

@dneise
Copy link
Member

dneise commented Apr 24, 2018

Thank you for your kind understanding...

* master: (60 commits)
  Add test that shows slicing breaks cam geom and fix it (cta-observatory#782)
  fix ctapipe build failure (cta-observatory#811)
  fix package name for yaml (should be pyyaml) (cta-observatory#810)
  Implement number of islands (cta-observatory#801)
  fixed ranges of cam-display so they correspond to fixed toymodel sims (cta-observatory#808)
  Fix unknown section example warning (cta-observatory#800)
  Fix timing parameters for case when there are negative values in image (cta-observatory#804)
  Update Timing Parameters (cta-observatory#799)
  speed up unit tests that use test_event fixture (cta-observatory#798)
  Add unit to h_max in HillasReconstructor (cta-observatory#797)
  Codacy code style improvements (cta-observatory#796)
  Minor changes: mostly deprecationwarning fixes (cta-observatory#787)
  Array plotting (cta-observatory#784)
  added a config file for github change-drafter plugin (cta-observatory#795)
  Simple HESS adaptations (cta-observatory#794)
  add test for sliced geometries for hillas calculation (cta-observatory#781)
  Impact intersection (cta-observatory#778)
  updated main documentation page (cta-observatory#792)
  Implement concentration image features (cta-observatory#791)
  Fix bad builds by changing channel name (missing pyqt package) (cta-observatory#793)
  ...

# Conflicts:
#	ctapipe/calib/camera/dl1.py
@watsonjj
Copy link
Contributor Author

watsonjj commented Nov 9, 2018

Can this PR be accepted yet?

The only failures are those present in the master.

@kosack
Copy link
Contributor

kosack commented Nov 9, 2018

yes, i think this is fine to merge once I get the tests working again. I'll do it shortly.

@kosack kosack merged commit 1efbd77 into cta-observatory:master Nov 9, 2018
@watsonjj watsonjj deleted the bokeh_plotter branch November 18, 2018 20:16
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

5 participants