main.c

Go to the documentation of this file.

/*! \file main.c
    \brief This chunk of source code controls the NST using functions from \e utilities.c ; function prototypes are in \e 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
 *
 *
 ******************************************************************************/

/*! \fn void load_basis();
 *  \brief Initialisation function which loads hard-coded values of projection and Pauli operators into memory.
 */

/*! \fn void init();
 *  \brief Main initialisation function that loads all the data into memory.
 */

/*! \fn void prepare_state(gsl_matrix_complex *density_matrix);
 *  \brief Responsible for directing which particular tomography to simulate and how. This is guided by the STATE global variable selector in header.h.
 */

/*! \fn void measure_state(gsl_matrix_complex *density_matrix, gsl_vector *n)
 *  \brief Simulates experimental measurements of a particular state, characterised by \e density_matrix, writes measurement outcomes into pre-allocated \e n.
 */

/*! \fn void apply_noise(gsl_vector *n);
 *  \brief Applies Poisson shot noise to measurements stored in \e n.
 */

/*! \fn void linear_estimate(gsl_vector *n, gsl_matrix_complex *density_linear);
 *  \brief Performs linear reconstruction of the underlying density matrix into \e density_linear from measurements in \e n.
 */

/*! \fn 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);
 *  \brief Performs MLE tomography on measurements stored in \e n, output in \e density_mle. See detailed docs for other arguments.
 *
 *  Detailed description: \n
 *  n - vector containing experimental measurements. \n
 *  density_physical - actual density matrix of the underlying state with state error applied to it, used to compute fidelity. \n
 *  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 \e density_quick, so we are in fact using the values of the Quick and Dirty routine as a starting point for MLE. \n
 *  density_mle - density matrix obtained using MLE is stored here. \n
 *  fidelities - stores fidelities between density_physical and density_mle at each step of the iteration as it progresses. \n
 *  times - vector of runtimes for each iteration. \n
 *  f_vals - values of the MLE function for each iteration. \n
 *  abs_gradient - absolute value of the gradient of the MLE function for each iteration. \n
 *  trial - number of trials to be performed.
 *
 *
 */

/*! \fn void quick_and_dirty(gsl_matrix_complex *density_linear, gsl_matrix_complex *density_quick);
 *  \brief Quick and Dirty Procedure as discussed in the paper, where eigenvalues of \e density_linear are corrected into \e density_quick.
 */

/*! \fn void quick_and_dirty_2(gsl_matrix_complex *density_linear, gsl_matrix_complex *density_quick);
 *  \brief Forced Purity routine as discussed in the paper, the naming is somewhat misleading, here \e density_quick is the resultant FP density matrix.
 */

/*! \fn double my_f(const gsl_vector * x, void * params);
 *  \brief MLE function used by GSL for optimisation, refer to GNU manual for this.
 */

/*! \fn void my_df(const gsl_vector * x, void * params, gsl_vector * g);
 *  \brief Analytic gradient of the MLE function, also used by GSL optimisation routine, refer to GNU manual for parameter meaning.
 */

/*! \fn  void my_fdf(const gsl_vector * x, void * params, double * f, gsl_vector * g);
 *  \brief 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.
 */


/*! \fn void fill_tangle_entropy_plane();
 *  \brief 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.
 */

/*! \fn void states_1_to_10();
 *  \brief Groups anything that has to do with values 1-10 of STATE global variable described in header.h.
 */

/*! \fn void find_N();
 *  \brief 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.
 */

/*! \fn int main(void);
 *  \brief Main function of the program.
 */

/*! \mainpage Numerical State Tomography (NST)
 *
 *
 *
 * \section intro_sec Introduction
 *
 * This software code is to accompany the paper "Quantum State Tomography: Is the 'best' the enemy of 'good enough'?"
 * All the numerics in the paper were performed using this source code. \n
 *
 * Click on the "Files" tab at the top left-hand corner to view the detailed documentation of the source code. \n
 *
 * Source code is available here: http://www.physics.utoronto.ca/~dfvj/publications/NST.tar.gz
 *
 * \section install_sec Installation
 *
 * \subsection step1 Step 1: Pre-requisites
 *
 * This software uses GNU Scientific Library version 1.9, which can be obtained here: ftp://ftp.gnu.org/gnu/gsl/gsl-1.9.tar.gz \n
 *
   We used GNU C compiler of various versions from the GNU compiler collection, but you're welcome to use any other compiler.
 *
 * Unpack and install the archive as per GSL Instructions, basic steps would be: \n
 * \e ./configure \e --prefix=/usr/local/gsl-1.9 \n
 * \e make \n
 * \e sudo \e make \e install \n
 *  last command "make install" has to be run as root \n
 *
 * You can install GSL anywhere you like, but for the sake of simplicity I'll refer to GSL's location as /usr/local/gsl-1/9.
 *
 * \subsection step2 Step 2: Modifying local system files
 *
 * If you installation location is different to /usr/local/gsl-1.9, then you also need to modify the Makefile, see the next section. \n
 *
 * If not, then you need to export an environment variable LD_LIBRARY_PATH as described here: http://www.gnu.org/software/gsl/manual/html_node/Shared-Libraries.html \n
 *
 * In short, add these lines to ~/.bashrc if you have BASH Shell:
 *
 * LD_LIBRARY_PATH=/usr/local/gsl-1.9/lib:$LD_LIBRARY_PATH \n
   export LD_LIBRARY_PATH \n
 
   Now restart the terminal, change directory to the unpacked NST source code and run "make". You may need to modify the Makefile (see next section). \n
 
   To change the parameters of the simulation edit the file "header.h" in the source code and recompile. Click on "Files" tab and read the documentation of "header.h" to set the parameters correctly.
 
   \subsection step3 Step 3: Modifying Makefile
 
   Locate the lines starting with LIBC and INCL and change the paths to point to the location of your GSL. \n
 
   The line COMPILER should also be set to the proper compiler, but most Unix/Linux distribution come with GCC installed by default.
 *
 * \section license_sec License
 *
 *  Using GSL -> www.gnu.org/software/gsl \n
 *  Copyright (C) 2007  Max S Kaznady \n
 *
 *  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; \n
 *
 *  AUTHOR \n
 *     Max S Kaznady,  [EMAIL PROTECTED] \n
 *     max.kaznady@gmail.com \n
 *     Summer, 2007 \n
 */

#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 <gsl/gsl_siman.h> /* annealing */

#include "header.h"

// /* set up parameters for this simulated annealing run */
// /* how many points do we try before stepping */
// #define N_TRIES 200
// /* how many iterations for each T? */
// #define ITERS_FIXED_T 1000
// /* max step size in random walk */
// #define STEP_SIZE 1.0
// /* Boltzmann constant */
// #define K 1.0
// /* initial temperature */
// #define T_INITIAL 0.008
// /* damping factor for temperature */
// #define MU_T 1.003
// #define T_MIN 2.0e-6

/* function prototypes */

void
load_basis();
void
init();
void
prepare_state(gsl_matrix_complex *density_matrix);
void
measure_state(gsl_matrix_complex *density_matrix, gsl_vector *n);
void
apply_noise(gsl_vector *n);
void
linear_estimate(gsl_vector *n, gsl_matrix_complex *density_linear);
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);
void
quick_and_dirty(gsl_matrix_complex *density_linear, gsl_matrix_complex *density_quick);
void
quick_and_dirty_2(gsl_matrix_complex *density_linear, gsl_matrix_complex *density_quick);

double
my_f(const gsl_vector * x, void * params);
void
my_df(const gsl_vector * x, void * params, gsl_vector * g);
void
my_fdf(const gsl_vector * x, void * params, double * f, gsl_vector * g);

void
fill_tangle_entropy_plane();
void
states_1_to_10();

void
find_N();

int
main(void) {

    //allocate global variables
    init();

    if (STATE == 11) {
        if (NQ == 2) {
            fill_tangle_entropy_plane();
        }
        else {
            fprintf(stderr, "Only defined for 2 qubits\n");
            exit(1);
        }
    } else if (STATE == 15) {
        find_N();
    } else { // select other states
        states_1_to_10();
    }

    return 0;
}


//Used to perform 2000000 trials of biased density matrices
void
fill_tangle_entropy_plane() {
    int numtrials = 99;
    int i, j, k;

    //pseudo-random density matrix
    gsl_matrix_complex *density_physical = gsl_matrix_complex_alloc(pow(2,NQ), pow(2,NQ));
    //its copy which keeps the state - theoretical density_physical
    gsl_matrix_complex *density_theoretical = gsl_matrix_complex_alloc(pow(2,NQ), pow(2,NQ));
    //measurement vector
    gsl_vector *n = gsl_vector_alloc(pow(4, NQ));

    //declare and allocate all matrices here
    gsl_matrix_complex *density_linear = gsl_matrix_complex_alloc(pow(2, NQ), pow(2, NQ));
    gsl_matrix_complex *density_quick = gsl_matrix_complex_alloc(pow(2,NQ), pow(2,NQ));
    gsl_matrix_complex *density_fp = gsl_matrix_complex_alloc(pow(2,NQ), pow(2,NQ));
    gsl_matrix_complex *density_mle = gsl_matrix_complex_alloc(pow(2, NQ), pow(2, NQ));
    gsl_matrix_complex *psi_outer = gsl_matrix_complex_alloc(pow(2, NQ), pow(2, NQ));
    gsl_vector_complex *psi_ket = gsl_vector_complex_alloc(pow(2, NQ));

    //record data analysis
    gsl_vector *f_mle, *f_quick, *f_fp,  *tau, *sl;
    //vectors that store fidelities and times for each iteration
    gsl_matrix *fidelities, *times, *f_vals, *abs_gradient;
    //pre-allocate vectors - used to store data
    //int num_records = 11000;
    int num_records = 2*MILLION;
    f_mle = gsl_vector_alloc(num_records);
    f_quick = gsl_vector_alloc(num_records);
    f_fp = gsl_vector_alloc(num_records);
    tau = gsl_vector_alloc(num_records);
    sl = gsl_vector_alloc(num_records);
    printf("HERE\n");
    //fake data structures below
    fidelities = gsl_matrix_alloc(1, MAXITER);
    times = gsl_matrix_alloc(1, MAXITER);
    f_vals = gsl_matrix_alloc(1, MAXITER);
    abs_gradient = gsl_matrix_alloc(1, MAXITER);

    gsl_vector_set_zero(f_mle);
    gsl_vector_set_zero(f_quick);
    gsl_vector_set_zero(f_fp);
    gsl_vector_set_zero(tau);
    gsl_vector_set_zero(sl);
    gsl_matrix_set_all(fidelities, GSL_NAN);
    gsl_matrix_set_all(times, GSL_NAN);
    gsl_matrix_set_all(f_vals, GSL_NAN);
    gsl_matrix_set_all(abs_gradient, GSL_NAN);

    //constants used to save density matrices at each iteration
    char output[MAXCHAR];
    char integer[MAXCHAR];

    int counter = 0;

    double eps_squared;
    double tangle_eps;

    //gsl_matrix_complex_set_zero(psi_outer);

    for (i=0; i<=numtrials; i++) { //control tangle
        tangle_eps = (double)i/numtrials/sqrt(2.0); //range 0 to Pi/2
        //theta = (double)i/numtrials;
        for (j=0; j<=numtrials; j++) { //control entropy
            for (k=0; k<=numtrials; k++) { //loop for each location
                random_density_matrix(density_theoretical);

                eps_squared = gsl_pow_2((double)j/numtrials);

                //restore Entangled State
                PSI_matrix(psi_outer, tangle_eps);

                gsl_matrix_complex_scale(density_theoretical, gsl_complex_rect(eps_squared, 0.0));
                gsl_matrix_complex_scale(psi_outer, gsl_complex_rect(1.0 - eps_squared, 0.0));
                gsl_matrix_complex_add(density_theoretical, psi_outer);

                //record position in H map
                gsl_vector_set(tau, counter, pow(concurrence(density_theoretical), 2));
                gsl_vector_set(sl, counter, entropy_linear(density_theoretical));
                printf("Tau: %e\n", gsl_vector_get(tau, counter));
                printf("Entropy Linear: %e\n", gsl_vector_get(sl, counter));

                //restore theoretical state
                gsl_matrix_complex_memcpy(density_physical, density_theoretical);

                //simulate imperfect state - now we have rho_physical
                printf("Measuring the state\n");
                measure_state(density_physical, n);

                //scale measurement vector here
                gsl_vector_scale(n, N);

                if (WRITE_DENSITY) {
                    //save imperfect state - we will try to reconstruct this state
                    printf("Writing data to disc - density_physical\n");
                    //save the result
                    strncpy(output, "density_physical_trial_", MAXCHAR);
                    snprintf(integer, MAXCHAR, "%d", counter);
                    strncat(output, integer, MAXCHAR);
                    write_data(output, density_physical);
                }

                //apply measurement noise
                printf("Applying experimental noise to data\n");
                apply_noise(n);

                printf("Performing linear reconstruction\n");
                //perform linear reconstruction
                linear_estimate(n, density_linear);

                if (WRITE_DENSITY) {
                    printf("Writing data to disc - density_linear\n");
                    strncpy(output, "density_linear_trial_", MAXCHAR);
                    snprintf(integer, MAXCHAR, "%d", counter);
                    strncat(output, integer, MAXCHAR);
                    write_data(output, density_linear);
                }

                //save density_linear
                gsl_matrix_complex_memcpy(density_fp, density_linear);

                //density_linear will get modified in the process
                printf("Performing Quick and Dirty\n");
                quick_and_dirty(density_linear, density_quick);

                if (WRITE_DENSITY) {
                    printf("Writing data to disc - density_quick\n");
                    //write_data("density_quick", density_quick);
                    strncpy(output, "density_quick_trial_", MAXCHAR);
                    snprintf(integer, MAXCHAR, "%d", counter);
                    strncat(output, integer, MAXCHAR);
                    write_data(output, density_quick);
                }

                //record fidelity
                gsl_vector_set(f_quick, counter, fidelity(density_quick, density_physical));
                printf("Fidelity Quick: %e\n", gsl_vector_get(f_quick, counter));

                //restore density_linear
                gsl_matrix_complex_memcpy(density_linear, density_fp);

                //density_linear will get modified in the process
                printf("Performing Forced Purity\n");
                quick_and_dirty_2(density_linear, density_fp);

                if (WRITE_DENSITY) {
                    printf("Writing data to disc - density_fp\n");
                    strncpy(output, "density_fp_trial_", MAXCHAR);
                    snprintf(integer, MAXCHAR, "%d", counter);
                    strncat(output, integer, MAXCHAR);
                    write_data(output, density_fp);
                }

                //record fidelity
                gsl_vector_set(f_fp, counter, fidelity(density_fp, density_physical));
                printf("Fidelity Forced Purity: %e\n", gsl_vector_get(f_fp, counter));

                printf("Performing MLE\n");
                MLE(n, density_physical, density_quick, density_mle, fidelities, times, f_vals, abs_gradient, 0);

                if (WRITE_DENSITY) {
                    printf("Writing data to disc - density_mle\n");
                    //last density_mle is written each time as "density_mle"
                    write_data("density_mle", density_mle);
                    strncpy(output, "density_mle_trial_", MAXCHAR);
                    snprintf(integer, MAXCHAR, "%d", counter);
                    strncat(output, integer, MAXCHAR);
                    write_data(output, density_mle);
                }

                gsl_vector_set(f_mle, counter, fidelity(density_mle, density_physical));
                printf("Fidelity MLE: %e\n", gsl_vector_get(f_mle, counter));

                counter++;
            }
        }
    }

//now perform MEMS state

    //numtrials = 99;
    numtrials = 999;

    gsl_matrix_complex *rho_mems = gsl_matrix_complex_alloc(4, 4);

    for (i = 0; i <= numtrials; i++) {

        //display_matrix(rho_mems);

        for (j = 0; j <= numtrials; j++) {

            //restore MEM State
            rho_MEMS(rho_mems, (double) i/numtrials);

            random_density_matrix(density_theoretical);

            //display_matrix(density_theoretical);

            eps_squared = pow((double)j/numtrials, 2);

            gsl_matrix_complex_scale(density_theoretical, gsl_complex_rect(eps_squared, 0));
            gsl_matrix_complex_scale(rho_mems, gsl_complex_rect(1 - eps_squared, 0));
            gsl_matrix_complex_add(density_theoretical, rho_mems);
            //display_matrix(rho_mems);

            //gsl_matrix_complex_memcpy(density_theoretical, rho_mems);

            //record position in H map
            gsl_vector_set(tau, counter, pow(concurrence(density_theoretical), 2));
            gsl_vector_set(sl, counter, entropy_linear(density_theoretical));
            printf("Tau: %e\n", gsl_vector_get(tau, counter));
            printf("Entropy Linear: %e\n", gsl_vector_get(sl, counter));

            //restore theoretical state
            gsl_matrix_complex_memcpy(density_physical, density_theoretical);

            //simulate imperfect state - now we have rho_physical
            printf("Measuring the state\n");
            measure_state(density_physical, n);

            //scale measurement vector here
            gsl_vector_scale(n, N);

            if (WRITE_DENSITY) {
                //save imperfect state - we will try to reconstruct this state
                printf("Writing data to disc - density_physical\n");
                //save the result
                strncpy(output, "density_physical_trial_", MAXCHAR);
                snprintf(integer, MAXCHAR, "%d", counter);
                strncat(output, integer, MAXCHAR);
                write_data(output, density_physical);
            }

            //apply measurement noise
            printf("Applying experimental noise to data\n");
            apply_noise(n);

            printf("Performing linear reconstruction\n");
            //perform linear reconstruction
            linear_estimate(n, density_linear);

            if (WRITE_DENSITY) {
                printf("Writing data to disc - density_linear\n");
                strncpy(output, "density_linear_trial_", MAXCHAR);
                snprintf(integer, MAXCHAR, "%d", counter);
                strncat(output, integer, MAXCHAR);
                write_data(output, density_linear);
            }

            //save density_linear
            gsl_matrix_complex_memcpy(density_fp, density_linear);

            //density_linear will get modified in the process
            printf("Performing Quick and Dirty\n");
            quick_and_dirty(density_linear, density_quick);

            if (WRITE_DENSITY) {
                printf("Writing data to disc - density_quick\n");
                //write_data("density_quick", density_quick);
                strncpy(output, "density_quick_trial_", MAXCHAR);
                snprintf(integer, MAXCHAR, "%d", counter);
                strncat(output, integer, MAXCHAR);
                write_data(output, density_quick);
            }

            //record fidelity
            gsl_vector_set(f_quick, counter, fidelity(density_quick, density_physical));
            printf("Fidelity Quick: %e\n", gsl_vector_get(f_quick, counter));

            //restore density_linear
            gsl_matrix_complex_memcpy(density_linear, density_fp);

            //density_linear will get modified in the process
            printf("Performing Forced Purity\n");
            quick_and_dirty_2(density_linear, density_fp);

            if (WRITE_DENSITY) {
                printf("Writing data to disc - density_fp\n");
                strncpy(output, "density_fp_trial_", MAXCHAR);
                snprintf(integer, MAXCHAR, "%d", counter);
                strncat(output, integer, MAXCHAR);
                write_data(output, density_fp);
            }

            //record fidelity
            gsl_vector_set(f_fp, counter, fidelity(density_fp, density_physical));
            printf("Fidelity Forced Purity: %e\n", gsl_vector_get(f_fp, counter));

            printf("Performing MLE\n");
            MLE(n, density_physical, density_quick, density_mle, fidelities, times, f_vals, abs_gradient, 0);

            if (WRITE_DENSITY) {
                printf("Writing data to disc - density_mle\n");
                //last density_mle is written each time as "density_mle"
                write_data("density_mle", density_mle);
                strncpy(output, "density_mle_trial_", MAXCHAR);
                snprintf(integer, MAXCHAR, "%d", counter);
                strncat(output, integer, MAXCHAR);
                write_data(output, density_mle);
            }

            gsl_vector_set(f_mle, counter, fidelity(density_mle, density_physical));
            printf("Fidelity MLE: %e\n", gsl_vector_get(f_mle, counter));

            counter++;

        }
    }


    printf("Writing data to disc - fidelities and times vectors\n");
    write_vector_real("f_mle", f_mle);
    write_vector_real("f_quick", f_quick);
    write_vector_real("f_fp", f_fp);
    write_vector_real("tau", tau);
    write_vector_real("entropy_linear", sl);
    //write_matrix_real("fidelities", fidelities);
    //write_matrix_real("times", times);
    //write_matrix_real("f_vals", f_vals);
    //write_matrix_real("abs_gradient", abs_gradient);

    //deallocate memory here:
    gsl_matrix_complex_free(density_physical);
    gsl_matrix_complex_free(density_theoretical);
    gsl_matrix_complex_free(density_linear);
    gsl_matrix_complex_free(density_quick);
    gsl_matrix_complex_free(density_mle);
    gsl_matrix_complex_free(psi_outer);
    gsl_vector_complex_free(psi_ket);
    gsl_matrix_complex_free(rho_mems);
    gsl_vector_free(n);
    gsl_matrix_free(fidelities);
    gsl_matrix_free(times);

}

void
find_N() {

    //pseudo-random density matrix
    gsl_matrix_complex *density_physical = gsl_matrix_complex_alloc(pow(2,NQ), pow(2,NQ));
    //its copy which keeps the state - theoretical density_physical
    gsl_matrix_complex *density_theoretical = gsl_matrix_complex_alloc(pow(2,NQ), pow(2,NQ));
    //measurement vector
    gsl_vector *n = gsl_vector_alloc(pow(4, NQ));
    gsl_vector *n_copy = gsl_vector_alloc(pow(4, NQ));

    //declare and allocate all matrices here
    gsl_matrix_complex *density_linear = gsl_matrix_complex_alloc(pow(2, NQ), pow(2, NQ));
    gsl_matrix_complex *density_quick = gsl_matrix_complex_alloc(pow(2,NQ), pow(2,NQ));
    gsl_matrix_complex *density_fp = gsl_matrix_complex_alloc(pow(2,NQ), pow(2,NQ));


    //constants used to save density matrices at each iteration
    char output[MAXCHAR];
    char integer[MAXCHAR];

    //record data analysis
    gsl_vector *f_quick, *f_fp, *tau, *sl;

    gsl_vector *N_vec;
    //pre-allocate vectors
    f_quick = gsl_vector_alloc(NUMTRIALS + 1);
    f_fp = gsl_vector_alloc(NUMTRIALS + 1);
    tau = gsl_vector_alloc(NUMTRIALS + 1);
    sl = gsl_vector_alloc(NUMTRIALS + 1);
    N_vec = gsl_vector_alloc(NUMTRIALS + 1);

    gsl_vector_set_zero(f_quick);
    gsl_vector_set_zero(f_fp);
    gsl_vector_set_zero(tau);
    gsl_vector_set_zero(sl);
    gsl_vector_set_zero(N_vec);

    int counter = 0, j;

    for (j = 0; j<=NUMTRIALS; j++)
    {
        // set the tangle value
        theta_handle = (double)j/NUMTRIALS*(1.0/sqrt(2));
        prepare_state(density_theoretical);

        //record position in H map
        gsl_vector_set(tau, counter, pow(concurrence(density_theoretical), 2));
        gsl_vector_set(sl, counter, entropy_linear(density_theoretical));
        printf("Tau: %e\n", gsl_vector_get(tau, counter));
        printf("Entropy Linear: %e\n", gsl_vector_get(sl, counter));

        //restore theoretical state
        gsl_matrix_complex_memcpy(density_physical, density_theoretical);

        //simulate imperfect state - now we have rho_physical
        printf("Measuring the state\n");
        measure_state(density_physical, n);

        gsl_vector_memcpy(n_copy, n);


        //to save time, skip the first N values, because they produce low
        //fidelities anyway - these were determined by trial and error
        if (NQ == 2)
            N = 10;
        else if (NQ == 3)
            N = 100;
        else if (NQ == 4)
            N = 1000;
        else if (NQ == 5)
            N = 5000;
        else if (NQ == 6)
            N = 10000;
        else if (NQ == 7)
            N = 130000;
        else if (NQ == 8)
            N = 5000000;

        if (WRITE_DENSITY) {
            //save imperfect state - we will try to reconstruct this state
            printf("Writing data to disc - density_physical\n");
            //save the result
            strncpy(output, "density_physical_trial_", MAXCHAR);
            snprintf(integer, MAXCHAR, "%d", counter);
            strncat(output, integer, MAXCHAR);
            write_data(output, density_physical);
        }
        while(1) {

            N++;
            //restore measurements
            gsl_vector_memcpy(n, n_copy);

            //scale measurement vector here
            gsl_vector_scale(n, N);

            //apply measurement noise
            //printf("Applying experimental noise to data\n");
            apply_noise(n);

            //printf("Performing linear reconstruction\n");
            //perform linear reconstruction
            linear_estimate(n, density_linear);

            //save density_linear
            //gsl_matrix_complex_memcpy(density_fp, density_linear);

            //printf("Writing data to disc - density_linear\n");
            //strncpy(output, "density_linear_trial_", MAXCHAR);
            //snprintf(integer, MAXCHAR, "%d", counter);
            //strncat(output, integer, MAXCHAR);
            //write_data(output, density_linear);

            //density_linear will get modified in the process
            //printf("Performing Quick and Dirty\n");
            //quick_and_dirty(density_linear, density_quick);

            //record fidelity
            //gsl_vector_set(f_quick, counter, fidelity(density_quick, density_physical));
            //printf("Fidelity Quick: %e\n", gsl_vector_get(f_quick, counter));

            //restore density_linear
            //gsl_matrix_complex_memcpy(density_linear, density_fp);

            //density_linear will get modified in the process
            //printf("Performing Forced Purity\n");
            quick_and_dirty_2(density_linear, density_fp);

            printf("N value is: %e\n", N);

            //record fidelity
            gsl_vector_set(f_fp, counter, fidelity(density_fp, density_physical));
            printf("Fidelity Forced Purity: %e\n", gsl_vector_get(f_fp, counter));

            //termination condition
            if (gsl_vector_get(f_fp, counter) >= 0.90)
                break;

        }
        gsl_vector_set(N_vec, counter, N);
        printf("N final value is: %e\n", N);

        //     if (WRITE_DENSITY) {
        //       printf("Writing data to disc - density_linear\n");
        //       strncpy(output, "density_linear_trial_", MAXCHAR);
        //       snprintf(integer, MAXCHAR, "%d", counter);
        //       strncat(output, integer, MAXCHAR);
        //       write_data(output, density_linear);
        //     }

        //if (WRITE_DENSITY) {
        //  printf("Writing data to disc - density_quick\n");
        //  //write_data("density_quick", density_quick);
        //  strncpy(output, "density_quick_trial_", MAXCHAR);
        //  snprintf(integer, MAXCHAR, "%d", counter);
        //  strncat(output, integer, MAXCHAR);
        //  write_data(output, density_quick);
        //}

        if (WRITE_DENSITY) {
            printf("Writing data to disc - density_fp\n");
            strncpy(output, "density_fp_trial_", MAXCHAR);
            snprintf(integer, MAXCHAR, "%d", counter);
            strncat(output, integer, MAXCHAR);
            write_data(output, density_fp);
        }

        write_vector_real("N_vec", N_vec);
        counter++;
    }

    //display_vector_real(f_mle);

    printf("Writing data to disc - fidelities and times vectors\n");
    write_vector_real("f_quick", f_quick);
    write_vector_real("f_fp", f_fp);
    write_vector_real("tau", tau);
    write_vector_real("entropy_linear", sl);
    write_vector_real("N_vec", N_vec);

    //deallocate memory here:
    gsl_matrix_complex_free(density_physical);
    gsl_matrix_complex_free(density_theoretical);
    gsl_matrix_complex_free(density_linear);
    gsl_matrix_complex_free(density_quick);
    gsl_vector_free(n);

}

void
states_1_to_10() {

    int i, range;
    //pseudo-random density matrix
    gsl_matrix_complex *density_physical = gsl_matrix_complex_alloc(pow(2,NQ), pow(2,NQ));
    //its copy which keeps the state - theoretical density_physical
    gsl_matrix_complex *density_theoretical = gsl_matrix_complex_alloc(pow(2,NQ), pow(2,NQ));
    //measurement vector
    gsl_vector *n = gsl_vector_alloc(pow(4, NQ));

    //declare and allocate all matrices here
    gsl_matrix_complex *density_linear = gsl_matrix_complex_alloc(pow(2, NQ), pow(2, NQ));
    gsl_matrix_complex *density_quick = gsl_matrix_complex_alloc(pow(2,NQ), pow(2,NQ));
    gsl_matrix_complex *density_fp = gsl_matrix_complex_alloc(pow(2,NQ), pow(2,NQ));
    gsl_matrix_complex *density_mle = gsl_matrix_complex_alloc(pow(2, NQ), pow(2, NQ));
    gsl_vector *result = gsl_vector_alloc(pow(4,NQ));

    //constants used to save density matrices at each iteration
    char output[MAXCHAR];
    char integer[MAXCHAR];

    //start with some state, unless it's a random state which gets changed each time
    if (STATE!=9) {
        printf("Preparing the state\n");
        prepare_state(density_physical);
        write_data("density_state", density_physical);
        //save the original state with no noise
        gsl_matrix_complex_memcpy(density_theoretical, density_physical);
    }

    //define range of values for Werner_handle (or tangle handle - state 12)
    if (STATE != 10 && STATE!=12)
        range = 0;
    else
        range = NUMTRIALS;

    //record data analysis
    gsl_vector *f_mle, *f_quick, *f_fp, *tau, *sl, *times_qd, *times_fp, *times_lin, *times_meas;
    //vectors that store fidelities and times for each iteration
    gsl_matrix *fidelities, *times;
    //number of iterations before the fidelity starts to decrease
    gsl_matrix *f_vals, *abs_gradient;
    //pre-allocate vectors
    f_mle = gsl_vector_alloc(NUMTRIALS * (range+1));
    f_quick = gsl_vector_alloc(NUMTRIALS * (range+1));
    f_fp = gsl_vector_alloc(NUMTRIALS * (range+1));
    tau = gsl_vector_alloc(NUMTRIALS * (range+1));
    sl = gsl_vector_alloc(NUMTRIALS * (range+1));
    //allocate runtime matrices
    fidelities = gsl_matrix_alloc(NUMTRIALS * (range+1), MAXITER);
    times = gsl_matrix_alloc(NUMTRIALS * (range+1), MAXITER);
    times_lin = gsl_vector_alloc(NUMTRIALS * (range+1));
    times_qd = gsl_vector_alloc(NUMTRIALS * (range+1));
    times_fp = gsl_vector_alloc(NUMTRIALS * (range+1));
    times_meas = gsl_vector_alloc(NUMTRIALS * (range+1));
    f_vals = gsl_matrix_alloc(NUMTRIALS * (range+1), MAXITER);
    abs_gradient = gsl_matrix_alloc(NUMTRIALS * (range+1), MAXITER);
    gsl_vector_set_zero(f_mle);
    gsl_vector_set_zero(f_quick);
    gsl_vector_set_zero(f_fp);
    gsl_vector_set_zero(tau);
    gsl_vector_set_zero(sl);
    gsl_vector_set_all(times_lin, GSL_NAN);
    gsl_vector_set_all(times_qd, GSL_NAN);
    gsl_vector_set_all(times_fp, GSL_NAN);
    gsl_vector_set_all(times_meas, GSL_NAN);
    gsl_matrix_set_all(fidelities, GSL_NAN);
    gsl_matrix_set_all(times, GSL_NAN);
    gsl_matrix_set_all(f_vals, GSL_NAN);
    gsl_matrix_set_all(abs_gradient, GSL_NAN);

    int counter = 0, j;

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

    for (j = 0; j<=range; j++)
    {
        if (STATE == 10 || STATE == 12) {
            //printf("range %d\n", range);
            //for a Werner state, vary the handle (state 10)
            Werner_handle = (double)j/range;
            //for a random state, generate another random state
            //or set the tangle value
            theta_handle = Werner_handle*(1.0/sqrt(2));
            //printf("theta handle %e\n", theta_handle);
            //printf("range %d\n", range);
            prepare_state(density_theoretical);
        }

        for (i=0; i<NUMTRIALS; i++) {

            //display_matrix(density_theoretical);

            //update progress
            printf("===========================\n");
            printf("Iteration percent complete: %e\n", (double)counter/(range+1)/NUMTRIALS*100);
            printf("===========================\n");

            if (STATE == 9 ) {
                //for a random state, generate another random state
                prepare_state(density_theoretical);
            }

            //record position in H map
            gsl_vector_set(tau, counter, pow(concurrence(density_theoretical), 2));
            gsl_vector_set(sl, counter, entropy_linear(density_theoretical));
            printf("Tau: %e\n", gsl_vector_get(tau, counter));
            printf("Entropy Linear: %e\n", gsl_vector_get(sl, counter));

            //restore theoretical state
            gsl_matrix_complex_memcpy(density_physical, density_theoretical);

            //simulate imperfect state - now we have rho_physical
            printf("Measuring the state\n");
            gettimeofday(&tpstart, NULL);
            measure_state(density_physical, n);
            gettimeofday(&tpend, NULL);
            timedif = (double) MILLION*(tpend.tv_sec - tpstart.tv_sec);  /* + tpend.tv_usec - tpstart.tv_usec; */
            gsl_vector_set(times_meas, counter, timedif);

            //scale measurement vector here
            gsl_vector_scale(n, N);

            if (WRITE_DENSITY) {
                //save imperfect state - we will try to reconstruct this state
                printf("Writing data to disc - density_physical\n");
                //save the result
                strncpy(output, "density_physical_trial_", MAXCHAR);
                snprintf(integer, MAXCHAR, "%d", counter);
                strncat(output, integer, MAXCHAR);
                write_data(output, density_physical);
            }

            //apply measurement noise
            printf("Applying experimental noise to data\n");
            apply_noise(n);

            printf("Performing linear reconstruction\n");
            //perform linear reconstruction
            gettimeofday(&tpstart, NULL);
            linear_estimate(n, density_linear);
            gettimeofday(&tpend, NULL);
            timedif = (double) MILLION*(tpend.tv_sec - tpstart.tv_sec) + tpend.tv_usec - tpstart.tv_usec;
            gsl_vector_set(times_lin, counter, timedif);

            if (WRITE_DENSITY) {
                printf("Writing data to disc - density_linear\n");
                strncpy(output, "density_linear_trial_", MAXCHAR);
                snprintf(integer, MAXCHAR, "%d", counter);
                strncat(output, integer, MAXCHAR);
                write_data(output, density_linear);
            }

            //save density_linear
            gsl_matrix_complex_memcpy(density_fp, density_linear);

            //density_linear will get modified in the process
            printf("Performing Quick and Dirty\n");
            gettimeofday(&tpstart, NULL);
            quick_and_dirty(density_linear, density_quick);
            gettimeofday(&tpend, NULL);
            timedif = (double) MILLION*(tpend.tv_sec - tpstart.tv_sec) + tpend.tv_usec - tpstart.tv_usec;
            gsl_vector_set(times_qd, counter, timedif);

            if (WRITE_DENSITY) {
                printf("Writing data to disc - density_quick\n");
                //write_data("density_quick", density_quick);
                strncpy(output, "density_quick_trial_", MAXCHAR);
                snprintf(integer, MAXCHAR, "%d", counter);
                strncat(output, integer, MAXCHAR);
                write_data(output, density_quick);
            }

            //record fidelity
            gsl_vector_set(f_quick, counter, fidelity(density_quick, density_physical));
            printf("Fidelity Quick: %e\n", gsl_vector_get(f_quick, counter));

            //restore density_linear
            gsl_matrix_complex_memcpy(density_linear, density_fp);

            //density_linear will get modified in the process
            printf("Performing Forced Purity\n");
            gettimeofday(&tpstart, NULL);
            quick_and_dirty_2(density_linear, density_fp);
            gettimeofday(&tpend, NULL);
            timedif = (double) MILLION*(tpend.tv_sec - tpstart.tv_sec) + tpend.tv_usec - tpstart.tv_usec;
            gsl_vector_set(times_fp, counter, timedif);

            if (WRITE_DENSITY) {
                printf("Writing data to disc - density_fp\n");
                strncpy(output, "density_fp_trial_", MAXCHAR);
                snprintf(integer, MAXCHAR, "%d", counter);
                strncat(output, integer, MAXCHAR);
                write_data(output, density_fp);
            }

            //record fidelity
            gsl_vector_set(f_fp, counter, fidelity(density_fp, density_physical));
            printf("Fidelity Forced Purity: %e\n", gsl_vector_get(f_fp, counter));

            if (DO_MLE) {
                printf("Performing MLE\n");
                MLE(n, density_physical, density_quick, density_mle, fidelities, times, f_vals, abs_gradient, counter);

                if (WRITE_DENSITY) {
                    printf("Writing data to disc - density_mle\n");
                    //last density_mle is written each time as "density_mle"
                    //write_data("density_mle", density_mle);
                    strncpy(output, "density_mle_trial_", MAXCHAR);
                    snprintf(integer, MAXCHAR, "%d", counter);
                    strncat(output, integer, MAXCHAR);
                    write_data(output, density_mle);
                }

                gsl_vector_set(f_mle, counter, fidelity(density_mle, density_physical));
                printf("Fidelity MLE: %e\n", gsl_vector_get(f_mle, counter));
            }

            counter++;
        }
    }

    //display_vector_real(f_mle);

    printf("Writing data to disc - fidelities and times vectors\n");
    write_vector_real("f_mle", f_mle);
    write_vector_real("f_quick", f_quick);
    write_vector_real("f_fp", f_fp);
    write_vector_real("tau", tau);
    write_vector_real("times_lin", times_lin);
    write_vector_real("times_meas", times_meas);
    write_matrix_real("fidelities", fidelities);
    write_matrix_real("times_mle", times);
    write_vector_real("times_qd", times_qd);
    write_vector_real("times_fp", times_fp);
    write_matrix_real("f_vals", f_vals);
    write_matrix_real("abs_gradient", abs_gradient);

    //deallocate memory here:
    gsl_matrix_complex_free(density_physical);
    gsl_matrix_complex_free(density_theoretical);
    gsl_matrix_complex_free(density_linear);
    gsl_matrix_complex_free(density_quick);
    gsl_matrix_complex_free(density_mle);
    gsl_vector_free(n);
    gsl_vector_free(result);
    //already deallocated
    gsl_matrix_free(fidelities);
    gsl_matrix_free(times);
}

void
quick_and_dirty(gsl_matrix_complex *density_linear, gsl_matrix_complex *density_quick) {

    int i;
    //allocate workspace
    gsl_eigen_hermv_workspace *w = gsl_eigen_hermv_alloc (pow(2,NQ));
    //allocate eigenvalue vector
    gsl_vector *eval = gsl_vector_alloc(pow(2,NQ));
    //allocate eigenvector matrix - mutually orthogonal, normalized to 1
    gsl_matrix_complex *evec = gsl_matrix_complex_alloc(pow(2,NQ), pow(2,NQ));

    //leave density_linear intact
    gsl_matrix_complex_memcpy(density_quick, density_linear);
    if ((gsl_eigen_hermv (density_quick, eval, evec, w)) != 0) {
        //could not find eigenvalues
        gsl_matrix_complex_memcpy(density_quick, density_linear);
        printf("Failed at finding eigenvalues and eigenvectors\n");
        printf("Quick and Dirty Estimate not available\n");
        return;
    } else {
        if (VERBOSE)
            printf("Succeeded in finding eigenvalues and eigenvectors\n");
    }

    //display_matrix(evec);
    //   printf("Original eigenvalues\n");
    //   for (i=0; i<pow(2,NQ); i++)
    //     printf("%e ", gsl_vector_get(eval, i));
    //   printf("\n");

    //replace negative eigenvalues by zero
    for (i=0; i<pow(2,NQ); i++)
        if (gsl_vector_get(eval, i)<0)
            gsl_vector_set(eval, i, 0);

    gsl_eigen_hermv_sort(eval, evec, GSL_EIGEN_SORT_ABS_ASC);

    //matrix with eigenvalues across the diagonal
    gsl_matrix_complex *eval_matrix = gsl_matrix_complex_alloc(pow(2,NQ), pow(2, NQ));
    gsl_matrix_complex_set_zero(eval_matrix);
    for (i=0; i<pow(2,NQ); i++)
        gsl_matrix_complex_set(eval_matrix, i, i, gsl_complex_rect(gsl_vector_get(eval, i), 0));


    //perform V*D*inv(V)
    gsl_matrix_complex *tmp = gsl_matrix_complex_alloc(pow(2,NQ), pow(2,NQ));
    gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, one, evec, eval_matrix, zero, tmp);
    gsl_blas_zgemm(CblasNoTrans, CblasConjTrans, one, tmp, evec, zero, density_quick);

    //re-normalize the trace
    gsl_matrix_complex_scale(density_quick, gsl_complex_rect(1.0/GSL_REAL(trace(density_quick)), 0));

    //deallocate memory
    gsl_matrix_complex_free(eval_matrix);
    gsl_eigen_hermv_free(w);
    gsl_vector_free(eval);
    gsl_matrix_complex_free(evec);
    gsl_matrix_complex_free(tmp);
}

void
quick_and_dirty_2(gsl_matrix_complex *density_linear, gsl_matrix_complex *density_quick) {

    int i;
    //allocate workspace
    gsl_eigen_hermv_workspace *w = gsl_eigen_hermv_alloc (pow(2,NQ));
    //allocate eigenvalue vector
    gsl_vector *eval = gsl_vector_alloc(pow(2,NQ));
    //allocate eigenvector matrix - mutually orthogonal, normalized to 1
    gsl_matrix_complex *evec = gsl_matrix_complex_alloc(pow(2,NQ), pow(2,NQ));

    //leave density_linear intact
    gsl_matrix_complex_memcpy(density_quick, density_linear);
    if ((gsl_eigen_hermv (density_quick, eval, evec, w)) != 0) {
        //could not find eigenvalues
        gsl_matrix_complex_memcpy(density_quick, density_linear);
        printf("Failed at finding eigenvalues and eigenvectors\n");
        printf("Quick and Dirty Estimate not available\n");
        return;
    } else {
        if (VERBOSE)
            printf("Succeeded in finding eigenvalues and eigenvectors\n");
    }

    //display_matrix(evec);
    gsl_eigen_hermv_sort(eval, evec, GSL_EIGEN_SORT_VAL_ASC);

    //   printf("Original eigenvalues\n");
    //   for (i=0; i<pow(2,NQ); i++)
    //     printf("%e ", gsl_vector_get(eval, i));
    //   printf("\n");

    for (i=0; i<pow(2,NQ); i++)
        if (i==pow(2,NQ)-1)
            gsl_vector_set(eval, i, 1.0);
        else
            gsl_vector_set(eval, i, 0.0);

    //   for (i=0; i<pow(2,NQ); i++)
    //     printf("%e ", gsl_vector_get(eval, i));
    //   printf("\n");

    //matrix with eigenvalues across the diagonal
    gsl_matrix_complex *eval_matrix = gsl_matrix_complex_alloc(pow(2,NQ), pow(2, NQ));
    gsl_matrix_complex_set_zero(eval_matrix);
    for (i=0; i<pow(2,NQ); i++)
        gsl_matrix_complex_set(eval_matrix, i, i, gsl_complex_rect(gsl_vector_get(eval, i), 0));


    //perform V*D*inv(V)
    gsl_matrix_complex *tmp = gsl_matrix_complex_alloc(pow(2,NQ), pow(2,NQ));
    gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, one, evec, eval_matrix, zero, tmp);
    gsl_blas_zgemm(CblasNoTrans, CblasConjTrans, one, tmp, evec, zero, density_quick);

    //re-normalize the trace
    gsl_matrix_complex_scale(density_quick, gsl_complex_rect(1.0/GSL_REAL(trace(density_quick)), 0));

    //deallocate memory
    gsl_matrix_complex_free(eval_matrix);
    gsl_eigen_hermv_free(w);
    gsl_vector_free(eval);
    gsl_matrix_complex_free(evec);
    gsl_matrix_complex_free(tmp);
}

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) {

    //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);

}


//here n are the parameters
double
my_f(const gsl_vector * x, void * n) {
    gsl_matrix_complex *rho_mat = gsl_matrix_complex_alloc(pow(2, NQ), pow(2, NQ));
    gsl_matrix_complex *T_mat = gsl_matrix_complex_alloc(pow(2, NQ), pow(2, NQ));
    make_T(x, T_mat);

    gsl_blas_zgemm(CblasConjTrans, CblasNoTrans, one, T_mat, T_mat, zero, rho_mat);
    //printf("TRACE: %e\n", GSL_REAL(trace(rho_mat)));
    gsl_matrix_complex_scale(rho_mat, gsl_complex_rect(1/GSL_REAL(trace(rho_mat)), 0));

    int i, nu;
    double ret=0;

    //init the counter
    int counter[NQ], old_counter[NQ];
    for (i=0; i<NQ-1; i++) {
        counter[i]=0;
        old_counter[i]=0;
    }
    counter[NQ-1]=-1;
    old_counter[NQ-1]=-1;

    //init mu_cell
    gsl_matrix_complex_memcpy(mu_cell[0], mu_array[0]);
    for (i=1; i<NQ; i++)
        kron(mu_cell[i-1], mu_array[0], mu_cell[i]);

    for (nu=0; nu<pow(4, NQ); nu++) {
        base_4_rep_inc(counter, nu);
        //generate new measurement operator
        mu_tensor(mu_cell, mu_array, old_counter, counter);
        //update counter
        for (i=0; i<NQ; i++)
            old_counter[i] = counter[i];

        //don't need T_mat anymore, use it as a temp matrix
        gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, one, mu_cell[NQ-1], rho_mat, zero, T_mat);
        //T_mat = mu_nu*rho
        //need: alpha_return = alpha_return + ((N*trace(mu_nu*rho)-n(nu+1))^2)/n(nu+1);

        ret+=pow(((double)N*GSL_REAL(trace(T_mat))-gsl_vector_get(n, nu)),2)/gsl_vector_get(n,nu);
    }
    ret*=(double)1/2;

    gsl_matrix_complex_free(T_mat);
    gsl_matrix_complex_free(rho_mat);

    return ret;
}

void //params == n
my_df(const gsl_vector * x_vec, void * n, gsl_vector * g) {
    gsl_matrix_complex *rho_mat = gsl_matrix_complex_alloc(pow(2, NQ), pow(2, NQ));
    gsl_matrix_complex *T_mat = gsl_matrix_complex_alloc(pow(2, NQ), pow(2, NQ));
    make_T(x_vec, T_mat);
    gsl_matrix_complex *T_dag_T = gsl_matrix_complex_alloc(pow(2, NQ), pow(2, NQ));
    gsl_blas_zgemm(CblasConjTrans, CblasNoTrans, one, T_mat, T_mat, zero, rho_mat);
    gsl_matrix_complex_memcpy(T_dag_T, rho_mat);
    //printf("TRACE: %e\n", GSL_REAL(trace(rho_mat)));
    double tr_T_dag_T = GSL_REAL(trace(T_dag_T));
    gsl_matrix_complex_scale(rho_mat, gsl_complex_rect(1/tr_T_dag_T,0));
    //display_matrix(T_dag_T);
    //display_matrix(rho_mat);

    int i, j;
    gsl_matrix *X, *Y;
    X = gsl_matrix_alloc(pow(2, NQ), pow(2, NQ));
    Y = gsl_matrix_alloc(pow(2, NQ), pow(2, NQ));
    for (i=0; i<pow(2, NQ); i++) {
        for (j=0; j<pow(2, NQ); j++) {
            gsl_matrix_set(X, i, j, GSL_REAL(gsl_matrix_complex_get(T_mat, i, j)));
            gsl_matrix_set(Y, i, j, GSL_IMAG(gsl_matrix_complex_get(T_mat, i, j)));
        }
    }

    gsl_matrix_complex *Answer = gsl_matrix_complex_alloc(pow(2, NQ), pow(2, NQ));

    gsl_matrix_complex_set_zero(Answer);

    int counter[NQ], old_counter[NQ];
    for (i=0; i<NQ-1; i++) {
        counter[i]=0;
        old_counter[i]=0;
    }
    counter[NQ-1]=-1;
    old_counter[NQ-1]=-1;

    //init mu_cell
    gsl_matrix_complex_memcpy(mu_cell[0], mu_array[0]);
    for (i=1; i<NQ; i++)
        kron(mu_cell[i-1], mu_array[0], mu_cell[i]);

    //temporary matrix used for multiplication
    gsl_matrix_complex *temp = gsl_matrix_complex_alloc(pow(2,NQ), pow(2, NQ));
    gsl_matrix_complex_set_zero(temp); //init just in case

    gsl_matrix *K, *L;
    K = gsl_matrix_alloc(pow(2, NQ), pow(2, NQ));
    L = gsl_matrix_alloc(pow(2, NQ), pow(2, NQ));

    gsl_matrix *XK, *YL, *XL, *YK;
    XK = gsl_matrix_alloc(pow(2, NQ), pow(2, NQ));
    YK = gsl_matrix_alloc(pow(2, NQ), pow(2, NQ));
    XL = gsl_matrix_alloc(pow(2, NQ), pow(2, NQ));
    YL = gsl_matrix_alloc(pow(2, NQ), pow(2, NQ));
    //gsl_matrix_complex *D_mat = gsl_matrix_complex_alloc(pow(2, NQ), pow(2, NQ));

    gsl_complex temp_num;
    double real, imag;
    double x, y, xk, xl, yl, yk;
    int nu;
    double tr_mu_rho, tr_mu_T_dag_T, scale;
    for (nu=0; nu<pow(4, NQ); nu++) {
        base_4_rep_inc(counter, nu);
        //generate new measurement operator
        mu_tensor(mu_cell, mu_array, old_counter, counter);
        for (i=0; i<NQ; i++)
            old_counter[i] = counter[i];

        gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, one, mu_cell[NQ-1], rho_mat, zero, temp);
        tr_mu_rho = GSL_REAL(trace(temp));
        gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, one, mu_cell[NQ-1], T_dag_T, zero, temp);
        tr_mu_T_dag_T = GSL_REAL(trace(temp));

        scale = ((double)N*tr_mu_rho - gsl_vector_get(n,nu))/gsl_vector_get(n,nu);

        for (i=0; i<pow(2, NQ); i++) {
            for (j=0; j<pow(2, NQ); j++) {
                gsl_matrix_set(K, i, j, GSL_REAL(gsl_matrix_complex_get(mu_cell[NQ-1], i, j)));
                gsl_matrix_set(L, i, j, GSL_IMAG(gsl_matrix_complex_get(mu_cell[NQ-1], i, j)));
            }
        }

        gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, X, K, 0.0, XK);
        gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, X, L, 0.0, XL);
        gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, Y, K, 0.0, YK);
        gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, Y, L, 0.0, YL);

        //fill in Answer:
        //only care about the lower diagonal, because that's where the answer is stored

        for (i=0; i<pow(2, NQ); i++) {
            for (j=0; j<i+1; j++) {
                x = gsl_matrix_get(X, i, j);
                y = gsl_matrix_get(Y, i, j);
                xk = gsl_matrix_get(XK, i, j);
                xl = gsl_matrix_get(XL, i, j);
                yk = gsl_matrix_get(YK, i, j);
                yl = gsl_matrix_get(YL, i, j);
                real = tr_T_dag_T*2*(xk - yl) - tr_mu_T_dag_T*2*x;
                imag = tr_T_dag_T*2*(xl + yk) - tr_mu_T_dag_T*2*y;
                //printf("%e\n", real);
                temp_num = gsl_complex_rect(real, imag);
                //temp_num *= scale*N;
                temp_num = gsl_complex_mul(temp_num, gsl_complex_rect(scale,0));
                gsl_matrix_complex_set(Answer, i, j,
                                       gsl_complex_add(gsl_matrix_complex_get(Answer,i,j),temp_num));

            }
        }

    }

    //printf("%e\n", 1/(tr_T_dag_T*tr_T_dag_T));

    gsl_matrix_complex_scale(Answer, gsl_complex_rect(N/(tr_T_dag_T*tr_T_dag_T), 0));
    //display_matrix(Answer);
    //now reseed gradient values from matrix into the gradient vector
    int row = 2;
    int col = 1;
    int prev_row = row;

    i=1;
    while (i<=pow(4,NQ)) {
        if (i<=pow(2, NQ)) {
            gsl_vector_set(g, i-1, GSL_REAL(gsl_matrix_complex_get(Answer, i-1,i-1)));
        } else {
            gsl_vector_set(g, i-1, GSL_REAL(gsl_matrix_complex_get(Answer, row-1,col-1)));
            i++;
            gsl_vector_set(g, i-1, GSL_IMAG(gsl_matrix_complex_get(Answer, row-1, col-1)));
            row++;
            col++;
            if (row>pow(2, NQ)) {
                prev_row++;
                row = prev_row;
                col=1;
            }
        }
        i++;
    }

    //deallocate memory
    gsl_matrix_complex_free(Answer);
    gsl_matrix_free(X);
    gsl_matrix_free(Y);
    gsl_matrix_free(K);
    gsl_matrix_free(L);
    //////////////////
    gsl_matrix_free(XK);
    gsl_matrix_free(YL);
    gsl_matrix_free(YK);
    gsl_matrix_free(XL);
    //gsl_matrix_free(D_mat);
    gsl_matrix_complex_free(T_mat);
    gsl_matrix_complex_free(rho_mat);
    gsl_matrix_complex_free(temp);
    gsl_matrix_complex_free(T_dag_T);


    //fin_diff_grad(x_vec, n, g);
}

void //params = n
my_fdf(const gsl_vector * x_vec, void *n, double * f, gsl_vector * g) {

    gsl_matrix_complex *rho_mat = gsl_matrix_complex_alloc(pow(2, NQ), pow(2, NQ));
    gsl_matrix_complex *T_mat = gsl_matrix_complex_alloc(pow(2, NQ), pow(2, NQ));
    make_T(x_vec, T_mat);
    gsl_matrix_complex *T_dag_T = gsl_matrix_complex_alloc(pow(2, NQ), pow(2, NQ));
    gsl_blas_zgemm(CblasConjTrans, CblasNoTrans, one, T_mat, T_mat, zero, rho_mat);
    gsl_matrix_complex_memcpy(T_dag_T, rho_mat);
    //printf("TRACE: %e\n", GSL_REAL(trace(rho_mat)));
    double tr_T_dag_T = GSL_REAL(trace(T_dag_T));
    gsl_matrix_complex_scale(rho_mat, gsl_complex_rect(1/tr_T_dag_T,0));
    //display_matrix(T_dag_T);
    //display_matrix(rho_mat);

    int i, j;
    gsl_matrix *X, *Y;
    X = gsl_matrix_alloc(pow(2, NQ), pow(2, NQ));
    Y = gsl_matrix_alloc(pow(2, NQ), pow(2, NQ));
    for (i=0; i<pow(2, NQ); i++) {
        for (j=0; j<pow(2, NQ); j++) {
            gsl_matrix_set(X, i, j, GSL_REAL(gsl_matrix_complex_get(T_mat, i, j)));
            gsl_matrix_set(Y, i, j, GSL_IMAG(gsl_matrix_complex_get(T_mat, i, j)));
        }
    }

    gsl_matrix_complex *Answer = gsl_matrix_complex_alloc(pow(2, NQ), pow(2, NQ));

    gsl_matrix_complex_set_zero(Answer);

    int counter[NQ], old_counter[NQ];
    for (i=0; i<NQ-1; i++) {
        counter[i]=0;
        old_counter[i]=0;
    }
    counter[NQ-1]=-1;
    old_counter[NQ-1]=-1;

    //init mu_cell
    gsl_matrix_complex_memcpy(mu_cell[0], mu_array[0]);
    for (i=1; i<NQ; i++)
        kron(mu_cell[i-1], mu_array[0], mu_cell[i]);

    //temporary matrix used for multiplication
    gsl_matrix_complex *temp = gsl_matrix_complex_alloc(pow(2,NQ), pow(2, NQ));
    gsl_matrix_complex_set_zero(temp); //init just in case

    gsl_matrix *K, *L;
    K = gsl_matrix_alloc(pow(2, NQ), pow(2, NQ));
    L = gsl_matrix_alloc(pow(2, NQ), pow(2, NQ));

    gsl_matrix *XK, *YL, *XL, *YK;
    XK = gsl_matrix_alloc(pow(2, NQ), pow(2, NQ));
    YK = gsl_matrix_alloc(pow(2, NQ), pow(2, NQ));
    XL = gsl_matrix_alloc(pow(2, NQ), pow(2, NQ));
    YL = gsl_matrix_alloc(pow(2, NQ), pow(2, NQ));
    //gsl_matrix_complex *D_mat = gsl_matrix_complex_alloc(pow(2, NQ), pow(2, NQ));

    gsl_complex temp_num;
    double real, imag;
    double x, y, xk, xl, yl, yk;
    int nu;
    double tr_mu_rho, tr_mu_T_dag_T, scale;
    (*f) = 0;
    for (nu=0; nu<pow(4, NQ); nu++) {
        base_4_rep_inc(counter, nu);
        //generate new measurement operator
        mu_tensor(mu_cell, mu_array, old_counter, counter);
        for (i=0; i<NQ; i++)
            old_counter[i] = counter[i];

        gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, one, mu_cell[NQ-1], rho_mat, zero, temp);
        tr_mu_rho = GSL_REAL(trace(temp));
        gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, one, mu_cell[NQ-1], T_dag_T, zero, temp);
        tr_mu_T_dag_T = GSL_REAL(trace(temp));

        scale = ((double)N*tr_mu_rho - gsl_vector_get(n,nu))/gsl_vector_get(n,nu);

        (*f) +=  pow((double)N*tr_mu_rho - gsl_vector_get(n,nu), 2)/gsl_vector_get(n,nu);

        for (i=0; i<pow(2, NQ); i++) {
            for (j=0; j<pow(2, NQ); j++) {
                gsl_matrix_set(K, i, j, GSL_REAL(gsl_matrix_complex_get(mu_cell[NQ-1], i, j)));
                gsl_matrix_set(L, i, j, GSL_IMAG(gsl_matrix_complex_get(mu_cell[NQ-1], i, j)));
            }
        }

        gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, X, K, 0.0, XK);
        gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, X, L, 0.0, XL);
        gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, Y, K, 0.0, YK);
        gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, Y, L, 0.0, YL);

        //fill in Answer:
        //only care about the lower diagonal, because that's where the answer is stored

        for (i=0; i<pow(2, NQ); i++) {
            for (j=0; j<i+1; j++) {
                x = gsl_matrix_get(X, i, j);
                y = gsl_matrix_get(Y, i, j);
                xk = gsl_matrix_get(XK, i, j);
                xl = gsl_matrix_get(XL, i, j);
                yk = gsl_matrix_get(YK, i, j);
                yl = gsl_matrix_get(YL, i, j);
                real = tr_T_dag_T*2*(xk - yl) - tr_mu_T_dag_T*2*x;
                imag = tr_T_dag_T*2*(xl + yk) - tr_mu_T_dag_T*2*y;
                //printf("%e\n", real);
                temp_num = gsl_complex_rect(real, imag);
                //temp_num *= scale*N;
                temp_num = gsl_complex_mul(temp_num, gsl_complex_rect(scale,0));
                gsl_matrix_complex_set(Answer, i, j,
                                       gsl_complex_add(gsl_matrix_complex_get(Answer,i,j),temp_num));

            }
        }

    }

    //printf("%e\n", 1/(tr_T_dag_T*tr_T_dag_T));

    (*f) *= (double) 1/2; //scale objective function

    gsl_matrix_complex_scale(Answer, gsl_complex_rect(N/(tr_T_dag_T*tr_T_dag_T), 0));
    //display_matrix(Answer);
    //now reseed gradient values from matrix into the gradient vector
    int row = 2;
    int col = 1;
    int prev_row = row;

    i=1;
    while (i<=pow(4,NQ)) {
        if (i<=pow(2, NQ)) {
            gsl_vector_set(g, i-1, GSL_REAL(gsl_matrix_complex_get(Answer, i-1,i-1)));
        } else {
            gsl_vector_set(g, i-1, GSL_REAL(gsl_matrix_complex_get(Answer, row-1,col-1)));
            i++;
            gsl_vector_set(g, i-1, GSL_IMAG(gsl_matrix_complex_get(Answer, row-1, col-1)));
            row++;
            col++;
            if (row>pow(2, NQ)) {
                prev_row++;
                row = prev_row;
                col=1;
            }
        }
        i++;
    }

    //deallocate memory
    gsl_matrix_complex_free(Answer);
    gsl_matrix_free(X);
    gsl_matrix_free(Y);
    gsl_matrix_free(K);
    gsl_matrix_free(L);
    //////////////////
    gsl_matrix_free(XK);
    gsl_matrix_free(YL);
    gsl_matrix_free(YK);
    gsl_matrix_free(XL);
    //gsl_matrix_free(D_mat);
    gsl_matrix_complex_free(T_mat);
    gsl_matrix_complex_free(rho_mat);
    gsl_matrix_complex_free(temp);
    gsl_matrix_complex_free(T_dag_T);


}

// void
// fin_diff_grad(const gsl_vector *x, gsl_vector *n, gsl_vector *g) {
//   int i;
//   double h = 0.0001;
//   gsl_vector *xl = gsl_vector_alloc(pow(4, NQ));
//   gsl_vector *xr = gsl_vector_alloc(pow(4, NQ));
//   for (i=0; i<pow(4, NQ); i++) {
//     gsl_vector_memcpy(xl, x);
//     gsl_vector_memcpy(xr, x);
//     gsl_vector_set(xl, i, gsl_vector_get(xl, i) - h);
//     gsl_vector_set(xr, i, gsl_vector_get(xr, i) + h);
//     gsl_vector_set(g, i, (my_f(xr, n) - my_f(xl, n))/(2*h));
//   }
//
//   gsl_vector_free(xl);
//   gsl_vector_free(xr);
//
// }

void
linear_estimate(gsl_vector *n, gsl_matrix_complex *density_linear) {
    //init inverse matrix
    gsl_matrix *gamma = gsl_matrix_alloc(4, 4);
    gsl_matrix_set(gamma, 0, 0, 1);
    gsl_matrix_set(gamma, 0, 1, 0);
    gsl_matrix_set(gamma, 0, 2, 0);
    gsl_matrix_set(gamma, 0, 3, 0);
    gsl_matrix_set(gamma, 1, 0, 1);
    gsl_matrix_set(gamma, 1, 1, 0);
    gsl_matrix_set(gamma, 1, 2, -1);
    gsl_matrix_set(gamma, 1, 3, 0);
    gsl_matrix_set(gamma, 2, 0, 1);
    gsl_matrix_set(gamma, 2, 1, 0);
    gsl_matrix_set(gamma, 2, 2, 0);
    gsl_matrix_set(gamma, 2, 3, -1);
    gsl_matrix_set(gamma, 3, 0, -1);
    gsl_matrix_set(gamma, 3, 1, 1);
    gsl_matrix_set(gamma, 3, 2, 0);
    gsl_matrix_set(gamma, 3, 3, 0);

    //display_matrix_real(gamma);

    gsl_matrix_complex_set_zero(density_linear);
    int base_4_nu[NQ], base_4_mu[NQ], base_4_nu_old[NQ], i;
    int nu, mu, bcount;
    double temp, B_inv_element;
    gsl_matrix_complex *temporary = gsl_matrix_complex_alloc(pow(2, NQ), pow(2, NQ));
    for (i=0; i<NQ-1; i++) {
        base_4_nu[i]=0;
        base_4_nu_old[i]=0;
        base_4_mu[i]=0;
    }
    base_4_nu[NQ-1] = -1;
    base_4_mu[NQ-1] = -1;
    base_4_nu_old[NQ-1] = -1;
//////////////////////////////////check code BELOW
    for (nu=0; nu<pow(4,NQ); nu++) {
        temp = 0;
        base_4_rep_inc(base_4_nu, nu);

        for (mu=0; mu<pow(4,NQ); mu++) {
            base_4_rep_inc(base_4_mu, mu);
            //       for (i=0; i<NQ; i++)
            //         printf("%d",base_4_mu[i]);
            //       printf("\n");
            B_inv_element = 1;
            for (bcount = 0; bcount<NQ; bcount++) {
                B_inv_element*=gsl_matrix_get(gamma, base_4_nu[bcount], base_4_mu[bcount]);
            }
            temp+=B_inv_element*gsl_vector_get(n, mu);
        }
        //printf("temp: %e\n", temp);
        //reset base_4_mu
        for (i=0; i<NQ-1; i++)
            base_4_mu[i] = 0;
        base_4_mu[NQ-1] = -1;
        //generate Gamma matrix
        sig_tensor(sigma_cell, sig_array, base_4_nu_old, base_4_nu);
        //printf("%d\n", nu);
        //display_matrix(sigma_cell[NQ-1]);
        for (i=0; i<NQ; i++)
            base_4_nu_old[i] = base_4_nu[i];
        //sigma_cell[NQ-1] now contains the matrix that we need
        gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, one, sigma_cell[NQ-1], sigma_cell[NQ-1], zero, temporary);

        //printf("%d\n", nu);
        //display_matrix(temporary);

        gsl_complex normalize;
        GSL_SET_COMPLEX(&normalize, sqrt(1/GSL_REAL(trace(temporary))), 0);
        gsl_matrix_complex_scale(sigma_cell[NQ-1], normalize);
        gsl_matrix_complex_scale(sigma_cell[NQ-1], gsl_complex_rect(temp, 0));
        //display_matrix(sigma_cell[NQ-1]);
        gsl_matrix_complex_add(density_linear, sigma_cell[NQ-1]);

    }

    //renormalize the matrix:
    gsl_complex tmp;
    GSL_SET_COMPLEX(&tmp, 1/gsl_vector_get(n, 0)/sqrt(pow(2,NQ)), 0);
    gsl_matrix_complex_scale(density_linear, tmp);
    //deallocate memory
    gsl_matrix_free(gamma);
    gsl_matrix_complex_free(temporary);

}

void
init() {
    //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
prepare_state(gsl_matrix_complex *density_random) {


    //init random matrix - used to apply experimental noise mask
    // The folloding creates a random density matrix "random"
    //gsl_matrix_complex *random = gsl_matrix_complex_alloc(pow(2,NQ), pow(2,NQ));
    //matrix_complex_random_init(random);
    //compute density_random = ctranspose(random)*random
    //as density_random = 1*ctranspose(random)*random + 0*density_random
    //gsl_blas_zgemm(CblasConjTrans, CblasNoTrans, one, random, random, zero, density_random);
    //gsl_matrix_complex_memcpy(random, density_random);
    //gsl_matrix_complex_scale(random, gsl_complex_rect(1/GSL_REAL(trace(density_random)), 0));

    //constants used to apply experimental purity noise
    //gsl_complex epsilon;
    //GSL_SET_COMPLEX(&epsilon, 0.1, 0);
    //gsl_complex one_min_epsilon;
    //GSL_SET_COMPLEX(&one_min_epsilon, 0.9, 0);

    if (STATE==1) { //GHZ state
        //density_random = (1-epsilon)*GHZ*ctranspose(GHZ) + epsilon*density_random;

        gsl_vector_complex *ghz_vec = gsl_vector_complex_alloc(pow(2, NQ));
        gsl_vector_complex_set_zero(ghz_vec);
        gsl_vector_complex_set(ghz_vec, 0, gsl_complex_rect(1/sqrt(2), 0));
        gsl_vector_complex_set(ghz_vec, pow(2, NQ)-1, gsl_complex_rect(1/sqrt(2), 0));
        //gsl_complex neg_eps = gsl_complex_negative(epsilon);

        //gsl_matrix_complex *outp = gsl_matrix_complex_alloc(pow(2, NQ), pow(2, NQ));
        outer_product(ghz_vec, density_random);

        //free memory
        gsl_vector_complex_free(ghz_vec);
    } else if (STATE==2) { //unentangled, pure state (unmixed)

        gsl_matrix_complex_set_zero(density_random);
        gsl_matrix_complex_set(density_random, 0, 0, gsl_complex_rect(1, 0));

    } else if (STATE==3 || STATE==4 || STATE ==5 || STATE == 10) { //Werner states
        gsl_matrix_complex *identity = gsl_matrix_complex_alloc(pow(2, NQ), pow(2, NQ));
        gsl_matrix_complex_set_zero(identity);
        int i;
        for (i=0; i<pow(2, NQ); i++)
            gsl_matrix_complex_set(identity, i, i, gsl_complex_rect(1/pow(2, NQ), 0));

        gsl_complex eps, one_min_eps;
        double value;
        if (STATE==3)
            value = 0.2;
        else if (STATE==4)
            value = 0.80475;
        else if (STATE==5)
            value = 0.99;
        else if (STATE == 10)
            value = Werner_handle;

        GSL_SET_COMPLEX(&eps, value, 0);
        GSL_SET_COMPLEX(&one_min_eps, 1.0-value, 0);

        //create GHZ matrix
        gsl_vector_complex *ghz_vec = gsl_vector_complex_alloc(pow(2, NQ));
        gsl_vector_complex_set_zero(ghz_vec);
        gsl_vector_complex_set(ghz_vec, 0, gsl_complex_rect(1/sqrt(2), 0));
        gsl_vector_complex_set(ghz_vec, pow(2, NQ)-1, gsl_complex_rect(1/sqrt(2), 0));
        outer_product(ghz_vec, density_random);
        gsl_vector_complex_free(ghz_vec);

        gsl_matrix_complex_scale(identity, one_min_eps);
        gsl_matrix_complex_scale(density_random, eps);

        gsl_matrix_complex_add(density_random, identity);

        gsl_matrix_complex_free(identity);
    } else if (STATE == 6) { //MEMS State - only available for 2 QB
        if (NQ != 2) {
            fprintf(stderr, "MEMS state is only defined for 2 Qubits\n");
            exit(1);
        }
        gsl_matrix_complex_set_zero(density_random);
        //     gsl_matrix_complex_set(density_random, 0, 0, gsl_complex_rect(1.0/3.0, 0));
        //     gsl_matrix_complex_set(density_random, 1, 1, gsl_complex_rect(1.0/3.0, 0));
        //     gsl_matrix_complex_set(density_random, 3, 3, gsl_complex_rect(1.0/3.0, 0));

        double gam = 0.74 + 0.0123;

        //     gsl_matrix_complex_set(density_random, 0, 3, gsl_complex_rect(gam/2.0, 0));
        //     gsl_matrix_complex_set(density_random, 3, 0, gsl_complex_rect(gam/2.0, 0));

        gsl_matrix_complex_set(density_random, 0, 0, gsl_complex_rect(gam/2.0, 0));
        gsl_matrix_complex_set(density_random, 0, 3, gsl_complex_rect(gam/2.0, 0));
        gsl_matrix_complex_set(density_random, 1, 1, gsl_complex_rect(1.0-gam, 0));
        gsl_matrix_complex_set(density_random, 3, 0, gsl_complex_rect(gam/2.0, 0));
        gsl_matrix_complex_set(density_random, 3, 3, gsl_complex_rect(gam/2.0, 0));

    } else if (STATE == 7) { //tau = 0; S_l = 1
        gsl_matrix_complex_set_zero(density_random);
        int i;
        for (i=0; i<pow(2, NQ); i++)
            gsl_matrix_complex_set(density_random, i, i, gsl_complex_rect(1/pow(2, NQ), 0));

    } else if (STATE == 8) { // W state |W> = 1/sqrt(NQ) * ( |00...01> + |00...10> + ... + |10...0> )

        gsl_vector_complex *w_state = gsl_vector_complex_alloc(pow(2, NQ));
        gsl_matrix_complex *middle = gsl_matrix_complex_alloc(1, 2);
        gsl_vector_complex_set_zero(w_state);
        // |0> = (1 0)'   |1> = (0 1)'
        int temp[NQ], i, j;
        for (i=0; i<NQ; i++)
            temp[i] = 0;


        //fake allocation
        gsl_matrix_complex *after = gsl_matrix_complex_alloc(1, 1);
        gsl_matrix_complex *before = gsl_matrix_complex_alloc(1, 1);
        for (i=0; i<NQ; i++) {
            temp[i] = 1;

            gsl_matrix_complex_free(before);
            gsl_matrix_complex *before = gsl_matrix_complex_alloc(1, 2);
            // init before
            if (temp[0] == 1) { //|1> = (0 1)'
                gsl_matrix_complex_set(before, 0, 0, gsl_complex_rect(0.0, 0.0));
                gsl_matrix_complex_set(before, 0, 1, gsl_complex_rect(1.0, 0.0));
            } else {  //|0> = (1 0)'
                gsl_matrix_complex_set(before, 0, 0, gsl_complex_rect(1.0, 0.0));
                gsl_matrix_complex_set(before, 0, 1, gsl_complex_rect(0.0, 0.0));
            }
            for (j=1; j<NQ; j++) {
                if (temp[j] == 1) { //|1> = (0 1)'
                    gsl_matrix_complex_set(middle, 0, 0, gsl_complex_rect(0.0, 0.0));
                    gsl_matrix_complex_set(middle, 0, 1, gsl_complex_rect(1.0, 0.0));
                } else {  //|0> = (1 0)'
                    gsl_matrix_complex_set(middle, 0, 0, gsl_complex_rect(1.0, 0.0));
                    gsl_matrix_complex_set(middle, 0, 1, gsl_complex_rect(0.0, 0.0));
                }
                gsl_matrix_complex_free(after);
                gsl_matrix_complex *after = gsl_matrix_complex_alloc(1, 2*(before->size2));
                kron(before,  middle, after);
                //save result in another vector
                gsl_matrix_complex_free(before);
                gsl_matrix_complex *before = gsl_matrix_complex_alloc(1, after->size2);
                gsl_matrix_complex_memcpy(before, after);
            }
            //       for (j=0; j<pow(2, NQ); j++)
            //         printf("%e ", GSL_REAL(gsl_matrix_complex_get(after, 0, j)));
            //       printf("\n");

            for (j=0; j<pow(2, NQ); j++)
                gsl_vector_complex_set(w_state, j, gsl_complex_add(gsl_vector_complex_get(w_state, j), gsl_matrix_complex_get(after, 0, j)));

            //reset the tensor track vector
            temp[i] = 0;
        }

        //rescale vector
        for (j=0; j<pow(2, NQ); j++)
            gsl_vector_complex_set(w_state, j, gsl_complex_rect(GSL_REAL(gsl_vector_complex_get(w_state, j))/(double)sqrt(NQ), 0.0));
        outer_product(w_state, density_random);

        //deallocate memory
        gsl_matrix_complex_free(middle);
        gsl_matrix_complex_free(before);
        gsl_matrix_complex_free(after);
        gsl_vector_complex_free(w_state);

    } else if (STATE == 9) { //random state
        random_density_matrix(density_random);
    } else if (STATE == 12 || STATE == 15) { //vary tangle value
        PSI_matrix(density_random, theta_handle);
    } else if (STATE == 14) { //middle tangle, pure state
        PSI_matrix(density_random, 1.0/sqrt(2.0)/2.0 + 0.029151) ;
    } else if (STATE == 16) { // W state
        W_state(density_random);
    }


    // if rho is the state matrix, random is a random density matrix, then
    // apply purity noise now, as
    //  density_random = (epsilon)*density_random + (1-epsilon)*random

    //scale 2 matrices, then add them
    //gsl_matrix_complex_scale(random, epsilon);
    //display_matrix(outp);

    //gsl_matrix_complex_scale(density_random, one_min_epsilon);

    //display_matrix(density_random);
    //gsl_matrix_complex_add(density_random, random);
    //gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, one_min_epsilon, outp, ones,
    //                                                  epsilon, density_random);


    //display_matrix(density_random);
    //printf("trace %e\n", GSL_REAL(trace(density_random)));
    //printf("trace %e\n", GSL_IMAG(trace(density_random)));
    //free memory
    //gsl_matrix_complex_free(random);


}

void
measure_state(gsl_matrix_complex *density_matrix, gsl_vector *n) {
    int counter[NQ], old_counter[NQ], i, nu;

    //simulate inability to create a perfect state experimentally
    gsl_complex epsilon = gsl_complex_rect(EPS, 0);
    gsl_complex one_min_epsilon = gsl_complex_rect(1.0-EPS, 0);
    gsl_matrix_complex *random = gsl_matrix_complex_alloc(pow(2,NQ), pow(2,NQ));
    gsl_matrix_complex *tmp = gsl_matrix_complex_alloc(pow(2,NQ), pow(2,NQ));
    matrix_complex_random_init(random);
    gsl_blas_zgemm(CblasConjTrans, CblasNoTrans, one, random, random, zero, tmp);
    gsl_matrix_complex_memcpy(random, tmp);
    //renormalize trace
    gsl_matrix_complex_scale(random, gsl_complex_rect(1/GSL_REAL(trace(tmp)), 0));
    //now random = ctranspose(R)*R/trace(ctranspose(R)*R)

    //density_random = (1-EPS)*rho_state + EPS*rho_random
    gsl_matrix_complex_scale(random, epsilon);
    gsl_matrix_complex_scale(density_matrix, one_min_epsilon);
    gsl_matrix_complex_add(density_matrix, random);
    gsl_matrix_complex_free(random);

    //enforce positive-definiteness:
    //gsl_matrix_complex_memcpy(tmp, density_matrix);
    //quick_and_dirty(tmp, density_matrix);

    //gsl_matrix_complex *tmp = gsl_matrix_complex_alloc(pow(2, NQ), pow(2, NQ));
    //initialize
    for (i=0; i<NQ-1; i++)
        counter[i]=0;
    counter[NQ-1]=-1;
    gsl_matrix_complex_memcpy(mu_cell[0], mu_array[0]);
    for (i=1; i<NQ; i++)
        kron(mu_cell[i-1], mu_array[0], mu_cell[i]);

    // start taking measurements
    for (nu=0; nu<pow(4, NQ); nu++) {

        //shift counter
        base_4_rep_inc(counter, nu);

        //generate new measurement operator
        mu_tensor(mu_cell, mu_array, old_counter, counter);
        //READ: mu_cell[NQ-1] will contain the latest entry, and it will get updated
        //the next time we run the code, so we can use it as a temporary variable

        //update counter
        for (i=0; i<NQ; i++)
            old_counter[i] = counter[i];

        //compute trace(mu_nu*density_random)
        gsl_blas_zgemm(CblasNoTrans, CblasNoTrans, one, mu_cell[NQ-1],
                       density_matrix, zero, tmp);

        //display_matrix(density_matrix);
        //printf("%e\n", GSL_REAL(trace(tmp)));
        gsl_vector_set(n, nu, GSL_REAL(trace(tmp)));
    }

    //do not scale here - N is a variable now, scale elsewhere
    //scale n = N*real(n);
    //gsl_vector_scale(n, N);

    //don't round - leave as is - gets scaled later

    //for (i=0; i<pow(4,NQ)-1; i++) {
    //  gsl_vector_set(n, i, round(gsl_vector_get(n, i)));
    //printf("%e\n", gsl_vector_get(n, i));
    //}

    //deallocate memory
    gsl_matrix_complex_free(tmp);
    //gsl_matrix_complex_free(random);
}


void
apply_noise(gsl_vector *n) {

    //temorary random number generator
    gsl_rng * r;
    if ((r = gsl_rng_alloc (gsl_rng_taus)) == NULL) {
        perror("gsl_rng_alloc");
        exit(1);
    }

    //set seed - random - use local time in microseconds
    struct timeval tv;
    struct timezone tz;
    struct tm *tm;
    gettimeofday(&tv, &tz);
    tm = localtime(&tv.tv_sec);
    if (VERBOSE)
        printf(" %d:%02d:%02d %d \n", tm->tm_hour, tm->tm_min, tm->tm_sec, (int)tv.tv_usec);
    //gsl_rng_set (r, time(NULL));
    gsl_rng_set (r, tv.tv_usec);

    //waste time - random seeds have to be spaced out, so slow down the simulation
    //   if (NQ<=2) {
    //     printf("Sleeping for 2 seconds\n");
    //     usleep(2000000);
    //   }
    //   if (NQ == 3) {
    //     printf("Sleeping for 1 second\n");
    //     usleep(1000000);
    //   }
    //higher numbers of qubits will automatically make the simulation slower,
    //giving greater spacings between random seeds


    //apply Poisson noise
    int i;
    for (i=0; i<pow(4, NQ); i++) {
        gsl_vector_set(n, i, gsl_ran_poisson(r, gsl_vector_get(n, i)));
    }

    //deallocate random number generator
    gsl_rng_free(r);
}

void
load_basis() {

    //Stokes Measurement Basis
    //gsl_matrix_complex *mu_array[4];
    mu_array[0] = gsl_matrix_complex_alloc(2, 2);
    mu_array[1] = gsl_matrix_complex_alloc(2, 2);
    mu_array[2] = gsl_matrix_complex_alloc(2, 2);
    mu_array[3] = gsl_matrix_complex_alloc(2, 2);

    gsl_complex half, zero, one, nhalf, halfc, nhalfc;
    GSL_SET_COMPLEX(&half, 0.5, 0);
    GSL_SET_COMPLEX(&zero, 0, 0);
    GSL_SET_COMPLEX(&one, 1, 0);
    GSL_SET_COMPLEX(&nhalf, -0.5, 0);
    GSL_SET_COMPLEX(&halfc, 0, 0.5);
    GSL_SET_COMPLEX(&nhalfc, 0, -0.5);
    //init mu0
    gsl_matrix_complex_set(mu_array[0], 0, 0, half);
    gsl_matrix_complex_set(mu_array[0], 0, 1, zero);
    gsl_matrix_complex_set(mu_array[0], 1, 0, zero);
    gsl_matrix_complex_set(mu_array[0], 1, 1, half);
    //init mu1
    gsl_matrix_complex_set(mu_array[1], 0, 0, one);
    gsl_matrix_complex_set(mu_array[1], 0, 1, zero);
    gsl_matrix_complex_set(mu_array[1], 1, 0, zero);
    gsl_matrix_complex_set(mu_array[1], 1, 1, zero);
    //init mu3
    gsl_matrix_complex_set(mu_array[2], 0, 0, half);
    gsl_matrix_complex_set(mu_array[2], 0, 1, nhalf);
    gsl_matrix_complex_set(mu_array[2], 1, 0, nhalf);
    gsl_matrix_complex_set(mu_array[2], 1, 1, half);
    //init mu4
    gsl_matrix_complex_set(mu_array[3], 0, 0, half);
    gsl_matrix_complex_set(mu_array[3], 0, 1, halfc);
    gsl_matrix_complex_set(mu_array[3], 1, 0, nhalfc);
    gsl_matrix_complex_set(mu_array[3], 1, 1, half);

    //   int i;
    //   for (i=0; i<4; i++) {
    //     display_matrix(mu_array[i]);
    //     printf("----------------------------------------\n");
    //   }

    sig_array[0] = gsl_matrix_complex_alloc(2, 2);
    sig_array[1] = gsl_matrix_complex_alloc(2, 2);
    sig_array[2] = gsl_matrix_complex_alloc(2, 2);
    sig_array[3] = gsl_matrix_complex_alloc(2, 2);
    //init sig0
    gsl_matrix_complex_set(sig_array[0], 0, 0, one);
    gsl_matrix_complex_set(sig_array[0], 0, 1, zero);
    gsl_matrix_complex_set(sig_array[0], 1, 0, zero);
    gsl_matrix_complex_set(sig_array[0], 1, 1, one);
    //init sig1
    gsl_matrix_complex_set(sig_array[1], 0, 0, zero);
    gsl_matrix_complex_set(sig_array[1], 0, 1, one);
    gsl_matrix_complex_set(sig_array[1], 1, 0, one);
    gsl_matrix_complex_set(sig_array[1], 1, 1, zero);
    gsl_complex imag;
    GSL_SET_COMPLEX(&imag, 0, 1);
    //init sig2
    gsl_matrix_complex_set(sig_array[2], 0, 0, zero);
    gsl_matrix_complex_set(sig_array[2], 0, 1, gsl_complex_negative(imag));
    gsl_matrix_complex_set(sig_array[2], 1, 0, imag);
    gsl_matrix_complex_set(sig_array[2], 1, 1, zero);
    //init sig3
    gsl_matrix_complex_set(sig_array[3], 0, 0, one);
    gsl_matrix_complex_set(sig_array[3], 0, 1, zero);
    gsl_matrix_complex_set(sig_array[3], 1, 0, zero);
    gsl_matrix_complex_set(sig_array[3], 1, 1, gsl_complex_negative(one));


    //return mu_array;
}

---