main.c File Reference

This chunk of source code controls the NST using functions from utilities.c ; function prototypes are in header.h. More...

#include <stdio.h>
#include <unistd.h>
#include <math.h>
#include <time.h>
#include <sys/time.h>
#include <string.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_multimin.h>
#include "header.h"

Go to the source code of this file.

Functions
void	apply_noise (gsl_vector *n)
	Applies Poisson shot noise to measurements stored in n.
void	fill_tangle_entropy_plane ()
	Routine which is used to fill the entire entropy-tangle plane for 2 qubits and perform all tomography routines on each state, this is described in detail in the paper.
void	find_N ()
	Locates the first lowest value of N (the number of times the experiment has to be performed in real life) for tomography to yield a sufficiently high fidelity value in practice.
void	init ()
	Main initialisation function that loads all the data into memory.
void	linear_estimate (gsl_vector *n, gsl_matrix_complex *density_linear)
	Performs linear reconstruction of the underlying density matrix into density_linear from measurements in n.
void	load_basis ()
	Initialisation function which loads hard-coded values of projection and Pauli operators into memory.
int	main (void)
	Main function of the program.
void	measure_state (gsl_matrix_complex *density_matrix, gsl_vector *n)
	Simulates experimental measurements of a particular state, characterised by density_matrix, writes measurement outcomes into pre-allocated n.
void	MLE (gsl_vector *n, gsl_matrix_complex *density_physical, gsl_matrix_complex *density_linear, gsl_matrix_complex *density_mle, gsl_matrix *fidelities, gsl_matrix *times, gsl_matrix *f_vals, gsl_matrix *abs_gradient, int trial)
	Performs MLE tomography on measurements stored in n, output in density_mle. See detailed docs for other arguments.
void	my_df (const gsl_vector *x, void *params, gsl_vector *g)
	Analytic gradient of the MLE function, also used by GSL optimisation routine, refer to GNU manual for parameter meaning.
double	my_f (const gsl_vector *x, void *params)
	MLE function used by GSL for optimisation, refer to GNU manual for this.
void	my_fdf (const gsl_vector *x, void *params, double *f, gsl_vector *g)
	MLE function and analytic gradient computed together, this saves time so that GSL doesn't have to call each function separately, refer to GSL manual for parameter meaning.
void	prepare_state (gsl_matrix_complex *density_matrix)
	Responsible for directing which particular tomography to simulate and how. This is guided by the STATE global variable selector in header.h.
void	quick_and_dirty (gsl_matrix_complex *density_linear, gsl_matrix_complex *density_quick)
	Quick and Dirty Procedure as discussed in the paper, where eigenvalues of density_linear are corrected into density_quick.
void	quick_and_dirty_2 (gsl_matrix_complex *density_linear, gsl_matrix_complex *density_quick)
	Forced Purity routine as discussed in the paper, the naming is somewhat misleading, here density_quick is the resultant FP density matrix.
void	states_1_to_10 ()
	Groups anything that has to do with values 1-10 of STATE global variable described in header.h.

---

Detailed Description

This chunk of source code controls the NST using functions from utilities.c ; function prototypes are in header.h.

Numerical State Tomography Using GSL -> www.gnu.org/software/gsl Copyright (C) 2007 Max S Kaznady

This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation;

AUTHOR Max S Kaznady, [EMAIL PROTECTED] max.kaznady@gmail.com Summer, 2007

Definition in file main.c.

---

Function Documentation

void init

(

)

Main initialisation function that loads all the data into memory.

scale parameter - number of times experiment is performed

Definition at line 1835 of file main.c.

References kron(), load_basis(), mu_array, mu_cell, N, NQ, one, ones, sig_array, sigma_cell, and zero.

Referenced by main().

       {
    //gsl_matrix_complex *mu_array = load_basis();
    load_basis();
    //init constants
    GSL_SET_COMPLEX(&zero, 0, 0);
    GSL_SET_COMPLEX(&one, 1, 0);
    ones = gsl_matrix_complex_alloc(pow(2, NQ), pow(2, NQ));
    gsl_matrix_complex_set_all(ones, one);
    //allocate mu_cell
    int i;
    mu_cell[0] = gsl_matrix_complex_alloc(2, 2);
    gsl_matrix_complex_memcpy(mu_cell[0], mu_array[0]);
    for (i=1; i<NQ; i++) {
        mu_cell[i] = gsl_matrix_complex_alloc((mu_cell[i-1]->size1)*2, (mu_cell[i-1]->size2)*2);
        //printf("entering kron\n");
        //printf("rows mu_cell[i-1] : rows %d\n", mu_cell[i-1]->size1);
        //printf("rows mu_array[0] : rows %d\n", mu_array[0]->size1);
        kron(mu_cell[i-1], mu_array[0], mu_cell[i]);
    }

    //initialize sigma_cell
    sigma_cell[0] = gsl_matrix_complex_alloc(2,2);
    gsl_matrix_complex_memcpy(sigma_cell[0], sig_array[0]);
    for (i=1; i<NQ; i++) {
        sigma_cell[i] = gsl_matrix_complex_alloc((sigma_cell[i-1]->size1)*2, (sigma_cell[i-1]->size2)*2);
        kron(sigma_cell[i-1], sig_array[0], sigma_cell[i]);
    }

    //! scale parameter - number of times experiment is performed
    N = 1000000;
}

void MLE	(	gsl_vector *	n,
		gsl_matrix_complex *	density_physical,
		gsl_matrix_complex *	density_linear,
		gsl_matrix_complex *	density_mle,
		gsl_matrix *	fidelities,
		gsl_matrix *	times,
		gsl_matrix *	f_vals,
		gsl_matrix *	abs_gradient,
		int	trial
	)

Performs MLE tomography on measurements stored in n, output in density_mle. See detailed docs for other arguments.

Detailed description:
n - vector containing experimental measurements.
density_physical - actual density matrix of the underlying state with state error applied to it, used to compute fidelity.
density_linear - entries of this matrix are used as an initial guess, the name is misleading, because the argument in the calling program is actually density_quick, so we are in fact using the values of the Quick and Dirty routine as a starting point for MLE.
density_mle - density matrix obtained using MLE is stored here.
fidelities - stores fidelities between density_physical and density_mle at each step of the iteration as it progresses.
times - vector of runtimes for each iteration.
f_vals - values of the MLE function for each iteration.
abs_gradient - absolute value of the gradient of the MLE function for each iteration.
trial - number of trials to be performed.

Definition at line 1166 of file main.c.

References abs_val(), CUT_ITER, fidelity(), KEEP_GOING, make_T(), MAX_FIDELITY, MAXITER, MILLION, my_df(), my_f(), my_fdf(), NQ, one, trace(), VERBOSE, and zero.

Referenced by fill_tangle_entropy_plane(), and states_1_to_10().

                                                                                                                                                                                                                                  {

    //first perform Cholesky Decomposition
    // FIX THIS! use a guess for now
    gsl_vector *x = gsl_vector_alloc(pow(4,NQ));
    //gsl_vector_set_all(x, -100);

    int counter, limiter, nu;
    counter = 1;
    for (limiter = 0; limiter < pow(2, NQ); limiter++) {
        for (nu=1+limiter; nu <= pow(2, NQ); nu++) {
            if (counter<=pow(2, NQ)) {
                gsl_vector_set(x, counter-1, GSL_REAL(gsl_matrix_complex_get(density_linear, nu-1, nu-1)));
            } else {
                gsl_vector_set(x, counter-1, GSL_REAL(gsl_matrix_complex_get(density_linear, nu-1, nu-1-limiter)));
                gsl_vector_set(x, counter, GSL_IMAG(gsl_matrix_complex_get(density_linear, nu-1, nu-1-limiter)));
                //nu++;
                counter++;

            }
            counter++;
        }
    }
    //display_vector_real(x);
    //display_matrix(density_linear);

    //   for (nu=0; nu<pow(4,NQ); nu++)
    //     gsl_vector_set(x, nu, abs(gsl_vector_get(x,nu)));

    //gsl_vector_set_all(x, -100);
    gsl_vector_set(x, 0, pow(4, -NQ));

    //initialize the minimizer
    const gsl_multimin_fdfminimizer_type *T;
    gsl_multimin_fdfminimizer *s;
    //T = gsl_multimin_fdfminimizer_conjugate_fr; //minimizer type
    //T = gsl_multimin_fdfminimizer_conjugate_pr;
    //T = gsl_multimin_fdfminimizer_steepest_descent;
    //T = gsl_multimin_fdfminimizer_vector_bfgs;
    T = gsl_multimin_fdfminimizer_vector_bfgs2;

    s = gsl_multimin_fdfminimizer_alloc(T, pow(4,NQ));
    gsl_multimin_function_fdf my_func;
    size_t iter = 0;
    int status;

    my_func.f = &my_f;
    my_func.df = &my_df;
    my_func.fdf = &my_fdf;
    //my_func.fdf = NULL;
    my_func.n = pow(4,NQ);
    my_func.params = n;

    //   Function: int gsl_multimin_fdfminimizer_set (gsl_multimin_fdfminimizer * s, gsl_multimin_function_fdf * fdf, const gsl_vector * x, double step_size, double tol)
    //
    //   This function initializes the minimizer s to minimize the function fdf starting from the initial point x. The size of the first trial step is given by step_size. The accuracy of the line minimization is specified by tol. The precise meaning of this parameter depends on the method used. Typically the line minimization is considered successful if the gradient of the function g is orthogonal to the current search direction p to a relative accuracy of tol, where dot(p,g) < tol |p| |g|. A tol value of 0.1 is suitable for most purposes, since line minimization only needs to be carried out approximately. Note that setting tol to zero will force the use of “exact” line-searches, which are extremely expensive. — Function: int gsl_multimin_fminimizer_set (gsl_multimin_fminimizer * s, gsl_multimin_function * f, const gsl_vector * x, const gsl_vector * step_size)

    gsl_multimin_fdfminimizer_set (s, &my_func, x, 0.01, 1e-4);
    //gsl_multimin_fdfminimizer_set (s, &my_func, x, 1e-15, 1e-15);
    //gsl_multimin_fdfminimizer_set (s, &my_func, x, 0.0001, 1e-5);
    //gsl_multimin_fdfminimizer_set (s, &my_func, x, 0.0001, 0);
    //gsl_multimin_fdfminimizer_set (s, &my_func, x, pow(1,-10^(NQ+4)), 1e-4);

    int truth=1;

    //temporary MLE matrix
    gsl_matrix_complex *temp_mle = gsl_matrix_complex_alloc(pow(2, NQ), pow(2, NQ));
    //temporary T matrix, used to construct temporary MLE matrix
    gsl_matrix_complex *T_mat = gsl_matrix_complex_alloc(pow(2, NQ), pow(2, NQ));

    //init a timer
    //struct timeval tpstart, tpend;
    //double timedif;

    struct timeval tpstart_iter, tpend_iter;
    double timedif_iter;

    double temp_fid;

    //allocate previous fidelity value - will be overwritten later
    double previous_fidelity = -999.0;
    double previous_f_value = -1.0; //for first iteration
    gsl_vector *previous_t = gsl_vector_alloc(pow(4, NQ));

    do {

        gettimeofday(&tpstart_iter, NULL);

        status = gsl_multimin_fdfminimizer_iterate (s);
        //printf("HERE\n");

        if (status)
            break;

        status = gsl_multimin_test_gradient (s->gradient, 1e-3);
        //status = gsl_multimin_test_gradient (s->gradient, 0);

        if (status == GSL_SUCCESS)
            printf ("Minimum found.\n");

        if (VERBOSE)
            printf("%d: f() = %7.3f \n", iter, s->f);

        //end timer
        gettimeofday(&tpend_iter, NULL);
        timedif_iter = (double) MILLION*(tpend_iter.tv_sec - tpstart_iter.tv_sec) + tpend_iter.tv_usec - tpstart_iter.tv_usec;

        //construct temporaty density matrix, check its fidelity
        make_T(s->x, T_mat);
        gsl_blas_zgemm(CblasConjTrans, CblasNoTrans, one, T_mat, T_mat, zero, temp_mle);
        gsl_matrix_complex_scale(temp_mle, gsl_complex_rect(1.0/GSL_REAL(trace(temp_mle)), 0));

        //calculate temporary fidelity
        temp_fid = fidelity(temp_mle, density_physical);
        //printf("fidelity : %e\n", temp_fid);

        //termination criteria
        if (temp_fid >= MAX_FIDELITY && KEEP_GOING==0)
            truth = 0;

        //stop minimizing if func value difference is less than 1
        if (CUT_ITER == 1) {
            if (iter>0) { //make sure this is not the first iteration
                if ((previous_f_value - s->f) <= 1.0) {
                    truth = 0;
                }
            }
        }
        //update function value
        previous_f_value = s->f;

        if (KEEP_GOING == 0) { //if we want to exit at highest fidelity
            //fidelity is starting to decrease
            //so stop the iteration, return to previous fidelity state
            if (previous_fidelity > temp_fid) {
                truth = 0;
                gsl_vector_memcpy(s->x, previous_t);
            }
        }

        //update fidelity and t_vector
        previous_fidelity = temp_fid;
        gsl_vector_memcpy(previous_t, s->x);

        //Record outcomes for this iteration
        //if (trial<NUMTRIALS) {
        gsl_matrix_set(fidelities, trial, iter, temp_fid);
        gsl_matrix_set(times, trial, iter, timedif_iter);
        gsl_matrix_set(f_vals, trial, iter, s->f);
        gsl_matrix_set(abs_gradient, trial, iter, abs_val(s->gradient));
        //}

        //another termination criteria
        if ((s->f) == GSL_POSINF) {
            printf("Either no point in using MLE or we have a bad guess\n");
            truth = 0;
        }

        //update iteration number
        iter++;

    } while (status == GSL_CONTINUE && truth && iter<MAXITER);

    if (!VERBOSE)
        printf("%d: f() = %7.3f \n", iter, s->f);

    //construct temporaty density matrix, check its fidelity
    make_T(s->x, T_mat);
    gsl_blas_zgemm(CblasConjTrans, CblasNoTrans, one, T_mat, T_mat, zero, density_mle);
    gsl_matrix_complex_scale(density_mle, gsl_complex_rect(1.0/GSL_REAL(trace(density_mle)), 0));

    //deallocate memory
    gsl_multimin_fdfminimizer_free(s);
    gsl_vector_free(x);
    gsl_vector_free(previous_t);
    gsl_matrix_complex_free(T_mat);
    gsl_matrix_complex_free(temp_mle);

}