source: codes/icosagcm/devel/Python/dynamico/parallel.py

from itertools import groupby
import numpy as np

def inverse_list(lst): return {j:i for i,j in enumerate(lst)}
def ordered_list(lst): 
        lst = sorted(zip( lst, range(len(lst)) ))
        lst, order = zip(*lst)
        order = inverse_list(order)
        order = [order[i] for i in range(len(lst))] # dict->list
        return lst,order  

from dynamico import getargs
log_master, log_world = getargs.getLogger(__name__)
INFO, DEBUG, ERROR = log_master.info, log_master.debug, log_world.error
INFO_ALL, DEBUG_ALL = log_world.info, log_world.debug

#----------------------- Parallel access to distributed array ------------------------------#

# The classes PDim, Get_Indices and PArrayX are used to read NetCDF arrays in parallel
# Each MPI process reads a contiguous chunk of the array once.
# Classes CstPArrayX emulate constant arrays that are not stored in a NetCDF file.
# With class LocArray1D, local data is copied from a user-provided array, rather than from a NetCDF file.
# Using the class Get_Indices, each MPI process can obtain data for an arbitrary list of indices
# Data is not re-read from file but scattered via MPI all2all

class PDim:
    def __init__(self, ncdim, comm):
        mpi_rank, mpi_size = comm.Get_rank(), comm.Get_size()
        self.ncdim, n = ncdim, len(ncdim)
        vtxdist = [(n*i)/mpi_size for i in range(mpi_size+1)]
        self.n, self.vtxdist = n, np.asarray(vtxdist, dtype=np.int32)
        self.comm, self.start, self.end = comm, vtxdist[mpi_rank], vtxdist[mpi_rank+1]
    def getter(self, index_list): return Get_Indices(self, index_list)
    def get(self, index_list, data):
        getter=self.getter(index_list)
        return getter(data)
class Get_Indices:
    def __init__(self, dim, index_list): # list MUST be a list of unique indices
        index_list, self.order = ordered_list(index_list)
        comm, vtxdist = dim.comm, dim.vtxdist
        self.comm, self.n, self.vtxdist = comm, len(index_list), vtxdist
        list_local = [ [ j-vtxdist[i] for j in index_list if j>=vtxdist[i] and j<vtxdist[i+1]] 
                      for i in range(len(vtxdist)-1) ]
        self.list_local, self.list_n = list_local, map(len, list_local)
        list_send = comm.alltoall(list_local) # may modify list_local !!
        self.list_send, self.dict = list_send, inverse_list(index_list)
    def get(self, get_data, put_data):
        data_send = [get_data(indices) for indices in self.list_send]
        data_recv = self.comm.alltoall(data_send)
        start=0
        for data,n in zip(data_recv,self.list_n):
            if n>0:
                end = start + n
                put_data(start,end,data)
                start=end
    def __call__(indices, self): # self is a PArrayND
        shape = list(self.data.shape)
        shape[0] = indices.n
        self.data_out = np.zeros(shape, dtype=self.data.dtype)
        indices.get(self.get_data, self.put_data )
#        return self.data_out
        return self.reorder(indices.order)

class PArrayND:
    def init_data(self, dim, data):
        self.dim, self.data = dim, np.asarray(data) # local chunk of data
    def max(self):
        max_loc = self.data.max()
        self.dim.comm.allreduce

class PArray1D(PArrayND):
    def __init__(self, dim, data): self.init_data(dim, data[dim.start:dim.end])
    def get_data(self, indices): return np.array(self.data[indices])
    def put_data(self, start, end, data): self.data_out[start:end]=data
    def reorder(self, order): return self.data_out[order]
class LocPArray1D(PArray1D):
    def __init__(self, dim, data): self.init_data(dim,data)
class CstPArray1D(PArray1D):
    def __init__(self, dim, dtype, val): 
        self.init_data(dim,np.full( (dim.end-dim.start,), val, dtype=dtype))

class PArray2D(PArrayND):
    def __init__(self, dim, data): self.init_data(dim, data[dim.start:dim.end,:])
    def get_data(self, indices): return self.data[indices,:]
    def put_data(self, start, end, data): self.data_out[start:end]=data
    def reorder(self, order): return self.data_out[order,:]

def PArray(dim, data): return {1:PArray1D, 2:PArray2D}[len(data.shape)](dim,data)

#-------------------------------------- Halo management -----------------------------------#

# Classes LocalDim, LocalArrayX and Halo_Xchange are used to 
# store local arrays with halos and update halos, respectively
# A LocalDim instance is created using a list of (global indices of) cells. The instance contains
# a lookup table associating a local index to a global index. 
# A Halo_Xchange instance is created first based on a LocalDim and a partitioning table giving MPI process owning each cell.
# This instance can the be used to create instances of Halo_Xchange. 
# LocalArray data can be initialized by reading from a PArray or a NumPy array containing local values. 
# The update() methods updates the halo using point-to-point MPI communications (send/recv) 

class LocalDim:
    def __init__(self, dim, cells):
        # dim : a PDim instance ; cells : a list of (global indices of) cells, in any order
        self.dim, self.loc2glob, self.glob2loc = dim, cells, inverse_list(list)
        
def dict2list(size, thedict, default):
    lst=[default]*size
    for rank,val in thedict.items() : lst[rank]=val
    return lst

class Halo_Xchange:
    def __init__(self, tag, dim, cells, part, reorder_cells=False): 
        # cells = global indices, part = rank owning each cell
        # if reorder_cells is True, cells will be reordered so that each halo is contiguous
        # reordering is not recommended if cells is a telescopic sum of increasing sets C0<C1<C2 ...
        comm = dim.comm
        mpi_size, mpi_rank = comm.Get_size(), comm.Get_rank()
        self.comm, self.tag, self.reorder_cell = comm, tag, reorder_cells
        # sort cells by rank, keep track of their global and local index
        recv = sorted(zip(part,cells,range(len(cells))))
        # group by MPI rank ; lst is a list of (rank,global,local) tuples so zip(*it) is a tuple of lists (rank,global,local)
        recv = { rank : zip(*lst) for rank,lst in groupby(recv,lambda x:x[0]) }
        # remove mpi_rank from dict and get list of own cells
        junk, own_global, own_local = recv.pop(mpi_rank)
        own_len = len(own_local)
        if reorder_cells :  own_local = range(own_len)
        self.own_len, self.own_local, self.get_own = own_len, list(own_local), Get_Indices(dim, own_global) 
        # figure out data we want to receive from other CPUs
        recv_rank = sorted(recv.keys())
        recv_glob = { rank : recv[rank][1] for rank in recv_rank} # global indices of cells to receive
        recv_len = { rank : len(recv_glob[rank]) for rank in recv_rank}
        if reorder_cells :
            recv_loc, start = {}, own_len
            for rank in recv_rank:
                end = start+recv_len[rank]
                recv_loc[rank] = range(start,end)
                start=end
            cells = sum([recv_glob[rank] for rank in recv_rank], own_global)
        else:
            recv_loc = { rank : recv[rank][2] for rank in recv_rank}
        self.recv_list = [(rank, list(recv_loc[rank])) for rank in recv_rank]
        self.cells, self.get_all = cells, Get_Indices(dim, cells)
        # now figure out the data we must send
        send=comm.alltoall(dict2list(mpi_size, recv_glob, []))
        send_len = [len(lst) for lst in send]
        send_rank = [rank for rank in range(mpi_size) if (send_len[rank]>0)] # CPUs asking for data
        DEBUG('send_rank %d %s'%(mpi_rank, send_rank))
        # asssociate local index to global index
        own_dict = { glob:loc for glob,loc in zip(own_global, own_local) }
        self.send_list = [(rank, [own_dict[i] for i in send[rank]]) for rank in send_rank]
    def set_dynamico_transfer(self,name):
        def list_to_args(lst):
            ranks, lst = ([], []) if lst == [] else zip(*lst) # list of tuples => tuple of lists
            lens = map(len, lst)
            return len(ranks), ranks, lens, sum(lens), sum(lst,[])
        index = { 'primal':1, 'edge':2, 'dual':3 }[name]
        self.index = index
        send_num, send_rank, send_len, send_size, send_list = list_to_args(self.send_list) 
        recv_num, recv_rank, recv_len, recv_size, recv_list = list_to_args(self.recv_list) 
        send_rank, send_len, send_list, recv_rank, recv_len, recv_list = [ np.asarray(x, dtype=np.int32)
            for x in send_rank, send_len, send_list, recv_rank, recv_len, recv_list ]
        send_list, recv_list = send_list+1, recv_list+1 # Fortran expects that indices start at 1
        self.dynamico_init_transfer(index, send_num, send_size, send_rank, send_len, send_list,
                                       recv_num, recv_size, recv_rank, recv_len, recv_list)
    def dynamico_init_transfer(*args) : pass # overriden in dev.parallel

class LocalArray: # a base class for arrays with halos
    def read_own(self, parray): self.put(self.halo.own_local, self.halo.get_own(parray))
    def read_all(self, parray): self.data = self.halo.get_all(parray)
    def update(self):
        halo=self.halo
        comm, tag = halo.comm, halo.tag
        for rank,loc in halo.send_list:
            comm.send(self.get(loc), dest=rank, tag=tag)
        for rank,loc in halo.recv_list:
            self.put(loc, comm.recv(source=rank, tag=tag) )
        
class LocalArray1(LocalArray): # a 1D array with halo
    def __init__(self, halo):
        self.halo, self.data = halo, np.zeros((len(halo.cells),))
    def get(self, cells): return np.asarray([self.data[i] for i in cells])
    def put(self,cells,data): self.data[cells]=data
class LocalArray2(LocalArray): # a 2D array with halo
    def __init__(self, halo, llm):
        self.halo, self.llm, self.data = halo, llm, np.array((halo.local_len,llm))
    def get(self, cells): return self.data[cells,:]
    def put(self,cells,data): self.data[cells,:]=data