diff --git a/pyqpanda-algorithm/example/QJEPG/demo_QJPEG.py b/pyqpanda-algorithm/example/QJEPG/demo_QJPEG.py new file mode 100644 index 0000000..db32f46 --- /dev/null +++ b/pyqpanda-algorithm/example/QJEPG/demo_QJPEG.py @@ -0,0 +1,82 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Demo: Quantum JPEG (QJPEG) image downsampling. + +The given image is down-sampled at several compression levels and compared +with classical box-averaging. + +Run: + python demo_QJPEG.py (default -P equals -S, no tiling) + python demo_QJPEG.py -P 256 -S 512 (large canvas with small patch, force tiling) +""" + +from pathlib import Path +from argparse import ArgumentParser +from time import perf_counter + +import numpy as np +import matplotlib.pyplot as plt + +import sys +sys.path.insert(0, str(Path(__file__).resolve().parents[2])) +from utils import load_img, psnr, box_downsample, QJPEG, IMG_DIR, OUT_DIR + + +def save_filename(args): + I: Path = args.input + S = args.img_size + P = args.patch_size + sfx = '' + sfx += f'_{P}' if S == P else f'_P{P}_S{S}' + if args.no_H_layer is False: + sfx += '_noH' + return f'{I.stem}-QJPEG{sfx}.png' + + +def run(args): + img = load_img(args.input, args.img_size) + print(f'Input img_shape: {img.shape}') + + qjpeg = QJPEG(patch_size=args.patch_size, wrap_H_layer=args.no_H_layer) + print(f'QJPEG patch_size: {qjpeg.patch_size}') + + results: dict[int, np.ndarray] = {} + for n_discard in (1, 2, 3): + ts_start = perf_counter() + out = qjpeg(img, n_discard=n_discard) + ref = box_downsample(img, 2 ** n_discard) + ts_end = perf_counter() + results[n_discard] = out + print(f'n_discard={n_discard}: {img.shape} -> {out.shape} | PSNR={psnr(out, ref):.2f} dB ({ts_end - ts_start:.3f}s)') + + N = 1 + len(results) + fig, axes = plt.subplots(1, N, figsize=(3 * N, 3)) + axes[0].imshow(img, vmin=0, vmax=1) + axes[0].set_title(f'original {img.shape[0]}x{img.shape[1]}') + axes[0].axis('off') + for ax, (nd, out) in zip(axes[1:], results.items()): + ax.imshow(np.clip(out, 0, 1), vmin=0, vmax=1) + ax.set_title(f'QJPEG {out.shape[0]}x{out.shape[1]}\nn_discard={nd}') + ax.axis('off') + fig.tight_layout() + OUT_DIR.mkdir(exist_ok=True) + fp = OUT_DIR / save_filename(args) + fig.savefig(fp, dpi=300, bbox_inches='tight') + print(f'Saved comparison figure to {fp}') + + +if __name__ == '__main__': + parser = ArgumentParser() + parser.add_argument('-I', '--input', default=IMG_DIR / 'Border Collie.webp', type=Path) + parser.add_argument('-S', '--img_size', default=256, type=int) + parser.add_argument('-P', '--patch_size', default=None, type=int) + parser.add_argument('-nH', '--no_H_layer', action='store_false', default='disable H layer, only experimental use') + args = parser.parse_args() + + args.patch_size = args.patch_size or args.img_size + + if args.img_size > 256: + print(f'>> [WARN] the first image size will exceed {args.img_size * 2**2}, might be VERY SLOW!! :(') + + run(args) diff --git a/pyqpanda-algorithm/example/QJEPG/demo_cyclic.py b/pyqpanda-algorithm/example/QJEPG/demo_cyclic.py new file mode 100644 index 0000000..39e59bc --- /dev/null +++ b/pyqpanda-algorithm/example/QJEPG/demo_cyclic.py @@ -0,0 +1,88 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Demo: QJPEG & unQJPEG forward-inverse cyclic. + +Down-samples the given image with QJPEG, then up-samples it back with the unQJPEG +quantum circuit (spectral zero-padding / sinc interpolation) and reports how well +the original resolution is recovered. Also demonstrates plain up-sampling of a +low-resolution image. + +Run: + python demo_cyclic.py + python demo_cyclic.py -S 512 +""" + +from pathlib import Path +from argparse import ArgumentParser +from time import perf_counter + +import numpy as np +import matplotlib.pyplot as plt + +import sys +sys.path.insert(0, str(Path(__file__).resolve().parents[2])) +from utils import load_img, psnr, bilinear_upsample, QJPEG, unQJPEG, IMG_DIR, OUT_DIR + + +def save_filename(args): + I: Path = args.input + S = args.size + sfx = '' + sfx += f'_{S}' + return f'{I.stem}-cyclic{sfx}.png' + + +def run(args): + size = args.size + half = size // 2 + orig = load_img(args.input, size) + print(f'Input img_shape: {orig.shape}') + + # 1) downsample to half resolution with QJPEG + ts_start = perf_counter() + lowres = QJPEG(patch_size=size)(orig, n_discard=1) + ts_end = perf_counter() + print(f'Downsampled with QJPEG: {orig.shape} -> {lowres.shape} ({ts_end - ts_start:.3f}s)') + + # 2) upsample back to full resolution with unQJPEG + ts_start = perf_counter() + recon = unQJPEG()(lowres, n_append=1) + ts_end = perf_counter() + psnr_recon = psnr(recon, orig) + print(f'Upsampled with unQJPEG: {lowres.shape} -> {recon.shape} | PSNR={psnr_recon:.2f} dB ({ts_end - ts_start:.3f}s)') + + # 3) bilinear upscaling of the same lowres image + ts_start = perf_counter() + bilinear = bilinear_upsample(lowres, 2) + ts_end = perf_counter() + psnr_bilinear = psnr(bilinear, orig) + print(f'Upsampled with BILINEAR: {lowres.shape} -> {bilinear.shape} | PSNR={psnr_bilinear:.2f} dB ({ts_end - ts_start:.3f}s)') + + imgs = [orig, np.clip(lowres, 0, 1), np.clip(recon, 0.0, 1.0), bilinear] + titles = [ + f'original {size}x{size}', + f'QJPEG {half}x{half}', + f'unQJPEG {size}x{size}\nPSNR={psnr_recon:.2f}', + f'bilinear {size}x{size}\nPSNR={psnr_bilinear:.2f}', + ] + + fig, axes = plt.subplots(1, 4, figsize=(12, 3)) + for ax, im, title in zip(axes, imgs, titles): + ax.imshow(im, vmin=0, vmax=1) + ax.set_title(title) + ax.axis('off') + fig.tight_layout() + OUT_DIR.mkdir(exist_ok=True) + fp = OUT_DIR / save_filename(args) + fig.savefig(fp, dpi=300, bbox_inches='tight') + print(f'Saved comparison figure to {fp}') + + +if __name__ == '__main__': + parser = ArgumentParser() + parser.add_argument('-I', '--input', default=IMG_DIR / 'Border Collie.webp', type=Path) + parser.add_argument('-S', '--size', default=256, type=int) + args = parser.parse_args() + + run(args) diff --git a/pyqpanda-algorithm/example/QJEPG/demo_unQJPEG.py b/pyqpanda-algorithm/example/QJEPG/demo_unQJPEG.py new file mode 100644 index 0000000..c82ceab --- /dev/null +++ b/pyqpanda-algorithm/example/QJEPG/demo_unQJPEG.py @@ -0,0 +1,79 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Demo: inverse Quantum JPEG (unQJPEG) image upsampling / super-resolution. + +The given image is up-sampled at several upscale levels and compared +with classical bilinear interpolation. + +Run: + python demo_unQJPEG.py (default -P equals -S, no tiling) + python demo_unQJPEG.py -P 32 -S 64 (large canvas with small patch, force tiling) +""" + +from pathlib import Path +from argparse import ArgumentParser +from time import perf_counter + +import numpy as np +import matplotlib.pyplot as plt + +import sys +sys.path.insert(0, str(Path(__file__).resolve().parents[2])) +from utils import load_img, psnr, bilinear_upsample, unQJPEG, IMG_DIR, OUT_DIR + + + +def save_filename(args): + I: Path = args.input + S = args.img_size + P = args.patch_size + sfx = '' + sfx += f'_{P}' if S == P else f'_P{P}_S{S}' + return f'{I.stem}-unQJPEG{sfx}.png' + + +def run(args): + img = load_img(args.input, args.img_size) + print(f'Input img_shape: {img.shape}') + + unqjpeg = unQJPEG(separated_QFT=args.no_separated_QFT) + print(f'patch_size: {args.patch_size}') + + results: dict[int, np.ndarray] = {} + for n_append in (1, 2, 3): + ts_start = perf_counter() + out = unqjpeg(img, n_append=n_append, patch_size=args.patch_size) + ref = bilinear_upsample(img, 2 ** n_append) + ts_end = perf_counter() + results[n_append] = out + print(f'n_append={n_append}: {img.shape} -> {out.shape} | PSNR={psnr(out, ref):.2f} dB ({ts_end - ts_start:.3f}s)') + + N = 1 + len(results) + fig, axes = plt.subplots(1, N, figsize=(3 * N, 3)) + axes[0].imshow(img, vmin=0, vmax=1) + axes[0].set_title(f'original {img.shape[0]}x{img.shape[1]}') + axes[0].axis('off') + for ax, (nd, out) in zip(axes[1:], results.items()): + ax.imshow(np.clip(out, 0, 1), vmin=0, vmax=1) + ax.set_title(f'unQJPEG {out.shape[0]}x{out.shape[1]}\nn_append={nd}') + ax.axis('off') + fig.tight_layout() + OUT_DIR.mkdir(exist_ok=True) + fp = OUT_DIR / save_filename(args) + fig.savefig(fp, dpi=300, bbox_inches='tight') + print(f'Saved comparison figure to {fp}') + + +if __name__ == '__main__': + parser = ArgumentParser() + parser.add_argument('-I', '--input', default=IMG_DIR / 'Border Collie.webp', type=Path) + parser.add_argument('-S', '--img_size', default=64, type=int) + parser.add_argument('-P', '--patch_size', default=None, type=int) + parser.add_argument('-nQFT', '--no_separated_QFT', action='store_false', help='disable separated QFT, only experimental use') + args = parser.parse_args() + + if args.img_size > 64: + print(f'>> [WARN] the last image size will exceed {args.img_size * 2**3}, could be VERY SLOW!! :(') + + run(args) diff --git a/pyqpanda-algorithm/example/QJEPG/img/Border Collie.webp b/pyqpanda-algorithm/example/QJEPG/img/Border Collie.webp new file mode 100644 index 0000000..30b7915 Binary files /dev/null and b/pyqpanda-algorithm/example/QJEPG/img/Border Collie.webp differ diff --git a/pyqpanda-algorithm/example/QJEPG/output/Border Collie-QJPEG_256.png b/pyqpanda-algorithm/example/QJEPG/output/Border Collie-QJPEG_256.png new file mode 100644 index 0000000..8695217 Binary files /dev/null and b/pyqpanda-algorithm/example/QJEPG/output/Border Collie-QJPEG_256.png differ diff --git a/pyqpanda-algorithm/example/QJEPG/output/Border Collie-cyclic_256.png b/pyqpanda-algorithm/example/QJEPG/output/Border Collie-cyclic_256.png new file mode 100644 index 0000000..0b370ad Binary files /dev/null and b/pyqpanda-algorithm/example/QJEPG/output/Border Collie-cyclic_256.png differ diff --git a/pyqpanda-algorithm/example/QJEPG/output/Border Collie-unQJPEG_PNone_S64.png b/pyqpanda-algorithm/example/QJEPG/output/Border Collie-unQJPEG_PNone_S64.png new file mode 100644 index 0000000..036ccf0 Binary files /dev/null and b/pyqpanda-algorithm/example/QJEPG/output/Border Collie-unQJPEG_PNone_S64.png differ diff --git a/pyqpanda-algorithm/example/QJEPG/utils.py b/pyqpanda-algorithm/example/QJEPG/utils.py new file mode 100644 index 0000000..2caac10 --- /dev/null +++ b/pyqpanda-algorithm/example/QJEPG/utils.py @@ -0,0 +1,52 @@ +import sys +from pathlib import Path + +import numpy as np +from numpy import ndarray +from PIL import Image +from PIL.Image import Resampling +import matplotlib ; matplotlib.use('Agg') + +sys.path.insert(0, str(Path(__file__).resolve().parents[2])) +from pyqpanda_alg.QJEPG import QJPEG, unQJPEG # keep + +BASE_PATH = Path(__file__).resolve().parent +IMG_DIR = BASE_PATH / 'img' +OUT_DIR = BASE_PATH / 'output' + + +def load_img(path: Path, size: int) -> ndarray: + img = Image.open(path).resize((size, size), resample=Resampling.LANCZOS) + return np.asarray(img, dtype=np.float64) / 255.0 # [H, W, C] / [H, W] + + +def psnr(a: ndarray, b: ndarray) -> float: + mse = np.mean((a - b) ** 2) + return float('inf') if mse == 0 else 10 * np.log10(1.0 / mse) + + +def nearest_neighbor_upsample(img: ndarray, factor: int = 2) -> ndarray: + return np.repeat(np.repeat(img, factor, axis=0), factor, axis=1) + + +def bilinear_upsample(img: ndarray, factor: int = 2) -> ndarray: + if len(img.shape) == 3: + H, W, C = img.shape + else: + H, W = img.shape + im = np.asarray(np.clip(img, 0, 1) * 255, dtype=np.uint8) + img = Image.fromarray(im).resize((W * factor, H * factor), resample=Resampling.BILINEAR) + hires = np.asarray(img, dtype=np.float32) / 255.0 + return hires + + +def box_downsample(img: ndarray, factor: int = 2) -> ndarray: + if len(img.shape) == 3: + H, W, C = img.shape + p = H // factor + lowres = img.reshape(p, factor, p, factor, C).mean(axis=(1, 3)) + else: + H, W = img.shape + p = H // factor + lowres = img.reshape(p, factor, p, factor).mean(axis=(1, 3)) + return lowres diff --git a/pyqpanda-algorithm/pyqpanda_alg/QJEPG/QJPEG.py b/pyqpanda-algorithm/pyqpanda_alg/QJEPG/QJPEG.py new file mode 100644 index 0000000..6b9db0c --- /dev/null +++ b/pyqpanda-algorithm/pyqpanda_alg/QJEPG/QJPEG.py @@ -0,0 +1,287 @@ +# 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 or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +import math +import numpy as np +from numpy import ndarray + +from pyqpanda3.core import CPUQVM, QProg, Encode, H + +from ..plugin import QFT + + +class QJPEG: + """ + Quantum JPEG (QJPEG) image down-sampling / compression. + + The QJPEG algorithm encodes an image patch into the amplitudes of a quantum + state (QPIE, Quantum Probability Image Encoding), applies the quantum Fourier + transform to move to the spatial-frequency domain, discards the high + spatial-frequency qubits, and transforms back. Reading out the measurement + probabilities of the remaining qubits yields a lower-resolution image, in the + same spirit as the classical JPEG that filters high-frequency components [1]. + + The image is first split into square ``patch_size`` x ``patch_size`` patches. + Each patch is processed independently by a ``2 * log2(patch_size)`` qubit + circuit and down-sampled by a factor of ``2 ** n_discard`` along each axis, so + a ``P x P`` patch becomes a ``(P / 2**n_discard) x (P / 2**n_discard)`` patch. + + Parameters + patch_size : ``int``\n + Side length of the square patches, must be a power of two. The circuit + uses ``2 * log2(patch_size)`` qubits, which is required to be in + ``[4, 20]`` (i.e. ``patch_size`` in ``[4, 1024]``). + wrap_H_layer : ``bool``, optional\n + Whether to wrap the transform with Hadamard layers, as proposed in [1]. + The Hadamards reduce the statistical fluctuations of each pixel while + preserving the contrast of the image, making the reconstruction faithful + to the classical box down-sampling. Strongly recommended (Default: True). + shots : ``int``, optional\n + If given, the output probabilities are sampled with this many + measurement shots (emulating finite-sampling noise on real hardware). + If ``None`` (Default), the exact probabilities are used. + seed : ``int``, optional\n + Random seed used when ``shots`` is set (Default: None). + + References + [1] Simone Roncallo, Lorenzo Maccone, Chiara Macchiavello. "Quantum JPEG". + AVS Quantum Science 5, 043803 (2023). arXiv:2306.09323. + + >>> import numpy as np + >>> from pyqpanda_alg.QJEPG import QJPEG + >>> img = np.random.rand(16, 16).astype(np.float64) + >>> qjpeg = QJPEG(patch_size=16) + >>> out = qjpeg(img, n_discard=1) # 16x16 -> 8x8 + >>> print(out.shape) + (8, 8) + """ + + def __init__(self, patch_size: int = 256, wrap_H_layer: bool = True, shots: int = None, seed: int = None): + assert patch_size & (patch_size - 1) == 0, '"patch_size" must be power of 2' + self.patch_size = patch_size + n_qubit = int(math.log2(patch_size)) * 2 + assert 4 <= n_qubit <= 20, f'invalid input "n_qubit": {n_qubit}, should be in range [4, 20]' + self.n_qubit = n_qubit # i.e. k_in/max_qubit + + self.wrap_H_layer = wrap_H_layer + self.shots = shots + self.seed = seed + # per-call state, set inside __call__ + self._grid: tuple[int, int] = None # (n_patch_H, n_patch_W) tiling of the current channel + self._scale = 1 # down-sampling factor per axis (2 ** n_discard) + + def __call__(self, img: ndarray, n_discard: int = 2) -> ndarray: + assert isinstance(img, ndarray) and img.dtype in [np.float32, np.float64], 'image must be a float ndarray' + assert 0 <= img.min() and img.max() <= 1.0, 'pixel value must be in [0, 1]' + orig_dtype = img.dtype + if orig_dtype != np.float64: # tackle precision issue + img = img.astype(np.float64) + if len(img.shape) == 3: + assert img.shape[-1] in [3, 4], "colored image must be in shape (H, W, C) with C in [3, 4]" + is_grey = False + else: + assert len(img.shape) == 2, "grey image must be in shape (H, W)" + is_grey = True + img = np.expand_dims(img, -1) + assert self.patch_size <= min(img.shape[:2]), '"patch_size" must be smaller than image size' + k = self.n_qubit // 2 # qubits per axis == log2(patch_size) + assert isinstance(n_discard, int) and 1 <= n_discard <= k - 1, f'"n_discard" must be an int in [1, {k - 1}] for patch_size={self.patch_size}' + self._scale = 2 ** n_discard + + H, W, C = img.shape + # pad to whole patches, [H, W, C] + img_pad, pads = self.pad(img) + # split channels, C * [H, W] + channels = [img_pad[..., c] for c in range(C)] + channels_processed: list[ndarray] = [] + # foreach C + for channel in channels: + # patchify, [B=(H*W/P**2), h=P, w=P] + patches = self.pachify(channel) + # vectorize, [B, D=P**2] + vectors, norm = self.vectorize(patches) + # QFT-discard-IQFT, [B, d=p**2] + qvectors = self.apply_QJPEG(vectors, n_discard) + # devectorize, [B, p, p] + devectors = self.devectorize(qvectors, norm) + # unpatchify, [Ho, Wo] + channel_unpatchify = self.unpachify(devectors) + channels_processed.append(channel_unpatchify) + # merge channel, [Ho, Wo, C] + img_merged = np.stack(channels_processed, axis=-1) + # trim pads, [Ho, Wo, C] + img_unpad = self.unpad(img_merged, pads) + # dtype back + img_unpad = img_unpad.astype(orig_dtype) + return img_unpad[..., 0] if is_grey else img_unpad + + def pad(self, img: ndarray): + """Zero-pad an ``[H, W, C]`` image so that H and W are whole multiples of + ``patch_size``, pasting the original content to the center. Returns the + padded image together with the ``(top, bottom, left, right)`` pad widths + (in original pixels), or ``None`` when no padding is needed.""" + H, W, C = img.shape + n_patch_H = math.ceil(H / self.patch_size) + n_patch_W = math.ceil(W / self.patch_size) + H_ex = n_patch_H * self.patch_size + W_ex = n_patch_W * self.patch_size + if H_ex == H and W_ex == W: + return img, None + pad_H = H_ex - H + pad_W = W_ex - W + pads = ( + math.floor(pad_H / 2), math.ceil(pad_H / 2), + math.floor(pad_W / 2), math.ceil(pad_W / 2), + ) + img_ex = np.zeros(shape=(H_ex, W_ex, C), dtype=img.dtype) + img_ex[pads[0]:pads[0] + H, pads[2]:pads[2] + W, :] = img + return img_ex, pads + + def unpad(self, img: ndarray, pads=None) -> ndarray: + """Trim the padding added by :meth:`pad` from an ``[H, W, C]`` image. The + pad widths are scaled down by the down-sampling factor to match the + processed (lower) resolution.""" + if not pads or all(e == 0 for e in pads): + return img + top, bottom, left, right = pads + s = self._scale + H, W, C = img.shape + t = round(top / s) + b = round(bottom / s) + l = round(left / s) + r = round(right / s) + return img[t:(H - b) if b else H, l:(W - r) if r else W, :] + + def pachify(self, img: ndarray) -> ndarray: + """Split a single-channel ``[H, W]`` image into non-overlapping + ``patch_size`` x ``patch_size`` patches, returned as ``[B, P, P]`` in + row-major patch order.""" + P = self.patch_size + H, W = img.shape + assert H % P == 0 and W % P == 0 + nH, nW = H // P, W // P + self._grid = (nH, nW) + return img.reshape(nH, P, nW, P).swapaxes(1, 2).reshape(nH * nW, P, P) + + def unpachify(self, patches: ndarray) -> ndarray: + """Re-assemble ``[B, p, p]`` (processed) patches into a single-channel + ``[Ho, Wo]`` image using the tiling recorded by :meth:`pachify`.""" + nH, nW = self._grid + B, p, _ = patches.shape + assert B == nH * nW + return patches.reshape(nH, nW, p, p).swapaxes(1, 2).reshape(nH * p, nW * p) + + def vectorize(self, patches: ndarray) -> tuple[ndarray, ndarray]: + """Vectorize ``[B, P, P]`` patches to amplitude vectors ``[B, P**2]``. + + Each patch is flattened (row-major) and encoded as a normalized quantum + state whose amplitudes are the square-roots of the intensities, i.e. + ``amplitude = sqrt(pixel / patch_sum)`` (so ``|amplitude|**2`` recovers the + normalized intensity). The per-patch intensity sum is returned in ``norm`` + so the absolute scale can be restored on decoding. Empty (all-zero) patches + map to a zero vector. + """ + B = patches.shape[0] + vect = patches.reshape(B, -1) + norm = vect.sum(axis=1) + amps = np.zeros_like(vect, dtype=patches.dtype) + nz = norm > 0 + amps[nz] = np.sqrt(vect[nz] / norm[nz, None]) + return amps, norm + + def devectorize(self, vectors: ndarray, norm: ndarray) -> ndarray: + """Decode output probability vectors ``[B, p**2]`` back into intensity + patches ``[B, p, p]``. The probabilities are rescaled by the stored patch + intensity ``norm`` and by the pixel-count ratio so that the mean intensity + of each patch is preserved by the down-sampling.""" + B, d = vectors.shape + p = int(round(math.sqrt(d))) + n_in = self.patch_size ** 2 + n_out = d + patches = vectors * norm[:, None] * (n_out / n_in) + return patches.reshape(B, p, p) + + def apply_QJPEG(self, vectors: ndarray, n_discard: int) -> ndarray: + """Run the QJPEG quantum circuit on each amplitude vector ``[B, P**2]`` and + return the output measurement-probability vectors ``[B, p**2]``. + + For every patch the circuit (i) amplitude-encodes the state, (ii) optionally + applies a Hadamard layer, (iii) applies the QFT over all ``n0`` qubits, + (iv) applies the inverse QFT over the lowest ``n1 = n0 - n_discard`` qubits, + (v) discards the ``n_discard`` high-frequency qubits at the top of each + (row / column) half, and (vi) reads the probabilities of the remaining + ``n2 = n0 - 2 * n_discard`` qubits. + """ + + ''' + LSB MSB + |--------------nq---------------| + |------ki-------|------ki-------| + |---ko---| drop |---ko---| drop | + ''' + nq = self.n_qubit # all qubits + ki = nq // 2 # half qubits (aka. k_in) + ko = ki - n_discard # half qubits kept (aka. k_out) + qv = list(range(nq)) + # drop the top(MSB)/last n_discard qubits of each half/axis + dropped = [i for p in range(2) for i in range(ki * p + ko, ki * p + ki)] + # qubits kept for measurement + measured = [i for i in qv if i not in dropped] + + B = vectors.shape[0] + d = 2 ** (nq - 2 * n_discard) # dim out + out = np.zeros((B, d), dtype=vectors.dtype) + rng = np.random.default_rng(self.seed) if self.shots else None + for b in range(B): + vec = vectors[b] + if not np.any(vec): continue + + prog = QProg(nq) + # data enc + enc = Encode() + enc.amplitude_encode(qv, list(vec)) + prog << enc.get_circuit() + # QFT + if self.wrap_H_layer: # 抵消QFT中的Hlayer + for qi in qv: + prog << H(qi) + prog << QFT(qv) + # 舍去列比特 (末尾之前ntilde个) + # iQFT + prog << QFT(qv[:nq - n_discard]).dagger() + if self.wrap_H_layer: # 抵消iQFT中的Hlayer + for qi in measured: + prog << H(qi) + # 舍去行比特 (中部之前ntilde个，仅保留measured) + + machine = CPUQVM() + machine.run(prog, 1) + prob_dict = machine.result().get_prob_dict(measured) + probs = self._prob_dict_to_vector(prob_dict, len(measured), d, vectors.dtype) + if self.shots: + probs = rng.multinomial(self.shots, probs) / self.shots + out[b] = probs + return out + + @staticmethod + def _prob_dict_to_vector(prob_dict: dict, n: int, dim: int, dtype=np.float64) -> ndarray: + """Convert a ``get_prob_dict`` result into a dense vector indexed so that + ``qubits[0]`` is the least-significant bit (matching the QPIE encoding).""" + vec = np.zeros(dim, dtype=dtype) + for key, val in prob_dict.items(): + idx = 0 + for k in range(n): + # key is big-endian over `qubits`: leftmost char == last qubit + if key[n - 1 - k] == '1': + idx |= (1 << k) + vec[idx] += val + return vec diff --git a/pyqpanda-algorithm/pyqpanda_alg/QJEPG/__init__.py b/pyqpanda-algorithm/pyqpanda_alg/QJEPG/__init__.py new file mode 100644 index 0000000..7bf0c57 --- /dev/null +++ b/pyqpanda-algorithm/pyqpanda_alg/QJEPG/__init__.py @@ -0,0 +1,18 @@ +''' +Adapted implementation of the essay "Quantum JPEG" (arXiv:2306.09323) from the official repo, +and the inversed form (i.e. unQJPEG) as we originally proposed :) + +The QJPEG algorithm uses the quantum Fourier transform to discard the high spatial-frequency +qubits of an image, downsampling it to a lower resolution. This allows one to capture, +compress, and send images even with limited quantum resources for storage and communication. +~ref: https://arxiv.org/abs/2306.09323 +~ref: https://github.com/simoneroncallo/quantum-jpeg +''' + +from .QJPEG import QJPEG +from .unQJPEG import unQJPEG + +__all__ = [ + 'QJPEG', + 'unQJPEG', +] diff --git a/pyqpanda-algorithm/pyqpanda_alg/QJEPG/unQJPEG.py b/pyqpanda-algorithm/pyqpanda_alg/QJEPG/unQJPEG.py new file mode 100644 index 0000000..c0af239 --- /dev/null +++ b/pyqpanda-algorithm/pyqpanda_alg/QJEPG/unQJPEG.py @@ -0,0 +1,287 @@ +# 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 or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +import math +import numpy as np +from numpy import ndarray + +from pyqpanda3.core import CPUQVM, QProg, Encode, H, CNOT + +from ..plugin import QFT + + +class unQJPEG: + """ + Inverse Quantum JPEG (unQJPEG) image up-sampling / super-resolution. + + unQJPEG is the inverse of :class:`~pyqpanda_alg.QJEPG.QJPEG`, originally + proposed within this library. Where QJPEG *discards* the high spatial-frequency + qubits to down-sample an image, unQJPEG *re-introduces* them as vacuum (zeros), + i.e. it zero-pads the spatial-frequency spectrum and transforms back, which is + the quantum analogue of sinc / Fourier interpolation. + + Each ``patch_size`` x ``patch_size`` patch is encoded into ``2 * log2(patch_size)`` + qubits (QPIE amplitude encoding). For each of the two image axes the quantum + Fourier transform maps the register to the frequency domain; ``n_append`` fresh + ``|0>`` qubits are inserted between the low- and high-frequency components (the + high-frequency band is padded with zeros); an inverse QFT of the enlarged + register returns to the spatial domain. Reading the probabilities yields an + up-sampled patch of side ``patch_size * 2 ** n_append``. + + The sign qubit (most-significant frequency bit) of each axis is duplicated onto + the appended qubits with ``CNOT`` gates so that the negative-frequency band is + mapped to the top of the enlarged register, preserving the Hermitian frequency + layout required for a real, interpolating reconstruction. + + Parameters + separated_QFT : ``bool``, optional\n + Whether the QFT operation should be axis-wise (Default: True). + shots : ``int``, optional\n + If given, the output probabilities are sampled with this many + measurement shots (emulating finite-sampling noise). If ``None`` + (Default), the exact probabilities are used. + seed : ``int``, optional\n + Random seed used when ``shots`` is set (Default: None). + + References + [1] Simone Roncallo, Lorenzo Maccone, Chiara Macchiavello. "Quantum JPEG". + AVS Quantum Science 5, 043803 (2023). arXiv:2306.09323. + + >>> import numpy as np + >>> from pyqpanda_alg.QJEPG import unQJPEG + >>> img = np.random.rand(8, 8).astype(np.float64) + >>> unqjpeg = unQJPEG() + >>> out = unqjpeg(img, n_append=1, patch_size=8) # 8x8 -> 16x16 + >>> print(out.shape) + (16, 16) + """ + + def __init__(self, separated_QFT: bool = True, shots: int = None, seed: int = None): + self.separated_QFT = separated_QFT + self.shots = shots + self.seed = seed + # per-call state, set inside __call__ + self.patch_size: int = None + self.n_qubit: int = None # aka. k_out/max_qubit + self._grid: tuple[int, int] = None # (n_patch_H, n_patch_W) tiling of the current channel + self._scale = 1 # up-sampling factor per axis (2 ** n_append) + + def __call__(self, img: ndarray, n_append: int = 1, patch_size: int = None) -> ndarray: + assert isinstance(img, ndarray) and img.dtype in [np.float32, np.float64], 'image must be a float array' + assert 0 <= img.min() and img.max() <= 1.0, 'pixel value must be in [0, 1]' + orig_dtype = img.dtype + if orig_dtype != np.float64: # tackle precision issue + img = img.astype(np.float64) + if len(img.shape) == 3: + assert img.shape[-1] in [3, 4], "colored image must be in shape (H, W, C) with C in [3, 4]" + is_grey = False + else: + assert len(img.shape) == 2, "grey image must be in shape (H, W)" + is_grey = True + img = np.expand_dims(img, -1) + assert isinstance(n_append, int) and n_append >= 1, '"n_append" must be a positive int' + if patch_size is None: + patch_size = min(img.shape[:2]) + assert patch_size & (patch_size - 1) == 0, '"patch_size" must be power of 2' + assert patch_size <= min(img.shape[:2]), '"patch_size" must be smaller than image size' + self.patch_size = patch_size + nq_init = int(math.log2(patch_size)) * 2 + assert 2 <= nq_init <= 18, f'invalid input "nq_init": {nq_init}, should be in range [2, 18]' + n_qubit = nq_init + n_append * 2 # qubits of the enlarged (output) register + assert 4 <= n_qubit <= 20, f'invalid output "n_qubit": {n_qubit}, should be in range [4, 20]; reduce patch_size or n_append' + self.n_qubit = n_qubit + self._scale = 2 ** n_append + + H, W, C = img.shape + # pad to whole patches, [H, W, C] + img_pad, pads = self.pad(img) + # split channels, C * [H, W] + channels = [img_pad[..., c] for c in range(C)] + channels_processed: list[ndarray] = [] + for channel in channels: + # patchify, [B, p, p] + patches = self.pachify(channel) + # vectorize, [B, d=p**2] + vectors, norm = self.vectorize(patches) + # QFT-append-IQFT (spectral zero-pad), [B, D=P**2] + qvectors = self.apply_unQJPEG(vectors, n_append) + # devectorize, [B, P, P] + devectors = self.devectorize(qvectors, norm) + # unpatchify, [Ho, Wo] + channel_unpatchify = self.unpachify(devectors) + channels_processed.append(channel_unpatchify) + # merge channel, [Ho, Wo, C] + img_merged = np.stack(channels_processed, axis=-1) + # trim pads, [Ho, Wo, C] + img_unpad = self.unpad(img_merged, pads) + # dtype back + img_unpad = img_unpad.astype(orig_dtype) + return img_unpad[..., 0] if is_grey else img_unpad + + def pad(self, img: ndarray): + """Zero-pad an ``[H, W, C]`` image so H and W are whole multiples of + ``patch_size``, pasting the content to the center. Returns the padded image + and the ``(top, bottom, left, right)`` pad widths, or ``None`` if unneeded.""" + H, W, C = img.shape + n_patch_H = math.ceil(H / self.patch_size) + n_patch_W = math.ceil(W / self.patch_size) + H_ex = n_patch_H * self.patch_size + W_ex = n_patch_W * self.patch_size + if H_ex == H and W_ex == W: + return img, None + pad_H = H_ex - H + pad_W = W_ex - W + pads = ( + math.floor(pad_H / 2), math.ceil(pad_H / 2), + math.floor(pad_W / 2), math.ceil(pad_W / 2), + ) + img_ex = np.zeros(shape=(H_ex, W_ex, C), dtype=img.dtype) + img_ex[pads[0]:pads[0] + H, pads[2]:pads[2] + W, :] = img + return img_ex, pads + + def unpad(self, img: ndarray, pads=None) -> ndarray: + """Trim the padding added by :meth:`pad`. The pad widths are scaled *up* by + the up-sampling factor to match the enlarged resolution.""" + if not pads or all(e == 0 for e in pads): + return img + top, bottom, left, right = pads + s = self._scale + H, W, C = img.shape + t = round(top * s) + b = round(bottom * s) + l = round(left * s) + r = round(right * s) + return img[t:(H - b) if b else H, l:(W - r) if r else W, :] + + def pachify(self, img: ndarray) -> ndarray: + """Split a single-channel ``[H, W]`` image into non-overlapping + ``patch_size`` x ``patch_size`` patches, returned as ``[B, p, p]``.""" + P = self.patch_size + H, W = img.shape + assert H % P == 0 and W % P == 0 + nH, nW = H // P, W // P + self._grid = (nH, nW) + return img.reshape(nH, P, nW, P).swapaxes(1, 2).reshape(nH * nW, P, P) + + def unpachify(self, patches: ndarray) -> ndarray: + """Re-assemble ``[B, P, P]`` (up-sampled) patches into a ``[Ho, Wo]`` image.""" + nH, nW = self._grid + B, P, _ = patches.shape + assert B == nH * nW + return patches.reshape(nH, nW, P, P).swapaxes(1, 2).reshape(nH * P, nW * P) + + def vectorize(self, patches: ndarray) -> tuple: + """Vectorize ``[B, p, p]`` patches to normalized amplitude vectors + ``[B, p**2]`` (``amplitude = sqrt(pixel / patch_sum)``) and return the + per-patch intensity ``norm``.""" + B = patches.shape[0] + vect = patches.reshape(B, -1) + norm = vect.sum(axis=1) + amps = np.zeros_like(vect, dtype=patches.dtype) + nz = norm > 0 + amps[nz] = np.sqrt(vect[nz] / norm[nz, None]) + return amps, norm + + def devectorize(self, vectors: ndarray, norm: object) -> ndarray: + """Decode output probability vectors ``[B, P**2]`` into intensity patches + ``[B, P, P]``, rescaled by ``norm`` and the pixel-count ratio so the mean + intensity of each patch is preserved by the up-sampling.""" + B, D = vectors.shape + P = int(round(math.sqrt(D))) + n_in = self.patch_size ** 2 + n_out = D + patches = vectors * norm[:, None] * (n_out / n_in) + return patches.reshape(B, P, P) + + def apply_unQJPEG(self, vectors: ndarray, n_append: int) -> ndarray: + """Run the unQJPEG (spectral zero-pad) circuit on each amplitude vector + ``[B, p**2]`` and return the up-sampled probability vectors ``[B, P**2]``, + with ``P = p * 2 ** n_append``.""" + + ''' + LSB MSB + |-------------nq--------------| + |------ko------|------ko------| + |---ki---| add |---ki---| add | + ''' + nq = self.n_qubit # all qubits + ko = nq // 2 # qubits per axis of the (enlarged) output (aka. k_in) + ki = ko - n_append # qubits per axis of the (init) input (aka. k_out) + qv = list(range(nq)) + # per-axis qubit layout of the enlarged register + col_init = list(range(0, ki)) + col_add = list(range(ki, ko)) + col_msb = col_init[-1] + row_init = list(range(ko, ko + ki)) + row_add = list(range(ko + ki, 2 * ko)) + row_msb = row_init[-1] + # data index bit order: col bits (LSB..MSB) then row bits (LSB..MSB) + enc_qubits = col_init + row_init + measured = qv + + B = vectors.shape[0] + D = 2 ** nq # dim out + out = np.zeros((B, D), dtype=vectors.dtype) + rng = np.random.default_rng(self.seed) if self.shots else None + for b in range(B): + vec = vectors[b] + if not np.any(vec): continue + + prog = QProg(nq) + # data enc + enc = Encode() + enc.amplitude_encode(enc_qubits, list(vec)) + prog << enc.get_circuit() + # QFT + if self.separated_QFT: + # per-axis QFT of the init registers + prog << QFT(col_init) + prog << QFT(row_init) + else: + # QFT over the whole init register + prog << QFT(enc_qubits) + # zero-pad the spectrum: copy the sign qubit onto the appended qubits + for qi in col_add: + prog << CNOT(col_msb, qi) + for qi in row_add: + prog << CNOT(row_msb, qi) + # iQFT + if self.separated_QFT: + # per-axis inverse QFT of the enlarged registers + prog << QFT(list(range(0, ko))).dagger() + prog << QFT(list(range(ko, 2 * ko))).dagger() + else: + # inverse QFT over the whole enlarged register + prog << QFT(qv).dagger() + + machine = CPUQVM() + machine.run(prog, 1) + prob_dict = machine.result().get_prob_dict(measured) + probs = self._prob_dict_to_vector(prob_dict, len(measured), D, vectors.dtype) + if self.shots: + probs = rng.multinomial(self.shots, probs) / self.shots + out[b] = probs + return out + + @staticmethod + def _prob_dict_to_vector(prob_dict: dict, n: int, dim: int, dtype=np.float64) -> ndarray: + """Convert a ``get_prob_dict`` result into a dense vector indexed so that + ``qubits[0]`` is the least-significant bit (matching the QPIE encoding).""" + vec = np.zeros(dim, dtype=dtype) + for key, val in prob_dict.items(): + idx = 0 + for k in range(n): + # key is big-endian over `qubits`: leftmost char == last qubit + if key[n - 1 - k] == '1': + idx |= (1 << k) + vec[idx] += val + return vec diff --git a/pyqpanda-algorithm/pyqpanda_alg/__init__.py b/pyqpanda-algorithm/pyqpanda_alg/__init__.py index 12d6808..61c004c 100644 --- a/pyqpanda-algorithm/pyqpanda_alg/__init__.py +++ b/pyqpanda-algorithm/pyqpanda_alg/__init__.py @@ -51,4 +51,4 @@ from . import Grover from . import QmRMR from . import QSEncode - +from . import QJEPG diff --git a/test/QJEPG/Test_QJPEG.py b/test/QJEPG/Test_QJPEG.py new file mode 100644 index 0000000..21fad47 --- /dev/null +++ b/test/QJEPG/Test_QJPEG.py @@ -0,0 +1,95 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +"""Tests for pyqpanda_alg.QJEPG.QJPEG (quantum JPEG down-sampling).""" + +import sys +from pathlib import Path + +import numpy as np +import pytest + +# make `pyqpanda_alg` importable from the nested package directory +sys.path.insert(0, str(Path(__file__).resolve().parents[2] / 'pyqpanda-algorithm')) +from pyqpanda_alg.QJEPG import QJPEG + + +def _smooth(size: int) -> np.ndarray: + x = np.linspace(0, 2 * np.pi, size) + xx, yy = np.meshgrid(x, x) + return (np.sin(xx) * np.cos(yy) * 0.4 + 0.5).astype(np.float64) + + +def _box_downsample(img: np.ndarray, factor: int) -> np.ndarray: + p = img.shape[0] // factor + return img.reshape(p, factor, p, factor).mean(axis=(1, 3)) + + +class TestQJPEG: + """QJPEG 量子图像下采样测试""" + + def test_output_shape(self): + img = _smooth(16) + out = QJPEG(patch_size=16)(img, n_discard=1) + assert isinstance(out, np.ndarray) + assert out.shape == (8, 8) + + @pytest.mark.parametrize("n_discard,expected", [(1, 8), (2, 4), (3, 2)]) + def test_downsample_factors(self, n_discard, expected): + img = _smooth(16) + out = QJPEG(patch_size=16)(img, n_discard=n_discard) + assert out.shape == (expected, expected) + + def test_matches_box_average(self): + """With the Hadamard layers, QJPEG reproduces classical box-averaging.""" + img = _smooth(16) + out = QJPEG(patch_size=16, wrap_H_layer=True)(img, n_discard=1) + box = _box_downsample(img, 2) + corr = np.corrcoef(out.ravel(), box.ravel())[0, 1] + assert corr > 0.99 + assert np.abs(out - box).mean() < 0.02 + + def test_tiling_multiple_patches(self): + """Image larger than patch_size is tiled and processed per patch.""" + img = _smooth(32) + out = QJPEG(patch_size=16)(img, n_discard=1) + assert out.shape == (16, 16) + + @pytest.mark.parametrize("channels", [3, 4]) # 3 = RGB, 4 = RGB+IR / RGBA + def test_multichannel_image(self, channels): + img = np.stack([_smooth(16) * (1.0 - 0.1 * c) for c in range(channels)], axis=-1) + out = QJPEG(patch_size=16)(img, n_discard=1) + assert out.shape == (8, 8, channels) + # channels are processed independently, each preserving its own mean + for c in range(channels): + assert abs(out[c].mean() - img[c].mean()) < 1e-4 + + def test_reject_unsupported_channel_count(self): + img = np.stack([_smooth(16)] * 5, axis=-1) # 5-channel is not supported + with pytest.raises(AssertionError): + QJPEG(patch_size=16)(img, n_discard=1) + + def test_padding_non_multiple_size(self): + img = _smooth(20) # not a multiple of patch_size=16 + out = QJPEG(patch_size=16)(img, n_discard=1) + # 20 padded to 32, down-sampled by 2, trimmed back to ~20/2 + assert out.shape == (10, 10) + + def test_shots_close_to_exact(self): + img = _smooth(8) + exact = QJPEG(patch_size=8)(img, n_discard=1) + sampled = QJPEG(patch_size=8, shots=100000, seed=0)(img, n_discard=1) + assert sampled.shape == exact.shape + assert np.abs(exact - sampled).max() < 0.05 + + def test_invalid_patch_size(self): + with pytest.raises(AssertionError): + QJPEG(patch_size=15) # not a power of two + + def test_invalid_n_discard(self): + img = _smooth(16) # patch_size 16 -> k = 4, valid n_discard in [1, 3] + with pytest.raises(AssertionError): + QJPEG(patch_size=16)(img, n_discard=4) + + +if __name__ == "__main__": + pytest.main([__file__, "-v", "-s"]) diff --git a/test/QJEPG/Test_unQJPEG.py b/test/QJEPG/Test_unQJPEG.py new file mode 100644 index 0000000..cece19d --- /dev/null +++ b/test/QJEPG/Test_unQJPEG.py @@ -0,0 +1,87 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +"""Tests for pyqpanda_alg.QJEPG.unQJPEG (inverse quantum JPEG up-sampling).""" + +import sys +from pathlib import Path + +import numpy as np +import pytest + +# make `pyqpanda_alg` importable from the nested package directory +sys.path.insert(0, str(Path(__file__).resolve().parents[2] / 'pyqpanda-algorithm')) +from pyqpanda_alg.QJEPG import QJPEG, unQJPEG + + +def _smooth(size: int) -> np.ndarray: + x = np.linspace(0, 2 * np.pi, size) + xx, yy = np.meshgrid(x, x) + return (np.sin(xx) * np.cos(yy) * 0.4 + 0.5).astype(np.float64) + + +def _box_downsample(img: np.ndarray, factor: int) -> np.ndarray: + p = img.shape[0] // factor + return img.reshape(p, factor, p, factor).mean(axis=(1, 3)) + + +class TestUnQJPEG: + """unQJPEG 量子图像上采样测试""" + + def test_output_shape(self): + img = _smooth(8) + out = unQJPEG()(img, n_append=1, patch_size=8) + assert isinstance(out, np.ndarray) + assert out.shape == (16, 16) + + @pytest.mark.parametrize("n_append,expected", [(1, 16), (2, 32)]) + def test_upsample_factors(self, n_append, expected): + img = _smooth(8) + out = unQJPEG()(img, n_append=n_append, patch_size=8) + assert out.shape == (expected, expected) + + def test_flat_stays_flat(self): + """A constant image must up-sample to the same constant, artifact-free.""" + flat = np.full((8, 8), 0.5, dtype=np.float64) + out = unQJPEG()(flat, n_append=1, patch_size=8) + assert out.shape == (16, 16) + assert np.allclose(out, 0.5, atol=1e-6) + + def test_non_negative_and_mean_preserving(self): + img = _smooth(8) + out = unQJPEG()(img, n_append=1, patch_size=8) + assert out.min() >= -1e-9 + # up-sampling preserves the overall mean intensity + assert abs(out.mean() - img.mean()) < 1e-6 + + def test_blockpool_recovers_original(self): + """Down-pooling the up-sampled image returns (close to) the input.""" + img = _smooth(8) + out = unQJPEG()(img, n_append=1, patch_size=8) + rec = _box_downsample(out, 2) + corr = np.corrcoef(rec.ravel(), img.ravel())[0, 1] + assert corr > 0.9 + + def test_round_trip_with_qjpeg(self): + """QJPEG down then unQJPEG up returns an image of the original size.""" + orig = _smooth(16) + low = QJPEG(patch_size=16)(orig, n_discard=1) + assert low.shape == (8, 8) + rec = unQJPEG()(low, n_append=1, patch_size=8) + assert rec.shape == (16, 16) + corr = np.corrcoef(rec.ravel(), orig.ravel())[0, 1] + assert corr > 0.9 + + @pytest.mark.parametrize("channels", [3, 4]) # 3 = RGB, 4 = RGB+IR / RGBA + def test_multichannel_image(self, channels): + img = np.stack([_smooth(8) * (1.0 - 0.1 * c) for c in range(channels)], axis=-1) + out = unQJPEG()(img, n_append=1, patch_size=8) + assert out.shape == (16, 16, channels) + + def test_invalid_patch_size(self): + img = _smooth(8) + with pytest.raises(AssertionError): + unQJPEG()(img, n_append=1, patch_size=6) # not a power of two + + +if __name__ == "__main__": + pytest.main([__file__, "-v", "-s"]) diff --git a/test/pytest.ini b/test/pytest.ini index 2887f64..ecdaaab 100644 --- a/test/pytest.ini +++ b/test/pytest.ini @@ -5,6 +5,7 @@ testpaths = QRAM QPCA QSVM + QJEPG python_files =