MEDCoupling,  multiprocessing
----------------------------

Cet exercice fait la supposition de Numpy Scipy sont correctement maĂŽtrisĂŠs, sinon voir :ref:`medcouplingnumpyptr`.
On va faire simuler un traitement un peu long (ici l'interpolation d'un maillage de 64000 cells avec un autre de 64000 cells).
On va faire le traitement d'abord en sÊquentiel puis en parallèle pour exploiter les coeurs de notre CPU.
Nous allons utiliser le module ``multiprocessing`` pour cela.

DĂŠbut de l'implĂŠmentation
~~~~~~~~~~~~~~~~~~~~~~~~~

Pour commencer l'exercice importer le module Python ``MEDCoupling``, ``MEDCouplingRemapper``, ``numpy``, ``scipy``, ``multiprocessing``
et ``datetime`` pour chronomĂŠtrer : ::

	import MEDCoupling as mc
	import MEDCouplingRemapper as mr
	import multiprocessing as mp
	from datetime import datetime
	from scipy.sparse import csr_matrix

CrĂŠer un maillage cartĂŠsien rĂŠgulier 3D avec 40 cells en X, Y et Z
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
CrĂŠons un maillage cartĂŠsien 3D m de pas rĂŠgulier entre 0. et 1. en X, Y et Z : ::

	nbCells=40
	arr=mc.DataArrayDouble(nbCells+1) ; arr.iota() ; arr/=nbCells
	m=mc.MEDCouplingCMesh() ; m.setCoords(arr,arr,arr)

Traduisons m en non structurĂŠ pour que le calcul soit plus long : ::

	m=m.buildUnstructured()

CrĂŠer une copie m2 de m translatĂŠe de la moitiĂŠ du pas en X, Y et Z ::

	m2=m.deepCopy()
	t=mc.DataArrayDouble(3)
	t[:]=1/(2*float(nbCells))
	m2.translate(t.getValues())

Calculer sĂŠquentiellement la matrice d'interpolation M de la projection entre m et m2 en P0P0
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

m sera considĂŠrĂŠ comme le maillage source et m2 sera considĂŠrĂŠ comme le maillage cible.
Profitons en pour chronomĂŠtrer le temps necessaire pour le traitement sĂŠquentiel.
Utilisons ``MEDCouplingRemapper`` pour cela. ::

	remap=mr.MEDCouplingRemapper()
	strt=datetime.now()
	assert(remap.prepare(m,m2,"P0P0")==1)
	print "time in sequential : %s"%(str(datetime.now()-strt))

Stockons la sparse matrix scipy dans ``matSeq``. ::

	matSeq=remap.getCrudeCSRMatrix()

Calculer cette même matrice M en parallèle avec multiprocessing.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Commencons par rĂŠcupĂŠrer le nombre de coeur de notre machine. ::

	nbProc=mp.cpu_count()

L'idĂŠe est de faire une mĂŠthode ``work`` prenant un tuple de longueur 2.
Le premier ĂŠlĂŠment du tuple est la partie du maillage ``m2`` considĂŠrĂŠe. Le 2eme ĂŠlĂŠment donne la correspondance entre les cells id de ``m2Part`` les cells id de ``m2``.

L'idĂŠe est d'interpoler ``m`` avec ``m2Part``.

On rÊcupèrera ensuite la matrice sparse ``myMat`` issue de ``m`` avec ``m2Part``.
Ensuite l'idĂŠe et de gĂŠnĂŠrer une matrice sparse ``mat2`` Ă  partir de ``myMat`` avec les ids globaux de ``m2``. ::

	def work(inp):
            m2Part,partToGlob=inp
	    myRemap=mr.MEDCouplingRemapper()
	    assert(myRemap.prepare(m,m2Part,"P0P0")==1)
	    myMat=myRemap.getCrudeCSRMatrix()
	    a=mc.DataArrayInt.Range(s.start,s.stop,s.step)
	    indptrnew=mc.DataArrayInt(m2.getNumberOfCells())
	    indptrnew.fillWithZero()
	    d=mc.DataArrayInt(myMat.indptr).deltaShiftIndex()
	    indptrnew[partToGlob]=d
	    indptrnew.computeOffsetsFull()
	    mat2=csr_matrix( (myMat.data,myMat.indices,indptrnew.toNumPyArray()), shape=(m2.getNumberOfCells(),m.getNumberOfCells()))
	    return mat2

Il s'agit dÊsormais de faire la liste des inputs à donner aux ``nbProc`` instances de ``work`` qui seront exÊcutÊs en parallèle.
Appelons cette liste python ``workToDo`` qui sera de longueur ``nbProc``.
On peut se faire aider de ``mc.DataArray.GetSlice``. ::

        workToDo=[]
        for i in xrange(nbProc):
              s=mc.DataArray.GetSlice(slice(0,m2.getNumberOfCells(),1),i,nbProc)
              part=m2[s]
              partToGlob=mc.DataArrayInt.Range(s.start,s.stop,s.step)
              workToDo.append((part,partToGlob))
              pass

On est maintenant prĂŞt pour faire travailler chacun des coeurs indĂŠpendamment. Pour ce faire, on crĂŠe un ``mp.Pool`` et on assigne Ă  chaque worker le travail ``work`` avec autant de worker que de coeurs. Et chronomĂŠtrons le tout ! ::


	strt=datetime.now()
	pool = mp.Pool()
	asyncResult = pool.map_async(work,workToDo)
	subMatrices = asyncResult.get()
	print "time in parallel (x%d) : %s"%(nbProc,str(datetime.now()-strt))

.. note:: A noter la magie ! On a transfĂŠrĂŠ entre le process maitre et les process esclave sans mĂŞme s'en rendre compte les maillages et les DataArrayInt contenus dans ``workToDo`` !
	  Merci Ă  la pickelisation des objets MEDCoupling :)

VĂŠrfication
~~~~~~~~~~~

VĂŠrifions que les matrices sont les mĂŞmes ! Sommons ``subMatrices`` (``matPar``) et regardons le nombre de non zĂŠros de la diffĂŠrence entre la ``matPar`` et ``matSeq``. ::

	matPar = sum(subMatrices)
	matDelta=matSeq-matPar
	assert(matDelta.nnz==0)