// Lemur OLAP library (c) 2003 National Research Council of Canada by Daniel Lemire, and Owen Kaser
 /**
 *  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 (version 2). 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 "normalizationscommon.h"

#include "normalutil.h"
#include "normalizations.h"
#include "pearsonnormalization.h"
#include "perlc-interface.h"
#include "slicesortingnormalization.h"
#include "greedyfrequencysort.h"
#include "multinormalizations.h"
#include "matching.h"
#include "cubestatistics.h"
#include "pngimage.h"
#include "perfectnormalization.h"
/*
 * Just dump info about Binary file.
 */
void dump(char * BinaryFileName) {
    cout << "[dump] Parsing " << BinaryFileName << "."<<endl;
    OwenParser<int,int> parser(BinaryFileName);
    cout << "[dump] Header is : \"" << parser.getHeader() << "\". "<<endl;
    vector<int> shape = parser.getShape();
    cout << "[dump] Number of dimensions is "<< shape.size() << "."<<endl;
    cout << "[dump] Numbers of attribute values per dimension follow... "<< endl;
    cout << "[dump] ";
    for(vector<int>::iterator i = shape.begin(); i != shape.end(); ++i) cout << *i << " ";
    cout <<endl;
    cout << "[dump] Done dumping info about file "<< BinaryFileName << "."<<endl;
}

void writePNG(string jointpng, float_table_type joint, const float min_proba = 0.0f, const float max_proba = 1.0f) {
    const int height = joint.size();
    const int width = joint[0].size();
    PNGImage image1(jointpng, width, height);
    image1.openFileWrite();
    image1.initWriteStructs();
    image1.writeHeader();
    if(max_proba - min_proba <= 0) {
        cerr << "Max = Min. Aborting Image Output." << endl;
        return;
    }
    png_byte * RGB = new png_byte[3* width * height];
    try {
        const int pixel_size = 3 * sizeof(png_byte);
        for (int i=0; i<height; i++) {
            for (int j=0; j<width; j++) {
                assert(joint[i][j] >= 0.0f);
                assert(joint[i][j] <= 1.0f);
                int value = (int) floor( (max_proba -joint[i][j])/(max_proba - min_proba) *256 );
                if( value >= 256 ) value = 255;
                if( value < 0 ) value = 0;
                png_byte pngvalue = (png_byte) value;
                RGB[i * pixel_size * width + j * pixel_size + 0] = pngvalue;  // red
                RGB[i * pixel_size * width + j * pixel_size + 1] = pngvalue;    // green
                RGB[i * pixel_size * width + j * pixel_size + 2] = pngvalue;    // blue
            }
        }
        image1.writeImage(RGB);
    } catch(...) {
        cout << "We are not doing well." <<endl;
        delete[] RGB;
    }
    image1.writeEnd();
    cout << "Written PNG image to " << jointpng << "." << endl;
}

void correlation(DataCube<int,int64>& DC) {
    const int d = DC.getShape().size();
    float_table_type correlation_matrix(d,row_type(d,0.0f));
    float_table_type percentage_matrix(d,row_type(d,0.0f));
    for(int dimension1 = 0; dimension1 < d; ++dimension1) {
        for(int dimension2 = 0; dimension2 < dimension1; ++dimension2) {
            float_table_type joint =
                CubeStatistics<int,int64>::jointDistribution(DC,dimension1,dimension2);
            const float correlation = CubeStatistics<int,int64>::correlation(joint);
            //const float max_proba = CubeStatistics<int,int64>::max(joint);
            //const float min_proba = CubeStatistics<int,int64>::min(joint);
            float_table_type separable =
                CubeStatistics<int,int64>::separableJointDistribution(joint);
            //const float max_proba_sep = CubeStatistics<int,int64>::max(separable);
            //const float min_proba_sep = CubeStatistics<int,int64>::min(separable);
            //const float max_p = max(max_proba,max_proba_sep);
            //const float min_p = min(min_proba,min_proba_sep);
            //cout << " Correlation between dim = "<< dimension1 << " and " <<
            // dimension2 <<
            //	" is " << correlation << " ("<<pow(correlation,2)*100<< "%)" << endl;
            //
            // Now, I'm about to print this off into images. For now, images
            // name are hardcoded. Not hard to change.
            //
            //writePNG("joint.png",joint,min_p,max_p);
            //writePNG("separate.png",separable,min_p,max_p );
            //
            //cout << "Now go and use an image viewer to check these distributions out! "
            //<<endl;
            correlation_matrix[dimension1][dimension2] = correlation;
        }
    }
    for(int dimension1 = 0; dimension1 < d; ++dimension1)
        correlation_matrix[dimension1][dimension1] = 1.0f;
    for(int dimension1 = 0; dimension1 < d; ++dimension1)
        for(int dimension2 = dimension1 + 1; dimension2 < d; ++dimension2)
            correlation_matrix[dimension1][dimension2]
            = correlation_matrix[dimension2][dimension1] ;
    for(int dimension1 = 0; dimension1 < d; ++dimension1)
        for(int dimension2 = 0; dimension2 < d; ++dimension2)
            percentage_matrix[dimension1][dimension2]
            = pow(correlation_matrix[dimension1][dimension2],2) * 100;
    cout << "Correlation matrix" << endl;
    for(int dimension1 = 0; dimension1 < d; ++dimension1)  {
        for(int dimension2 = 0; dimension2 < d; ++dimension2) {
            cout << setw(7) << setprecision(2)
            <<  correlation_matrix[dimension2][dimension1]  << " " ;
        }
        cout << endl;
    }
    cout << "Percentages matrix" << endl;
    for(int dimension1 = 0; dimension1 < d; ++dimension1)  {
        for(int dimension2 = 0; dimension2 < d; ++dimension2) {
            cout << setw(7) ;
            if(percentage_matrix[dimension2][dimension1] < 1.0f)
                cout	<< setprecision(1)  ;
            else if (percentage_matrix[dimension2][dimension1] < 10.0f)
                cout	<< setprecision(2)  ;
            else cout	<< setprecision(3)  ;
            cout	<< percentage_matrix[dimension2][dimension1] ;
            cout	<< "% " ;
        }
        cout << endl;
    }
    //Recall that correlation may not be significant for small number of attribute
    //values.
}

void displayStatistics(DataCube<int,int64>&	DC) {
    const vector<int> shape = DC.getShape();
    const uint d = shape.size();
    FrequencySort<int,int64,GreaterComparator<pair<uint64,int> > > fs;
    int64 allocatedcells = -1;
    // this will dump the the percentage of allocated cells in each attribute values,
    // the point is to see to what extend the cube is "spread out" or concentrated
    for(uint dim = 0; dim < d ; ++dim) {
        cout << "[stats] dim = " << dim << " size = " << shape[dim] << " : \n";
        deque<pair<int64,int> > sortedfreqs =  fs.sortedFrequencyHistogram(DC, dim);
        // first compute sum
        int64 sum = 0;
        for(deque<pair<int64,int> >::const_iterator it = sortedfreqs.begin(); it != sortedfreqs.end(); ++it)
            sum += it->first;
        if(allocatedcells < 0) allocatedcells = (int64) sum;
        else assert((int64) allocatedcells == sum);// sanity
        float cumul = 0.0f;
        for(deque<pair<int64,int> >::const_iterator it = sortedfreqs.begin(); it != sortedfreqs.end(); ++it) {
            float percent =  it->first / (float) allocatedcells * 100;
            cout << "\t" << percent << "% ";
            cumul += percent;
            cout << "( cumul = " << cumul << "%)"<<endl;
        }
        assert(cumul > 99); assert (cumul < 101);
        cout << endl;
    }
    cout << "[stats] allocated cells = " << allocatedcells << " density = " <<
    allocatedcells/(float)DC.getVolume() << endl;

}

string tempFileName(const vector<int>& GroupBys) {
    stringstream strs;
    for(vector<int>::const_iterator i = GroupBys.begin(); i != GroupBys.end(); ++i) strs << *i << "_";
    string temp;
    strs >> temp;
    temp = "DCTemp_"+temp+".bin"; // we will use this for our data cube
    return temp;
}


int main(int argc, char * argv[]) {
    srand( (unsigned)time( NULL ) ); // whenever we need random data
    // TODO: use getopt of any such variant
    cout << " Normalization tool. (c) 2003 NRC/CNRC by Daniel Lemire, Owen Kaser " << endl;
    // we display usage information
    cout << " Usage: " << argv[0] << " BinaryFileName -groupby n1,n2,n3 " << endl;
    cout << " or  " << argv[0] << " BinaryFileName -dumpinfo " << endl;
    //
    if(argc < 2) return 0;
    // Initialize variables
    vector<int> GroupBys;
    char * BinaryFileName = NULL;
    char * normFileName = NULL;
    bool stat = false; // whether to display some statistics about the cube
    bool fast = false; // omit all but the slow ones, useful for debugging
    bool verbose = false;
    bool timing = false;
    bool toy = false;
    bool run_optimal = false;// whether to try to find the optimal normalization
    bool preScramble = false;
    enum{ NONE, UP_EVEN, UP_SQRT } roundupcube = UP_SQRT;
    int64 system_memory = numeric_limits<int64>::max();
    vector<vector<int> > ToyData;
    //
    // next follow our ugly parsing of the command line, to be replaced by getopt
    //
    if(argv[1][0] != '-')	BinaryFileName = argv[1];
    for(int k = 1; k < argc ; ++k) {
        if(strcmp(argv[k],"-groupby")==0) {
            if(k+1 == argc) {cerr << "Missing groupbys" << endl; return 8;}
            char * groupbys = argv[++k];
            string number;
            for( int i = 0; groupbys[i] != 0; ++i) {
                if(groupbys[i] == ',') {
                    GroupBys.push_back(atoi(number.c_str()));
                    number.clear();
                } else {
                    number += groupbys[i];
                }
            }
            GroupBys.push_back(atoi(number.c_str()));
        }
        else if (strcmp(argv[k],"-dumpinfo")==0){
            dump(BinaryFileName);
        }
        else if (strcmp(argv[k],"-normfile")==0) {
            normFileName = argv[++k];
        }
        else if (strcmp(argv[k],"-stat")==0) {
            stat = true;
        }
        else if (strcmp(argv[k],"-fast")==0) {
            fast = true;
        }
        else if (strcmp(argv[k],"-optimal")==0) {
            run_optimal  = true;
        } 
        else if (strcmp(argv[k],"-toy")==0) {
            if(k+1 == argc) {cerr << "Missing toy data" << endl; return 9;}
            toy = true;
            ToyData.push_back(vector<int>());
            // expect what to follow to have the format  1 2 / 3 4 // where // means "end"
            //  and to be a matrix...
            do {
                char * content = argv[++k];
                if(strcmp(content, "//") == 0)
                    break;
                if(strcmp(content, "/") == 0) {
                    ToyData.push_back(vector<int>()); continue;
                }
                if( strcmp(content,"?") == 0 ) {// random number
                    const float random_number =  rand() /(float) RAND_MAX;
                    if(random_number > 0.5f)
                        ToyData[ToyData.size()-1].push_back(1);
                    else
                        ToyData[ToyData.size()-1].push_back(0);
                    continue;
                }
                //cout << " must be a number = " << content << endl;
                ToyData[ToyData.size()-1].push_back(atoi(content));
            } while (k + 1 < argc );
            GroupBys.resize(0);
            GroupBys.push_back(0);
            GroupBys.push_back(1);
        }
        else if (strcmp(argv[k],"-noroundup")==0) {
            roundupcube = NONE;
        }
        else if (strcmp(argv[k],"-roundup_even")==0) {
            roundupcube = UP_EVEN;
        }
        else if (strcmp(argv[k],"-roundup_sqrt")==0) {
            roundupcube = UP_SQRT;
        }
        else if (strcmp(argv[k],"-memorylimit")==0) {  // system memory, in kilobytes (suffices ignored)
            ++k;
            if(k == argc) {cerr << "Missing memory limit" << endl; return 10;}
            system_memory = atoll(argv[k]);
        }
        else if (strcmp(argv[k],"-groupbits")==0) {  // packed binary groupby info
            ++k;
            if(k == argc) { cerr << "Missing groupbits argument" << endl; return 11; }
            if (GroupBys.size() != 0) {
                cerr << "[Warning]: have both -groupby and -groupbits" << endl;
                GroupBys.resize(0);
            }
            unsigned long gbits = atoll(argv[k]);
            bitset<32> bv(gbits);
            for (int pos=0; pos < 32; pos++)
                if (bv.test(pos)) GroupBys.push_back(pos);
        }
        else if (strcmp(argv[k],"-prescramble")==0) {
            preScramble = true;
        }
        else if (argv[k][0] == 'v') {
            cout << "Option not recognized = "<< argv[k] << endl;
        }
    }
    if(toy) cout << "toy mode" << endl;
    if(!toy && (GroupBys.size() == 0)) {
        cout << "No groupby specified. Exiting. "<< endl;
        return 1;
    }
    if(!toy && (BinaryFileName == NULL)) {
        cout << "No file specified. Exiting. " << endl;
        return 2;
    }
    //
    // end of the parsing
    //
    if(verbose) {
        cout << "[info] Group by : ";
        for(vector<int>::iterator i = GroupBys.begin(); i != GroupBys.end(); ++i) {
            cout << *i << " ";
        }
        cout <<endl;
    }
    // next we contruct a temporary file name
    string temp = tempFileName(GroupBys);
    vector<int> DCShape(GroupBys.size());
    // parsing the data
    OwenParser<int,int64> parser(BinaryFileName);
    if(verbose) cout << "[info] Header is : \"" << parser.getHeader() << "\". "<<endl;
    vector<int> shape = parser.getShape();
    if(toy) {
        shape.resize(2);
        shape[0] = ToyData.size();
        shape[1] = ToyData[0].size();
    }
    cout << "[info] Your data cube will have dimensions... "<<endl;
    cout << "[info]  ";
    for(uint i = 0; i < GroupBys.size(); ++i) {
        DCShape[i] = shape[GroupBys[i]];
        cout << DCShape[i]<< " ";
    }
    cout << endl;
    // done parsing the data
    clock_t start, finish;
    double NombreDeSecondes = 0.0;// that's NumberOfSeconds in French to test your bilinguism -DL

    // some useful chunkings
    vector<int> chunkedByTwos(DCShape.size(),2);
    vector<int> chunkedByFours(DCShape.size(),4);
    vector<int> chunkedBySqrt = ChunkedDataCube<int,int64>::sqrtChunkShapes(DCShape);
    vector<int> volumeTwoChunking(DCShape.size(),1); volumeTwoChunking[0]=2;
    if (toy) { // hack, clean up
        chunkedByTwos = vector<int>(2,2);
        chunkedByFours = vector<int>(2,4);
        volumeTwoChunking = vector<int>(2,1);
        volumeTwoChunking[0]=2;
        chunkedBySqrt = vector<int>(2,3);  // probably a lie
    }


    //
    // First I have to option of rounding up the dimensions of the cube.
    // I got worried, rightly so I think, that too many of the chunks have
    // relatively small size. That's because our cube is so small. It tends
    // to throw off everything.
    //
    if(!toy) {
        switch (roundupcube) {
        case UP_EVEN:
            cout <<"[info] rounding up cube so that all dimensions are divisible by 2!" <<endl;
            cout <<"[info]  ";
            for(uint i = 0; i < GroupBys.size(); ++i) {
                DCShape[i] += (DCShape[i] % 2);  // definitely safe if DCShape[i] is positive
                cout << DCShape[i]<< " ";
            }
            cout <<endl;
            break;

        case UP_SQRT:
            // MUST be done AFTER chunkedBySqrt calculated, unless you want
            // to contemplate their interaction...this is not nice.

            cout <<"[info] rounding up cube so that all dimensions are divisible by sqrt-based chunking" <<endl;
            cout <<"[info]  ";
            for(uint i = 0; i < GroupBys.size(); ++i) {
                int chunk_span = chunkedBySqrt[i];
                DCShape[i] = (DCShape[i] + chunk_span - 1) / chunk_span * chunk_span;
                cout << DCShape[i]<< " ";
            }
            cout <<endl;
            break;

        case NONE:
        default:
            break;
        }
    }
    if(toy) {
        // replace the data cube by a toy data cube
        vector<int> smallDCShape(2);
        smallDCShape[0] = ToyData.size();
        smallDCShape[1] = ToyData[0].size();
        DCShape = smallDCShape;
    }
    RAMCube<int,int64>	DC (DCShape);
    if(toy) {
        DC.open();
        DC.fillWithZeroes(); // no longer needed thanks to Sean's patch!
        for(uint row = 0; row < ToyData.size(); ++ row) {
            for(uint col = 0; col < ToyData[row].size(); ++col) {
                assert(ToyData[row].size() == (uint) DCShape[1]);
                DC.put(ToyData[row][col],row,col);
            }
        }
    } else {
        if (DC.getApproxDenseStorageSize() > 0.25 * system_memory * 1024) {
            cout << "Cube storage requires  " <<
            (DC.getApproxDenseStorageSize()/(1024*1024)) <<
            "MB and this exceeds specified limit.  Exiting." << endl;
            return 12;
        }

        if (preScramble) {
            RAMCube<int,int64> tempDC(DCShape);
            tempDC.open();
            //tempDC.fillWithZeroes(); // no longer needed thanks to Sean's patch
            parser.fill(tempDC,GroupBys);
            vector<vector<int> > scrambler = PermutationUtil::randomPermutation(DCShape);
            DC.open();
            //DC.fillWithZeroes(); // no longer needed thanks to Sean's patch!
            NormalUtil<int,int64>::copyTo(tempDC, DC,  scrambler);
            tempDC.close();
        }
        else {
            DC.open();
            //DC.fillWithZeroes(); // no longer needed
            parser.fill(DC,GroupBys);
        }
    }
    NombreDeSecondes =  (double)(finish - start) / CLOCKS_PER_SEC;
    if(verbose) cout << "[info] done building data cube." <<endl;
    if(timing) cout << "[time] It took " << NombreDeSecondes << " s."<<endl;
    // Data cube should be alright now.
    //
    if(stat) {
      displayStatistics(DC);
    }
    const bool corr = true;
    if(corr) {	// this is the attempt at checking whether the distribution is separable
      correlation(DC);
    }

    //
    // Ok, enough fooling around!
    //
    
    /*
     * First we define the normalizations we will use
     * If you add a new normalization, put it here.
     */
    deque<Normalization<int, int64> *> n;
    n.push_back( new Normalization<int,int64> ());
    n.push_back( new SliceSort<int,int64>());// leave it in position 1
    if(run_optimal) {// expensive!!!
      n.push_back( new PerfectNormalization<int,int64> (chunkedByTwos)); 
    }

    //    n.push_back( new GreedyFrequencySort<int,int64> ());
    //
    // MatchingNorm appears to be broken. 
    //
    // Normalization Weighted Matching iterated over all dimensions
    //  normalizer: graphs.cpp:119: 
    //  std::vector<int, std::allocator<int> > 
    //  GraphForMatching::readRothbergOutput(const char*) const: Assertion `adam && eve' failed.
                //  
                // Owen: hmmm.  I put it back and ran a few tests that it survived.  It is important to
                // have rounding to an even number of values.  If you discover a test case, let me know
                // how to replicate the bug...  if /tmp/Rothberg.in isn't too big, send it to me, also.
    //
    if(!toy && roundupcube == UP_EVEN) n.push_back( new MatchingNorm<int,int64> (MatchingNorm<int,int64>::ALL_DIMS) );
    //
    //
    //
    //    n.push_back( (MatchingNorm<int,int64> (1) * MatchingNorm<int,int64> ( 0 )).clone());
    //    n.push_back( (MatchingNorm<int,int64> (2) * MatchingNorm<int,int64> (1) * MatchingNorm<int,int64> (0)).clone());

    // n.push_back( new GreedyIterSort<int,int64> ());
    if(!fast) {
        if(normFileName != 0)
            n.push_back( new NormalizationReader<int,int64> (normFileName));
        n.push_back(new IteratedSliceCluster<int,int64>(chunkedByTwos));
        n.push_back(new IteratedSliceCluster<int,int64>(chunkedBySqrt));
        if (*min_element(DCShape.begin(),DCShape.end()) >= 4)
            n.push_back(new IteratedSliceCluster<int,int64>(chunkedByFours));
        n.push_back(( IteratedSliceCluster<int,int64> (chunkedByTwos) *
                      FrequencySort<int,int64,LessComparator<pair<uint64,int> >  >()).clone());
    }
#if 0
    n.push_back(new RandomNormalization<int, int64>() );
#endif

    /*
     *  next we actually compute the normalizations
     *
     *  This takes a while.
     */
    vector<vector<int> > normal[n.size()];
    for (uint i = 0; i < n.size(); ++i) {
        cout << "Normalization " + n[i]->getTextName() << endl;
        start = clock();
        normal[i] = n[i]->computeNormal(DC);
        finish = clock();
        NombreDeSecondes =  (double)(finish - start) / CLOCKS_PER_SEC;
        if(timing) cout << "[time] It took " << NombreDeSecondes << " s."<<endl;
        NormalUtil<int,int64>::printSmall2d(DC,normal[i],chunkedBySqrt);
    }
    /*
     *  Ok, The normalizations have been computed.
     *  
     *  last, we display the results
     */
    int maxm = min(8,*min_element(DCShape.begin(),DCShape.end()));
    int maxmPlusSpecials = maxm+2;  // forgive me!
    for(int m = 2; m <= maxmPlusSpecials ; ++m) {
        vector<int> ChunkShape;
        if (m > maxm) {  // yuck
            cout << "----------------------" << endl;
            // process irregular chunks here...

            switch (m - maxm) {
            case 1:
                ChunkShape = chunkedBySqrt;
                cout << "Using sqrt-based chunking: cube shape (in chunks) approximates chunk shape (in cells)"
                << endl;
                break;
            case 2:
                ChunkShape = volumeTwoChunking;
                cout << "Using 2x1x1...x1 chunking" << endl;
                break;
            }
        }
        else {
            ChunkShape = vector<int>(DC.getShape().size(), m);
            cout << "------------------------" << endl;
            cout << "Using regular chunks of size " << m << endl;
        }
        vector<pair<int, string> > scores;
        for (uint i = 0; i < n.size(); ++i) {
            // this is really ugly code, but I want nicer formatting of the
            // output so I can read it even when I'm tired.
            if(verbose) cout << "Computing cost for "<<n[i]->getTextName() << endl;
            const int cost = HOLAPUtil<int,int64>::cost(DC,ChunkShape,normal[i],verbose);
            scores.push_back(pair<int,string>(cost,n[i]->getTextName()) );
        }
        sort(scores.begin(), scores.end());
        pair<int,string> last = scores[scores.size()-1];
        cout << endl;
        for(uint algo = 0; algo < scores.size(); ++algo) {
            cout << setw(3) << (algo+1)<< "=" <<
            setw(6) << setprecision(3) << scores[algo].first << " ("<<
            setw(6) << setprecision(3) <<
            (scores[algo].first /(float) last.first * 100) <<"%)" <<
            "- "<<scores[algo].second << endl;
        }
        cout << endl;
    }
    /*
     *
     * We display the independence product. It is a bound on how
     * bad FrequencySort can be.
     */
    const pair<float,float> IndependenceInfo = CubeStatistics<int,int64>::independenceInfo(DC);   
    cout << "Bound on frequency sort lack of optimality = " << setw(7) << setprecision(7)  
      << IndependenceInfo.first << " (useless measure?) " << endl;
    cout << " Independence sum = "   << setw(7) << setprecision(7) << IndependenceInfo.second
    << " (1 = independent, 0 = dependent) "<< endl;
    /*
     *
     * That's it for the independence product.
     */
    bool dualnormalization = false;// use it if you want... I'm turning it off for now
    if(dualnormalization) {
        DualNormalizationScheme<int,int64> dns;
        DualNormalization dn = dns.computeNormal(DC);
        int i = 1; //position of FrequencySort or equivalent
        for(int m = 2; m <= maxm ; ++m) {
            vector<int> ChunkShape(DC.getShape().size(), m) ;
            cout << "m = "<< m << endl;
            float fs = HOLAPUtil<int,int64>::cost(DC,ChunkShape,normal[i],false);
            float dfs = DualHOLAPUtil<int,int64>::cost(DC,ChunkShape,dn,false);
            cout << " Frequency Sort = " <<  (int) fs
            << " (" << (fs / max(fs,dfs) * 100)<< "%)" <<  endl;
            cout << " Dual Frequency Sort = " << (int)  dfs
            << " (" << (dfs / max(fs,dfs) * 100)<< "%)" <<   endl;
        }
    }
    //
    // We must clean up!
    //
    DC.close();
    //
    for (uint i = 0; i < n.size(); ++i) delete n[i];
}


