/*****************************************************************
 * author Daniel Lemire 
 * (c) NRC 
 ************************************************************
 * Created March 23rd 2003: Convenient routine for toy experiments on
 *								efficient storage for data cubes. Work done by D. Lemire
 *								and O. 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.
 */

using namespace std;

#include <cmath>
#include <cassert>
#include <deque>
#include <vector>
#include <algorithm>
#include <map>
#include <iostream>

// it appears that false and true are defined in C++, at least under gcc
//bool false = 0;
//bool true = 1;

/*
* Check to see if value appears in the list permutation. There may be a built-in way
* to do this in STL. In fact, there is, it is "find", but for some reason I chose not
* to use it. A mistake?
*
* Also, this method is a double negative as it returns false if the value is found.
*/
bool doesntContain(deque<short>& permutation, short value) {
  for( deque<short>::iterator iter = permutation.begin(); iter != permutation.end() ; iter++) 
    if(*iter == value) return false;// ok, it is bad form to have a double negative, but allow me this one
  return true;
}

/* 
* For high performance purposes where memory is not an issue, it is better to compute the permutations
* once and for all, store the object and reuse it at will.
* Now, I should really discard some of these permutations based on the fact that I use block encoded
* data cubes, but then, this call would depend on m and I'm afraid it would bring added complexity
*/
void permutations(const short maximum, deque<short> & CurrentPermutation, deque<deque<short> > & PermutationSet) {
  for(short i = 0 ; i < maximum ; ++i) 
    if( doesntContain(CurrentPermutation, i) ) {
      deque<short> NewPermutation(CurrentPermutation);
      NewPermutation.push_back(i);
      if(NewPermutation.size() == maximum) 
        PermutationSet.push_back(NewPermutation);
      else
        permutations(maximum, NewPermutation, PermutationSet);      
    }
}

/*
* a bit ugly, but it appears to run quite fast.
*/
inline void normalize(const vector<vector<bool> >& array, vector<vector<bool> >& buffer, deque<short> & rows,  deque<short> & cols) {
  const int N = buffer.size();
  int r = 0, c = 0;
  for( deque<short>::iterator row = rows.begin(); /*row != rows.end()*/ r< N; ++row, ++r) {
    const vector<bool> & CurrentRowInArray = array[r];
    vector<bool> & CurrentRowInBuffer = buffer[*row];
    c = 0;
    for( deque<short>::iterator col = cols.begin(); /*col != cols.end()*/ c < N; ++col, ++c) {
        CurrentRowInBuffer[*col] = CurrentRowInArray[c];
    }
  }  
}



/*
* I don't think that Shannon entropy was the greatest of all ideas. I still leave this
* stuff here for historians to worry about when this work will win the nobel prize.
*
* Compiler will probably refuse to take this inline, try anyhow.
*/
inline float entropy(const vector<vector<bool> >& buffer,  int m) {
  const int N = buffer.size();
  int NumberOfSymbols = 1 << (m*m);
  int Samples = (N/m) * (N/m);
  map<int, float> p;
  for (int r = 0; r < N; r+= m) {
    for (int c = 0; c < N; c+= m) {
      int Symbol = 0;
      for(int i = 0; i < m ; ++i) {
        const vector<bool>& row = buffer[r+i];
        for(int j = 0; j < m ; ++j) {
          if(row[c+j])  Symbol += 1 << (i*m + j);
        }
      }
      map<int,float>::iterator iter = p.find(Symbol);
      if(iter == p.end()) 
        p.insert(pair<int,float>(Symbol, 1.0f / Samples));
      else 
        iter->second += 1.0f / Samples; 
    }
  }
  float entropy = 0.0f;
  for(map<int,float>::iterator iter = p.begin(); iter != p.end(); ++iter) 
    entropy -= iter->second * log(iter->second);
  entropy /= log((double)Samples);//log(NumberOfSymbols); here we normalized entropy
  return entropy;
}


/* This measure the efficiency as I defined in the paper. The effeciency is
* provably between 1 and 1/density. 1 being very good efficiency and 1/density
* being very poor efficiency.
*
* Compiler will probably refuse to take this inline, try anyhow.
*/
inline float holap(const vector<vector<bool> >& buffer,  int m) {
  const int N = buffer.size();
  int NumberOfSymbols = 1 << (m*m);
  int Samples = (N/m) * (N/m);
  float HOLAP = 0.0f;
  int TotalNonZero = 0;
  const int dimension = 2;
  for (int r = 0; r < N; r+= m) {
    for (int c = 0; c < N; c+= m) {
      int NonZero = 0;
      for(int i = 0; i < m ; ++i) {
        const vector<bool>& row = buffer[r+i];
        for(int j = 0; j < m ; ++j) {
          if(row[c+j]) ++NonZero; 
        }
      }
      TotalNonZero += NonZero;
      HOLAP += NonZero * dimension > m*m ?  m*m  :  dimension * NonZero ;
    }
  }
  if(TotalNonZero != N*N)
    return HOLAP  / TotalNonZero;// this is the definition I gave in the paper
  return 1.0f;
}


void print(const vector<vector<bool> > & array) {
  for(vector<vector<bool> >::const_iterator row = array.begin(); row != array.end(); ++row) { 
    for(vector<bool>::const_iterator value = row->begin(); value != row->end(); ++value) {
      cout << *value << " " ;
    }
    cout << endl;
  }
}


vector<vector<bool> > randomMatrix(const short N) {
  vector<vector<bool> > answer(N,vector<bool>(N,false));
  for(int i = 0 ; i < N * N / 2; ++i) {// the matrix is half and half
    //cout << "i = "<<i<<endl;
    while (true) {// this is stupid and inefficient
      int position =  (int)( ( (double) rand() * N * N ) / RAND_MAX );
      //cout << "Position = " << position<< " value = "<< answer[position/N][position%N]  << endl;
      if	(answer[position/N][position%N] == false) {
        answer[position/N][position%N] = true;
        break;
      }
    }
  }
  return answer;
}

inline void exchange(vector<vector<bool> >& x, const short Pos1, const short Pos2, const short Size) {
  for(int i = 0; i  < Size ; ++i) {
    vector<bool> old = x[i+Pos1];
    x[i+Pos1] = x[i+Pos2];
    x[i+Pos2] = old;// better hope that things don't overlap
  }
}

inline void unexchange(vector<vector<bool> >& x, const short Pos1, const short Pos2, const short Size) {
  exchange(x,Pos2,Pos1,Size);
}


