diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/core/noise.c | 392 |
1 files changed, 392 insertions, 0 deletions
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; +} |