/**
 * Description not yet available.
 * \param
 */
dmatrix sin(const dmatrix& m)
{
  ivector cmin(m.rowmin(),m.rowmax());
  ivector cmax(m.rowmin(),m.rowmax());
  int i;
  for (i=m.rowmin();i<=m.rowmax();i++)
  {
    cmin(i)=m(i).indexmin();
    cmax(i)=m(i).indexmax();
  }
  dmatrix tmp(m.rowmin(),m.rowmax(),cmin,cmax);
  for (i=m.rowmin();i<=m.rowmax();i++)
  {
    tmp(i)=sin(m(i));
  }
  return tmp;
}

/**
 * Description not yet available.
 * \param
 */
dmatrix elem_prod(const dmatrix& m, const dmatrix& m2)
{
  ivector cmin(m.rowmin(),m.rowmax());
  ivector cmax(m.rowmin(),m.rowmax());
  int i;
  for (i=m.rowmin();i<=m.rowmax();i++)
  {
    cmin(i)=m(i).indexmin();
    cmax(i)=m(i).indexmax();
  }
  dmatrix tmp(m.rowmin(),m.rowmax(),cmin,cmax);
  for (i=m.rowmin();i<=m.rowmax();i++)
  {
    tmp(i)=elem_prod(m(i),m2(i));
  }
  return tmp;
}

/**
Return the sum of the diagonal of a square matrix mat.

\param mat is a square scalar matrix
*/
double trace(const dmatrix& mat)
{
  if (mat.colmin() != mat.rowmin() || mat.colmax() != mat.rowmax() )
  {
    cerr << "Matrix not square in trace\n";
    ad_exit(1);
  }

  double sum = 0.0;
  for (int i = mat.colmin(); i <= mat.colmax(); ++i)
  {
    sum += mat.elem(i, i);
  }
  return sum;
}

 double clusteringValidity::getAverageToCentroidDiameter(const dmatrix& m1) const {
   dvector a(m1.columns());
   int i,j;
   l2Distance<double> dist;
   double distance=0.0;
   for (i=0; i<m1.rows(); i++) {
     a.add(m1.getRow(i));
   }
   a.divide(m1.rows());
   for (j=0; j< m1.rows(); j++) {
     distance+=dist.apply(a,m1.getRow(j));
   }
   if (m1.rows()>0) {
     return (2*distance/(double)m1.rows());
   } else {
     return 2*distance;
   }
   
 }

 double clusteringValidity::getMaximumDistance(const dmatrix& m1,
                                               const dmatrix& m2) const {
   int i,j;
   dmatrix distances(m1.rows(),m2.rows());
   l2Distance<double> dist;
   for (i=0; i<m1.rows(); i++) {
     for (j=0; j<m2.rows(); j++) {
       distances[i][j]=dist.apply(m1.getRow(i),m2.getRow(j));
     }
   }
   return distances.maximum();
 }

 double clusteringValidity::getAverageDistance(const dmatrix& m1,
                                               const dmatrix& m2) const {
   double distance=0.0;
   int i,j;
   l2Distance<double> dist;
   for (i=0; i<m1.rows(); i++) {
     for (j=0; j<m2.rows(); j++) {
       distance+=dist.apply(m1.getRow(i),m2.getRow(j));
     }
   }
   distance=distance/((double)m1.rows()*(double)m2.rows());
   return distance;
 }

 double clusteringValidity::getAverageDiameter(const dmatrix& m1) const {
   double distance=0.0;
   int j,k;
   l2Distance<double> dist;
   for (j=0; j<m1.rows(); j++) {
     for (k=0; k<m1.rows(); k++) {
       distance+=dist.apply(m1.getRow(j),
                            m1.getRow(k));
     }
   }
   if (m1.rows()>1) {
     return (distance/((double)m1.rows()*
                       (double)(m1.rows()-1)));
   } else {
     return distance;
   }
 }

/**
 * Description not yet available.
 * \param
 */
dvariable ghk(const dvar_vector& lower,const dvar_vector& upper,
  const dvar_matrix& Sigma, const dmatrix& eps)
{
  RETURN_ARRAYS_INCREMENT();
  int n=lower.indexmax();
  int m=eps.indexmax();
  dvariable ssum=0.0;
  dvar_matrix ch=choleski_decomp(Sigma);
  dvar_vector l(1,n);
  dvar_vector u(1,n);

  for (int k=1;k<=m;k++)
  {
    dvariable weight=1.0;
    l=lower;
    u=upper;
    for (int j=1;j<=n;j++)
    {
      l(j)/=ch(j,j);
      u(j)/=ch(j,j);
      dvariable Phiu=cumd_norm(u(j));
      dvariable Phil=cumd_norm(l(j));
      weight*=Phiu-Phil;
      dvariable eta=inv_cumd_norm((Phiu-Phil)*eps(k,j)+Phil+1.e-30);
      for (int i=j+1;i<=n;i++)
      {
        dvariable tmp=ch(i,j)*eta;
        l(i)-=tmp;
        u(i)-=tmp;
      }
    }
    ssum+=weight;
  }
  RETURN_ARRAYS_DECREMENT();
  return ssum/m;
}

/**
 * Description not yet available.
 * \param
 */
dmatrix symmetrize(const dmatrix& m)
{
    if (m.rowmin() != m.colmin() || m.rowmax() != m.colmax() )
    {
        cerr << " Non square matrix passed to dmatrix symmetrize\n";
        ad_exit(1);
    }
    int rmin=m.rowmin();
    int rmax=m.rowmax();

    dmatrix s(rmin,rmax,rmin,rmax);
    for (int i=rmin; i<=rmax; i++)
    {
        s(i,i)=m(i,i);

        for (int j=rmin; j<i; j++)
        {
            s(i,j)=(m(i,j)+m(j,i))/2.;
            s(j,i)=s(i,j);
        }
    }
    return s;
}

dvariable mult_likelihood(const dmatrix &o, const dvar_matrix &p, dvar_matrix &nu, 
                          const dvariable &log_vn)
{

	// kludge to ensure observed and predicted matrixes are the same size
	if(o.colsize()!=p.colsize() || o.rowsize()!=p.rowsize())
	{
		cerr<<"Error in multivariate_t_likelihood, observed and predicted matrixes"
		" are not the same size\n";
		ad_exit(1);
	}
	dvariable vn = mfexp(log_vn);
	dvariable ff = 0.0;
	int r1 = o.rowmin();
	int r2 = o.rowmax();
	int c1 = o.colmin();
	int c2 = o.colmax();

	for(int i = r1; i <= r2; i++ )
	{
		dvar_vector sobs = vn * o(i)/sum(o(i));  //scale observed numbers by effective sample size.
		ff -= gammln(vn);
		for(int j = c1; j <= c2; j++ )
		{
			if( value(sobs(j)) > 0.0 )
				ff += gammln(sobs(j));
		}
		ff -= sobs * log(TINY + p(i));
		dvar_vector o1=o(i)/sum(o(i));
		dvar_vector p1=p(i)/sum(p(i));
		nu(i) = elem_div(o1-p1,sqrt(elem_prod(p1,1.-p1)/vn));


	}
	// exit(1);
	return ff;
}

		/**
		   calculate Pseudo-Inverse using SVD(Singular Value Decomposition)
		   by lapack library DGESVD (_a can be non-square matrix)
		*/
		int calcPseudoInverse(const dmatrix &_a, dmatrix &_a_pseu, double _sv_ratio)
		{
				int i, j, k;
				char jobu  = 'A';
				char jobvt = 'A';
				int m = (int)_a.rows();
				int n = (int)_a.cols();
				int max_mn = max(m,n);
				int min_mn = min(m,n);

				dmatrix a(m,n);
				a = _a;

				int lda = m;
				double *s = new double[max_mn];
				int ldu = m;
				double *u = new double[ldu*m];
				int ldvt = n;
				double *vt = new double[ldvt*n];
				int lwork = max(3*min_mn+max_mn, 5*min_mn);     // for CLAPACK ver.2 & ver.3
				double *work = new double[lwork];
				int info;

				for(i = 0; i < max_mn; i++) s[i] = 0.0;
		   
				dgesvd_(&jobu, &jobvt, &m, &n, &(a(0,0)), &lda, s, u, &ldu, vt, &ldvt, work,
						&lwork, &info);


				double smin, smax=0.0;
				for (j = 0; j < min_mn; j++) if (s[j] > smax) smax = s[j];
				smin = smax*_sv_ratio; 			// default _sv_ratio is 1.0e-3
				for (j = 0; j < min_mn; j++) if (s[j] < smin) s[j] = 0.0;

				//------------ calculate pseudo inverse   pinv(A) = V*S^(-1)*U^(T)
				// S^(-1)*U^(T)
				for (j = 0; j < m; j++){
						if (s[j]){
								for (i = 0; i < m; i++) u[j*m+i] /= s[j];
						}
						else {
								for (i = 0; i < m; i++) u[j*m+i] = 0.0;
						}
				}

				// V * (S^(-1)*U^(T)) 
				_a_pseu.resize(n,m);
				for(j = 0; j < n; j++){
						for(i = 0; i < m; i++){
								_a_pseu(j,i) = 0.0;
								for(k = 0; k < min_mn; k++){
										if(s[k]) _a_pseu(j,i) += vt[j*n+k] * u[k*m+i];
								}
						}
				}

				delete [] work;
				delete [] vt;
				delete [] s;
				delete [] u;

				return info;
		}

inline void orthogonalize(int N__,
                          int n__,
                          std::vector<wave_functions*> wfs__,
                          int idx_bra__,
                          int idx_ket__,
                          dmatrix<T>& o__,
                          wave_functions& tmp__)
{
    PROFILE("sddk::wave_functions::orthogonalize");

    auto pu = wfs__[0]->pu();
        
    /* project out the old subspace:
     * |\tilda phi_new> = |phi_new> - |phi_old><phi_old|phi_new> */
    if (N__ > 0) {
        inner(*wfs__[idx_bra__], 0, N__, *wfs__[idx_ket__], N__, n__, 0.0, o__, 0, 0);
        transform(pu, -1.0, wfs__, 0, N__, o__, 0, 0, 1.0, wfs__, N__, n__);
    }

    /* orthogonalize new n__ x n__ block */
    inner(*wfs__[idx_bra__], N__, n__, *wfs__[idx_ket__], N__, n__, 0.0, o__, 0, 0);

    /* single MPI rank */
    if (o__.blacs_grid().comm().size() == 1) {
        bool use_magma{false};
        #if defined(__GPU) && defined(__MAGMA)
        if (pu == GPU) {
            use_magma = true;
        }
        #endif

        if (use_magma) {
            #ifdef __GPU
            /* Cholesky factorization */
            if (int info = linalg<GPU>::potrf(n__, o__.template at<GPU>(), o__.ld())) {
                std::stringstream s;
                s << "error in GPU factorization, info = " << info;
                TERMINATE(s);
            }
            /* inversion of triangular matrix */
            if (linalg<GPU>::trtri(n__, o__.template at<GPU>(), o__.ld())) {
                TERMINATE("error in inversion");
            }
            #endif
        } else { /* CPU version */
            //check_hermitian("OVLP", o__, n__);
            //o__.serialize("overlap.dat", n__);
            /* Cholesky factorization */
            if (int info = linalg<CPU>::potrf(n__, &o__(0, 0), o__.ld())) {
                std::stringstream s;
                s << "error in factorization, info = " << info << std::endl
                  << "number of existing states: " << N__ << std::endl
                  << "number of new states: " << n__ << std::endl
                  << "number of wave_functions: " << wfs__.size() << std::endl
                  << "idx_bra: " << idx_bra__ << " " << "idx_ket:" << idx_ket__;
                TERMINATE(s);
            }
            /* inversion of triangular matrix */
            if (linalg<CPU>::trtri(n__, &o__(0, 0), o__.ld())) {
                TERMINATE("error in inversion");
            }
            if (pu == GPU) {
                #ifdef __GPU
                acc::copyin(o__.template at<GPU>(), o__.ld(), o__.template at<CPU>(), o__.ld(), n__, n__);
                #endif
            }
        }

        /* CPU version */
        if (pu == CPU) {
            /* multiplication by triangular matrix */
            for (auto& e: wfs__) {
                /* wave functions are complex, transformation matrix is complex */
                if (std::is_same<T, double_complex>::value) {
                    linalg<CPU>::trmm('R', 'U', 'N', e->pw_coeffs().num_rows_loc(), n__, double_complex(1, 0),
                                      reinterpret_cast<double_complex*>(o__.template at<CPU>()), o__.ld(),
                                      e->pw_coeffs().prime().at<CPU>(0, N__), e->pw_coeffs().prime().ld());

                    if (e->has_mt() && e->mt_coeffs().num_rows_loc()) {
                        linalg<CPU>::trmm('R', 'U', 'N', e->mt_coeffs().num_rows_loc(), n__, double_complex(1, 0),
                                          reinterpret_cast<double_complex*>(o__.template at<CPU>()), o__.ld(),
                                          e->mt_coeffs().prime().at<CPU>(0, N__), e->mt_coeffs().prime().ld());
                    }
                }
                /* wave functions are real (psi(G) = psi^{*}(-G)), transformation matrix is real */
                if (std::is_same<T, double>::value) {
                    linalg<CPU>::trmm('R', 'U', 'N', 2 * e->pw_coeffs().num_rows_loc(), n__, 1.0,
                                      reinterpret_cast<double*>(o__.template at<CPU>()), o__.ld(),
                                      reinterpret_cast<double*>(e->pw_coeffs().prime().at<CPU>(0, N__)), 2 * e->pw_coeffs().prime().ld());

                    if (e->has_mt() && e->mt_coeffs().num_rows_loc()) {
                        linalg<CPU>::trmm('R', 'U', 'N', 2 * e->mt_coeffs().num_rows_loc(), n__, 1.0,
                                          reinterpret_cast<double*>(o__.template at<CPU>()), o__.ld(),
                                          reinterpret_cast<double*>(e->mt_coeffs().prime().at<CPU>(0, N__)), 2 * e->mt_coeffs().prime().ld());
                    }
                }
            }
        }
        #ifdef __GPU
        if (pu == GPU) {
//.........这里部分代码省略.........

/**
 * Description not yet available.
 * \param
 */
dvar_matrix operator*(const dvar_matrix& m1, const dmatrix& cm2)
 {
   if (m1.colmin() != cm2.rowmin() || m1.colmax() != cm2.rowmax())
   {
     cerr << " Incompatible array bounds in "
     "dmatrix operator*(const dvar_matrix& x, const dmatrix& m)\n";
     ad_exit(21);
   }
   dmatrix cm1=value(m1);
   //dmatrix cm2=value(m2);
   dmatrix tmp(m1.rowmin(),m1.rowmax(), cm2.colmin(), cm2.colmax());
#ifdef OPT_LIB
   const size_t rowsize = (size_t)cm2.rowsize();
#else
   const int _rowsize = cm2.rowsize();
   assert(_rowsize > 0);
   const size_t rowsize = (size_t)_rowsize;
#endif
   try
   {
     double* temp_col = new double[rowsize];
     temp_col-=cm2.rowmin();
     for (int j=cm2.colmin(); j<=cm2.colmax(); j++)
     {
       for (int k=cm2.rowmin(); k<=cm2.rowmax(); k++)
       {
         temp_col[k] = cm2.elem(k,j);
       }
       for (int i=cm1.rowmin(); i<=cm1.rowmax(); i++)
       {
         double sum=0.0;
         dvector& temp_row = cm1(i);
         for (int k=cm1.colmin(); k<=cm1.colmax(); k++)
         {
           sum+=temp_row(k) * (temp_col[k]);
           // sum+=temp_row(k) * cm2(k,j);
         }
         tmp(i,j)=sum;
       }
     }
     temp_col+=cm2.rowmin();
     delete [] temp_col;
     temp_col = 0;
   }
   catch (std::bad_alloc& e)
   {
     cerr << "Error[" << __FILE__ << ':' << __LINE__
          << "]: Unable to allocate array.\n";
     //ad_exit(21);
     throw e;
   }
   dvar_matrix vtmp=nograd_assign(tmp);
   save_identifier_string("TEST1");
   //m1.save_dvar_matrix_value();
   m1.save_dvar_matrix_position();
   cm2.save_dmatrix_value();
   cm2.save_dmatrix_position();
   vtmp.save_dvar_matrix_position();
   save_identifier_string("TEST6");
   gradient_structure::GRAD_STACK1->
            set_gradient_stack(dmcm_prod);
   return vtmp;
 }

 double clusteringValidity::getAverageInterpointDistance(const dmatrix& m1,
                                                         const dmatrix& m2) const {
   l2Distance<double> dist;
   int i;
   dvector a(m1.columns());
   dvector b(m2.columns());
   for (i=0; i<m1.rows();i++) {
     a.add(m1.getRow(i));
   }
   a.divide(m1.rows()); // centroid 1
   for (i=0; i<m2.rows();i++) {
     b.add(m2.getRow(i));
   }
   b.divide(m2.rows()); // centroid 2
   double distance=0.0;
   for (i=0; i<m1.rows(); i++) {
     distance+=dist.apply(m1.getRow(i),a);
   }
   for (i=0; i<m2.rows(); i++) {
     distance+=dist.apply(m2.getRow(i),b);
   }
   return (distance/(m1.rows()+m2.rows()));
 }

  bool MLP::trainSteepestSequential(const dmatrix& data,
                                    const ivector& internalIds) {

    const parameters& param = getParameters();
    char buffer[256];
    bool abort = false;
    scramble<int> scrambler;
    int i,j,k;
    double tmpError;
    ivector idx;
    idx.resize(data.rows(),0,false,false);
    for (i=0;i<idx.size();++i) {
      idx.at(i)=i;
    }

    if (param.momentum > 0) {
      // with momentum
      dvector grad,delta(weights.size(),0.0);

      for (i=0; !abort && (i<param.maxNumberOfEpochs); ++i) {
        scrambler.apply(idx); // present the pattern in a random sequence
        totalError = 0;
        for (j=0;j<idx.size();++j) {
          k=idx.at(j);
          calcGradient(data.getRow(k),internalIds.at(k),grad);
          computeActualError(internalIds.at(k),tmpError);
          totalError+=tmpError;
          delta.addScaled(param.learnrate,grad,param.momentum,delta);
          weights.add(delta);
        }

        // update progress info object
        if (validProgressObject()) {
          sprintf(buffer,"Error=%f",totalError/errorNorm);
          getProgressObject().step(buffer);
          abort = abort || (totalError/errorNorm <= param.stopError);
          abort = abort || getProgressObject().breakRequested();
        }
      }
    } else {
      // without momentum
      ivector idx;
      idx.resize(data.rows(),0,false,false);
      dvector grad;

      int i,j,k;
      double tmpError;
      for (i=0;i<idx.size();++i) {
        idx.at(i)=i;
      }
      for (i=0; !abort && (i<param.maxNumberOfEpochs); ++i) {
        scrambler.apply(idx); // present the pattern in a random sequence
        totalError = 0;
        for (j=0;j<idx.size();++j) {
          k=idx.at(j);
          calcGradient(data.getRow(k),internalIds.at(k),grad);
          computeActualError(internalIds.at(k),tmpError);
          totalError+=tmpError;
          weights.addScaled(param.learnrate,grad);
        }

        // update progress info object
        if (validProgressObject()) {
          sprintf(buffer,"Error=%f",totalError/errorNorm);
          getProgressObject().step(buffer);
          abort = abort || (totalError/errorNorm <= param.stopError);
          abort = abort || getProgressObject().breakRequested();
        }
      }
    }
    return true;
  }

void Body::calcCMJacobian(Link *base, dmatrix &J)
{
    // prepare subm, submwc
    JointPathPtr jp;
    if (base){
        jp = getJointPath(rootLink(), base);
        Link *skip = jp->joint(0);
        skip->subm = rootLink()->m;
        skip->submwc = rootLink()->m*rootLink()->wc;
        Link *l = rootLink()->child;
        if (l){
            if (l != skip) {
                l->calcSubMassCM();
                skip->subm += l->subm;
                skip->submwc += l->submwc;
            }
            l = l->sibling;
            while(l){
                if (l != skip){
                    l->calcSubMassCM();
                    skip->subm += l->subm;
                    skip->submwc += l->submwc;
                }
                l = l->sibling;
            }
        }
        
        // assuming there is no branch between base and root
        for (int i=1; i<jp->numJoints(); i++){
            l = jp->joint(i);
            l->subm = l->parent->m + l->parent->subm;
            l->submwc = l->parent->m*l->parent->wc + l->parent->submwc;
        }
        
        J.resize(3, numJoints());
    }else{
        rootLink()->calcSubMassCM();
        J.resize(3, numJoints()+6);
    }
    
    // compute Jacobian
    std::vector<int> sgn(numJoints(), 1);
    if (jp) {
        for (int i=0; i<jp->numJoints(); i++) sgn[jp->joint(i)->jointId] = -1;
    }
    
    for (int i=0; i<numJoints(); i++){
        Link *j = joint(i);
        switch(j->jointType){
        case Link::ROTATIONAL_JOINT:
        {
            Vector3 omega(sgn[j->jointId]*j->R*j->a);
            Vector3 arm((j->submwc-j->subm*j->p)/totalMass_);
            Vector3 dp(omega.cross(arm));
            J.col(j->jointId) = dp;
            break;
        }
        default:
            std::cerr << "calcCMJacobian() : unsupported jointType("
                      << j->jointType << std::endl;
        }
    }
    if (!base){
        int c = numJoints();
        J(0, c  ) = 1.0; J(0, c+1) = 0.0; J(0, c+2) = 0.0;
        J(1, c  ) = 0.0; J(1, c+1) = 1.0; J(1, c+2) = 0.0;
        J(2, c  ) = 0.0; J(2, c+1) = 0.0; J(2, c+2) = 1.0;

        Vector3 dp(rootLink()->submwc/totalMass_ - rootLink()->p);
        J(0, c+3) =    0.0; J(0, c+4) =  dp(2); J(0, c+5) = -dp(1);
        J(1, c+3) = -dp(2); J(1, c+4) =    0.0; J(1, c+5) =  dp(0);
        J(2, c+3) =  dp(1); J(2, c+4) = -dp(0); J(2, c+5) =    0.0;
    }
}

  // implements the Fuzzy C Means algorithm
  bool fuzzyCMeans::train(const dmatrix& data) {

    bool ok=true;
    int t=0;
    // create the distance functor according to the paramter norm
    distanceFunctor<double>* distFunc = 0;
    switch (getParameters().norm)  {
      case parameters::L1:
        distFunc = new l1Distance<double>;
        break;
      case parameters::L2:
        distFunc = new l2Distance<double>;
        break;
      default:
        break;
    }
    int nbOfClusters=getParameters().nbOfClusters;
    int nbOfPoints=data.rows();
    if(nbOfClusters>nbOfPoints) {
      setStatusString("more Clusters than points");
      ok = false;
    }
    double q=getParameters().fuzzifier;
    if (q<=1) {
      setStatusString("q has to be bigger than 1");
      ok = false;
    }
    // select some points of the given data to initialise the centroids
    selectRandomPoints(data,nbOfClusters,centroids);
    // initialize variables
    centroids.resize(nbOfClusters,data.columns(),0.0);
    dmatrix memberships(nbOfPoints, nbOfClusters, 0.0);
    double terminationCriterion=0;
    double newDistance;
    dvector newCenter(data.columns());
    dvector currentPoint(data.columns());
    dmatrix newCentroids(nbOfClusters,data.columns(),0.0);
    double sumOfMemberships=0;
    double membership=0;
    double dist1;
    double dist2;
    int i,j,k,m;
    do {
        // calculate new memberships
      memberships.fill(0.0);  //  clear old memberships
      for (i=0; i<nbOfPoints; i++) {
        for (j=0; j<nbOfClusters; j++) {
          newDistance=0;
          dist1=distFunc->apply(data.getRow(i),
                                centroids.getRow(j));
          for (k=0; k<nbOfClusters; k++) {
            dist2=distFunc->apply(data.getRow(i),
                                  centroids.getRow(k));
       // if distance is 0, normal calculation of membership is not possible.
            if (dist2!=0) {
              newDistance+=pow((dist1/dist2),(1/(q-1)));
            }
          }
      // if point and centroid are equal
          if (newDistance!=0)
            memberships.at(i,j)=1/newDistance;
          else {
            dvector row(memberships.columns(),0.0);
            memberships.setRow(i,row);
            memberships.at(i,j)=1;
            break;
          }
        }
      }
      t++;  // counts the iterations

     // calculate new centroids based on modified memberships
      for (m=0; m<nbOfClusters; m++) {
        newCenter.fill(0.0);
        sumOfMemberships=0;
        for (i=0; i<nbOfPoints; i++) {
          currentPoint=data.getRow(i);
          membership=pow(memberships.at(i,m),q);
          sumOfMemberships+=membership;
          currentPoint.multiply(membership);
          newCenter.add(currentPoint);
        }
        newCenter.divide(sumOfMemberships);
        newCentroids.setRow(m,newCenter);
      }
      terminationCriterion=distFunc->apply(centroids,newCentroids);
      centroids=newCentroids;
    }
    // the termination criterions
    while ( (terminationCriterion>getParameters().epsilon)
            && (t<getParameters().maxIterations));

    int nbClusters = nbOfClusters;
    //Put the id information into the result object
    //Each cluster has the id of its position in the matrix
    ivector tids(nbClusters);
    for (i=0; i<nbClusters; i++) {
      tids.at(i)=i;
    }
//.........这里部分代码省略.........

		/**
		   solve linear equation using LU decomposition
		   by lapack library DGESVX (_a must be square matrix)
		*/
		int solveLinearEquationLU(const dmatrix &_a, const dvector &_b, dvector &_x)
		{
				assert(_a.cols() == _a.rows() && _a.cols() == _b.size() );

				int n = (int)_a.cols();
				int nrhs = 1;

				int lda = n;

				std::vector<int> ipiv(n);

				int ldb = n;

				int info;

				// compute the solution
#ifndef USE_CLAPACK_INTERFACE
  				char fact      = 'N';
				char transpose = 'N';

				double *af = new double[n*n];

				int ldaf = n;

				char equed = 'N';

				double *r = new double[n];
				double *c = new double[n];

				int ldx = n;

				double rcond;

				double *ferr = new double[nrhs];
				double *berr = new double[nrhs];
				double *work = new double[4*n];

				int *iwork = new int[n];

			    _x.resize(n);			// memory allocation for the return vector
				dgesvx_(&fact, &transpose, &n, &nrhs, const_cast<double *>(&(_a(0,0))), &lda, af, &ldaf, &(ipiv[0]),
						&equed, r, c, const_cast<double *>(&(_b(0))), &ldb, &(_x(0)), &ldx, &rcond,
						ferr, berr, work, iwork, &info);

				delete [] iwork;
				delete [] work;
				delete [] berr;
				delete [] ferr;
				delete [] c;
				delete [] r;

				delete [] af;
#else
				_x = _b;
				info = clapack_dgesv(CblasColMajor,
									 n, nrhs, const_cast<double *>(&(a(0,0))), lda, &(ipiv[0]),
									 &(_x(0)), ldb);
#endif

				return info;
		}

 bool normModHubertStat::apply(const std::vector<dmatrix>& 
                               clusteredData, double& index, 
                               const dmatrix& centroids) const {
   index =0.0;
   int nbClusters=clusteredData.size();
   l2Distance<double> dist;
   int nbPoints=0;
   int i,j,k,l;
   // count the points that are in clusteredData
   for (i=0; i<nbClusters; i++) {
     nbPoints+=clusteredData[i].rows();
   }
   // x is the distance matrix. It has in the i-th rows and j-th column
   // the distance between the i-th and j-th point
   // y is an other distance matrix. It has in the i-th row and j-th column
   // the distance of the centroids of the clusters they belong to.
   dmatrix x(nbPoints,nbPoints);
   dmatrix y(nbPoints,nbPoints);
   int row=0;
   int col=0;
   for (i=0; i<nbClusters; i++) {
     for (j=0; j<clusteredData[i].rows(); j++) {
       for (k=0; k<nbClusters; k++) {
         for (l=0; l<clusteredData[k].rows(); l++) {
           if (col>row) {
             y.at(row,col)=dist.apply(centroids.getRow(i),
                                      centroids.getRow(k));
             x.at(row,col)=dist.apply(clusteredData[i].getRow(j),
                                      clusteredData[k].getRow(l));
           }
           if (col<nbPoints-1) col++;
           else { 
             col=0;
             row++;
           }  
         }
       }
     }
   }
   double m=0.5*(nbPoints*(nbPoints-1));
   double meanX=x.sumOfElements()/m;
   double meanY=y.sumOfElements()/m;
   double tmp1=meanX*meanX;
   double tmp2=meanY*meanY;
   double varianzX=0.0;
   double varianzY=0.0;
   for (i=0; i<nbPoints; i++) {
     for (j=i+1; j<nbPoints; j++) {
       varianzX+=x[i][j]*x[i][j]-tmp1;
       varianzY+=y[i][j]*y[i][j]-tmp2;
     }
   }
   varianzX=varianzX/m;
   varianzY=varianzY/m;
   varianzX=sqrt(varianzX);
   varianzY=sqrt(varianzY);
   double varianz=varianzX*varianzY;
   for (i=0; i<nbPoints; i++) {
     for (j=i+1; j<nbPoints; j++) {
       index+=(x[i][j]-meanX)*(y[i][j]-meanY);
     }
   }
   index=index/(m*varianz); 
 
   return true;
 }

  bool dunnIndex::apply(const std::vector<dmatrix>& clusteredData,
                        double& index, const dmatrix& centroids) const {
    int nbClusters=clusteredData.size();
    double denominator=0.0;
    int i,j;
    l2Distance<double> dist;
    dvector diameters(nbClusters);
    parameters param;
    param=getParameters();
    // pointer to the function which implements the measure according to the 
    // parameters
    double (lti::clusteringValidity::*diamFunc)(const dmatrix&) const;
    
#ifdef _LTI_MSC_DOT_NET_2003
    // nasty bug in this version of the .NET compiler
#define QUALIFIER
#else 
#define QUALIFIER &lti::dunnIndex::
#endif

    switch (param.diameterMeasure) {
      case parameters::Standard:
        diamFunc=QUALIFIER getStandardDiameter;
        break;
      case parameters::Average:
        diamFunc=QUALIFIER getAverageDiameter;
        break;
      case parameters::Centroid:
        diamFunc=QUALIFIER getAverageToCentroidDiameter;
        break;
      default:
        diamFunc=QUALIFIER getStandardDiameter;
        setStatusString("Unknown diameterMeasure in clusteringValidity\n");
        return false;
    }
    // compute all diameters of all clusters
    for (i=0; i<nbClusters; i++) {
      diameters[i]=(this->*diamFunc)(clusteredData[i]);
    }
    denominator=diameters.maximum();
    // pointer to the function which calculates the distance of the functions
    // a pointer to a function is used, because the function will be called 
    // many times later
    double (lti::clusteringValidity::*distFunc)
      (const dmatrix&,const dmatrix&) const ;
    // set pointer to function which is set by the parameter distanceMeasure
    switch (param.distanceMeasure) {
      case parameters::Minimum:
        distFunc=QUALIFIER getMinimumDistance;
        break;
      case parameters::Maximum:
        distFunc=QUALIFIER getMaximumDistance;
        break;
      case parameters::Mean:
        distFunc=QUALIFIER getAverageDistance;
        break;
      case parameters::Centroids:
        distFunc=QUALIFIER getCentroidDistance;
        break;
      case parameters::Interpoint:
        distFunc=QUALIFIER getAverageInterpointDistance;
        break;
      default:
        distFunc=QUALIFIER getAverageDistance;
        setStatusString("Unknown distanceMeasure in clusteringValidity\n");
        return false;
    }
    // compute the distances of all clusters to each other
    int counter=0;
    dvector distanceVector(static_cast<int>(.5*(nbClusters*(nbClusters-1))));
    for (i=0; i<nbClusters; i++) {
      for (j=i+1; j<nbClusters; j++) {
        if (distFunc==QUALIFIER getCentroidDistance) {
          distanceVector[counter]=dist.apply(centroids.getRow(i),
                                             centroids.getRow(j));
        } else {
          distanceVector[counter]=(this->*distFunc)(clusteredData[i],
                                                    clusteredData[j]);
        }
        counter++;
      }
    }
    distanceVector.divide(denominator);
    index=distanceVector.minimum();

    return true;
  }