src/qcggs.cpp

/***************************************************************************
 *   Quantum Construct (qC++)                                                                  *
 *   The Quantum Physics Computational Library                                         *
 *   Copyright (C) 2005 by Shekhar S. Chandra                                      *
 *   Shekhar.Chandra@sci.monash,edu.au                                     *
 *                                                                         *
 *   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; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

#include "qcggs.h"

QCGGS::QCGGS()
{
    twoCompBEC = false;
}

QCGGS::~QCGGS()
{
    //dtor
}

void QCGGS::setNoOfImages(int number)
{
    noOfImages = number;
    GSPhases.resizeAndPreserve(noOfImages);
    GSPhases2.resizeAndPreserve(noOfImages);
    Images.resizeAndPreserve(noOfImages);
    Images2.resizeAndPreserve(noOfImages);
}

void QCGGS::setSystem(QCQuantumSystem<complexType,2> &system)
{
    ptrSystem = &system;
    ptrField = system.accessData();
    int rows = system.getNoOfRows(), cols = system.getNoOfColumns();
        for(int j = 0; j < noOfImages; j ++)
                GSPhases(j).setSize(rows,cols);
        for(int j = 0; j < noOfImages; j ++)
                Images(j).setSize(rows,cols);
    twoCompBEC = false;
}

void QCGGS::setSystem(QC2ComponentBEC<2> &system)
{
    ptrSystem = &system;
    ptrField = system.accessComp1();
    ptrField2 = system.accessComp2();
    int rows = system.getNoOfRows(), cols = system.getNoOfColumns();

        for(int j = 0; j < noOfImages; j ++)
        {
                GSPhases(j).setSize(rows,cols);
                GSPhases2(j).setSize(rows,cols);
        }
        for(int j = 0; j < noOfImages; j ++)
        {
                Images(j).setSize(rows,cols);
                Images2(j).setSize(rows,cols);
        }
        twoCompBEC = true;
}

void QCGGS::setGuess(parameterType mag, parameterType phase)
{
    if(ptrField != NULL)
    {
        if(!twoCompBEC)
        {
            for(int k = 0; k < ptrField->rows(); k ++)
                for(int j = 0; j < ptrField->cols(); j ++)
                    (*ptrField)(j,k) = mag*exp(i*phase);
        }
        else
        {
            for(int k = 0; k < ptrField->rows(); k ++)
                for(int j = 0; j < ptrField->cols(); j ++)
                    (*ptrField)(j,k) = (*ptrField2)(j,k) = mag*exp(i*phase);
        }
    }
    else
        exit(0);
}

void QCGGS::loadImages(string extension)
{
    string filename;
    char digit;

    for(int j = 0; j < noOfImages; j ++)
        {
            digit = static_cast<char> (j + 48);
            if(!twoCompBEC)
            {
            filename = filePrefix + digit + extension;
            Images(j).readFromFile(filename.c_str());
            }
            else
            {
            filename = filePrefix + '1' + digit + extension;
            Images(j).readFromFile(filename.c_str());
            filename = filePrefix + '2' + digit + extension;
            Images2(j).readFromFile(filename.c_str());
            }
        }
}

void QCGGS::writeResults(string extension)
{
    QCPhase Phase;
    string filename;
    char digit;

    for(int j = 0; j < noOfImages; j ++)
    {
        digit = static_cast<char> (j + 48);
        if(!twoCompBEC)
            {
            filename = filePrefix + digit + extension;
            Phase.outputPhaseFunction(GSPhases(j),filename.c_str());
            }
            else
            {
            filename = filePrefix + '1' + digit + extension;
            Phase.outputPhaseFunction(GSPhases(j),filename.c_str());
            filename = filePrefix + '2' + digit + extension;
            Phase.outputPhaseFunction(GSPhases2(j),filename.c_str());
            }
    }
}

void QCGGS::writeOriginals(string extension)
{
    QCPhase Phase;
    string filename;
    char digit;

    for(int j = 0; j < noOfImages; j ++)
    {
        digit = static_cast<char> (j + 48);
        if(!twoCompBEC)
            {
            filename = filePrefix + digit + extension;
            Phase.outputPhaseFunction(Images(j),filename.c_str());
            }
            else
            {
            filename = filePrefix + '1' + digit + extension;
            Phase.outputPhaseFunction(Images(j),filename.c_str());
            filename = filePrefix + '2' + digit + extension;
            Phase.outputPhaseFunction(Images2(j),filename.c_str());
            }
    }
}

void QCGGS::copyField(QCComplexField<2>& dest, QCComplexField<2>& source)
{
        for(int k = 0; k < source.noOfColumns(); k ++)
                for(int j = 0; j < source.noOfRows(); j ++)
                        dest.assignElement(j,k,source.retrieveElement(j,k));
}

void QCGGS::copyField(QCComplexField<2>& dest, Array<complexType,2> &source)
{
        for(int k = 0; k < source.cols(); k ++)
                for(int j = 0; j < source.rows(); j ++)
                        dest.assignElement(j,k,source(j,k));
}

void QCGGS::updateMagnitude(QCComplexField<2>& image, QCComplexField<2> &result)
{
        double currentPhase = 0.0;
        int buffer = 0/*fieldSize/5*/;

        for(int k = buffer; k < result.noOfColumns()-buffer; k ++)
                for(int j = buffer; j < result.noOfRows()-buffer; j ++)
                {
                        currentPhase = atan2(result.retrieveElement(j,k).imag(),result.retrieveElement(j,k).real());
                        result.assignElement(j,k,sqrt(norm(image.retrieveElement(j,k)))*exp(i*currentPhase));
                }
}

void QCGGS::updateMagnitude(QCComplexField<2>& image, Array<complexType,2>& result)
{
        double currentPhase = 0.0;
        int buffer = 0/*fieldSize/5*/;

        for(int k = buffer; k < result.cols()-buffer; k ++)
                for(int j = buffer; j < result.rows()-buffer; j ++)
                {
                        currentPhase = atan2(result(j,k).imag(),result(j,k).real());
                        result(j,k) = sqrt(norm(image.retrieveElement(j,k)))*exp(i*currentPhase);
                }
}

parameterType QCGGS::RMSError(QCComplexField<2>& original, QCComplexField<2> &result)
{
        parameterType sigma = 0.0, totalNorm = 0.0, tmp = 0.0;
        int frame = result.noOfRows()/10;

        for(int k = frame; k < original.noOfColumns()-frame; k ++)
                for(int j = frame; j < original.noOfRows()-frame; j ++)
                {
                        tmp = sqrt(norm(original.retrieveElement(j,k))) - sqrt(norm(result.retrieveElement(j,k)));
                        sigma = sigma + tmp*tmp;
                        totalNorm = totalNorm + norm(original.retrieveElement(j,k));
                }

        sigma = sqrt(sigma/totalNorm);
        return sigma;
}

parameterType QCGGS::RMSError(QCComplexField<2>& original, Array<complexType,2> &result)
{
        parameterType sigma = 0.0, totalNorm = 0.0, tmp = 0.0;
        int frame = result.rows()/10;

        for(int k = frame; k < original.noOfColumns()-frame; k ++)
                for(int j = frame; j < original.noOfRows()-frame; j ++)
                {
                        tmp = sqrt(norm(original.retrieveElement(j,k))) - sqrt(norm(result(j,k)));
                        sigma = sigma + tmp*tmp;
                        totalNorm = totalNorm + norm(original.retrieveElement(j,k));
                }

        sigma = sqrt(sigma/totalNorm);
        return sigma;
}

parameterType QCGGS::iterate(QCBEC<2> &bec)
{
    //Forward Part of the Algorithm
    for(int k = 0; k < noOfImages-1; k ++)
    {
        updateMagnitude(Images(k),*ptrField); //Write magnitude to result field
        copyField(GSPhases(k),*ptrField); //Restore intermediate phase result
        for(int j = 0; j < maxTimeSteps; j ++)
            bec.evolve(true); //Evolve to next Image
    }
    rmsError = RMSError(Images(noOfImages-1),*ptrField);
    //Reverse Part of the Algorithm
    for(int k = noOfImages-1; k > 0; k --)
    {
        updateMagnitude(Images(k),*ptrField); //Write magnitude to result field
        copyField(GSPhases(k),*ptrField); //Restore intermediate phase result
        for(int j = 0; j < maxTimeSteps; j ++)
            bec.evolve(false); //Evolve to previous Image
    }
        iterationNo ++;
        return rmsError;
}

parameterType QCGGS::iterate(QC2ComponentBEC<2> &bec)
{
    //Forward Part of the Algorithm
    for(int k = 0; k < noOfImages-1; k ++)
    {
        updateMagnitude(Images(k),*ptrField); //Write magnitude to result field
        updateMagnitude(Images2(k),*ptrField2);
        copyField(GSPhases(k),*ptrField); //Restore intermediate phase result
        copyField(GSPhases2(k),*ptrField2);
        for(int j = 0; j < maxTimeSteps; j ++)
            bec.evolve(true); //Evolve to next Image
    }
    rmsError = RMSError(Images(noOfImages-1),*ptrField);
    //Reverse Part of the Algorithm
    for(int k = noOfImages-1; k > 0; k --)
    {
        updateMagnitude(Images(k),*ptrField); //Write magnitude to result field
        updateMagnitude(Images2(k),*ptrField2);
        copyField(GSPhases(k),*ptrField); //Restore intermediate phase result
        copyField(GSPhases2(k),*ptrField2);
        for(int j = 0; j < maxTimeSteps; j ++)
            bec.evolve(false); //Evolve to previous Image
    }
        iterationNo ++;
        return rmsError;
}

---