cluster.cc

模糊聚类的算法实现程序

int 
FCM_Clustering::Berechne_Dist (DVektor & Dist, DVektor & Datum,
			       DVektorArray & Cluster)
{
  int Result = 0, i;

  for (i = 0; i < Cluster.Lese_Groesse (); i++) {
    Dist[i] = Die_Norm.d_quadr (Datum, Cluster[i]);
    if (Dist[i] == 0)
      Result++;			/* Wieviele sind 0 ? */
  }

  return (Result);
};

void 
FCM_Clustering::Berechne_possib_U (DVektorArray & Daten,
				   DVektorArray & Cluster, DMatrix & U)
{
  int i, j;

  for (i = 0; i < Cluster.Lese_Groesse (); i++) {	/* alle Cluster */
    for (j = 0; j < Daten.Lese_Groesse (); j++) {
      U[i][j] = 1.0 / (1.0 + pow (Die_Norm.d_quadr (Daten[j], Cluster[i]) / eta[i],
			      1.0 / (M - 1.0)));
    }
  }
};

void 
FCM_Clustering::Mache_Possibilistisch (DVektorArray & Daten,
				       DVektorArray & Cluster)
{
  int l;
  DMatrix U[2];
  char Test = FALSE;

  U[0] = Zugehoerigkeit;
  U[1].DMatrix (U[0].Lese_Dim_m (), U[0].Lese_Dim_n (), NULL);

  l = 0;

  if (eta.Lese_Dim () == 0)	/* Am Anfang */
    eta = Berechne_eta (Daten, Cluster, U[l]);

  if (Possib_Zaehler <= 2) {
    while ((Test == FALSE) && (Iterationen < Max_Iterationen)) {
      Iterationen++;
      /*    v_i's berechnen    */
      Cluster = Berechne_v (U[l], Daten);
      l = (l + 1) % 2;
/*************************/
      /*    U_neu berechnen    */
/*************************/
      Berechne_possib_U (Daten, Cluster, U[l]);

      if (Die_Norm.d (U[0], U[1]) < Delta_U) {
	Test = TRUE;
	Possib_Zaehler++;
	if (Possib_Zaehler <= 2)
	  eta = Berechne_eta (Daten, Cluster, U[l]);
	else {
	  Fertig = TRUE;
	}
      }
    };
  }				/* if(Possib_Zaehler <= 2) */
  Zugehoerigkeit = U[l];
};

void 
FCM_Clustering::Berechne_Form (DVektorArray & Daten,
			       DVektorArray & Cluster)
{
  int DatenNr;
  double Max;

  Radien.Setze_Dim (Cluster.Lese_Groesse ());
  Radien.NullVektor ();

  for (DatenNr = 0; DatenNr < Daten.Lese_Groesse (); DatenNr++) {	/* alle Daten */
    Max = Cluster[Zuordnung[DatenNr]].Abstand (Daten[DatenNr],
					       Die_sind_relevant);
    if (Max > Radien[Zuordnung[DatenNr]])
      Radien[Zuordnung[DatenNr]] = Max;
  }
};

/**********************************************************/
/**********************************************************/
GK_Clustering::GK_Clustering (Guete_Typ Guete, int Soll_Iterationen,
	       char Possib)
{
  sprintf (Name, "%s", "GK");
};

GK_Clustering::~GK_Clustering ()
{
};

Cluster_Typ 
GK_Clustering::Lese_Typ ()
{
  return (GK);
};

DMatrixArray 
GK_Clustering::Lese_parallel_Form (IVektor & Dimensionen,
				   DVektorArray & Daten)
return Result (Cluster.Lese_Groesse (), Dimensionen.Lese_Dim (),
	       Dimensionen.Lese_Dim (), NULL);
{
  int ClusterNr;

  if ((Dimensionen.Lese_Dim () == 2) &&
      (Dimensionen[0] != Dimensionen[1])) {
    for (ClusterNr = 0; ClusterNr < Cluster.Lese_Groesse (); ClusterNr++) {
      Result[ClusterNr][0][0] = Laengen[ClusterNr][Dimensionen[0]];
      Result[ClusterNr][1][1] = Laengen[ClusterNr][Dimensionen[1]];
    }
  }
/************** Jetzt 3D *******************/
  else {			/* (Dimensionen.Lese_Dim == 3) */
    if ((Dimensionen[0] != Dimensionen[1]) &&
	(Dimensionen[1] != Dimensionen[2]) &&
	(Dimensionen[0] != Dimensionen[2])) {
      for (ClusterNr = 0; ClusterNr < Cluster.Lese_Groesse (); ClusterNr++) {
	Result[ClusterNr][0][0] = Laengen[ClusterNr][Dimensionen[0]];
	Result[ClusterNr][1][1] = Laengen[ClusterNr][Dimensionen[1]];
	Result[ClusterNr][2][2] = Laengen[ClusterNr][Dimensionen[2]];
      }
    }
  }
};

DMatrixArray 
GK_Clustering::Lese_Form (IVektor & Dim,
			  DVektorArray & Daten)
return Result (Cluster.Lese_Groesse (), Dim.Lese_Dim (),
	       Dim.Lese_Dim (), NULL);
{
  int ClusterNr, DatenNr, DimNr, i;
  double Winkel, a, b, Nullstellen[2], Faktor[2], Verhaeltnis, Max;

/************* 2D-Fall *********************/
  if ((Dim.Lese_Dim () == 2) && (Dim[0] != Dim[1]) &&
      Die_sind_relevant[Dim[0]] &&
      Die_sind_relevant[Dim[1]]) {
    DMatrix KS (2, 2, NULL), KS2 (2, 2, NULL);
    DVektor Rotiert (2, 0, NULL);
    char Welche_drinnen;

    for (ClusterNr = 0; ClusterNr < Cluster.Lese_Groesse (); ClusterNr++) {
      Welche_drinnen = FALSE;
      for (i = 0; i < Zuordnung.Lese_Dim (); i++) {
	if (Zuordnung[i] == ClusterNr) {
	  Welche_drinnen = TRUE;
	  break;
	}
      }
      if (Welche_drinnen) {
	if (Covar[ClusterNr][Dim[0]][Dim[0]] == Covar[ClusterNr][Dim[1]][Dim[1]]) {	/*Sonderfall */
	  Winkel = 0;
	  Result[ClusterNr][0][0] = 1.0;	/* x -   */
	  Result[ClusterNr][1][0] = 0.0;	/* Achse */
	  Result[ClusterNr][0][1] = 0.0;	/* y -   */
	  Result[ClusterNr][1][1] = 1.0;	/* Achse */
	} else {
	  Winkel = (atan (2 * Covar[ClusterNr][Dim[0]][Dim[1]] /
			  (Covar[ClusterNr][Dim[0]][Dim[0]] - Covar[ClusterNr][Dim[1]][Dim[1]])) / 2);
	  if (Winkel < 0)
	    Winkel += M_PI_2;
	  Result[ClusterNr][0][0] = cos (Winkel);
	  Result[ClusterNr][1][0] = sin (Winkel);
	  Result[ClusterNr][0][1] = -Result[ClusterNr][1][0];
	  Result[ClusterNr][1][1] = Result[ClusterNr][0][0];
	}
	/* Die Achsenverhaeltnisse der Ellipsen feststellen */
	a = -(Covar[ClusterNr][Dim[0]][Dim[0]] + Covar[ClusterNr][Dim[1]][Dim[1]]);
	b = (Covar[ClusterNr][Dim[0]][Dim[0]] * Covar[ClusterNr][Dim[1]][Dim[1]] -
	Covar[ClusterNr][Dim[0]][Dim[1]] * Covar[ClusterNr][Dim[0]][Dim[1]]);
	Nullstellen[0] = -sqrt (a * a / 4 - b) - a / 2.0;
	Nullstellen[1] = sqrt (a * a / 4 - b) - a / 2.0;

	/* lange zu kurzer Seite */
	Verhaeltnis = sqrt (fabs (Nullstellen[1] / Nullstellen[0]));	/* >= 1 */

	/* Jetzt den Faktor herausfinden */
	KS = Result[ClusterNr].Gauss_Inverse ();
	KS2 = KS;
	KS.Zeile_Multiplizieren (0, 1.0 / Verhaeltnis);
	KS2.Zeile_Multiplizieren (1, 1.0 / Verhaeltnis);

	Faktor[0] = Faktor[1] = 0.0;
	for (DatenNr = 0; DatenNr < Daten.Lese_Groesse (); DatenNr++) {
	  if (ClusterNr == Zuordnung[DatenNr]) {
	    Rotiert[0] = KS[0][0] * (Daten[DatenNr][Dim[0]] - Cluster[ClusterNr][Dim[0]]) +
	      KS[0][1] * (Daten[DatenNr][Dim[1]] - Cluster[ClusterNr][Dim[1]]);
	    Rotiert[1] = KS[1][0] * (Daten[DatenNr][Dim[0]] - Cluster[ClusterNr][Dim[0]]) +
	      KS[1][1] * (Daten[DatenNr][Dim[1]] - Cluster[ClusterNr][Dim[1]]);
	    Max = Rotiert * Rotiert;
	    if (Max > Faktor[0])
	      Faktor[0] = Max;
	    Rotiert[0] = KS2[0][0] * (Daten[DatenNr][Dim[0]] - Cluster[ClusterNr][Dim[0]]) +
	      KS2[0][1] * (Daten[DatenNr][Dim[1]] - Cluster[ClusterNr][Dim[1]]);
	    Rotiert[1] = KS2[1][0] * (Daten[DatenNr][Dim[0]] - Cluster[ClusterNr][Dim[0]]) +
	      KS2[1][1] * (Daten[DatenNr][Dim[1]] - Cluster[ClusterNr][Dim[1]]);
	    Max = Rotiert * Rotiert;
	    if (Max > Faktor[1])
	      Faktor[1] = Max;
	  }
	}
	/* Jetzt noch in die Form eintragen */
	if (Faktor[0] < Faktor[1]) {	/* Dann nicht tauschen */
	  Result[ClusterNr].Spalte_Multiplizieren (
					 0, Verhaeltnis * sqrt (Faktor[0]));
	  Result[ClusterNr].Spalte_Multiplizieren (1, sqrt (Faktor[0]));
	} else {
	  Result[ClusterNr].Spalte_Multiplizieren (
					 1, Verhaeltnis * sqrt (Faktor[1]));
	  Result[ClusterNr].Spalte_Multiplizieren (0, sqrt (Faktor[1]));
	}
      }
    }
  }
/************** Jetzt 3D *******************/
  else if (Dim.Lese_Dim () == 3) {
    if ((Dim[0] != Dim[1]) && (Dim[1] != Dim[2]) && (Dim[0] != Dim[2])) {
      DMatrix KS (3, 3, NULL);
      DVektorArray Rotiert (3, Daten.Lese_Groesse (), NULL);
      double Achse[3];
      char Welche_drinnen;

      /* Zunaechst die Verhaletnisse der Achsen herausfinden */
      for (ClusterNr = 0; ClusterNr < Cluster.Lese_Groesse (); ClusterNr++) {
	Welche_drinnen = FALSE;
	for (i = 0; i < Zuordnung.Lese_Dim (); i++) {
	  if (Zuordnung[i] == ClusterNr) {
	    Welche_drinnen = TRUE;
	    break;
	  }
	}
	if (Welche_drinnen) {
	  Achse[0] = 0.0;
	  Achse[1] = 0.0;
	  Achse[2] = 0.0;
	  Result[ClusterNr] = (Covar[ClusterNr].Projiziere (Dim)).Eigenmatrix ();
	  KS = Result[ClusterNr].Gauss_Inverse ();
	  for (DatenNr = 0; DatenNr < Daten.Lese_Groesse (); DatenNr++) {
	    if (Zuordnung[DatenNr] == ClusterNr) {
	      for (DimNr = 0; DimNr < 3; DimNr++) {
		Rotiert[DatenNr][DimNr] =
		  KS[DimNr][0] * (Daten[DatenNr][Dim[0]] - Cluster[ClusterNr][Dim[0]]) +
		  KS[DimNr][1] * (Daten[DatenNr][Dim[1]] - Cluster[ClusterNr][Dim[1]]) +
		  KS[DimNr][2] * (Daten[DatenNr][Dim[2]] - Cluster[ClusterNr][Dim[2]]);
		Achse[DimNr] += Zugehoerigkeit[ClusterNr][DatenNr] *
		  Rotiert[DatenNr][DimNr] * Rotiert[DatenNr][DimNr];
	      }
	    }
	  }			/* for(DatenNr=0 ; DatenNr < Daten.Lese_Groesse() ; DatenNr++) */
	  /* Jetzt stehen die Achsenverhaeltnisse der Ellipsen fest */
	  /* Jetzt den Faktor herausfinden, indem jeder Punkt getestet wird */
	  Faktor[0] = 0.0;
	  for (DatenNr = 0; DatenNr < Daten.Lese_Groesse (); DatenNr++) {
	    if (Zuordnung[DatenNr] == ClusterNr) {
	      Max = 0.0;
	      for (DimNr = 0; DimNr < 3; DimNr++)
		Max += Rotiert[DatenNr][DimNr] * Rotiert[DatenNr][DimNr] /
		  Achse[DimNr];
	      if (Max > Faktor[0])
		Faktor[0] = Max;
	    }
	  }
	  Result[ClusterNr].Spalte_Multiplizieren (0, sqrt (Achse[0] * Faktor[0]));
	  Result[ClusterNr].Spalte_Multiplizieren (1, sqrt (Achse[1] * Faktor[0]));
	  Result[ClusterNr].Spalte_Multiplizieren (2, sqrt (Achse[2] * Faktor[0]));
	}
      }
    }
  }
};

void 
GK_Clustering::Mein_Speichern (FILE * File)
{
  AMatrizen.Speichern (File);
  if (Art == OPTIMAL)
    ABuffer.Speichern (File);
};

int 
GK_Clustering::Mein_Laden (FILE * File)
{
  if (!AMatrizen.Laden (File))
    return (0);
  if (Art == OPTIMAL) {
    if (!ABuffer.Laden (File))
      return (0);
  }
  return (1);
};

DVektorArray 
GK_Clustering::Mein_Clustern (DVektorArray & Daten,
			      DVektorArray & Vorgabe_Protos)
return Result (Daten.Lese_Dim (), Vorgabe_Protos.Lese_Groesse (), NULL);
{
  int i, l, Clusteranzahl = Vorgabe_Protos.Lese_Groesse ();
  DMatrix U[2];

  /* Initialisierungen */
  if (Zugehoerigkeit.Lese_Dim_m () != Vorgabe_Protos.Lese_Groesse ()) {
    AMatrizen.Setze_Dim (Daten.Lese_Dim (), Daten.Lese_Dim ());
    AMatrizen.Setze_Groesse (Clusteranzahl);
    if (rho.Lese_Dim () != Clusteranzahl) {	/* Dann eben alle auf 1 setzen */
      rho.Setze_Dim (Clusteranzahl);
      for (i = 0; i < Clusteranzahl; i++) {
	rho[i] = 1;
	AMatrizen[i].Einheitsmatrix ();
      }
    }
    U[0] = DMatrix (Clusteranzahl, Daten.Lese_Groesse (), NULL);
    Berechne_U (Daten, Vorgabe_Protos, U[0]);
  } else
    U[0] = Zugehoerigkeit;
  U[1] = DMatrix (Clusteranzahl, Daten.Lese_Groesse (), NULL);

/********************/
  /* Schleifen-Anfang */
/********************/
  l = 0;
  while ((ClusterFertig != TRUE) && (Iterationen < Max_Iterationen)) {
    Iterationen++;
/***********************/
    /* II. v_i's berechnen */
/***********************/
    Result = Berechne_v (U[l], Daten);
/************************/
    /* III. A[] bestimmen  */
/************************/
    Berechne_Ai (Daten, U[l], Result);
    if (Lese_Fehlerstatus () != KEINFEHLER) {
      ClusterFertig = TRUE;
      Reset_Fehlerstatus ();
      Result. ~ DVektorArray ();
      return (Result);
    }
/***********************/
    /* IV. U_neu bestimmen */
/***********************/
    l = (l + 1) % 2;
    Berechne_U (Daten, Result, U[l]);
/*************************/
    /* IV. ||U[l-1]-U[l]||   */
/*************************/
    if (Die_Norm.d (U[0], U[1]) <= Delta_U) {
      ClusterFertig = TRUE;
    }
  };

  Zugehoerigkeit = U[l];
};

int 
GK_Clustering::Berechne_Dist (DVektor & Dist, DVektor & Datum,
			      DVektorArray & Cluster)
{
  int Result = 0, i;

  for (i = 0; i < Cluster.Lese_Groesse (); i++) {
    Dist[i] = Die_A_Norm.d_quadr (AMatrizen[i], Datum, Cluster[i]);
    if (Dist[i] == 0)
      Result++;			/* Wieviele sind 0 ? */
  }

  return (Result);
};

void 
GK_Clustering::Berechne_possib_U (DVektorArray & Daten,
				  DVektorArray & Cluster, DMatrix & U)
{
  int i, j;

  for (i = 0; i < Cluster.Lese_Groesse (); i++) {	/* alle Cluster */
    for (j = 0; j < Daten.Lese_Groesse (); j++) {
      U.Setze_i_j (i, j, 1.0 / (1.0 + pow (Die_A_Norm.d_quadr (AMatrizen[i],
					     Daten[j], Cluster[i]) / eta[i],
				       1.0 / (M - 1.0))));
    }
  }
};

void 
GK_Clustering::Mache_Possibilistisch (DVektorArray & Daten,
				      DVektorArray & Cluster)
{
  int l;
  DMatrix U[2];
  char Test = FALSE;

  U[0] = Zugehoerigkeit;
  U[1].DMatrix (U[0].Lese_Dim_m (), U[0].Lese_Dim_n (), NULL);

  l = 0;
  if (eta.Lese_Dim () == 0.0) {	/* Am Anfang */
    Berechne_Ai (Daten, U[l], Cluster);
    eta = Berechne_eta (Daten, Cluster, U[l]);
  }
  if (Possib_Zaehler <= 2) {
    while ((Test == FALSE) && (Iterationen < Max_Iterationen)) {
      Iterationen++;
      /*    v_i's berechnen    */
      Cluster = Berechne_v (U[l], Daten);
      Berechne_Ai (Daten, U[l], Cluster);
      l = (l + 1) % 2;
/*************************/
      /*    U_neu berechnen    */
/*************************/
      Berechne_possib_U (Daten, Cluster, U[l]);
      if (Die_Norm.d (U[0], U[1]) < Delta_U) {
	Test = TRUE;
	Possib_Zaehler++;
	if (Possib_Zaehler <= 2) {
	  Berechne_Ai (Daten, U[l], Cluster);
	  eta = Berechne_eta (Daten, Cluster, U[l]);
	} else {
	  Fertig = TRUE;
	}
      }
    };
  }				/* if(Possib_Zaehler <= 2) */
  Zugehoerigkeit = U[l];
};

void 
GK_Clustering::Berechne_Ai (DVektorArray & Daten, DMatrix & U,
			    DVektorArray & Cluster)
{
  int i;
  double Buffer;

  Berechne_Covarianz (Daten, U, AMatrizen, Cluster);
  Reset_Fehlerstatus ();
  for (i = 0; i < Cluster.Lese_Groesse (); i++) {
    Buffer = pow (rho[i] * AMatrizen[i].Determinante (), 1.0 / Daten.Lese_Dim ());
    AMatrizen[i] = AMatrizen[i].Gauss_Inverse ();
    if (Lese_Fehlerstatus () != KEINFEHLER)
      break;
    AMatrizen[i] *= Buffer;
  }				/* for(i=0 ; i < Cluster.Lese_Groesse() ; i++) */
}

DVektor 
GK_Clustering::Berechne_eta (DVektorArray & Daten,
			     DVektorArray & Cluster, DMatrix & U)
return Result (Cluster.Lese_Groesse (), 0, NULL);
{
  int i, k;
  double Potenz, Ni, Buffer;

  for (i = 0; i < Cluster.Lese_Groesse (); i++) {
    Buffer = 0;
    Ni = 0;
    for (k = 0; k < Daten.Lese_Groesse (); k++) {
      Potenz = pow (U[i][k], M);
      Buffer += Potenz * Die_A_Norm.d_quadr (AMatrizen[i], Daten[k], Cluster[i]);
      Ni += Potenz;
    }
    Result[i] = (K / Ni) * Buffer;
  }
}

void 
GK_Clustering::parallel_Form_Berechnen (DVektorArray & Daten,
					DVektorArray & Cluster)
{
  int ClusterNr, Dim, DatenNr;
  double Max;
  DVektorArray Achsen (Daten.Lese_Dim (), Cluster.Lese_Groesse (), NULL);
  DVektor Faktor (Cluster.Lese_Groesse (), 0, NULL);

  Laengen.Setze_Dim (Daten.Lese_Dim ());
  Laengen.Setze_Groesse (Cluster.Lese_Groesse ());

  for (DatenNr = 0; DatenNr < Daten.Lese_Groesse (); DatenNr++) {
    for (Dim = 0; Dim < Daten.Lese_Dim (); Dim++)
      Achsen[Zuordnung[DatenNr]][Dim] +=
	(Daten[DatenNr][Dim] - Cluster[Zuordnung[DatenNr]][Dim]) *
	(Daten[DatenNr][Dim] - Cluster[Zuordnung[DatenNr]][Dim]);
  }				/* for(DatenNr=0 ; DatenNr < Daten.Lese_Groesse() ; DatenNr++) */
  /* Jetzt stehen die Achsenverhaeltnisse der Ellipsen fest */

  /* Jetzt den Faktor herausfinden, indem jeder Punkt getestet wird */
  for (DatenNr = 0; DatenNr < Daten.Lese_Groesse (); DatenNr++) {
    Max = 0;
    for (Dim = 0; Dim < Daten.Lese_Dim (); Dim++) {
      if (Achsen[Zuordnung[DatenNr]][Dim] != 0)
	Max += (Daten[DatenNr][Dim] - Cluster[Zuordnung[DatenNr]][Dim]) *
	  (Daten[DatenNr][Dim] - Cluster[Zuordnung[DatenNr]][Dim]) /
	  Achsen[Zuordnung[DatenNr]][Dim];
    }
    if (Max > Faktor[Zuordnung[DatenNr]])
      Faktor[Zuordnung[DatenNr]] = Max;