Require a new set of python wrappers that enables user to:

• Initialise a geometry and field map
• Build a set of transfer maps
• Fire beam ellipses through the transfer maps
• Fire particles through the transfer maps
• Return particles, <some> transfer maps, beam ellipse as a thin wrapper to a numpy.matrix

• Modify the geometry
• Modify the field map

This is designed to make e.g. #1292 easier to do

Some general comments:

• Wrapper code will be placed in py_cpp/optics directory
  • May need to modify build script to make sure it builds in child directory
• We use numpy for all matrices

Going through the requirements:

Initialise a geometry and field map¶

• We can already initialise a geometry and field map from python using globals module

Build a set of transfer maps - 2 days¶

• As the public interface of all *OpticsModel classes is the same, propose implementing a single OpticsModel wrapper that takes a string argument indicating which base class to use
• This simply does init - and calls relevant underlying OpticsModel class. No need for anything else.
• Extension: fire particles using native G4 integration; propose adding a "transfer map" that enables this functionality

Fire beam ellipses through the transfer maps - 2 days¶

• Implement a CovarianceMatrix wrapper that wraps following functions
  • explicit CovarianceMatrix(const Matrix<double>& matrix);
  • static const CovarianceMatrix CreateFromPennParameters(...)
  • static const CovarianceMatrix CreateFromTwissParameters(...)
  • Note as an alternative to wrapping the last two, xboa has functions that will return a numpy matrix based on Penn/Twiss; may be easier to use this

Fire particles through the transfer maps - 2 days¶

• Implement a PhaseSpaceVector wrapper that takes x, px, y, py, z, pz, t, energy
• Extension: Implement a PhaseSpaceVector wrapper that takes an xboa.Hit as control. Means extracting phase space information from the hit

Modify the geometry, Modify the field map - 1 week¶

• Implement a wrapper for MiceModules that enables one to interface with the MiceModules objects; Json may be my friend here
  • Either wrap the SetPropertyDoubleThis/SetPropertyIntThis/etc functions (preferred)
  • Or make a json interface
• Modify G4 geometry needs investigating; I recall writing down somewhere how to do this; google is my friend
• Update() function that rebuilds BTFieldConstructor (and geometry)

Total time estimate: 2 weeks of reasonable work

1. Build a set of transfer maps

First pass done in bzr+ssh://bazaar.launchpad.net/~chris-rogers/maus/1293/ src/py_cpp/optics/OpticsModel.hh; can initialise an OpticsModel class. Things I havent done:

• Added the various Transport(...) classes
• Hardcoded to LinearApproximationOpticsModel
• Placed in optics_model python module; probably when I am done I will move everything into an optics python module

So what I have done is figured out/reminded myself how to wrap a C++ class in python and turn it into a python class. Total time to do this was about 5 hours

Added wrapper for PyCovarianceMatrix. This uses numpy matrices for the matrix work.

Things I have done:

• Module and class initialisation
• Wrappers for various initialisation functions

Things I havent done:

• Added a Transport(...) method Still needs z_start and z_end or equivalent
• Added accessors to matrix elements; could go the whole hog and make it an iterable
• Added print oeprator
• Added equality operators
• Class docstring needs init doc added
• Style test will fail
• Unit test needs extension when I have accessors

Total time to do this was about 10 hours.

I spent another 5 or 6 hours on this, just adding the Transport(CovarianceMatrix) method. The reason this took a bit longer was

• I had to remember how to make C API functions callable from another python module
• Transport(CovarianceMatrix, start, end) was throwing a segmentation fault when no planes were loaded (=> no transfer maps) and it took me a bit of time to isolate the problem. There is a NULL dereference in this situation. At the moment I am wrapping Transport(CovarianceMatrix, end) instead. I will have to find a better way but it will do for now.

Note I modified OpticsModel to hard coded PolynomialOpticsModel

I did some work on the Covariance Matrix class and also added PhaseSpaceVector class and Transport wrapper in Optics. I updated #note-4 with some strike throughs for things I have done. To do in PhaseSpaceVector:

Add accessors to matrix elements
Add print operator
Unit test extensions
Docstring checks esp init
Style tests will no doubt fail
Modify Transport so that it takes z start and z end

Next up - geometry interface

Another ~2.5 hours, doing some tidying around PyPhaseSpaceVector (updated note #6 with changes) pushed as r957 in bzr+ssh://bazaar.launchpad.net/~chris-rogers/maus/1293/

more jobs:

• optics_model unit tests should just use a drift.
• add an integration test that duplicates the optics calculations done using current optics test C++ code.
• Because docstrings can't be put on init functions, init should be dumb code that does nothing and we should instead define static C functions like "new_..." or "create_..." which do anything more complicated.

Another ~ 3 hours

Cleaned up init functions so they don't take any parameters. Moved stuff to create_blah functions (where I can put a docstring).

Added more tests on covariance matrix and phase space vector. PhaseSpaceVector at least isn't acting as expected in a drift. I need to dig a little deeper.

Now looking at the geometry side. Making a wrapper for MiceModule - done the header and init function. Not quite sure how Optics handles it's geometry, it probably isn't taking stuff from Globals (and that is a bug).

Plan for PyMiceModule

• init fuction
• str function
• get_name
• set_property_value
• get_property_value
• get_children (child MiceModules)
• set_children
• get_parent
• set_parent
• globals interface
  • get_mice_modules
  • set_mice_modules - which reinitialises field maps and G4 geometry as appropriate

Another 5 hours - partial implementation of MiceModule wrapper

Another 3 hours or so - some things in #8 commented. Another job:

• Check memory usage carefully.

Another two hours - I started implementing set_parent but after thinking about it for a while I realised that this will cause problems. Because there is no reverse memory management in the tree (child does not own memory of parent) there is no way to keep memory under control. Say we have a tree

A
|--\
|  |
B  C

What happens if we take a pointer to B but then delete A? Only way is to forbid B from accessing A (at least on python side). If A is deleted then B is safe because it knows nothing about A on python side.

Probably means that similar problem exists on C++ side and is a lurking evil... which is not exposed right now...

This is probably just another case of me sticking my nose into things I don't understand, but I don't see why this is a problem. If A is deleted then so shouldn't B and C? Generally you don't want orphaned tree nodes.

Say we make a python reference to B, call it py_B? Remember B has no right to ownership of A.

• If py_B is a shallow copy then when we delete B the python reference disappears -> memory fault
• If py_B is a deep copy then when we delete B, py_B is okay; but when we delete A, py_B's pointer up to A is no longer valid -> memory fault

Rogers, Chris wrote:

Say we make a python reference to B, call it py_B? Remember B has no right to ownership of A.

Ah, you're talking about external references/copies to tree nodes...

• If py_B is a shallow copy then when we delete B the python reference disappears -> memory fault

Right, but this is an issue regardless. If you have a shallow copy of any of the nodes and delete it or a connected one, you potentially have problems.

• If py_B is a deep copy then when we delete B, py_B is okay; but when we delete A, py_B's pointer up to A is no longer valid -> memory fault

Yep. The issue isn't really the data structure having too many linkages. The issue is your mistrust in whoever is using the data structure. The only way around this if you don't trust the users is to make sure that "illegal" references to the nodes can't be made by hiding the node objects behind an interface of some sort. If there is only one actor using the tree at a time (i.e. the interface code or some core code that does pre-prescribed actions), then it's simply the singleton actor's job not to do something stupid.

Right - but the point of Python is it's a scripting tool/open to dumb users so should never expose a possible segmentation fault or anything. For example, we could have py_B_1 and py_B_2 both of which have pointers to A (which could be a python object). What happens if we delete py_B_1? py_B_2? What happens if A was a python object already (py_A) and user deletes A anyway? When is A deleted? There may be a way through, but it is a little tricky. I would rather not modify MiceModule code also.

Absolutely! When I said "hiding the node objects behind an interface", I meant using data encapsulation within Python. In other words, create a class that allows traversal of the tree and access to node data but never hands the user a node reference. I realize python doesn't have private data members so this isn't a complete guarantee of safety, but it's the best one can do in these regards. At least you can say "not supported" if someone circumvents the interface and tries to manipulate the internal data themselves.

Another 12 hours or so

• Change Construct(...) to make a trivial geometry (i.e. 1 mm^3 Galactic cube)
  • Move main construction into a Reset function to reduce testing/maintenance over head
• Make a ResetFields(...) function that has the job of resetting the fields - and incidentally sets the main fields
• Make a ResetGeometry(...) function that has the job of resetting the geometry - and incidentally sets the main geometry
• Break into sub functions and test each of the blocks in AddDaughter - that is
  • Set up top level G4 volume
  • Set up special G4Detector volumes
  • Set up "normal" volumes
  • Set up "None" volumes (recursively check for child volumes in None type and throw if found)
  • SetUserLimits
  • SetVisAttributes
  • BuildSensitiveDetector
  • Recurse through the MiceModule tree

Tested everything and checked that it works - mostly by e.g. firing particles and checking that they do what is expected. Test is a little slow to run but meh. ~ 8 hours.

Now need to spend a day running valgrind and sorting through memory leaks. This may address the online recon memory leak #1328, but should at least make it easier to see what is going wrong...

I think I have done all of the memory cleanup from Geant4 I can manage. Remaining leaks are:

• There is something in the depths of Geant4 physics list which appears to be unreachable memory, initialised in MAUS::MAUSPhysicsList::ConstructProcess()
• Geant4 allocates memory in the optics routines of the Materials list (gamma energy vs refractive index lookup) that is unreachable

• SciFiPlane, CkovMirror, KLFiber, KLGlue all take ownership of memory for the Geant4 volumes (physical, logical) but need method to hand this memory to Geant4 which should properly own memory for G4 volumes. So currently these objects cannot be deleted.

• Worryingly when MAUS::MAUSGeant4Manager::Tracking is called some darkness happens, again in the physics table, that means Geant4 generates unreachable memory.
• Some G4NavigationHistory thing, whatever that means

Another 12 hours or so of tidying up... styles tests pass, unit tests all pass. Checking integration tests now...

Peter, would you like to do a code review? Up to you - it's in branch

bzr+ssh://bazaar.launchpad.net/~chris-rogers/maus/1293/

Another 8 hours or so of tidying up. Now passes on test server. I haven't done:

• an example usage
• integration test (e.g. can I reproduce results from e.g. existing optics mapping test)
• any "top level" documentation e.g. user documentation on maus_cpp python library - but there are docstrings on everything
• should all the optics stuff go into an optics submodule - e.g. maus_cpp.optics.covariance_matrix.CovarianceMatrix?

Total time spent ~ 66 hours; 66 / 8 => a little over 8 man days. About as expected, 2 weeks of decent work (but Rogers has lots of distractions so really it took more like 4 weeks real time). In fact, the break down wasn't quite right - the python API stuff took a bit less time than the Geant4 and MiceModules stuff.

Sure, I'll take a look and post any comments here.

Thanks. Just be aware, to get it to build you have to run

${MAUS_ROOT_DIR}/third_party/bash/39numpy.bash

I need numpy development libraries (.h files) in order to build the C API against numpy.

I can't get it to build. It continues to complain about not being able to find numpy even after I install it separately as you indicated. I tried building the third party software manually and just made a mess out of things. I don't feel like spending any more time on this. The only issue I found skimming through the code was the following:

py_cpp/optics:

• PyOpticsModel.hh does not seem to have a function declaration for transport_phase_space_vector().

Ah, sorry about that - it did build on the test server, not sure why it isn't building for you. Thanks for trying anyway.