MongoDB database handles in parallel constructs

[pavlis]

I am not 100% sure of this, but I think I've figure out why my parallel version of pwstack (part of pwmig I've been testing using mspass) is failing with an obscure error about not being unable to pickle a PyCapsule. My current working hypothesis is that the error is being created by a somewhat hidden (from a rookie perspective anyway) attempt to serialize a MongoDB Database class.

Let me break this algorithm down in stages so you can see how this is happening. The algorithm first creates a list of strings that are queries for MongoDB find:

allqueries=list()
    for sid in source_id_list:
        for rids in staids:
            # build_wfquery returns a dict with lat, lon,
            # i, j, and a query string  That is a simple
            # construct so don't think it will be a bottleneck.
            # May have been better done with a dataframe
            #q=dask.delayed(build_wfquery)(sid,rids)
            q=build_wfquery(sid,rids)
            # debug
            #print(q['ix1'],q['ix2'],q['fold'])
            allqueries.append(q)

Pretty stock python code.

Next I take that list of strings and turn it into a dask bag:

mybag=dask.bag.from_sequence(allqueries)

which is also pretty standard stuff.

I think the problem is right here with the next line:

mybag = mybag.map(lambda q : read_ensembles(db,q,control))

db is a mspasspy.db.Database handle and, I believe, is not serializable. I am pretty sure it needs to be serializable since it appears in the arg list to read_ensemble. Can you, @wangyinz , confirm that hypothesis?

The reason I'm not positive about this issue is that the error chain tracking this really "leads down a rathole" to again us a cliche. The error message python posts is totally obscure:

TypeError: cannot pickle 'PyCapsule' object

Furthermore, when I followed the call sequence with spyder I got into a very confusing chain of results I (preserved here](pavlis/parallel_pwmig#2) (the rathole in the cliche above). A simple summary of what I found was the data structure dask passes around is pretty opaque and finding the problem from the internals is a "needle in a haystack problem", which I guess is yet another cliche.

If, as I suspect strongly my hypothesis is true, there are some followup issues we need to deal with:

1. The parallel processing section of our user's manual needs an extensive discussion emphasizing strongly to warn the user not to implicitly serialize a Database class as I did.
2. The User's Manual also needs guidance on how one should interact with the database in in any parallel constructs.
3. I put this in discussion instead of issues because there are some design holes this process has revealed. I'll expand on them once you, @wangyinz , confirm my memory about this limitation of the Database class. One thing we will need for sure, however, is a general way to replace my (apparently incorrect) application of read_ensemble above.

[wangyinz]

hmmm.... what is this read_ensembles function? Probably this is something you defined in parallel_pwmig. Anyway, you are correct that the database handle is not serializable. However, we've run into this before and found that we need to be able to serialize it for the read_distributed_data function, so I have implemented the needed method to make the handler serializable here:

mspass/python/mspasspy/db/database.py

Lines 193 to 205 in 3fe0c52

	def __getstate__(self):
	ret = self.__dict__.copy()
	ret["_Database__client"] = self.client.__repr__()
	return ret
	def __setstate__(self, data):
	# somewhat weird that this import is requiired here but it won't
	# work without it. Not sure how the symbol MongoClient is required
	# here but it is - ignore if a lint like ide says MongoClient is not used
	from pymongo import MongoClient
	data["_Database__client"] = eval(data["_Database__client"])
	self.__dict__.update(data)

Therefore, I am not quite sure where does the PyCapsule come from. btw, if it is indeed from the database handler, I remember the error message is different, and it should relate to the Client object instead of PyCapsule. The only place I found that used py::capsule before is here:

mspass/cxx/python/utility/utility_py.cc

Lines 201 to 209 in 85c8f46

	/* The use of capsule to avoid copy is found at
	https://github.com/pybind/pybind11/issues/1042 */
	.def_property_readonly("unit_vector", [](const SphericalCoordinate& self) {
	double* v = SphericalToUnitVector(self);
	// auto capsule = py::capsule(v, [](void *v) { delete[] reinterpret_cast<double*>(v); });
	// return py::array(3, v, capsule);
	vector<double> unit_vector(v, v + 3);
	delete[] v;
	return pybind11::module::import("numpy").attr("array")(unit_vector);

where you can clearly see that the relevant code is already commented out with a new numpy implementation. I guess the only way to sort out the problem is to have me run and debug the code. We attempted this before but was blocked by a separate issue. Maybe it is time to reinitiate that process.

[pavlis]

Seems my hypothesis may well be false then. As noted I was guessing because I couldn't untangle the messy data structure dask was passing through the scheduler. I learned a lot about how dask works by doing that process with spyder, but it was complicated enough I got lost in the woods, stepped back and looked at what was being done, and made the conjecture that Database was to blame even though I couldn't actually even see any references to Database in the data I was viewing with spyder. It seems my conjecture was likely wrong. Based on what you are saying I will dig back into this before punting it to you.

FYI, I did learn a key thing in working on this to make the "PyCapsule" reference a bit less confusing. One source I read (no longer can easily find it) described it as a way to push around an opaque pointer. In C that approach is abused way too often and used even in stdlib for things malloc where a void* pointer has to be cast into whatever the code wants it to be. A source, by the way, of a large number of software bugs. What is mysterious is what is being miscast in this code into a PyCapsule? I certainly didn't sort it out yesterday, but let me have another go at it before I punt it to you.

Putting aside the needs for fixing pwstack/pwmig, this still raises design and documentation issues for MsPASS. We really really need to have a clear set of documentation on the whole serialization issue. Seems to me that is a huge weakness in both dask and spark to not have a common tool for testing the feasibility of serializing a random python object. I suspect strongly it is because the dominant use of these packages are for things like dataframes where what is atomic is very rigidly defined relative to a bag/rdd where the contents can be anything. Our use, I think, is particularly fringe because we have all these custom classes bound to python with pybind11. Not sure what is realistic to support in MsPASS, but an absolute mininum is more extensive documentation on parallelizing a workflow. The one draft page we have is far far from sufficient, especially since the parallel thing is one of the key innovations we are advertising for this package.

[pavlis]

New hypothesis after painstaking stepping through dask scheduler code. It seems likely your guess that read_ensembles was the problem child. If my guess is correct, we do have an issue with obspy. This will get little verbose because I'm going to paste in a bunch of code to make it clear what I think the issue is here. This is the read_ensembles function that I'm now suspecting may be the problem:

def read_ensembles(db,querydata,control):
    """
    Constructs a query from dict created by build_wfquery, runs it
    on the wf_Seismogram collection of db, and then calls read_ensemble_data
    on the cursor returned by MongoDB.  It the sets the ensemble
    metadata for lat, lon, ix1, and ix2 before returning the ensembles

    We need to pass the control class to access the stack_count_cutoff attribute.
    """
    # don't even issue a query if the fold is too low
    fold=querydata['fold']
    if fold<=control.stack_count_cutoff:
        d=SeismogramEnsemble()
    else:
        query=querydata['query']
        n=db.wf_Seismogram.count_documents(query)
        #debug
        #print(query,n)
        if n==0:
            # This shouldn't ever really be executed unless stack_count_cutoff
            # is 0 or the query is botched
            d=SeismogramEnsemble()
        else:
            cursor=db.wf_Seismogram.find(query)
            # Note control.data_tag can be a None type here - see
            # control object constructor
            d=db.read_ensemble_data(cursor,collection='wf_Seismogram',
                                    normalize=['source','site'],
                                    data_tag=control.data_tag)
        if len(d.member) > 0:
            # When the ensemble is not empty we have to compute the
            # slowness vector of the incident wavefield using source \
            # coordinates and the pseudostation location.  This section
            # depends on a feature of the reader that it the normalize
            # parameter causes all members to have source coordinates loaded
            # this small function copies source metadata from the first
            # live member
            srcdata = get_source_metadata(d)
            if len(srcdata)==0:
                d.kill()
                d.elog.log_error("pwstack",
                        "Ensembled read from database has no live members",
                        ErrorSeverity.Invalid)
            else:
                # post the source metadata to the ensemble metadata
                for k in srcdata:
                    d[k] = srcdata[k]
                # a bit of a weird way to fetch the pseudostation
                # coordinates but the only way without using additional arts
                pslat=querydata['lat']
                pslon=querydata['lon']
                georesult=gps2dist_azimuth(srcdata['source_lat'],
                                        srcdata['source_lon'],pslat,pslon)
                # obspy's function we just called returns distance in m in element 0 of a tuple
                # their travel time calculator it is degrees so we need this conversion
                Rearth=6378.164
                dist=kilometers2degrees(georesult[0]/1000.0)
                # We pass the model object through control because I think
                # there is a nontrivial overhead in creating it
                arrivals=control.model.get_travel_times(
                    source_depth_in_km=srcdata['source_depth'],
                    distance_in_degree=dist,phase_list=['P']
                    )
                ray_param=arrivals[0].ray_param
                umag=ray_param/Rearth    # need slowness in s/km but ray_param is s/radian
                baz=georesult[2]   # The obspy function seems to return back azimuth
                az=baz+180.0
                # az can be > 360 here but all known trig function algs handl this automatically
                ux=umag*math.sin(math.radians(az))
                uy=umag*math.cos(math.radians(az))
                d.put('ux0', ux)
                d.put('uy0', uy)
                # We intentionally save the same attributes with
                # different keys.  One has units of degrees (the two
                # with the pseudostation tag) and the others have
                # units of radians.   The quantities in radian are
                # required for the C++ pwstack_ensemble function.
                d.put('pseudostation_lat',pslat)
                d.put('pseudostation_lon',pslon)
                d.put('lat0',math.radians(pslat))
                d.put('lon0',math.radians(pslon))
                d.put('ix1',querydata['ix1'])
                d.put('ix2',querydata['ix2'])
                d.put('gridname',control.pseudostation_gridname)
    return d

I'm most suspicious of this line in that functions:

                arrivals=control.model.get_travel_times(
                    source_depth_in_km=srcdata['source_depth'],
                    distance_in_degree=dist,phase_list=['P']
                    )

The problem, I think, is that the taup travel time calculator is passed to this function through the python class received as the symbol control. control.model is a taup calculator handle created like this:

self.model = TauPyModel(model=modelname)

Here is dir of it:

IPdb [3]: dir(self.model)
['__class__', '__delattr__', '__dict__', '__dir__', '__doc__', '__eq__', '__format__', '__ge__', '__getattribute__', '__gt__', '__hash__', '__init__', '__init_subclass__', '__le__', '__lt__', '__module__', '__ne__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__slotnames__', '__str__', '__subclasshook__', '__weakref__', 'get_pierce_points', 'get_pierce_points_geo', 'get_ray_paths', 'get_ray_paths_geo', 'get_travel_times', 'get_travel_times_geo', 'model', 'planet_flattening', 'verbose']

Because control is an argument received by read_ensembles that I'm trying to run through a map call I'm high suspicious this is causing the PyCapsule error. The reason I suspect that function is that when I've used that function before I often get mysterious warnings from TauPyModel about initialization issues that suggest strongly it has some odd state dependency. Since my last tests hinted strongly that read_ensembles was the problem and that is the only part of the function that seems prone to any odd behavior I think that is a reasonable hypothesis.

I'm going to attempt a workaround of this by divorcing the travel time calculations from dask and inserting the travel times as metadata into the input data. Then I can remove the travel time functions from read_ensembles to test the hypothesis that the travel time calculator is causing this error.

If my current hypothesis proves true that is a VERY important warning to our users as travel time calculations are a very very common requirement for seismology algortihms

[pavlis]

It is starting to look like my hypothesis about obspy's travel time calculator being a problem is an issue. It turned out my original idea was not going to be easy to do and would not have been very definitive anyway. I realized it was easier to write a pure test program that did more or less what the large function was doing that I suspected was trouble. So, I wrote this little test program I would like both of you to see if you can run:

from obspy.taup import TauPyModel
from obspy.geodetics import gps2dist_azimuth,kilometers2degrees
import numpy as np
import math

import dask.bag as db
import dask

def ttcompute(dist,model):
    arrivals=model.get_travel_times(source_depth_in_km=10.0,distance_in_degree=dist,
                                    phase_list=['P'])
    ray_param=arrivals[0].ray_param
    Rearth=6378.164
    umag=ray_param/Rearth    # need slowness in s/km but ray_param is s/radian
    baz=0.0
    az=baz+180.0
    # az can be > 360 here but all known trig function algs handl this automatically
    ux=umag*math.sin(math.radians(az))
    uy=umag*math.cos(math.radians(az))
    return [ux,uy]

# This one is ran to completetion without any issues.
dask.config.set(scheduler='processes',num_workers=2) 
#dask.config.set(scheduler='synchronous')

model = TauPyModel(model='ak135')

# Drive dask with a sequence of epicentral distances
distance = np.linspace(1.0,90.0,1000)
#for d in distance:
#    x=ttcompute(d,model)
#    print(x)
bg = db.from_sequence(distance)
bg = bg.map(ttcompute,model)
x=bg.compute()
print(x)

When I run this as above I get this chain of errors:

runfile('/home/pavlis/data/copy_pwmigtest/test_dask_taup.py', wdir='/home/pavlis/data/copy_pwmigtest')
Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "/home/pavlis/anaconda3/lib/python3.8/multiprocessing/spawn.py", line 116, in spawn_main
    exitcode = _main(fd, parent_sentinel)
  File "/home/pavlis/anaconda3/lib/python3.8/multiprocessing/spawn.py", line 125, in _main
    prepare(preparation_data)
  File "/home/pavlis/anaconda3/lib/python3.8/multiprocessing/spawn.py", line 236, in prepare
    _fixup_main_from_path(data['init_main_from_path'])
  File "/home/pavlis/anaconda3/lib/python3.8/multiprocessing/spawn.py", line 287, in _fixup_main_from_path
    main_content = runpy.run_path(main_path,
  File "/home/pavlis/anaconda3/lib/python3.8/runpy.py", line 265, in run_path
    return _run_module_code(code, init_globals, run_name,
  File "/home/pavlis/anaconda3/lib/python3.8/runpy.py", line 97, in _run_module_code
    _run_code(code, mod_globals, init_globals,
  File "/home/pavlis/anaconda3/lib/python3.8/runpy.py", line 87, in _run_code
    exec(code, run_globals)
  File "/home/pavlis/data/copy_pwmigtest/test_dask_taup.py", line 44, in <module>
    x=bg.compute()
  File "/home/pavlis/anaconda3/lib/python3.8/site-packages/dask/base.py", line 288, in compute
    (result,) = compute(self, traverse=False, **kwargs)
  File "/home/pavlis/anaconda3/lib/python3.8/site-packages/dask/base.py", line 570, in compute
    results = schedule(dsk, keys, **kwargs)
  File "/home/pavlis/anaconda3/lib/python3.8/site-packages/dask/multiprocessing.py", line 219, in get
    result = get_async(
  File "/home/pavlis/anaconda3/lib/python3.8/site-packages/dask/local.py", line 505, in get_async
    fire_tasks(chunksize)
  File "/home/pavlis/anaconda3/lib/python3.8/site-packages/dask/local.py", line 500, in fire_tasks
    fut = submit(batch_execute_tasks, each_args)
  File "/home/pavlis/anaconda3/lib/python3.8/concurrent/futures/process.py", line 645, in submit
    self._start_queue_management_thread()
  File "/home/pavlis/anaconda3/lib/python3.8/concurrent/futures/process.py", line 584, in _start_queue_management_thread
    self._adjust_process_count()
  File "/home/pavlis/anaconda3/lib/python3.8/concurrent/futures/process.py", line 608, in _adjust_process_count
    p.start()
  File "/home/pavlis/anaconda3/lib/python3.8/multiprocessing/process.py", line 121, in start
    self._popen = self._Popen(self)
  File "/home/pavlis/anaconda3/lib/python3.8/multiprocessing/context.py", line 284, in _Popen
    return Popen(process_obj)
  File "/home/pavlis/anaconda3/lib/python3.8/multiprocessing/popen_spawn_posix.py", line 32, in __init__
    super().__init__(process_obj)
  File "/home/pavlis/anaconda3/lib/python3.8/multiprocessing/popen_fork.py", line 19, in __init__
    self._launch(process_obj)
  File "/home/pavlis/anaconda3/lib/python3.8/multiprocessing/popen_spawn_posix.py", line 42, in _launch
    prep_data = spawn.get_preparation_data(process_obj._name)
  File "/home/pavlis/anaconda3/lib/python3.8/multiprocessing/spawn.py", line 154, in get_preparation_data
    _check_not_importing_main()
  File "/home/pavlis/anaconda3/lib/python3.8/multiprocessing/spawn.py", line 134, in _check_not_importing_main
    raise RuntimeError('''
RuntimeError: 
        An attempt has been made to start a new process before the
        current process has finished its bootstrapping phase.

        This probably means that you are not using fork to start your
        child processes and you have forgotten to use the proper idiom
        in the main module:

            if __name__ == '__main__':
                freeze_support()
                ...

        The "freeze_support()" line can be omitted if the program
        is not going to be frozen to produce an executable.
Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "/home/pavlis/anaconda3/lib/python3.8/multiprocessing/spawn.py", line 116, in spawn_main
    exitcode = _main(fd, parent_sentinel)
  File "/home/pavlis/anaconda3/lib/python3.8/multiprocessing/spawn.py", line 125, in _main
    prepare(preparation_data)
  File "/home/pavlis/anaconda3/lib/python3.8/multiprocessing/spawn.py", line 236, in prepare
    _fixup_main_from_path(data['init_main_from_path'])
  File "/home/pavlis/anaconda3/lib/python3.8/multiprocessing/spawn.py", line 287, in _fixup_main_from_path
    main_content = runpy.run_path(main_path,
  File "/home/pavlis/anaconda3/lib/python3.8/runpy.py", line 265, in run_path
    return _run_module_code(code, init_globals, run_name,
  File "/home/pavlis/anaconda3/lib/python3.8/runpy.py", line 97, in _run_module_code
    _run_code(code, mod_globals, init_globals,
  File "/home/pavlis/anaconda3/lib/python3.8/runpy.py", line 87, in _run_code
    exec(code, run_globals)
  File "/home/pavlis/data/copy_pwmigtest/test_dask_taup.py", line 44, in <module>
    x=bg.compute()
  File "/home/pavlis/anaconda3/lib/python3.8/site-packages/dask/base.py", line 288, in compute
    (result,) = compute(self, traverse=False, **kwargs)
  File "/home/pavlis/anaconda3/lib/python3.8/site-packages/dask/base.py", line 570, in compute
    results = schedule(dsk, keys, **kwargs)
  File "/home/pavlis/anaconda3/lib/python3.8/site-packages/dask/multiprocessing.py", line 219, in get
    result = get_async(
  File "/home/pavlis/anaconda3/lib/python3.8/site-packages/dask/local.py", line 505, in get_async
    fire_tasks(chunksize)
  File "/home/pavlis/anaconda3/lib/python3.8/site-packages/dask/local.py", line 500, in fire_tasks
    fut = submit(batch_execute_tasks, each_args)
  File "/home/pavlis/anaconda3/lib/python3.8/concurrent/futures/process.py", line 645, in submit
    self._start_queue_management_thread()
  File "/home/pavlis/anaconda3/lib/python3.8/concurrent/futures/process.py", line 584, in _start_queue_management_thread
    self._adjust_process_count()
  File "/home/pavlis/anaconda3/lib/python3.8/concurrent/futures/process.py", line 608, in _adjust_process_count
    p.start()
  File "/home/pavlis/anaconda3/lib/python3.8/multiprocessing/process.py", line 121, in start
    self._popen = self._Popen(self)
  File "/home/pavlis/anaconda3/lib/python3.8/multiprocessing/context.py", line 284, in _Popen
    return Popen(process_obj)
  File "/home/pavlis/anaconda3/lib/python3.8/multiprocessing/popen_spawn_posix.py", line 32, in __init__
    super().__init__(process_obj)
  File "/home/pavlis/anaconda3/lib/python3.8/multiprocessing/popen_fork.py", line 19, in __init__
    self._launch(process_obj)
  File "/home/pavlis/anaconda3/lib/python3.8/multiprocessing/popen_spawn_posix.py", line 42, in _launch
    prep_data = spawn.get_preparation_data(process_obj._name)
  File "/home/pavlis/anaconda3/lib/python3.8/multiprocessing/spawn.py", line 154, in get_preparation_data
    _check_not_importing_main()
  File "/home/pavlis/anaconda3/lib/python3.8/multiprocessing/spawn.py", line 134, in _check_not_importing_main
    raise RuntimeError('''
RuntimeError: 
        An attempt has been made to start a new process before the
        current process has finished its bootstrapping phase.

        This probably means that you are not using fork to start your
        child processes and you have forgotten to use the proper idiom
        in the main module:

            if __name__ == '__main__':
                freeze_support()
                ...

        The "freeze_support()" line can be omitted if the program
        is not going to be frozen to produce an executable.
Traceback (most recent call last):

  File "/home/pavlis/data/copy_pwmigtest/test_dask_taup.py", line 44, in <module>
    x=bg.compute()

  File "/home/pavlis/anaconda3/lib/python3.8/site-packages/dask/base.py", line 288, in compute
    (result,) = compute(self, traverse=False, **kwargs)

  File "/home/pavlis/anaconda3/lib/python3.8/site-packages/dask/base.py", line 570, in compute
    results = schedule(dsk, keys, **kwargs)

  File "/home/pavlis/anaconda3/lib/python3.8/site-packages/dask/multiprocessing.py", line 219, in get
    result = get_async(

  File "/home/pavlis/anaconda3/lib/python3.8/site-packages/dask/local.py", line 506, in get_async
    for key, res_info, failed in queue_get(queue).result():

  File "/home/pavlis/anaconda3/lib/python3.8/concurrent/futures/_base.py", line 432, in result
    return self.__get_result()

  File "/home/pavlis/anaconda3/lib/python3.8/concurrent/futures/_base.py", line 388, in __get_result
    raise self._exception

BrokenProcessPool: A process in the process pool was terminated abruptly while the future was running or pending.

Note if you change the dask.config.set line to use scheduler='syncrhonous' it should run (did for me). Output is meaningless but is correct. Note you will also get the same answer if you run the for loop version.

Now what happens here is a completely different problem from what I was getting with the pwstack program, but seems to suggest an issue with the TauPyModel class.

Could you guys verify this behaves the same for you? I ran this from my local standalone machine not from the docker container. We need to rule out the issue being something with my local version of dask. (No reason to think that but this is not our standard container.)

[pavlis]

As recommended in our last meeting I tried running that test script above with our current container. I did not get the error above and it ran correctly. Hence, my hypothesis of taup being an issue was most definitely wrong.

Seems I must have an installation problem with my local copy of dask. It is conceivable the reason I'm getting the original error trying to get my pwstack program to run is related to that. Need some help sorting that out. Perhaps a topic for our next group meeting.