GSL Routine als Klassenfunktion

[Klaus82]

Hi,

angestoßen durch meinen anderen Thread über Funktionen habe ich nun versucht eine GSL Routine zur Integration als Klassenfunktion zu programmieren. Das Beispiel habe ich ganz gut nachvollziehen können und habe es auch geschafft meine Funktion integrieren zu lassen.

Der Knackpunkt bei der GSL Vorlage ist wohl, dass eine Funktion f definiert wird, in der die Funktionsvorschrift steht. Später wird diese mittels F.function & f übergeben, sodass gsl_integration_qags( & F, ...) loslegen kann.

Also dachte ich mir, dass ich schlau bin und so eine Struktur f als Klassenfunktion schreibe:

double IonizationCsDistr::f(double x, void* params)
{
	double T = *(double*) params;
	double f = ionization(T,x);
	return f;
}

Dabei ist ionization (T,x) eben meine zu integrierende Funktion, die ich über den operator() in einer Klasse definiert habe.

Insgesamt sieht das ganze dann so aus - ich habe gleich versucht eine Initialisierungsliste einzubauen, hoffe das ist so richtig:

// ionization_cs_distr.h
#ifndef IONIZATION_CS_DISTR_H_
#define IONIZATION_CS_DISTR_H_

#include <math.h>
#include <gsl/gsl_integration.h>

#include "/home/klaus/Cpp/Water/Ionization/ionization_cs.h"

class IonizationCsDistr{
public:
	IonizationCsDistr(const double & I) : ionization(I), r_I(I) {}
	double operator() (double T, double W);

private:
	IonizationCS ionization;
	double f(double x, void* params);
	const double & r_I;
};

double IonizationCsDistr::operator() (double T, double W)
{
	gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000);	

	double result, error;
	double a = 0;								// lower boundary
	double b = 0.5 * (T - r_I);	// max energy
	gsl_function F;
	F.function = & f;
	F.params = & T;	

	double result_max = gsl_integration_qags (&F, a, b, 0, 1e-7, 1000, w, & result, & error);	

	gsl_integration_qags (&F, a, W, 0, 1e-7, 1000, w, & result, & error);

	return result / result_max;

}

double IonizationCsDistr::f(double x, void* params)
{
	double T = *(double*) params;
	double f = ionization(T,x);
	return f;
}

#endif
// End of File

Wenn ich diese header Datei dann aber in meine Hauptdatei include knirscht es beim kompilieren:

/home/klaus/Cpp/Water/Ionization/ionization_cs_distr.h:29: error: ISO C++ forbids taking the address
of an unqualified or parenthesized non-static member function to form a pointer to member function.
Say ‘&IonizationCsDistr::f’

/home/klaus/Cpp/Water/Ionization/ionization_cs_distr.h:29: error: cannot convert ‘double
(IonizationCsDistr::*)(double, void*)’ to ‘double (*)(double, void*)’ in assignment

Ich habe dann versucht wie wild die Klassenfunktion umzubennen oder die Übergabe F.function & f mittels der Angabe der Fehlermeldung Say &IonizationCsDistr::f.

Aber ich würde ja gerne verstehen was da knirscht. Ich schätze mal, dass es darum geht, dass die Klassenfunktion nicht Anforderungen erfüllt, welche F.function & erwartet? Also der Prototype?

Gruß,
Klaus.

[Klaus82]

Hi,

ich habe noch im Netz gesucht und bin dann wiederum hier im Forum fündig geworden: Link.

Dabei scheint die Antwort des Users auf seine Problematik folgendes zu sein:

Jetzt fang ich doch wieder an die RefCon-Pointer runterzucasten und auf den daraus hervorgehenden Objekten rumzumachen.

Ist das auch die Antwort auf meine Problematik? Immerhin sieht die Situation ähnlich aus, dass nämlich laut Fehlermeldung

‘double (*)(double, void*)’

gewünscht ist, ich allerdings mit

‘double (IonizationCsDistr::*)(double, void*)’

Ich habe also noch die Umgebung der Klasse IonizationCsDistr:: mit drin.

Und was meint der User dann mit RefCon-Pointer runtercasten ? Ich google mal...

Gruß,
Klaus.

[Klaus82]

Argh!!

Also ich jetzt im Netz gefunden, dass man solch eine Konstruktion möglichst vermeiden soll, also versuche ich jetzt eine Funktion vor der Klasse zu definieren und dieser Funktion dann über void* param die Klasse zu übergeben, doch jetzt kriege ich die Fehlermeldung

‘void*’ is not a pointer-to-object type

Ja was denn nun? Ich dachte das tolle an void wäre gerade, dass kein Datentyp vorher deklariert sein muss? Oder ist ein Datentyp jetzt wieder was anderes als ein Objekt?

Ich dachte ich könnte dann folgendes machen.

// ionization_cs_distr.h
#ifndef IONIZATION_CS_DISTR_H_
#define IONIZATION_CS_DISTR_H_

#include <math.h>
#include <gsl/gsl_integration.h>

#include "/home/klaus/Cpp/Water/Ionization/ionization_cs.h"

double f(double x, void* params)
{
	double T = params->T;
	double f = params->ionization(T,x);
	return f;
}

class IonizationCsDistr{
public:
	IonizationCsDistr(const double & I) : ionization(I), r_I(I) {}
	double operator() (double T, double W);

private:
	IonizationCS ionization;
	const double & r_I;
};

double IonizationCsDistr::operator() (double T, double W)
{
	gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000);	

	double result, error;
	double a = 0;								// lower boundary
	double b = 0.5 * (T - r_I);	// max energy
	gsl_function F;
	F.function = & f;
	F.params = & ionization;	

	double result_max = gsl_integration_qags (&F, a, b, 0, 1e-7, 1000, w, & result, & error);	

	gsl_integration_qags (&F, a, W, 0, 1e-7, 1000, w, & result, & error);

	return result / result_max;

}
#endif
// End of File

Gruß,
Klaus.

[Klaus82]

So,

ich denke ich habe das Typcasten jetzt hinbekommen, ich habe jetzt:

// ionization_cs_distr.h
#ifndef IONIZATION_CS_DISTR_H_
#define IONIZATION_CS_DISTR_H_

#include <math.h>
#include <gsl/gsl_integration.h>

#include "/home/klaus/Cpp/Water/Ionization/ionization_cs.h"

double f(double x, void* params)
{
	double T = ((IonizationCS*)params) -> T ;
	double f = ((IonizationCS*)params) -> ionization(T,x);
	return f;
}

class IonizationCsDistr{
public:
	IonizationCsDistr(const double & I) : ionization(I), r_I(I) {}
	double operator() (double T, double W);

private:
	IonizationCS ionization;
	const double & r_I;
};

double IonizationCsDistr::operator() (double T, double W)
{
	gsl_integration_workspace * w = gsl_integration_workspace_alloc (1000);	

	double result, error;
	double a = 0;								// lower boundary
	double b = 0.5 * (T - r_I);	// max energy
	gsl_function F;
	F.function = & f;
	F.params = & ionization;	

	double result_max = gsl_integration_qags (&F, a, b, 0, 1e-7, 1000, w, & result, & error);	

	gsl_integration_qags (&F, a, W, 0, 1e-7, 1000, w, & result, & error);

	return result / result_max;

}
#endif
// End of File

Jetzt bekomme ich nur die Fehlermeldung

ionization_cs_distr.h:13: error: ‘ionization’ was not declared in this scope
ionization_cs_distr.h:13: error: ‘(IonizationCS*)params’ cannot be used as a function

Jo, weil mir immer noch nicht klar ist, wie ich operator() mittels -> aurufen? Ich finde es schön, dass mir gezeigt wurde, dass ich Klassen quasie auch als Funktionen aufrufen kann, doch wenn ich einen Pointer habe, dann kann operator() nicht aufrufen, weil es keine Memberfunktion ist? Es wäre eine Art Memberoperator, aber wie kriege ich das hin?

Hier noch die ionization class

// ionization_cs.h

#ifndef IONIZATION_CS_H
#define IONIZATION_CS_H

class IonizationCS
{
public:
	IonizationCS(double ionization_potential);
	IonizationCS(double T, double ionization_potential);
	double operator() (double T);
	double operator() (double T, double W);
	double T;
private:
	const double A1, A2, B1, B2, a0, Ry, N;
	double S, F1, F2, I;
};

// constructor
inline IonizationCS::IonizationCS(double ionization_potential):
	A1(0.94), A2(1.13), B1(2.3), B2(22.0), // Liquid
//	A1(1.31), A2(0.37), B1(0.0), B2(0.0),	// Inner shells
        a0(5.3e-9),                            // Bohrradius [a0] = cm
	Ry(13.6),                              // Rydberg energy [Ry] = eV
	N(2),				       // number of electrons per shell
	I(ionization_potential)
{
	S = 4 * M_PI * N * pow(a0,2)*pow(Ry/I,2);
}

inline IonizationCS::IonizationCS(double T, double ionization_potential):
	T(T),
	A1(0.94), A2(1.13), B1(2.3), B2(22.0), // Liquid
//	A1(1.31), A2(0.37), B1(0.0), B2(0.0),	// Inner shells
	a0(5.3e-9),                             // Bohrradius [a0] = cm
	Ry(13.6),                               // Rydberg energy [Ry] = eV
	N(2),					// number of electrons per shell
	I(ionization_potential)

{
	S = 4 * M_PI * N * pow(a0,2)*pow(Ry/I,2);
}

double IonizationCS::operator() (double T)
{
	const double t = T/I;
	F1 = A1 * log(t)/(t + B1);
	F2 = A2 / (t + B2);
	return //
	1e16 * // for better plotting and fitting because of the little numbers	
	S * F1 * (t-1)/(2*pow(t,2)) * (t + 1 - (4 * sqrt(t))/(3 + t))	//
	+ S * F2 * (1 - 1/t - log(t)/(t+1));	
}

double IonizationCS::operator() (double T, double W)
{
	const double t = T/I;
	const double w = W/I;
	F1 = A1 * log(t)/(t + B1);
	F2 = A2 / (t + B2);	
	return 
	1e16 * // for better plotting and fitting because of the little numbers
	0.5 * (// according to the Rudd Paper
	S*F1*(0.5 - 1/(2*pow(w+1,2)) - 1/(2*pow(t,2)) //
	+ 1/(2*pow(t-w,2)) //
	+ 2*(t-2*w-1)/(pow(t+1,2)*sqrt(w+1)*sqrt(t-w)) //
	- 2*(t-1) / ( pow(t+1,2)*sqrt(t) ) )
	+
	S*F2*( 1 - 1/(w+1) - 1/t + 1/(t-w) //
	+ 1/(w+1) * log((t-w)/(t*(w+1))) ));
}

#endif

Ich bin doch jetzt soooo kurz davor, oder nicht?

Gruß,
Klaus.

[Klaus82]

So,

jetzt habe ich es: Derefernzieren ist das Zauberwort:

double f(double x, void* params)
{
	double T = ((IonizationCS*)params) -> T ;
	double f = (*(IonizationCS*)params)(T,x);
	return f;
}

Auch wenn ich mir jetzt gerne anhöre was für ein Spaghetti Programmierer ich bin - es funzt!

Gruß,
Klaus.

[seldon]

Wenn T eine Membervariable von IonizationCS ist, warum braucht operator() dann T als Parameter?

Falls T nur in IonizationCS steht, um die Callback-Schnittstelle befriedigen zu können, dann wäre es im Zweifel sinnvoller, etwas wie

struct IonizationCallbackData {
  double T;
  Ionization *ionization;
};

double callback_f(double x, void *context_raw) {
  IonizationCallbackData *context = static_cast<IonizationCallbackData*>(context_raw);

  return (*context->ionization)(context->T, x);
}

zu schreiben.

[Klaus82]

Hm,

ich glaube diese Schnittstelle hilft mir bei meinem Vorhaben alles in eine Klasse zu packen.

Im Folgenden will ich die Klassenfunktion

double cs(double I, double T, double W)

nach W integrieren. Und die Integration bis zur oberen Grenze W soll als Funktion in die Klassenfunktion

cs_distr(double I, double T, double W)

Also wenn ich das ganze in die GSL reinpacken möchte, dann benötige ich anhand des gegebenen Beispiels grob folgende Funktion:

double f(double W, void * params)
{
  return cs(I,T,W);
}

Wobei mir jetzt nicht ganz klar ist, wie ich über Params die beiden Variablen I und T weitergeben soll.
Ich würde mir ja gerne lokal in der Klassenfunktion eine Struktur params und eine Funktion f definieren, aber das darf ich ja nicht.
Ich könnte das auch außerhalb machen, doch benötige ich für spätere Anwendungen der GSL die Bezeichnungen noch.
Oder ich fange an das ganze params 1, params 2 und die Funktionen f,g,h ... zu nennen.

Gruß,
Klaus.

// ionization.cs - class to coup with all processes belonging to ionization
#ifndef IONIZATION_H_
#define IONIZATION_H_

#include <math.h>
#include <gsl/gsl_integration.h>

using namespace std;

class IONIZATION
{
public:
	IONIZATION();
	double cs(double I, double T, double W);
	double cs_distr(double I, double T, double W);
private:
	const double A1,A2,B1,B2,a0,Ry,N;
	gsl_function F;
//	F.function & f;
//	F.params = &

};

IONIZATION::IONIZATION() : 
A1(0.94), A2(1.13), B1(2.3), B2(22.0),a0(5.3e-9),
Ry(13.6),N(2) {}

double IONIZATION::cs(double I, double T, double W)
{
	if(T < I)
	{
		cout << "IONIZATION IMPOSSIBLE!" << endl;
		return 0;
	}
	else
	{
		const double t = T/I;
		const double w = W/I;
		double S = 4 * M_PI * N * pow(a0,2)*pow(Ry/I,2);		
		double F1 = A1 * log(t)/(t + B1);
		double F2 = A2 / (t + B2);	
		return 
		1e16 * // for better plotting and fitting because of the little numbers
		0.5 * (// according to the Rudd Paper
		S*F1*(0.5 - 1/(2*pow(w+1,2)) - 1/(2*pow(t,2)) //
		+ 1/(2*pow(t-w,2)) //
		+ 2*(t-2*w-1)/(pow(t+1,2)*sqrt(w+1)*sqrt(t-w)) //
		- 2*(t-1) / ( pow(t+1,2)*sqrt(t) ) )
		+
		S*F2*( 1 - 1/(w+1) - 1/t + 1/(t-w) //
		+ 1/(w+1) * log((t-w)/(t*(w+1))) ));
	}
}

double IONIZATION::cs_distr(double I, double T, double W)
{

}

#endif

[Klaus82]

Also ich denke meine Funktion müsste ungefähr so aussehen:

double f(double x, void * params)
{
	return ((IONIZATION*)params)->cs(I,T,x);
}

Nur wie kriege ich dann die Information über I und T da rein?

Ist das der Grund für diese Callback Schnittstelle?

Gruß,
Klaus.

[Klaus82]

Also ich habe mir jetzt eine Struktur gebastelt als Vermittler, aber ich verstehe den Fehler der Kompilierung nicht - ich denke es hängt irgendwie mit der Vorwärtsdeklaration zusammen.

Hier der Codeschnippsel:

// ionization.cs - class to coup with all processes belonging to ionization
#ifndef IONIZATION_H_
#define IONIZATION_H_

using namespace std;

class IONIZATION;

struct NECESSARY_EVIL;

double f(double x, void * params);

class IONIZATION
{
public:
	double cs_distr(double I, double T, double W);
// ..
};

double IONIZATION::cs_distr(double I, double T, double W)
{
	NECESSARY_EVIL ne(I,T,this);
  // ...

}

struct NECESSARY_EVIL
{
	NECESSARY_EVIL(double I, double T, IONIZATION * p_io) : I(I), T(T), p_io(p_io) {}
	const double I,T;
	IONIZATION * p_io;
};

double f(double x, void * params)
{
	double I = ((NECESSARY_EVIL*)params)->I;
	double T = ((NECESSARY_EVIL*)params)->T;
	return ((NECESSARY_EVIL*)params)->p_io->cs(I,T,x);
}

#endif

Und die Fehlermeldung ist jetzt eben:

Header_Files/ionization.h:118: error: variable ‘NECESSARY_EVIL ne’ has initializer but incomplete type

Aber rein vom Syntax ist die doch sauber initialisiert, oder nicht?

Gruß,
Klaus.