+/*********************************************************************************************************************************/
+double LinearAlgebra::calcKendall(vector<double>& x, vector<double>& y, double& sig){
+ try {
+ if (x.size() != y.size()) { m->mothurOut("[ERROR]: vector size mismatch."); m->mothurOutEndLine(); return 0.0; }
+
+ //format data
+ vector<spearmanRank> xscores;
+ for (int i = 0; i < x.size(); i++) {
+ spearmanRank member(toString(i), x[i]);
+ xscores.push_back(member);
+ }
+
+ //sort xscores
+ sort(xscores.begin(), xscores.end(), compareSpearman);
+
+ //convert scores to ranks of x
+ vector<spearmanRank*> ties;
+ int rankTotal = 0;
+ for (int j = 0; j < xscores.size(); j++) {
+ rankTotal += (j+1);
+ ties.push_back(&(xscores[j]));
+
+ if (j != xscores.size()-1) { // you are not the last so you can look ahead
+ if (xscores[j].score != xscores[j+1].score) { // you are done with ties, rank them and continue
+ for (int k = 0; k < ties.size(); k++) {
+ float thisrank = rankTotal / (float) ties.size();
+ (*ties[k]).score = thisrank;
+ }
+ ties.clear();
+ rankTotal = 0;
+ }
+ }else { // you are the last one
+ for (int k = 0; k < ties.size(); k++) {
+ float thisrank = rankTotal / (float) ties.size();
+ (*ties[k]).score = thisrank;
+ }
+ }
+ }
+
+
+ //format data
+ vector<spearmanRank> yscores;
+ for (int j = 0; j < y.size(); j++) {
+ spearmanRank member(toString(j), y[j]);
+ yscores.push_back(member);
+ }
+
+ //sort yscores
+ sort(yscores.begin(), yscores.end(), compareSpearman);
+
+ //convert to ranks
+ map<string, float> rank;
+ vector<spearmanRank> yties;
+ rankTotal = 0;
+ for (int j = 0; j < yscores.size(); j++) {
+ rankTotal += (j+1);
+ yties.push_back(yscores[j]);
+
+ if (j != yscores.size()-1) { // you are not the last so you can look ahead
+ if (yscores[j].score != yscores[j+1].score) { // you are done with ties, rank them and continue
+ for (int k = 0; k < yties.size(); k++) {
+ float thisrank = rankTotal / (float) yties.size();
+ rank[yties[k].name] = thisrank;
+ }
+ yties.clear();
+ rankTotal = 0;
+ }
+ }else { // you are the last one
+ for (int k = 0; k < yties.size(); k++) {
+ float thisrank = rankTotal / (float) yties.size();
+ rank[yties[k].name] = thisrank;
+ }
+ }
+ }
+
+
+ int numCoor = 0;
+ int numDisCoor = 0;
+
+ //associate x and y
+ vector<spearmanRank> otus;
+ for (int l = 0; l < xscores.size(); l++) {
+ spearmanRank member(xscores[l].name, rank[xscores[l].name]);
+ otus.push_back(member);
+ }
+
+ int count = 0;
+ for (int l = 0; l < xscores.size(); l++) {
+
+ int numWithHigherRank = 0;
+ int numWithLowerRank = 0;
+ float thisrank = otus[l].score;
+
+ for (int u = l+1; u < xscores.size(); u++) {
+ if (otus[u].score > thisrank) { numWithHigherRank++; }
+ else if (otus[u].score < thisrank) { numWithLowerRank++; }
+ count++;
+ }
+
+ numCoor += numWithHigherRank;
+ numDisCoor += numWithLowerRank;
+ }
+
+ double p = (numCoor - numDisCoor) / (float) count;
+
+ sig = calcKendallSig(x.size(), p);
+
+ return p;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "LinearAlgebra", "calcKendall");
+ exit(1);
+ }
+}
+double LinearAlgebra::ran0(int& idum)
+{
+ const int IA=16807,IM=2147483647,IQ=127773;
+ const int IR=2836,MASK=123459876;
+ const double AM=1.0/double(IM);
+ int k;
+ double ans;
+
+ idum ^= MASK;
+ k=idum/IQ;
+ idum=IA*(idum-k*IQ)-IR*k;
+ if (idum < 0) idum += IM;
+ ans=AM*idum;
+ idum ^= MASK;
+ return ans;
+}
+
+double LinearAlgebra::ran1(int &idum)
+{
+ const int IA=16807,IM=2147483647,IQ=127773,IR=2836,NTAB=32;
+ const int NDIV=(1+(IM-1)/NTAB);
+ const double EPS=3.0e-16,AM=1.0/IM,RNMX=(1.0-EPS);
+ static int iy=0;
+ static vector<int> iv(NTAB);
+ int j,k;
+ double temp;
+
+ if (idum <= 0 || !iy) {
+ if (-idum < 1) idum=1;
+ else idum = -idum;
+ for (j=NTAB+7;j>=0;j--) {
+ k=idum/IQ;
+ idum=IA*(idum-k*IQ)-IR*k;
+ if (idum < 0) idum += IM;
+ if (j < NTAB) iv[j] = idum;
+ }
+ iy=iv[0];
+ }
+ k=idum/IQ;
+ idum=IA*(idum-k*IQ)-IR*k;
+ if (idum < 0) idum += IM;
+ j=iy/NDIV;
+ iy=iv[j];
+ iv[j] = idum;
+ if ((temp=AM*iy) > RNMX) return RNMX;
+ else return temp;
+}
+
+double LinearAlgebra::ran2(int &idum)
+{
+ const int IM1=2147483563,IM2=2147483399;
+ const int IA1=40014,IA2=40692,IQ1=53668,IQ2=52774;
+ const int IR1=12211,IR2=3791,NTAB=32,IMM1=IM1-1;
+ const int NDIV=1+IMM1/NTAB;
+ const double EPS=3.0e-16,RNMX=1.0-EPS,AM=1.0/double(IM1);
+ static int idum2=123456789,iy=0;
+ static vector<int> iv(NTAB);
+ int j,k;
+ double temp;
+
+ if (idum <= 0) {
+ idum=(idum==0 ? 1 : -idum);
+ idum2=idum;
+ for (j=NTAB+7;j>=0;j--) {
+ k=idum/IQ1;
+ idum=IA1*(idum-k*IQ1)-k*IR1;
+ if (idum < 0) idum += IM1;
+ if (j < NTAB) iv[j] = idum;
+ }
+ iy=iv[0];
+ }
+ k=idum/IQ1;
+ idum=IA1*(idum-k*IQ1)-k*IR1;
+ if (idum < 0) idum += IM1;
+ k=idum2/IQ2;
+ idum2=IA2*(idum2-k*IQ2)-k*IR2;
+ if (idum2 < 0) idum2 += IM2;
+ j=iy/NDIV;
+ iy=iv[j]-idum2;
+ iv[j] = idum;
+ if (iy < 1) iy += IMM1;
+ if ((temp=AM*iy) > RNMX) return RNMX;
+ else return temp;
+}
+
+double LinearAlgebra::ran3(int &idum)
+{
+ static int inext,inextp;
+ static int iff=0;
+ const int MBIG=1000000000,MSEED=161803398,MZ=0;
+ const double FAC=(1.0/MBIG);
+ static vector<int> ma(56);
+ int i,ii,k,mj,mk;
+
+ if (idum < 0 || iff == 0) {
+ iff=1;
+ mj=labs(MSEED-labs(idum));
+ mj %= MBIG;
+ ma[55]=mj;
+ mk=1;
+ for (i=1;i<=54;i++) {
+ ii=(21*i) % 55;
+ ma[ii]=mk;
+ mk=mj-mk;
+ if (mk < int(MZ)) mk += MBIG;
+ mj=ma[ii];
+ }
+ for (k=0;k<4;k++)
+ for (i=1;i<=55;i++) {
+ ma[i] -= ma[1+(i+30) % 55];
+ if (ma[i] < int(MZ)) ma[i] += MBIG;
+ }
+ inext=0;
+ inextp=31;
+ idum=1;
+ }
+ if (++inext == 56) inext=1;
+ if (++inextp == 56) inextp=1;
+ mj=ma[inext]-ma[inextp];
+ if (mj < int(MZ)) mj += MBIG;
+ ma[inext]=mj;
+ return mj*FAC;
+}
+
+double LinearAlgebra::ran4(int &idum)
+{
+#if defined(vax) || defined(_vax_) || defined(__vax__) || defined(VAX)
+ static const unsigned long jflone = 0x00004080;
+ static const unsigned long jflmsk = 0xffff007f;
+#else
+ static const unsigned long jflone = 0x3f800000;
+ static const unsigned long jflmsk = 0x007fffff;
+#endif
+ unsigned long irword,itemp,lword;
+ static int idums = 0;
+
+ if (idum < 0) {
+ idums = -idum;
+ idum=1;
+ }
+ irword=idum;
+ lword=idums;
+ psdes(lword,irword);
+ itemp=jflone | (jflmsk & irword);
+ ++idum;
+ return (*(float *)&itemp)-1.0;
+}
+
+void LinearAlgebra::psdes(unsigned long &lword, unsigned long &irword)
+{
+ const int NITER=4;
+ static const unsigned long c1[NITER]={
+ 0xbaa96887L, 0x1e17d32cL, 0x03bcdc3cL, 0x0f33d1b2L};
+ static const unsigned long c2[NITER]={
+ 0x4b0f3b58L, 0xe874f0c3L, 0x6955c5a6L, 0x55a7ca46L};
+ unsigned long i,ia,ib,iswap,itmph=0,itmpl=0;
+
+ for (i=0;i<NITER;i++) {
+ ia=(iswap=irword) ^ c1[i];
+ itmpl = ia & 0xffff;
+ itmph = ia >> 16;
+ ib=itmpl*itmpl+ ~(itmph*itmph);
+ irword=lword ^ (((ia = (ib >> 16) |
+ ((ib & 0xffff) << 16)) ^ c2[i])+itmpl*itmph);
+ lword=iswap;
+ }
+}
+/*********************************************************************************************************************************/
+double LinearAlgebra::calcKendallSig(double n, double r){
+ try {
+
+ double sig = 0.0;
+ double svar=(4.0*n+10.0)/(9.0*n*(n-1.0));
+ double z= r/sqrt(svar);
+ sig=erfcc(fabs(z)/1.4142136);
+
+ if (isnan(sig) || isinf(sig)) { sig = 0.0; }
+
+ return sig;
+ }
+ catch(exception& e) {
+ m->errorOut(e, "LinearAlgebra", "calcKendallSig");
+ exit(1);
+ }
+}