/*
 * solver_bndy_isFocusBndy.cpp
 *
 *  Created on: Feb 6, 2014
 *      Author: chuan
 */

#include "../misc/misc_interpl.h"
#include "solver.h"


bool needByteSwap(ifstream& inFileStream)
{
   int32_t v;
   inFileStream.read(reinterpret_cast<char*>(&v),sizeof(int32_t));
   inFileStream.seekg(0, inFileStream.beg);
   bool swap = (v < 0 || v >= 65536) ? true : false;
   return swap;
}

long getFileLength(ifstream& inFileStream)
{
    inFileStream.seekg (0, inFileStream.end);
    long length = inFileStream.tellg();
    inFileStream.seekg (0, inFileStream.beg);
    return length;
}


void swap4Bytes(char *in)
{
    char temp;
    temp = *(in + 3);
    *(in + 3) = *(in + 0);
    *(in + 0) = temp;
    temp = *(in + 2);
    *(in + 2) = *(in + 1);
    *(in + 1) = temp;
}

int32_t readInt32FromStream(ifstream& inFileStream, bool swap)
{
    char bytes[4];
    inFileStream.read(bytes, 4);
    if(swap) swap4Bytes(bytes);
    int32_t value = *reinterpret_cast<int32_t*>(bytes);
    return value;
}

float readFloatFromStream(ifstream& inFileStream, bool swap)
{
    char bytes[4];
    inFileStream.read(bytes, 4);
    if(swap) swap4Bytes(bytes);
    float value = *reinterpret_cast<float*>(bytes);
    return value;
}


bool CDelphiSolver::isFocusBndy(delphi_real *** phimap)
{
    int32_t recordhead, recordfoot;
    char charhead[4], charfoot[4], chartemp[4];
    char toblbl[20], label[10], title[60], botlbl[16];
    delphi_real fScale1, goff, gChildBoxCenter, goff1;
    SGrid<delphi_real> fgBoxCenter1, fgXYZ, fcXYZ, fgTemp;
    delphi_integer iGrid1, isgrid;
    int i, ix, iy, iz, iout, flag_s, flag_f; //flag_s: start; flag_f: finish.
    delphi_real *** phimap_pre;
    string strLine, strSubLine;
    SGrid < delphi_real > oldmid_1;

    size_t MAXWIDTH = 45;
#if !defined(MCCE) && !defined(PRIME)

    if (debug_solver) cout << "#### in isFocusBndy: phiintype: " << phiintype << endl;

    if (phiintype == 0) 
    {   // CODE to read old .phi (delphi/c++) not standard fortran compatible
        // Thus commented out, new code writes and reads fortran compatible
        // standard .phi format, to be consistent with implementation of .phi in
        // all other tools e.g. visualizers
        // SHAILESH: JAN 14, 2020
        /*
        ifstream ifPhiiFileStream;
        ifPhiiFileStream.open(strPhiiFile.c_str());
        if (!ifPhiiFileStream.is_open()) 
        {
            throw CUnknownPhiiFile(strPhiiFile);
            return false;
        }

        cout << "\n";
        cout << " ---------------- Focussing boundary condition ----------------- " << endl;
        cout << left << setw(MAXWIDTH) << " Read potential from file " << " : ";
        cout << strPhiiFile << endl << endl;

        ifPhiiFileStream.read(toblbl, 21);
        ifPhiiFileStream.read(label, 10);
        ifPhiiFileStream.read(title, 54);

        for (iz = 0; iz < iGrid; iz++) 
        {
            for (iy = 0; iy < iGrid; iy++) 
            {
                for (ix = 0; ix < iGrid; ix++) 
                {
                    ifPhiiFileStream.read(reinterpret_cast<char*>(&phimap[iz][iy][ix]), sizeof(delphi_real));
                }
            }
        }

        ifPhiiFileStream.read(botlbl, 16);

        ifPhiiFileStream.read(reinterpret_cast<char*>(&fScale1), sizeof(fScale1));
        ifPhiiFileStream.read(reinterpret_cast<char*>(&fgBoxCenter1), sizeof(fgBoxCenter1));
        ifPhiiFileStream.read(reinterpret_cast<char*>(&iGrid1), sizeof(iGrid1));

        ifPhiiFileStream.close();
        */

        ifstream ifPhiiFileStream;
        ifPhiiFileStream.open(strPhiiFile.c_str());
        if (!ifPhiiFileStream.is_open())
        {
            throw CUnknownPhiiFile(strPhiiFile);
            return false;
        }

        cout << "\n";
        cout << " ---------------- Focussing boundary condition ----------------- " << endl;
        cout << left << setw(MAXWIDTH) << " Read potential from file " << " : ";
        cout << strPhiiFile << endl << endl;
        long fileLength = getFileLength(ifPhiiFileStream);
        bool swap = needByteSwap(ifPhiiFileStream);
        ifPhiiFileStream.seekg(fileLength-8, ifPhiiFileStream.beg);
        iGrid1 = readInt32FromStream(ifPhiiFileStream, swap);
        ifPhiiFileStream.seekg(0, ifPhiiFileStream.beg);
        //cout << "fileLength=" << fileLength << "  iGrid1=" << iGrid1 << endl;

        recordhead = readInt32FromStream(ifPhiiFileStream, swap);
        ifPhiiFileStream.read(toblbl, 20);
        //cout << "swap=" << swap << " toplbl='" << toblbl << "' record-head=" << recordhead << endl;
        recordfoot = readInt32FromStream(ifPhiiFileStream, swap);
        if(recordhead != recordfoot)
        {
           stringstream ss;
           ss << "toplbl record head="<< recordhead<<" and foot="<<recordfoot<<" are unequal in " << strPhiiFile << endl;
           throw CInvalidRecordRead(ss.str());
        }
        recordhead = readInt32FromStream(ifPhiiFileStream, swap);
        ifPhiiFileStream.read(label, 10);
        ifPhiiFileStream.read(title, 60);
        recordfoot = readInt32FromStream(ifPhiiFileStream, swap);
        if(recordhead != recordfoot)
        {
           stringstream ss;
           ss << "toplbl record head="<< recordhead<<" and foot="<<recordfoot<<" are unequal in " << strPhiiFile << endl;
           throw CInvalidRecordRead(ss.str());
        }

        recordhead = readInt32FromStream(ifPhiiFileStream, swap);
        int32_t iGridRecord = (int32_t) cbrt(recordhead/4);
        if(iGridRecord != iGrid1)
        {
            stringstream ss;
            ss << "Phimap Grid data record gridsize= " << iGridRecord << " not equal to value shown in file=" << iGrid1 << endl;
            throw CInvalidRecordRead(ss.str());
        }
        //--------- allocate memory for previous phimap -------
        phimap_pre_v.assign(iGrid1 * iGrid1 * iGrid1, 0.0);
        phimap_pre = pdc->getKey_Ptr < delphi_real > ("phimap_pre", iGrid1, iGrid1, iGrid1); // pointer to 3D phimap
        for (iz = 0; iz < iGrid1; iz++)
        {
            for (iy = 0; iy < iGrid1; iy++)
            {
                for (ix = 0; ix < iGrid1; ix++)
                {
                    //ifPhiiFileStream.read(reinterpret_cast<char*>(&gridpointpotential), sizeof(gridpointpotential));
                    // phimap[iz][iy][ix] = readFloatFromStream(ifPhiiFileStream, swap);
                    phimap_pre[iz][iy][ix] = readFloatFromStream(ifPhiiFileStream, swap);
                }
            }
        }
        recordfoot = readInt32FromStream(ifPhiiFileStream, swap);
        if(recordhead != recordfoot)
        {
           stringstream ss;
           ss << "griddata record head="<< recordhead<<" and foot="<<recordfoot<<" are unequal in " << strPhiiFile << endl;
           throw CInvalidRecordRead(ss.str());
        }

        recordhead = readInt32FromStream(ifPhiiFileStream, swap);
        ifPhiiFileStream.read(botlbl, 16);
        recordfoot = readInt32FromStream(ifPhiiFileStream, swap);
        if(recordhead != recordfoot)
        {
           stringstream ss;
           ss << "botlbl record head="<< recordhead<<" and foot="<<recordfoot<<" are unequal in " << strPhiiFile << endl;
           throw CInvalidRecordRead(ss.str());
        }

        recordhead = readInt32FromStream(ifPhiiFileStream, swap);
        fScale1 = readFloatFromStream(ifPhiiFileStream, swap);
        fgBoxCenter1.nX = readFloatFromStream(ifPhiiFileStream, swap);
        fgBoxCenter1.nY = readFloatFromStream(ifPhiiFileStream, swap);
        fgBoxCenter1.nZ = readFloatFromStream(ifPhiiFileStream, swap);
        iGrid1 = readInt32FromStream(ifPhiiFileStream, swap);
        recordfoot = readInt32FromStream(ifPhiiFileStream, swap);
        if(recordhead != recordfoot)
        {
           stringstream ss;
           ss << "bottom boxparams record head="<< recordhead<<" and foot="<<recordfoot<<" are unequal in " << strPhiiFile << endl;
           throw CInvalidRecordRead(ss.str());
        }

        ifPhiiFileStream.close();

    } 
    else if (phiintype == 1) // cube format
    {
        ifstream ifPhiiFileStream;
        ifPhiiFileStream.open(strPhiiFile.c_str());
        if (!ifPhiiFileStream.is_open()) 
        {
            throw CUnknownPhiiFile(strPhiiFile);
            return false;
        }

        cout << "\n";
        cout << " ---------------- Focussing boundary condition ----------------- " << endl;
        cout << left << setw(MAXWIDTH) << " Read potential from file " << " : ";
        cout << strPhiiFile << endl << endl;


        // first line
        getline(ifPhiiFileStream, strLine);

        strSubLine  = strLine.substr(0, 10);
        fScale1     = atof(strSubLine.c_str());

        strSubLine  = strLine.substr(10, 6);
        iGrid1      = atoi(strSubLine.c_str());

        strSubLine  = strLine.substr(16, 10);
        oldmid_1.nX = atof(strSubLine.c_str());

        strSubLine  = strLine.substr(26, 10);
        oldmid_1.nY = atof(strSubLine.c_str());

        strSubLine  = strLine.substr(36, 10);
        oldmid_1.nZ = atof(strSubLine.c_str());

        fgBoxCenter1 = oldmid_1;

		/*==============================================================================
		 * THIS SEGMENT OF THE CODE WAS WRITTEN TO GET THE CENTER OF THE PARENT PHIMAP
		 * FROM THE USER VIA ACENTER. HOWEVER, THIS RESULTED INTO A NEW ISSUE WHERE THE
		 * PARENT POTENTIALS WERE ACCIDENTALLY TRANSLATED.
		 *==============================================================================
		 */
		/*
        // Warn if the center fo the parent run does not coincide with the box center
        // provided by the user for the focussing run. One provided by the user will be used.
        if ( fgBoxCenter1 != fgAcenter ) 
        {
            CAcenterMismatch warn_acenter_mismatch(fgBoxCenter1, fgAcenter);
            fgBoxCenter1 = fgAcenter;
        }
		*/
        
		cout << left << setw(MAXWIDTH) << " Scale of the parent run " << " : " << fScale1 << endl;
        cout << left << setw(MAXWIDTH) << " Grid size of the parent run" << " : "  << iGrid1 << endl;
        cout << left << setw(MAXWIDTH) << " Box center of the parent run" << " : " << fgBoxCenter1.nX  
             << " " << fgBoxCenter1.nY    << " " << fgBoxCenter1.nZ  << endl;
       

        cout << left << setw(MAXWIDTH) << " Scale of the child run" << " : " << fScale << endl;
        cout << left << setw(MAXWIDTH) << " Grid size of the child run" << " : " << iGrid  << endl; 
        cout << left << setw(MAXWIDTH) << " Box center of the child run" << " : " << fgBoxCenter.nX << " " 
             << fgBoxCenter.nY << " " << fgBoxCenter.nZ << endl;
        

        //--------- allocate memory for previous phimap -------
        phimap_pre_v.assign(iGrid1 * iGrid1 * iGrid1, 0.0);
        phimap_pre = pdc->getKey_Ptr < delphi_real > ("phimap_pre", iGrid1, iGrid1, iGrid1); // pointer to 3D phimap


        // skip the 2nd line like every cube file parser
        getline(ifPhiiFileStream, strLine);

        // 3rd line : fetch number of atoms and origin (corner)
        getline(ifPhiiFileStream, strLine);
        
        // num of atoms (fake but important)    
        strSubLine = strLine.substr(0,5);  
        int num_atoms = atoi(strSubLine.c_str());
        
        // 4th line : read iGrid1 and scale 
        // --------- REMOVE TO AVOID REDUNDANT READ/SKIPPING --------- 
		// getline(ifPhiiFileStream, strLine);

        // skip lines 4-7
        for (i = 1; i <= (3 + num_atoms); i++ )
            getline(ifPhiiFileStream, strLine);
        
       
        // start reading the phi values
        flag_s = 0;
        for (iz = 0; iz < iGrid1; iz++) 
        {
            for (iy = 0; iy < iGrid1; iy++) 
            {
                flag_s = 0;
                while (flag_s < iGrid1) 
                {
                    getline(ifPhiiFileStream, strLine);
                    if (ifPhiiFileStream.eof()) break;
                    
                    if (flag_s + 5 > iGrid1 - 1) 
                            flag_f = iGrid1 - 1;
                    else
                        flag_f = flag_s + 5;
                    
                    i = 0; //for ith number in a line
                    for (ix = flag_s; ix <= flag_f; ix++) 
                    {
                        strSubLine = strLine.substr(13 * i, 13);
                        i++; //for ith number in a line
                        //cout << "strSubLine: " << strSubLine << endl;
                        //phimap_pre_v[iz*iGrid1*iGrid1+iy*iGrid1+ix] = atof( strSubLine.c_str() );
                        phimap_pre[iz][iy][ix] = atof(strSubLine.c_str());
                        //cout << phimap_pre_v[iz*iGrid1*iGrid1+iy*iGrid1+ix] << " ";
                    }
                    flag_s = flag_s + 6;
                    //cout << endl;
                }
            }
        }

        ifPhiiFileStream.close();

    } 
    else 
    {
        cout << "Focusing is not CUBE or PHI" << endl;
        exit(0);
    }
#endif

#ifdef MCCE
    phimap_pre_v = pmcce->phimap;
    fScale1 = pmcce->scale1;
    fgBoxCenter1 = pmcce->oldmid1;
    iGrid1 = pmcce->igrid1;

    phimap_pre = pdc->getKey_Ptr<delphi_real>("phimap_pre",iGrid1,iGrid1,iGrid1);
#endif

#ifdef PRIME
    phimap_pre_v = pPrime->phimap;
    fScale1 = pPrime->scale1;
    fgBoxCenter1 = pPrime->oldmid1;
    iGrid1 = pPrime->igrid1;

    phimap_pre = pdc->getKey_Ptr<delphi_real>("phimap_pre",iGrid1,iGrid1,iGrid1);
#endif

    if (fZero > abs(fScale1 - fScale)) 
    {
#ifdef VERBOSE      
        cout << "scales are the same.\n";
        cout << "therefore assuming this to be a continuence\n";
#endif      
    } 
    else 
    {
#ifdef VERBOSE          
        cout << "\n";
        cout << " focussing potential map:\n";
        if (debug_solver) cout << title << endl;
        cout << "original scale (grids/A)      : " << fScale1 << endl;
        cout << "object centre at (A) : " << fgBoxCenter1.nX << " "
        << fgBoxCenter1.nY << " " << fgBoxCenter1.nZ << endl;
#endif  
        //----- check to see that new grid lies within old one that is going to provide bc's
        iout = 0;
        
		// child run
		goff = (iGrid + 1.0) / 2.0;
        
		// parent run
		goff1 = (iGrid1 + 1.0) / 2.0;

        for (iz = 0; iz < iGrid; iz += iGrid - 1) 
        {
            for (iy = 0; iy < iGrid; iy += iGrid - 1) 
            {
                for (ix = 0; ix < iGrid; ix += iGrid - 1) 
                {
                    fgXYZ.nX = (delphi_real)(ix + 1);
                    fgXYZ.nY = (delphi_real)(iy + 1);
                    fgXYZ.nZ = (delphi_real)(iz + 1);

                    //for each new grid corner, calculate old grid coords
                    fcXYZ = (fgXYZ - goff) / fScale + fgBoxCenter;
                    fgTemp = (fcXYZ - fgBoxCenter1) * fScale1 + goff1;

                    if (optORLE < delphi_real > (fgTemp, 1.0) || optORGE < delphi_real > (fgTemp, (delphi_real) iGrid1)) iout = 1;
                }
            }
        }

        if (0 != iout) throw COutsideFocus(fScale1, fgBoxCenter1, fScale, fgBoxCenter);

        /*
         * for each boundary point
         *    convert to delphi_real coordinates
         *    convert to old grid coordinates
         *    interpolate potential
         * note that can use same potential array for boundaries since old potentials at boundary are not used for new ones
         */

        //----- save new grid size, and set temporarily to 65
        isgrid = iGrid;               // isgrid is the grid size of the child run
        iGrid = iGrid1;               // iGrid is now the grid size of the parent run
        gChildBoxCenter = (isgrid + 1.0) / 2.0;  // midpoint of child run's box
#ifdef VERBOSE      
        cout << "pulling boundary values out of old potential map...\n";
#endif
        //cout << "phimap_pre[0]: " << phimap_pre[0][0][0] << endl;
        //cout << "phimap_pre[1][1][2]: " << phimap_pre[0][0][1] << endl;

        for (iz = 0; iz < isgrid; iz++) 
        {
            for (iy = 0; iy < isgrid; iy++) 
            {
                for (ix = 0; ix < isgrid; ix += isgrid - 1) 
                {
                    fgXYZ.nX = (delphi_real)(ix + 1);
                    fgXYZ.nY = (delphi_real)(iy + 1);
                    fgXYZ.nZ = (delphi_real)(iz + 1);
                    fcXYZ    = (fgXYZ - gChildBoxCenter) / fScale + fgBoxCenter;
                    fgTemp   = (fcXYZ - fgBoxCenter1) * fScale1 + goff1;
					
                    //for each new grid side, calculate old grid coords find potential
                    //phimap[iz][iy][ix] = interpl(iGrid,phimap_pre,fgTemp);
                    phimap[ix][iy][iz] = interpl(iGrid, phimap_pre, fgTemp);
                }
            }
        }

        for (iz = 0; iz < isgrid; iz++) 
        {
            for (iy = 0; iy < isgrid; iy += isgrid - 1) 
            {
                for (ix = 0; ix < isgrid; ix++) 
                {
                    fgXYZ.nX = (delphi_real)(ix + 1);
                    fgXYZ.nY = (delphi_real)(iy + 1);
                    fgXYZ.nZ = (delphi_real)(iz + 1);
                    fcXYZ    = (fgXYZ - gChildBoxCenter) / fScale + fgBoxCenter;
                    fgTemp   = (fcXYZ - fgBoxCenter1) * fScale1 + goff1;

                    //for each new grid side, calculate old grid coords find potential
                    //phimap[iz][iy][ix] = interpl(iGrid,phimap_pre,fgTemp);
                    phimap[ix][iy][iz] = interpl(iGrid, phimap_pre, fgTemp);
                }
            }
        }

        for (iz = 0; iz < isgrid; iz += isgrid - 1) 
        {
            for (iy = 0; iy < isgrid; iy++) 
            {
                for (ix = 0; ix < isgrid; ix++) 
                {
                    fgXYZ.nX = (delphi_real)(ix + 1);
                    fgXYZ.nY = (delphi_real)(iy + 1);
                    fgXYZ.nZ = (delphi_real)(iz + 1);
                    fcXYZ    = (fgXYZ - gChildBoxCenter) / fScale + fgBoxCenter;
                    fgTemp   = (fcXYZ - fgBoxCenter1) * fScale1 + goff1;

                    //for each new grid side, calculate old grid coords find potential
                    //phimap[iz][iy][ix] = interpl(iGrid,phimap_pre,fgTemp);
                    phimap[ix][iy][iz] = interpl(iGrid, phimap_pre, fgTemp);
                }
            }
        }

        iGrid = isgrid;
    }

    if (debug_solver) 
    {
        cout << "phimap[0][0][0]: " << phimap[0][0][0] << endl;
        cout << "phimap[1][1][1]: " << phimap[1][1][1] << endl;
        cout << "phimap[1][1][0]: " << phimap[1][1][0] << endl;
        cout << "phimap[1][0][0]: " << phimap[1][0][0] << endl;

        cout << "iGrid: " << iGrid << endl;
        cout << "phimap[iGrid-1][iGrid-1][iGrid-1]: " << phimap[iGrid - 1][iGrid - 1][iGrid - 1] << endl;
        cout << "phimap[1][1][1]: " << phimap[1][1][1] << endl;
        cout << "phimap[1][1][iGrid-1]: " << phimap[1][1][iGrid - 1] << endl;
        cout << "phimap[1][iGrid-1][iGrid-1]: " << phimap[1][iGrid - 1][iGrid - 1] << endl;
    }

    /** remove phimap_pre pointers*/
    for (int i = 0; i != iGrid1; ++i) 
    {
        delete[] phimap_pre[i];
    }
    delete[] phimap_pre;

    //exit(0);

    return true;
}
