// Hilbert class implementation file

#include <stdafx.h>
#include "resource.h"

#include <math.h>
#include <stdlib.h>
#include <assert.h>
#include <iomanip.h>
#include <time.h>

#include "hilbert.h"

// constructors

CHilbert::CHilbert(int nqubits)	
{
	n_qbs = nqubits;
	n_bs = (int)pow(2, n_qbs);
	
	stateamp = new CComplex[n_bs];	// default complex constructor called 
	assert(stateamp != 0);			// i.e all set to 0+0i
}

CHilbert::CHilbert(int nqubits, int basis_index)	
{
	n_qbs = nqubits;
	n_bs = (int)pow(2, n_qbs);
	
	stateamp = new CComplex[n_bs];	// default complex constructor called 
	assert(stateamp != 0);			// i.e all set to 0+0i

	// set to initial state
	CComplex* comp = new CComplex(1,0);
	if(basis_index < n_bs)	// check if out of range
		stateamp[basis_index] = *comp;
	else
		cout << "Hilbert initialisation problem!" << endl;

	// tidy up
	delete comp;
}

// destructor

CHilbert::~CHilbert()
{	
	delete[] stateamp;	// reclaim the memory
}

// accessor functions

int CHilbert::numbits(void)
{
	return n_qbs;
}

int CHilbert::numbs(void)
{
	return n_bs;
}

CComplex CHilbert::retamp(int index)
{
	return stateamp[index];
}

// set function

void CHilbert::setamp(int basis_index, CComplex amp)
{
	stateamp[basis_index] = amp;
}

// functions

void CHilbert::normalise(void)
{
	double norm = 0;
	double sqrtnorm = 0;

	for(int i = 0; i < n_bs; i++)
		norm += stateamp[i].magsqrd();	// add up current probabilities

	// check for zero norm -> obvious problems
	if(norm == 0)
		cout << "Overall prob = 0!" << endl;
	
	sqrtnorm = sqrt(norm);

	for(i = 0; i < n_bs; i++)
		stateamp[i] /= sqrtnorm;	// normalise

}

void CHilbert::reverse(void)
{
	CHilbert* state_vec_tmp = new CHilbert(n_qbs);	// temp state_vecert 
	*state_vec_tmp = *this;							// make copy

	int revindex;

	CBstate* basis = new CBstate(n_qbs);

	for(int i = 0; i < n_bs; i++)
	{
		basis->setbasis(i);
		revindex = basis->rev_index();	// reversed basis state

		state_vec_tmp->setamp(i, this->retamp(revindex));	
	}

	// tidy up
	delete basis;
	*this = *state_vec_tmp; 
	delete state_vec_tmp;
}

CHilbert* CHilbert::subspaceprojection(int n, int* list)
{
	// subspace projection onto an n-qubit sub space defined
	// by the n-qubit indices in list 

	CHilbert* substate_vec = new CHilbert(n);	// subspace state_vector 

	CBstate* basis = new CBstate(n_qbs);	// basis state of entire system

	int subindex;

	for(int i = 0; i < n_bs; i++)
	{
		basis->setbasis(i);	// set basis state of entire system at this index
		subindex = basis->small_index(n, list);	// index of sub space basis state 
												// defined by qubit indices in list

		// add to current probability for this sub space basis state  
		substate_vec->setamp(subindex, substate_vec->stateamp[subindex] + this->stateamp[i]);
	}

	// tidy up
	delete basis;
	return substate_vec;
}

int CHilbert::measure(int sub_n_qbs, int* qbits, THREADPARAMS* ptp)
{
	// (MFC) get activity progress object ptr
	CProgressCtrl *progress_ctrl_ptr = (CProgressCtrl*) ptp->pView->GetDlgItem(IDC_ACTIVITY_PROGRESS);
	progress_ctrl_ptr->SetPos(0);

	int n_bs_sub = (int)pow(2, sub_n_qbs);	// number of basis states spanning
									// the sub space being measured
	int value = -1;		// if value not set, -1 signals an error

	// random number in range 0 -> 1
	double randnum = (double)rand() / RAND_MAX;
	
	double accumprob = 0;	// accumulative probability

	CHilbert* state_vec_tmp = new CHilbert(n_qbs);	// temporary state_vector 
	
	// basis state of the entire space
	CBstate* fullbasis = new CBstate(n_qbs);	
	int sub_fullindex;		// index to measured subspace basis state 
							// defined by the qubit indices of the measured space
							// e.g. if entire space basis state = |10110>
							// sub space basis state defined by qubits 1 and 3
							// is |01>
	
	// basis state of the measured space 
	// defined by the qubit indices in qubits list
	CBstate* subbasis = new CBstate(sub_n_qbs);	
	int subindex;	// index to the measured subspace basis state								

	// loop through all basis states of the measured subspace
	for(int j = 0; j < n_bs_sub; j++) 
	{
		subbasis->setbasis(j);	// j'th basis state
		subindex = subbasis->indexnumber();	// index to the j'th basis state
											// just j! Here for clarity!? 

		// (MFC)
		progress_ctrl_ptr->SetPos((int)(100*(float)j/n_bs_sub));

		// sum the probabilities for the j'th measured subspace basis
		// by stepping through basis states of entire space and looking 
		// for matches 
		for(int i = 0; i < n_bs; i++)
		{
			fullbasis->setbasis(i);	// i'th basis state of the entire space 
			// get the index corresponding to the measured basis state
			sub_fullindex = fullbasis->small_index(sub_n_qbs, qbits);

			if(subindex == sub_fullindex)	// got match
				accumprob += stateamp[i].magsqrd();
		}
		// if random exceed -> entire space is restricted to 
		// this particular measured subspace basis state
		if(accumprob >= randnum)	
		{
			// do the preojection, restricting system to the
			// measured subspace basis state
			for(int k = 0; k < n_bs; k++) 
			{
				fullbasis->setbasis(k);
				sub_fullindex = fullbasis->small_index(sub_n_qbs, qbits);
				
				if(subindex == sub_fullindex)
					state_vec_tmp->setamp(k, this->retamp(k));
			}

			*this = *state_vec_tmp;
			this->normalise();		// ensure unity measurement probability
			value = j;				// your measured value

			break;	// break outer loop
		}
	}

	// tidy up
	delete subbasis;
	delete fullbasis;
	delete state_vec_tmp;

	// (MFC)
	progress_ctrl_ptr->SetPos(100);

	if(value >= 0)
		return value;
	else
	{
		cout << "Doh! accumulative probability = " << accumprob << endl;
		return value;	// value not yet set
	}
}

// modula exponentiation function application

void CHilbert::modexp(int n_qbs_reg1, int n_qbs_reg2, int x, int n)
{
	// |reg2|reg1>
	// |F(A)|A>

	// number of basis states spanning subspace defined by reg1
	int reg1_n_bs = (int)pow(2, n_qbs_reg1);	

	CComplex* comp = new CComplex(1.0);

	// loop around basis states of the subspace
	// defined by reg1 -> i.e. all values of A
	for(int i = 0; i < reg1_n_bs; i++)
	{
		// set amplitude at basis state which corresponds 
		// to the appropriate A and corresponding F(A)   
		int bigbasisindex = CModexp::modexp(x, i, n) * (int)pow(2,n_qbs_reg1) + i;
		setamp(bigbasisindex,*comp);
	}

	this->normalise();	// ensure unity measurement probability
						// (comp magnitude = 1)

	// tidy up
	delete comp;
}

ostream& operator <<(ostream& out, CHilbert& state_vec)
{
	out << setiosflags(ios::fixed) << setprecision(4);

	out << " ";
	for(int i = (state_vec.numbits() - 1); i >= 0; i--)
	{
		if(i < 10)
			out << i;
		else
			out << (char)((int)'A' + i - 10);
	}
	out << endl;

	CBstate* basis = new CBstate(state_vec.n_qbs);

	for(i = 0; i < state_vec.n_bs; i++)
	{
		basis->setbasis(i);
		out << *basis << ' ';

		out << state_vec.stateamp[i];
		out << ' ' << state_vec.stateamp[i].magsqrd() << ' ';// << endl;
		for(int j = 0; j < (int)(state_vec.stateamp[i].magsqrd()*100.0); j++)
			out << '*';

		out << endl;
	}

	delete basis;
	return out;
}

CHilbert &CHilbert::operator =(CHilbert& rhs)
{	
	for(int i = 0; i < n_bs; i++)
		this->setamp(i, rhs.stateamp[i]);

	this->n_qbs = rhs.n_qbs;
	this->n_bs = rhs.n_bs;

	return *this;
}