/* truecluster/src/matchindex.c
 * Truecluster matching heuristic
 * (C) 2000-2007 Dr. Jens Oehlschlägel
 *
 * 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

void R_tc_matchindex(
  int *observed   // integer input matrix
, int *rn_        // no. of rows
, int *cn_        // no. of cols
, int *rpos      // integer output rows
, int *cpos      // integer output cols
)
{
  int rn= *rn_, cn= *cn_;
  int rn1 = rn;
  int r, c, i, imax, nmax, ri, ci;
  int tmpint;
  double tmpdouble;
  double tsum, maxl, maxu;
  int maxo;
  double expected, res, ssres;
  Rboolean used_random=FALSE;
  double *rsum = Calloc(sizeof(double)*rn, double);
  double *csum = Calloc(sizeof(double)*cn, double);

  while(rn>1 && cn>1){

    maxl  = R_NegInf;
    maxu  = R_NegInf;
    maxo = -1;

    // calculate row sums and total sum
    tsum = 0;
    for (r=0; r<rn; r++){
      tmpdouble = 0;
      for (c=0; c<cn; c++){
        tmpdouble += observed[rpos[r]+cpos[c]*rn1];
      }
      rsum[r] = tmpdouble;
      tsum += tmpdouble;
    }

    // calculate col sums divided by total sum
    for (c=0; c<cn; c++){
      tmpdouble = 0;
      for (r=0; r<rn; r++){
        tmpdouble += observed[rpos[r]+cpos[c]*rn1];
      }
      csum[c] = tmpdouble / tsum;
    }

    // calculate expected
    for (r=0; r<rn; r++){
      for (c=0; c<cn; c++){
        i = rpos[r]+cpos[c]*rn1;
        expected = rsum[r] * csum[c];
        res = observed[i] - expected;
        if (expected==0)
          ssres = 0;
        else if (res<0)
          ssres = -1.0 * res*res / expected;
        else
          ssres = res*res / expected;
        //Rprintf("r=%d c=%d maxl=%f ssres=%f maxu=%f maxo=%d obs=%d\n", r, c, maxl, ssres, maxu, maxo, observed[i]);
        if ((ssres>maxu) || (ssres>maxl && observed[i]>maxo)){
          nmax = 1;
          imax = r + c*rn;
          maxl = ssres - ssres*RELTOL;
          maxu = ssres + ssres*RELTOL;
          maxo = observed[i];
          //Rprintf("NEWMAX maxl=%f ssres=%f maxu=%f \n", maxl, ssres, maxu);
        }else if (ssres>maxl){
          nmax++;
          if (!used_random) {
            GetRNGstate(); used_random = TRUE;
          }
          tmpdouble = (1.0/nmax);
          //Rprintf("EQUAL p=%f\n", tmpdouble);
          if (unif_rand()>tmpdouble)
            imax = r + c*rn;
        }
      }
    }

    ri = imax % rn;
    ci = imax / rn;
    //Rprintf("rn=%d cn=%d imax=%d ri=%d ci=%d\n", rn, cn, imax, ri, ci);
    tmpint   = rpos[rn-1];
    rpos[rn-1] = rpos[ri];
    rpos[ri] = tmpint;
    tmpint   = cpos[cn-1];
    cpos[cn-1] = cpos[ci];
    cpos[ci] = tmpint;

    rn--;
    cn--;
  }
  if(used_random) {
    PutRNGstate();
  }
  Free(rsum);
  Free(csum);

}
