diff --git a/dimod/generators/__init__.py b/dimod/generators/__init__.py index f5f61d8da..1335befbe 100644 --- a/dimod/generators/__init__.py +++ b/dimod/generators/__init__.py @@ -13,15 +13,16 @@ # limitations under the License. from dimod.generators.anti_crossing import * +from dimod.generators.binpacking import * from dimod.generators.chimera import * from dimod.generators.constraints import * from dimod.generators.fcl import * from dimod.generators.gates import * from dimod.generators.graph import * from dimod.generators.integer import * -from dimod.generators.random import * from dimod.generators.knapsack import * -from dimod.generators.binpacking import * +from dimod.generators.magic_square import * from dimod.generators.multi_knapsack import * +from dimod.generators.random import * from dimod.generators.satisfiability import * -from dimod.generators.magic_square import * \ No newline at end of file +from dimod.generators.wireless import * \ No newline at end of file diff --git a/dimod/generators/wireless.py b/dimod/generators/wireless.py new file mode 100644 index 000000000..6fc5188d1 --- /dev/null +++ b/dimod/generators/wireless.py @@ -0,0 +1,985 @@ +# Copyright 2023 D-Wave Systems Inc. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express 93or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +from collections import namedtuple +from functools import wraps +from itertools import product +from typing import Literal, Optional, Tuple, Union + +import numpy as np + +import dimod + +__all__ = ['mimo', 'coordinated_multipoint', ] + +mod_params = namedtuple("mod_params", ["bits_per_transmitter", + "amplitudes", + "transmitters_per_bit", + "number_of_amps", + "bits_per_symbol"]) +mod_config = { + "BPSK": mod_params(1, 1, 1, 1, 1), + "QPSK": mod_params(2, 1, 2, 1, 1), + "16QAM": mod_params(4, 2, 4, 2, 2), + "64QAM": mod_params(6, 4, 6, 3, 3), + "256QAM": mod_params(8, 8, 8, 5, 4)} + +def _quadratic_form(y, F): + """Convert :math:`O(v) = ||y - F v||^2` to sparse quadratic form. + + Constructs coefficients for the form + :math:`O(v) = v^{\dagger} J v - 2 \Re(h^{\dagger} v) + \\text{offset}`. + + Args: + y: Received signal as a NumPy column vector of complex or real values. + + F: Wireless channel as an :math:`i \\times j` NumPy matrix of complex + values, where :math:`i` rows correspond to :math:`y_i` receivers + and :math:`j` columns correspond to :math:`v_i` transmitted symbols. + + Returns: + Three tuple of offset, as a real scalar, linear biases :math:`h`, as a dense + real vector, and quadratic interactions, :math:`J`, as a dense real symmetric + matrix. + """ + + if len(F.shape) != 2 or F.shape[0] != y.shape[0]: + raise ValueError("F should have shape (n, m) for some m, n " + "and n should equal y.shape[1];" + f" given: {F.shape}, n={y.shape[1]}") + + offset = np.matmul(y.imag.T, y.imag) + np.matmul(y.real.T, y.real) + h = -2 * np.matmul(F.T.conj(), y) + J = np.matmul(F.T.conj(), F) + + return offset, h, J + +def _real_quadratic_form(h, J, modulation=None): + """Separate real and imaginary parts of quadratic form. + + Unwraps objective function on complex variables as an objective function of + concatenated real variables, first the real and then the imaginary part. + + Args: + h: Linear biases as a dense NumPy vector. + + J: Quadratic interactions as a dense symmetric matrix. + + modulation: Modulation. Supported values are 'BPSK', 'QPSK', '16QAM', + '64QAM', and '256QAM'. + + Returns: + Two-tuple of linear biases, :math:`h`, as a NumPy real vector with any + imaginary part following the real part, and quadratic interactions, + :math:`J`, as a real matrix with any imaginary part moved to above and + below the diagonal. + """ + # Here, for BPSK F-induced complex parts of h and J are discarded: + # Given y = F x + nu, for independent and identically distributed channel F + # and complex noise nu, the system of equations defined by the real part is + # sufficient to define the canonical decoding problem. + # In essence, rotate y to the eigenbasis of F, throw away the orthogonal noise + # (the complex problem as the real part with suitable adjustment factor 2 to + # signal to noise ratio: F^{-1}*y = I*x + F^{-1}*nu) + # JR: revisit and prove + + if modulation and modulation not in mod_config: + raise ValueError(f"Unsupported modulation: {modulation}") + + if modulation != 'BPSK' and (np.iscomplex(h).any() or np.iscomplex(J).any()): + hR = np.concatenate((h.real, h.imag), axis=0) + JR = np.concatenate((np.concatenate((J.real, J.imag), axis=0), + np.concatenate((J.imag.T, J.real), axis=0)), + axis=1) + return hR, JR + return h.real, J.real + +def _amplitude_modulated_quadratic_form(h, J, modulation="BPSK"): + """Amplitude-modulate the quadratic form. + + Updates bias amplitudes for quadrature amplitude modulation. + + Args: + h: Linear biases as a NumPy vector. + + J: Quadratic interactions as a matrix. + + modulation: Modulation. Supported values are non-quadrature modulation + 'BPSK' (the default) and quadrature modulations 'QPSK', '16QAM', + '64QAM', and '256QAM'. + + Returns: + Two-tuple of amplitude-modulated linear biases, :math:`h`, as a NumPy + vector and amplitude-modulated quadratic interactions, :math:`J`, as + a matrix. + """ + if modulation not in mod_config: + raise ValueError(f"Unsupported modulation: {modulation}") + + amps = 2 ** np.arange(mod_config[modulation].number_of_amps) + hA = np.kron(amps[:, np.newaxis], h) + JA = np.kron(np.kron(amps[:, np.newaxis], amps[np.newaxis, :]), J) + return hA, JA + +def _symbols_to_bits(symbols, modulation="BPSK"): + """Convert quadrature amplitude modulated (QAM) symbols to bits. + + Args: + symbols: Transmitted symbols as a NumPy column vector. + + modulation: Modulation. Supported values are the default non-quadrature + modulation, binary phase-shift keying (BPSK, or 2-QAM), and + quadrature modulations 'QPSK', '16QAM', '64QAM', and '256QAM'. + + Returns: + Bits as a NumPy array. + """ + if modulation not in mod_config: + raise ValueError(f"Unsupported modulation: {modulation}") + + if modulation == 'BPSK': + return symbols.copy() + + if modulation == 'QPSK': + return np.concatenate((symbols.real, symbols.imag)) + + bits_per_real_symbol = mod_config[modulation].bits_per_symbol + + # A map from integer parts to real is clearest (and sufficiently performant), + # generalizes to Gray coding more easily as well: + + symb_to_bits = {np.sum([x * 2**xI for xI, x in enumerate(bits)]): bits + for bits in product(*bits_per_real_symbol * [(-1, 1)])} + bits = np.concatenate( + [np.concatenate(([symb_to_bits[symb][prec] for symb in symbols.real.flatten()], + [symb_to_bits[symb][prec] for symb in symbols.imag.flatten()])) + for prec in range(bits_per_real_symbol)]) + + if len(symbols.shape) > 2: + raise ValueError(f"`symbols` should be 1 or 2 dimensional but is shape {symbols.shape}") + + if symbols.ndim == 1: # If symbols shaped as vector, return as vector + bits.reshape((len(bits), )) + + return bits + +def _yF_to_hJ(y, F, modulation="BPSK"): + """Convert :math:`O(v) = ||y - F v||^2` to modulated quadratic form. + + Constructs coefficients for the form + :math:`O(v) = v^{\dagger} J v - 2 \Re(h^{\dagger} v) + \\text{offset}`. + + Args: + y: Received symbols as a NumPy column vector of complex or real values. + + F: Wireless channel as an :math:`i \\times j` NumPy matrix of complex + values, where :math:`i` rows correspond to :math:`y_i` receivers + and :math:`j` columns correspond to :math:`v_i` transmitted symbols. + + modulation: Modulation. Supported values are the default non-quadrature + modulation, 'BPSK', and quadrature modulations 'QPSK', '16QAM', + '64QAM', and '256QAM'. + + Returns: + Three tuple of amplitude-modulated linear biases :math:`h`, as a NumPy + vector, amplitude-modulated quadratic interactions, :math:`J`, as a + matrix, and offset as a real scalar. + """ + offset, h, J = _quadratic_form(y, F) # Conversion to quadratic form + + # Separate real and imaginary parts of quadratic form: + h, J = _real_quadratic_form(h, J, modulation) + + # Amplitude-modulate the biases in the quadratic form: + h, J = _amplitude_modulated_quadratic_form(h, J, modulation) + + return h, J, offset + +def _bits_to_symbols(bits, + modulation="BPSK", + num_transmitters=None): + """Convert bits to modulated symbols. + + Args: + bits: Bits as a NumPy array. + + modulation: Modulation. Supported values are the default non-quadrature + modulation, 'BPSK', and quadrature modulations 'QPSK', '16QAM', + '64QAM', and '256QAM'. + + Returns: + Transmitted symbols as a NumPy vector. + """ + if modulation not in mod_config: + raise ValueError(f"Unsupported modulation: {modulation}") + + num_bits = len(bits) + + if num_transmitters is None: + num_transmitters = num_bits // mod_config[modulation].transmitters_per_bit + + if num_transmitters == num_bits: + symbols = bits + else: + num_amps, rem = divmod(len(bits), (2*num_transmitters)) + if num_amps > 64: + raise ValueError('Complex encoding is limited to 64 bits in' + 'real and imaginary parts; `num_transmitters` is' + 'too small') + if rem != 0: + raise ValueError('number of bits must be divisible by `num_transmitters` ' + 'for modulation schemes') + + bitsR = np.reshape(bits, (num_amps, 2 * num_transmitters)) + amps = 2 ** np.arange(0, num_amps)[:, np.newaxis] + + symbols = np.sum(amps*bitsR[:, :num_transmitters], axis=0) \ + + 1j*np.sum(amps*bitsR[:, num_transmitters:], axis=0) + + return symbols + +def _lattice_to_attenuation_matrix(lattice, + transmitters_per_node=1, + receivers_per_node=1, + neighbor_root_attenuation=1): + """Generate an attenuation matrix from a given lattice. + + The attenuation matrix, a NumPy :class:`~numpy.ndarray` matrix with a row + for each receiver and column for each transmitter, specifies the expected + root-power of transmission between integer-indexed transmitters and + receivers. + + It sets uniform transmission of power for on-site transmitter-receiver pairs + and uniform transmission from transmitters to receivers up to graph distance + 1. + """ + # Developer note: this code should exploit sparsity, and track the label map. + # The code could be generalized to account for asymmetric transmission patterns, + # or real-valued spatial structure. + + if any('num_transmitters' in lattice.nodes[n] for n in lattice.nodes) or \ + any('num_receivers' in lattice.nodes[n] for n in lattice.nodes): + + node_to_transmitters = {} #Integer labels of transmitters at node + node_to_receivers = {} #Integer labels of receivers at node + t_ind = 0 + r_ind = 0 + for n in lattice.nodes: + num = transmitters_per_node + if 'num_transmitters' in lattice.nodes[n]: + num = lattice.nodes[n]['num_transmitters'] + node_to_transmitters[n] = list(range(t_ind, t_ind + num)) + t_ind = t_ind + num + + num = receivers_per_node + if 'num_receivers' in lattice.nodes[n]: + num = lattice.nodes[n]['num_receivers'] + node_to_receivers[n] = list(range(r_ind, r_ind + num)) + r_ind = r_ind + num + + A = np.zeros(shape=(r_ind, t_ind)) + for n0 in lattice.nodes: + root_receivers = node_to_receivers[n0] + for r in root_receivers: + for t in node_to_transmitters[n0]: + A[r, t] = 1 + for neigh in lattice.neighbors(n0): + for t in node_to_transmitters[neigh]: + A[r, t] = neighbor_root_attenuation + else: + A = np.identity(lattice.number_of_nodes()) + # Uniform case: + node_to_int = {n:idx for idx, n in enumerate(lattice.nodes())} + for n0 in lattice.nodes: + root = node_to_int[n0] + for neigh in lattice.neighbors(n0): + A[node_to_int[neigh], root] = neighbor_root_attenuation + A = np.tile(A, (receivers_per_node, transmitters_per_node)) + + node_to_receivers = {n: [node_to_int[n] + i*len(node_to_int) for + i in range(receivers_per_node)] for n in node_to_int} + + node_to_transmitters = {n: [node_to_int[n] + i*len(node_to_int) for + i in range(transmitters_per_node)] for n in node_to_int} + + return A, node_to_transmitters, node_to_receivers + +def create_channel(num_receivers: int = 1, + num_transmitters: int = 1, + F_distribution: Optional[Tuple[str, str]] = None, + random_state: Optional[Union[int, np.random.mtrand.RandomState]] = None, + attenuation_matrix: Optional[np.ndarray] = None) -> Tuple[np.ndarray, float]: + """Create a channel model. + + Args: + num_receivers: Number of receivers. + + num_transmitters: Number of transmitters. + + F_distribution: Distribution for the channel. Supported values are: + + * First value: ``normal`` and ``binary``. + * Second value: ``real`` and ``complex``. + + random_state: Seed for a random state or a random state. + + attenuation_matrix: Root of the power associated with each + transmitter-receiver channel. + + Returns: + Two-tuple of channel and channel power, where the channel is an + :math:`i \times j` matrix with :math:`i` rows corresponding to the + receivers and :math:`j` columns to the transmitters, and channel power + is a number. + """ + + if num_receivers < 1 or num_transmitters < 1: + raise ValueError('At least one receiver and one transmitter are required.') + + channel_power = num_transmitters + + random_state = np.random.default_rng(random_state) + + if F_distribution is None: + F_distribution = ('normal', 'complex') + elif type(F_distribution) is not tuple or len(F_distribution) != 2: + raise ValueError('F_distribution should be a tuple of strings or None') + + if F_distribution[0] == 'normal': + if F_distribution[1] == 'real': + F = random_state.normal(0, 1, size=(num_receivers, num_transmitters)) + else: + channel_power = 2 * num_transmitters + F = random_state.normal(0, 1, size=(num_receivers, num_transmitters)) + \ + 1j*random_state.normal(0, 1, size=(num_receivers, num_transmitters)) + elif F_distribution[0] == 'binary': + if F_distribution[1] == 'real': + F = (1 - 2*random_state.integers(2, size=(num_receivers, num_transmitters))) + else: + channel_power = 2*num_transmitters #For integer precision purposes: + F = ((1 - 2*random_state.integers(2, size=(num_receivers, num_transmitters))) + + 1j*(1 - 2*random_state.integers(2, size=(num_receivers, num_transmitters)))) + + if attenuation_matrix is not None: + if np.iscomplex(attenuation_matrix).any(): + raise ValueError('attenuation_matrix must not have complex values') + F = F * attenuation_matrix #Dense format for now, this is slow. + channel_power *= np.mean(np.sum(attenuation_matrix * attenuation_matrix, + axis=0)) / num_receivers + + return F, channel_power + +def _constellation_properties(modulation): + """Return bits per symbol, amplitudes, and mean power for QAM constellation. + + Constellation mean power makes the standard assumption that symbols are + sampled uniformly at random for the signal. + """ + if modulation not in mod_config: + raise ValueError('Unsupported modulation method') + + amps = 1 + 2*np.arange(mod_config[modulation].amplitudes) + constellation_mean_power = 1 if modulation == 'BPSK' else 2 * np.mean(amps*amps) + + return mod_config[modulation].bits_per_transmitter, amps, constellation_mean_power + +def _create_transmitted_symbols(num_transmitters, + amps=[-1, 1], + quadrature=True, + random_state=None): + """Generate symbols. + + Symbols are generated uniformly at random as a function of the quadrature + and amplitude modulation. + + The power per symbol is not normalized, it is proportional to :math:`N_t*sig2`, + where :math:`sig2 = [1, 2, 10, 42]` for BPSK, QPSK, 16QAM and 64QAM respectively. + + The complex and real-valued parts of all constellations are integer. + + Args: + num_transmitters: Number of transmitters. + + amps: Amplitudes as an iterable. + + quadrature: Quadrature (True) or only phase-shift keying such as BPSK (False). + + random_state: Seed for a random state or a random state. + + Returns: + + Symbols, as a column vector of length ``num_transmitters``. + + """ + + if np.iscomplex(amps).any(): + raise ValueError('Amplitudes cannot have complex values') + + if any(np.modf(amps)[0]): + raise ValueError('Amplitudes must have integer values') + + random_state = np.random.default_rng(random_state) + + if not quadrature: + transmitted_symbols = random_state.choice(amps, size=(num_transmitters, 1)) + else: + transmitted_symbols = random_state.choice(amps, size=(num_transmitters, 1)) \ + + 1j*random_state.choice(amps, size=(num_transmitters, 1)) + + return transmitted_symbols + +def _create_signal(F, + transmitted_symbols=None, + channel_noise=None, + SNRb=float('Inf'), + modulation='BPSK', + channel_power=None, + random_state=None): + """Simulate a transmission signal. + + Generates random transmitted symbols and optionally noise, math:`y = F v + n`, + where the channel, :math:`F`, is assumed to consist of independent and + identically distributed (i.i.d) elements such that + :math:`F\dagger*F = N_r I[N_t]*cp` where :math:`I` is the identity matrix + and :math:`cp` the channel power; the transmitted symbols, :math:`v` or + ``transmitted_symbols``, are assumed to consist of i.i.d unscaled + constellations elements (integer valued in real and complex parts), and + mean constellation power dictates a rescaling relative to + :math:`E[v v\dagger] = I[Nt]`. + + ``channel_noise`` is assumed, or created, to be suitably scaled: + :math:`N_0 I[Nt] = SNRb `. + + Args: + F: Wireless channel as an :math:`i \times j` matrix of complex values, + where :math:`i` rows correspond to :math:`y_i` receivers and :math:`j` + columns correspond to :math:`v_i` transmitted symbols. + + transmitted_symbols: Transmitted symbols as a column vector. + + channel_noise: Channel noise as a column vector. + + SNRb: Signal-to-noise ratio. + + modulation: Modulation. Supported values are 'BPSK', 'QPSK', '16QAM', + '64QAM', and '256QAM'. + + channel_power: Channel power. By default, equal to the number + of transmitters. + + random_state: Seed for a random state or a random state. + + Returns: + Four-tuple of received signals (``y``), transmitted symbols (``v``), + channel noise, and random_state, where ``y`` is a column vector of length + equal to the rows of ``F``. + """ + + num_receivers = F.shape[0] + num_transmitters = F.shape[1] + + random_state = np.random.default_rng(random_state) + + bits_per_transmitter, amps, constellation_mean_power = _constellation_properties(modulation) + + if transmitted_symbols is not None: + if modulation == 'BPSK' and np.iscomplex(transmitted_symbols).any(): + raise ValueError(f"BPSK transmitted signals must be real") + if modulation != 'BPSK' and np.isreal(transmitted_symbols).any(): + raise ValueError(f"Quadrature transmitted signals must be complex") + else: + quadrature = modulation != 'BPSK' + transmitted_symbols = _create_transmitted_symbols( + num_transmitters, amps=amps, quadrature=quadrature, random_state=random_state) + + if SNRb <= 0: + raise ValueError(f"signal-to-noise ratio must be positive. SNRb={SNRb}") + + if SNRb == float('Inf'): + y = np.matmul(F, transmitted_symbols) + elif channel_noise is not None: + y = channel_noise + np.matmul(F, transmitted_symbols) + else: + # Energy_per_bit: + if channel_power == None: + #Assume proportional to num_transmitters; i.e., every channel + # component is RMSE 1 and 1 bit + channel_power = num_transmitters + + #Eb is the same for QPSK and BPSK + # Eb/N0 = SNRb (N0 = 2 sigma^2, the one-sided PSD ~ kB T at antenna) + # SNRb and Eb, together imply N0 + Eb = channel_power * constellation_mean_power / bits_per_transmitter + N0 = Eb / SNRb + # Noise is complex by definition, hence 1/2 power in real & complex parts + sigma = np.sqrt(N0 / 2) + + # Channel noise of covariance N0*I_{NR}. Noise is complex by definition, + # although for real channel and symbols we need only worry about + # real part: + channel_noise = sigma*(random_state.normal(0, 1, size=(num_receivers, 1)) \ + + 1j*random_state.normal(0, 1, size=(num_receivers, 1))) + if modulation == 'BPSK' and np.isreal(F).all(): + channel_noise = channel_noise.real + + y = channel_noise + np.matmul(F, transmitted_symbols) + + return y, transmitted_symbols, channel_noise, random_state + +def mimo(modulation: Literal["BPSK", "QPSK", "16QAM", "64QAM", "256QAM"] = "BPSK", + y: Union[np.array, None] = None, + F: Union[np.array, None] = None, + *, + transmitted_symbols: Union[np.array, None] = None, + channel_noise: Union[np.array, None] = None, + num_transmitters: int = None, + num_receivers: int = None, + SNRb: float = float('Inf'), + seed: Union[None, int, np.random.RandomState] = None, + F_distribution: Union[None, tuple] = None, + attenuation_matrix = None) -> dimod.BinaryQuadraticModel: + """Generate a multi-input multiple-output (MIMO) channel-decoding problem. + + In radio networks, `MIMO `_ is a method + of increasing link capacity by using multiple transmission and receiving + antennas to exploit multipath propagation. + + Users transmit complex-valued symbols over a random channel, :math:`F`, + subject to additive white Gaussian noise. Given the received signal, + :math:`y`, the log likelihood of a given symbol set, :math:`v`, is + :math:`MLE = argmin || y - F v ||_2`. When :math:`v` is encoded as + a linear sum of bits, the optimization problem is a binary quadratic model. + + Depending on its parameters, this function can model code division multiple + access (CDMA) _[#T02, #R20], 5G communication networks _[#Prince], or + other problems. + + Args: + modulation: Constellation (symbol set) users can transmit. Symbols are + assumed to be transmitted with equal probability. Supported values + are: + + * 'BPSK' + + Binary Phase Shift Keying. Transmitted symbols are :math:`+1, -1`; + no encoding is required. A real-valued channel is assumed. + + * 'QPSK' + + Quadrature Phase Shift Keying. Transmitted symbols are + :math:`1+1j, 1-1j, -1+1j, -1-1j` normalized by + :math:`\\frac{1}{\\sqrt{2}}`. Bits are encoded as a real vector + concatenated with an imaginary vector. + + * '16QAM' + + Each user is assumed to select independently from 16 symbols. + The transmitted symbol is a complex value that can be encoded + by two bits in the imaginary part and two bits in the real + part. Highest precision real and imaginary bit vectors are + concatenated to lower precision bit vectors. + + * '64QAM' + + A QPSK symbol set is generated and symbols are further amplitude + modulated by an independently and uniformly distributed random + amount from :math:`[1, 3]`. + + * '256QAM' + + A QPSK symbol set is generated and symbols are further amplitude + modulated by an independently and uniformly distributed random + amount from :math:`[1, 3, 5]`. + + y: Complex- or real-valued received signal, as a NumPy array. If + ``None``, generated from other arguments. + + F: Complex- or real-valued channel, as a NumPy array. If ``None``, + generated from other arguments. Note that for correct interpretation + of SNRb, channel power should be normalized to ``num_transmitters``. + + num_transmitters: Number of users. Each user transmits one symbol per + frame. + + num_receivers: Number of receivers of a channel. Must be consistent + with the length of any provided signal, ``len(y)``. + + SNRb: Signal-to-noise ratio per bit used to generate the noisy signal + when ``y`` is not provided. If ``float('Inf')``, no noise is + added. + + :math:`SNR_b = E_b/N_0`, where :math:`E_b` is energy per bit, + and :math:`N_0` is the one-sided power-spectral density. :math:`N_0` + is typically :math:`k_B T` at the receiver. To convert units of + :math:`dB` to :math:`SNR_b` use :math:`SNRb=10^{SNR_b[decibels]/10}`. + + seed: Random seed, as an integer, or state, as a + :class:`numpy.random.RandomState` instance. + + transmitted_symbols: + Set of symbols transmitted. Used in combination with ``F`` to + generate the received signal, :math:`y`. The number of transmitted + symbols must be consistent with ``F``. + + For BPSK and QPSK modulations, statistics of the ensemble do not + depend on the choice: all choices are equivalent. By default, + symbols are chosen for all users as :math:`1` or :math:`1 + 1j`, + respectively. Note that for correct analysis by some solvers, applying + spin-reversal transforms may be necessary. + + For QAM modulations such as 16QAM, amplitude randomness affects + likelihood in a non-trivial way. By default, symbols are chosen + independently and identically distributed from the constellations. + + channel_noise: Channel noise as a NumPy array of complex values. Must + be consistent with the number of receivers. + + F_distribution: + Zero-mean, variance-one distribution, in tuple form + ``(distribution, type)``, used to generate each element in ``F`` + when ``F`` is not provided . Supported values are: + + * ``'normal'`` or ``'binary'`` for the distribution + * ``'real'`` or ``'complex'`` for the type + + For large numbers of receivers and transmitters, statistical + properties of the likelihood are weakly dependent on the + distribution. Choosing ``'binary'`` allows for integer-valued + Hamiltonians while ``'normal'`` is a more typical model. The channel + can be real or complex; in many cases this represents a superficial + distinction up to rescaling. For real-valued symbols (BPSK) the + default is ``('normal', 'real')``; otherwise, the default is + ``('normal', 'complex')``. + + attenuation_matrix: + Root of the power associated with each transmitter-receiver channel; + use for sparse and structured codes. + + Returns: + Binary quadratic model defining the log-likelihood function. + + Example: + + This example generates an instance of a CDMA problem in the high-load + regime, near a first-order phase transition. + + >>> num_transmitters = 64 + >>> transmitters_per_receiver = 1.5 + >>> SNRb = 5 + >>> bqm = dimod.generators.mimo(modulation='BPSK', + ... num_transmitters = 64, + ... num_receivers = round(num_transmitters / transmitters_per_receiver), + ... SNRb=SNRb, + ... F_distribution = ('binary', 'real')) + + .. [#T02] T. Tanaka IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 48, NO. 11, NOVEMBER 2002 + + .. [#R20] J. Raymond, N. Ndiaye, G. Rayaprolu and A. D. King, + "Improving performance of logical qubits by parameter tuning and topology compensation," + 2020 IEEE International Conference on Quantum Computing and Engineering (QCE), + Denver, CO, USA, 2020, pp. 295-305, doi: 10.1109/QCE49297.2020.00044. + + .. [#Prince] Various (https://paws.princeton.edu/) + """ + + random_state = np.random.default_rng(seed) + + if y is not None: + if len(y.shape) == 1: + y = y.reshape(y.shape[0], 1) + elif len(y.shape) != 2 or y.shape[1] != 1: + raise ValueError(f"y should have shape (n, 1) or (n,) for some n; given: {y.shape}") + + if F is None: + + if num_transmitters: + if num_transmitters <= 0: + raise ValueError('Configured number of transmitters must be positive') + elif transmitted_symbols: + num_transmitters = len(transmitted_symbols) + else: + ValueError('`num_transmitters` is not specified and cannot' + 'be inferred from `F` or `transmitted_symbols` (both None)') + + if num_receivers: + if num_receivers <= 0: + raise ValueError('Configured number of receivers must be positive') + elif y: + num_receivers = y.shape[0] + elif channel_noise: + num_receivers = channel_noise.shape[0] + else: + raise ValueError('`num_receivers` is not specified and cannot' + 'be inferred from `F`, `y` or `channel_noise` (all None)') + + F, channel_power = create_channel(num_receivers=num_receivers, + num_transmitters=num_transmitters, + F_distribution=F_distribution, + random_state=random_state, + attenuation_matrix=attenuation_matrix) + else: + channel_power = num_transmitters + + if y is None: + y, _, _, _ = _create_signal(F, + transmitted_symbols=transmitted_symbols, + channel_noise=channel_noise, + SNRb=SNRb, + modulation=modulation, + channel_power=channel_power, + random_state=random_state) + + h, J, offset = _yF_to_hJ(y, F, modulation) + + return dimod.BQM(h[:, 0], J, 'SPIN', offset=offset) + +def coordinated_multipoint(lattice: 'networkx.Graph', + modulation: Literal["BPSK", "QPSK", "16QAM", "64QAM", "256QAM"] = "BPSK", + y: Optional[np.array] = None, + F: Union[np.array, None] = None, + *, + transmitted_symbols: Union[np.array, None] = None, + channel_noise: Union[np.array, None] = None, + SNRb: float = float('Inf'), + seed: Union[None, int, np.random.RandomState] = None, + F_distribution: Union[None, str] = None) -> dimod.BinaryQuadraticModel: + """Generate a coordinated multi-point (CoMP) decoding problem. + + In `coordinated multipoint (CoMP) `_ + neighboring cellular base stations coordinate transmissions and jointly + process received signals. + + Users transmit complex-valued symbols over a random channel, :math:`F`, + subject to additive white Gaussian noise. Given the received signal, + :math:`y`, the log likelihood of a given symbol set, :math:`v`, is + :math:`MLE = argmin || y - F v ||_2`. When :math:`v` is encoded as + a linear sum of bits, the optimization problem is a binary quadratic model. + + Args: + lattice: Geometry, as a :class:`networkx.Graph`, defining + the set of nearest-neighbor base stations. + + Each base station has ``num_receivers`` receivers and + ``num_transmitters`` local transmitters, set as either attributes + of the graph or as per-node values. Transmitters from neighboring + base stations are also received. + + modulation: Constellation (symbol set) users can transmit. Symbols are + assumed to be transmitted with equal probability. Supported values + are: + + * 'BPSK' + Binary Phase Shift Keying. Transmitted symbols are :math:`+1, -1`; + no encoding is required. A real-valued channel is assumed. + * 'QPSK' + Quadrature Phase Shift Keying. Transmitted symbols are + :math:`1+1j, 1-1j, -1+1j, -1-1j` normalized by + :math:`\\frac{1}{\\sqrt{2}}`. Bits are encoded as a real vector + concatenated with an imaginary vector. + * '16QAM' + Each user is assumed to select independently from 16 symbols. + The transmitted symbol is a complex value that can be encoded + by two bits in the imaginary part and two bits in the real + part. Highest precision real and imaginary bit vectors are + concatenated to lower precision bit vectors. + * '64QAM' + A QPSK symbol set is generated and symbols are further amplitude + modulated by an independently and uniformly distributed random + amount from :math:`[1, 3]`. + * '256QAM' + A QPSK symbol set is generated and symbols are further amplitude + modulated by an independently and uniformly distributed random + amount from :math:`[1, 3, 5]`. + + y: Complex- or real-valued received signal, as a NumPy array. If + ``None``, generated from other arguments. + + F: Transmission channel. Currently not supported and must be ``None``. + + transmitted_symbols: Set of symbols transmitted. Used in combination + with ``F`` to generate the received signal, :math:`y`. The number + of transmitted symbols must be consistent with ``F``. + + For BPSK and QPSK modulations, statistics of the ensemble do not + depend on the choice: all choices are equivalent. By default, + symbols are chosen for all users as :math:`1` or :math:`1 + 1j`, + respectively. Note that for correct analysis by some solvers, + applying spin-reversal transforms may be necessary. + + For QAM modulations such as 16QAM, amplitude randomness affects + likelihood in a non-trivial way. By default, symbols are chosen + independently and identically distributed from the constellations. + + channel_noise: Channel noise as a complex value. + + SNRb: Signal-to-noise ratio per bit, :math:`SNRb=10^{SNR_b[decibels]/10}`, + used to generate the noisy signal when ``y`` is not provided. + If ``float('Inf')``, no noise is added. + + seed: Random seed, as an integer, or state, as a + :class:`numpy.random.RandomState` instance. + + F_distribution: Zero-mean, variance-one distribution, in tuple form + ``(distribution, type)``, used to generate each element in ``F`` + when ``F`` is not provided . Supported values are: + + * ``'normal'`` or ``'binary'`` for the distribution + * ``'real'`` or ``'complex'`` for the type + + For large numbers of receivers and transmitters, statistical + properties of the likelihood are weakly dependent on the + distribution. Choosing ``'binary'`` allows for integer-valued + Hamiltonians while ``'normal'`` is a more typical model. The channel + can be real or complex; in many cases this represents a superficial + distinction up to rescaling. For real-valued symbols (BPSK) the + default is ``('normal', 'real')``; otherwise, the default is + ``('normal', 'complex')``. + + Returns: + bqm: Binary quadratic model defining the log-likelihood function. + + Example: + + Generate an instance of a CDMA problem in the high-load regime, near a + first-order phase transition: + + .. doctest:: # TODO: reconsider example/default graph + :skipif: True + + >>> import networkx as nx + >>> G = nx.complete_graph(4) + >>> nx.set_node_attributes(G, values={n:2*n+1 for n in G.nodes()}, name='num_transmitters') + >>> nx.set_node_attributes(G, values={n:2 for n in G.nodes()}, name='num_receivers') + >>> transmitted_symbols = np.random.choice([1, -1], + ... size=(sum(nx.get_node_attributes(G, "num_transmitters").values()), 1)) + >>> bqm = dimod.generators.coordinated_multipoint(G, + ... modulation='BPSK', + ... transmitted_symbols=transmitted_symbols, + ... SNRb=5, + ... F_distribution = ('binary', 'real')) + + """ + + if not hasattr(lattice, 'edges') or not hasattr(lattice, 'nodes'): # not nx.Graph: + raise ValueError('Lattice must be a :class:`networkx.Graph`') + + attenuation_matrix, _, _ = _lattice_to_attenuation_matrix(lattice) + + num_receivers, num_transmitters = attenuation_matrix.shape + + bqm = mimo( + modulation=modulation, + y=y, + F=F, + transmitted_symbols=transmitted_symbols, + channel_noise=channel_noise, + num_transmitters=num_transmitters, + num_receivers=num_receivers, + SNRb=SNRb, + seed=seed, + F_distribution=F_distribution, + attenuation_matrix=attenuation_matrix) + + return bqm + +# Linear-filter functions. These are not used for bit-encoding MIMO problems +# and are maintained here for user convenience + +def linear_filter(F, method='zero_forcing', SNRoverNt=float('Inf'), PoverNt=1): + """Construct a linear filter for estimating transmitted signals. + + Following the conventions of MacKay\ [#Mackay]_, a filter is constructed + for independent and identically distributed Gaussian noise at power spectral + density :math:`N_0`: + + :math:`N_0 I[N_r] = E[n n^{\dagger}]` + + For independent and identically distributed (i.i.d), zero mean, transmitted symbols, + + :math:`P_c I[N_t] = E[v v^{\dagger}]` + + where :math:`P_{c}`, the constellation's mean power, is equal to + :math:`(1, 2, 10, 42)N_t` for BPSK, QPSK, 16QAM, 64QAM respectively. + + For an i.i.d channel, + + :math:`N_r N_t = E_F[T_r[F F^{\dagger}]] \qquad \Rightarrow \qquad E[||F_{\mu, i}||^2] = 1` + + Symbols are assumed to be normalized: + + :math:`\\frac{SNR}{N_t} = \\frac{P_c}{N_0}` + + :math:`SNR_b = \\frac{SNR}{N_t B_c}` + + where :math:`B_c` is bit per symbol, equal to :math:`(1, 2, 4, 8)` + for BPSK, QPSK, 16QAM, 64QAM respectively. + + Typical use case: set :math:`\\frac{SNR}{N_t} = SNR_b`. + + .. [#Mackay] Matthew R. McKay, Iain B. Collings, Antonia M. Tulino. + "Achievable sum rate of MIMO MMSE receivers: A general analytic framework" + IEEE Transactions on Information Theory, February 2010 + arXiv:0903.0666 [cs.IT] + + Reference: + + https://www.youtube.com/watch?v=U3qjVgX2poM + """ + if method not in ['zero_forcing', 'matched_filter', 'MMSE']: + raise ValueError('Unsupported filtering method') + + if method == 'zero_forcing': + # Moore-Penrose pseudo inverse + return np.linalg.pinv(F) + + Nr, Nt = F.shape + + if method == 'matched_filter': # F = root(Nt/P) Fcompconj + return F.conj().T / np.sqrt(PoverNt) + + # method == 'MMSE': + return np.matmul( + F.conj().T, + np.linalg.pinv(np.matmul(F, F.conj().T) + np.identity(Nr)/SNRoverNt) + ) / np.sqrt(PoverNt) + +def filter_marginal_estimator(x: np.array, modulation: str): + """Map filter output to valid symbols. + + Takes the continuous filter output and maps each estimated symbol to + the nearest valid constellation value. + """ + if modulation is not None: + if modulation == 'BPSK' or modulation == 'QPSK': + max_abs = 1 + elif modulation == '16QAM': + max_abs = 3 + elif modulation == '64QAM': + max_abs = 7 + elif modulation == '128QAM': + max_abs = 15 + else: + raise ValueError('Unknown modulation') + + #Real part (nearest): + x_R = 2*np.round((x.real - 1)/2) + 1 + x_R = np.where(x_R < -max_abs, -max_abs, x_R) + x_R = np.where(x_R > max_abs, max_abs, x_R) + + if modulation != 'BPSK': + x_I = 2*np.round((x.imag - 1)/2) + 1 + x_I = np.where(x_I <- max_abs, -max_abs, x_I) + x_I = np.where(x_I > max_abs, max_abs, x_I) + return x_R + 1j*x_I + return x_R diff --git a/docs/reference/generators.rst b/docs/reference/generators.rst index 146994007..c69601d23 100644 --- a/docs/reference/generators.rst +++ b/docs/reference/generators.rst @@ -43,9 +43,11 @@ Optimization .. autosummary:: :toctree: generated/ + coordinated_multipoint independent_set maximum_independent_set maximum_weight_independent_set + mimo random_bin_packing random_knapsack random_multi_knapsack diff --git a/releasenotes/notes/feature-mimo-59c3d6d02335061c.yaml b/releasenotes/notes/feature-mimo-59c3d6d02335061c.yaml new file mode 100644 index 000000000..90df73f1c --- /dev/null +++ b/releasenotes/notes/feature-mimo-59c3d6d02335061c.yaml @@ -0,0 +1,8 @@ +--- +features: + - | + Add ``mimo()`` function for generating a multi-input + multiple-output (MIMO) channel-decoding problem. + - | + Add ``coordinated_multipoint()`` function for generating a coordinated multi-point + (CoMP) decoding problem. diff --git a/tests/test_generators.py b/tests/test_generators.py index 9b33b0069..2cedea26f 100644 --- a/tests/test_generators.py +++ b/tests/test_generators.py @@ -27,7 +27,6 @@ else: _networkx = True - class TestRandomGNMRandomBQM(unittest.TestCase): def test_bias_generator(self): def gen(n): @@ -1003,7 +1002,6 @@ def integers(self, *args, **kwargs): for trial in range(10): with self.subTest(trial=trial): bqm = dimod.generators.random_2in4sat(4, 1, seed=seed) - # in the ground state they should not be all equal ss = dimod.ExactSolver().sample(bqm) self.assertEqual(sum(ss.first.sample.values()), 0) @@ -1129,3 +1127,566 @@ def test_constraints_squares(self): self.assertEqual(term, -2) else: self.assertEqual(term, 24) + +@unittest.skipUnless(_networkx, "no networkx installed") +class TestMIMO(unittest.TestCase): + + def setUp(self): + self.rng = np.random.default_rng(1) + self.symbols_bpsk = np.asarray([[-1, 1]]) + self.symbols_qam = lambda a: np.array([[complex(i, j)] \ + for i in range(-a, a + 1, 2) for j in range(-a, a + 1, 2)]) + + def _make_honeycomb(L: int): + """Generate 2L by 2L triangular lattice. + + The generated lattice has open boundaries and cut corners to make a hexagon. + + Args: + L: Length of lattice. + + Returns: + :class:`networkx.Graph`. + """ + G = nx.Graph() + + G.add_edges_from([((x, y), (x, y + 1)) for x in range(2*L + 1) for y in range(2*L)]) + G.add_edges_from([((x, y), (x + 1, y)) for x in range(2*L) for y in range(2*L + 1)]) + G.add_edges_from([((x, y), (x + 1, y +1 )) for x in range(2*L) for y in range(2*L)]) + G.remove_nodes_from([(i, j) for j in range(L) for i in range(L + 1 + j, 2*L + 1) ]) + G.remove_nodes_from([(i, j) for i in range(L) for j in range(L + 1 + i, 2*L + 1)]) + + return G + + self._make_honeycomb = lambda L: _make_honeycomb(L) + + def _effective_fields(self, bqm): + num_var = bqm.num_variables + effFields = np.zeros(num_var) + for key in bqm.quadratic: + effFields[key[0]] += bqm.adj[key[0]][key[1]] + effFields[key[1]] += bqm.adj[key[0]][key[1]] + for key in bqm.linear: + effFields[key] += bqm.linear[key] + return effFields + + def test_filter_marginal_estimators(self): + # Tested but so far this function is unused + fme = dimod.generators.wireless.filter_marginal_estimator + + filtered_signal = self.rng.random(20) + np.arange(-20, 20, 2) + estimated_source = fme(filtered_signal, 'BPSK') + self.assertTrue(0 == len(set(estimated_source).difference(np.arange(-1, 3, 2)))) + self.assertTrue(np.all(estimated_source[:-1] <= estimated_source[1:])) + + filtered_signal = filtered_signal + 1j*(-self.rng.random(20) + np.arange(20, -20, -2)) + + for modulation in ['QPSK','16QAM','64QAM']: + estimated_source = fme(filtered_signal, modulation=modulation) + self.assertTrue(np.all(np.flip(estimated_source.real) == estimated_source.imag)) + + def test_linear_filter(self): + + Nt = 5 + Nr = 7 + # linear_filter(F, method='zero_forcing', PoverNt=1, SNRoverNt = 1) + F = self.rng.normal(size=(Nr,Nt)) + 1j*self.rng.normal(size=(Nr,Nt)) + Fsimple = np.identity(Nt) # Nt=Nr + + #BPSK, real channel: + transmitted_symbolsQAM = dimod.generators.wireless._create_transmitted_symbols(Nt, + amps=[-3, -1, 1, 3], quadrature=True) + + y = np.matmul(F, transmitted_symbolsQAM) + + # Defaults + W = dimod.generators.wireless.linear_filter(F=F) + self.assertEqual(W.shape,(Nt,Nr)) + + # Check arguments: + W = dimod.generators.wireless.linear_filter(F=F, + method='matched_filter', PoverNt=0.5, SNRoverNt=1.2) + self.assertEqual(W.shape,(Nt,Nr)) + + # Over-constrained noiseless channel by default, zero_forcing and MMSE are perfect: + for method in ['zero_forcing', 'MMSE']: + W = dimod.generators.wireless.linear_filter(F=F, method=method) + reconstructed_symbols = np.matmul(W,y) + self.assertTrue(np.all(np.abs(reconstructed_symbols - transmitted_symbolsQAM) < 1e-8)) + + # matched_filter and MMSE (non-zero noise) are erroneous given interfered signal: + W = dimod.generators.wireless.linear_filter(F=F, method='MMSE', PoverNt=0.5, SNRoverNt=1) + reconstructed_symbols = np.matmul(W, y) + self.assertTrue(np.all(np.abs(reconstructed_symbols - transmitted_symbolsQAM) > 1e-8)) + + def test_quadratic_forms(self): + # Quadratic form must evaluate to match original objective: + num_var = 3 + num_receivers = 5 + F = self.rng.normal(0, 1, size=(num_receivers, num_var)) + \ + 1j*self.rng.normal(0, 1, size=(num_receivers, num_var)) + y = self.rng.normal(0, 1, size=(num_receivers, 1)) + \ + 1j*self.rng.normal(0, 1, size=(num_receivers, 1)) + + # Random test case: + vUnwrap = self.rng.normal(0, 1, size=(2*num_var, 1)) + v = vUnwrap[:num_var, :] + 1j*vUnwrap[num_var:, :] + vec = y - np.matmul(F, v) + val1 = np.matmul(vec.T.conj(), vec) + + # Check complex quadratic form + k, h, J = dimod.generators.wireless._quadratic_form(y, F) + val2 = np.matmul(v.T.conj(), np.matmul(J, v)) + (np.matmul(h.T.conj(), v)).real + k + self.assertLess(abs(val1 - val2), 1e-8) + + # Check unwrapped complex quadratic form: + h, J = dimod.generators.wireless._real_quadratic_form(h, J) + val3 = np.matmul(vUnwrap.T, np.matmul(J, vUnwrap)) + np.matmul(h.T, vUnwrap) + k + self.assertLess(abs(val1 - val3), 1e-8) + + # Check zero energy for y generated from F: + y = np.matmul(F, v) + k, h, J = dimod.generators.wireless._quadratic_form(y, F) + val2 = np.matmul(v.T.conj(), np.matmul(J, v)) + (np.matmul(h.T.conj(), v)).real + k + self.assertLess(abs(val2), 1e-8) + h, J = dimod.generators.wireless._real_quadratic_form(h, J) + val3 = np.matmul(vUnwrap.T, np.matmul(J, vUnwrap)) + np.matmul(h.T, vUnwrap) + k + self.assertLess(abs(val3), 1e-8) + + def test_real_quadratic_form(self): + h_in, J_in = np.array([1, 1]), np.array([2]) + h, J = dimod.generators.wireless._real_quadratic_form(h_in, J_in) + self.assertTrue(np.array_equal(h_in, h)) + self.assertTrue(np.array_equal(J_in, J)) + + h_in, J_in = np.array([1+1j, 1-1j]), np.array([[0, 2], [0, 0]]) + h, J = dimod.generators.wireless._real_quadratic_form(h_in, J_in) + self.assertTrue(len(h) == 2*len(h_in)) + self.assertTrue(J.shape == (4, 4)) + + h, J = dimod.generators.wireless._real_quadratic_form(h_in, J_in, 'BPSK') + self.assertTrue(np.array_equal(np.real(h_in), h)) + self.assertTrue(np.array_equal(J_in, J)) + + def test_amplitude_modulated_quadratic_form(self): + num_var = 3 + h = self.rng.random(size=(num_var, 1)) + J = self.rng.random(size=(num_var, num_var)) + mods = ['BPSK', 'QPSK', '16QAM', '64QAM'] + mod_pref = [1, 1, 2, 3] + for offset in [0]: + for modI, modulation in enumerate(mods): + hO, JO = dimod.generators.wireless._amplitude_modulated_quadratic_form(h, + J, modulation=modulation) + self.assertEqual(hO.shape[0], num_var*mod_pref[modI]) + self.assertEqual(JO.shape[0], hO.shape[0]) + self.assertEqual(JO.shape[0], JO.shape[1]) + + max_val = 2**mod_pref[modI]-1 + self.assertLess(abs(max_val*np.sum(h)-np.sum(hO)), 1e-8) + self.assertLess(abs(max_val*max_val*np.sum(J)-np.sum(JO)), 1e-8) + #self.assertEqual(h.shape[0], num_var*mod_pref[modI]) + #self.assertLess(abs(bqm.offset-np.sum(np.diag(J))), 1e-8) + + def test_yF_to_hJ(self): + F = np.array([[0, 1], [1, 1]]) + + y = np.ones(2) + h_bpsk, J_bpsk, o_bpsk = dimod.generators.wireless._yF_to_hJ(y, F, 'BPSK') + h_qpsk, J_qpsk, o_qpsk = dimod.generators.wireless._yF_to_hJ(y, F, 'QPSK') + self.assertTrue(np.array_equal(h_bpsk, h_qpsk)) + self.assertTrue(np.array_equal(J_bpsk, J_qpsk)) + self.assertTrue(np.array_equal(o_bpsk, o_qpsk)) + + y = np.array([1, -1+1j]) + h_bpsk, J_bpsk, o_bpsk = dimod.generators.wireless._yF_to_hJ(y, F, 'BPSK') + h_qpsk, J_qpsk, o_qpsk = dimod.generators.wireless._yF_to_hJ(y, F, 'QPSK') + h_16, J_16, o_16 = dimod.generators.wireless._yF_to_hJ(y, F, '16QAM') + self.assertFalse(np.array_equal(h_bpsk, h_qpsk)) + self.assertFalse(np.array_equal(J_bpsk, J_qpsk)) + self.assertFalse(np.array_equal(h_qpsk, h_16)) + self.assertFalse(np.array_equal(J_qpsk, J_16)) + self.assertTrue(np.array_equal(h_bpsk, h_qpsk[:, :2])) + self.assertTrue(np.array_equal(J_bpsk, J_qpsk[:2, :2])) + self.assertTrue(np.array_equal(o_bpsk, o_qpsk)) + self.assertTrue(np.array_equal(h_qpsk, h_16[:1, :])) + self.assertTrue(np.array_equal(J_qpsk, J_16[:4, :4])) + self.assertTrue(np.array_equal(o_qpsk, o_16)) + + def test_bits_to_symbols(self): + symbols = dimod.generators.wireless._bits_to_symbols(self.symbols_bpsk, + modulation='BPSK') + self.assertTrue(np.array_equal(self.symbols_bpsk, symbols)) + + symbols = dimod.generators.wireless._bits_to_symbols(self.symbols_bpsk, + modulation='BPSK', num_transmitters=1) + self.assertTrue(np.array_equal(self.symbols_bpsk, symbols)) + + symbols = dimod.generators.wireless._bits_to_symbols(self.symbols_qam(1), + modulation='QPSK') + self.assertEqual(len(symbols), 2) + + symbols = dimod.generators.wireless._bits_to_symbols(self.symbols_qam(1), + modulation='16QAM') + self.assertEqual(len(symbols), 1) + + with self.assertRaises(ValueError): + bits = dimod.generators.wireless._bits_to_symbols(self.symbols_qam(1), + modulation='QPSK', num_transmitters=3) + + + def test_symbols_to_bits(self): + # Standard symbol cases (2D input): + bits = dimod.generators.wireless._symbols_to_bits(self.symbols_bpsk, + modulation='BPSK') + self.assertEqual(bits.sum(), 0) + self.assertTrue(bits.ndim, 2) + + bits = dimod.generators.wireless._symbols_to_bits(self.symbols_qam(1), + modulation='QPSK') + self.assertEqual(bits[:len(bits//2)].sum(), 0) + self.assertEqual(bits[len(bits//2):].sum(), 0) + self.assertTrue(bits.ndim, 2) + + bits = dimod.generators.wireless._symbols_to_bits(self.symbols_qam(3), + modulation='16QAM') + self.assertEqual(bits[:len(bits//2)].sum(), 0) + self.assertEqual(bits[len(bits//2):].sum(), 0) + + bits = dimod.generators.wireless._symbols_to_bits(self.symbols_qam(5), + modulation='64QAM') + self.assertEqual(bits[:len(bits//2)].sum(), 0) + self.assertEqual(bits[len(bits//2):].sum(), 0) + + # Standard symbol cases (1D input): + bits = dimod.generators.wireless._symbols_to_bits( + self.symbols_qam(1).reshape(4,), + modulation='QPSK') + self.assertTrue(bits.ndim, 1) + self.assertEqual(bits[:len(bits//2)].sum(), 0) + self.assertEqual(bits[len(bits//2):].sum(), 0) + + # Unsupported input + with self.assertRaises(ValueError): + bits = dimod.generators.wireless._symbols_to_bits(self.symbols_bpsk, + modulation='unsupported') + + def test_BPSK_symbol_coding(self): + #This is simply read in read out. + num_bits = 5 + bits = self.rng.choice([-1, 1], size=num_bits) + symbols = dimod.generators.wireless._bits_to_symbols(bits=bits, modulation='BPSK') + self.assertTrue(np.all(bits == symbols)) + bits = dimod.generators.wireless._symbols_to_bits(symbols=bits, modulation='BPSK') + self.assertTrue(np.all(bits == symbols)) + + def test_constellation_properties(self): + _cp = dimod.generators.wireless._constellation_properties + self.assertEqual(_cp("QPSK")[0], 2) + self.assertEqual(sum(_cp("16QAM")[1]), 4) + self.assertEqual(_cp("64QAM")[2], 42.0) + with self.assertRaises(ValueError): + bits_per_transmitter, amps, constellation_mean_power = _cp("dummy") + + def test_create_transmitted_symbols(self): + _cts = dimod.generators.wireless._create_transmitted_symbols + self.assertTrue(_cts(1, amps=[-1, 1], quadrature=False)[0][0] in [-1, 1]) + self.assertTrue(_cts(1, amps=[-1, 1])[0][0].real in [-1, 1]) + self.assertTrue(_cts(1, amps=[-1, 1])[0][0].imag in [-1, 1]) + self.assertEqual(len(_cts(5, amps=[-1, 1])), 5) + self.assertTrue(np.isin(_cts(20, amps=[-1, -3, 1, 3]).real, [-1, -3, 1, 3]).all()) + self.assertTrue(np.isin(_cts(20, amps=[-1, -3, 1, 3]).imag, [-1, -3, 1, 3]).all()) + with self.assertRaises(ValueError): + transmitted_symbols = _cts(1, amps=[-1.1, 1], quadrature=False) + with self.assertRaises(ValueError): + transmitted_symbols = _cts(1, amps=np.array([-1, 1.1]), quadrature=False) + with self.assertRaises(ValueError): + transmitted_symbols = _cts(1, amps=np.array([-1, 1+1j])) + + def test_complex_symbol_coding(self): + num_symbols = 5 + mod_pref = [1, 2, 3] + mods = ['QPSK', '16QAM', '64QAM'] + + for modI, mod in enumerate(mods): + num_bits = 2*num_symbols*mod_pref[modI] + max_symb = 2**mod_pref[modI]-1 + + #uniform encoding (max bits = max amplitude symbols): + bits = np.ones(num_bits) + symbols = max_symb*np.ones(num_symbols) + 1j*max_symb*np.ones(num_symbols) + symbols_enc = dimod.generators.wireless._bits_to_symbols(bits=bits, modulation=mod) + self.assertTrue(np.all(symbols_enc == symbols )) + bits_enc = dimod.generators.wireless._symbols_to_bits(symbols=symbols, modulation=mod) + self.assertTrue(np.all(bits_enc == bits)) + + #random encoding: + bits = self.rng.choice([-1, 1], size=num_bits) + symbols_enc = dimod.generators.wireless._bits_to_symbols(bits=bits, modulation=mod) + bits_enc = dimod.generators.wireless._symbols_to_bits(symbols=symbols_enc, modulation=mod) + self.assertTrue(np.all(bits_enc == bits)) + + def test_mimo(self): + for num_transmitters, num_receivers in [(1, 1), (5, 1), (1, 3), (11, 7)]: + F = self.rng.normal(0, 1, size=(num_receivers, num_transmitters)) + \ + 1j*self.rng.normal(0, 1, size=(num_receivers, num_transmitters)) + y = self.rng.normal(0, 1, size=(num_receivers, 1)) + \ + 1j*self.rng.normal(0, 1, size=(num_receivers, 1)) + bqm = dimod.generators.wireless.mimo(modulation='QPSK', y=y, F=F) + + mod_pref = [1, 1, 2, 3] + mods = ['BPSK', 'QPSK', '16QAM', '64QAM'] + for modI, modulation in enumerate(mods): + bqm = dimod.generators.wireless.mimo(modulation=modulation, + num_transmitters=num_transmitters, num_receivers=num_receivers) + if modulation == 'BPSK': + constellation = [-1, 1] + dtype = np.float64 + else: + max_val = 2**mod_pref[modI] - 1 + dtype = np.complex128 + # All 1 bit encoding (max symbol in constellation) + constellation = [real_part + 1j*imag_part + for real_part in range(-max_val, max_val+1, 2) + for imag_part in range(-max_val, max_val+1, 2)] + + F_simple = np.ones(shape=(num_receivers, num_transmitters), dtype=dtype) + transmitted_symbols_max = np.ones(shape=(num_transmitters, 1), + dtype=dtype)*constellation[-1] + transmitted_symbols_random = self.rng.choice(constellation, + size=(num_transmitters, 1)) + transmitted_bits_random = dimod.generators.wireless._symbols_to_bits( + symbols=transmitted_symbols_random.flatten(), modulation=modulation) + + #Trivial channel (F_simple), machine numbers + bqm = dimod.generators.wireless.mimo(modulation=modulation, + F=F_simple, transmitted_symbols=transmitted_symbols_max, + SNRb=float('Inf')) + + ef = self._effective_fields(bqm) + self.assertLessEqual(np.max(ef), 0) + self.assertLessEqual(abs(bqm.energy((np.ones(bqm.num_variables), + np.arange(bqm.num_variables)))), 1e-10) + + #Random channel, potential precision + bqm = dimod.generators.wireless.mimo(modulation=modulation, + num_transmitters=num_transmitters, num_receivers=num_receivers, + transmitted_symbols=transmitted_symbols_max, + SNRb=float('Inf')) + ef=self._effective_fields(bqm) + self.assertLessEqual(np.max(ef), 0) + self.assertLess(abs(bqm.energy((np.ones(bqm.num_variables), + np.arange(bqm.num_variables)))), 1e-8) + + # Add noise, check that offset is positive (random, scales as num_var/SNRb) + bqm = dimod.generators.wireless.mimo(modulation=modulation, + num_transmitters=num_transmitters, num_receivers=num_receivers, + transmitted_symbols=transmitted_symbols_max, SNRb=1) + self.assertLess(0, abs(bqm.energy((np.ones(bqm.num_variables), + np.arange(bqm.num_variables))))) + + # Random transmission, should match bit encoding. Bit-encoded energy should be minimal + bqm = dimod.generators.wireless.mimo(modulation=modulation, + num_transmitters=num_transmitters, num_receivers=num_receivers, + transmitted_symbols=transmitted_symbols_random, SNRb=float('Inf')) + self.assertLess(abs(bqm.energy((transmitted_bits_random, + np.arange(bqm.num_variables)))), 1e-8) + + def create_channel(self): + # Test some defaults + c, cp = dimod.generators.wireless.create_channel()[0] + self.assertEqual(cp, 2) + self.assertEqual(c.shape, (1, 1)) + + c, cp = dimod.generators.wireless.create_channel(5, 5, + F_distribution=("normal", "real")) + self.assertTrue(np.isin(c, [-1, 1]).all()) + self.assertEqual(cp, 5) + + c, cp = dimod.generators.wireless.create_channel(5, 5, + F_distribution=("binary", "complex")) + self.assertTrue(np.isin(c, [-1-1j, -1+1j, 1-1j, 1+1j]).all()) + self.assertEqual(cp, 10) + + n_trans = 40 + c, cp = dimod.generators.wireless.create_channel(30, n_trans, + F_distribution=("normal", "real")) + self.assertLess(c.mean(), 0.2) + self.assertLess(c.std(), 1.3) + self.assertGreater(c.std(), 0.7) + self.assertEqual(cp, n_trans) + + c, cp = dimod.generators.wireless.create_channel(30, n_trans, + F_distribution=("normal", "complex")) + self.assertLess(c.mean().complex, 0.2) + self.assertLess(c.real.std(), 1.3) + self.assertGreater(c.real.std(), 0.7) + self.assertEqual(cp, 2*n_trans) + + c, cp = dimod.generators.wireless.create_channel(5, 5, + F_distribution=("binary", "real"), + attenuation_matrix=np.array([[1, 2], [3, 4]])) + self.assertLess(c.ptp(), 8) + self.assertEqual(cp, 30) + + def test_create_signal(self): + # Only required parameters + got, sent, noise, _ = dimod.generators.wireless._create_signal(F=np.array([[1]])) + self.assertEqual(got, sent) + self.assertTrue(all(np.isreal(got))) + self.assertIsNone(noise) + + got, sent, _, __ = dimod.generators.wireless._create_signal(F=np.array([[-1]])) + self.assertEqual(got, -sent) + + got, sent, noise, _ = dimod.generators.wireless._create_signal(F=np.array([[1], [1]])) + self.assertEqual(got.shape, (2, 1)) + self.assertEqual(sent.shape, (1, 1)) + self.assertIsNone(noise) + + got, sent, _, __ = dimod.generators.wireless._create_signal(F=np.array([[1, 1]])) + self.assertEqual(got.shape, (1, 1)) + self.assertEqual(sent.shape, (2, 1)) + + # Optional parameters + got, sent, _, __ = dimod.generators.wireless._create_signal(F=np.array([[1]]), modulation="QPSK") + self.assertTrue(all(np.iscomplex(got))) + self.assertTrue(all(np.iscomplex(sent))) + self.assertEqual(got.shape, (1, 1)) + self.assertEqual(got, sent) + + got, sent, _, __ = dimod.generators.wireless._create_signal(F=np.array([[1]]), + transmitted_symbols=np.array([[1]])) + self.assertEqual(got, sent) + self.assertEqual(got[0][0], 1) + + with self.assertRaises(ValueError): # Complex symbols for BPSK + a, b, c, d = dimod.generators.wireless._create_signal(F=np.array([[1]]), + transmitted_symbols=np.array([[1+1j]])) + + with self.assertRaises(ValueError): # Non-complex symbols for non-BPSK + a, b, c, d = dimod.generators.wireless._create_signal(F=np.array([[1]]), + transmitted_symbols=np.array([[1]]), modulation="QPSK") + + noise = 0.2+0.3j + got, sent, _, __ = dimod.generators.wireless._create_signal(F=np.array([[1]]), + transmitted_symbols=np.array([[1]]), channel_noise=noise) + self.assertEqual(got, sent) + got, sent, _, __ = dimod.generators.wireless._create_signal(F=np.array([[1]]), + transmitted_symbols=np.array([[1]]), channel_noise=noise, SNRb=10 ) + self.assertEqual(got, sent + noise) + got, sent, _, __ = dimod.generators.wireless._create_signal(F=np.array([[1]]), + transmitted_symbols=np.array([[1]]), SNRb=10 ) + self.assertNotEqual(got, sent) + + def test_coordinated_multipoint(self): + bqm = dimod.generators.wireless.coordinated_multipoint(lattice=nx.complete_graph(1), + modulation='BPSK') + lattice = self._make_honeycomb(1) + bqm = dimod.generators.wireless.coordinated_multipoint(lattice=lattice) + num_var = lattice.number_of_nodes() + self.assertEqual(num_var,bqm.num_variables) + self.assertEqual(21,bqm.num_interactions) + # Transmitted symbols are 1 by default + lattice = self._make_honeycomb(2) + bqm = dimod.generators.wireless.coordinated_multipoint(lattice=lattice, + modulation='BPSK', SNRb=float('Inf')) + self.assertLess(abs(bqm.energy((np.ones(bqm.num_variables), bqm.variables))), 1e-10) + + def test_attenuation_matrix(self): + #Check that attenuation matches the matrix + lattice=nx.Graph() + num_var = 10 + lattice.add_nodes_from(n for n in range(num_var)) + + A,_,_ = dimod.generators.wireless._lattice_to_attenuation_matrix(lattice) + self.assertFalse(np.any(A-np.identity(num_var))) + + for t_per_node in range(1,3): + for r_per_node in range(1,3): + A,_,_ = dimod.generators.wireless._lattice_to_attenuation_matrix( + lattice, + transmitters_per_node=t_per_node, + receivers_per_node=r_per_node, + neighbor_root_attenuation=self.rng.random()) + self.assertFalse(np.any(A - np.tile(np.identity(num_var), (r_per_node, t_per_node)))) + + for ea in range(2): + lattice.add_edge(ea, ea+1) + neighbor_root_attenuation = self.rng.random() + A,_,_ = dimod.generators.wireless._lattice_to_attenuation_matrix( + lattice, neighbor_root_attenuation=2) + self.assertFalse(np.any(A - A.transpose())) + self.assertTrue(all(A[eap, eap + 1]==2 for eap in range(ea + 1))) + + ## Check num_transmitters and num_receivers override: + nx.set_node_attributes(lattice, values=3, name="num_transmitters") + nx.set_node_attributes(lattice, values=1, name="num_receivers") + A,_,_ = dimod.generators.wireless._lattice_to_attenuation_matrix( + lattice, + transmitters_per_node=2, + receivers_per_node=2) + self.assertEqual(A.shape, (num_var, 3*num_var)) + nx.set_node_attributes(lattice, values=0, name="num_receivers") + nx.set_node_attributes(lattice, values={0:2, 3:1}, name="num_receivers") + nx.set_node_attributes(lattice, values=0, name="num_transmitters") + nx.set_node_attributes(lattice, values={i:1 for i in [0,1,2,4]}, name="num_transmitters") + + # t/r2 -- t -- t r t #We can assume the ntr and ntt arguments. + Acorrect = np.array([[1, 2, 0, 0], [1, 2, 0, 0], [0, 0, 0, 0]]) + A,_,_ = dimod.generators.wireless._lattice_to_attenuation_matrix( + lattice, neighbor_root_attenuation=2) + self.assertFalse(np.any(A - Acorrect)) + + def test_noise_scale(self): + # After applying use_offset, the expected energy is the sum of noise terms. + # (num_transmitters/SNRb)*sum_{mu=1}^{num_receivers} nu_mu^2 , where =1 under default channels + # We can do a randomized test (for practicl purpose, I fix the seed to avoid rare outliers): + for num_transmitters in [256]: + for SNRb in [0.1]: #[0.1,10] + for mods in [('BPSK', 1, 1, 1),('64QAM', 2, 42, 6)]: #,('QPSK',2,2,2),('16QAM',2,10,4)]: + mod, channel_power_per_transmitter, constellation_mean_power, bits_per_transmitter = mods + for num_receivers in [num_transmitters*4]: #[num_transmitters//4,num_transmitters]: + EoverN = (channel_power_per_transmitter * \ + constellation_mean_power/bits_per_transmitter/SNRb) * \ + num_transmitters * num_receivers + if mod=='BPSK': + EoverN *= 2 #Real part only + for seed in range(1): + bqm0 = dimod.generators.wireless.mimo(modulation=mod, + num_transmitters=num_transmitters, num_receivers=num_receivers, + seed=seed) + bqm = dimod.generators.wireless.mimo(modulation=mod, + num_transmitters=num_transmitters, num_receivers=num_receivers, + SNRb=SNRb, seed=seed) + + #E[n^2] constructed from offsets correctly: + scale_n = (bqm.offset - bqm0.offset)/EoverN + self.assertGreater(1.5, scale_n) + self.assertLess(0.5, scale_n) + #scale_n_alt = np.sum(abs(n)**2,axis=0)/EoverN) + + for num_transmitter_block in [2]: #[1,2]: + lattice_size = num_transmitters//num_transmitter_block + for num_receiver_block in [1]: #[1,2]: + # Similar applies for COMP, up to boundary conditions. Choose a symmetric lattice: + num_receiversT = lattice_size * num_receiver_block + num_transmittersT = lattice_size * num_transmitter_block + EoverN = (channel_power_per_transmitter * \ + constellation_mean_power/bits_per_transmitter/SNRb) * \ + num_transmittersT * num_receiversT + + if mod=='BPSK': + EoverN *= 2 #Real part only + lattice = nx.Graph() + lattice.add_edges_from((i, (i + 1)%lattice_size) for i in + range(num_transmitters//num_transmitter_block)) + for seed in range(1): + bqm = dimod.generators.wireless.coordinated_multipoint(lattice=lattice, + modulation=mod, SNRb=SNRb) + bqm0 = dimod.generators.wireless.coordinated_multipoint(lattice=lattice, + modulation=mod) + scale_n = (bqm.offset - bqm0.offset)/EoverN + self.assertGreater(1.5, scale_n) + #self.assertLess(0.5, scale_n) +