+# @author: Ismael Boureima
+import argparse
+import os, sys,time
+from .toolz import log, amber, blue, green, purple, red
+import cupy as cp
+try:
+ from cupy.cuda import nccl
+ NCCL_backend = 1
+except:
+ NCCL_backend = 0
+
+try:
+ import asyncio
+ import ucp
+ UCX_backend = 1
+except:
+ UCXL_backend = 0
+
+import numpy as np
+from pyDNMFk.communicators import NCCLComm
+from mpi4py import MPI
+from pyDNMFk.get_reduce_comm_tree import getReduceCommTree
+
+
+[docs]class GetTopology():
+ r"""
+ Constructor for the GetTopology class.
+
+ This class is used for identifying the system topology in the context of MPI.
+ It determines the global and local ranks, identifies the host system, and
+ determines GPU availability and assignments.
+
+ Args:
+ - comm (MPI Comm): The MPI communicator. Defaults to MPI.COMM_WORLD.
+ - gmasterRank (int): Rank of the global master. Defaults to 0.
+ """
+ def __init__(self, comm=None, gmasterRank=0):
+ if comm is None: comm = MPI.COMM_WORLD
+ rank = comm.Get_rank()
+ size = comm.Get_size()
+ myHost = os.uname()[1]
+ GMASTER = (rank==gmasterRank)
+ nGPUs_local = cp.cuda.runtime.getDeviceCount()
+ self.comm, self.rank, self.size, self.myHost, self.GMASTER, self.nGPU_local = comm, rank, size, myHost, GMASTER, nGPUs_local
+
+ if GMASTER:
+ hosts, local_size, members = [myHost], {myHost:1}, {myHost:[gmasterRank]}
+ for r in range(size):
+ if r != gmasterRank:
+ rh = comm.recv(source=r, tag=10*r)
+ if rh not in hosts:
+ hosts.append(rh)
+ local_size[rh]=0
+ members[rh] = []
+ local_size[rh] += 1
+ if r not in members[rh]: members[rh].append(r)
+ else:
+ hosts, local_size = None, None
+ comm.send(myHost, dest=gmasterRank, tag=10*rank)
+ comm.barrier()
+
+ if GMASTER:
+ #if rank == 0: print("[+] OK020")
+ global_uID = None
+ if NCCL_backend:
+ #if rank == 0: print("[+] OK021")
+ try:
+ global_uID = nccl.get_unique_id()
+ #if rank == 0: print("[+] OK022")
+ except:
+ if rank == 0: print("[!] UNABLE to get NCCL GLOBAL_UID")
+ nHosts = len(hosts)
+ stride = int(size/nHosts)
+ MPI_TOPO = {'GMASTER':rank, 'hosts':hosts, 'nHosts':nHosts, 'stride':stride, 'nccl_GUID':global_uID, 'local_size':local_size, 'members':members}
+ MPI_TOPO['LMASTER'] = {}
+ #if rank == 0: print("[+] OK023")
+ for host in hosts:
+ idx = hosts.index(host)
+ MPI_TOPO['LMASTER'][host]=idx
+ else:
+ MPI_TOPO = None
+ comm.barrier()
+
+ MPI_TOPO = comm.bcast(MPI_TOPO, root=0)
+ #if rank == 0: print("[+] OK04")
+ nHosts = MPI_TOPO['nHosts']
+ stride = MPI_TOPO['stride']
+ hosts = MPI_TOPO['hosts']
+ local_size = MPI_TOPO['local_size']
+ idx = hosts.index(myHost)
+ lrank = rank // nHosts
+ myMaster = MPI_TOPO['LMASTER'][myHost]
+ LMASTER = (lrank==0)
+ self.topology, self.nHost, self.hosts, self.host_idx, self.mpi_stride = MPI_TOPO, nHosts, hosts, idx, stride
+ self.local_size, self.lrank, self.myMaster, self.LMASTER = local_size, lrank, myMaster, LMASTER
+
+ if GMASTER:
+ gIDs = {myHost:np.arange(nGPUs_local, dtype=np.int32)}
+ nGPUs = nGPUs_local
+ lastID = gIDs[myHost][-1]
+ for r in list(MPI_TOPO['LMASTER'].keys()):
+ if MPI_TOPO['LMASTER'][r] != rank:
+ ng = comm.recv(source=MPI_TOPO['LMASTER'][r], tag=100*MPI_TOPO['LMASTER'][r])
+ nGPUs += ng
+ gIDs[r] = np.arange(ng, dtype=np.int32)+lastID+1
+ lastID = gIDs[r][-1]
+ nGPUs = min(nGPUs, size)
+ gIDs['nGPUs'] = min(nGPUs, size)
+ elif ( LMASTER and (not GMASTER) ):
+ gIDs = None
+ ng = None
+ comm.send(nGPUs_local, dest=MPI_TOPO['GMASTER'], tag=100*rank)
+ else:
+ ng = None
+ gIDs = None
+ comm.barrier()
+ #if rank == 0: print("[+] OK06")
+ gIDs = comm.bcast(gIDs, root=MPI_TOPO['GMASTER'])
+ #if rank == 0: print("[+] OK07")
+ nGPUs = gIDs['nGPUs']
+ del gIDs['nGPUs']
+ MPI_TOPO['gIDs'] = gIDs
+ MPI_TOPO['nGPUs'] = nGPUs
+ if lrank < nGPUs_local:
+ self.gID = MPI_TOPO['gIDs'][myHost][lrank]
+
+[docs] def showTopology(self):
+ """
+ Method to display the detected topology for debugging purposes.
+ """
+ if self.GMASTER:
+ role = f"<{red('GMASTER')}>"
+ elif self.LMASTER and (not self.GMASTER):
+ role = f"<{blue('LMASTER')}>"
+ else:
+ role = f"<{amber('WORKER')}>"
+ log(msg = f"{role} on host {self.myHost}", rank=self.rank, lrank=self.lrank)
+ self.comm.barrier()
+
+[docs] def getReduceCommTree(self, VRBZ=False):
+ """
+ Retrieve a communication tree for efficient reduce operations.
+
+ Args:
+ - VRBZ (bool): Verbosity level for debug outputs. Defaults to False.
+ """
+ x = np.arange(self.topology['nGPUs'])
+ self.reduceCommTree = getReduceCommTree(x=x, VRBZ=VRBZ)
+
+
+
+
+[docs]class superCommunicator():
+ r"""
+ Constructor for the superCommunicator class.
+
+ This class acts as a wrapper around the MPI and NCCL communication protocols,
+ enabling simplified and efficient communication setup based on the detected
+ system topology.
+
+ Args:
+ - comm (MPI Comm): The MPI communicator.
+ - gmasterRank (int): Rank of the global master. Defaults to 0.
+ """
+ def __init__(self, comm, gmasterRank=0):
+ self.globalComm = comm
+ self.topology = GetTopology(comm=comm, gmasterRank=gmasterRank)
+ self.comm = {'GLOBAL':comm}
+ self.comms = ['GLOBAL']
+ self.rank = self.topology.rank
+ self.GMASTER = self.topology.GMASTER
+ self.LMASTER = self.topology.LMASTER
+ self.myMaster = self.topology.myMaster
+ self.inContext = {'GLOBAL':True}
+ self.contexts = list(self.inContext.keys())
+
+[docs] def createComm(self, root, members,name, backEnd='NCCL'):
+ """
+ Create a new communication context based on the provided parameters.
+
+ Args:
+ - root (int): Root rank for this communicator.
+ - members (list): List of member ranks included in this communicator.
+ - name (str): Name for this communicator.
+ - backEnd (str): Communication backend (e.g., 'NCCL'). Defaults to 'NCCL'.
+ """
+ if name in self.comms:
+ if self.GMASTER:
+ log(msg = f"[{red('!')} ]Communcator {name} existAlready", rank=self.rank, lrank=self.lrank)
+ return
+ self.comm[name] = NCCL_COMM(comm=self.globalComm, root=root, members=members, name=name)
+ if self.rank in members:
+ self.inContext[name]=True
+ else:
+ self.inContext[name]=False
+ self.contexts = list(self.inContext.keys())
+
+
+[docs]class NCCL_COMM():
+ r"""
+ Constructor for the NCCL_COMM class.
+
+ This class encapsulates NCCL-specific communication functions, enabling
+ efficient GPU-to-GPU communications.
+
+ Args:
+ - comm (MPI Comm): The MPI communicator.
+ - root (int): Root rank for this communicator.
+ - members (list): List of member ranks included in this communicator.
+ - name (str): Name for this communicator.
+ """
+ def __init__(self, comm, root, members, name):
+ if comm.rank in members:
+ self.rank = comm.rank
+ self.lrank = members.index(comm.rank)
+ self.size = len(members)
+ self.ROOT = root
+ uID = None
+ self.MASTER = (comm.rank==root)
+ if self.MASTER:
+ uID = nccl.get_unique_id()
+ print(f"[+] <LMASTER>[{comm.rank}]: uID = {uID[:8]}")
+ for m in members:
+ if m != root: comm.send(uID, dest=m, tag=1000*m)
+ else:
+ uID = comm.recv(source=root, tag=1000*self.rank)
+ print(f"[+] <WORKER>[{comm.rank}]: uID = {uID[:8]}")
+ #comm.barrier()
+ if uID is not None:
+ print(f" -> Building NCCL comm [{len(members)}, uID:{uID[:8]}, lrank:{self.lrank}]")
+ comm.barrier()
+ ncclcomm = nccl.NcclCommunicator.initAll(members)
+ if self.MASTER: log(msg = f"[{green('+')}] NCCL comm {name} built OK", rank=self.rank, lrank=self.lrank)
+ else:
+ ncclcomm = None
+ self.comm = ncclcomm
+
+# @author: Ismael Boureima
+import sys
+from mpi4py import MPI
+
+# Attempt to import the cupy library for GPU array operations.
+try:
+ import cupy as cp
+except ImportError:
+ print("Unable to import Cupy.", file=sys.stderr)
+ cp = None
+
+# Attempt to import NCCL for GPU-to-GPU communication support from cupy.
+try:
+ from cupy.cuda import nccl
+except ImportError:
+ import sys
+ print("NCCL is not supported.", file=sys.stderr)
+ nccl = None
+
+
+[docs]def get_NCCL_unique_id():
+ """
+ Fetch a unique ID for NCCL group communication.
+
+ Returns:
+ - Unique ID if NCCL is available, otherwise None.
+ """
+ if nccl is not None: return nccl.get_unique_id()
+
+
+
+[docs]class MPIComm():
+ """
+ Class to handle basic MPI communication.
+
+ This class acts as a simple wrapper around basic MPI communication
+ functions for simplification and clarity.
+ """
+ def __init__(self):
+ self.comm = MPI.COMM_WORLD
+ self.size = self.comm.Get_size()
+ self.rank = self.comm.Get_rank()
+
+[docs] def Allreduce(self, send_arr, op=MPI.SUM, stream=None):
+ """
+ Perform an all-reduce operation across all MPI processes.
+
+ Args:
+ send_arr (array-like): Data to be reduced.
+ op (MPI.Op): The reduction operation to be applied. Default is MPI.SUM.
+ stream: CUDA stream for async operations (currently unused).
+
+ Returns:
+ array-like: The reduced array.
+ """
+
+ return self.comm.allreduce(send_arr)
+
+[docs] def Reduce(self, send_arr, op=MPI.SUM, root=0, stream=None):
+ """
+ Perform a reduce operation over MPI processes.
+
+ Args:
+ send_arr (array-like): Data to be reduced.
+ op (MPI.Op): The reduction operation to be applied. Default is MPI.SUM.
+ root (int): Rank of the root process. Default is 0.
+ stream: CUDA stream for async operations (currently unused).
+ """
+ self.comm.reduce(send_arr, op, root=root)
+
+[docs] def Bcast(self, arr, root=0, stream=None):
+ """
+ Broadcast data to all MPI processes.
+
+ Args:
+ arr (array-like): Data to be broadcasted.
+ root (int): Rank of the root process initiating the broadcast. Default is 0.
+ stream: CUDA stream for async operations (currently unused).
+ """
+ self.comm.bcast(arr, root=root)
+
+
+
+
+[docs]class NCCLComm(nccl.NcclCommunicator):
+ """
+ Class to handle NCCL GPU communication.
+
+ This class acts as a wrapper around NCCL's communication functions
+ optimized for GPU-to-GPU communications.
+ """
+ def __init__(self, ndev, commId, rank):
+ """
+ Args:
+ ndev (int): Number of devices (GPUs) involved in the communication.
+ commId (ncclUniqueId): A unique NCCL ID for group communication.
+ rank (int): Rank of the current process within the group.
+ """
+ super(NCCLComm, self).__init__(ndev, commId, rank)
+ self.comm = nccl.NcclCommunicator
+ self.size = ndev
+ self.rank = rank
+ self.commId = commId
+
+
+[docs] def get_unique_id(self):
+ """
+ Returns:
+ ncclUniqueId: The unique NCCL ID for this communicator.
+ """
+ return self.commId
+
+[docs] def get_NCCL_count_dtype(self, arr):
+ """
+ Determine the data count and data type for NCCL communication based on array dtype.
+ """
+ if arr.dtype == cp.complex64:
+ return arr.size*2, nccl.NCCL_FLOAT32
+ elif arr.dtype == cp.complex128:
+ return arr.size*2, nccl.NCCL_FLOAT64
+ elif arr.dtype == cp.float32:
+ return arr.size, nccl.NCCL_FLOAT32
+ elif arr.dtype == cp.float64:
+ return arr.size, nccl.NCCL_FLOAT64
+ else:
+ raise ValueError("This dtype is not supported by NCCL.")
+
+[docs] def Allreduce(self, send_arr, recv_arr, op=nccl.NCCL_SUM, stream=None):
+ """
+ Perform an all-reduce operation across all GPU processes using NCCL.
+ """
+ sendbuf = send_arr.__cuda_array_interface__['data'][0]
+ recvbuf = recv_arr.__cuda_array_interface__['data'][0]
+ count, datatype = self.get_NCCL_count_dtype(send_arr)
+ if stream is None:
+ stream = cp.cuda.Stream.null.ptr
+ else:
+ stream = stream.ptr
+ super(NCCLComm, self).allReduce(sendbuf, recvbuf, count, datatype, op, stream)
+
+[docs] def Reduce(self, send_arr, recv_arr, op=nccl.NCCL_SUM, root=0, stream=None):
+ """
+ Perform a reduce operation across all GPU processes using NCCL.
+ """
+ sendbuf = send_arr.__cuda_array_interface__['data'][0]
+ recvbuf = recv_arr.__cuda_array_interface__['data'][0]
+ count, datatype = self.get_NCCL_count_dtype(send_arr)
+ if stream is None:
+ stream = cp.cuda.Stream.null.ptr
+ else:
+ stream = stream.ptr
+ super(NCCLComm, self).reduce(sendbuf, recvbuf, count, datatype, op, root, stream)
+
+[docs] def Bcast(self, arr, root=0, stream=None):
+ """
+ Perform broadcast operation from root to all GPU processes using NCCL.
+ """
+ buff = arr.__cuda_array_interface__['data'][0]
+ count, datatype = self.get_NCCL_count_dtype(arr)
+ if stream is None:
+ stream = cp.cuda.Stream.null.ptr
+ else:
+ stream = stream.ptr
+ super(NCCLComm, self).bcast(buff, count, datatype, root, stream)
+
+
+
+
+
+[docs]class SUPER_COMM():
+ """
+ A unified communicator that supports both MPI (for CPUs) and NCCL (for GPUs).
+ """
+ def __init__(self, ndev=None, commId=None, rank=None, backend="MPI"):
+ """
+ Args:
+ ndev (int, optional): Number of devices (GPUs) involved in the communication.
+ commId (ncclUniqueId, optional): A unique NCCL ID for group communication.
+ rank (int, optional): Rank of the current process within the group.
+ backend (str): Communication backend to use. Can be "MPI" or "NCCL".
+ """
+ self.backend = backend.upper()
+ if self.backend in ["MPI","CPU","DEFAULT"]:
+ self.comm = MPI_COMM()
+ elif self.backend in ["NCCL","GPU","CUDA","BEST","PERFORMANCE"]:
+ assert ndev is not None, "[!][COMM ERROR] SUPER_COMM NCCL backend requires valid ndev"
+ assert commId is not None, "[!][COMM ERROR] SUPER_COMM NCCL backend requires valid commID"
+ assert rank is not None, "[!][COMM ERROR] SUPER_COMM NCCL backend requires valid rank"
+ self.comm = NCCL_COMM(ndev=ndev, commId=commId, rank=rank)
+ #self.comm = NCCLComm(ndev=ndev, commId=commId, rank=rank)
+ else:
+ raise "[!][COMM ERROR] Backend {} not supported".format(self.backend)
+ self.size = self.comm.size
+ self.rank = self.comm.rank
+
[docs]def init(arg):
- """ Global variables declaration here. The variables declared within this function in this file are shared across other files and functions during import."""
+ """ Global variables declaration here. The variables declared within this function in this file are shared across other files and functions during import."""
global time, flag
time = {}
flag = 0
+# @author: Ismael Boureima
+# Necessary imports for the module
+import os, time
+import numpy as np
+import numpy
+from queue import Queue
+from .toolz import log, blue, green, red, amber
+
+# Try importing cupy modules for GPU support
+try:
+ import cupy as cp
+ import cupy, cupyx
+except:
+ raise Exception(f"[{red('!!')}] Unable to import Cupy")
+
+# Try importing necessary cupy backend modules
+try:
+ from cupy.cuda import nccl
+except:
+ raise Exception(f"[{red('!!')}] Unable to import NCCL Cupy backend")
+
+# Importing specific functions from cupy to work with
+from cupy import asarray as _asarray
+from cupy import asnumpy as _asnumpy
+from cupy import divide as _divide
+from cupy import matmul as _matmul
+from cupy import multiply as _multiply
+from cupy import zeros as _zeros
+
+# Additional imports for the module
+from .communicators import NCCLComm
+from .comm_utils import GetTopology
+from mpi4py import MPI
+from .cupyCuSPARSELib import spMM as _spMM
+from .cupyCuSPARSELib import spMat, spRandMat, spMM
+from .cupy_utils import pin_memory as _pin_mem
+from scipy.sparse import coo_matrix, csc_matrix, csr_matrix, dia_matrix, issparse
+import scipy
+from .utils import data_operations
+
+
+# Defining some constants for the module
+PRECISIONS = {'float32':np.float32, 'float64':np.float64, 'int16':np.int16, 'int32':np.int32,'int64':np.int64}
+WORKLOADS = {'nmf':'NMF', 'nmfk':'NMFk', 'bench':'BENCHMARK' }
+NMF_NORMS = {'fro':'FROBENUS', 'kl':'KL'}
+NMF_METHODS = {'mu':'Multiplicative Update'}
+NMF_INITS = {'rand':'RANDOM', 'random':'RANDOM', 'nnsvd':'Non-Negative SVD', 'svd':'SVD'}
+
+# Utility function to raise an exception for unsupported operations
+[docs]def noSupportFor(X):
+ raise Exception(f"[{red('!!')}] {amber(X)} is not yet supported on GPU")
+
+
+# The main class to perform NMF on GPU
+[docs]class cudaNMF():
+
+ def __init__(self, A_ij,k, params, factors=None, save_factors=False):
+ r"""
+ Initialize the cudaNMF class.
+
+ Parameters
+ ----------
+ A_ij : array-like
+ The matrix to be factorized.
+ k : int
+ Rank of the factorization.
+ params : dict
+ Configuration parameters for NMF.
+ factors : array-like, optional
+ Initial guess for the factors.
+ save_factors : bool, optional
+ Whether to save factors or not.
+
+ """
+ self.A_ij = A_ij
+ self.params = params
+ self.k = k
+
+ self.factors = factors
+ self.save_factors = save_factors
+ self.VRBZ = params.verbose
+ try:
+ self.comm = params.comm1
+ except:
+ print(f"{149*'#'}")
+ print("[{}] --[{} comm {}]--".format(red("!"),blue('MPI'), red("NOT FOUND")))
+ self.comm = MPI.COMM_WORLD
+ self.rank = self.comm.rank
+ if self.rank == 0:
+ print(f"{149*'#'}")
+ log("[{}] --[{} comm built {}]--".format(green("+"),blue('MPI'), green("OK")))
+ self.rank = self.comm.rank
+ self.eps = np.finfo(self.A_ij.dtype).eps
+ self.params.eps = self.eps
+ self.precision = self.params.precision
+ self.topo_dim = self.params.topo_dim
+ self.local_m,self.local_n = self.A_ij.shape
+ self.p_r, self.p_c = self.params.p_r, self.params.p_c
+ #[1] Identify distributed system's topology
+ self.NCCL_COMM_BUILT = False
+ self.identify_distributed_system_compute_topology()
+ #[2] Identify Global compute grid parameters
+ self.data_op = data_operations(self.A_ij, self.params)
+ self.params = self.data_op.params
+ self.GPU_GRID_BUILT = False
+ #[3] Identify A matrix partition method
+ self.find_optimat_partition_strategy()
+ #[4] Build GPU communication COMM
+ self.checkGPUContext()
+ self.buildSubCommunicators()
+ METHOD = f"NMF_{self.GPU_PARTITION}_cuda"
+ #[5] Check format (dense/sparse) of A matrix
+ self.SPARSE_BUFF_PINNED_MEM_PARAMS_CONFIGURED = False
+ self.check_A_matrix_format()
+ if self.A_IS_SPARSE:
+ METHOD = f"sparse_{METHOD}"
+ else:
+ METHOD = f"dense_{METHOD}"
+ self.NMFMethod = METHOD
+ #[6] Identifiy arrays that will be cached on GPU
+ self.A_IS_BATCHED = self.params.A_is_batched
+ self.A_is_cached,self.W_is_cached,self.H_is_cached = None, None, None
+ try:
+ gpu_arrays = params.GPU_arrays
+ except:
+ if self.rank == 0:
+ self.log(msg = f"[!][CACHING] GPU arrays unspecified")
+ self.log(msg = f"[!][CACHING] Assuming A, W, and H are cached on GPU")
+ gpu_arrays = []
+ self.gpu_arrays = []
+ self.automatic_gpu_array_selection = False
+ for ARR in gpu_arrays:
+ arr = ARR.lower()
+ if arr in ['auto', 'optim']:
+ self.automatic_gpu_array_selection = True
+ else:
+ self.gpu_arrays.append(arr)
+ if self.automatic_gpu_array_selection:
+ if self.rank == 0: self.log(msg = f"[{red('!')}][{amber('COMING SOON')}] Automatic/Optimal GPU array selection is {red('Not yet supported')}")
+ #[7] Build GPU MEM grid
+ self.nBatch = self.params.nBatch
+ self.batchQeueSize = self.params.batchQeueSize
+ self.qs = self.batchQeueSize
+ self.SQ = None
+ self.MAXSTREAMS = self.params.MAXSTREAMS
+ if not self.GPU_GRID_BUILT: self.get_gpu_grid_params()
+ self.BATCHING_AXIS = None
+ self.COL_BATCHING_AXIS = None
+ self.ROW_BATCHING_AXIS = None
+ self.find_optimal_baching_axis()
+ self.set_batching_params()
+ if self.A_IS_BATCHED:
+ self.NMFMethod = "{}_BATCHED".format(self.NMFMethod)
+ else:
+ self.cache_A_on_device()
+ self.X_per_initialized = False
+ self.show_A_info()
+ #[8] Prepare events
+ self.GB = 1024**3
+ self.TB = 1024*self.GB
+ self.dt = {'NMF':0.0, 'W_up':0.0, 'H_up':0.0,
+ "allRed_XHT":0.0, "allRed_WHHT":0.0,
+ 'allRed_WTX':0.0, 'allRed_WTWH':0.0,
+ 'allRed_WTW':0.0, 'allRed_HHT':0.0,
+ 'H2D_A':0.0, 'H2D_H':0.0, 'H2D_W':0.0,
+ 'D2H_A':0.0, 'D2H_H':0.0, 'D2H_W':0.0,
+ 'units':'[ms]'}
+ self.events = {}
+ self.events['start'] = cp.cuda.Event()
+ self.events['end'] = cp.cuda.Event()
+ self.events['h2d_s'] = cp.cuda.Event()
+ self.events['h2d_e'] = cp.cuda.Event()
+ self.events['d2h_s'] = cp.cuda.Event()
+ self.events['d2h_e'] = cp.cuda.Event()
+ self.events['reduce_s'] = cp.cuda.Event()
+ self.events['reduce_e'] = cp.cuda.Event()
+ self.events['nmf_start'] = cp.cuda.Event()
+ self.events['nmf_end'] = cp.cuda.Event()
+ #[9] Build random generation initial state seeds # This is used so that cofactors are initialized
+ # locally and therefor avoiding broadcat communications
+ if self.rank==0:
+ iter_seeds = np.random.randint(100*self.params.end_k, size=2*self.params.end_k)
+ else:
+ iter_seeds = None
+ iter_seeds=self.comm.bcast(iter_seeds, root=0)
+ self.iter_seeds = iter_seeds
+ #[10] Configure Pinned Memory management Pool
+ #self.PINNEDMEMPOOL = cp.cuda.PinnedMemoryPool()
+ #cp.cuda.set_pinned_memory_allocator(self.PINNEDMEMPOOL.malloc)
+ if self.A_IS_BATCHED:
+ self.MEMPOOL = cp.get_default_memory_pool()
+ self.PINNEDMEMPOOL = cp.cuda.PinnedMemoryPool()
+ cp.cuda.set_pinned_memory_allocator(self.PINNEDMEMPOOL.malloc)
+ else:
+ self.MEMPOOL = cp.get_default_memory_pool()
+ self.PINNEDMEMPOOL = cp.get_default_pinned_memory_pool()
+ self.showMemStats()
+ #[11] Build Qeues of CUDA STREAMS for Async copies management
+ self.STREAMS = {}
+ self.stream = []
+ self.FREE_STREAM = {}
+ self.REDUCE_STREAM = cp.cuda.stream.Stream()
+ self.MAXSTREAMS = min(self.nBatch, self.MAXSTREAMS)
+ for n in range(self.MAXSTREAMS):
+ self.stream.append("s%04d" %n)
+ self.STREAMS[self.stream[-1]] = cp.cuda.stream.Stream()
+ self.FREE_STREAM[self.stream[-1]] = cp.cuda.Event()
+
+
+[docs] def log(self,msg):
+ r"""
+ Prints a log message.
+
+ Parameters
+ ----------
+ msg : str
+ The message to be logged.
+
+ """
+ log(msg = msg, rank=self.rank, lrank=self.lrank)
+
+
+[docs] def isCupyObject(self,x):
+ r"""
+ Checks if the input object is a Cupy object.
+
+ Parameters
+ ----------
+ x : object
+ The object to be checked.
+
+ Returns
+ -------
+ bool
+ True if x is a Cupy object, False otherwise.
+
+ """
+ return type(x) in [cp._core.core.ndarray]
+
+
+[docs] def getArrayType(self,x):
+ r"""
+ Returns the type of the array (CPU/GPU, dense/sparse, etc.).
+
+ Parameters
+ ----------
+ x : array-like
+ The array whose type needs to be determined.
+
+ Returns
+ -------
+ str
+ Type of the array (e.g., "CPU_DENSE", "GPU_SPARSE").
+
+ """
+ T = type(x)
+ if T in [list]:
+ return "LIST[{}]".format(self.getObjectType(x[0]))
+ elif T in [np.ndarray]:
+ return "CPU_DENSE"
+ elif T in [cp._core.core.ndarray]:
+ return "GPU_DENSE"
+ elif T in [scipy.sparse.coo.coo_matrix, scipy.sparse.csr.csr_matrix, scipy.sparse.csc.csc_matrix]:
+ return "CPU_SPARSE"
+ elif T in [cupyx.scipy.sparse.coo.coo_matrix, cupyx.scipy.sparse.csr.csr_matrix, cupyx.scipy.sparse.csc.csc_matrix]:
+ return "GPU_SPARSE"
+ else:
+ return "UNKNOWN"
+
+
+[docs] def identify_distributed_system_compute_topology(self):
+ r"""
+ Identifies the distributed system's compute topology and sets up the necessary
+ attributes related to topology.
+ """
+ self.topology = GetTopology(self.comm)
+ if self.VRBZ : self.show_topology()
+ if self.rank == 0:
+ global_uID = nccl.get_unique_id()
+ else:
+ global_uID = None
+ self.global_uID = self.comm.bcast(global_uID, root=0)
+ self.comm.barrier()
+ self.myHost = self.topology.myHost
+ self.lrank = self.topology.lrank
+ self.nGPU_local = self.topology.nGPU_local
+ self.nGPUs = self.topology.topology['nGPUs']
+
+
+[docs] def compute_global_dim(self): ##[IDB]: This is needed here
+ r"""Computes global dimensions m and n from given chunk sizes for any grid configuration"""
+ if self.p_r != 1 and self.p_c == 1:
+ self.global_n = self.local_n
+ self.global_m = self.comm.allreduce(self.local_m)
+ elif self.p_c != 1 and self.p_r == 1:
+ self.global_n = self.comm.allreduce(self.local_n)
+ self.global_m = self.local_m
+ else:
+ if self.rank % self.p_c == 0:
+ self.global_m = self.local_m
+ else:
+ self.global_m = 0
+ self.global_m = self.comm.allreduce(self.global_m)
+ if self.rank // self.p_c == 0:
+ self.global_n = self.local_n
+ else:
+ self.global_n = 0
+ self.global_n = self.comm.allreduce(self.global_n)
+ self.comm.barrier()
+
+
+[docs] def checkGPUContext(self):
+ r"""
+ Determines if the current instance is operating within a GPU context and sets up
+ GPU-related attributes.
+ """
+ if self.lrank < self.nGPU_local:
+ self.IN_GPU_CONTEXT = True
+ self.lID = self.lrank # Local GPU ID
+ self.gID = self.topology.gID # Global GPU ID
+ cp.cuda.Device(self.lrank).use()
+ #self.device = cp.cuda.Device()
+ self.cudaDevice = cp.cuda.runtime.getDevice()
+ if self.VRBZ: self.log(msg=f"I am using <Global Device{self.gID}> | <Local Device{self.cudaDevice}> on {self.myHost}")
+ else:
+ self.IN_GPU_CONTEXT = False
+ self.lID = None
+ self.gID = None
+ if self.VRBZ: self.log(msg = f"I am not using CUDA")
+ self.comm.barrier()
+ self.nccl_group_size = self.comm.allreduce(int(self.IN_GPU_CONTEXT), op=MPI.SUM)
+ if self.rank == 0: self.log(msg = f"[+][NCCL] GROUP size = {self.nccl_group_size}")
+
+
+[docs] def check_A_matrix_format(self):
+ r"""
+ Checks the format of matrix A, decides if it's sparse or dense, and sets the
+ necessary attributes related to the matrix's format.
+ """
+ if type(self.A_ij) in [scipy.sparse.coo.coo_matrix, scipy.sparse.csr.csr_matrix, scipy.sparse.csc.csc_matrix]:
+ if self.params.densify_A: # Dense A and Array format
+ self.A_ij = self.A_ij.toarray() # Converting to dense for now (Fix Not Efficient)
+ self.A_in_sparse_format = False
+ self.A_IS_SPARSE = False
+ else: # Sparse A and sparse format
+ self.A_in_sparse_format = True
+ self.A_IS_SPARSE = True
+ self.A_nnz = self.A_ij.nnz
+ self.A_density = self.A_nnz/(self.global_m*self.global_n)
+ self.MM = _spMM
+ else:
+ self.A_in_sparse_format = False
+ if not self.A_in_sparse_format:
+ self.A_nnz = np.count_nonzero(self.A_ij)
+ self.A_nnz = self.comm.allreduce(self.A_nnz)
+ self.A_density = self.A_nnz/(self.global_m*self.global_n)
+ if self.A_density > self.params.sparsity_tresh: # Dense A and Array format
+ #self.NMFMethod = 'dense_{}'.format(self.NMFMethod)
+ self.A_IS_SPARSE = False
+ self.MM = _matmul
+ else: # Dense A and sparse format
+ self.A_IS_SPARSE = True
+ if self.params.sparsify_A:
+ self.A_ij = csr_matrix(self.A_ij)
+ self.A_in_sparse_format = True
+ self.MM = _spMM
+ else:
+ self.A_IS_SPARSE = False
+ self.MM = _matmul
+ self.A_IS_DENSE = not self.A_IS_SPARSE
+
+
+[docs] def show_A_info(self):
+ r"""
+ Logs and displays information about matrix A, including dimensions, density,
+ and other related attributes.
+ """
+ if self.rank == 0:
+ self.log(f"[{green('+')}] [ {blue('Local')} ] A_ij[{amber(self.local_m)}, {amber(self.local_n)}] ({red(round(np.prod(self.A_ij.shape)*4/(1024**3), 4) )} {red('GB')}) built {green('OK')}")
+ self.log(f"[{green('+')}] [{blue('Distributed')}] A_ij[{amber(self.global_m)}, {amber(self.global_n)}] ({red( round(self.global_m*self.global_n*4/(1024**3), 4 ) )} {red('GB')}) built {green('OK')}")
+ self.log(f"[{green('+')}] Density = {self.A_density}")
+ if self.COL_PARTITION:
+ self.log(f"[{green('+')}] Nr = nBatch = {self.nBatch} -> Batch Size = m/nBatch = {self.batch_size}")
+ else:
+ self.log(f"[{green('+')}] Nc = nBatch = {self.nBatch} -> Batch Size = n/nBatch = {self.batch_size}")
+ if self.last_chunk > 0 : self.log(f"[{amber('!')}] Uneven split -> Last Batch Size = {self.last_chunk}")
+ #self.log(f"[{green('+')}] Nc = nGPUs = {self.grid_Nc} -> J = n/Nc = {self.grid_J}")
+ self.log(f"[{green('+')}] Batch Qeue size = {self.batchQeueSize}")
+ print(f"{149*'#'}")
+ self.log(f"[{green('+')}] {blue('WORKLOAD')} : {green(WORKLOADS[self.params.work_load.lower()])}")
+ self.log(f"[{green('+')}] {blue('NMF METHOD')} : {green(NMF_METHODS[self.params.method.lower()])}")
+ self.log(f"[{green('+')}] {blue('NMF NORM')} : {green(NMF_NORMS[self.params.norm.lower()])}")
+ self.log(f"[{green('+')}] {blue('NMF INIT')} : {green(NMF_INITS[self.params.init.lower()])}")
+ self.log(f"[{green('+')}] {blue('Number of perturbation')} : {green(str(self.params.perturbations))}")
+ self.log(f"[{green('+')}] {blue('Iterations/Perturbation')} : {green( str(self.params.itr)) }")
+ self.log(f"[{green('+')}] {blue('Range k')} : {green(str(self.params.start_k))} < k < {amber(str(self.params.end_k))}")
+ self.log(f"[{green('+')}] {blue('Delta k')} : {green(str(self.params.step))}")
+ self.log(f"[{green('+')}] {blue('Sill Threshold')} : {green(str(self.params.sill_thr)) }")
+ self.log(f"[{green('+')}] {blue('Noise')} : {green(str(self.params.noise_var)) }")
+ self.log(f"[{green('+')}] {blue('GPU Arrays')} : {green(self.gpu_arrays)}")
+ print(f"{149*'#'}")
+ print(f"{149*'#'}")
+ self.sub_comm.barrier()
+
+
+[docs] def show_topology(self):
+ r"""
+ Displays information about the system's topology.
+ """
+ self.topology.showTopology()
+ self.comm.barrier()
+ log(f"[+] topo_dim = {self.params.topo_dim}", rank=self.rank)
+ log(f"[+] myHost = {self.topology.myHost}", rank=self.rank)
+ log(f"[+] lrank = {self.topology.lrank}", rank=self.rank)
+ log(f"[+] nGPU_local = {self.topology.nGPU_local}", rank=self.rank)
+ log(f"[+] nGPUs = {self.topology.topology['nGPUs']}", rank=self.rank)
+ self.comm.barrier()
+
+
+[docs] def buildSubCommunicators(self, rebuildNCCLCOM=False):
+ r"""
+ Args:
+ rebuildNCCLCOM (bool): Flag indicating if NCCL Communicator should be rebuilt.
+
+ Builds sub-communicators for the distributed system. Optionally rebuilds the NCCL communicator.
+ """
+ ...
+ if self.IN_GPU_CONTEXT:
+ color = 55
+ else:
+ color = 77
+ self.color = color
+ self.comm.barrier()
+ self.sub_comm = self.comm.Split(color, self.rank)
+ self.comm.barrier()
+ if ((not self.NCCL_COMM_BUILT) or (rebuildNCCLCOM)): self.buildNCCLSubCommunicator()
+
+
+[docs] def buildNCCLSubCommunicator(self):
+ r"""
+ Builds the NCCL sub-communicator for GPU context. Sets the related attributes.
+ """
+ if self.IN_GPU_CONTEXT:
+ if self.rank == 0:
+ if self.VRBZ: self.log(f"[+] Broadcasting NCCL global_UID")
+ global_uID = nccl.get_unique_id()
+ else:
+ global_uID = None
+ self.global_uID = self.comm.bcast(global_uID, root=0)
+ if self.VRBZ: self.log(f"[+] NCCL global_UID received OK")
+ self.comm.barrier()
+ self.nccl_comm = NCCLComm(ndev=self.nccl_group_size, commId=self.global_uID, rank=self.rank)
+ else:
+ self.nccl_comm = None
+ self.comm.barrier()
+ self.NCCL_COMM_BUILT = True
+ if self.rank == 0:
+ print(f"{149*'#'}")
+ self.log(f"[{green('+')}] --[{blue('NCCL')} comm built {green('OK')}]--")
+ print(f"{149*'#'}")
+ self.comm.barrier()
+
+
+[docs] def get_gpu_grid_params(self):
+ r"""
+ Fetches or calculates GPU grid parameters based on the current setup.
+ Sets the related attributes.
+ """
+ if self.GPU_GRID_BUILT: return
+ self.grid_loc_m, self.grid_loc_n = self.local_m, self.local_n
+ self.grid_glob_m, self.grid_glob_n = self.global_m, self.global_n
+ if self.params.init.lower() in ['nnsvd','svd','singular']:
+ if self.COL_PARTITION: # n>m
+ self.grid_L = int(self.grid_glob_m/self.nGPUs)
+ self.grid_q = int(self.grid_L/self.batch_size)
+ else:
+ self.grid_L = int(self.grid_glob_n/self.nGPUs)
+ self.grid_q = int(self.grid_L/self.batch_size)
+
+
+
+[docs] def cache_A_on_device(self):
+ r"""
+ Caches matrix A on the device for faster access during GPU operations.
+ """
+ if not self.GPU_GRID_BUILT: self.get_gpu_grid_params()
+ if self.IN_GPU_CONTEXT:
+ self.A_d = []
+ for b in range(self.nBatch):
+ if self.COL_PARTITION:
+ i0, i1 = b*self.batch_size, (b+1)*self.batch_size
+ if self.A_IS_SPARSE:
+ self.A_d.append(spMat( self.A_ij[i0:i1, :], Format='csr'))
+ else:
+ self.A_d.append(_asarray(self.A_ij[i0:i1, :], dtype=self.precision))
+ else:
+ j0, j1 = b*self.batch_size, (b+1)*self.batch_size
+ if self.A_IS_SPARSE:
+ self.A_d.append(spMat( self.A_ij[:, j0:j1], Format='csr'))
+ else:
+ self.A_d.append(_asarray(self.A_ij[:, j0:j1], dtype=self.precision))
+ else:
+ self.A_d = None
+ self.comm.barrier()
+ self.A_is_cached = True
+
+
+[docs] def cache_H_on_device(self): # W.I.P
+ r"""
+ This method sets the flag indicating that the matrix H is cached on the device.
+ """
+ self.H_is_cached = True
+
+
+[docs] def cache_W_on_device(self): # W.I,P
+ r"""
+ This method sets the flag indicating that the matrix W is cached on the device.
+ """
+ self.W_is_cached = True
+
+
+
+[docs] def find_optimat_partition_strategy(self):
+ r"""
+ This method determines the best GPU partitioning strategy based on the provided parameters.
+ It supports partitioning in 'row', 'column' or 'auto' modes.
+ """
+ # Attempt to set partition type from parameters; default to 'auto' on failure.
+ try:
+ partition_type = self.params.gpu_partition.lower()
+ except:
+ partition_type = 'auto'
+
+ # Check the partition type and set the GPU_PARTITION attribute accordingly.
+ # Multiple alias names are supported for each partitioning strategy.
+ if partition_type in ['row', 'row1', '1row' ,'row_1d', '1d_row']:
+ self.GPU_PARTITION = 'row_1d'
+ # ... (similar checks for other partition types) ...
+ elif partition_type in ['col', 'col1', '1col' ,'col_1d', '1d_col', 'column', 'column_1d', '1d_column']:
+ self.GPU_PARTITION = 'col_1d'
+ # If 'auto' mode, optimize for grid partition.
+ elif partition_type in ['auto', 'optim', 'optimal']:
+ # ... (logic to automatically choose the partition type) ...
+ self.GPU_PARTITION = 'auto'
+ if self.rank == 0:
+ log(msg = f"[{amber('!')}][PARTION]: GPU_PARTITION is set to {amber(self.GPU_PARTITION)}", rank=self.rank, lrank=self.lrank)
+ log(msg = f"[{amber('.')}][PARTION]: Optimizing for grid partion ...", rank=self.rank, lrank=self.lrank)
+ if self.global_m >= self.global_n:
+ self.GPU_PARTITION = 'row_1d'
+ else:
+ self.GPU_PARTITION = 'col_1d'
+ if self.rank == 0: log(msg = f"[{green('+')}][PARTION]: GPU_PARTITION is set to {green(self.GPU_PARTITION)}", rank=self.rank, lrank=self.lrank)
+ else:
+ raise Exception(f"[!!] GPU grid partition {red(params.gpu_partition)} is not supported. Try {green('auto')}, {green('column')}, {green('row')}")
+ # Set ROW and COL partition flags based on GPU_PARTITION.
+ if self.GPU_PARTITION in ['row_1d']:
+ self.COL_PARTITION = True
+ else:
+ self.COL_PARTITION = False
+ self.ROW_PARTITION = not self.COL_PARTITION
+
+
+
+[docs] def find_optimal_baching_axis(self):
+ try:
+ self.BATCHING_AXIS = self.params.batching_axis.lower()
+ except:
+ self.BATCHING_AXIS = 'auto'
+ if self.BATCHING_AXIS in ['row', 'row1', '1row' ,'row_1d', '1d_row']:
+ self.BATCHING_AXIS = 'row_1d'
+ elif self.BATCHING_AXIS in ['col', 'col1', '1col' ,'col_1d', '1d_col', 'column', 'column_1d', '1d_column']:
+ self.BATCHING_AXIS = 'col_1d'
+ elif self.BATCHING_AXIS in ['auto', 'optim', 'optimal']:
+ self.BATCHING_AXIS = 'auto'
+ if self.rank == 0:
+ log(msg = f"[{amber('!')}][BATCHING]: STRATEGY is set to {amber(self.GPU_PARTITION)}", rank=self.rank, lrank=self.lrank)
+ log(msg = f"[{amber('.')}][BATCHING]: Optimizing for BATCHING STRATEGY...", rank=self.rank, lrank=self.lrank)
+ if self.grid_loc_m >= self.grid_loc_n:
+ self.BATCHING_AXIS = 'row_1d'
+ else:
+ self.BATCHING_AXIS = 'col_1d'
+ if self.rank == 0: log(msg = f"[{green('+')}][BATCHING]: STRATEGY is set to {green(self.BATCHING_AXIS)}", rank=self.rank, lrank=self.lrank)
+ else:
+ raise Exception(f"[!!] BATCHING STRATEGY {red(params.BATCHING_AXIS)} is not supported. Try {green('auto')}, {green('column')}, {green('row')}")
+ if self.BATCHING_AXIS in ['row_1d']:
+ self.ROW_BATCHING_AXIS = True
+ else:
+ self.ROW_BATCHING_AXIS = False
+ self.COL_BATCHING_AXIS = not self.ROW_BATCHING_AXIS
+
+
+[docs] def set_batching_params(self):
+ if self.ROW_BATCHING_AXIS: # Row-wise batchingn
+ self.batch_size = self.local_m//self.nBatch
+ self.last_chunk = self.local_m%self.nBatch
+ self.grid_I = self.batch_size # min(self.params.IMAX, self.grid_glob_m)
+ self.params.grid_I = self.grid_I
+ self.grid_Nr, self.grid_Nc = self.nBatch, 1 # int(self.grid_loc_m/self.grid_I), int(1*self.nGPUs)
+ self.grid_J = self.local_n # int(self.grid_loc_n/self.grid_Nc)
+ else: # Column-wise batching
+ self.batch_size = self.local_n//self.nBatch
+ self.last_chunk = self.local_n%self.nBatch
+ self.grid_J = self.batch_size # min(self.params.JMAX, self.grid_loc_n)
+ self.params.grid_J = self.grid_J
+ self.grid_Nr, self.grid_Nc = 1,self.nBatch # 1, int(self.grid_loc_n/self.grid_J)
+ self.grid_I = self.local_m # int(self.grid_loc_m/self.grid_Nr)
+
+
+[docs] def get_managed_stream_queue(self):
+ """
+ Initializes and returns a dictionary for managing stream queues for GPU operations.
+ """
+ SQ={}
+ SQ['Queue'] = Queue(maxsize = self.qs)
+ SQ['READY'] = {}
+ #[1] En-queue working streams and set their states to ready to work
+ for j in range(self.qs):
+ i = j+1
+ SQ['READY'][i] = cp.cuda.stream.Stream()
+ SQ['Queue'].put(i)
+ SQ['REDUCE'] = cp.cuda.stream.Stream()
+ SQ['COMPLETED'] = []
+ return SQ
+
+
+[docs] def configure_sparse_buffers_pinned_mem_params(self):
+ r"""
+ Configures the parameters for pinned memory buffers when dealing with sparse matrices.
+ """
+ if self.SPARSE_BUFF_PINNED_MEM_PARAMS_CONFIGURED: return
+ assert self.BATCHING_AXIS is not None, "[!] PINNED MEM BUFFERS are used only for BATCHING"
+ assert self.A_IS_SPARSE, "[!] SPARSE PINNED MEM BUFFERS are used only for BATCHING A in SPARSE FORMAT"
+ self.buff_idx, self.buff_ptr = [], []
+ self.sparse_vect_size = {'dat':[],'idx':[],'ptr':[]}
+ self.sparseData_max_vect_size, self.sparseIdx_max_vect_size, self.sparsePtr_max_vect_size = 0,0,0
+ for b in range(self.nBatch):
+ #if self.ROW_PARTITION:
+ if self.ROW_BATCHING_AXIS:
+ i0, i1 = b*self.batch_size, (b+1)*self.batch_size
+ self.sparse_vect_size['dat'].append(len(self.A_ij[i0:i1, :].data))
+ self.sparse_vect_size['idx'].append(len(self.A_ij[i0:i1, :].indices))
+ self.sparse_vect_size['ptr'].append(len(self.A_ij[i0:i1, :].indptr))
+ self.sparseData_max_vect_size = max(self.sparseData_max_vect_size, len(self.A_ij[i0:i1, :].data))
+ self.sparseIdx_max_vect_size = max(self.sparseIdx_max_vect_size, len(self.A_ij[i0:i1, :].indices))
+ self.sparsePtr_max_vect_size = max(self.sparsePtr_max_vect_size, len(self.A_ij[i0:i1, :].indptr))
+ else: # Col Partition
+ j0, j1 = b*self.batch_size, (b+1)*self.batch_size
+ self.sparse_vect_size['dat'].append(len(self.A_ij[:, j0:j1].data))
+ self.sparse_vect_size['idx'].append(len(self.A_ij[:, j0:j1].indices))
+ self.sparse_vect_size['ptr'].append(len(self.A_ij[:, j0:j1].indptr))
+ self.sparseData_max_vect_size = max(self.sparseData_max_vect_size, len(self.A_ij[:, j0:j1].data))
+ self.sparseIdx_max_vect_size = max(self.sparseIdx_max_vect_size, len(self.A_ij[:, j0:j1].indices))
+ self.sparsePtr_max_vect_size = max(self.sparsePtr_max_vect_size, len(self.A_ij[:, j0:j1].indptr))
+ #print("[b%04d] lenData, lenIdx, lenPtr = {},{},{}".format(self.sparseData_max_vect_size, self.sparseIdx_max_vect_size, self.sparsePtr_max_vect_size) %b)
+ self.SPARSE_BUFF_PINNED_MEM_PARAMS_CONFIGURED = True
+
+
+
+
+[docs] def allocate_sparse_batch_buffers(self):
+ r"""
+ Allocate memory for sparse batch buffers.
+
+ This method configures memory parameters if not already set. Based on the specified
+ batching axis, it initializes the appropriate buffers to store the data.
+ """
+ # Check if memory parameters for sparse buffers are configured
+ if not self.SPARSE_BUFF_PINNED_MEM_PARAMS_CONFIGURED:
+ self.configure_sparse_buffers_pinned_mem_params()
+
+ # Allocate memory based on column batching axis
+ if self.COL_BATCHING_AXIS:
+ self.H_d, self.X_d = [], []
+ for n in range(self.batchQeueSize):
+ self.H_d.append(cp.empty((self.k, self.batch_size), dtype=self.A_ij.dtype))
+ self.X_d.append(cp.empty(self.sparseData_max_vect_size, dtype=self.A_ij.dtype))
+ self.buff_idx.append(cp.empty(self.sparseIdx_max_vect_size, dtype=cp.int32))
+ self.buff_ptr.append(cp.empty(self.sparsePtr_max_vect_size, dtype=cp.int32))
+ else: # ROW_BATCHING_AXIS
+ self.W_d, self.X_d = [], []
+ for n in range(self.batchQeueSize):
+ self.W_d.append(cp.empty((self.batch_size, self.k), dtype=self.A_ij.dtype))
+ self.X_d.append(cp.empty(self.sparseData_max_vect_size, dtype=self.A_ij.dtype))
+ self.buff_idx.append(cp.empty(self.sparseIdx_max_vect_size, dtype=cp.int32))
+ self.buff_ptr.append(cp.empty(self.sparsePtr_max_vect_size, dtype=cp.int32))
+
+
+[docs] def allocate_dense_batch_buffers(self):
+ r"""
+ Allocate memory for dense batch buffers.
+
+ Based on the specified batching axis, this method initializes the appropriate
+ buffers to store dense data.
+ """
+
+ # Allocate memory based on column batching axis
+ if self.COL_BATCHING_AXIS:
+ self.H_d, self.X_d = [], []
+ for n in range(self.batchQeueSize):
+ self.H_d.append(cp.empty((self.k, self.batch_size), dtype=self.A_ij.dtype))
+ #self.X_d.append(cp.empty((self.grid_loc_m, self.batch_size), dtype=self.A_ij.dtype))
+ self.X_d.append(cp.empty((self.grid_glob_m, self.batch_size), dtype=self.A_ij.dtype))
+ else: # ROW_BATCHING_AXIS
+ self.W_d, self.X_d = [], []
+ for n in range(self.batchQeueSize):
+ self.W_d.append(cp.empty((self.batch_size, self.k), dtype=self.A_ij.dtype))
+ #self.X_d.append(cp.empty((self.batch_size, self.grid_loc_n), dtype=self.A_ij.dtype))
+ self.X_d.append(cp.empty((self.batch_size, self.grid_glob_n), dtype=self.A_ij.dtype))
+
+
+[docs] def allocate_gpu_batch_buffers(self):
+ r"""
+ Allocate memory for batch buffers on GPU.
+
+ Depending on whether the data is sparse or dense, it calls the appropriate
+ method to initialize the buffers.
+ """
+ # Check if the data is batched
+ if self.A_IS_BATCHED:
+ if self.A_IS_SPARSE:
+ self.allocate_sparse_batch_buffers()
+ else:
+ self.allocate_dense_batch_buffers()
+
+[docs] def showMemStats(self, msg=None):
+ r"""
+ Display memory statistics.
+
+ It logs information about used and available memory, including details about
+ the pinned memory pool.
+
+ Parameters:
+ msg (str): Optional message to include in the log.
+ """
+ # Only process with rank 0 displays the memory stats
+ if self.rank == 0:
+ if msg is not None: log(msg = f"[{amber('!')}][MEM INFO] {msg}", rank=self.rank, lrank=self.lrank)
+ # Log info based on whether the data is batched or not
+ if self.A_IS_BATCHED:
+ if msg is not None: log(msg = f"[{amber('!')}][MEM INFO] {msg}", rank=self.rank, lrank=self.lrank)
+ log(msg = f"[{amber('!')}][PINNEDMEM]: FREE BLOCKS = { amber(self.PINNEDMEMPOOL.n_free_blocks()) }", rank=self.rank, lrank=self.lrank)
+ else:
+ log(msg = f"[{amber('!')}][MEMPOOL]: USED = { amber(self.MEMPOOL.used_bytes()/self.GB)} GB || TOTAL = { red(self.MEMPOOL.total_bytes()/self.GB) } GB", rank=self.rank, lrank=self.lrank)
+ log(msg = f"[{amber('!')}][PINNEDMEM]: FREE BLOCKS = { amber(self.PINNEDMEMPOOL.n_free_blocks()) }", rank=self.rank, lrank=self.lrank)
+
+[docs] def sampleA(self, noise_var, method, seed=None):
+ r"""
+ Samples matrix A based on a specified method and noise variance.
+
+ Parameters:
+ - noise_var (float): Variance of the noise.
+ - method (str): Sampling method, either 'uniform' or 'poisson'.
+ - seed (int, optional): Random seed for reproducibility.
+
+ Raises:
+ - Exception if an unsupported sampling method is provided.
+ """
+ self.noise_var = noise_var
+ self.seed = seed
+ if self.seed != None:
+ cp.random.seed(self.seed)
+ self.sampling_method = method
+ if self.sampling_method == 'uniform':
+ self.randM()
+ elif self.sampling_method == 'poisson':
+ self.poisson()
+
+
+[docs] def randM(self):
+ r"""
+ Perturbs the elements of X by multiplying them with a uniform random number in the range (1-epsilon, 1+epsilon).
+ """
+
+ # Check if X_per has been initialized, if not, initialize the arrays
+ if not self.X_per_initialized: self.X_per, self.X_idx, self.X_ptr = [], [], []
+ #self.showMemStats(msg = " RandM() 0")
+
+ # Free up memory blocks in the memory pool
+ self.MEMPOOL.free_all_blocks()
+ #self.showMemStats(msg = " RandM() MEMPOOL Freed OK")
+
+ # Get the max available memory in bytes
+ MAXMEM = self.MEMPOOL.total_bytes()*1.0
+ if self.A_IS_BATCHED: #CPU Array
+ for b in range(self.nBatch):
+ if self.A_IS_SPARSE:
+ if self.ROW_BATCHING_AXIS:
+ i0, i1 = b*self.batch_size, (b+1)*self.batch_size
+ X = self.A_ij[i0:i1, :]
+ else: # Col Axis
+ j0, j1 = b*self.batch_size, (b+1)*self.batch_size
+ X = self.A_ij[:, j0:j1]
+ M = np.random.random_sample(X.data.shape).astype(self.A_ij.dtype)
+ M = 2.0*self.noise_var * M + self.noise_var + 1.0
+ X.data *= M
+ if not self.X_per_initialized: # Initialize sparse array pinned mem buffers
+ self.X_per.append(_pin_mem(np.append(X.data,np.zeros(self.sparseData_max_vect_size-len(X.data), dtype=X.data.dtype))))
+ self.X_idx.append(_pin_mem(np.append(X.indices,np.zeros(self.sparseIdx_max_vect_size-len(X.indices), dtype=X.indices.dtype))))
+ self.X_ptr.append(_pin_mem(np.append(X.indptr,np.zeros(self.sparsePtr_max_vect_size-len(X.indptr), dtype=X.indptr.dtype))))
+ else:
+ #print(f"Before : type(X_per) = {self.X_per[b].dtype}")
+ self.X_per[b] = np.append(X.data,np.zeros(self.sparseData_max_vect_size-len(X.data), dtype=X.data.dtype)).astype(self.A_ij.dtype)
+ self.X_idx[b] = np.append(X.indices,np.zeros(self.sparseIdx_max_vect_size-len(X.indices), dtype=X.indices.dtype))
+ self.X_ptr[b] = np.append(X.indptr,np.zeros(self.sparsePtr_max_vect_size-len(X.indptr), dtype=X.indptr.dtype))
+ #print(f"Afterpe : type(X_per) = {self.X_per[b].dtype}")
+
+ else: # A in Dense
+ if self.ROW_BATCHING_AXIS:
+ i0, i1 = b*self.batch_size, (b+1)*self.batch_size
+ M = np.random.random_sample(self.A_ij[i0:i1, :].shape).astype(self.A_ij.dtype)
+ M = 2.0*self.noise_var * M + self.noise_var + 1.0
+ if not self.X_per_initialized:
+ self.X_per.append(_pin_mem(np.multiply(self.A_ij[i0:i1, :], M)))
+ else:
+ self.X_per[b] = np.multiply(self.A_ij[i0:i1, :], M)
+ else: # Row Axis
+ j0, j1 = b*self.batch_size, (b+1)*self.batch_size
+ M = np.random.random_sample(self.A_ij[:, j0:j1].shape).astype(self.A_ij.dtype)
+ M = 2.0*self.noise_var * M + self.noise_var + 1.0
+ if not self.X_per_initialized:
+ self.X_per.append(_pin_mem(np.multiply(self.A_ij[:, j0:j1], M)))
+ else:
+ self.X_per[b] = np.multiply(self.A_ij[:, j0:j1], M)
+
+ else : #GPU Array
+ for b in range(self.nBatch):
+ if self.A_IS_SPARSE:
+ X = self.A_d[b]*1.0
+ M = cp.random.random_sample(self.A_d[b].data.shape, dtype=self.A_d[b].dtype)
+ M = 2.0*self.noise_var * M + self.noise_var + 1.0
+ X.data *= M
+ #self.X_per.append(X)
+ if not self.X_per_initialized:
+ self.X_per.append(X)
+ else:
+ self.X_per[b] = X
+ else: # A in Dense
+ self.showMemStats(msg = " RandM() BEFORE X_per init b%04d" %b)
+ MAXMEM = max(MAXMEM, self.MEMPOOL.total_bytes()*1.0)
+ if not self.X_per_initialized:
+ self.X_per.append(cp.random.random_sample(self.A_d[b].shape, dtype=self.A_d[b].dtype))
+ self.showMemStats(msg = " RandM() AFTER X_per init 1 b%04d" %b)
+ MAXMEM = max(MAXMEM, self.MEMPOOL.total_bytes()*1.0)
+ else:
+ self.X_per[b] = cp.random.random_sample(self.A_d[b].shape, dtype=self.A_d[b].dtype)
+ self.showMemStats(msg = " RandM() AFTER X_per set 1 b%04d" %b)
+ MAXMEM = max(MAXMEM, self.MEMPOOL.total_bytes()*1.0)
+ self.X_per[b] = 2.0*self.noise_var
+ self.X_per[b] += self.noise_var + 1.0
+ self.X_per[b] = cp.multiply(self.A_d[b], self.X_per[b])
+ self.showMemStats(msg = " RandM() X_per init OK b%04d" %b)
+ if self.rank == 0: log(msg = f"[{amber('!')}][MEMPOOL]: randM(): MEM Peak = {blue(MAXMEM/self.GB)} GB", rank=self.rank, lrank=self.lrank)
+ self.MEMPOOL.free_all_blocks()
+
+ if not self.X_per_initialized: self.X_per_initialized = True
+
+
+[docs] def init_factors(self):
+ r"""
+ Initializes NMF factor matrices W and H based on the selected initialization method.
+ """
+ # Random initialization
+ if self.params.init == 'rand':
+ if self.topo_dim=='2d':
+ #raise Exception(f"[{red('!!')}] 2D topology is not yet supported on GPU")
+ noSupport("2D topology")
+ elif self.topo_dim in ['1d']:
+ rng_h = np.random.RandomState(self.iter_seeds[self.k])
+ rng_d = cp.random.RandomState(self.iter_seeds[self.k])
+ if self.ROW_BATCHING_AXIS:
+ self.H_d = cp.random.rand(self.k, self.grid_glob_n).astype(self.A_ij.dtype) # H [k, J]
+ self.H_is_cached = True
+ if self.A_IS_BATCHED: # W is distributed, H is replicate on Hostd
+ self.W_h = []
+ for n in range(self.nBatch):
+ self.W_h.append(_pin_mem(rng_h.rand(self.batch_size, self.k).astype(self.A_ij.dtype)))
+ self.W_is_cached = False
+ else:
+ self.W_d = rng_d.rand(self.grid_glob_m, self.k).astype(self.A_ij.dtype) # W_h [m, k]
+ self.W_is_cached = True
+ else: # COL_BATCHING_AXIS:
+ self.W_d = cp.random.rand(self.grid_glob_m, self.k).astype(self.A_ij.dtype) # W_d [I, k]
+ self.W_is_cached = True
+ if self.A_IS_BATCHED: # H is distributed, W is replicate on Hostd
+ self.H_h = []
+ for n in range(self.nBatch):
+ self.H_h.append(_pin_mem(rng_h.rand(self.k, self.batch_size).astype(self.A_ij.dtype)))
+ self.H_is_cached = False
+ else:
+ self.H_d = rng_d.rand(self.k, self.grid_glob_n).astype(self.A_ij.dtype) # W_h [m, k]
+ self.H_is_cached = True
+ # NNSVD-based initialization
+ elif self.params.init in ['nnsvd','svd','singular']:
+ if self.topo_dim == '1d':
+ self.AA = []
+ if self.A_IS_BATCHED:
+ for i in range(self.nGPUs):
+ self.AA.append(np.zeros(self.grid_L, self.grid_L))
+ else:
+ for i in range(self.nGPUs):
+ self.AA.append(_zeros((self.grid_L, self.grid_L)))
+ self.svdSoFar = [[1,0,0]]
+ self.W_i, self.H_j = self.nnsvd(flag=1, verbose=0)
+ elif self.topo_dim == '2d':
+ raise Exception('NNSVD init only available for 1D topology, please try with 1d topo.')
+ else:
+ # Raise an error if an unsupported initialization method is chosen
+ raise Exception(f"[{red('!!')}] Only Random and nnsvd init is supported on GPU")
+
+
+
+
+[docs] def normalize_features(self):
+ r"""
+ Normalizes the NMF factor matrices W and H.
+ """
+ # If A is batched
+ if self.A_IS_BATCHED: # W.I.P
+ if self.COL_PARTITION: # W.I.P
+ pass
+ else:
+ Wall_norm = self.W_d.sum(axis=0, keepdims=True)
+ if self.params.p_r != 1: self.nccl_comm.Allreduce(Wall_norm, Wall_norm, op=nccl.NCCL_SUM, stream=None)
+ self.W_d/= Wall_norm+ self.eps
+ self.H_h *= cp.asnumpy(Wall_norm.T)
+
+ else: # If A is not batched
+ Wall_norm = self.W_d.sum(axis=0, keepdims=True)
+ if self.topo_dim == '2d':
+ raise Exception(f"[{red('!!')}] 2D topology is not yet supported on GPU")
+ elif self.topo_dim == '1d':
+ if self.params.p_r != 1:
+ self.nccl_comm.Allreduce(Wall_norm, Wall_norm, op=nccl.NCCL_SUM, stream=None)
+ else:
+ raise Exception(f"[{red('!!')}] Topology {self.topo_dim}is not yet supported")
+ self.W_d /= Wall_norm+ self.eps
+ self.H_d *= Wall_norm.T
+
+
+[docs] def relative_err(self):
+ r"""
+ Computes the relative reconstruction error of the NMF decomposition.
+ """
+
+ # Check topology and raise an error for unsupported topologies
+ if self.topo_dim == '2d': raise Exception(f"[{red('!!')}] 2D topology is not yet supported on GPU")
+ err = 0.0
+ # Get the max available memory in bytes
+ MAXMEM = self.MEMPOOL.total_bytes()*1.0
+
+ # Compute the error for batched CPU arrays
+ if self.A_IS_BATCHED: #CPU Array
+ for b in range(self.nBatch):
+ if self.ROW_BATCHING_AXIS: # ROW AXIS
+ self.H_h = cp.asnumpy(self.H_d) # Download H
+ if self.A_IS_SPARSE:
+ err += np.square( self.X_per[b].toarray() - (self.W_h[b] @ self.H_h) ).sum()
+ else:
+ err += np.square( self.X_per[b] - (self.W_h[b] @ self.H_h) ).sum()
+ else: # COL AXIS
+ self.W_h = cp.asnumpy(self.W_d) # Downlod W
+ if self.A_IS_SPARSE:
+ #err += np.square( self.X_per[b].toarray() - (self.W_h @ self.H_h[b]) ).sum()
+ err += 0.0
+ else:
+ err += np.square( self.X_per[b] - (self.W_h @ self.H_h[b]) ).sum()
+ self.glob_norm_A = self.dist_norm(self.A_ij)
+ else: # Compute the error for GPU arrays
+ #self.showMemStats(msg = " relative_err() STARTING BATCHES")
+ MAXMEM = max(MAXMEM, self.MEMPOOL.total_bytes()*1.0)
+ for b in range(self.nBatch):
+ if self.ROW_BATCHING_AXIS: # Row Partition
+ i0, i1 = b*self.grid_I, (b+1)*self.grid_I
+ if self.A_IS_SPARSE:
+ err += cp.square( self.A_d[b].toarray() - (_matmul(self.W_d[i0:i1], self.H_d)) ).sum()
+ else:
+
+ err += cp.square( self.A_d[b] - (_matmul(self.W_d[i0:i1], self.H_d)) ).sum()
+ MAXMEM = max(MAXMEM, self.MEMPOOL.total_bytes()*1.0)
+ self.MEMPOOL.free_all_blocks()
+
+ else: # Col Partition
+ j0, j1 = b*self.grid_J, (b+1)*self.grid_J
+ if self.A_IS_SPARSE:
+ err += cp.square( self.A_d[b].toarray() - (_matmul(self.W_d, self.H_d[:, j0:j1])) ).sum()
+ else:
+ err += cp.square( self.A_d[b] - (_matmul(self.W_d, self.H_d[:, j0:j1])) ).sum()
+ err = err.get()
+ self.glob_norm_A = self.dist_norm(self.A_d)
+ #self.glob_norm_A = self.dist_norm(self.X_per)
+ #del err
+ err= self.sub_comm.allreduce(err)
+ self.glob_norm_err = np.sqrt(err)
+ self.recon_err = self.glob_norm_err / self.glob_norm_A
+ if self.rank == 0: log(msg = f"[{amber('!')}][MEMPOOL]: relative_err(): MEM Peak = {blue(MAXMEM/self.GB)} GB", rank=self.rank, lrank=self.lrank)
+
+[docs] def dist_norm(self, X, proc=-1, norm='fro', axis=None):
+ r"""Computes the distributed norm of an array.
+
+ Args:
+ - X (list/array): The input array or list of arrays.
+ - proc (int, optional): Processor ID. Default is -1.
+ - norm (str, optional): Type of matrix norm. Default is 'fro' (Frobenius norm).
+ - axis (optional): Axis for norm computation. Default is None.
+
+ Returns:
+ - float: The computed distributed norm.
+ """
+
+ # If X is a list, this is a batched computation.
+ if type(X) in [list]: # BACHED
+ array_type = self.getArrayType(X[0])
+ N = len(X)
+ err = 0.0
+ try:
+ device, density = array_type.split('_')
+ except:
+ raise Exception('[!][dist_norm()] UNABLE TO IDENTIFY ARRAY OBJECT TYPE {}!!'.format(array_type))
+ if device in ['CPU']:
+ if density in ['DENSE']:
+ for n in range(N): err += np.square(X[n]).sum()
+ elif density in ['SPARSE']:
+ for n in range(N): err += np.square(X[n].data).sum()
+ else:
+ raise Exception('[!][dist_norm()] UNABLE TO IDENTIFY ARRAY OBJECT TYPE {}!!'.format(array_type))
+ elif device in ['GPU']:
+ if density in ['DENSE']:
+ for n in range(N): err += cp.square(X[n]).sum()
+ elif density in ['SPARSE']:
+ for n in range(N): err += cp.square(X[n].data).sum()
+ else:
+ raise Exception('[!][dist_norm()] UNABLE TO IDENTIFY ARRAY OBJECT TYPE {}!!'.format(array_type))
+ err = err.get()
+
+ else:
+ raise Exception('[!][dist_norm()] UNABLE TO IDENTIFY ARRAY OBJECT TYPE {}!!'.format(array_type))
+ else: # If X is not a list, this is a non-batched computation.
+ array_type = self.getArrayType(X)
+ try:
+ device, density = array_type.split('_')
+ except:
+ raise Exception('[!][dist_norm()] UNABLE TO IDENTIFY ARRAY OBJECT TYPE {}!!'.format(array_type))
+ if device in ['CPU']:
+ if density in ['DENSE']:
+ err=np.square(X).sum()
+ elif density in ['SPARSE']:
+ err=np.square(X.data).sum()
+ elif device in ['GPU']:
+ if density in ['DENSE']:
+ err = cp.square(X).sum()
+ elif density in ['SPARSE']:
+ err = cp.square(X.data).sum()
+ else:
+ raise Exception('[!][dist_norm()] UNABLE TO IDENTIFY ARRAY OBJECT TYPE {}!!'.format(array_type))
+ err = err.get()
+ else:
+ raise Exception('[!][dist_norm()] UNABLE TO IDENTIFY ARRAY OBJECT TYPE {}!!'.format(array_type))
+
+ err = self.sub_comm.allreduce(err)
+ return np.sqrt(err)
+
+
+[docs] def fit(self, factors=None, save_factors=False, rtrn=['W','H']):
+ r"""Fits the NMF model based on provided parameters.
+
+ Args:
+ - factors (list, optional): Initial factor matrices. Default is None.
+ - save_factors (bool, optional): Flag to decide whether to save factor matrices. Default is False.
+ - rtrn (list, optional): List of matrix names to return. Default is ['W','H'].
+
+ Returns:
+ - tuple: Factor matrices and meta information.
+ """
+ if self.params.norm.upper() == 'FRO':
+ if self.params.method.upper() == 'MU':
+ self.Fro_MU_update(factors=factors, save_factors=save_factors, rtrn=rtrn)
+ else:
+ raise Exception('Not a valid method: Choose (mu)')
+ elif self.params.norm.upper() == 'KL':
+ if self.params.method.upper() == 'MU':
+ self.KL_MU_update(factors=factors, save_factors=save_factors, rtrn=rtrn)
+ else:
+ raise Exception('Not a valid method: Choose (mu)')
+ else:
+ raise Exception('Not a valid norm: Choose (fro/kl)')
+ #return self.W_i, self.H_j
+ metha = {}
+ metha['dt'] = self.dt
+ metha['err'] = self.recon_err
+ if self.A_IS_BATCHED:
+ if self.ROW_BATCHING_AXIS:
+ return self.W_h, cp.asnumpy(self.H_d), metha
+ else:
+ return cp.asnumpy(self.W_d), self.H_h, metha
+ #return self.W_h, self.H_h, self.recon_err
+ else:
+ return cp.asnumpy(self.W_d), cp.asnumpy(self.H_d), metha #self.recon_err.get()
+
+
+[docs] def Fro_MU_update(self, factors=None, save_factors=False, rtrn=['W','H']):
+ r"""Performs updates for NMF using the Frobenius norm and MU optimization method.
+
+ Args:
+ - factors (list, optional): Initial factor matrices. Default is None.
+ - save_factors (bool, optional): Flag to decide whether to save factor matrices. Default is False.
+ - rtrn (list, optional): List of matrix names to return. Default is ['W','H'].
+
+ Returns:
+ - None
+ """
+ # Memory management and initialization steps.
+ self.showMemStats(msg = " Fro_MU_update() 0")
+ self.MEMPOOL.free_all_blocks()
+ #self.showMemStats(msg = " Fro_MU_update() MEMPOOL Freed OK")
+ self.PINNEDMEMPOOL.free_all_blocks()
+ # If the matrix is batched, perform the updates using the appropriate method based on partitioning axis.
+ if self.A_IS_BATCHED: # BATCHED CPU -> GPU
+ if self.COL_BATCHING_AXIS:
+ self.FroNMF_1D_row(factors=factors, save_factors=save_factors, rtrn=rtrn)
+ elif self.ROW_BATCHING_AXIS:
+ #self.FroNMF_1D_col(factors=factors, save_factors=save_factors, rtrn=rtrn)
+ self.FroNMF_1D_row_batched(factors=factors, save_factors=save_factors, rtrn=rtrn)
+ else:
+ raise Exception(f"[{red('!!')}] Grid Partition UNDEFINED!! ")
+ else: # LOCAL BATCHING on GPU Mem
+ if self.COL_BATCHING_AXIS:
+ self.FroNMF_1D_col_partion(factors=factors, save_factors=save_factors, rtrn=rtrn)
+ elif self.ROW_BATCHING_AXIS:
+ self.FroNMF_1D_row_partion(factors=factors, save_factors=save_factors, rtrn=rtrn)
+ else:
+ raise Exception(f"[{red('!!')}] Grid Partition UNDEFINED!! ")
+
+
+[docs] def FroNMF_1D_col_partion(self,factors=None, save_factors=False, rtrn=False): #W.I.P
+ r"""Performs NMF on 1D column partitioned data using the Frobenius norm.
+
+ Args:
+ - factors (list, optional): Initial factor matrices. Default is None.
+ - save_factors (bool, optional): Flag to decide whether to save factor matrices. Default is False.
+ - rtrn (bool, optional): Flag to decide whether to return factor matrices. Default is False.
+
+ Returns:
+ - None
+ """
+ #[1] Initialize perturbed cofactors
+ self.init_factors()
+ #[2] Initialiaze different accumulators
+ self.events['nmf_start'].record()
+ self.dt['W_up'], self.dt['H_up'], self.dt['allRed_XHT'], self.dt['allRed_WHHT'] = 0.0, 0.0, 0.0, 0.0
+ for i in range(self.params.itr):
+ #[3] Initialize iteration Accumulators
+ dt_w, dt_h, dt_reduce1, dt_reduce2 = 0.0, 0.0, 0.0, 0.0
+ WT = self.W_d.T # [m, k].T -> WT [k, m]
+ WTW = _matmul(WT, self.W_d) # [k, m] @ [m, k] -> WTW [k, k]
+ WHHT = _zeros((self.grid_glob_m, self.k), dtype=self.A_ij.dtype) # WHHT [m, k] # ACCUMULATOR for W@H@H.T
+ XHT = _zeros((self.grid_glob_m, self.k), dtype=self.A_ij.dtype) # XHT [m, k] # ACCUMULATOR for X@H.T
+ #[4] Strart local batching
+ for b in range(self.nBatch):
+ j0, j1 = b*self.batch_size, (b+1)*self.batch_size
+ #[5]//////////////////////////// (H update)> ////////////////////////////////////////////////////////
+ self.events['start'].record()
+ WTWH = _matmul(WTW, self.H_d[:, j0:j1]) + self.eps # [k, k] @ [k, b] -> WTWH [k, b]
+ if self.A_IS_SPARSE:
+ WTX = self.MM(self.X_per[p], self.W_d, transpA=True).T # ([m, b].T @ [m, k]).T -> [b, k].T -> WTX [k, b]
+ else:
+ WTX = _matmul(WT, self.X_per[b]) # [k, m] @ [m, b] -> XHT [k, b]
+ WTX = _multiply(WTX, self.H_d[:, j0:j1]) # [k, b] * [k, b] -> HWTX [k, b]
+ self.H_d[:, j0:j1] = _divide(WTX, WTWH) # [k, b] / [k, b] -> H [k, b]
+ self.events['end'].record()
+ self.events['end'].synchronize()
+ #del WTWH, WTX
+ #self.MEMPOOL.free_all_blocks()
+ dt_h += cp.cuda.get_elapsed_time(self.events['start'], self.events['end']) # [ms]
+ #[6]//////////////////////////// (W update)> ////////////////////////////////////////////////////////
+ if self.params.W_update:
+ self.events['start'].record()
+ HT = self.H_d[:, j0:j1].T # [k, b].T -> HT [b, k]
+ HHT = _matmul(self.H_d[:, j0:j1], HT) # [k, b] @ [b, k] -> HTH [k, k]
+ WHHT += _matmul(self.W_d, HHT) # [m, k] @ [k, k] -> WHHT [m, k]
+ if self.A_IS_SPARSE:
+ XHT += self.MM(self.X_per[p], HT, transpA=False) # [b, n].T @ [b, k] -> [n, k].T -> [k, n]
+ else:
+ XHT += _matmul(self.self.X_per[b], HT) # [m, b] @ [b, k] -> xht [m, k]
+ self.events['end'].record()
+ self.events['end'].synchronize()
+ #del HT, HHT
+ #self.MEMPOOL.free_all_blocks()
+ dt_w += cp.cuda.get_elapsed_time(self.events['start'], self.events['end']) # [ms]
+ #[7] Aggregate Global WHHT and XHT
+ self.events['start'].record()
+ if self.nGPUs > 1:
+ self.events['reduce_s'].record()
+ nccl_comm.Allreduce(WHHT, WHHT, op=nccl.NCCL_SUM, stream=None)
+ self.events['reduce_e'].record()
+ self.events['reduce_e'].synchronize()
+ dt_reduce1 += cp.cuda.get_elapsed_time(self.events['reduce_s'], self.events['reduce_e'])
+ self.events['reduce_s'].record()
+ nccl_comm.Allreduce(XHT, XHT, op=nccl.NCCL_SUM, stream=None)
+ self.events['reduce_e'].record()
+ self.events['reduce_e'].synchronize()
+ dt_reduce2 += cp.cuda.get_elapsed_time(self.events['reduce_s'], self.events['reduce_e'])
+ #[8] Finalize W-Update
+ WHHT = WHHT + self.eps
+ XHT = _multiply(self.W_d, XHT) # [m, k] * [m, k] -> [m, k]
+ self.W_d = _divide(XHT, WHHT) # [m, k] / [m, k] -> [m, k]
+ self.events['end'].record()
+ self.events['end'].synchronize()
+ dt_w += cp.cuda.get_elapsed_time(self.events['start'], self.events['end']) # [ms]
+ #[9] Accumulate iteration performance metric
+ self.dt['W_up'] += dt_w*1.0
+ self.dt['H_up'] += dt_h*1.0
+ self.dt['allRed_WHHT'] += dt_reduce1*1.0
+ self.dt['allRed_XHT'] += dt_reduce2*1.0
+ #[10]
+ if i % 10 == 0:
+ #if self.A_is_Large:
+ # self.H_h = np.maximum(self.H_h, self.eps)
+ #else:
+ # self.H_d = cp.maximum(self.H_d, self.eps)
+ self.H_d = cp.maximum(self.H_d, self.eps)
+ self.W_d = cp.maximum(self.W_d, self.eps)
+
+ #[10] Compute Average performance metric
+ self.dt['W_up'] /= float(self.params.itr)
+ self.dt['H_up'] /= float(self.params.itr)
+ self.dt['allRed_XHT'] /= float(self.params.itr)
+ self.dt['allRed_WHHT'] /= float(self.params.itr)
+ self.events['nmf_end'].record()
+ self.events['nmf_end'].synchronize()
+ self.dt['NMF'] = cp.cuda.get_elapsed_time(self.events['nmf_start'], self.events['nmf_end']) # [ms]
+ #[11] Calculate relative error
+ self.relative_err()
+
+
+
+[docs] def FroNMF_1D_row_partion(self,factors=None, save_factors=False, rtrn=False):
+ r"""
+ This method performs Frobenius Norm Non-negative Matrix Factorization (FroNMF)
+ on a 1D row partitioned matrix.
+
+ Parameters:
+ - factors: Initial guess for the factors (Optional)
+ - save_factors: Boolean, whether to save the computed factors or not (Default is False)
+ - rtrn: Boolean, whether to return the factors (Default is False)
+ """
+
+ # [1] Monitor initial memory usage
+ self.showMemStats(msg = " FroNMF_1D_row_partion() 0")
+
+ # [2] Initialization of perturbed cofactors
+ self.init_factors()
+ #self.showMemStats(msg = " FroNMF_1D_row_partion() init_factors OK")
+ # [3] Initialize various accumulators for performance measurement and computation
+
+ self.events['nmf_start'].record()
+ self.dt['W_up'], self.dt['H_up'], self.dt['allRed_WTX'], self.dt['allRed_WTWH'] = 0.0, 0.0, 0.0, 0.0
+
+ # Main iterative loop for NMF computation
+ for i in range(self.params.itr):
+ # Initialize accumulators for this iteration
+ dt_w, dt_h, dt_reduce1, dt_reduce2 = 0.0, 0.0, 0.0, 0.0
+
+ # [4] Setup for W matrix update
+ if self.params.W_update:
+ HT = self.H_d.T # HT [n, k]
+ HHT = _matmul(self.H_d, HT) # [k, n] @ [n, k] -> HHT [k, k]
+
+ # Initialize accumulators for batch computations
+ WTWH = _zeros((self.k, self.grid_glob_n), dtype=self.A_ij.dtype) # WTWH [k, n] # ACCUMULATOR for WT@W@H
+ WTX = _zeros((self.k, self.grid_glob_n), dtype=self.A_ij.dtype) # WTX [k, n] # ACCUMULATOR for WT@A
+
+ # [5] Batching loop for local computations
+ for b in range(self.nBatch):
+ # Determine current batch range
+ i0, i1 = b*self.batch_size, (b+1)*self.batch_size
+ #self.showMemStats(msg = " FroNMF_1D_row_partion() 2 b%04d" %b)
+ #[6]//////////////////////////// (W update)> ////////////////////////////////////////////////////////
+ if self.params.W_update:
+ self.events['start'].record()
+ WHHT = _matmul(self.W_d[i0:i1,:], HHT) + self.eps # [b, k] @ [k, k] -> WHHT [b, k]
+ #self.showMemStats(msg = " FroNMF_1D_row_partion() 3 b%04d" %b)
+ if self.A_IS_SPARSE:
+ XHT = self.MM(self.X_per[b], HT, transpA=False) # [b, n] @ [n, k] -> XHT [b, k]
+ else:
+ XHT = _matmul(self.X_per[b], HT) # [b, n] @ [n, k] -> XHT [b, k]
+ #self.showMemStats(msg = " FroNMF_1D_row_partion() 4 b%04d" %b)
+ XHT = _multiply(self.W_d[i0:i1,:], XHT) # [b, k] * [b, k] -> WXHT [b, k]
+ #self.showMemStats(msg = " FroNMF_1D_row_partion() 5 b%04d" %b)
+ self.W_d[i0:i1,:] = _divide(XHT, WHHT) # [b, k] / [b, k] -> W [b, k]
+ #self.showMemStats(msg = " FroNMF_1D_row_partion() W-up b%04d OK" %b)
+ self.events['end'].record()
+ self.events['end'].synchronize()
+ del WHHT, XHT
+ self.MEMPOOL.free_all_blocks()
+ dt_w += cp.cuda.get_elapsed_time(self.events['start'], self.events['end']) # [ms]
+ #[7]//////////////////////////// (H update)> ////////////////////////////////////////////////////////
+ #self.showMemStats(msg = " FroNMF_1D_row_partion() 7")
+ self.events['start'].record()
+ WT = self.W_d[i0:i1,:].T # [b, k].T -> WT [k, b]
+ #self.showMemStats(msg = " FroNMF_1D_row_partion() 8")
+ WTW = _matmul(WT, self.W_d[i0:i1,:]) # [k, b] @ [b, k] -> WTW [k, k]
+ #self.showMemStats(msg = " FroNMF_1D_row_partion() 9")
+ #print(f"[!!!][iter %04d | batch %04d] shape(WTW)={WTW.shape} || shape(H_d) = {self.H_d.shape}" %(i, b))
+ WTWH += _matmul(WTW, self.H_d) # [k, k] @ [k, n] -> WTWH [k, n]
+ #self.showMemStats(msg = " FroNMF_1D_row_partion() 10 ")
+ if self.A_IS_SPARSE:
+ WTX += self.MM(self.X_per[b], self.W_d[i0:i1,:], transpA=True).T # [b, n].T @ [b, k] -> [n, k].T -> [k, n]
+ else:
+ WTX += _matmul(WT, self.X_per[b]) # [k, b] @ [b, n] -> WTX [k, n]
+ #self.showMemStats(msg = " FroNMF_1D_row_partion() H-up B%04d OK" %b)
+ self.events['end'].record()
+ self.events['end'].synchronize()
+ # Free memory from intermediate variables
+ del WT, WTW
+ self.MEMPOOL.free_all_blocks()
+ dt_h += cp.cuda.get_elapsed_time(self.events['start'], self.events['end']) # [ms]
+
+ # Free memory from intermediate variable
+ del HHT
+ self.MEMPOOL.free_all_blocks()
+ # [8] Aggregate global WTWH and WTX values (used for multi-GPU setup)
+ self.events['start'].record()
+ if self.nGPUs > 1:
+ # Perform Allreduce operation on WTWH and WTX
+ # This is to aggregate results from all GPUs
+
+ self.events['reduce_s'].record()
+ self.nccl_comm.Allreduce(WTWH, WTWH, op=nccl.NCCL_SUM, stream=None)
+ self.events['reduce_e'].record()
+ self.events['reduce_e'].synchronize()
+ dt_reduce1 += cp.cuda.get_elapsed_time(self.events['reduce_s'], self.events['reduce_e'])
+ self.events['reduce_s'].record()
+ self.nccl_comm.Allreduce(WTX, WTX, op=nccl.NCCL_SUM, stream=None)
+ self.events['reduce_e'].record()
+ self.events['reduce_e'].synchronize()
+ dt_reduce2 += cp.cuda.get_elapsed_time(self.events['reduce_s'], self.events['reduce_e'])
+ # [9] Update H with global results
+ WTWH = WTWH + self.eps
+ #self.showMemStats(msg = " FroNMF_1D_row_partion() 12")
+ WTX = _multiply(self.H_d, WTX) # [k, n] * [k, n] -> [k, n]
+ #self.showMemStats(msg = " FroNMF_1D_row_partion() 13")
+ self.H_d = _divide(WTX, WTWH) # [k, n] / [k, n] -> [k, n]
+ #self.showMemStats(msg = " FroNMF_1D_row_partion() NMF OK")
+ self.events['end'].record()
+ self.events['end'].synchronize()
+ dt_h += cp.cuda.get_elapsed_time(self.events['start'], self.events['end']) # [ms]
+ # Accumulate iteration performance metric
+ self.dt['W_up'] += dt_w*1.0
+ self.dt['H_up'] += dt_h*1.0
+ self.dt['allRed_WTX'] += dt_reduce1*1.0
+ self.dt['allRed_WTWH'] += dt_reduce2*1.0
+ #[10] Clipping to ensure non-negativity after every 10 iterations
+ if i % 10 == 0:
+ #if self.A_is_Large:
+ # self.H_h = np.maximum(self.H_h, self.eps)
+ #else:
+ # self.H_d = cp.maximum(self.H_d, self.eps)
+ self.H_d = cp.maximum(self.H_d, self.eps)
+ self.W_d = cp.maximum(self.W_d, self.eps)
+ self.showMemStats(msg = " FroNMF_1D_row_partion() NMF OK")
+ # [11] Compute Average Performance metrics over all iterations
+ del WTWH, WTX
+ self.MEMPOOL.free_all_blocks()
+ self.showMemStats(msg = " FroNMF_1D_row_partion() ALL DELETIONS DONE")
+ self.dt['W_up'] /= float(self.params.itr)
+ self.dt['H_up'] /= float(self.params.itr)
+ self.dt['allRed_WTX'] /= float(self.params.itr)
+ self.dt['allRed_WTWH'] /= float(self.params.itr)
+ self.events['nmf_end'].record()
+ self.events['nmf_end'].synchronize()
+ self.dt['NMF'] = cp.cuda.get_elapsed_time(self.events['nmf_start'], self.events['nmf_end']) # [ms]
+ # [12] Calculate relative error
+ self.showMemStats(msg = " FroNMF_1D_row_partion() BEFORE ERROR Eval")
+ if self.A_IS_SPARSE:
+ self.recon_err = 0.0
+ else:
+ self.relative_err()
+ self.showMemStats(msg = " FroNMF_1D_row_partion() AFTER ERROR Eval")
+ # Free all memory blocks
+ self.MEMPOOL.free_all_blocks()
+ self.showMemStats(msg = " FroNMF_1D_row_partion() LAST CALL MEM POOL FReed OK")
+
+
+
+[docs] def FroNMF_1D_row_batched(self,factors=None, save_factors=False, rtrn=False): #W.I.P
+ r"""
+ Computes the Frobenius Norm Non-negative Matrix Factorization (FroNMF) for the 1D row of the data matrix
+ using a batched approach.
+
+ Args:
+ - factors: Not used in this function but could represent the initial factors for matrix factorization.
+ - save_factors: A boolean flag to determine if factors should be saved.
+ - rtrn: A boolean flag to determine if the function should return some values (not implemented yet).
+
+ Returns:
+ - None (but updates internal attributes of the object and potentially could return some values based on `rtrn` flag).
+ """
+ #[1] Initialize perturbed cofactors
+ self.init_factors()
+ #[2] Initialize accumulators and record the start event for NMF
+ self.events['nmf_start'].record()
+ self.dt['W_up'], self.dt['H_up'], self.dt['allRed_WTX'], self.dt['allRed_WTWH'] = 0.0, 0.0, 0.0, 0.0
+ self.dt['H2D_A'], self.dt['H2D_W'], self.dt['D2H_W'] = 0.0, 0.0, 0.0
+ SQ = self.get_managed_stream_queue() # Get Managed stream Queue
+
+ # Loop for the number of iterations (NMF optimization steps)
+ for i in range(self.params.itr):
+ #[3] Initialize iteration Accumulators
+ dt_w, dt_h, dt_reduce1, dt_reduce2 = 0.0, 0.0, 0.0, 0.0
+ if self.params.W_update:
+ HT = self.H_d.T # HT [n, k]
+ HHT = _matmul(self.H_d, HT) # [k, n] @ [n, k] -> HHT [k, k]
+ #WTWH = _zeros((self.k, self.grid_glob_n), dtype=self.A_ij.dtype) # WTWH [k, n] # ACCUMULATOR for WT@W@H
+ ##WTX = _zeros((self.k, self.grid_glob_n), dtype=self.A_ij.dtype) # WTX [k, n] # ACCUMULATOR for WT@A
+ WTX, WTWH = [], []
+ H2D_A, H2D_W, D2H_W, W_up, H_up, AR_1, AR_2 = [], [], [], [], [], [], []
+ #[4] Start processing data in batches
+ for b in range(self.nBatch):
+ st_key = SQ['Queue'].get() # Get stream from the queue. If queue is empty, wait until a stream is available.
+ stream = SQ['READY'][st_key]
+ i0, i1 = b*self.batch_size, (b+1)*self.batch_size
+ #strm_idx = b % self.MAXSTREAMS
+ #strm_name = self.stream[strm_idx]
+ #stream = self.STREAMS[strm_name]
+ with stream:
+ Q_idx = b % self.batchQeueSize
+ #stream.wait_event(self.FREE_STREAM[strm_name])
+ self.events['h2d_s'].record()
+ self.W_d[Q_idx].set(self.W_h[b])
+ self.events['h2d_e'].record()
+ self.events['h2d_e'].synchronize()
+ H2D_W.append( cp.cuda.get_elapsed_time(self.events['h2d_s'], self.events['h2d_e']) ) # [ms]
+ self.events['h2d_s'].record()
+ if self.A_IS_SPARSE: # CopyH2D_Async Sparse array
+ #print(f"[!!] [itr %04d][b %02d] type(X_per[b])= {self.X_per[b].dtype}" %(i,b))
+ self.X_d[Q_idx].set(self.X_per[b].astype(self.A_ij.dtype))
+ self.buff_idx[Q_idx].set(self.X_idx[b])
+ self.buff_ptr[Q_idx].set(self.X_ptr[b])
+ else:
+ self.X_d[Q_idx].set(self.X_per[b]) # CopyH2D_Async
+ self.events['h2d_e'].record()
+ self.events['h2d_e'].synchronize()
+ #print(f"[!!] [itr %04d][b %02d] H2D_A = {cp.cuda.get_elapsed_time(self.events['h2d_s'], self.events['h2d_e'])}" %(i,b))
+ H2D_A.append( cp.cuda.get_elapsed_time(self.events['h2d_s'], self.events['h2d_e']) ) # [ms]
+
+ #[5]//////////////////////////// (W update)> ////////////////////////////////////////////////////////
+ if self.params.W_update:
+ self.events['start'].record()
+ WHHT = _matmul(self.W_d[Q_idx], HHT) + self.eps # [b, k] @ [k, k] -> WHHT [b, k]
+ if self.A_IS_SPARSE:
+ #XHT = self.MM(self.X_per[p], HT, transpA=False) # [b, n] @ [n, k] -> XHT [b, k]
+ #print(f"[!!] BS = {self.batch_size} || vect_size_dat = {self.sparse_vect_size['dat'][b]}, vect_size_idx = {self.sparse_vect_size['idx'][b]}, vect_size_ptr = {self.sparse_vect_size['ptr'][b]},")
+ X_d = cupyx.scipy.sparse.csr_matrix((self.X_d[Q_idx][:self.sparse_vect_size['dat'][b]],
+ self.buff_idx[Q_idx][:self.sparse_vect_size['idx'][b]],
+ self.buff_ptr[Q_idx][:self.sparse_vect_size['ptr'][b]]), shape=(self.batch_size, self.grid_loc_n))
+ #print(f"[!!] shape(X_d) ={X_d.shape} || shape(HT) = {HT.shape} || BS = {self.batch_size} || vect_size_dat = {self.sparse_vect_size['dat'][b]}")
+ #XHT = self.MM(X_d, HT, transpA=False) # [b, n] @ [n, k] -> XHT [s, k]
+ XHT = X_d @ HT
+ else:
+ #XHT = _matmul(self.X_per[b], HT) # [b, n] @ [n, k] -> XHT [b, k]
+ XHT = self.MM(self.X_d[Q_idx], HT) # [b, n] @ [n, k] -> XHT [s, k]
+ #XHT = _multiply(self.W_d[i0:i1,:], XHT) # [b, k] * [b, k] -> WXHT [b, k]
+ #self.W_d[i0:i1,:] = _divide(XHT, WHHT) # [b, k] / [b, k] -> W [b, k]
+ self.W_d[Q_idx] = _divide( _multiply(self.W_d[Q_idx], XHT), _matmul(self.W_d[Q_idx],HHT)+self.params.eps) # [bs,k] <- ([bs,k]*[bs,k]) / ([bs,k]@[k, k] + epsi)
+ #WT = self.W_d[Q_idx].T # [k, b]
+ self.events['d2h_s'].record()
+ self.W_d[Q_idx].get(out=self.W_h[b]) # CopyD2H_Async W_d -> W_h
+ self.events['d2h_e'].record()
+ self.events['d2h_e'].synchronize()
+ D2H_W.append( cp.cuda.get_elapsed_time(self.events['d2h_s'], self.events['d2h_e']) ) # [ms]
+ self.events['end'].record()
+ self.events['end'].synchronize()
+ #del WHHT, XHT
+ #self.MEMPOOL.free_all_blocks()
+ W_up.append(cp.cuda.get_elapsed_time(self.events['start'], self.events['end']) ) # [ms]
+ #[6]//////////////////////////// (H update)> ////////////////////////////////////////////////////////
+ self.events['start'].record()
+ #WT = self.W_d[i0:i1,:].T # [b, k].T -> WT [k, b]
+ WT = self.W_d[Q_idx].T # [b, k].T -> WT [k, b]
+ #WTW = _matmul(WT, self.W_d[i0:i1,:]) # [k, b] @ [b, k] -> WTW [k, k]
+ WTW = _matmul(WT, self.W_d[Q_idx]) # [k, b] @ [b, k] -> WTW [k, k]
+ #print(f"[!!!][iter %04d | batch %04d] shape(WTW)={WTW.shape} || shape(H_d) = {self.H_d.shape}" %(i, b))
+ #WTWH += _matmul(WTW, self.H_d) # [k, k] @ [k, n] -> WTWH [k, n]
+ WTWH.append(_matmul(WTW, self.H_d)) # [k, k] @ [k, n] -> WTWH [k, n]
+ if self.A_IS_SPARSE:
+ #WTX += self.MM(X_d, self.W_d[Q_idx], transpA=True).T # [b, n].T @ [b, k] -> [n, k].T -> [k, n]
+ WTX.append(self.MM(X_d, self.W_d[Q_idx], transpA=True).T) # [b, n].T @ [b, k] -> [n, k].T -> [k, n]
+ else:
+ #WTX += _matmul(WT, self.X_d[Q_idx]) # [k, b] @ [b, n] -> WTX [k, n]
+ WTX.append(_matmul(WT, self.X_d[Q_idx])) # [k, b] @ [b, n] -> WTX [k, n]
+ self.events['end'].record()
+ self.events['end'].synchronize()
+ #del WT, WTW
+ #self.MEMPOOL.free_all_blocks()
+ H_up.append(cp.cuda.get_elapsed_time(self.events['start'], self.events['end'])) # [ms]
+ completed_event = stream.record()
+ SQ['COMPLETED'].append(completed_event)
+ SQ['Queue'].put(st_key)
+ # Make sure all streams are done
+ for b in range(len(SQ['COMPLETED'])):
+ SQ['REDUCE'].wait_event(SQ['COMPLETED'][b])
+ cp.cuda.Device().synchronize()
+
+ # Local reduces
+ WTWH, WTX = sum(WTWH), sum(WTX)
+ #device.synchronize() # Wait for Local reduce to finish
+ #[7] GLOBAL REDUCE: Aggregate Global WTWH and WTX
+ self.events['start'].record()
+ if self.nGPUs > 1:
+ self.events['reduce_s'].record()
+ self.nccl_comm.Allreduce(WTWH, WTWH, op=nccl.NCCL_SUM, stream=None)
+ self.events['reduce_e'].record()
+ self.events['reduce_e'].synchronize()
+ dt_reduce1 += cp.cuda.get_elapsed_time(self.events['reduce_s'], self.events['reduce_e'])
+ self.events['reduce_s'].record()
+ self.nccl_comm.Allreduce(WTX, WTX, op=nccl.NCCL_SUM, stream=None)
+ self.events['reduce_e'].record()
+ self.events['reduce_e'].synchronize()
+ dt_reduce2 += cp.cuda.get_elapsed_time(self.events['reduce_s'], self.events['reduce_e'])
+ #[8] Perform H-Update
+ WTWH = WTWH + self.eps
+ WTX = _multiply(self.H_d, WTX) # [k, n] * [k, n] -> [k, n]
+ self.H_d = _divide(WTX, WTWH) # [k, n] / [k, n] -> [k, n]
+ self.events['end'].record()
+ self.events['end'].synchronize()
+ dt_h += cp.cuda.get_elapsed_time(self.events['start'], self.events['end']) # [ms]
+ #[9] Accumulate iteration performance metric
+ #print(f"[!!] [itr %04d] H2D_A = {sum(H2D_A)}" %i)
+ self.dt['W_up'] += dt_w*1.0
+ self.dt['H_up'] += dt_h*1.0
+ self.dt['allRed_WTX'] += dt_reduce1*1.0
+ self.dt['allRed_WTWH'] += dt_reduce2*1.0
+ self.dt['H2D_A'] += sum(H2D_A)
+ self.dt['H2D_W'] += sum(H2D_W)
+ self.dt['D2H_W'] += sum(D2H_W)
+ #[10]
+ #if i % 2000 == 0:
+ if False:
+ self.H_d = cp.maximum(self.H_d, self.eps)
+ self.W_d = cp.maximum(self.W_d, self.eps)
+ #self.relative_err()
+ #print(f"[+] [iter%08d] err = {self.recon_err}" %i)
+
+ #[10] Compute Average performance metric
+ self.dt['W_up'] /= float(self.params.itr)
+ self.dt['H_up'] /= float(self.params.itr)
+ self.dt['allRed_WTX'] /= float(self.params.itr)
+ self.dt['allRed_WTWH'] /= float(self.params.itr)
+ self.dt['H2D_A'] /= float(self.params.itr)
+ self.dt['H2D_W'] /= float(self.params.itr)
+ self.dt['D2H_W'] /= float(self.params.itr)
+ print(f"[!!!!] H2D(A) = {self.dt['H2D_A']} | H2D(W) = {self.dt['H2D_W']} | D2H(W) = {self.dt['D2H_W']}")
+ if False: #i % 10 == 0:
+ if self.A_IS_BATCHED:
+ self.H_h = np.maximum(self.H_h, self.eps)
+ else:
+ self.H_d = cp.maximum(self.H_d, self.eps)
+ self.W_d = cp.maximum(self.W_d, self.eps)
+ if i == self.params.itr - 1:
+ self.events['nmf_end'].record()
+ self.events['nmf_end'].synchronize()
+ self.dt['NMF'] = cp.cuda.get_elapsed_time(self.events['nmf_start'], self.events['nmf_end'])
+ #self.W_d, self.H_d = self.normalize_features(self.W_d, self.H_d)
+ #self.normalize_features()
+ if self.A_IS_SPARSE:
+ self.recon_err = 0.0
+ else:
+ self.relative_err()
+ #self.recon_err = 0.0
+ else:
+ print(f"[!] NOT OOOOOOOOK i = {i} | self.params.itr = {self.params.itr }")
+ del SQ
+
+
+
+
+[docs] def FroNMF_1D_col(self,factors=None, save_factors=False, rtrn=False):
+ r"""
+ Performs Frobenius Non-negative Matrix Factorization (NMF) using
+ 1D column-based partitioning with GPU acceleration.
+
+ Parameters:
+ - factors : (Optional) Initial values for matrix factors.
+ - save_factors : (Optional) If True, saves factorized matrices.
+ - rtrn : (Optional) If True, returns factorized matrices.
+
+ Returns:
+ - Factors (if rtrn=True).
+ """
+ print("[!!!] in NMF_ROW_BACHED")
+ self.eps = np.finfo(self.A_ij.dtype).eps
+ self.params.eps = self.eps
+ self.init_factors()
+ xht, W, X, WTW, WTX = 0.0, 0.0, 0.0, 0.0, 0.0
+ NMF_avg, H_avg, W_avg, red0_avg, red1_avg = 0.0, 0.0, 0.0, 0.0, 0.0
+ GLOBAL_NORMX = 0.0
+ self.events['nmf_start'].record()
+ for i in range(self.params.itr):
+ dt_wt,dt_wtw,dt_wtx = 0, 0, 0
+ dt_HHT, dt_WHHT, dt_xht, dt_wxht, dt_reduce0,dt_reduce1 = 0, 0, 0, 0,0, 0
+ if self.A_IS_BATCHED:
+ if self.params.W_update:
+ #//////////////////////////// (W update)> ////////////////////////////////////////////////////////
+ self.events['start'].record()
+ HT = self.H_d.T
+ HHT = _matmul(self.H_d, HT) # [k, n]@[n,k] -> HHT [k,k]
+ WTX, WTWH = [], []
+ for p in range(self.nBatch):
+ #print("[!!!!!!][FroNMF_1D_col][5] len(self.H_d) = {}".format(len(self.H_d)))
+ strm_idx = p % self.MAXSTREAMS
+ strm_name = self.stream[strm_idx]
+ stream = self.STREAMS[strm_name]
+ with stream:
+ Q_idx = strm_idx % self.batchQeueSize
+ stream.wait_event(self.FREE_STREAM[strm_name])
+ self.W_d[Q_idx].set(self.W_h[p]) # CopyH2D_Async W_h -> >W_d [bs, k]
+ WT = self.W_d[Q_idx].T # [k, bs]
+ #print("[+] shape(W_d) = {} | shape(self.W_h) = {}".format(self.W_d[Q_idx].shape, self.W_h[p].shape))
+ #print("[+] shape(H_d) = {} | shape(self.H_h) = {}".format(self.H_d.shape, HT.shape))
+ if self.A_IS_SPARSE:
+ self.X_d[Q_idx].set(self.X_per[p])
+ #print("[+] shape(X_d) = {} | shape(self.X_h) = {}".format(self.X_d[Q_idx].shape, self.X_per[p].shape))
+ self.buff_idx[Q_idx].set(self.X_idx[p])
+ self.buff_ptr[Q_idx].set(self.X_ptr[p])
+ X_d = cupyx.scipy.sparse.csr_matrix((self.X_d[Q_idx][:self.sparse_vect_size['dat'][p]],
+ self.buff_idx[Q_idx][:self.sparse_vect_size['idx'][p]],
+ self.buff_ptr[Q_idx][:self.sparse_vect_size['ptr'][p]]), shape=(self.grid_loc_m, self.batch_size))
+ XHT = self.MM(X_d, HT, transpA=False) # [bs, n] @ [n, k] -> XHT [bs, k]
+ else:
+ self.X_d[Q_idx].set(self.X_per[p]) # CopyH2D_Async
+ #print("[+] shape(X_d) = {} | shape(self.X_h) = {}".format(self.X_d[Q_idx].shape, self.X_per[p].shape))
+ XHT = self.MM(self.X_d[Q_idx], HT)
+ self.W_d[Q_idx] = _divide( _multiply(self.W_d[Q_idx], XHT), _matmul(self.W_d[Q_idx],HHT)+self.params.eps) # [bs,k] <- ([bs,k]*[bs,k]) / ([bs,k]@[k, k] + epsi)
+ self.W_d[Q_idx].get(out=self.W_h[p]) # CopyD2H_Async W_d -> W_h
+ #del HT
+ #self.MEMPOOL.free_all_blocks()
+ #//////////////////////////// (H update)> ////////////////////////////////////////////////////////
+ if self.A_IS_SPARSE:
+ wtx = _matmul( WT, X_d) # [k, bs] @ [bs, n] -> [k, n]
+ else:
+ wtx = _matmul( WT, self.X_d[Q_idx]) # [k, bs] @ [bs, n] -> [k, n]
+ #wtx = _multiply(H, wtx)
+ WTX.append(wtx)
+ wtw = _matmul(WT, self.W_d[Q_idx]) # [k, bs] @ [bs, k] -> [k, k]
+ #print("[+] shape(wtw) = {} | shape(self.H_d) = {}".format(wtw.shape, self.H_d.shape))
+ wtwh = _matmul(wtw, self.H_d) # [k, k] @ [k, n] -> [k, n]
+ WTWH.append(wtwh)
+ self.FREE_STREAM[strm_name] = stream.record()
+ # Block the `REDUCE_STREAM` until all events occur. This does not block host.
+ # This is not required when reduction is performed in the default (Stream.null)
+ # stream unless streams are created with `non_blocking=True` flag.
+ for p in range(self.batchQeueSize):
+ strm_idx = p % self.MAXSTREAMS
+ strm_name = self.stream[strm_idx]
+ self.REDUCE_STREAM.wait_event(self.FREE_STREAM[strm_name])
+ # Local reduces
+ WTWH, WTX = sum(WTWH), sum(WTX)
+ #device.synchronize() # Wait for Local reduce to finish
+ # Global reduce
+ self.nccl_comm.Allreduce(WTWH, WTWH, op=nccl.NCCL_SUM, stream=None)
+ self.sub_comm.barrier()
+ self.nccl_comm.Allreduce(WTX, WTX , op=nccl.NCCL_SUM, stream=None)
+ self.sub_comm.barrier()
+ self.H_d = _divide( _multiply(self.H_d, WTX), WTWH+self.params.eps )
+ else: # GPU MATRIX
+ for p in range(self.grid_Nc):
+ j0, j1 = p*self.grid_J, (p+1)*self.grid_J
+ #print(f"[?] shape(self.X_per[p]) = {self.X_per[p].shape}")
+ #WTX = _matmul(WT, self.X_per[:,j0:j1]) # [k, I] @ [I, J] -> [k, J]
+ if self.A_IS_SPARSE:
+ WTX = self.MM(self.X_per[p], self.W_d, transpA=True).T # [I, J].T @ [I, k] -> [J, k].T -> [k, J]
+ else:
+ WTX = self.MM(WT, self.X_per[p]) # [k, I] @ [I, J] -> [k, J]
+ WTWH = _matmul(WTW, self.H_d[:,j0:j1]) # [k, k] @ [k, J] -> [k, J]
+ WTWH += self.eps # [k, J]
+ self.events['reduce_s'].record()
+ if self.grid_Nr > 1:
+ self.nccl_comm.Allreduce(WTX, WTX, op=nccl.NCCL_SUM, stream=None)
+ self.events['reduce_e'].record()
+ self.events['reduce_e'].synchronize()
+ dt_reduce1 += cp.cuda.get_elapsed_time(self.events['reduce_s'], self.events['reduce_e'])
+ WTX = _multiply(WTX, self.H_d[:,j0:j1]) # [k, J] * [k, J] -> [k, J]
+ #WTX = _multiply(WTX, h) # [k, J] * [k, J] -> [k, J]
+ self.H_d[:, j0:j1] = WTX/WTWH # _divide(WTX, WTWH) # [k, J] / [k, J] -> [k, J]
+ #del WTWH, WTX
+ #self.MEMPOOL.free_all_blocks()
+
+ self.events['end'].record()
+ self.events['end'].synchronize()
+ del WTWH, WTX
+
+ if False: #i % 10 == 0:
+ if self.A_IS_BATCHED:
+ self.H_h = np.maximum(self.H_h, self.eps)
+ else:
+ self.H_d = cp.maximum(self.H_d, self.eps)
+ self.W_d = cp.maximum(self.W_d, self.eps)
+ if i == self.params.itr - 1:
+ self.events['nmf_end'].record()
+ self.events['nmf_end'].synchronize()
+ self.dt['NMF'] = cp.cuda.get_elapsed_time(self.events['nmf_start'], self.events['nmf_end'])
+ #self.W_d, self.H_d = self.normalize_features(self.W_d, self.H_d)
+ #self.normalize_features()
+ if self.A_IS_SPARSE:
+ NORM = cp.trace(_matmul(WTW, HHT))
+ NORM -= 2.0*cp.trace(_matmul(WT, XHT))
+ if (self.grid_Nr > 1):
+ nccl_comm.Allreduce(NORM, NORM, op=nccl.NCCL_SUM,stream=None)
+ GLOBAL_NORMX = GLOBAL_NORMX + NORM
+ self.recon_err = GLOBAL_NORMX.get()
+ else:
+ self.relative_err()
+
+
+[docs] def FroNMF_1D_row(self,factors=None, save_factors=False, rtrn=False): # [W.I.P][B.U.G]
+ r"""
+ Perform 1D Matrix Factorization using the Frobenius Norm Minimization (FroNMF).
+
+ This function updates W and H matrices for matrix factorization, where
+ the original matrix A is approximated as a product of W and H matrices.
+
+ Parameters:
+ - factors (optional): External factors for initialization (if any). Not currently utilized.
+ - save_factors (bool, optional): Flag to decide if intermediate factors should be saved.
+ - rtrn (bool, optional): Flag to decide if factors should be returned. Not currently utilized."""
+
+ ##Attributes Updated:
+ #- eps: Machine epsilon for the datatype of A_ij to avoid division by zero.
+ #- params.eps: Machine epsilon parameter.
+ #- Various time tracking events and intermediate matrices for computation.
+
+ #Notes:
+ #- This function uses CUDA operations and relies heavily on streams and events
+ # to synchronize and parallelize operations for better performance.
+ #- The function handles both batched and sparse data but has some
+ # placeholders which suggest further extensions or checks are planned (e.g., for KL updates).
+ #- The 'W_update' parameter from params decides if W should be updated in the current iteration.
+
+ #[Work in Progress - Some features are still being developed.]
+ #[Bug Warning - There might be unresolved bugs in the code.]
+
+ self.eps = np.finfo(self.A_ij.dtype).eps
+ self.params.eps = self.eps
+ self.init_factors()
+ xht, W, X, WTW, WTX = 0.0, 0.0, 0.0, 0.0, 0.0
+ NMF_avg, H_avg, W_avg, red0_avg, red1_avg = 0.0, 0.0, 0.0, 0.0, 0.0
+ #self.prune = var_init(self.params, 'prune', default=True)
+ #if self.prune: self.prune_all()
+ self.events['nmf_start'].record()
+ j0, j1 = int(self.gID*self.grid_J), int( (self.gID+1)*self.grid_J )
+ for i in range(self.params.itr):
+ dt_wt,dt_wtw,dt_wtx = 0, 0, 0
+ dt_HHT, dt_WHHT, dt_xht, dt_wxht, dt_reduce0,dt_reduce1 = 0, 0, 0, 0,0, 0
+ #//////////////////////////// (W update)> ////////////////////////////////////////////////////////
+ if self.params.W_update:
+ self.events['start'].record()
+ HT = self.H_d.T # [J, k]
+ HHT = _matmul(self.H_d, HT) # Local # [k, I] @ [I, k] -> [k, k]
+ #self.comm.barrier()
+ self.sub_comm.barrier()
+ self.events['reduce_s'].record()
+ self.nccl_comm.Allreduce(HHT, HHT, op=nccl.NCCL_SUM, stream=None)
+ self.events['reduce_e'].record()
+ self.events['reduce_e'].synchronize()
+ dt_reduce0 = cp.cuda.get_elapsed_time(self.events['reduce_s'], self.events['reduce_e'])
+ self.dt['allRed_HHT'] = dt_reduce0*1.0
+ if self.A_IS_BATCHED:
+ for p in range(self.grid_Nr):
+ i0, i1 = p*self.grid_I, (p+1)*self.grid_I
+ #print(f"[+]!!!!!!!!!! p = {p}, i0-i1 = {i0}-{i1}")
+ W = _asarray(self.W_h[i0:i1, :],dtype=self.A_ij.dtype) # [I, k], In stream 1, create event w_cp
+ X = _asarray(self.X_per[i0:i1,:],dtype=self.A_ij.dtype) # [I, J], In stream 2, create event x_cp
+ WHHT = _matmul(W, HHT) # [I, k] @ [k, k] -> [I, k], wait w_cp == True, in stream0
+ WHHT += self.eps # [I, k]
+ XHT = _matmul(X, HT) # [I, J] @ [J, k] -> [I, k], wait x_cp == True, in stream0
+ XHT = W*XHT # [I, k] * [I, k] -> [I, k]
+ self.events['reduce_s'].record()
+ if self.grid_Nc > 1:
+ self.nccl_comm.Allreduce(XHT, XHT, op=nccl.NCCL_SUM, stream=None)
+ self.events['reduce_e'].record()
+ self.events['reduce_e'].synchronize()
+ dt_reduce1 += cp.cuda.get_elapsed_time(self.events['reduce_s'], self.events['reduce_e'])
+ self.dt['allRed_XHT'] = dt_reduce1*1.0
+ W = XHT/WHHT # _divide(XHT, WHHT) # [I, k] / [I, k] -> [I, k]
+ self.W_h[i0:i1, :] = _asnumpy(W) # Local host W Update. Note W is the same everywhere
+ del XHT, WHHT
+ else:
+ for p in range(self.grid_Nr):
+ i0, i1 = p*self.grid_I, (p+1)*self.grid_I
+ WHHT = _matmul(self.W_d[i0:i1, :], HHT) # [I, k] @ [k, k] -> [I, k]
+ WHHT += self.eps # [I, k]
+ if self.A_IS_SPARSE:
+ XHT = self.MM(self.X_per[p], HT, transpA=False) # [I, J] @ [J, k] -> [I, k]
+ else:
+ #XHT = _matmul(self.X_per[i0:i1, :], HT) # [I, J] @ [J, k] -> [I, k]
+ XHT = self.MM(self.X_per[p], HT) # [I, J] @ [J, k] -> [I, k]
+ #print(f"[+] p = {p}, [i0:i1] = [{i0}:{i1}]: SHAPE(W_d) = {self.W_d[i0:i1, :].shape}, SHAPE(XHT)= {XHT.shape}")
+ XHT = self.W_d[i0:i1, :]*XHT # [I, k] * [I, k] -> [I, k]
+ self.events['reduce_s'].record()
+ if self.grid_Nc > 1:
+ self.nccl_comm.Allreduce(XHT, XHT, op=nccl.NCCL_SUM, stream=None)
+ self.events['reduce_e'].record()
+ self.events['reduce_e'].synchronize()
+ dt_reduce1 += cp.cuda.get_elapsed_time(self.events['reduce_s'], self.events['reduce_e'])
+ self.dt['allRed_XHT'] = dt_reduce1*1.0
+ self.W_d[i0:i1, :] = XHT/WHHT # _divide(XHT, WHHT) # [I, k] / [I, k] -> [I, k]
+ del XHT, WHHT
+ self.events['end'].record()
+ self.events['end'].synchronize()
+ del HT
+ dt_w = cp.cuda.get_elapsed_time(self.events['start'],self.events['end']) # [ms]
+ self.dt['W_up'] = dt_w*1.0
+ #//////////////////////////// (H update)> ////////////////////////////////////////////////////////
+ self.events['start'].record()
+ WTW = _zeros((self.k, self.k), dtype=self.A_ij.dtype) # [k, k]
+ WTX = _zeros((self.k, self.grid_J), dtype=self.A_ij.dtype) # [k, J]
+ if self.A_IS_BATCHED:
+ for p in range(self.grid_Nr):
+ i0, i1 = p*self.grid_I, (p+1)*self.grid_I
+ W = _asarray(self.W_h[i0:i1, :],dtype=self.A_ij.dtype) # [I, k], In stream 1, create event w_cp
+ X = _asarray(self.X_per[i0:i1,:],dtype=self.A_ij.dtype) # [I, J], In stream 2, create event x_cp
+ WT = W.T # [k, I], wait w_cp == True, in stream0
+ WTW += _matmul(WT, W) # [k, I] @ [I, k] -> [k, k], in stream 0
+ WTX += _matmul(WT, X) # [k, I] @ [I, J] -> [k, J], wait x_cp == True, in stream0
+ del WT , W, X
+ else:
+ for p in range(self.grid_Nr):
+ i0, i1 = p*self.grid_I, (p+1)*self.grid_I
+ #W = self.W_d[i0:i1, :] # [I, k]
+ WT = self.W_d[i0:i1, :].T # [k, I]
+ WTW += _matmul(WT, self.W_d[i0:i1, :]) # [k, I] @ [I, k] -> [k, k]
+ if self.A_IS_SPARSE:
+ WTX += self.MM(self.X_per[p], self.W_d[i0:i1,:], transpA=True).T # [I, J].T @ [I, k] -> [J, k].T -> [k, J]
+ else:
+ #WTX += self.MM(WT, self.X_per[i0:i1,:]) # [k, I] @ [I, J] -> [k, J]
+ WTX += self.MM(WT, self.X_per[p])
+ del WT #,W
+ WTWH = _matmul(WTW,self.H_d) # [k, k] @ [k, J] -> [k, J]
+ del WTW
+ WTWH += self.eps
+ WTX = _multiply(WTX,self.H_d) # [k, J] * [k, J] -> [k, J]
+ #H = _divide(WTX, WTW) # [k, J] / [k, J] -> [k, J]
+ self.H_d = _divide(WTX, WTWH) # [k, J] / [k, J] -> [k, J]
+ self.events['end'].record()
+ self.events['end'].synchronize()
+ del WTX #, WTW
+ dt_h = cp.cuda.get_elapsed_time(self.events['start'], self.events['end']) # [ms]
+ self.dt['H_up'] = dt_h*1.0
+ if self.A_IS_BATCHED:
+ noSupportFor("KL update for LARGE A_ij")
+ else:
+ WTXWH = _zeros((self.k, self.grid_J), dtype=self.A_ij.dtype) # WTXWH [k, J]
+ for p in range(self.grid_Nr):
+ i0, i1 = p*self.grid_I, (p+1)*self.grid_I
+ XWH = _matmul(self.W_d[i0:i1, :], self.H_d) # XWH = W@H : [I, k] @ [k, J] -> [I, J]
+ if self.A_IS_SPARSE:
+ noSupportFor(" SPARSE A_ij")
+ else:
+ XWH = self.X_per[p]/(XWH + self.eps) # XWH = X/(XWH+e) : [I, J] / [I, J] + eps -> [I, J]
+ WTXWH += _matmul(self.W_d[i0:i1, :].T, XWH) # WTXWH = W.T@XWH@HT : [k, I] @ [I, J] -> [k, J]
+ self.H_d = self.H_d*WTXWH/X1[:, None]
+ del WTXWH,XWH
+
+ self.relative_err()
+
+
+# @author: Ismael Boureima
+import cupy, cupyx, numpy
+import scipy
+from .toolz import amber, blue, green, red
+from scipy.sparse import coo_matrix, csc_matrix, csr_matrix, dia_matrix, issparse
+
+[docs]def spMat(A, Format='coo', shape=None, dtype=None, copy=False):
+ r"""
+ Converts a given matrix A into a specified sparse format using CuPy.
+
+ Parameters:
+ -----------
+ A : ndarray or sparse matrix
+ Input matrix to be converted.
+ Format : {'coo', 'csr', 'csc', 'dia'}
+ Desired sparse format. Default is 'coo'.
+ shape : tuple, optional
+ Desired shape for the sparse matrix.
+ dtype : data-type, optional
+ Data type of the result.
+ copy : bool, optional
+ If true, it guarantees the input data A is not modified. Default is False.
+
+ Returns:
+ --------
+ Sparse matrix in the desired format and library (cupy).
+ """
+ fmt = Format.lower()
+ assert fmt in ['coo', 'csr', 'csc', 'dia'], f"[{red('!')}] Format must be in [{green('coo')}, {green('crs')}, {green('csc')}, {green('dia')}] "
+ isNdArr = False
+ if type(A) in [numpy.ndarray]: isNdArr=True
+ if fmt in ['coo']:
+ if isNdArr:
+ return cupyx.scipy.sparse.coo_matrix(coo_matrix(A), shape=shape, dtype=dtype, copy=copy)
+ else:
+ return cupyx.scipy.sparse.coo_matrix(A, shape=shape, dtype=dtype, copy=copy)
+ elif fmt in ['csc']:
+ if isNdArr:
+ return cupyx.scipy.sparse.csc_matrix(coo_matrix(A), shape=shape, dtype=dtype, copy=copy)
+ else:
+ return cupyx.scipy.sparse.csc_matrix(A, shape=shape, dtype=dtype, copy=copy)
+ elif fmt in ['csr']:
+ if isNdArr:
+ return cupyx.scipy.sparse.csr_matrix(coo_matrix(A), shape=shape, dtype=dtype, copy=copy)
+ else:
+ return cupyx.scipy.sparse.csr_matrix(A, shape=shape, dtype=dtype, copy=copy)
+ elif fmt in ['dia']:
+ if isNdArr:
+ return cupyx.scipy.sparse.dia_matrix(coo_matrix(A), shape=shape, dtype=dtype, copy=copy)
+ else:
+ return cupyx.scipy.sparse.dia_matrix(A, shape=shape, dtype=dtype, copy=copy)
+ else:
+ print(f"[!] sparse matrix format '{fmt}' is not supported")
+
+
+
+[docs]def spCu2Sc(A, copy=False):
+ r"""
+ Converts a CuPy sparse matrix to a SciPy sparse matrix.
+
+ Parameters:
+ -----------
+ A : cupyx sparse matrix
+ Input matrix to be converted.
+ copy : bool, optional
+ If true, it guarantees the input data A is not modified. Default is False.
+
+ Returns:
+ --------
+ Sparse matrix in the desired format and library (scipy).
+ """
+ fmt = A.format.lower()
+ assert fmt in ['coo', 'csr', 'csc', 'dia'], f"[{red('!')}] Format must be in [{green('coo')}, {green('crs')}, {green('csc')}, {green('dia')}] "
+ if fmt in ['coo']:
+ return scipy.sparse.coo_matrix((A.data.get(),A.indices.get(),A.indptr.get()),shape=A.shape,dtype=A.dtype)
+ elif fmt in ['csc']:
+ return scipy.sparse.csc_matrix((A.data.get(),A.indices.get(),A.indptr.get()),shape=A.shape,dtype=A.dtype)
+ elif fmt in ['csr']:
+ return scipy.sparse.csr_matrix((A.data.get(),A.indices.get(),A.indptr.get()),shape=A.shape,dtype=A.dtype)
+ elif fmt in ['dia']:
+ return scipy.sparse.dia_matrix((A.data.get(),A.indices.get(),A.indptr.get()),shape=A.shape,dtype=A.dtype)
+ else:
+ print(f"[!] sparse matrix format '{fmt}' is not supported")
+
+[docs]def spSc2Cu(A, copy=False):
+ r"""
+ Converts a SciPy sparse matrix to a CuPy sparse matrix.
+
+ Parameters:
+ -----------
+ A : scipy sparse matrix
+ Input matrix to be converted.
+ copy : bool, optional
+ If true, it guarantees the input data A is not modified. Default is False.
+
+ Returns:
+ --------
+ Sparse matrix in the desired format and library (cupy).
+ """
+ fmt = A.format.lower()
+ assert fmt in ['coo', 'csr', 'csc', 'dia'], f"[{red('!')}] Format must be in [{green('coo')}, {green('crs')}, {green('csc')}, {green('dia')}] "
+ if fmt in ['coo']:
+ return cupyx.scipy.sparse.coo_matrix((cp.array(A.data), cp.array(A.indices), cp.array(A.indptr)), shape=A.shape, dtype=A.dtype)
+ elif fmt in ['csc']:
+ return cupyx.scipy.sparse.csc_matrix((cp.array(A.data), cp.array(A.indices), cp.array(A.indptr)), shape=A.shape, dtype=A.dtype)
+ elif fmt in ['csr']:
+ return cupyx.scipy.sparse.csr_matrix((cp.array(A.data), cp.array(A.indices), cp.array(A.indptr)), shape=A.shape, dtype=A.dtype)
+ elif fmt in ['dia']:
+ return cupyx.scipy.sparse.dia_matrix((cp.array(A.data), cp.array(A.indices), cp.array(A.indptr)), shape=A.shape, dtype=A.dtype)
+ else:
+ print(f"[!] sparse matrix format '{fmt}' is not supported")
+
+
+[docs]def spRandMat(nr, nc, density, Format='coo', dtype=cupy.float32):
+ r"""
+ Generates a random sparse matrix using CuPy.
+
+ Parameters:
+ -----------
+ nr : int
+ Number of rows.
+ nc : int
+ Number of columns.
+ density : float
+ Desired density for the sparse matrix. Values between 0 and 1.
+ Format : {'coo', 'csr', 'csc', 'dia'}, optional
+ Desired sparse format. Default is 'coo'.
+ dtype : data-type, optional
+ Data type of the result. Default is cupy.float32.
+
+ Returns:
+ --------
+ Random sparse matrix in the desired format and library (cupy).
+ """
+ fmt = Format.lower()
+ return cupyx.scipy.sparse.random(nr, nc, density=density, format=fmt, dtype=dtype)
+
+
+[docs]def spMM(A, B, alpha=1.0, beta=0.0, transpA=False, transpB=False):
+ r"""
+ Performs matrix multiplication of sparse matrix A and dense matrix B using CuPy.
+
+ Parameters:
+ -----------
+ A : cupyx sparse matrix
+ Left sparse matrix.
+ B : ndarray
+ Right dense matrix.
+ alpha : float, optional
+ Scalar multiplier for the product of A and B. Default is 1.0.
+ beta : float, optional
+ Scalar multiplier for the initial matrix C (if provided). Default is 0.0.
+ transpA : bool, optional
+ If True, transpose matrix A before multiplication. Default is False.
+ transpB : bool, optional
+ If True, transpose matrix B before multiplication. Default is False.
+
+ Returns:
+ --------
+ Resultant matrix after multiplication.
+ """
+ assert cupy.cusparse.check_availability('spmm'), "[!] spmm is not available"
+ #if transpA : A = A.T
+ #if transpB : B = B.T
+ if not A.has_canonical_format: A.sum_duplicates()
+ B = cupy.array(B, order='f')
+ return cupy.cusparse.spmm( A, B, alpha=alpha, beta=beta, transa=transpA, transb=transpB)
+
+
+[docs]def spMM_with_C(A, B, C , alpha=1.0, beta=0.0, transpA=False, transpB=False):
+ r"""
+ Performs matrix multiplication of sparse matrix A and dense matrix B, and adds it to matrix C using CuPy.
+
+ Parameters:
+ -----------
+ A : cupyx sparse matrix
+ Left sparse matrix.
+ B : ndarray
+ Right dense matrix.
+ C : ndarray
+ Dense matrix to which the result is added.
+ alpha : float, optional
+ Scalar multiplier for the product of A and B. Default is 1.0.
+ beta : float, optional
+ Scalar multiplier for the initial matrix C. Default is 0.0.
+ transpA : bool, optional
+ If True, transpose matrix A before multiplication. Default is False.
+ transpB : bool, optional
+ If True, transpose matrix B before multiplication. Default is False.
+
+ Returns:
+ --------
+ Resultant matrix after multiplication and addition to C.
+ """
+ assert cupy.cusparse.check_availability('spmm'), "[!] spmm is not available"
+ #if transpA : A = A.T
+ #if transpB : B = B.T
+ if not A.has_canonical_format: A.sum_duplicates()
+ B = cupy.array(B, order='f')
+ C = cupy.array(C, order='f')
+ return cupy.cusparse.spmm( A, B, c=C, alpha=alpha, beta=beta, transa=transpA, transb=transpB)
+
+
+# @author: Ismael Boureima
+
+TPB = 32
+
+import cupy as cp
+import numpy as np
+#from sgemm import sgemm
+from .toolz import amber, blue, green, purple, red
+
+# Initializing a pinned memory pool in CuPy to speed up host-device memory transfers.
+# However, the following two lines are commented out and thus inactive.
+#pinned_memory_pool = cupy.cuda.PinnedMemoryPool()
+#cupy.cuda.set_pinned_memory_allocator(pinned_memory_pool.malloc)
+#mem = cp.cuda.alloc_pinned_memory(array.nbytes)
+
+[docs]def pin_memory(array):
+ """Allocate memory in the pinned (or "page-locked") area. Data in this memory can be transferred to the GPU faster."""
+ mem = cp.cuda.alloc_pinned_memory(array.nbytes)
+ ret = np.frombuffer(mem, array.dtype, array.size).reshape(array.shape)
+ ret[...] = array
+ return ret
+
+[docs]def read_code(code_filename, params):
+ """Read a code file and prepend it with CUDA kernel parameter definitions."""
+ with open(code_filename, 'r') as f:
+ code = f.read()
+ for k, v in params.items():
+ code = '#define ' + k + ' ' + str(v) + '\n' + code
+ return code
+
+
+[docs]def benchmark(func, args, n_run=1):
+ """Benchmark a given function with provided arguments over n_run runs. Return a list of execution times."""
+ times = []
+ for _ in range(n_run):
+ start = cp.cuda.Event()
+ end = cp.cuda.Event()
+ start.record()
+ func(*args)
+ end.record()
+ end.synchronize()
+ times.append(cp.cuda.get_elapsed_time(start, end)) # milliseconds
+ return times
+
+[docs]def timeExecution(func):
+ """Decorator to time the execution of a function and print it."""
+ def wrapper(*args, **kwargs):
+ start, end = cp.cuda.Event(), cp.cuda.Event()
+ start.record()
+ out = func(*args, **kwargs)
+ end.record()
+ end.synchronize()
+ t = cp.cuda.get_elapsed_time(start, end) # milliseconds
+ print(f"[INFO]: {blue(func.__name__)} ran in {red(round(t,4))} [ms]")
+ return out
+ return wrapper
+
+@timeExecution
+def _asarray(x):
+ """Convert the input to a CuPy array."""
+ return cp.asarray(x, dtype=x.dtype)
+
+@timeExecution
+def _zeros(*args, **kwargs):
+ """Return a new CuPy array of given shape and type, filled with zeros."""
+ return cp.zeros(*args, **kwargs)
+
+@timeExecution
+def _asnumpy(x):
+ """Convert a CuPy array to a NumPy array."""
+ return cp.asnumpy(x)
+
+@timeExecution
+def _divide(a,b):
+ """Divide two CuPy arrays element-wise."""
+ return cp.divide(a,b, out=None)
+
+@timeExecution
+def _matmul(a,b):
+ """Matrix multiplication of two CuPy arrays."""
+ return cp.matmul(a,b, out=None)
+
+@timeExecution
+def _multiply(a,b):
+ """Multiply two CuPy arrays element-wise."""
+ return cp.multiply(a,b)
+
+
+
+
+
# @Author: Manish Bhattarai, Erik Skau
import argparse
@@ -160,7 +235,7 @@ Source code for pyDNMFk.data_generator
[docs]def parser():
- r"""
+ r"""
Reads the input arguments from the user and parses the parameters to the data generator module.
"""
parser = argparse.ArgumentParser(description='Data generator arguments')
@@ -175,7 +250,7 @@ Source code for pyDNMFk.data_generator
[docs]class data_generator():
- r"""
+ r"""
Generates synthetic data in distributed manner where each MPI process generates a chunk from the data parallelly.
The W matrix is generated with gaussian distribution whereas the H matrix is random.
@@ -211,7 +286,7 @@ Source code for pyDNMFk.data_generator
# self.factor = k
[docs] def gauss_matrix_generator(self, n, k):
- r"""
+ r"""
Construct a matrix of dimensions n by k where the ith column is a Gaussian kernel corresponding to approximately N(i*n/k, 0.01*n^2)
Parameters
@@ -230,28 +305,28 @@ Source code for pyDNMFk.data_generator
offset = n / k / 2 - 0.5
noverk = n / k
- coeff = -k / (.01 * n ** 2)
+ coeff = -k / (.01 * n ** 2)
return lambda i, j: np.exp(coeff * (i - (j * noverk + offset)) ** 2)
[docs] def determine_block_index_range_asymm(self):
- '''Determines the start and end indices for the Data block for each rank'''
+ '''Determines the start and end indices for the Data block for each rank'''
chunk_ind = np.unravel_index(self.rank, self.pgrid)
start_inds = [i * (n // k) + min(i, n % k) for n, k, i in zip(self.shape, self.pgrid, chunk_ind)]
end_inds = [(i + 1) * (n // k) + min((i + 1), n % k) - 1 for n, k, i in zip(self.shape, self.pgrid, chunk_ind)]
return start_inds, end_inds
[docs] def determine_block_shape_asymm(self):
- '''Determines the shape for the Data block for each rank'''
+ '''Determines the shape for the Data block for each rank'''
start_inds, end_inds = self.determine_block_index_range_asymm()
return [(j - i + 1) for (i, j) in zip(start_inds, end_inds)]
[docs] def random_matrix_generator(self, n, k, seed):
- '''Generator for random matric with given seed'''
+ '''Generator for random matric with given seed'''
np.random.seed(seed)
return np.random.rand(n, k)
[docs] def dist_fromfunction(self, func, shape, pgrid, *args, unravel_index=np.unravel_index, **kwargs):
- """
+ """
produces X_{i,j} = func(i,j) in a distributed manner, so that each processor has an array_split section of X according to the grid.
"""
grid_index = unravel_index()
@@ -260,7 +335,7 @@ Source code for pyDNMFk.data_generator
return np.fromfunction(lambda *x: func(*[a + b for a, b in zip(x, start_index)]), block_shape, *args, **kwargs)
[docs] def unravel_column(self):
- '''finds the column rank for 2d grid'''
+ '''finds the column rank for 2d grid'''
def wrapper(*args, **kwargs):
row, col = np.unravel_index(self.rank, self.pgrid)
@@ -269,29 +344,29 @@ Source code for pyDNMFk.data_generator
return wrapper
[docs] def unravel_row(self): # ,ind, shape):
- '''finds the row rank for 2d grid'''
+ '''finds the row rank for 2d grid'''
row, col = np.unravel_index(self.rank, self.pgrid)
return (row // self.pgrid[0], col)
[docs] def create_folder_dir(self, fpath):
- '''Create a folder if doesn't exist'''
+ '''Create a folder if doesn't exist'''
try:
os.mkdir(fpath)
except:
pass
[docs] def generate_factors_data(self):
- """Generates the chunk of factors W,H and data X for each MPI process"""
+ """Generates the chunk of factors W,H and data X for each MPI process"""
W_gen = self.dist_fromfunction(self.gauss_matrix_generator(self.m, self.k), (self.m, self.k), (self.p_r, 1),
- unravel_index=self.unravel_column())
+ unravel_index=self.unravel_column()).astype(np.float32)
H_gen = self.random_matrix_generator(self.k, self.determine_block_shape_asymm()[1],
- self.unravel_row()[1])
+ self.unravel_row()[1]).astype(np.float32)
X_gen = W_gen @ H_gen
print('For rank=', self.rank, ' dimensions of W,H and X are ', W_gen.shape, H_gen.shape, X_gen.shape)
return W_gen, H_gen, X_gen
[docs] def fit(self):
- '''generates and save factors'''
+ '''generates and save factors'''
W_gen, H_gen, X_gen = self.generate_factors_data()
self.create_folder_dir(self.fpath)
self.create_folder_dir(self.fpath + 'W_factors')
@@ -311,47 +386,73 @@ Source code for pyDNMFk.data_generator
data_gen.fit()
-
-
-
-
+
+
+
+
-
-
-
+
+
+
-
-
+
+
+
-
+
+
\ No newline at end of file
diff --git a/docs/_modules/pyDNMFk/data_io.html b/docs/_modules/pyDNMFk/data_io.html
index 7ae59f4..5e0166d 100644
--- a/docs/_modules/pyDNMFk/data_io.html
+++ b/docs/_modules/pyDNMFk/data_io.html
@@ -1,160 +1,235 @@
-
-
-
+
+
+
+
+
+
+
+ pyDNMFk.data_io — pyDNMFk 0.0.1 documentation
-
- pyDNMFk.data_io — pyDNMFk 1.0.0 documentation
-
+
-
-
+
+
+
+
+
+
+
+
+
+
+
+
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-
+ Skip to main content
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
-
+
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
+
-
-
- - »
-
- - Module code »
-
- - pyDNMFk.data_io
-
-
- -
-
-
-
-
+
+
+
+
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Source code for pyDNMFk.data_io
# @author: Manish Bhattarai
import glob
import os
-
+import scipy
import h5py
import pandas as pd
from scipy.io import loadmat
@@ -163,7 +238,7 @@ Source code for pyDNMFk.data_io
[docs]class data_read():
- r"""Class for reading data.
+ r"""Class for reading data.
Parameters
----------
@@ -189,6 +264,7 @@ Source code for pyDNMFk.data_io
self.fname = args.fname
self.comm = args.comm1
self.rank = self.comm.rank
+ self.precision = args.precision if args.precision else 'float32'
self.data = 0
if self.ftype == 'folder':
self.file_path = self.fpath + self.fname + str(self.comm.rank) + '.npy'
@@ -197,28 +273,34 @@ Source code for pyDNMFk.data_io
[docs] @comm_timing()
def read(self):
- r"""Data read function"""
+ r"""Data read function"""
return self.read_dat()
[docs] @comm_timing()
def read_file_npy(self):
- r"""Numpy data read function"""
+ r"""Numpy data read function"""
self.data = np.load(self.file_path)
[docs] @comm_timing()
def read_file_csv(self):
- r"""CSV data read function"""
+ r"""CSV data read function"""
self.data = pd.read_csv(self.file_path, header=None).values
[docs] @comm_timing()
def read_file_mat(self):
- r"""mat file read function"""
+ r"""mat file read function"""
self.data = loadmat(self.file_path)['X']
+[docs] @comm_timing()
+ def read_file_npz(self):
+ r"""mat file read function"""
+ self.data = scipy.sparse.load_npz(self.file_path)
+ #self.data = loadmat(self.file_path)['X']
+
[docs] @comm_timing()
def data_partition(
self):
- r"""
+ r"""
This function divides the input matrix into chunks as specified by grid configuration.
Return n array of shape (nrows_i, ncols_i) where i is the index of each chunk.
@@ -233,16 +315,19 @@ Source code for pyDNMFk.data_io
[docs] @comm_timing()
def save_data_to_file(self, fpath):
- r"""This function saves the splitted data to numpy array indexed with chunk number"""
+ r"""This function saves the splitted data to numpy array indexed with chunk number"""
fname = fpath + 'A_' + self.comm.rank + '.npy'
np.save(fname, self.data)
[docs] @comm_timing()
def read_dat(self):
- r"""Function for reading the data and split into chunks to be reach by each MPI rank"""
+ r"""Function for reading the data and split into chunks to be reach by each MPI rank"""
if self.ftype == 'npy':
self.read_file_npy()
self.data_partition()
+ elif self.ftype == "npz":
+ self.read_file_npz()
+ self.data_partition()
elif self.ftype == 'csv' or self.ftype == 'txt':
self.read_file_csv()
self.data_partition()
@@ -251,11 +336,11 @@ Source code for pyDNMFk.data_io
self.data_partition()
if self.ftype == 'folder':
self.read_file_npy()
- return self.data
+ return self.data.astype(self.precision)
[docs]class split_files_save():
- r"""Rank 0 based data read, split and save"""
+ r"""Rank 0 based data read, split and save"""
@comm_timing()
def __init__(self, data, pgrid, fpath):
@@ -267,14 +352,14 @@ Source code for pyDNMFk.data_io
[docs] @comm_timing()
def split_files(self):
- r"""Compute the index range for each block and partition the data as per the chunk"""
+ r"""Compute the index range for each block and partition the data as per the chunk"""
dtr_blk_idx = [determine_block_params(rank, self.pgrid, self.data.shape).determine_block_index_range_asymm() for
rank in range(np.product(self.pgrid))]
self.split = [self.data[i[0][0]:i[1][0] + 1, i[0][1]:i[1][1] + 1] for i in dtr_blk_idx]
[docs] @comm_timing()
def save_data_to_file(self):
- r"""Function to save the chunks into numpy files"""
+ r"""Function to save the chunks into numpy files"""
s = 0
self.split = self.split_files()
for i in range(self.p_r * self.p_c):
@@ -286,7 +371,7 @@ Source code for pyDNMFk.data_io
[docs]class data_write():
- r"""Class for writing data/results.
+ r"""Class for writing data/results.
Parameters
----------
@@ -310,7 +395,7 @@ Source code for pyDNMFk.data_io
[docs] @comm_timing()
def create_folder_dir(self, fpath):
- r"""Create directory if not present"""
+ r"""Create directory if not present"""
try:
os.mkdir(fpath)
except:
@@ -318,7 +403,7 @@ Source code for pyDNMFk.data_io
[docs] @comm_timing()
def save_factors(self, factors, reg=False):
- r"""Save the W and H factors for each MPI process"""
+ r"""Save the W and H factors for each MPI process"""
self.create_folder_dir(self.fpath)
if reg == True:
W_factors_pth = self.fpath + 'W_reg_factors/'
@@ -342,7 +427,7 @@ Source code for pyDNMFk.data_io
[docs] @comm_timing()
def save_cluster_results(self, params):
- r"""Save cluster results to a h5 file with rank 0"""
+ r"""Save cluster results to a h5 file with rank 0"""
if self.rank == 0:
with h5py.File(self.fpath + 'results.h5', 'w') as hf:
hf.create_dataset('clusterSilhouetteCoefficients', data=params['clusterSilhouetteCoefficients'])
@@ -350,11 +435,12 @@ Source code for pyDNMFk.data_io
hf.create_dataset('L_err', data=params['L_err'])
hf.create_dataset('L_errDist', data=params['L_errDist'])
hf.create_dataset('avgErr', data=params['avgErr'])
- hf.create_dataset('ErrTol', data=params['recon_err'])
+ hf.create_dataset('ErrTol', data=params['recon_err'])
+ hf.create_dataset('AIC', data=params['AIC'])
[docs]class read_factors():
- r"""Class for reading saved factors.
+ r"""Class for reading saved factors.
Parameters
----------
@@ -367,24 +453,24 @@ Source code for pyDNMFk.data_io
@comm_timing()
def __init__(self, factors_path, pgrid):
self.factors_path = factors_path
- self.W_path = self.factors_path + 'W/*'
- self.H_path = self.factors_path + 'H/*'
+ self.W_path = self.factors_path + 'W_reg_factors/*'
+ self.H_path = self.factors_path + 'H_reg_factors/*'
self.p_grid = pgrid
self.load_factors()
[docs] @comm_timing()
def custom_read_npy(self, fpath):
- r"""Read numpy files"""
+ r"""Read numpy files"""
data = np.load(fpath)
return data
[docs] @comm_timing()
def read_factor(self, fpath):
- """Read factors as chunks and stack them"""
+ """Read factors as chunks and stack them"""
files = glob.glob(fpath)
data = []
if len(files) == 1:
- data = self.custom_read_npy(files)
+ data = self.custom_read_npy(files[0])
else:
for file in np.sort(files):
data.append(self.custom_read_npy(file))
@@ -392,7 +478,7 @@ Source code for pyDNMFk.data_io
[docs] @comm_timing()
def load_factors(self):
- r"""Load the final stacked factors for visualization"""
+ r"""Load the final stacked factors for visualization"""
W_data, ct_W = self.read_factor(self.W_path)
H_data, ct_H = self.read_factor(self.H_path)
if ct_W > 1: W_data = np.vstack((W_data))
@@ -405,47 +491,73 @@ Source code for pyDNMFk.data_io
return W_data, H_data
-
-
-
-
+
+
+
+
-
-
-
+
+
+
-
-
+
+
+
-
+
+
\ No newline at end of file
diff --git a/docs/_modules/pyDNMFk/dist_clustering.html b/docs/_modules/pyDNMFk/dist_clustering.html
index 401f520..bb71af9 100644
--- a/docs/_modules/pyDNMFk/dist_clustering.html
+++ b/docs/_modules/pyDNMFk/dist_clustering.html
@@ -1,162 +1,237 @@
-
-
-
+
+
+
+
+
+
+
+ pyDNMFk.dist_clustering — pyDNMFk 0.0.1 documentation
-
- pyDNMFk.dist_clustering — pyDNMFk 1.0.0 documentation
-
+
-
-
+
+
+
+
+
+
+
+
+
+
+
+
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-
+ Skip to main content
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
-
+
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-
-
-
- - »
-
- - Module code »
-
- - pyDNMFk.dist_clustering
-
-
- -
-
-
-
-
+
+
+
+
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Source code for pyDNMFk.dist_clustering
# @Author: Manish Bhattarai, Erik Skau
from .utils import *
[docs]class custom_clustering():
- r"""
+ r"""
Greedy algorithm to approximate a quadratic assignment problem to cluster vectors. Given p groups of k vectors, construct k clusters, each cluster containing a single vector from each of the p groups. This clustering approximation uses cos distances and mean centroids.
Parameters
@@ -182,7 +257,7 @@ Source code for pyDNMFk.dist_clustering
[docs] @comm_timing()
def normalize_by_W(self):
- r'''Normalize the factors W and H'''
+ r'''Normalize the factors W and H'''
Wall_norm = (self.W_all * self.W_all).sum(axis=0)
if self.p_r != 1:
Wall_norm = self.comm1.allreduce(Wall_norm)
@@ -193,7 +268,7 @@ Source code for pyDNMFk.dist_clustering
[docs] @comm_timing()
def mad(self, data, flag=1, axis=-1):
- r'''Compute the median/mean absolute deviation'''
+ r'''Compute the median/mean absolute deviation'''
if flag == 1: # the median absolute deviation
return np.nanmedian(np.absolute(data - np.nanmedian(data, axis=axis, keepdims=True)), axis=axis)
else: # flag = 0 the mean absolute deviation
@@ -202,7 +277,7 @@ Source code for pyDNMFk.dist_clustering
[docs] @comm_timing()
def change_order(self, tens):
- r'''change the order of features'''
+ r'''change the order of features'''
ans = list(range(len(tens)))
for p in tens:
ans[p[0]] = p[1]
@@ -210,7 +285,7 @@ Source code for pyDNMFk.dist_clustering
[docs] @comm_timing()
def greedy_lsa(self, A):
- r"""Return the permutation order"""
+ r"""Return the permutation order"""
X = A.copy()
pairs = []
for i in range(X.shape[0]):
@@ -223,7 +298,7 @@ Source code for pyDNMFk.dist_clustering
[docs] @comm_timing()
def dist_feature_ordering(self, centroids, W_sub):
- r'''return the features in proper order'''
+ r'''return the features in proper order'''
k = W_sub.shape[1]
dist = centroids.T @ W_sub
if self.p_r != 1:
@@ -235,7 +310,7 @@ Source code for pyDNMFk.dist_clustering
[docs] @comm_timing()
def dist_custom_clustering(self, centroids=None, vb=0):
- """
+ """
Performs the distributed custom clustering
Parameters
@@ -281,7 +356,7 @@ Source code for pyDNMFk.dist_clustering
[docs] @comm_timing()
def dist_silhouettes(self):
- """
+ """
Computes the cosine distances silhouettes of a distributed clustering of vectors.
Returns
@@ -314,7 +389,7 @@ Source code for pyDNMFk.dist_clustering
[docs] @comm_timing()
def fit(self):
- r"""
+ r"""
Calls the sub routines to perform distributed custom clustering and compute silhouettes
Returns
@@ -341,47 +416,73 @@ Source code for pyDNMFk.dist_clustering
return result
-
-
-
-
+
+
+
+
-
-
-
+
+
+
-
-
+
+
+
-
+
+
\ No newline at end of file
diff --git a/docs/_modules/pyDNMFk/dist_comm.html b/docs/_modules/pyDNMFk/dist_comm.html
index ead97cd..8f512f1 100644
--- a/docs/_modules/pyDNMFk/dist_comm.html
+++ b/docs/_modules/pyDNMFk/dist_comm.html
@@ -1,159 +1,234 @@
-
-
-
+
+
+
+
+
+
+
+ pyDNMFk.dist_comm — pyDNMFk 0.0.1 documentation
-
- pyDNMFk.dist_comm — pyDNMFk 1.0.0 documentation
-
+
-
-
+
+
+
+
+
+
+
+
+
+
+
+
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-
+ Skip to main content
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
-
+
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-
-
-
- - »
-
- - Module code »
-
- - pyDNMFk.dist_comm
-
-
- -
-
-
-
-
+
+
+
+
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Source code for pyDNMFk.dist_comm
# @author: Manish Bhattarai
[docs]class MPI_comm():
- """Initialization of MPI communicator to construct the cartesian topology and sub communicators
+ """Initialization of MPI communicator to construct the cartesian topology and sub communicators
Parameters
----------
@@ -176,7 +251,7 @@ Source code for pyDNMFk.dist_comm
self.coord2d = self.cartesian2d.Get_coords(self.rank)
[docs] def cart_1d_row(self):
- """
+ """
Constructs a cartesian row communicator through construction of a sub communicator across rows
Returns
@@ -190,7 +265,7 @@ Source code for pyDNMFk.dist_comm
return self.cartesian1d_row
[docs] def cart_1d_column(self):
- """
+ """
Constructs a cartesian column communicator through construction of a sub communicator across columns
Returns
@@ -204,52 +279,78 @@ Source code for pyDNMFk.dist_comm
return self.cartesian1d_column
[docs] def Free(self):
- """ Frees the sub communicators"""
+ """ Frees the sub communicators"""
self.cart_1d_row().Free()
self.cart_1d_column().Free()
-
-
-
-
+
+
+
+
-
-
-
+
+
+
-
-
+
+
+
-
+
+
\ No newline at end of file
diff --git a/docs/_modules/pyDNMFk/dist_nmf.html b/docs/_modules/pyDNMFk/dist_nmf.html
index 8adb823..3f941f4 100644
--- a/docs/_modules/pyDNMFk/dist_nmf.html
+++ b/docs/_modules/pyDNMFk/dist_nmf.html
@@ -1,155 +1,230 @@
-
-
-
+
+
+
+
+
+
+
+ pyDNMFk.dist_nmf — pyDNMFk 0.0.1 documentation
-
- pyDNMFk.dist_nmf — pyDNMFk 1.0.0 documentation
-
+
-
-
+
+
+
+
+
+
+
+
+
+
+
+
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-
+ Skip to main content
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
-
+
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
+
-
-
- - »
-
- - Module code »
-
- - pyDNMFk.dist_nmf
-
-
- -
-
-
-
-
+
+
+
+
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Source code for pyDNMFk.dist_nmf
# @author: Manish Bhattarai
from numpy import matlib
@@ -158,7 +233,7 @@ Source code for pyDNMFk.dist_nmf
[docs]class nmf_algorithms_2D():
- """
+ """
Performs the distributed NMF operation along 2D cartesian grid
Parameters
@@ -217,7 +292,7 @@ Source code for pyDNMFk.dist_nmf
self.local_H_n = self.H_ij.shape[1]
[docs] def update(self):
- """Performs 1 step Update for factors W and H based on NMF method and corresponding norm minimization
+ """Performs 1 step Update for factors W and H based on NMF method and corresponding norm minimization
Returns
-------
@@ -247,7 +322,7 @@ Source code for pyDNMFk.dist_nmf
[docs] @comm_timing()
def global_gram(self, A):
- r""" Distributed gram computation
+ r""" Distributed gram computation
Computes the global gram operation of matrix A
.. math:: A^TA
@@ -271,7 +346,7 @@ Source code for pyDNMFk.dist_nmf
[docs] @comm_timing()
def global_mm(self, A, B):
- r""" Distributed matrix multiplication
+ r""" Distributed matrix multiplication
Computes the global matrix multiplication of matrix A and B
.. math:: AB
@@ -292,12 +367,12 @@ Source code for pyDNMFk.dist_nmf
self.comm1.barrier()
return AB_glob
- '''Functions for Fro MU NMF update'''
+ '''Functions for Fro MU NMF update'''
[docs] @comm_timing()
def ATW_glob(self):
- r""" Distributed computation of W^TA
+ r""" Distributed computation of W^TA
Computes the global matrix multiplication of matrix W and A
.. math:: W^TA
@@ -327,7 +402,7 @@ Source code for pyDNMFk.dist_nmf
[docs] @comm_timing()
def AH_glob(self, H_ij=None):
- r""" Distributed computation of AH^T
+ r""" Distributed computation of AH^T
Computes the global matrix multiplication of matrix A and H
.. math:: AH^T
@@ -359,7 +434,7 @@ Source code for pyDNMFk.dist_nmf
[docs] def Fro_MU_update_H(self):
- r"""
+ r"""
Frobenius norm based multiplicative update of H parameter
Function computes updated H parameter for each mpi rank
@@ -379,7 +454,7 @@ Source code for pyDNMFk.dist_nmf
[docs] def Fro_MU_update_W(self):
- r"""
+ r"""
Frobenius norm based multiplicative update of W parameter
Function computes updated H parameter for each mpi rank
@@ -398,7 +473,7 @@ Source code for pyDNMFk.dist_nmf
self.W_ij *= AH / WHTH
[docs] def Fro_MU_update(self, W_update=True):
- r"""
+ r"""
Frobenius norm based multiplicative update of W and H parameter
Function computes updated W and H parameter for each mpi rank
@@ -415,11 +490,11 @@ Source code for pyDNMFk.dist_nmf
self.Fro_MU_update_W()
self.Fro_MU_update_H()
- '''Functions for KL MU NMF update'''
+ '''Functions for KL MU NMF update'''
[docs] @comm_timing()
def gather_W_H(self, gW=True, gH=True):
- r"""
+ r"""
Gathers W and H factors across cartesian groups
i.e H_ij -> H_j if gH=True and W_ij -> W_i and gW=True
@@ -445,7 +520,7 @@ Source code for pyDNMFk.dist_nmf
[docs] @comm_timing()
def WTU_glob(self):
- r""" Distributed computation of W^TU
+ r""" Distributed computation of W^TU
Computes the global matrix multiplication of matrix W and U for KL
.. math:: W^TU
@@ -472,7 +547,7 @@ Source code for pyDNMFk.dist_nmf
[docs] @comm_timing()
def UHT_glob(self):
- r""" Distributed computation of UH^T
+ r""" Distributed computation of UH^T
Computes the global matrix multiplication of matrix W and U for KL
.. math:: UH^T
@@ -502,7 +577,7 @@ Source code for pyDNMFk.dist_nmf
return tmp
[docs] def KL_MU_update_W(self):
- r"""
+ r"""
KL divergence based multiplicative update of W parameter
Function computes updated W parameter for each mpi rank
@@ -522,7 +597,7 @@ Source code for pyDNMFk.dist_nmf
self.W_ij *= sk / (X2 + self.eps)
[docs] def KL_MU_update_H(self):
- r"""
+ r"""
Frobenius norm based multiplicative update of H parameter
Function computes updated H parameter for each mpi rank
@@ -542,7 +617,7 @@ Source code for pyDNMFk.dist_nmf
self.H_ij *= ks / (X1 + self.eps)
[docs] def KL_MU_update(self, W_update=True):
- r"""
+ r"""
KL divergence based multiplicative update of W and H parameter
Function computes updated W and H parameter for each mpi rank
@@ -559,10 +634,10 @@ Source code for pyDNMFk.dist_nmf
self.KL_MU_update_W()
self.KL_MU_update_H()
- '''Functions for FRO HALS NMF update'''
+ '''Functions for FRO HALS NMF update'''
[docs] def FRO_HALS_update_W(self):
- r"""
+ r"""
Frobenius norm minimization based HALS update of W parameter
Function computes updated W parameter for each mpi rank
@@ -585,7 +660,7 @@ Source code for pyDNMFk.dist_nmf
self.W_ij[:, kk] /= ss
[docs] def FRO_HALS_update_H(self):
- r"""
+ r"""
Frobenius norm minimization based HALS update of H parameter
Function computes updated H parameter for each mpi rank
@@ -605,7 +680,7 @@ Source code for pyDNMFk.dist_nmf
self.H_ij[kk, :] = np.maximum(temp_vec, self.eps)
[docs] def FRO_HALS_update(self, W_update=True):
- r"""
+ r"""
Frobenius norm minimization based HALS update of W and H parameter
Function computes updated W and H parameter for each mpi rank
@@ -622,11 +697,11 @@ Source code for pyDNMFk.dist_nmf
self.FRO_HALS_update_W()
self.FRO_HALS_update_H()
- '''Functions for FRO BCD NMF update'''
+ '''Functions for FRO BCD NMF update'''
[docs] @comm_timing()
def globalSqNorm(self, comm, X):
- """ Calc global squared norm of any matrix"""
+ """ Calc global squared norm of any matrix"""
normX = np.linalg.norm(X)
sqnormX = normX * normX
Xnorm = self.comm1.allreduce(sqnormX, op=MPI.SUM)
@@ -634,7 +709,7 @@ Source code for pyDNMFk.dist_nmf
[docs] @comm_timing()
def initWandH(self):
- """ Initialize the parameters for BCD updates"""
+ """ Initialize the parameters for BCD updates"""
# global Frobenius norm u calculated before
Xnorm = self.globalSqNorm(self.comm1, self.A_ij)
@@ -654,7 +729,7 @@ Source code for pyDNMFk.dist_nmf
return Wm, Hm, HHT, AHT, W_old, H_old, obj_old, Xnorm
[docs] def FRO_BCD_update(self, W_update=True, itr=1000):
- r"""
+ r"""
Frobenius norm minimization based BCD update of W and H parameter
Function computes updated W and H parameter for each mpi rank
@@ -680,7 +755,7 @@ Source code for pyDNMFk.dist_nmf
# relerr1 = relerr2 = []
# Iterate for max iterations:
for i in range(max_iters):
- """
+ """
W update
"""
HHTnorm_old = HHTnorm #
@@ -694,7 +769,7 @@ Source code for pyDNMFk.dist_nmf
self.comm1.barrier()
self.W_ij /= globalWSum
WTW = self.global_gram(self.W_ij)
- """
+ """
H update
"""
WTWnorm_old = WTWnorm
@@ -733,7 +808,7 @@ Source code for pyDNMFk.dist_nmf
[docs]class nmf_algorithms_1D():
- """
+ """
Performs the distributed NMF operation along 1D cartesian grid
Parameters
@@ -785,7 +860,7 @@ Source code for pyDNMFk.dist_nmf
self.local_H_n = self.H_j.shape[1]
[docs] def update(self):
- """Performs 1 step Update for factors W and H based on NMF method and corresponding norm minimization
+ """Performs 1 step Update for factors W and H based on NMF method and corresponding norm minimization
Returns
-------
@@ -814,7 +889,7 @@ Source code for pyDNMFk.dist_nmf
[docs] @comm_timing()
def global_gram(self, A, p=1):
- r""" Distributed gram computation
+ r""" Distributed gram computation
Computes the global gram operation of matrix A
.. math:: A^TA
@@ -839,7 +914,7 @@ Source code for pyDNMFk.dist_nmf
[docs] @comm_timing()
def global_mm(self, A, B, p=-1):
- r""" Distributed matrix multiplication
+ r""" Distributed matrix multiplication
Computes the global matrix multiplication of matrix A and B
.. math:: AB
@@ -855,7 +930,8 @@ Source code for pyDNMFk.dist_nmf
AB_glob : ndarray
"""
- AB_loc = np.matmul(A, B)
+ #AB_loc = np.matmul(A, B)
+ AB_loc = A@B
if p != 1:
AB_glob = self.comm1.allreduce(AB_loc, op=MPI.SUM)
self.comm1.barrier()
@@ -863,11 +939,11 @@ Source code for pyDNMFk.dist_nmf
AB_glob = AB_loc
return AB_glob
- '''Functions for Fro MU NMF update'''
+ '''Functions for Fro MU NMF update'''
[docs] @comm_timing()
def Fro_MU_update_W(self):
- r"""
+ r"""
Frobenius norm based multiplicative update of W parameter
Function computes updated H parameter for each mpi rank
@@ -879,15 +955,16 @@ Source code for pyDNMFk.dist_nmf
-------
self.W_i : ndarray
"""
+ HH_T = self.global_gram(self.H_j.T, p=self.p_c)
+ AH = self.global_mm(self.A_ij, self.H_j.T, p=self.p_c)
+ #WHTH = np.matmul(self.W_i, HH_T) + self.eps
+ WHTH = (self.W_i @ HH_T) + self.eps
+ self.W_i *= AH / WHTH
- W_TW = self.global_gram(self.W_i, p=self.p_r)
- AtW = self.global_mm(self.W_i.T, self.A_ij, p=self.p_r)
- HWtW = np.matmul(self.H_j.T, W_TW) + self.eps
- self.H_j *= AtW / HWtW.T
[docs] @comm_timing()
def Fro_MU_update_H(self):
- r"""
+ r"""
Frobenius norm based multiplicative update of H parameter
Function computes updated H parameter for each mpi rank
@@ -898,15 +975,21 @@ Source code for pyDNMFk.dist_nmf
Returns
-------
self.H_j : ndarray"""
+ W_TW = self.global_gram(self.W_i, p=self.p_r)
+ #print(f"[+] WTW = {W_TW}")
+ AtW = self.global_mm(self.W_i.T, self.A_ij, p=self.p_r)
+ #print(f"[+] ATW = {AtW}")
+ #HWtW = np.matmul(self.H_j.T, W_TW) + self.eps
+ HWtW = (self.H_j.T @ W_TW) + self.eps
+ #print(f"[+] HWTW = {HWtW.T}")
+ h_up = self.H_j*AtW / HWtW.T
+ #print(f"H_UPDATE = {h_up}")
+ self.H_j *= AtW / HWtW.T
- HH_T = self.global_gram(self.H_j.T, p=self.p_c)
- AH = self.global_mm(self.A_ij, self.H_j.T, p=self.p_c)
- WHTH = np.matmul(self.W_i, HH_T) + self.eps
- self.W_i *= AH / WHTH
[docs] @comm_timing()
def Fro_MU_update(self, W_update=True):
- r"""
+ r"""
Frobenius norm based multiplicative update of W and H parameter
Function computes updated W and H parameter for each mpi rank
@@ -923,11 +1006,11 @@ Source code for pyDNMFk.dist_nmf
self.Fro_MU_update_W()
self.Fro_MU_update_H()
- '''Functions for KL MU NMF update'''
+ '''Functions for KL MU NMF update'''
[docs] @comm_timing()
def sum_along_axis(self, X, p=1, axis=0):
- r"""
+ r"""
Performs sum of the matrix along given axis
Parameters
@@ -955,7 +1038,7 @@ Source code for pyDNMFk.dist_nmf
[docs] @comm_timing()
def glob_UX(self, axis):
- """Perform a global operation UX for W and H update with KL"""
+ """Perform a global operation UX for W and H update with KL"""
UX = self.A_ij / (self.W_i @ self.H_j + self.eps)
if axis == 1:
UX = self.global_mm(self.W_i.T, UX, p=self.p_r)
@@ -964,7 +1047,7 @@ Source code for pyDNMFk.dist_nmf
return UX
[docs] def KL_MU_update_W(self):
- r"""
+ r"""
KL divergence based multiplicative update of W parameter
Function computes updated W parameter for each mpi rank
@@ -983,7 +1066,7 @@ Source code for pyDNMFk.dist_nmf
self.W_i *= sk / (X2 + self.eps)
[docs] def KL_MU_update_H(self):
- r"""
+ r"""
KL divergence based multiplicative update of H parameter
Function computes updated H parameter for each mpi rank
@@ -1002,7 +1085,7 @@ Source code for pyDNMFk.dist_nmf
self.H_j *= sk / (X2 + self.eps)
[docs] def KL_MU_update(self, W_update=True):
- r"""
+ r"""
KL divergence based multiplicative update of W and H parameter
Function computes updated W and H parameter for each mpi rank
@@ -1021,10 +1104,10 @@ Source code for pyDNMFk.dist_nmf
self.KL_MU_update_W()
self.KL_MU_update_H()
- '''Functions for FRO HALS NMF update'''
+ '''Functions for FRO HALS NMF update'''
[docs] def FRO_HALS_update_W(self):
- r"""
+ r"""
Frobenius norm minimization based HALS update of W parameter
Function computes updated W parameter for each mpi rank
@@ -1046,7 +1129,7 @@ Source code for pyDNMFk.dist_nmf
self.W_i[:, kk] /= ss
[docs] def FRO_HALS_update_H(self):
- r"""
+ r"""
Frobenius norm minimization based HALS update of H parameter
Function computes updated H parameter for each mpi rank
@@ -1067,7 +1150,7 @@ Source code for pyDNMFk.dist_nmf
# self.H_j = self.H_j.T
[docs] def FRO_HALS_update(self, W_update=True):
- r"""
+ r"""
Frobenius norm minimizatio based HALS update of W and H parameter
Function computes updated W and H parameter for each mpi rank
@@ -1087,11 +1170,11 @@ Source code for pyDNMFk.dist_nmf
self.FRO_HALS_update_H()
- '''Functions for FRO BCD NMF update'''
+ '''Functions for FRO BCD NMF update'''
[docs] @comm_timing()
def globalSqNorm(self, X, p=-1):
- """ Calc global squared norm of any matrix"""
+ """ Calc global squared norm of any matrix"""
normX = np.linalg.norm(X)
sqnormX = normX * normX
if p != 1:
@@ -1103,7 +1186,7 @@ Source code for pyDNMFk.dist_nmf
[docs] @comm_timing()
def initWandH(self):
- """ Initialize the parameters for BCD updates"""
+ """ Initialize the parameters for BCD updates"""
# global Frobenius norm u calculated before
Xnorm = self.globalSqNorm(self.A_ij)
@@ -1122,7 +1205,7 @@ Source code for pyDNMFk.dist_nmf
return Wm, Hm, HHT, AHT, W_old, H_old, obj_old, Xnorm
[docs] def FRO_BCD_update(self, W_update=True, itr=1000):
- r"""
+ r"""
Frobenius norm minimization based BCD update of W and H parameter
Function computes updated W and H parameter for each mpi rank
@@ -1147,7 +1230,7 @@ Source code for pyDNMFk.dist_nmf
# relerr1 = relerr2 = []
# Iterate for max iterations:
for i in range(max_iters):
- """
+ """
W update
"""
HHTnorm_old = HHTnorm #
@@ -1163,7 +1246,7 @@ Source code for pyDNMFk.dist_nmf
if self.p_r == 1: globalWSum = localWsum
self.W_i /= globalWSum
WTW = self.global_gram(self.W_i, p=self.p_r)
- """
+ """
H update
"""
WTWnorm_old = WTWnorm
@@ -1200,47 +1283,73 @@ Source code for pyDNMFk.dist_nmf
obj_old = obj
-
-
-
-
+
+
+
+
-
-
-
+
+
+
-
-
+
+
+
-
+
+
\ No newline at end of file
diff --git a/docs/_modules/pyDNMFk/dist_svd.html b/docs/_modules/pyDNMFk/dist_svd.html
index c98ad74..4fe57ee 100644
--- a/docs/_modules/pyDNMFk/dist_svd.html
+++ b/docs/_modules/pyDNMFk/dist_svd.html
@@ -1,155 +1,230 @@
-
-
-
+
+
+
+
+
+
+
+ pyDNMFk.dist_svd — pyDNMFk 0.0.1 documentation
-
- pyDNMFk.dist_svd — pyDNMFk 1.0.0 documentation
-
+
-
-
+
+
+
+
+
+
+
+
+
+
+
+
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-
+ Skip to main content
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
-
+
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-
-
-
- - »
-
- - Module code »
-
- - pyDNMFk.dist_svd
-
-
- -
-
-
-
-
+
+
+
+
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Source code for pyDNMFk.dist_svd
# @Author: Gopinath Chennupati, Manish Bhattarai, Erik Skau
from datetime import datetime
@@ -160,7 +235,7 @@ Source code for pyDNMFk.dist_svd
[docs]class DistSVD():
- r"""
+ r"""
Distributed Computation of SVD along 1D distribution of the data. Only U or V is distributed based on data size.
Parameters
@@ -220,7 +295,7 @@ Source code for pyDNMFk.dist_svd
[docs] @comm_timing()
def normalize_by_W(self, Wall, Hall, comm1):
- """Normalize the factors W and H"""
+ """Normalize the factors W and H"""
Wall_norm = Wall.sum(axis=0, keepdims=True)
if self.proc_rows != 1:
Wall_norm = comm1.allreduce(Wall_norm, op=MPI.SUM)
@@ -232,7 +307,7 @@ Source code for pyDNMFk.dist_svd
[docs] @comm_timing()
def randomUnitVector(self, d):
- """
+ """
Construnct a rondom unit vector
"""
unnormalized = [normalvariate(0, 1) for _ in range(d)]
@@ -241,14 +316,14 @@ Source code for pyDNMFk.dist_svd
[docs] @comm_timing()
def globalGram(self, X, Y):
- """Compute the global gram betwee X and Y"""
+ """Compute the global gram betwee X and Y"""
B = X @ Y
B = self.grid_comm.allreduce(B)
return B
[docs] @comm_timing()
def svd1D(self):
- """
+ """
One dimensional SVD
"""
buf = self.randomUnitVector(min(self.globalm, self.globaln))
@@ -291,7 +366,7 @@ Source code for pyDNMFk.dist_svd
[docs] @comm_timing()
def calc_norm(self, vec):
- """Compute the norm of vector"""
+ """Compute the norm of vector"""
partial_sq_sum = sum(vec * vec)
global_sq_sum = self.grid_comm.allreduce(partial_sq_sum)
# norm = torch.sqrt(global_sq_sum)
@@ -300,7 +375,7 @@ Source code for pyDNMFk.dist_svd
[docs] @comm_timing()
def svd(self):
- """
+ """
Computes the SVD for a given matrix
Returns
@@ -340,7 +415,7 @@ Source code for pyDNMFk.dist_svd
[docs] @comm_timing()
def rel_error(self, U, S, V):
- """Computes the relative error between the reconstructed data with factors vs original data"""
+ """Computes the relative error between the reconstructed data with factors vs original data"""
X_recon = U @ S @ V
err_num = np.sum((self.A - X_recon) ** 2)
norm_deno = np.sum(self.A ** 2)
@@ -351,7 +426,7 @@ Source code for pyDNMFk.dist_svd
[docs] @comm_timing()
def nnsvd(self, flag=1, verbose=1):
- r"""
+ r"""
Computes the distributed Non-Negative SVD(NNSVD) components from the computed SVD factors.
Parameters
@@ -420,47 +495,73 @@ Source code for pyDNMFk.dist_svd
return self.normalize_by_W(W, H, self.grid_comm)
-
-
-
-
+
+
+
+
-
-
-
+
+
+
-
-
+
+
+
-
+
+
\ No newline at end of file
diff --git a/docs/_modules/pyDNMFk/plot_results.html b/docs/_modules/pyDNMFk/plot_results.html
index ffd30d7..61f54c0 100644
--- a/docs/_modules/pyDNMFk/plot_results.html
+++ b/docs/_modules/pyDNMFk/plot_results.html
@@ -1,164 +1,239 @@
-
-
-
+
+
+
+
+
+
+
+ pyDNMFk.plot_results — pyDNMFk 0.0.1 documentation
-
- pyDNMFk.plot_results — pyDNMFk 1.0.0 documentation
-
+
-
-
+
+
+
+
+
+
+
+
+
+
+
+
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-
+ Skip to main content
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
-
+
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
+
-
-
- - »
-
- - Module code »
-
- - pyDNMFk.plot_results
-
-
- -
-
-
-
-
+
+
+
+
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Source code for pyDNMFk.plot_results
# @author: Manish Bhattarai
import matplotlib
from matplotlib import pyplot as plt
from .data_io import *
-
+from h5py import File
[docs]def plot_err(err):
- """Plots the relative error for NMF decomposition as a function of number of iterations"""
+ """Plots the relative error for NMF decomposition as a function of number of iterations"""
idx = np.linspace(1, len(err), len(err))
plt.plot(idx, err)
plt.xlabel('Iterations')
@@ -169,7 +244,7 @@ Source code for pyDNMFk.plot_results
[docs]def read_plot_factors(factors_path, pgrid):
- """Reads the factors W and H and Plots them"""
+ """Reads the factors W and H and Plots them"""
W, H = read_factors(factors_path, pgrid)
plot_W(W)
plt.savefig(factors_path + 'W.png')
@@ -178,7 +253,7 @@ Source code for pyDNMFk.plot_results
[docs]def plot_W(W):
- """Reads a factor and plots into subplots for each component"""
+ """Reads a factor and plots into subplots for each component"""
m, k = W.shape
params = {'legend.fontsize': 60,
@@ -215,10 +290,10 @@ Source code for pyDNMFk.plot_results
plt.show()
-[docs]def plot_results(startProcess, endProcess, RECON, RECON1, SILL_MIN, out_put, name):
- """Plots the relative error and Silhouette results for estimation of k"""
+[docs]def plot_results(startProcess, endProcess, stepProcess,RECON, RECON1, SILL_MIN, out_put, name):
+ """Plots the relative error and Silhouette results for estimation of k"""
######################################## Plotting ####################################################
- t = range(startProcess, endProcess + 1)
+ t = range(startProcess, endProcess + 1,stepProcess)
fig, ax1 = plt.subplots(num=None, figsize=(10, 6), dpi=300, facecolor='w', edgecolor='k')
title = 'Num'
color = 'tab:red'
@@ -235,9 +310,7 @@ Source code for pyDNMFk.plot_results
# manipulate the y-axis values into percentage
vals = ax1.get_yticks()
ax1.set_yticklabels(['{:,.0%}'.format(x) for x in vals])
-
# ax1.legend(loc=0)
-
ax2 = ax1.twinx() # instantiate a second axes that shares the same x-axis
color = 'tab:blue'
ax2.set_ylabel('Minimum Stability', color=color) # we already handled the x-label with ax1
@@ -246,19 +319,61 @@ Source code for pyDNMFk.plot_results
# ax2.legend(loc=1)
fig.tight_layout() # otherwise the right y-label is slightly clipped
# plt.show()
-
# added these three lines
lns = lns1 + lns2 + lns3
labs = [l.get_label() for l in lns]
ax1.legend(lns, labs, loc=0)
-
plt.savefig(out_put + '/' + name + '_selection_plot.pdf')
-
plt.close()
+[docs]def plot_results_fpath(params):
+ """Plots the relative error and Silhouette results for estimation of k from given folder location"""
+ ######################################## Plotting ####################################################
+ t = range(params.start_k,params.end_k + 1,params.step)
+ fig, ax1 = plt.subplots(num=None, figsize=(10, 6), dpi=300, facecolor='w', edgecolor='k')
+ title = 'Num'
+ color = 'tab:red'
+ ax1.set_xlabel('Total Signatures')
+ ax1.set_ylabel('Mean L2 %', color=color)
+ ax1.set_title(title)
+ RECON = []
+ RECON1 = []
+ SILL_MIN = []
+ for k in t:
+ results_paths = params.results_path + str(k) + '/'
+ data = File(results_paths+'/results.h5','r')
+ RECON.append(np.array(data['L_errDist']))
+ RECON1.append(np.array(data['avgErr']))
+ SILL_MIN.append(round(np.min(np.array(data['clusterSilhouetteCoefficients'])), 2))
+
+ lns1 = ax1.plot(t, RECON, marker='o', linestyle=':', color=color, label='Mean L2 %')
+ lns3 = ax1.plot(t, RECON1, marker='X', linestyle=':', color='tab:green', label="Relative error %")
+ ax1.tick_params(axis='y', labelcolor=color)
+ ax1.xaxis.set_ticks(np.arange(min(t), max(t) + 1, 1))
+ # ax1.axvspan(shadow_start, shadow_end, alpha=0.20, color='#ADD8E6')
+ # ax1.axvspan(shadow_alternative_start, shadow_alternative_end, alpha=0.20, color='#696969')
+ # manipulate the y-axis values into percentage
+ vals = ax1.get_yticks()
+ ax1.set_yticklabels(['{:,.0%}'.format(x) for x in vals])
+ # ax1.legend(loc=0)
+ ax2 = ax1.twinx() # instantiate a second axes that shares the same x-axis
+ color = 'tab:blue'
+ ax2.set_ylabel('Minimum Stability', color=color) # we already handled the x-label with ax1
+ lns2 = ax2.plot(t, SILL_MIN, marker='s', linestyle="-.", color=color, label='Minimum Stability')
+ ax2.tick_params(axis='y', labelcolor=color)
+ # ax2.legend(loc=1)
+ fig.tight_layout() # otherwise the right y-label is slightly clipped
+ # plt.show()
+ # added these three lines
+ lns = lns1 + lns2 + lns3
+ labs = [l.get_label() for l in lns]
+ ax1.legend(lns, labs, loc=0)
+ plt.savefig(params.results_path + '/' + params.fname + '_selection_plot.pdf')
+ plt.close()
+
[docs]def box_plot(dat, respath):
- """Plots the boxplot from the given data and saves the results"""
+ """Plots the boxplot from the given data and saves the results"""
dat.plot.bar()
plt.xlabel('operation')
plt.ylabel('timing(sec)')
@@ -268,7 +383,7 @@ Source code for pyDNMFk.plot_results
[docs]def timing_stats(fpath):
- """Reads the timing stats dictionary from the stored file and parses the data. """
+ """Reads the timing stats dictionary from the stored file and parses the data. """
import copy
data = pd.read_csv(fpath).iloc[0, 1:]
breakdown_level_2 = {'init': ['__init__', 'init_factors'],
@@ -286,7 +401,7 @@ Source code for pyDNMFk.plot_results
'compute': 'fit'}
results = {}
- ''''Data parsing'''
+ ''''Data parsing'''
breakdown_level_1_dat = copy.deepcopy(breakdown_level_1)
breakdown_level_2_dat = copy.deepcopy(breakdown_level_2)
@@ -315,7 +430,7 @@ Source code for pyDNMFk.plot_results
[docs]def plot_timing_stats(fpath, respath):
- ''' Plots the timing stats for the MPI operation.
+ ''' Plots the timing stats for the MPI operation.
fpath: Stats data path
respath: Path to save graph'''
res1, res2 = timing_stats(fpath)
@@ -327,47 +442,73 @@ Source code for pyDNMFk.plot_results
box_plot(pd.DataFrame([tmp]).loc[0, :], respath)
-
-
-
-
+
+
+
+
-
-
-
+
+
+
-
-
+
+
+
-
+
+
\ No newline at end of file
diff --git a/docs/_modules/pyDNMFk/pyDNMF.html b/docs/_modules/pyDNMFk/pyDNMF.html
index 56fb7fe..f66ff00 100644
--- a/docs/_modules/pyDNMFk/pyDNMF.html
+++ b/docs/_modules/pyDNMFk/pyDNMF.html
@@ -1,155 +1,230 @@
-
-
-
+
+
+
+
+
+
+
+ pyDNMFk.pyDNMF — pyDNMFk 0.0.1 documentation
-
- pyDNMFk.pyDNMF — pyDNMFk 1.0.0 documentation
-
+
-
-
+
+
+
+
+
+
+
+
+
+
+
+
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-
+ Skip to main content
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
-
+
-
-
-
-
+
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
+
-
-
- - »
-
- - Module code »
-
- - pyDNMFk.pyDNMF
-
-
- -
-
-
-
-
+
+
+
+
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Source code for pyDNMFk.pyDNMF
# @author: Manish Bhattarai
@@ -157,10 +232,12 @@ Source code for pyDNMFk.pyDNMF
from .dist_nmf import *
from .dist_svd import *
from .utils import *
-
+import scipy
+from scipy.linalg import norm as scipy_norm
+from scipy.sparse.linalg import norm as scipy_sparse_norm
[docs]class PyNMF():
- r"""
+ r"""
Performs the distributed NMF decomposition of given matrix X into factors W and H
Parameters
@@ -219,9 +296,12 @@ Source code for pyDNMFk.pyDNMF
self.params.eps = self.eps
self.norm = var_init(self.params,'norm',default='kl')
self.method = var_init(self.params,'method',default='mu')
+ self.prune = var_init(self.params, 'prune', default=True)
self.save_factors = save_factors
self.params.itr = var_init(self.params,'itr',default=5000)
self.itr = self.params.itr
+ self.W_start, self.W_end = 0,0
+ self.H_start, self.H_end = 0, 0
try:
self.W_update = self.params.W_update
except:
@@ -231,7 +311,9 @@ Source code for pyDNMFk.pyDNMF
self.topo = '2d'
else:
self.topo = '1d'
- self.compute_global_dim()
+ self.params.topo = self.topo
+ self.data_op = data_operations(self.A_ij, self.params)
+ self.params = self.data_op.params
if factors is not None:
if self.topo == '1d':
self.W_i = factors[0].astype(self.A_ij.dtype)
@@ -241,45 +323,23 @@ Source code for pyDNMFk.pyDNMF
self.H_ij = factors[1].astype(self.A_ij.dtype)
else:
self.init_factors()
+ if self.prune:
+ if self.topo=='2d': self.A_ij,self.W_ij,self.H_ij = self.data_op.prune_all(self.W_ij,self.H_ij)
+ elif self.topo == '1d': self.A_ij,self.W_i,self.H_j = self.data_op.prune_all(self.W_i, self.H_j)
+
-[docs] @comm_timing()
- def compute_global_dim(self):
- """Computes global dimensions m and n from given chunk sizes for any grid configuration"""
- self.loc_m, self.loc_n = self.A_ij.shape
- if self.p_r != 1 and self.p_c == 1:
- self.params.n = self.loc_n
- self.params.m = self.comm1.allreduce(self.loc_m)
- elif self.p_c != 1 and self.p_r == 1:
- self.params.n = self.comm1.allreduce(self.loc_n)
- self.params.m = self.loc_m
- else:
- if self.rank % self.p_c == 0:
- self.params.m = self.loc_m
- else:
- self.params.m = 0
- self.params.m = self.comm1.allreduce(self.params.m)
- if self.rank // self.p_c == 0:
- self.params.n = self.loc_n
- else:
- self.params.n = 0
- self.params.n = self.comm1.allreduce(self.params.n)
- self.comm1.barrier()
- # if self.rank == 0: print('Data dimensions=(', self.params.m, self.params.n, ')')
[docs] @comm_timing()
def init_factors(self):
- """Initializes NMF factors with rand/nnsvd method"""
+ """Initializes NMF factors with rand/nnsvd method"""
if self.init == 'rand':
- if self.topo == '2d':
- dtr_blk_m = determine_block_params(self.cart_1d_column, (self.p_c, 1), (self.A_ij.shape[0], self.k))
- m_loc = dtr_blk_m.determine_block_shape_asymm()[0]
- dtr_blk_n = determine_block_params(self.cart_1d_row, (1, self.p_r), (self.k, self.A_ij.shape[1]))
- n_loc = dtr_blk_n.determine_block_shape_asymm()[1]
- self.W_ij = np.random.rand(m_loc, self.k).astype(self.A_ij.dtype)
- self.H_ij = np.random.rand(self.k, n_loc).astype(self.A_ij.dtype)
+ if self.topo=='2d':
+ self.W_ij = np.random.rand(self.params.m_loc, self.k).astype(self.A_ij.dtype)
+ self.H_ij = np.random.rand(self.k, self.params.n_loc).astype(self.A_ij.dtype)
elif self.topo == '1d':
+ #np.random.seed(0)
if self.p_c == 1:
self.W_i = np.random.rand(self.m_loc, self.k).astype(self.A_ij.dtype)
if self.rank == 0:
@@ -302,9 +362,14 @@ Source code for pyDNMFk.pyDNMF
elif self.topo == '2d':
raise Exception('NNSVD init only available for 1D topology, please try with 1d topo.')
+ #if self.rank == 0:
+ # print("=============== [iter %04d] =================" %-1)
+ # print(f'W_gpu = {self.W_i}')
+ # print(f'H_gpu = {self.H_j}')
+
[docs] @comm_timing()
def fit(self):
- r"""
+ r"""
Calls the sub routines to perform distributed NMF decomposition with initialization for a given norm minimization and update method
Returns
@@ -317,6 +382,7 @@ Source code for pyDNMFk.pyDNMF
Reconstruction error for NMF decomposition
"""
for i in range(self.itr):
+ #if self.rank == 0: print("=============== [iter %04d] =================" %i)
if self.method.lower() == 'bcd': i = self.itr - 1
if self.topo == '2d':
self.W_ij, self.H_ij = nmf_algorithms_2D(self.A_ij, self.W_ij, self.H_ij, params=self.params).update()
@@ -326,41 +392,47 @@ Source code for pyDNMFk.pyDNMF
if i == self.itr - 1:
self.W_ij, self.H_ij = self.normalize_features(self.W_ij, self.H_ij)
self.relative_err()
- if self.verbose == True:
- if self.rank == 0: print('relative error is:', self.recon_err)
+ #if self.verbose == True:
+ # #if self.rank == 0: print('relative error is:', self.recon_err)
if self.save_factors:
data_write(self.params).save_factors([self.W_ij, self.H_ij])
self.comm.Free()
+ if self.prune: self.W_ij,self.H_ij = self.data_op.unprune_factors(self.W_ij, self.H_ij)
return self.W_ij, self.H_ij, self.recon_err
elif self.topo == '1d':
self.W_i, self.H_j = nmf_algorithms_1D(self.A_ij, self.W_i, self.H_j, params=self.params).update()
+ #if self.rank == 0:
+ # print(f'W_gpu = {self.W_i}')
+ # print(f'H_gpu = {self.H_j}')
if i % 10 == 0:
self.H_j = np.maximum(self.H_j, self.eps)
self.W_i = np.maximum(self.W_i, self.eps)
if i == self.itr - 1:
self.W_i, self.H_j = self.normalize_features(self.W_i, self.H_j)
self.relative_err()
- if self.verbose == True:
- if self.rank == 0: print('\nrelative error is:', self.recon_err)
+ #if self.verbose == True: if self.rank == 0: print('\nrelative error is:', self.recon_err)
+ # if self.rank == 0: print('\nrelative error is:', self.recon_err)
if self.save_factors:
data_write(self.params).save_factors([self.W_i, self.H_j])
+ if self.prune:
+ self.W_i, self.H_j = self.data_op.unprune_factors(self.W_i, self.H_j)
return self.W_i, self.H_j, self.recon_err
[docs] @comm_timing()
def normalize_features(self, Wall, Hall):
- """Normalizes features Wall and Hall"""
- Wall_norm = Wall.sum(axis=0, keepdims=True) + self.eps
+ """Normalizes features Wall and Hall"""
+ Wall_norm = Wall.sum(axis=0, keepdims=True)
if self.topo == '2d':
Wall_norm = self.comm1.allreduce(Wall_norm, op=MPI.SUM)
elif self.topo == '1d':
if self.p_r != 1: Wall_norm = self.comm1.allreduce(Wall_norm, op=MPI.SUM)
- Wall /= Wall_norm
+ Wall /= Wall_norm+ self.eps
Hall *= Wall_norm.T
return Wall, Hall
[docs] @comm_timing()
def cart_2d_collect_factors(self):
- """Collects factors along each sub communicators"""
+ """Collects factors along each sub communicators"""
self.H_j = self.cart_1d_row.allgather(self.H_ij)
self.H_j = np.hstack((self.H_j))
self.W_i = self.cart_1d_column.allgather(self.W_ij)
@@ -368,7 +440,7 @@ Source code for pyDNMFk.pyDNMF
[docs] @comm_timing()
def relative_err(self):
- """Computes the relative error for NMF decomposition"""
+ """Computes the relative error for NMF decomposition"""
if self.topo == '2d': self.cart_2d_collect_factors()
self.glob_norm_err = self.dist_norm(self.A_ij - self.W_i @ self.H_j)
self.glob_norm_A = self.dist_norm(self.A_ij)
@@ -376,15 +448,20 @@ Source code for pyDNMFk.pyDNMF
[docs] @comm_timing()
def dist_norm(self, X, proc=-1, norm='fro', axis=None):
- """Computes the distributed norm"""
- nm = np.linalg.norm(X, axis=axis, ord=norm)
+ """Computes the distributed norm"""
+ if type(X) in [numpy.matrix, numpy.ndarray]:
+ nm = np.linalg.norm(X, axis=axis, ord=norm)
+ elif type(X) in [scipy.sparse.coo.coo_matrix, scipy.sparse.csr.csr_matrix, scipy.sparse.csc.csc_matrix]:
+ nm = scipy_sparse_norm(X, ord=norm, axis=axis)
+ else:
+ raise Exception("[!!] type(X): {} is not understaood".format(type(X)))
if proc != 1:
nm = self.comm1.allreduce(nm ** 2)
return np.sqrt(nm)
[docs] @comm_timing()
def column_err(self):
- """Computes the distributed column wise norm"""
+ """Computes the distributed column wise norm"""
dtr_blk = determine_block_params(self.comm1, (self.p_r, self.p_c), (self.params.m, self.params.n))
dtr_blk_idx = dtr_blk.determine_block_index_range_asymm()
dtr_blk_shp = dtr_blk.determine_block_shape_asymm()
@@ -394,8 +471,12 @@ Source code for pyDNMFk.pyDNMF
L_errDist_deno = np.zeros(self.n_loc)
Arecon = self.W_i @ self.H_j
for q in range(self.A_ij.shape[1]):
- L_errDist_num[q] = np.sum((self.A_ij[:, q] - Arecon[:, q]) ** 2)
- L_errDist_deno[q] = np.sum(self.A_ij[:, q] ** 2)
+ if scipy.sparse.issparse(self.A_ij):
+ L_errDist_num[q] = (np.square(self.A_ij[:, q] - Arecon[:, q][:,None])).sum()
+ L_errDist_deno[q] = (self.A_ij[:, q].data ** 2).sum()
+ else:
+ L_errDist_num[q] = ((self.A_ij[:, q] - Arecon[:, q]) ** 2).sum()
+ L_errDist_deno[q] = (self.A_ij[:, q] ** 2).sum()
col_err_num[dtr_blk_idx[0][1]:dtr_blk_idx[0][1] + dtr_blk_shp[1]] = L_errDist_num
col_err_deno[dtr_blk_idx[0][1]:dtr_blk_idx[0][1] + dtr_blk_shp[1]] = L_errDist_deno
col_err_num = self.comm1.allreduce(col_err_num)
@@ -404,47 +485,73 @@ Source code for pyDNMFk.pyDNMF
return col_err
-
-
-
-
+
+
+
+
-
-
-
+
+
+
-
-
+
+
+
-
+
+
\ No newline at end of file
diff --git a/docs/_modules/pyDNMFk/pyDNMFk.html b/docs/_modules/pyDNMFk/pyDNMFk.html
index 5e4ab66..9bdb298 100644
--- a/docs/_modules/pyDNMFk/pyDNMFk.html
+++ b/docs/_modules/pyDNMFk/pyDNMFk.html
@@ -1,155 +1,230 @@
-
-
-
+
+
+
+
+
+
+
+ pyDNMFk.pyDNMFk — pyDNMFk 0.0.1 documentation
-
- pyDNMFk.pyDNMFk — pyDNMFk 1.0.0 documentation
-
+
-
-
+
+
+
+
+
+
+
+
+
+
+
+
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-
+ Skip to main content
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
-
+
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
+
-
-
- - »
-
- - Module code »
-
- - pyDNMFk.pyDNMFk
-
-
- -
-
-
-
-
+
+
+
+
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
Source code for pyDNMFk.pyDNMFk
#@author: Manish Bhattarai
from scipy.stats import wilcoxon
@@ -158,8 +233,32 @@ Source code for pyDNMFk.pyDNMFk
from .pyDNMF import *
from .plot_results import *
+import cupy as cp
+from cupy import asarray as _asarray
+from cupy import asnumpy as _asnumpy
+from cupy import divide as _divide
+from cupy import matmul as _matmul
+from cupy import multiply as _multiply
+from cupy import zeros as _zeros
+from .cupyCuSPARSELib import spMM as _spMM
+from .cupyCuSPARSELib import spMat, spRandMat, spMM
+from scipy.sparse import coo_matrix, csc_matrix, csr_matrix, dia_matrix, issparse
+import scipy
+#from .NMF_utils import sparse_NMF_generator
+from .cudaNMF import cudaNMF
+#from sgemm import sgemm as _sgemm
+
+from cupy.cuda import nccl
+from .communicators import NCCLComm
+from .toolz import log, blue, green, red, amber
+from .comm_utils import GetTopology
+
+import time
+
+
+
[docs]class sample():
- """
+ """
Generates perturbed version of data based on sampling distribution.
Parameters
@@ -188,7 +287,7 @@ Source code for pyDNMFk.pyDNMFk
[docs] @comm_timing()
def randM(self):
- """
+ """
Multiplies each element of X by a uniform random number in (1-epsilon, 1+epsilon).
"""
@@ -198,13 +297,13 @@ Source code for pyDNMFk.pyDNMFk
[docs] @comm_timing()
def poisson(self):
- """Resamples each element of a matrix from a Poisson distribution with the mean set by that element. Y_{i,j} = Poisson(X_{i,j}"""
+ """Resamples each element of a matrix from a Poisson distribution with the mean set by that element. Y_{i,j} = Poisson(X_{i,j}"""
self.X_per = np.random.poisson(self.X).astype(self.X.dtype)
[docs] @comm_timing()
def fit(self):
- r"""
+ r"""
Calls the sub routines to perform resampling on data
Returns
@@ -221,7 +320,7 @@ Source code for pyDNMFk.pyDNMFk
[docs]class PyNMFk():
- r"""
+ r"""
Performs the distributed NMF decomposition with custom clustering for estimating hidden factors k
Parameters
@@ -278,10 +377,21 @@ Source code for pyDNMFk.pyDNMFk