How to use the mpi4py.MPI.COMM_WORLD.Allreduce function in mpi4py

def _mpi_average(self, x):
        buf = np.zeros_like(x)
        MPI.COMM_WORLD.Allreduce(x, buf, op=MPI.SUM)
        buf /= MPI.COMM_WORLD.Get_size()
        return buf

def update(self, x):
        x = x.astype('float64')
        n = int(np.prod(self.shape))
        totalvec = np.zeros(n*2+1, 'float64')
        addvec = np.concatenate([x.sum(axis=0).ravel(), np.square(x).sum(axis=0).ravel(), np.array([len(x)],dtype='float64')])
        MPI.COMM_WORLD.Allreduce(addvec, totalvec, op=MPI.SUM)
        self.incfiltparams(totalvec[0:n].reshape(self.shape), totalvec[n:2*n].reshape(self.shape), totalvec[2*n])

def allreduce(*args, **kwargs):
    return MPI.COMM_WORLD.Allreduce(*args, **kwargs)

def update(self, x):
        x = x.astype('float64')
        n = int(np.prod(self.shape))
        totalvec = np.zeros(n*2+1, 'float64')
        addvec = np.concatenate([x.sum(axis=0).ravel(), np.square(x).sum(axis=0).ravel(), np.array([len(x)],dtype='float64')])
        if MPI is not None:
            MPI.COMM_WORLD.Allreduce(addvec, totalvec, op=MPI.SUM)
        self.incfiltparams(totalvec[0:n].reshape(self.shape), totalvec[n:2*n].reshape(self.shape), totalvec[2*n])

def _meshAdvector(self, tectonicX, tectonicY, timer):
        """
        Advect the mesh horizontally and interpolate mesh information
        """

        # Move horizontally coordinates
        XYZ = np.zeros((self.gpoints,3))
        XYZ[:,0] = self.gcoords[:,0] + tectonicX*timer
        XYZ[:,1] = self.gcoords[:,1] + tectonicY*timer

        # Get mesh variables to interpolate
        # Elevation
        elev = np.zeros(self.gpoints)
        elev.fill(-1.e8)
        elev[self.natural2local] = self.hLocal.getArray().copy()
        MPI.COMM_WORLD.Allreduce(MPI.IN_PLACE, elev, op=MPI.MAX)
        interpolator = CloughTocher2DInterpolator(XYZ[:,:2],elev)
        nelev = interpolator(self.lcoords[:,:2])
        id_NaNs = np.isnan(nelev)
        nelev[id_NaNs] = 0.

        # Erosion/deposition
        erodep = np.zeros(self.gpoints)
        erodep.fill(-1.e8)
        erodep[self.natural2local] = self.cumEDLocal.getArray().copy()
        MPI.COMM_WORLD.Allreduce(MPI.IN_PLACE, erodep, op=MPI.MAX)
        interpolator = CloughTocher2DInterpolator(XYZ[:,:2],erodep)
        nerodep = interpolator(self.lcoords[:,:2])
        nerodep[id_NaNs] = 0.

        # Soil thickness
        if self.Ksed > 0.:

def allmean(x):
        assert isinstance(x, np.ndarray)
        out = np.empty_like(x)
        MPI.COMM_WORLD.Allreduce(x, out, op=MPI.SUM)
        out /= nworkers
        return out

def allmean(x):
        assert isinstance(x, np.ndarray)
        if MPI is not None:
            out = np.empty_like(x)
            MPI.COMM_WORLD.Allreduce(x, out, op=MPI.SUM)
            out /= nworkers
        else:
            out = np.copy(x)

        return out

def mpi_moments(x, axis=0):
    x = np.asarray(x, dtype='float64')
    newshape = list(x.shape)
    newshape.pop(axis)
    n = np.prod(newshape,dtype=int)
    totalvec = np.zeros(n*2+1, 'float64')
    addvec = np.concatenate([x.sum(axis=axis).ravel(), 
        np.square(x).sum(axis=axis).ravel(), 
        np.array([x.shape[axis]],dtype='float64')])
    MPI.COMM_WORLD.Allreduce(addvec, totalvec, op=MPI.SUM)
    sum = totalvec[:n]
    sumsq = totalvec[n:2*n]
    count = totalvec[2*n]
    if count == 0:
        mean = np.empty(newshape); mean[:] = np.nan
        std = np.empty(newshape); std[:] = np.nan
    else:
        mean = sum/count
        std = np.sqrt(np.maximum(sumsq/count - np.square(mean),0))
    return mean, std, count

def initAMRData(self,gridData):
        ngrids=len(gridData['gridParam'])
        local2global_tmp=np.zeros((ngrids,),'i')
        local2global=np.zeros((ngrids,),'i')
        for i in range(len(gridData['qParam'])):
            local2global_tmp[gridData['qParam'][i][0]]=i
        MPI.COMM_WORLD.Allreduce(local2global_tmp,local2global)
        #
        # these variables are from high-order FE off-body implementation
        # have to see if p=0 defaults work
        #
        nf=3
        qstride=1
        qnodein=0.0
        qnodeinC=arrayToDblPtr(np.array([qnodein,0.0],'d'))
        qnodesize=1

        idata=np.zeros((ngrids*11),'i')
        rdata=np.zeros((ngrids*6),'d')
        for i in range(ngrids):
            m=i*11
            idata[m]   =gridData['gridParam'][i][0] # global id of patch
            idata[m+1] =gridData['gridParam'][i][1] # level number of patch

def allmean(x):
        assert isinstance(x, np.ndarray)
        out = np.empty_like(x)
        MPI.COMM_WORLD.Allreduce(x, out, op=MPI.SUM)
        out /= nworkers
        return out