/*-----------------------------------------------------------*
 *
 * File : STD.c
 * Author :  NTM
 * Created : 16/07/06
 *
 *
 * Copyright (C) Nicolas Thierry-Mieg, 2006.
 *
 *
 * This file is part of InterPool, written by 
 * Nicolas Thierry-Mieg (CNRS, France) Nicolas.Thierry-Mieg@imag.fr
 *
 * InterPool 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.
 *
 * InterPool 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 InterPool; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 *
 *------------------------------------------------------------*/

#include <stdio.h> // printf
#include <stdlib.h> // exit
#include <math.h> // log

#include <string.h> /* memset */
#include "types.h" /* MOT */
#include "fonc.h" /* firstDivisor */
#include "pools.h" /* pool, getpoolvalue, setpoolvalue */
#include "STD.h"


/************************************************************************
 ******************* LOCAL FUNCTIONS ************************************
 ************************************************************************/

/////////////// DECLARATIONS ///////////////


//build one layer of STD pools, and place it (where it belongs) in tabpool
static void buildlayer(MOT* tabpool, int n, int q, int layer);


//save the set of pools tabpool in file named designFile
// format for the file is identical to buildPools in pools.h.
static void writeDesign(char* designFile, MOT* tabpool, int n, int nbPools);


//////////////// BODIES ///////////////////


///////////////////////////////////////////////////////////////
// buildlayer: builds one layer of STD pools, and places it 
// (where it belongs) in tabpool.
// In practice we only set to 11 the bits corresponding to variables
// present in each pool.
// Therefore we assume that tabpool was initialized to 0's
static void buildlayer(MOT* tabpool, int n, int q, int layer)
{
  /* compression of q rel to n */
  int comp = compression(q,n) ;
  
  //make sure layer<q 
  if (layer>=q)
    {
      fprintf(stderr,"in buildlayere, layer >= q! dying\n");
      exit(1);
    }

  // for each variable, place it in the correct pool of this layer
  {
    int variable;
    for (variable=0;variable<n;variable++)
      { 
	//pool number where variable should go
	int poolnb=0;
     
	int c;
	for(c=0;c<=comp;c++)
	  {
	    //calculate poolnb (according to STD definition)
	    poolnb=poolnb+(int)(pow(layer,c))*(int)(variable/(pow(q,c)));
	  }
	poolnb=poolnb % q;
	
	// find the memory area for this pool
	{
	  int absPoolNum = layer*q + poolnb ;
	  MOT* poolForVar=pool(tabpool,absPoolNum,n);
	
	  // set the bits coding for variable in pool to 11
	  setpoolvalue(poolForVar,variable,11);
	}
      }
  }
}

///////////////////////////////////////////////////////////////
// writeDesign:
// designFile is a string, it holds the filename (including
// path) where the pool from tabpool must be saved.
// format in the file is compatible with readDesign.
static void writeDesign(char* designFile, MOT* tabpool, int n, int nbPools)
{
  FILE* file;

  /* make sure the file doesn't already exist */
  file = fopen(designFile, "r");
  if (file != NULL)
    {
      fprintf(stderr,"In writeDesign: cannot write to %s, file already exists\n",designFile);
      exit(1);
    }
  
  /* open file for writing */
  file = fopen(designFile, "w");
  if (file == NULL)
    {
      fprintf(stderr,"In writeDesign: cannot open %s for writing\n",designFile);
      exit(1);
    }

  {
    int absPoolNum ;
    for (absPoolNum=0; absPoolNum < nbPools; absPoolNum++)
      {
	MOT* thispool = pool(tabpool, absPoolNum, n) ;
	
	/* for each var, see if it in thispool */
	int var ;
	for(var=0; var<n; var++)
	  if(getpoolvalue(thispool, var) == 11)
	    fprintf(file,"%d:",var) ;
	/* add newline */
	fprintf(file, "\n") ;
      }
  }
  fclose(file);
}


/************************************************************************
 ******************* EXPORTED FUNCTIONS *********************************
 ************************************************************************/

/////////////////////////////////////////////////////////////////////
// compression: calculate the compression value of q relative to n, ie
// the smallest c such that q^(c+1) >= n.
int compression(int q, int n)
{
  /* sanity check (make sure caller didn't swap params) */
  if (q > n)
    {
      fprintf(stderr, 
	      "in compression: q>n (values: %d, %d). Did you swap the params?\n",
	      q, n) ;
      exit(1) ;
    }
  int c = 0 ;
  int qpowc = q ;
  while (qpowc < n)
    {
      qpowc *= q ;
      c++ ;
    }
  return c ;
}


///////////////////////////////////////////////////////////////
// optimalq:
// given n and nbPosVars, returns the "optimal" value to use
// for q when using STD.
//
// Choice of q:
// For best results, we want approx. 50% of positive exps. The following holds:
// P[pool containing F fishes is positive] = 1 - (1-fracposvars)**F
// Therefore 50% positive exps is equiv to:
//   F = -1/(log2(1-fracposvars))   (see the demo in code3BestTau, in the 
//                                   perl simulator)
// However, here we use STD to generate pools, in which 
//   F = int(n/q) +- 1, which can be written approximately as F=n/q.
// Therefore using STD for fishes, 50% exps are positive if q = -n*(log2(1-fracposvars)).
int optimalq(int n, int nbPosVars)
{
  double fracposvars = (double) nbPosVars / n ; // since nbPosVars is cast to double, n is autocast
  double mylog = log(1-fracposvars) / M_LN2 ; /* log2(X) = log(X)/log(2) , M_LN2 is in math.h */
  double bestq = -1 * n * mylog ;
  return ((int)(bestq + 0.5)) ;
  // add 0.5 and cast to int returns the nearest int
}


/*!
  \brief return the average pool size when using STD.
*/
int poolSizeSTD(int n, int nbPools)
{
  return(n / calculateQ(nbPools)) ;
}


/*!
  \brief calculate q (STD parameter): it is the largest prime
  factor of nbPools.
*/
int calculateQ(int nbPools)
{
  int tmpNP = nbPools ;
  int divisor ;
  while ((divisor = firstDivisor(tmpNP)) != 0)
    tmpNP /= divisor ;
  /* now tmpNP is the largest prime factor of nbPools: this is q */
  return (tmpNP) ;
}


/*!
  \brief calculate k (STD parameter).
*/
int calculateK(int nbPools)
{
  int q = calculateQ(nbPools) ;
  /* k == nbPools / q */
  int k = nbPools / q ;
  return(k) ;
}



//////////////////////////////////////////////////////////////////
// STD: build the pools STD(n;q;k) and save them to a
// file in the interPool design format.
// File will be named STD.n<n>.q<q>.k<k> and placed in DESIGN_DIR.
void STD(int n,int q,int k)
{

  //------------------------------------------------
  // make sure the constraints for STD are respected
  //------------------------------------------------
  if (q >= n)
    {
      fprintf(stderr, "In STD, q >= n, error\n");
      exit(1);
    }
  
  else if (k > q)
    {
      fprintf(stderr, "In STD, error k>q\n");
      exit (1);
    }
  
  else if (firstDivisor(q) != 0)
    {
      fprintf(stderr, "In STD, error: q is not prime\n");
      exit(1);
    }
  

  //--------------------------------------------------
  // build the set of pools.
  //--------------------------------------------------

  {
    /* filename: string for storing the complete filename (including path) */
    char filename[DESIGNFILELENGTH];
    int numchars ;
    FILE* file ;



    numchars = snprintf(filename,DESIGNFILELENGTH,"%s/STD.n%d.q%d.k%d",DESIGN_DIR,n,q,k);
    // make sure file name was not too long
    if ((numchars < 0) || (numchars >= DESIGNFILELENGTH))
      {
	fprintf(stderr, "in STD, error building filename: too long?\n");
	exit(1);
      }
    
    // try to open file
    file = fopen(filename, "r");
    
    if(file != NULL)
      /* file already exists, do nothing */
      fclose(file) ;
    else
      {
	// file doesn't exist, create the pools and save them to file
	int nbPools = q*k ;

	/* allocate the necessary memory for tabpool */
	int tabpoolsize = nbmotparpool(n)*nbPools ; // number of MOTs in tabpool

	MOT* tabpool=(MOT *)malloc(tabpoolsize * sizeof(MOT));
	if(tabpool==NULL)
	  {
	    fprintf(stderr,"in STD, not enough memory for tabpool\n");
	    exit(1);
	  }
	// initialize tabpool with 0's (we will later only change the required values to 11)
	memset(tabpool, 0, tabpoolsize * sizeof(MOT));

	int layer ;
	for(layer=0; layer < k; layer++)
	  buildlayer(tabpool,n,q,layer);

	// save the pools
	writeDesign(filename,tabpool,n,nbPools);

	free(tabpool) ;
      }
  }
}



///////////////////////////////////////////////////////////////
// buildqs: 
// returns a sorted array (allocated here) of at most nbofqs values for
// the parameter q of STD.
// these values are all prime numbers, >= qmin and <= qmax; or ZERO.
// They are selected so that the interval centered on qopt 
// and which contains them all is as small as possible.
// If the interval has some values smaller than qmin, it is
// shifted towards larger q values.
// Similarly, if it has some values larger than qmax it is shifted left.
// Finally if {qmin,..,qmax} doesn't contain at least nbofqs different primes,
// then the returned array contains all the possible primes preceded by zeros.
int* buildqs(int qopt, int nbofqs, int qmin, int qmax)
{
  /* allocate memory */
  int* allqs = (int*)malloc(nbofqs * sizeof(int)) ;
  if (allqs==NULL)
    {
      fprintf(stderr, "in buildqs, no more memory for allqs\n") ;
      exit(1) ;
    }
  /* initialize with zeros */
  memset(allqs, 0, nbofqs * sizeof(int)) ;

  /* sanity check: qmin must be positive */
  if (qmin < 0)
    {
      fprintf(stderr, "in buildqs, qmin is negative (value %d)\n", qmin) ;
      exit(1) ;
    }
  /* ... and smaller than qmax */
  if (qmin > qmax)
    {
      fprintf(stderr, "in buildqs, qmin > qmax (values %d %d)\n", qmin, qmax) ;
      exit(1) ;
    }

  int index = 0 ; // number of q values already saved

  // we want qopt to be odd, increment if it's even
  if ((qopt%2)==0)
    qopt++ ;

  /* keep qopt if it's prime and >= qmin  and <= qmax */
  /* use 2 ifs instead of &&, to make sure the arg to firstDivisor is >0 */
  if ((qopt >= qmin) && (qopt <= qmax))
    if (firstDivisor(qopt)==0)
      allqs[index++] = qopt ;

  // now start filling allqs with other primes
  int smallq = qopt ; // will go down from qopt
  int bigq = qopt ; // will go up from qopt

  while(index < nbofqs)
    {
      if ((bigq >= qmax) && (smallq <= qmin))
	{
	  /* can't generate anymore usable primes, remaining slots will stay at 0 */
	  break ;
	}

      bigq += 2 ; // +2 because we want to stay odd
      if ((bigq >= qmin) && (bigq <= qmax))
	if (firstDivisor(bigq)==0)
	  {
	    allqs[index++] = bigq ;
            if (index == nbofqs) // this was the last needed q
	      break ;
	  }
      

      smallq -= 2 ;
      if ((smallq >= qmin) && (smallq <= qmax))
	if (firstDivisor(smallq)==0)
	  allqs[index++] = smallq ;
    }

  // now sort the list, and return it
  qsort((void*)allqs, nbofqs, sizeof(int), compInts) ;
  return(allqs) ;
}


