/*********************************************************
** Copyright (c) 2005
** University of Washington
** Licensed under the terms set forth by University of
** Washington. If you did not sign such a license, you
** are using this software/code illegally and you do not
** have permission to use, modify, or redistribute
** this or any files in this software package.
**
** File: OdpThread.cpp  $Revision: 253 $
**
**********************************************************/
#include <iostream>
#include "edge.h"
#include "OdpThread.h"
#include "ProgressThread.h"
#include "MathUtility.h"

/*
 * Function that does the actual work
 */
int COdpThread::Go(THREAD_ID_T tid)
{
  double *pDen = new double[m_observedResultSize];
  double current = 0;
// the following arithmetic should be double to avoid overflow
// gah - changed to agree with change in odp3
//  double max = (double)m_m * (m_m + m_m0) + m_m0*(m_m0+m_m)*m_B;
// gah - odp?t compute the appropriate value for the specific odp routine.
//       this allows fine tuning of the timing calculation for each.
  double max = odp1t(m_m, m_m, m_n, m_k)+odp2t(m_m, m_m0, m_n, m_k)+
          m_B*(odp3t(m_m0, m_m, m_n, m_k)+odp4t(m_m0, m_m0, m_n, m_k));

  /* Compute the denominator */
  odp2(m_m, m_m0, m_n, m_k, m_pXX, m_pXK, m_pMu, m_pSigma0, m_pNGrp, 
    m_pNii, pDen, UpdateProgressThreadProgress, 
    this, &current, &max);

  /* Compute scores for observed statistics */
  odp1(m_m, m_m, m_n, m_k, m_pXX, m_pXK, m_pMu, m_pSigma, m_pNGrp, 
    m_pNii, m_pObservedResult, UpdateProgressThreadProgress, 
    this, &current, &max);

  DivideVector(m_pObservedResult, pDen, m_observedResultSize);

  /* Compute scores for null statistics */
// gah - try with Sigma instead of SigmaBS and m instead of m0
// gah - also Mu instead of MuBS, to match with the brca.r version.
  nullODP(m_m, m_m0, m_n, m_k, m_B, m_pRes, m_pMu, m_pSigma, m_pSigma0BS, 
    m_pVV, m_pID, m_pNii, m_pNGrp, 
    m_pNullResult, UpdateProgressThreadProgress, this, &current, &max);

  delete [] pDen;

#ifdef EDGE_DEBUG
//  cerr << max << " " << (max-current) << "\n";
  fprintf(stderr,"max=%20.0f max-current=%20.0f\n",max,max-current);
#endif
  return 0;
}

/* 
 * Function: CopyVector
 *
 * Copies an array
 */
double* COdpThread::CopyVector(double* vec, int len)
{
  return CMathUtility::CopyVector(vec, len);
}

int* COdpThread::CopyVector(int* vec, int len)
{
  return CMathUtility::CopyVector(vec, len);
}

/*
 * Frees resources associated with a COdpThread object.
 */
COdpThread::~COdpThread()
{
  delete [] m_pXX;
  delete [] m_pXK;
  delete [] m_pMu;
  delete [] m_pSigma;
  delete [] m_pSigma0;
  delete [] m_pMuBS;
  delete [] m_pSigmaBS;
  delete [] m_pSigma0BS;
  delete [] m_pID;
  delete [] m_pNGrp;
  delete [] m_pRes;
  delete [] m_pVV;
  delete [] m_pNii;
  delete [] m_pObservedResult;
  delete [] m_pNullResult;
}

/* 
 * Function:  DivideVector
 *
 * Divides each element of pNum by the corresponding 
 * member of pDenom. Both vectors must have length len.
 *
 */
void COdpThread::DivideVector(double *pNum, double *pDenom, int len)
{
  for(int i = 0; i < len; i++)
    {
      pNum[i] = pNum[i]/pDenom[i];
    }
}

