add: density noise backend

4 files changed, 463 insertions, 0 deletions

diff --git a/include/core/noise.h b/include/core/noise.h
new file mode 100644
index 0000000..6232d8e
--- /dev/null
+++ b/include/core/noise.h
@@ -0,0 +1,57 @@
+#pragma once
+
+#include <stdbool.h>
+#include <stddef.h>
+
+#include "maths/complex.h"
+#include "maths/matrix.h"
+
+struct PsiKrausOperator
+{
+	struct PsiMatrix matrix;
+	const char *name;
+};
+
+struct PsiKrausOperator psi_new_kraus_operator(const char *name, struct PsiMatrix matrix);
+void psi_free_kraus_operator(struct PsiKrausOperator *op);
+
+struct PsiNoiseChannel
+{
+	const char *name;
+	struct PsiKrausOperator *operators;
+	size_t operator_count;
+	size_t num_qubits;
+};
+
+struct PsiNoiseChannel psi_new_noise_channel(const char *name, const struct PsiKrausOperator *operators, size_t count, size_t num_qubits);
+void psi_free_noise_channel(struct PsiNoiseChannel *channel);
+
+struct PsiNoiseChannel psi_depolarising_channel(double p);
+struct PsiNoiseChannel psi_amplitude_damping_channel(double gamma);
+struct PsiNoiseChannel psi_phase_damping_channel(double gamma);
+struct PsiNoiseChannel psi_bit_flip_channel(double p);
+struct PsiNoiseChannel psi_phase_flip_channel(double p);
+struct PsiNoiseChannel psi_bit_phase_flip_channel(double p);
+struct PsiNoiseChannel psi_generalised_amplitude_damping_channel(double p, double gamma);
+
+struct PsiDensityMatrix
+{
+	struct PsiComplex *data;
+	size_t dim;
+	size_t num_qubits;
+};
+
+struct PsiDensityMatrix psi_new_density_matrix(size_t num_qubits);
+struct PsiDensityMatrix psi_new_density_matrix_from_state(const struct PsiComplex *state, size_t len);
+void psi_free_density_matrix(struct PsiDensityMatrix *dm);
+
+struct PsiComplex psi_get_density_matrix(struct PsiDensityMatrix dm, size_t row, size_t col);
+void psi_set_density_matrix(struct PsiDensityMatrix *dm, size_t row, size_t col, struct PsiComplex value);
+struct PsiComplex psi_trace_density_matrix(struct PsiDensityMatrix dm);
+double psi_purity_density_matrix(struct PsiDensityMatrix dm);
+bool psi_is_pure_density_matrix(struct PsiDensityMatrix dm, double tolerance);
+void psi_density_matrix_probabilities(struct PsiDensityMatrix dm, double *out);
+void psi_apply_unitary_density_matrix(struct PsiDensityMatrix *dm, struct PsiMatrix gate, const size_t *targets, size_t target_count);
+void psi_apply_noise_channel(struct PsiDensityMatrix *dm, struct PsiNoiseChannel channel, size_t target);
+double psi_measure_probability_density_matrix(struct PsiDensityMatrix dm, size_t qubit, size_t outcome);
+double psi_fidelity_density_matrix(struct PsiDensityMatrix dm, const struct PsiComplex *state);
diff --git a/include/psi.h b/include/psi.h
index 45cf8a1..28f2edd 100644
--- a/include/psi.h
+++ b/include/psi.h
@@ -18,3 +18,4 @@ const char *psi_version(void);
 #include "core/circuit.h"
 #include "core/kernel.h"
 #include "core/runtime.h"
+#include "core/noise.h"
diff --git a/src/core/noise.c b/src/core/noise.c
new file mode 100644
index 0000000..9aaf1d9
--- /dev/null
+++ b/src/core/noise.c
@@ -0,0 +1,392 @@
+#include "core/noise.h"
+
+#include <assert.h>
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+
+struct PsiKrausOperator psi_new_kraus_operator(const char *name, struct PsiMatrix matrix)
+{
+	return (struct PsiKrausOperator)
+	{
+		matrix,
+		name,
+	};
+}
+
+void psi_free_kraus_operator(struct PsiKrausOperator *op)
+{
+	psi_free_matrix(&op->matrix);
+}
+
+struct PsiNoiseChannel psi_new_noise_channel(const char *name, const struct PsiKrausOperator *operators, size_t count, size_t num_qubits)
+{
+	struct PsiKrausOperator *owned = malloc(count * sizeof(struct PsiKrausOperator));
+	assert(owned != NULL || count == 0);
+
+	if (count > 0)
+		memcpy(owned, operators, count * sizeof(struct PsiKrausOperator));
+
+	return (struct PsiNoiseChannel)
+	{
+		name,
+		owned,
+		count,
+		num_qubits,
+	};
+}
+
+void psi_free_noise_channel(struct PsiNoiseChannel *channel)
+{
+	for (size_t i = 0; i < channel->operator_count; i++)
+		psi_free_matrix(&channel->operators[i].matrix);
+
+	free(channel->operators);
+	channel->operators = NULL;
+	channel->operator_count = 0;
+}
+
+struct PsiNoiseChannel psi_depolarising_channel(double p)
+{
+	double sqrt_1_p = sqrt(1.0 - p);
+	double sqrt_p3 = sqrt(p / 3.0);
+
+	struct PsiKrausOperator ops[] = {
+		psi_new_kraus_operator("K0", psi_matrix(2, 2,
+			psi_new_complex(sqrt_1_p, 0.0), psi_new_complex(0.0, 0.0),
+			psi_new_complex(0.0, 0.0), psi_new_complex(sqrt_1_p, 0.0))),
+		psi_new_kraus_operator("K1(X)", psi_matrix(2, 2,
+			psi_new_complex(0.0, 0.0), psi_new_complex(sqrt_p3, 0.0),
+			psi_new_complex(sqrt_p3, 0.0), psi_new_complex(0.0, 0.0))),
+		psi_new_kraus_operator("K2(Y)", psi_matrix(2, 2,
+			psi_new_complex(0.0, 0.0), psi_new_complex(0.0, -sqrt_p3),
+			psi_new_complex(0.0, sqrt_p3), psi_new_complex(0.0, 0.0))),
+		psi_new_kraus_operator("K3(Z)", psi_matrix(2, 2,
+			psi_new_complex(sqrt_p3, 0.0), psi_new_complex(0.0, 0.0),
+			psi_new_complex(0.0, 0.0), psi_new_complex(-sqrt_p3, 0.0))),
+	};
+
+	return psi_new_noise_channel("Depolarising", ops, 4, 1);
+}
+
+struct PsiNoiseChannel psi_amplitude_damping_channel(double gamma)
+{
+	double sqrt_gamma = sqrt(gamma);
+	double sqrt_1_gamma = sqrt(1.0 - gamma);
+
+	struct PsiKrausOperator ops[] = {
+		psi_new_kraus_operator("K0", psi_matrix(2, 2,
+			psi_new_complex(1.0, 0.0), psi_new_complex(0.0, 0.0),
+			psi_new_complex(0.0, 0.0), psi_new_complex(sqrt_1_gamma, 0.0))),
+		psi_new_kraus_operator("K1", psi_matrix(2, 2,
+			psi_new_complex(0.0, 0.0), psi_new_complex(sqrt_gamma, 0.0),
+			psi_new_complex(0.0, 0.0), psi_new_complex(0.0, 0.0))),
+	};
+
+	return psi_new_noise_channel("AmplitudeDamping", ops, 2, 1);
+}
+
+struct PsiNoiseChannel psi_phase_damping_channel(double gamma)
+{
+	double sqrt_gamma = sqrt(gamma);
+	double sqrt_1_gamma = sqrt(1.0 - gamma);
+
+	struct PsiKrausOperator ops[] = {
+		psi_new_kraus_operator("K0", psi_matrix(2, 2,
+			psi_new_complex(1.0, 0.0), psi_new_complex(0.0, 0.0),
+			psi_new_complex(0.0, 0.0), psi_new_complex(sqrt_1_gamma, 0.0))),
+		psi_new_kraus_operator("K1", psi_matrix(2, 2,
+			psi_new_complex(0.0, 0.0), psi_new_complex(0.0, 0.0),
+			psi_new_complex(0.0, 0.0), psi_new_complex(sqrt_gamma, 0.0))),
+	};
+
+	return psi_new_noise_channel("PhaseDamping", ops, 2, 1);
+}
+
+struct PsiNoiseChannel psi_bit_flip_channel(double p)
+{
+	double sqrt_1_p = sqrt(1.0 - p);
+	double sqrt_p = sqrt(p);
+
+	struct PsiKrausOperator ops[] = {
+		psi_new_kraus_operator("K0(I)", psi_matrix(2, 2,
+			psi_new_complex(sqrt_1_p, 0.0), psi_new_complex(0.0, 0.0),
+			psi_new_complex(0.0, 0.0), psi_new_complex(sqrt_1_p, 0.0))),
+		psi_new_kraus_operator("K1(X)", psi_matrix(2, 2,
+			psi_new_complex(0.0, 0.0), psi_new_complex(sqrt_p, 0.0),
+			psi_new_complex(sqrt_p, 0.0), psi_new_complex(0.0, 0.0))),
+	};
+
+	return psi_new_noise_channel("BitFlip", ops, 2, 1);
+}
+
+struct PsiNoiseChannel psi_phase_flip_channel(double p)
+{
+	double sqrt_1_p = sqrt(1.0 - p);
+	double sqrt_p = sqrt(p);
+
+	struct PsiKrausOperator ops[] = {
+		psi_new_kraus_operator("K0(I)", psi_matrix(2, 2,
+			psi_new_complex(sqrt_1_p, 0.0), psi_new_complex(0.0, 0.0),
+			psi_new_complex(0.0, 0.0), psi_new_complex(sqrt_1_p, 0.0))),
+		psi_new_kraus_operator("K1(Z)", psi_matrix(2, 2,
+			psi_new_complex(sqrt_p, 0.0), psi_new_complex(0.0, 0.0),
+			psi_new_complex(0.0, 0.0), psi_new_complex(-sqrt_p, 0.0))),
+	};
+
+	return psi_new_noise_channel("PhaseFlip", ops, 2, 1);
+}
+
+struct PsiNoiseChannel psi_bit_phase_flip_channel(double p)
+{
+	double sqrt_1_p = sqrt(1.0 - p);
+	double sqrt_p = sqrt(p);
+
+	struct PsiKrausOperator ops[] = {
+		psi_new_kraus_operator("K0(I)", psi_matrix(2, 2,
+			psi_new_complex(sqrt_1_p, 0.0), psi_new_complex(0.0, 0.0),
+			psi_new_complex(0.0, 0.0), psi_new_complex(sqrt_1_p, 0.0))),
+		psi_new_kraus_operator("K1(Y)", psi_matrix(2, 2,
+			psi_new_complex(0.0, 0.0), psi_new_complex(0.0, -sqrt_p),
+			psi_new_complex(0.0, sqrt_p), psi_new_complex(0.0, 0.0))),
+	};
+
+	return psi_new_noise_channel("BitPhaseFlip", ops, 2, 1);
+}
+
+struct PsiNoiseChannel psi_generalised_amplitude_damping_channel(double p, double gamma)
+{
+	double sqrt_p = sqrt(p);
+	double sqrt_1_p = sqrt(1.0 - p);
+	double sqrt_gamma = sqrt(gamma);
+	double sqrt_1_gamma = sqrt(1.0 - gamma);
+
+	struct PsiKrausOperator ops[] = {
+		psi_new_kraus_operator("K0", psi_matrix(2, 2,
+			psi_new_complex(sqrt_p, 0.0), psi_new_complex(0.0, 0.0),
+			psi_new_complex(0.0, 0.0), psi_new_complex(sqrt_p * sqrt_1_gamma, 0.0))),
+		psi_new_kraus_operator("K1", psi_matrix(2, 2,
+			psi_new_complex(0.0, 0.0), psi_new_complex(sqrt_p * sqrt_gamma, 0.0),
+			psi_new_complex(0.0, 0.0), psi_new_complex(0.0, 0.0))),
+		psi_new_kraus_operator("K2", psi_matrix(2, 2,
+			psi_new_complex(sqrt_1_p * sqrt_1_gamma, 0.0), psi_new_complex(0.0, 0.0),
+			psi_new_complex(0.0, 0.0), psi_new_complex(sqrt_1_p, 0.0))),
+		psi_new_kraus_operator("K3", psi_matrix(2, 2,
+			psi_new_complex(0.0, 0.0), psi_new_complex(0.0, 0.0),
+			psi_new_complex(sqrt_1_p * sqrt_gamma, 0.0), psi_new_complex(0.0, 0.0))),
+	};
+
+	return psi_new_noise_channel("GeneralisedAmplitudeDamping", ops, 4, 1);
+}
+
+struct PsiDensityMatrix psi_new_density_matrix(size_t num_qubits)
+{
+	size_t dim = (size_t)1 << num_qubits;
+	struct PsiComplex *data = calloc(dim * dim, sizeof(struct PsiComplex));
+	assert(data != NULL);
+	data[0] = psi_new_complex(1.0, 0.0);
+
+	return (struct PsiDensityMatrix)
+	{
+		data,
+		dim,
+		num_qubits,
+	};
+}
+
+struct PsiDensityMatrix psi_new_density_matrix_from_state(const struct PsiComplex *state, size_t len)
+{
+	size_t dim = len;
+	size_t num_qubits = 0;
+	while (((size_t)1 << num_qubits) < dim)
+		num_qubits++;
+
+	struct PsiComplex *data = malloc(dim * dim * sizeof(struct PsiComplex));
+	assert(data != NULL);
+
+	for (size_t i = 0; i < dim; i++)
+		for (size_t j = 0; j < dim; j++)
+			data[i * dim + j] = psi_mul_complex(state[i], psi_conjugate_complex(state[j]));
+
+	return (struct PsiDensityMatrix)
+	{
+		data,
+		dim,
+		num_qubits,
+	};
+}
+
+void psi_free_density_matrix(struct PsiDensityMatrix *dm)
+{
+	free(dm->data);
+	dm->data = NULL;
+	dm->dim = 0;
+	dm->num_qubits = 0;
+}
+
+struct PsiComplex psi_get_density_matrix(struct PsiDensityMatrix dm, size_t row, size_t col)
+{
+	assert(row < dm.dim && col < dm.dim);
+	return dm.data[row * dm.dim + col];
+}
+
+void psi_set_density_matrix(struct PsiDensityMatrix *dm, size_t row, size_t col, struct PsiComplex value)
+{
+	assert(row < dm->dim && col < dm->dim);
+	dm->data[row * dm->dim + col] = value;
+}
+
+struct PsiComplex psi_trace_density_matrix(struct PsiDensityMatrix dm)
+{
+	struct PsiComplex sum = psi_new_complex(0.0, 0.0);
+	for (size_t i = 0; i < dm.dim; i++)
+		sum = psi_add_complex(sum, dm.data[i * dm.dim + i]);
+
+	return sum;
+}
+
+double psi_purity_density_matrix(struct PsiDensityMatrix dm)
+{
+	struct PsiComplex sum = psi_new_complex(0.0, 0.0);
+	for (size_t i = 0; i < dm.dim; i++)
+		for (size_t j = 0; j < dm.dim; j++)
+			sum = psi_add_complex(sum, psi_mul_complex(dm.data[i * dm.dim + j], dm.data[j * dm.dim + i]));
+
+	return sum.real;
+}
+
+bool psi_is_pure_density_matrix(struct PsiDensityMatrix dm, double tolerance)
+{
+	return fabs(psi_purity_density_matrix(dm) - 1.0) < tolerance;
+}
+
+void psi_density_matrix_probabilities(struct PsiDensityMatrix dm, double *out)
+{
+	for (size_t i = 0; i < dm.dim; i++)
+		out[i] = dm.data[i * dm.dim + i].real;
+}
+
+void psi_apply_unitary_density_matrix(struct PsiDensityMatrix *dm, struct PsiMatrix gate, const size_t *targets, size_t target_count)
+{
+	size_t g = target_count;
+	size_t gate_dim = (size_t)1 << g;
+	size_t dim = dm->dim;
+
+	size_t *target_bits = malloc(g * sizeof(size_t));
+	assert(target_bits != NULL || g == 0);
+	for (size_t t = 0; t < g; t++)
+		target_bits[t] = dm->num_qubits - 1 - targets[t];
+
+	size_t non_target_mask = dim - 1;
+	for (size_t t = 0; t < g; t++)
+		non_target_mask &= ~((size_t)1 << target_bits[t]);
+
+	struct PsiComplex *new_data = calloc(dim * dim, sizeof(struct PsiComplex));
+	assert(new_data != NULL);
+
+	for (size_t i = 0; i < dim; i++)
+		for (size_t j = 0; j < dim; j++)
+		{
+			struct PsiComplex sum = psi_new_complex(0.0, 0.0);
+
+			for (size_t k = 0; k < gate_dim; k++)
+				for (size_t l = 0; l < gate_dim; l++)
+				{
+					size_t src_i = i & non_target_mask;
+					size_t src_j = j & non_target_mask;
+
+					for (size_t idx = 0; idx < g; idx++)
+					{
+						if ((k >> (g - 1 - idx)) & 1)
+							src_i |= (size_t)1 << target_bits[idx];
+						if ((l >> (g - 1 - idx)) & 1)
+							src_j |= (size_t)1 << target_bits[idx];
+					}
+
+					size_t tgt_i = 0;
+					size_t tgt_j = 0;
+					for (size_t idx = 0; idx < g; idx++)
+					{
+						if ((i >> target_bits[idx]) & 1)
+							tgt_i |= (size_t)1 << (g - 1 - idx);
+						if ((j >> target_bits[idx]) & 1)
+							tgt_j |= (size_t)1 << (g - 1 - idx);
+					}
+
+					struct PsiComplex u_ik = gate.data[tgt_i * gate_dim + k];
+					struct PsiComplex u_jl_dag = psi_conjugate_complex(gate.data[tgt_j * gate_dim + l]);
+					struct PsiComplex rho_kl = dm->data[src_i * dim + src_j];
+
+					sum = psi_add_complex(sum, psi_mul_complex(psi_mul_complex(u_ik, rho_kl), u_jl_dag));
+				}
+
+			new_data[i * dim + j] = sum;
+		}
+
+	free(target_bits);
+	free(dm->data);
+	dm->data = new_data;
+}
+
+void psi_apply_noise_channel(struct PsiDensityMatrix *dm, struct PsiNoiseChannel channel, size_t target)
+{
+	assert(channel.num_qubits == 1);
+
+	size_t dim = dm->dim;
+	size_t target_bit = dm->num_qubits - 1 - target;
+
+	struct PsiComplex *new_data = calloc(dim * dim, sizeof(struct PsiComplex));
+	assert(new_data != NULL);
+
+	for (size_t op = 0; op < channel.operator_count; op++)
+	{
+		struct PsiMatrix k = channel.operators[op].matrix;
+
+		for (size_t i = 0; i < dim; i++)
+			for (size_t j = 0; j < dim; j++)
+			{
+				size_t i_target = (i >> target_bit) & 1;
+				size_t j_target = (j >> target_bit) & 1;
+
+				for (size_t ki = 0; ki < 2; ki++)
+					for (size_t kj = 0; kj < 2; kj++)
+					{
+						size_t src_i = (i & ~((size_t)1 << target_bit)) | (ki << target_bit);
+						size_t src_j = (j & ~((size_t)1 << target_bit)) | (kj << target_bit);
+
+						struct PsiComplex k_elem = k.data[i_target * 2 + ki];
+						struct PsiComplex k_dag_elem = psi_conjugate_complex(k.data[j_target * 2 + kj]);
+						struct PsiComplex rho_elem = dm->data[src_i * dim + src_j];
+
+						struct PsiComplex term = psi_mul_complex(psi_mul_complex(k_elem, rho_elem), k_dag_elem);
+						new_data[i * dim + j] = psi_add_complex(new_data[i * dim + j], term);
+					}
+			}
+	}
+
+	free(dm->data);
+	dm->data = new_data;
+}
+
+double psi_measure_probability_density_matrix(struct PsiDensityMatrix dm, size_t qubit, size_t outcome)
+{
+	size_t target_bit = dm.num_qubits - 1 - qubit;
+	double prob = 0.0;
+
+	for (size_t i = 0; i < dm.dim; i++)
+		if (((i >> target_bit) & 1) == outcome)
+			prob += dm.data[i * dm.dim + i].real;
+
+	return prob;
+}
+
+double psi_fidelity_density_matrix(struct PsiDensityMatrix dm, const struct PsiComplex *state)
+{
+	struct PsiComplex sum = psi_new_complex(0.0, 0.0);
+
+	for (size_t i = 0; i < dm.dim; i++)
+		for (size_t j = 0; j < dm.dim; j++)
+			sum = psi_add_complex(sum, psi_mul_complex(psi_mul_complex(psi_conjugate_complex(state[i]), dm.data[i * dm.dim + j]), state[j]));
+
+	return sum.real;
+}
diff --git a/tester/main.c b/tester/main.c
index 346b16a..858be2e 100644
--- a/tester/main.c
+++ b/tester/main.c
@@ -215,5 +215,18 @@ int main(void)
 		psi_format_probability(probs[7], pbuf, sizeof pbuf));
 
 	psi_free_quantum_circuit(&ghz);
+
+	struct PsiDensityMatrix dm = psi_new_density_matrix(1);
+	printf("dm      = purity %g (pure? %d)\n", psi_purity_density_matrix(dm), psi_is_pure_density_matrix(dm, 1e-10));
+
+	struct PsiNoiseChannel bf = psi_bit_flip_channel(0.25);
+	psi_apply_noise_channel(&dm, bf, 0);
+
+	double dprobs[2];
+	psi_density_matrix_probabilities(dm, dprobs);
+	printf("noisy   = p [%g, %g], purity %g\n", dprobs[0], dprobs[1], psi_purity_density_matrix(dm));
+
+	psi_free_noise_channel(&bf);
+	psi_free_density_matrix(&dm);
 }