/* truecluster/src/combinatorics.c
 * Some combinatorial functions
 * (C) 2000-2007 Dr. Jens Oehlschlägel
 * R_tc_combination_noreplace
 * R_tc_combination_replace
 * R_tc_permutation_noreplace
 * R_tc_permutation_replace
 *
 * 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 <R_ext/Arith.h>        /* NA handling */
#include <Rmath.h>              /* fmin2 */
#include <R_ext/Random.h>       /* ..RNGstate */

#include <R_ext/Applic.h>       /* NA handling */

#include <R_ext/Utils.h>      /* sort routines et al */

#define RELTOL 1e-5

// create matrix with all combinations (without replacement)
void R_tc_combination_noreplace(
  int *M    // K x C output matrix in {0..N-1}
, int *N_   // N input number of elements to combine
, int *K_   // K input number of elements per combination
)
{
    int  N = *N_, K = *K_;
    int  j,k,c=0;
    int *i;
    i = Calloc(sizeof(int)*K, int);
    for (k = K-1; k>=0; k--){
        i[k] = N - (K-k);
    }
    k = 0;
    while(1){
      for(j=0;j<K;j++){
        //Rprintf("c=%d k=%d i=%d\n", j, i[j]);
        M[c+j] = i[j];
      }
      //Rprintf("\n");
      while(i[k]==k){
        k++;
        if (k==K)
          break;
      }
      if (k==K)
        break;
      i[k]--;
      while(k>0){
        k--;
        i[k] = i[k+1] - 1;
      }
      c += K;
    }
    Free(i);
}

// create matrix with all combinations (with replacement)
void R_tc_combination_replace(
  int *M    // K x C output matrix in {0..N-1}
, int *N_   // N input number of elements to combine
, int *K_   // K input number of elements per combination
)
{
    int  N = *N_, K = *K_;
    int  j,k,c=0;
    int *i;
    i = Calloc(sizeof(int)*K, int);
    for (k = K-1; k>=0; k--){
        i[k] = N - 1;
    }
    k = 0;
    while(1){
      for(j=0;j<K;j++){
        //Rprintf("c=%d k=%d i=%d\n", j, i[j]);
        M[c+j] = i[j];
      }
      //Rprintf("\n");
      while(i[k]==0){
        k++;
        if (k==K)
          break;
      }
      if (k==K)
        break;
      i[k]--;
      while(k>0){
        k--;
        i[k] = i[k+1];
      }
      c += K;
    }
    Free(i);
}


void R_tc_permutation_noreplace_sub(
  int *M    // K x C output matrix in {0..N-1}
, int *pool
, int *i
, int *c
, int N    // N input number of elements to permute
, int K    // K input number of elements per permutation
, int k     // current k
)
{
  int j,l,swapper;
  int lastk = K-1;
  for (j=k;j<N;j++){
    i[k] = pool[j];
    if (k==lastk){
      for(l=0;l<K;l++){
        M[(*c)+l] = i[l];
      }
      (*c)+=K;
    }else{
      swapper = pool[k];
      pool[k] = pool[j];
      pool[j] = swapper;
      R_tc_permutation_noreplace_sub(M,pool,i,c,N,K,k+1);
      swapper = pool[k];
      pool[k] = pool[j];
      pool[j] = swapper;
    }
  }
}

// create matrix with all permutations (without replacement)
void R_tc_permutation_noreplace(
  int *M    // K x C output matrix in {0..N-1}
, int *N_   // N input number of elements to permute
, int *K_   // K input number of elements per permutation
)
{
    int  N = *N_, K = *K_;
    int  k,c=0;
    int *i, *pool;
    i = Calloc(sizeof(int)*K, int);
    pool = Calloc(sizeof(int)*N, int);
    for (k = 0; k<N; k++){
        pool[k] = k;
    }
    R_tc_permutation_noreplace_sub(M,pool,i,&c,N,K,0);
    Free(pool);
    Free(i);
}


/* speedup not efficient, not used
void R_tc_permutation_noreplace_tryspeedup(
  int *M    // K x C output matrix in {0..N-1}
, int *N_   // N input number of elements to permute
, int *K_   // K input number of elements per permutation
)
{
    int  N = *N_, K = *K_;
    int  nComb=1,nPerm=1,nTotal=1;    // Number of combinations, subpermutations, total permutations
    int  k,p,c;
    int  cTotal,cComb,cPerm;          // Matrix row 0 pointer, incremented by K
    int *i, *pool, *Mperm;

    for (k=N;k>(N-K);k--)
      nTotal *= k;
    for (k=K;k>0;k--)
      nPerm *= k;
    nComb = nTotal / nPerm;

    //Rprintf("nComb=%d  nPerm=%d  nTotal=%d\n",nComb,nPerm,nTotal);

    // determine nComb combinations
    R_tc_combination_noreplace(M,N_,K_);
    cTotal=nComb*K;

    // permutate each combinations (nPerm)
    i = Calloc(sizeof(int)*K, int);
    pool = Calloc(sizeof(int)*K, int);
    Mperm = Calloc(sizeof(int)*K*nPerm, int);
    for (k = 0; k<K; k++){
        pool[k] = k;
    }
    cPerm=0;
    R_tc_permutation_noreplace_sub(Mperm,pool,i,&cPerm,K,K,0); // KxK not NxK
    // we skip the first permutation (the self-permutation)
    cPerm=K;
    for (p=1;p<nPerm;p++){
      cComb=0;
      for (c=0;c<nComb;c++){
        for (k=0;k<K;k++)
          M[cTotal+k] = M[cComb+Mperm[cPerm+k]];
        cTotal += K;
        cComb  += K;
      }
      cPerm  += K;
    }
    Free(Mperm);
    Free(pool);
    Free(i);
}
*/



// create matrix with all permutations (with replacement)
void R_tc_permutation_replace(
  int *M    // K x C output matrix in {0..N-1}
, int *N_   // N input number of elements to permute
, int *K_   // K input number of elements per permutation
)
{
    int  N = *N_, K = *K_;
    int  j,k,c=0;
    int *i;
    i = Calloc(sizeof(int)*K, int);
    for (k = K-1; k>=0; k--){
        i[k] = N - 1;
    }
    k = 0;
    while(1){
      for(j=0;j<K;j++){
        //Rprintf("c=%d k=%d i=%d\n", j, i[j]);
        M[c+j] = i[j];
      }
      //Rprintf("\n");
      while(i[k]==0){
        k++;
        if (k==K)
          break;
      }
      if (k==K)
        break;
      i[k]--;
      while(k>0){
        k--;
        i[k] = N - 1;
      }
      c += K;
    }
    Free(i);
}


