regel.cc

UNIX下使用的基于模糊Ｃ聚类分析的原码

/********************************************************/
/*   filename: regel.cc                                 */
/*                                                      */
/********************************************************/
/* programmed by: Oliver Wagner                         */
/* last change: 22-03-95                                */
/********************************************************/
#include "regel.h"

/********************************************************/
/*          Alles zur Klasse Trapez_Parameter           */
/********************************************************/
Trapez_Parameter::Trapez_Parameter ()
{
  /* hier muss noch was rein */

  So = Trapez;
};

Trapez_Parameter::Trapez_Parameter (const Trapez_Parameter & Das)
{
  So = Trapez;
  Daten[0] = Das.Daten[0];
  Daten[1] = Das.Daten[1];
  Daten[2] = Das.Daten[2];
  Daten[3] = Das.Daten[3];
  Steigung[0] = Das.Steigung[0];
  Steigung[1] = Das.Steigung[1];
};

Trapez_Parameter::Trapez_Parameter (const Parameter & Das)
{
  if (Das.So == Trapez) {
    So = Trapez;
    Daten[0] = ((Trapez_Parameter &) Das).Daten[0];
    Daten[1] = ((Trapez_Parameter &) Das).Daten[1];
    Daten[2] = ((Trapez_Parameter &) Das).Daten[2];
    Daten[3] = ((Trapez_Parameter &) Das).Daten[3];
    Steigung[0] = ((Trapez_Parameter &) Das).Steigung[0];
    Steigung[1] = ((Trapez_Parameter &) Das).Steigung[1];
  } else
    Fehlermeldung ("Trapez_Parameter::Trapez_Parameter(const Parameter& Das)", KOMPATIBILITAETSFEHLER);
};

void 
Trapez_Parameter::print ()
{
  printf ("\nTrapez_Parameter: %f 0, %f 1, %f 1, %f 0\n", Daten[0], Daten[1],
	  Daten[2], Daten[3]);
};

void 
Trapez_Parameter::Setze_Parameter (Parameter & Damit)
{
  int i;

  for (i = 0; i < 4; i++)
    Daten[i] = Damit.Lese_xParameter (i);
  if (Daten[0] != Daten[1])
    Steigung[0] = 1 / (Daten[1] - Daten[0]);
  if (Daten[2] != Daten[3])
    Steigung[1] = 1 / (Daten[2] - Daten[3]);
};

void 
Trapez_Parameter::Setze_xParameter (int Index, double Damit)
{
  Daten[Index] = Damit;
  if (Index < 2) {		/* 1 oder 0 */
    if (Daten[0] != Daten[1])
      Steigung[0] = 1 / (Daten[1] - Daten[0]);
  } else {
    if (Daten[2] != Daten[3])
      Steigung[1] = 1 / (Daten[2] - Daten[3]);
  }
};

double 
Trapez_Parameter::Lese_yParameter (int i)
{
  if ((i == 0) || (i == 3))
    return (0);
  else
    return (1);
};


Parameter & Trapez_Parameter::operator = (Parameter & Das) {
  int i;

  if (Das.So == Trapez) {
    for (i = 0; i < 4; i++)
      Daten[i] = Das.Lese_xParameter (i);
    Steigung[0] = ((Trapez_Parameter &) Das).Steigung[0];
    Steigung[1] = ((Trapez_Parameter &) Das).Steigung[1];
  } else
    Fehlermeldung ("Trapez_Parameter::operator=()", KOMPATIBILITAETSFEHLER);
};

char 
Trapez_Parameter::Finde_Parameter (DVektorArray & Daraus,
				   double Z_Grenze)
{
  double ptest[4], f, Fehler, Fehlerbuffer;
  int i, k, Index = 0, Schritt;
  DVektorArray Kurve;

  Kurve = Daraus.Quicksort (0);
  Kurve = Kurve.Mache_Konvex (1, Z_Grenze);

  if (Kurve.Lese_Groesse () == 0)
    return (FALSE);		/* Das ging so nicht */
  /* Initialisierung der Parameter */
  Daten[0] = Kurve.Lese_Min (0);
  Daten[3] = Kurve.Lese_Max (0);
//  Daten[1]=Daten[0]+(Daten[3]-Daten[0])/3;
  //  Daten[2]=Daten[1]+(Daten[3]-Daten[0])/3;
  Daten[1] = Daten[2] = Kurve[Kurve.Lese_MaxNummer (1)][0];

  f = PARAMETER_F_PROZENT * (Daten[3] - Daten[0]);

  /* jetzt alle Parameter nacheinander aendern */
  for (Schritt = 0; Schritt < PARAMETER_SCHRITTE; Schritt++) {
    for (k = 0; k < 4; k++) {
      /* zunaechst die moeglichen anderen Parameter errechnen */
      ptest[0] = Daten[k] - f;
      if (k != 0) {		/* 1, 2 oder 3 */
	if (ptest[0] < Daten[k - 1])
	  ptest[0] = Daten[k - 1];
      } else {
	if (ptest[0] < Kurve.Lese_Min (0))
	  ptest[0] = Kurve.Lese_Min (0);
      }
      ptest[1] = Daten[k];
      ptest[2] = Daten[k] + f;
      if (k != 3) {		/* 0, 1 oder 2 */
	if (ptest[2] > Daten[k + 1])
	  ptest[2] = Daten[k + 1];
      } else {
	if (ptest[3] < Kurve.Lese_Max (0))
	  ptest[3] = Kurve.Lese_Max (0);
      }
      /* Jetzt die moeglichen Parameter ausprobieren */
      for (i = 0, Fehler = MAXDOUBLE; i < 3; i++) {
	Daten[k] = ptest[i];
	Fehlerbuffer = Parameter_Fehler (Kurve);
	if (Fehlerbuffer < Fehler) {
	  Index = i;
	  Fehler = Fehlerbuffer;
	}
      }
      Daten[k] = ptest[Index];
    }
  }

  if (Daten[0] != Daten[1])
    Steigung[0] = 1 / (Daten[1] - Daten[0]);
  if (Daten[2] != Daten[3])
    Steigung[1] = 1 / (Daten[2] - Daten[3]);

  return (TRUE);
};

void 
Trapez_Parameter::Speichern (FILE * File)
{
  fprintf (File, "%lf %lf %lf %lf\n", Daten[0], Daten[1], Daten[2], Daten[3]);
};

double 
Trapez_Parameter::Berechne_Schwerpunkt (double Hoehe, DVektor & x,
					DVektor & y)
{
  double Integral_fx = 0, Integral_f = 0, Result;

  Berechne_Form (Hoehe, x, y);

  /* erstes Segment */
  if (Hoehe != 0) {
    if (x[1] != x[0]) {		/* nicht senkrecht */
      Integral_fx += Steigung[0] * ((pow (x[1], 3) - pow (x[0], 3)) / 3 -
				x[0] * (pow (x[1], 2) - pow (x[0], 2)) / 2);
      Integral_f += Hoehe * (x[1] - x[0]) / 2;
    }
    /* letztes Segment */
    if (x[3] != x[2]) {
      Integral_fx += Steigung[1] * ((pow (x[3], 3) - pow (x[2], 3)) / 3) +
	((Hoehe - x[2] * Steigung[1]) * (pow (x[3], 2) - pow (x[2], 2)) / 2);
      Integral_f += Hoehe * (x[3] - x[2]) / 2;
    }
    /* mittleres Segment */
    Integral_fx += Hoehe * (pow (x[2], 2) - pow (x[1], 2)) / 2;
    Integral_f += (x[2] - x[1]) * Hoehe;

    if (Integral_f == 0) {	/* Falls nur "Ausgabe-Peak" */
      Result = x[0];
    } else {
      Result = Integral_fx / Integral_f;
    }
  } else {
    Result = x[0] + (x[3] - x[0]) / 2;
  }
  return (Result);
};

double 
Trapez_Parameter::Berechne_Hoehe (double Eingabe)
{
  double Result;

  if ((Eingabe >= Daten[0]) && (Eingabe <= Daten[3])) {		/* im Trapez */
    if (Eingabe < Daten[1]) {
      Result = (Eingabe - Daten[0]) / (Daten[1] - Daten[0]);
    } else if (Eingabe > Daten[2]) {
      Result = (Eingabe - Daten[3]) / (Daten[2] - Daten[3]);
    } else
      Result = 1;
  } else
    Result = 0;

  return (Result);
};

void 
Trapez_Parameter::Berechne_Form (double Hoehe, DVektor & x,
				 DVektor & y)
{
  x.Setze_Dim (4);
  y.Setze_Dim (4);

  x[0] = Daten[0];
  x[3] = Daten[3];
  y[0] = y[3] = 0;
  y[1] = y[2] = Hoehe;

  if (Hoehe != 0) {
    /* erstes Segment */
    if (Daten[1] != Daten[0]) {	/* nicht senkrecht */
      x[1] = Daten[0] + Hoehe / Steigung[0];
    } else
      x[1] = Daten[1];
    /* letztes Segment */
    if (Daten[3] != Daten[2]) {
      x[2] = Daten[3] + Hoehe / Steigung[1];
    } else
      x[2] = Daten[3];
    /* mittleres Segment */
  } else {
    x[1] = Daten[0];
    x[2] = Daten[3];
  }
};


double 
Trapez_Parameter::Parameter_Fehler (DVektorArray & Damit)
{
  double Result = 0;
  int i = 0;

  while ((i < Damit.Lese_Groesse ()) && (Damit[i][0] < Daten[0]))
    Result += pow (Damit[i++][1], 2);
  while ((i < Damit.Lese_Groesse ()) && (Damit[i][0] < Daten[1])) {
    Result += pow (Damit[i][1] - (Damit[i][0] - Daten[0]) /
		   (Daten[1] - Daten[0]), 2);
    i++;
  }
  while ((i < Damit.Lese_Groesse ()) && (Damit[i][0] < Daten[2]))
    Result += pow (1 - Damit[i++][1], 2);
  while ((i < Damit.Lese_Groesse ()) && (Damit[i][0] < Daten[3])) {
    Result += pow (Damit.Lese_Vektor (i, 1) - (Damit[i][0] - Daten[2]) /
		   (Daten[3] - Daten[2]), 2);
    i++;
  }
  while (i < Damit.Lese_Groesse ())
    Result += pow (Damit[i++][1], 2);

  return (Result);
};

Trapez_Parameter 
Trapez_einlesen (FILE * File)
return Result;
{
  if (4 != fscanf (File, "%lf %lf %lf %lf", &(Result.Daten[0]), &(Result.Daten[1]),
		   &(Result.Daten[2]), &(Result.Daten[3])))
    Fehlermeldung ("Trapez_Parameter Trapez_einlesen", DATEIFORMATFEHLER);

  if (Lese_Fehlerstatus () == KEINFEHLER) {
    if (Result.Daten[0] != Result.Daten[1])
      Result.Steigung[0] = 1 / (Result.Daten[1] - Result.Daten[0]);
    if (Result.Daten[2] != Result.Daten[3])
      Result.Steigung[1] = 1 / (Result.Daten[2] - Result.Daten[3]);
  }
};
/********************************************************/
/*          Alles zur Klasse Konvex_Parameter           */
/********************************************************/
Konvex_Parameter::Konvex_Parameter ()
{
  /* hier muss noch was rein */

  So = Konvex;
};

Konvex_Parameter::Konvex_Parameter (const Konvex_Parameter & Das)
{
  So = Konvex;
  Daten = ((Konvex_Parameter &) Das).Daten;
};

Konvex_Parameter::Konvex_Parameter (const Parameter & Das)
{
  if (Das.So == Konvex) {
    So = Konvex;
    Daten = ((Konvex_Parameter &) Das).Daten;
  } else
    Fehlermeldung ("Konvex_Parameter::Konvex_Parameter(const Parameter& Das)", KOMPATIBILITAETSFEHLER);
};

void 
Konvex_Parameter::print ()
{
  int i;

  printf ("\nKonvex_Parameter:");

  for (i = 0; i < Daten.Lese_Dim (); i++)
    printf (" %f %f,", Daten[0][i], Daten[1][i]);
  printf ("\n");
};

void 
Konvex_Parameter::Setze_Parameter (Parameter & Damit)
{
  int i;

  for (i = 0; i < Damit.Lese_Parameter_Anzahl (); i++)
    Daten[0][i] = Damit.Lese_xParameter (i);
};

void 
Konvex_Parameter::Setze_xParameter (int Index, double Damit)
{
  Daten[0][Index] = Damit;
};

void 
Konvex_Parameter::Setze_yParameter (int Index, double Damit)
{
  Daten[1][Index] = Damit;
};

Parameter & Konvex_Parameter::operator = (Parameter & Das) {
  if (Das.So == Konvex) {
    Daten = ((Konvex_Parameter &) Das).Daten;
  } else
    Fehlermeldung ("Konvex_Parameter::operator=()", KOMPATIBILITAETSFEHLER);
};

char 
Konvex_Parameter::Finde_Parameter (DVektorArray & Daraus,
				   double Z_Grenze)
{
  int i;

  DVektorArray Buffer;

  Buffer = Daraus.Quicksort (0);
  Buffer = Buffer.Mache_Konvex (1, Z_Grenze);

  if (Buffer.Lese_Groesse () == 0)
    return (FALSE);
  Daten. ~ DVektorArray ();
  Daten.DVektorArray (Buffer.Lese_Groesse (), 2, NULL);

  for (i = 0; i < Buffer.Lese_Groesse (); i++) {
    Daten[0][i] = Buffer[i][0];
    Daten[1][i] = Buffer[i][1];
  }
  return (TRUE);
};

void 
Konvex_Parameter::Speichern (FILE * File)
{
  int i;

  fprintf (File, "%d ", Daten.Lese_Dim ());
  for (i = 0; i < Daten.Lese_Dim (); i++)
    fprintf (File, "%lf %lf ", Daten[0][i], Daten[1][i]);
  fprintf (File, "\n");
};

double 
Konvex_Parameter::Berechne_Schwerpunkt (double Hoehe, DVektor & x,
					DVektor & y)
{
  int SegmentNr;
  double Integral_fx = 0, Integral_f = 0, Result, Steigung;

  Berechne_Form (Hoehe, x, y);

  if ((Hoehe == 0) || (x.Lese_Dim () < 2)) {
    if (x.Lese_Dim () < 1)	/* ganz ganz selten */
      Result = x[0];
    else
      Result = (x[x.Lese_Dim () - 1] - x[0]) / 2;
  } else {
    for (SegmentNr = 0; SegmentNr < (x.Lese_Dim () - 1); SegmentNr++) {
      if (x[SegmentNr + 1] != x[SegmentNr]) {
	Steigung = (y[SegmentNr + 1] - y[SegmentNr]) /
	  (x[SegmentNr + 1] - x[SegmentNr]);
	Integral_fx += Steigung * ((pow (x[SegmentNr + 1], 3) - pow (x[SegmentNr], 3)) / 3) +
	  ((y[SegmentNr] - x[SegmentNr] * Steigung) * (pow (x[SegmentNr + 1], 2) - pow (x[SegmentNr], 2)) / 2);
	// (y2-y1)(x2-x1)/2 + y1(x2-x1) =
	// (x2-x1)(y2+y1)/2
	Integral_f += (y[SegmentNr] + y[SegmentNr + 1]) *
	  (x[SegmentNr + 1] - x[SegmentNr]) / 2;
      }
    }
    if (Integral_f == 0) {	/* Falls nur "Ausgabe-Peak" */
      Result = Daten[0][0];
    } else {
      Result = Integral_fx / Integral_f;
    }
  }

  return (Result);
};

double 
Konvex_Parameter::Berechne_Hoehe (double Eingabe)
{
  double Result;
  int Groesse = Daten.Lese_Dim (), i;

  if ((Eingabe >= Daten[0][0]) && (Eingabe <= Daten[0][Groesse - 1])) {		/* jetzt den Bereich suchen, in dem Eingabe liegt */
    for (i = 0; i < Groesse - 1; i++) {
      if ((Eingabe >= Daten[0][i]) && (Eingabe <= Daten[0][i + 1]))
	break;
    }

    if (Daten[0][i] != Daten[0][i + 1]) {
      if (Daten[1][i] != Daten[1][i + 1])
	Result = Daten[1][i] + (Daten[1][i + 1] - Daten[1][i]) *
	  (Eingabe - Daten[0][i]) / (Daten[0][i + 1] - Daten[0][i]);
      else
	Result = Daten[1][i];
    } else {
      if (Daten[1][i + 1] > Daten[1][i])
	Result = Daten[1][i + 1];
      else
	Result = Daten[1][i];
    }
    return (Result);
  } else
    return (0);
};

void 
Konvex_Parameter::Berechne_Form (double Hoehe, DVektor & x,
				 DVektor & y)
{
  DVektor xtemp (Daten.Lese_Dim () + 1, 0, NULL), ytemp (Daten.Lese_Dim () + 1, 0, NULL);
  int i, soviel = 0;

  /* Suchen, bis oberhalb Hoehe */
  for (i = 0; i < Daten.Lese_Dim (); i++) {
    soviel++;
    if (Daten[1][i] < Hoehe) {
      xtemp[i] = Daten[0][i];
      ytemp[i] = Daten[1][i];
    } else {
      ytemp[i] = Hoehe;
      if (i == 0) {		/* gleich der erste oberhalb */
	xtemp[0] = Daten[0][0];
      }
      /* sonst den Punkt erst noch berechnen */ 
      else {