Quantum.cpp

/*
Quantum.cpp
Written by Walter O. Krawec
 
Copyright (c) 2016 Walter O. Krawec
 
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
 
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
 
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
 
#include "Quantum.h"
 
namespace quantum
{
void Basis::orthoganalize()
{
    if(B.empty()) return;
    std::vector < std::vector <Complex> > u;
    u.push_back(B[0]);
    for(unsigned int i=1; i<B.size(); ++i)
    {
        u.push_back(B[i]);
        Complex top, bottom, p, c;
        for(unsigned int k=0; k<i; ++k)
        {
            top = innerProduct(B[i], u[k]);
            p = innerProduct(u[k], u[k]);
            bottom = algebra::ComplexHelper::invert(p);
            p = algebra::ComplexHelper::mult(top, bottom);
            for(unsigned int j=0; j<B[i].size(); ++j)
            {
                c = algebra::ComplexHelper::mult(p, u[k][j]);
                u[i][j] = u[i][j] - c;
            }
        }
    }
 
    for(unsigned int i=0; i<B.size(); ++i)
    {
        double n = norm(u[i]);
        if(n > 0.00000001)
            n = 1.0f/n;
        B[i] = u[i];
        for(unsigned int j=0; j<B[i].size(); ++j)
        {
          B[i][j] *= n;
        }
    }
}
 
Complex KetBra::trace(std::vector <Basis>& basis)
{
    unsigned int basisIndex = basis.size()+1;
    for(unsigned int i=0; i<basis.size(); ++i)
    {
        if(!basis[i].B.empty())
            basisIndex = i;
    }
 
    Complex tr=1, zero = 0;
 
    for(unsigned int i=0; i<ket.size(); ++i)
    {
        //if(basis[i].B.empty() && ket[i] != bra[i])
        //  return zero;
        if(i != ket.size()-1 && ket[i] != bra[i])
            return zero;
        if(i == ket.size()-1)
        {
            Complex sum=0;
 
            for(unsigned int j=0; j<basis[basis.size()-1].B[ket[i]].size(); ++j)
            {
                //std::cout << "<" << bra[i] << "|" << ket[i] << ">";
                    Complex b = conj(basis[basis.size()-1].B[ket[i]][j]);
 
                Complex b2 = algebra::ComplexHelper::mult(basis[basis.size()-1].B[bra[i]][j], b);
                //b = algebra::ComplexHelper::mult(tr, b2);
                sum += b2;
            }
            Complex b;
            b = algebra::ComplexHelper::mult(tr, sum);
            tr = b;
        }
    }
 
    Complex c = algebra::ComplexHelper::mult(p, tr);
    tr = c;
    //if(tr.imag > 0.00001 || tr.imag < -0.00001) throw (std::string("trace imaginary!"));
    return tr;
}
Complex KetBra::trace()
{
 
    Complex tr=1, zero = 0;
 
    for(unsigned int i=0; i<ket.size(); ++i)
    {
        if(ket[i] != bra[i])
            return zero;
    }
 
    Complex c = algebra::ComplexHelper::mult(p, tr);
    tr = c;
    return tr;
}
 
 
void DensityList::applyOp(int space)
{
    // assumes space+1 is set to |0>
    // will map |00> -> |00> + |11> + |22> + ... + |n-1n-1>
    // |10> -> |0n> + ...
 
    // assumes dim(space+1) >= dim(space)^2
 
    std::list <KetBra>::iterator Iter;
    std::list <KetBra> newDensity;
    for(Iter=density.begin(); Iter != density.end(); ++Iter)
    {
        for(int i=0; i<dimension[space]; ++i)
        {
            for(int j=0; j<dimension[space]; ++j)
            {
                KetBra kb;
                kb = *Iter;
                kb.ket[space] = i;
                kb.ket[space+1] = Iter->ket[space]*dimension[space]+i;
                kb.bra[space] = j;
                kb.bra[space+1] = Iter->bra[space]*dimension[space]+j;
                newDensity.push_back(kb);
            }
        }
    }
 
    density = newDensity;
}
 
void DensityList::conditionalSet(std::vector <int>& targetSpace, std::vector <int>& targetVal, std::vector <int>& condSpace, std::vector <int>& condVal)
{
    std::list <KetBra>::iterator Iter;
    for(Iter = density.begin(); Iter != density.end(); ++Iter)
    {
        bool val = true;
        for(int i=0; i<(int)condSpace.size(); ++i)
        {
            if(Iter->ket[condSpace[i]] != condVal[i])
                val = false;
        }
        if(val)
        {
            for(int i=0; i<(int)targetSpace.size(); ++i)
                Iter->ket[targetSpace[i]] = targetVal[i];
        }
 
        val = true;
        for(int i=0; i<(int)condSpace.size(); ++i)
        {
            if(Iter->bra[condSpace[i]] != condVal[i])
                val = false;
        }
        if(val)
        {
            for(int i=0; i<(int)targetSpace.size(); ++i)
                Iter->bra[targetSpace[i]] = targetVal[i];
        }
    }
}
 
void DensityList::applyOp(int startSpace, int endSpace)
{
    // maps |0..0> -> |0..00> + |0..11> + ... etc. etc.
    // first compute N = dim(input space)
    int N = 1;
    for(int i=startSpace; i<endSpace; ++i)
        N *= dimension[i];
 
    std::list <KetBra>::iterator Iter;
    std::list <KetBra> newDensity;
    for(Iter=density.begin(); Iter != density.end(); ++Iter)
    {
        // first compute |index> where |i> lives in input_space:
        int indexKet=0, indexBra = 0, dim=1;
        for(int i=startSpace; i<endSpace; ++i)
        {
            indexKet += dim*Iter->ket[i];
            indexBra += dim*Iter->bra[i];
            dim *= dimension[i];
        }
        std::vector <int> k,b;
        k.resize(endSpace-startSpace);
        b.resize(endSpace-startSpace);
        for(int i=0; i<endSpace-startSpace; ++i)
        {
            k[i] = 0;
            b[i] = 0;
        }
        for(int i=0; i<N; ++i)
        {
            for(int a=0; a<endSpace-startSpace; ++a)
                b[a] = 0;
 
            for(int j=0; j<N; ++j)
            {
                KetBra kb;
                kb = *Iter;
                for(int m=startSpace; m<endSpace; ++m)
                {
                    kb.ket[m] = k[m-startSpace];
                    kb.bra[m] = b[m-startSpace];
                }
                kb.ket[endSpace] = indexKet*N+i;
                kb.bra[endSpace] = indexBra*N+j;
                newDensity.push_back(kb);
 
                for(int m=0; m<endSpace-startSpace; ++m)
                {
                    b[m] ++;
                    if(b[m] >= dimension[m+startSpace])
                        b[m] = 0;
                    else
                        break;
                }
            }
            for(int m=0; m<endSpace-startSpace; ++m)
            {
                k[m] ++;
                if(k[m] >= dimension[m+startSpace])
                    k[m] = 0;
                else
                    break;
            }
        }
    }
 
    density = newDensity;
}
 
void DensityList::trace(int space, DensityList& output)
{
    std::list <KetBra>::iterator Iter;
    output = *this;
    output.density.clear();
 
    for(Iter=density.begin(); Iter != density.end(); ++Iter)
    {
        if(Iter->ket[space] == Iter->bra[space])
        {
            KetBra kb = *Iter;
            kb.ket[space] = -1;
            kb.bra[space] = -1;
            output.density.push_back(kb);
        }
    }
}
 
void DensityList::trace(int space)
{
    DensityList output;
    std::list <KetBra>::iterator Iter;
    output = *this;
    output.density.clear();
 
    for(Iter=density.begin(); Iter != density.end(); ++Iter)
    {
        if(Iter->ket[space] == Iter->bra[space])
        {
            KetBra kb = *Iter;
            kb.ket[space] = -1;
            kb.bra[space] = -1;
            output.density.push_back(kb);
        }
    }
    density = output.density;
}
 
void DensityList::applyOp(int space1, int space2, Basis& op)
{
    std::list <KetBra>::iterator Iter;
    std::list <KetBra> newDensity;
 
    for(Iter=density.begin(); Iter != density.end(); ++Iter)
    {
        std::vector < std::vector <int> > newKets, newBras;
        std::vector <Complex> pk, pb;
 
        // first ket side:
        // compute index:
        int index = Iter->ket[space1] * dimension[space2] + Iter->ket[space2];
        // now expand:
        int c1=0, c2=0;
        for(unsigned int i=0; i<op.B[index].size(); ++i)
        {
            std::vector <int> k = Iter->ket;
            k[space1] = c1;
            k[space2] = c2;
            pk.push_back(op.B[index][i]);
            newKets.push_back(k);
            ++c2;
            if(c2 >= dimension[space2])
            {
                c2 = 0;
                ++c1;
            }
        }
 
        // now bra side:
        index = Iter->bra[space1] * dimension[space2] + Iter->bra[space2];
        // now expand:
        c1=0;
        c2=0;
        for(unsigned int i=0; i<op.B[index].size(); ++i)
        {
            std::vector <int> b = Iter->bra;
            b[space1] = c1;
            b[space2] = c2;
            Complex c = conj(op.B[index][i]);
 
            pb.push_back(c);
            newBras.push_back(b);
            ++c2;
            if(c2 >= dimension[space2])
            {
                c2 = 0;
                ++c1;
            }
        }
 
        // now combine:
        for(unsigned int i=0; i<newKets.size(); ++i)
        {
            for(unsigned int j=0; j<newBras.size(); ++j)
            {
                KetBra kb;
                kb.ket = newKets[i];
                kb.bra = newBras[j];
                Complex c;
                c = algebra::ComplexHelper::mult(pk[i], pb[j]);
                kb.p = algebra::ComplexHelper::mult(Iter->p, c);
                if(creal(kb.p) != 0 || cimag(kb.p) != 0)
                    newDensity.push_back(kb);
            }
        }
    }
 
    density = newDensity;
}
 
void DensityList::applyOp(int space1, int space2, int space3, Basis& op)
{
    std::list <KetBra>::iterator Iter;
    std::list <KetBra> newDensity;
 
    for(Iter=density.begin(); Iter != density.end(); ++Iter)
    {
        std::vector < std::vector <int> > newKets, newBras;
        std::vector <Complex> pk, pb;
 
        // first ket side:
        // compute index:
        int index = Iter->ket[space1] * dimension[space2] * dimension[space3] + Iter->ket[space2] * dimension[space3] + Iter->ket[space3];
        // now expand:
        int c1=0, c2=0, c3=0;
        for(unsigned int i=0; i<op.B[index].size(); ++i)
        {
            std::vector <int> k = Iter->ket;
            k[space1] = c1;
            k[space2] = c2;
            k[space3] = c3;
            pk.push_back(op.B[index][i]);
            newKets.push_back(k);
            ++c3;
            if(c3 >= dimension[space3])
            {
                c3 = 0;
                ++c2;
                if(c2 >= dimension[space2])
                {
                    c2 = 0;
                    ++c1;
                }
            }
        }
 
        // now bra side:
        index = Iter->bra[space1] * dimension[space2] * dimension[space3] + Iter->bra[space2] * dimension[space3] + Iter->bra[space3];
        // now expand:
        c1=0;
        c2=0;
        c3=0;
        for(unsigned int i=0; i<op.B[index].size(); ++i)
        {
            std::vector <int> b = Iter->bra;
            b[space1] = c1;
            b[space2] = c2;
            b[space3] = c3;
            Complex c = conj(op.B[index][i]);
 
            pb.push_back(c);
            newBras.push_back(b);
            ++c3;
            if(c3 >= dimension[space3])
            {
                c3 = 0;
                ++c2;
                if(c2 >= dimension[space2])
                {
                    c2 = 0;
                    ++c1;
                }
            }
        }
 
        // now combine:
        for(unsigned int i=0; i<newKets.size(); ++i)
        {
            for(unsigned int j=0; j<newBras.size(); ++j)
            {
                KetBra kb;
                kb.ket = newKets[i];
                kb.bra = newBras[j];
                Complex c;
                c = algebra::ComplexHelper::mult(pk[i], pb[j]);
                kb.p = algebra::ComplexHelper::mult(Iter->p, c);
 
                if(fabs(creal(kb.p)) >= 0.0000001 || fabs(cimag(kb.p)) >= 0.0000001)
                    newDensity.push_back(kb);
            }
        }
    }
 
    density = newDensity;
}
 
void DensityList::applyConditionalOp(int condSpace, int cond, int space, algebra::mat& U)
{
    std::list <KetBra>::iterator Iter;
 
    std::list <KetBra> newDensity;
    Complex one=1;
 
    for(Iter=density.begin(); Iter != density.end(); ++Iter)
    {
        std::vector < std::vector <int> > ket, bra;
        std::vector <Complex> pk, pb;
        if(Iter->ket[condSpace] == cond)
        {
            std::vector <int> k = Iter->ket;
            for(int r=0; r<U.getRows(); ++r)
            {
                k[space] = r;
                    Complex c = U(r, Iter->ket[space]);
                pk.push_back(c);
                ket.push_back(k);
            }
        }
        else
        {
            ket.push_back(Iter->ket);
            pk.push_back(one);
        }
 
        if(Iter->bra[condSpace] == cond)
        {
            std::vector <int> b = Iter->bra;
            for(int r=0; r<U.getRows(); ++r)
            {
                b[space] = r;
                Complex c = conj(U(r, Iter->bra[space]));
 
                pb.push_back(c);
                bra.push_back(b);
            }
        }
        else
        {
            bra.push_back(Iter->bra);
            pb.push_back(one);
        }
 
        for(unsigned int i=0; i<ket.size(); ++i)
        {
            for(unsigned int j=0; j<bra.size(); ++j)
            {
                KetBra kb;
                kb.ket = ket[i];
                kb.bra = bra[j];
                Complex c;
                c = algebra::ComplexHelper::mult(pk[i], pb[j]);
                kb.p = algebra::ComplexHelper::mult(Iter->p, c);
                if(creal(kb.p) != 0 || cimag(kb.p) != 0)
                    newDensity.push_back(kb);
            }
        }
    }
 
    density = newDensity;
}
 
void DensityList::applyOp(int space, algebra::mat& U)
{
    std::list <KetBra>::iterator Iter;
 
    std::list <KetBra> newDensity;
    Complex one=1;
    for(Iter=density.begin(); Iter != density.end(); ++Iter)
    {
        std::vector < std::vector <int> > ket, bra;
        std::vector <Complex> pk, pb;
        //if(Iter->ket[condSpace] == cond)
        {
            std::vector <int> k = Iter->ket;
            for(int r=0; r<U.getRows(); ++r)
            {
                k[space] = r;
                Complex c = U(r, Iter->ket[space]);
                pk.push_back(c);
                ket.push_back(k);
            }
        }
        /*else
        {
            ket.push_back(Iter->ket);
            pk.push_back(one);
        }*/
 
        //if(Iter->bra[condSpace] == cond)
        {
            std::vector <int> b = Iter->bra;
            for(int r=0; r<U.getRows(); ++r)
            {
                b[space] = r;
                Complex c = conj(U(r, Iter->bra[space]));
 
                pb.push_back(c);
                bra.push_back(b);
            }
        }
        /*else
        {
            bra.push_back(Iter->bra);
            pb.push_back(one);
        }*/
 
        for(unsigned int i=0; i<ket.size(); ++i)
        {
            for(unsigned int j=0; j<bra.size(); ++j)
            {
                KetBra kb;
                kb.ket = ket[i];
                kb.bra = bra[j];
                Complex c;
                c = algebra::ComplexHelper::mult(pk[i], pb[j]);
                kb.p = algebra::ComplexHelper::mult(Iter->p, c);
                if(creal(kb.p) != 0 || cimag(kb.p) != 0)
                    newDensity.push_back(kb);
            }
        }
    }
 
    density = newDensity;
}
 
// WK: Added for GA+prac-QKD project - very much a hack specific to this project
// assumes transitSpace is the index for the T wire; T+1 is the E guess wire (must be dim 2)
// and T+2 is arbitrary dimension (Eaux)
void DensityList::applyOp(algebra::mat& U, int transitSpace)
{
    std::list <KetBra>::iterator Iter;
 
    std::list <KetBra> newDensity;
    Complex one=1;
    for(Iter=density.begin(); Iter != density.end(); ++Iter)
    {
        std::vector < std::vector <int> > ket, bra;
        std::vector <Complex> pk, pb;
        {
          std::vector <int> k = Iter->ket;
          int n = dimension[transitSpace+2];
          int a = Iter->ket[transitSpace+2];
          int g = Iter->ket[transitSpace+1];
          int t = Iter->ket[transitSpace];
          int col = a+g*n+2*n*t;
 
          for(int a2 = 0; a2<n; ++a2){
            for(int g2=0; g2<2; ++g2){
              for(int t2=0; t2<2; ++t2){
            int row = a2+g2*n+2*n*t2;
            Complex c = U(row, col);
            k[transitSpace] = t2;
            k[transitSpace+1] = g2;
            k[transitSpace+2] = a2;
            pk.push_back(c);
            ket.push_back(k);
              }
            }
          }
        }
 
        {
          std::vector <int> b = Iter->bra;
          int n = dimension[transitSpace+2];
          int a = Iter->bra[transitSpace+2];
          int g = Iter->bra[transitSpace+1];
          int t = Iter->bra[transitSpace];
          int col = a+g*n+2*n*t;
 
          for(int a2 = 0; a2<n; ++a2){
            for(int g2=0; g2<2; ++g2){
              for(int t2=0; t2<2; ++t2){
            int row = a2+g2*n+2*n*t2;
            Complex c = conj(U(row, col));
            b[transitSpace] = t2;
            b[transitSpace+1] = g2;
            b[transitSpace+2] = a2;
            pb.push_back(c);
            bra.push_back(b);
              }
            }
          }
        }
 
        for(unsigned int i=0; i<ket.size(); ++i)
        {
            for(unsigned int j=0; j<bra.size(); ++j)
            {
                KetBra kb;
                kb.ket = ket[i];
                kb.bra = bra[j];
                Complex c;
                c = algebra::ComplexHelper::mult(pk[i], pb[j]);
                kb.p = algebra::ComplexHelper::mult(Iter->p, c);
                if(creal(kb.p) != 0 || cimag(kb.p) != 0)
                    newDensity.push_back(kb);
            }
        }
    }
 
    density = newDensity;
}
 
 
void DensityList::project(int space, int state, DensityList& output)
{
    std::list <KetBra>::iterator Iter;
    output = *this;
    output.density.clear();
 
    for(Iter=density.begin(); Iter != density.end(); ++Iter)
    {
        if(Iter->ket[space] == Iter->bra[space] && Iter->ket[space] == state)
        {
            KetBra kb = *Iter;
            output.density.push_back(kb);
        }
    }
}
 
void DensityList::project(int space, int state)
{
    std::list <KetBra>::iterator Iter;
    std::list <KetBra> newDensity;
 
    for(Iter=density.begin(); Iter != density.end(); ++Iter)
    {
        if(Iter->ket[space] == Iter->bra[space] && Iter->ket[space] == state)
        {
            KetBra kb = *Iter;
            newDensity.push_back(kb);
        }
    }
 
    density = newDensity;
}
 
void DensityList::applyHadamard(int space)
{
  algebra::mat H;
  H.create(2);
  H(0,0) = 1.0/sqrt(2.0);
  H(0,1) = 1.0/sqrt(2.0);
  H(1,0) = 1.0/sqrt(2.0);
  H(1,1) = -1.0/sqrt(2.0);
  applyOp(space, H);
}
 
void DensityList::applyRotation(int space, double rotation)
{
  // this might be more efficient if we just go through list...
  algebra::mat R;
  R.create(2);
  R(0,0) = 1;
  R(0,1) = 0;
  R(1,0) = 0;
  R(1,1) = cexp(rotation);
  applyOp(space, R);
}
 
void DensityList::applySU2(int space, double p, double theta, double psi)
{
  if(p > 1)
    p = 1;
  if(p < 0)
    p = 0;
  // this might be more efficient if we just go through list...
  algebra::mat R;
  R.create(2);
  R(0,0) = sqrt(p)*cexp(I*theta);
  R(0,1) = -sqrt(1-p)*cexp(-I*psi);
  R(1,0) = sqrt(1-p)*cexp(I*psi);
  R(1,1) = sqrt(p)*cexp(-I*theta);
 
  //R.print(std::cout);
  applyOp(space, R);
}
 
 
void DensityList::applyHadamard(int target, int control)
{
  algebra::mat H;
  H.create(2);
  H(0,0) = 1.0/sqrt(2.0);
  H(0,1) = 1.0/sqrt(2.0);
  H(1,0) = 1.0/sqrt(2.0);
  H(1,1) = -1.0/sqrt(2.0);
  if(control >= 0)
    applyConditionalOp(control, 1, target, H);
  else
    applyOp(target, H);
}
 
void DensityList::applyPiOver8(int target, int control)
{
  algebra::mat H;
  H.create(2);
  H(0,0) = 1.0;
  H(0,1) = 0;
  H(1,0) = 0;
  H(1,1) = cos(3.1415/4.0) + I*sin(3.1415/4.0);
  if(control >= 0)
    applyConditionalOp(control, 1, target, H);
  else
    applyOp(target, H);
}
 
 
void DensityList::applySU2(int target, int control, double p, double theta, double psi)
{
  if(p > 1)
    p = 1;
  if(p < 0)
    p = 0;
  // this might be more efficient if we just go through list...
  algebra::mat R;
  R.create(2);
  R(0,0) = sqrt(p)*cexp(I*theta);
  R(0,1) = -sqrt(1-p)*cexp(-I*psi);
  R(1,0) = sqrt(1-p)*cexp(I*psi);
  R(1,1) = sqrt(p)*cexp(-I*theta);
 
  //R.print(std::cout);
  if(control >= 0)
    applyConditionalOp(control, 1, target, R);
  else
    applyOp(target, R);
}
 
void DensityList::applyNOT(int target, int control)
{
  if(control >= 0)
    CNOT(control, target);
  else
    NOT(target);
}
 
 
 
 
void DensityList::CNOT(int controlSpace, int targetSpace)
{
  if(controlSpace == targetSpace)
    return;
  std::list <KetBra>::iterator Iter;
  //std::list <KetBra> newDensity;
  
  for(Iter=density.begin(); Iter != density.end(); ++Iter)
    {
      if(Iter->ket[controlSpace] == 1)
    Iter->ket[targetSpace] = 1-Iter->ket[targetSpace];
      if(Iter->bra[controlSpace] == 1)
    Iter->bra[targetSpace] = 1-Iter->bra[targetSpace];
    }
}
 
void DensityList::NOT(int space)
{
    std::list <KetBra>::iterator Iter;
    //std::list <KetBra> newDensity;
 
    for(Iter=density.begin(); Iter != density.end(); ++Iter)
    {
        Iter->ket[space] = 1-Iter->ket[space];
        Iter->bra[space] = 1-Iter->bra[space];
    }
}
 
void DensityList::splitAndSave(double p1, double p2, int space)
{
    std::list <KetBra>::iterator Iter;
    std::list <KetBra> newDensity;
 
    for(Iter = density.begin(); Iter != density.end(); ++Iter)
    {
        KetBra kb;
        kb = *Iter;
        kb.p *= p1;
 
        newDensity.push_back(kb);
        kb.ket[space] = 1;
        kb.bra[space] = 1;
        newDensity.push_back(kb);
    }
 
    density = newDensity;
}
void DensityList::swap(int space1, int space2)
{
    std::list <KetBra>::iterator Iter;
 
    for(Iter = density.begin(); Iter != density.end(); ++Iter)
    {
        int temp = Iter->ket[space1];
        Iter->ket[space1] = Iter->ket[space2];
        Iter->ket[space2] = temp;
 
        temp = Iter->bra[space1];
        Iter->bra[space1] = Iter->bra[space2];
        Iter->bra[space2] = temp;
    }
}
 
double DensityList::entropy()
{
  algebra::mat M;
  computeDensityMat(M);
  return M.entropy();
}
 
double DensityList::ShannonEntropy()
{
  algebra::mat M;
  computeDensityMat(M);
  double sum = 0;
  for(int i=0; i<M.getRows(); ++i){
    if(fabs(cimag(M(i,i))) > 0.0000001){
      std::cout << "Density List Error: Negative eigenvalue.\n";
      throw std::string("Error");
    }
    double e = creal( M(i,i) );
    sum -= e*SafeLog(e);
  }
  return sum;
}
 
void DensityList::measureAndSave(int measureSpace, int saveSpace)
{
    // assumes dim[saveSpace] >= dim[measureSpace]
    std::list <KetBra>::iterator Iter;
    std::list <KetBra> newDensity;
 
    for(Iter = density.begin(); Iter != density.end(); ++Iter)
    {
        if(Iter->ket[measureSpace] == Iter->bra[measureSpace])
        {
            KetBra kb = *Iter;
            kb.ket[saveSpace] = Iter->ket[measureSpace];
            kb.bra[saveSpace] = kb.ket[saveSpace];
            newDensity.push_back(kb);
        }
    }
 
    density = newDensity;
}
 
void DensityList::measure(int measureSpace)
{
    // assumes dim[saveSpace] >= dim[measureSpace]
    std::list <KetBra>::iterator Iter;
    std::list <KetBra> newDensity;
 
    for(Iter = density.begin(); Iter != density.end(); ++Iter)
    {
        if(Iter->ket[measureSpace] == Iter->bra[measureSpace])
        {
            KetBra kb = *Iter;
            newDensity.push_back(kb);
        }
    }
 
    density = newDensity;
}
 
double DensityList::calculatePr(int measureSpace, int state)
{
    // first |state><state| rho |state><state|:
    std::list <KetBra>::iterator Iter;
    std::list <KetBra> tempDensity;
 
    for(Iter = density.begin(); Iter != density.end(); ++Iter)
    {
        if(Iter->ket[measureSpace] == state && Iter->ket[measureSpace] == Iter->bra[measureSpace])
        {
            KetBra kb = *Iter;
            tempDensity.push_back(kb);
        }
    }
 
    if(tempDensity.empty())
        return 0.0;
 
    // now calculate tr(tempDensity)
    Complex p=0, c;
 
    for(Iter = tempDensity.begin(); Iter != tempDensity.end(); ++Iter)
    {
        c = Iter->trace();
        p += c;
    }
 
    // now calculate tr(density)
    Complex p2=0;
 
    for(Iter = density.begin(); Iter != density.end(); ++Iter)
    {
        c = Iter->trace();
        p2 += c;
    }
 
    c = algebra::ComplexHelper::invert(p2);
    p2 = algebra::ComplexHelper::mult(c, p);
 
    if(cimag(p2) > 0.00000001 || cimag(p2) < -0.00000001)
        throw std::string("Trace error.");
    return creal(p2);
}
 
double DensityList::trace()
{
  std::list <KetBra>::iterator Iter;
  Complex p2=0;
 
  for(Iter = density.begin(); Iter != density.end(); ++Iter)
    p2 += Iter->trace();;
  
  if(cimag(p2) > 0.00000001 || cimag(p2) < -0.00000001)
    throw std::string("Trace error.");
  return creal(p2);
}
 
double DensityList::calculatePr(int space1, int state1, int space2, int state2)
{
    // first |state><state| rho |state><state|:
    std::list <KetBra>::iterator Iter;
    std::list <KetBra> tempDensity;
 
    for(Iter = density.begin(); Iter != density.end(); ++Iter)
    {
        if(Iter->ket[space1] == state1 && Iter->ket[space1] == Iter->bra[space1] &&
            Iter->ket[space2] == state2 && Iter->ket[space2] == Iter->bra[space2])
        {
            KetBra kb = *Iter;
            tempDensity.push_back(kb);
        }
    }
 
    if(tempDensity.empty())
        return 0.0;
 
    // now calculate tr(tempDensity)
    Complex p=0, c;
 
    for(Iter = tempDensity.begin(); Iter != tempDensity.end(); ++Iter)
    {
      c = Iter->trace();
      p += c;
    }
 
    // now calculate tr(density)
    Complex p2=0;
 
    for(Iter = density.begin(); Iter != density.end(); ++Iter)
    {
      c = Iter->trace();
      p2 += c;
    }
 
    c = algebra::ComplexHelper::invert(p2);
    p2 = algebra::ComplexHelper::mult(c, p);
 
    if(cimag(p2) > 0.00000001 || cimag(p2) < -0.00000001)
        throw std::string("Trace error.");
    return creal(p2);
}
 
double DensityList::calculatePr(const std::vector <int>& spaceList, const std::vector <int>& stateList)
{
    // first |state><state| rho |state><state|:
    std::list <KetBra>::iterator Iter;
    std::list <KetBra> tempDensity;
 
    for(Iter = density.begin(); Iter != density.end(); ++Iter)
    {
        bool measure = true;
        for(unsigned int i=0; i<spaceList.size(); ++i)
        {
            if(Iter->ket[spaceList[i]] != stateList[i])
            {
                measure = false;
                break;
            }
        }
        if(measure)
        {
            KetBra kb = *Iter;
            tempDensity.push_back(kb);
        }
    }
 
    if(tempDensity.empty())
        return 0.0;
 
    // now calculate tr(tempDensity)
    Complex p=0, c;
 
    for(Iter = tempDensity.begin(); Iter != tempDensity.end(); ++Iter)
    {
        c = Iter->trace();
        p += c;
    }
 
    // now calculate tr(density)
    Complex p2=0;
 
    for(Iter = density.begin(); Iter != density.end(); ++Iter)
    {
        c = Iter->trace();
        p2 += c;
    }
 
    c = algebra::ComplexHelper::invert(p2);
    p2 = algebra::ComplexHelper::mult(c, p);
 
    if(cimag(p2) > 0.00000001 || cimag(p2) < -0.00000001)
        throw std::string("Trace error.");
    return creal(p2);
}
 
void DensityList::computeDensityMatFast(algebra::mat& output)
{
    int n = 1;
    for(unsigned int i=0; i<density.begin()->ket.size(); ++i)
    {
        if(density.begin()->ket[i] != -1)
            n *= dimension[i];
    }
 
    int basisIndex = dimension.size()-1;
    // basisIndex = last index not traced out!
    for(basisIndex=dimension.size()-1; basisIndex>=0; --basisIndex){
      if(density.begin()->ket[basisIndex] != -1)
        break;
    }
    output.create(n);
    output.zero();
 
    //int count = density.size();
    //int counter = 0;
 
    std::list <KetBra>::iterator Iter;
    for(Iter = density.begin(); Iter != density.end(); ++Iter)
    {
        //if(count > 5000 && counter%500 == 0)
        //  std::cout << "\t" << counter << "/" << count << "\n";
        //++counter;
        // change basis to be |i>|e><j|<e|
        int indexI=0, indexJ=0;
        for(unsigned int i=0; i<Iter->ket.size(); ++i)
        {
            if(Iter->ket[i] == -1 || i == basisIndex) continue;
            indexI = indexI*dimension[i] + Iter->ket[i];
            indexJ = indexJ*dimension[i] + Iter->bra[i];
        }
 
        int startRow = indexI*dimension[basisIndex];
        int startCol = indexJ*dimension[basisIndex];
        int i = Iter->ket[basisIndex];
        int j = Iter->bra[basisIndex];
 
        output(startRow+i, startCol+j) += Iter->p;
    }
 
    Complex tr = output.trace();
    Complex tr2 = algebra::ComplexHelper::invert(tr);
    output.multiplyScalar(tr2);
}
 
void DensityList::computeDensityMat(algebra::mat& output)
{
    // first determine size:
    int n = 1;
    for(unsigned int i=0; i<density.begin()->ket.size(); ++i)
    {
        if(density.begin()->ket[i] != -1)
            n *= dimension[i];
    }
 
 
    output.create(n);
    output.zero();
 
    std::list <KetBra>::iterator Iter;
    algebra::mat k, b, t, t2, temp, temp2;
 
    int count = density.size();
    int counter = 0;
 
    //double smallRe = fabs(density.begin()->p.real);
    for(Iter = density.begin(); Iter != density.end(); ++Iter)
    {
        //if(count > 5000 && counter%500 == 0)
        //  std::cout << "\t" << counter << "/" << count << "\n";
        //if(fabs(Iter->p.real) < smallRe)
        //  smallRe = fabs(Iter->p.real);
        ++counter;
        k.create(1, 1);
        b.create(1, 1);
        k(0,0) = 1.0;
        b(0,0) = 1.0;
        for(unsigned int i=0; i<Iter->ket.size(); ++i)
        {
            if(Iter->ket[i] == -1) continue;
 
            t.create(dimension[i], 1);
            t2.create(1, dimension[i]);
            t.zero();
            t(Iter->ket[i], 0) = 1.0;
 
            k.tensor(&t, &temp);
            k = temp;
 
            t2.zero();
            t2(0, Iter->bra[i]) = 1.0;
 
            b.tensor(&t2, &temp);
            b = temp;
 
        }
 
        k.multiply(b, temp);
        temp.multiplyScalar(Iter->p);
 
        output.add(&temp, &temp2);
        output = temp2;
    }
 
    Complex tr = output.trace();
    Complex tr2 = algebra::ComplexHelper::invert(tr);
    output.multiplyScalar(tr2);
 
    //std::cout << "\t\tsmallRe = " << smallRe << "\n";
}
 
}