C++ MyParser::EvalRemoveSingularity方法代码示例

void NonLinearFit::calculateFitCurveData(double *X, double *Y)
{
	MyParser parser;
	for (int i=0; i<d_p; i++)
		parser.DefineVar(d_param_names[i].ascii(), &d_results[i]);

	QMapIterator<QString, double> i(d_constants);
 	while (i.hasNext()) {
     	i.next();
		parser.DefineConst(i.key().ascii(), i.value());
 	}

	double x;
	parser.DefineVar("x", &x);
	parser.SetExpr(d_formula.ascii());

	if (d_gen_function){
		double X0 = d_x[0];
		double step = (d_x[d_n - 1] - X0)/(d_points - 1);
		for (int i=0; i<d_points; i++){
		    x = X0 + i*step;
			X[i] = x;
			Y[i] = parser.EvalRemoveSingularity(&x, false);
		}
	} else {
		for (int i=0; i<d_points; i++) {
		    x = d_x[i];
			X[i] = x;
			Y[i] = parser.EvalRemoveSingularity(&x, false);
		}
	}
}

void NonLinearFit::removeDataSingularities()
{
	MyParser parser;
	for (int i = 0; i < d_p; i++){
		double param = gsl_vector_get(d_param_init, i);
		parser.DefineVar(d_param_names[i].ascii(), &param);
	}

	QMapIterator<QString, double> it(d_constants);
 	while (it.hasNext()) {
     	it.next();
		parser.DefineConst(it.key().ascii(), it.value());
 	}

	double xvar;
	parser.DefineVar("x", &xvar);
	parser.SetExpr(d_formula.ascii());

    for (int i = 0; i < d_n; i++){
    	xvar = d_x[i];
    	try {
			parser.EvalRemoveSingularity(&xvar);
    	} catch(MyParser::Pole){
			QApplication::restoreOverrideCursor();
			QMessageBox::critical((ApplicationWindow *)parent(), QObject::tr("QtiPlot"),
			QObject::tr("Ignored data point at x = %1.").arg(xvar));

    		removePole(i);
    	}
    }
}

double NonLinearFit::eval(double *par, double x)
{
	MyParser parser;
	for (int i=0; i<d_p; i++)
		parser.DefineVar(d_param_names[i].ascii(), &par[i]);

	QMapIterator<QString, double> i(d_constants);
 	while (i.hasNext()) {
     	i.next();
		parser.DefineConst(i.key().ascii(), i.value());
 	}

	parser.DefineVar("x", &x);
	parser.SetExpr(d_formula.ascii());
    return parser.EvalRemoveSingularity(&x, false);
}

double user_d(const gsl_vector * x, void *params) {
    int n = ((struct FitData *)params)->n;
    int p = ((struct FitData *)params)->p;
    double *X = ((struct FitData *)params)->X;
    double *Y = ((struct FitData *)params)->Y;
    double *sigma = ((struct FitData *)params)->sigma;

	NonLinearFit *fitter = (NonLinearFit *)((struct FitData *) params)->fitter;
	const char *function = fitter->formula().ascii();
	QStringList parNames = fitter->parameterNames();

    double val=0;
    MyParser parser;
    try {
        double *parameters = new double[p];
        double xvar;
        parser.DefineVar("x", &xvar);
        for (int i=0; i < p; i++) {
            parameters[i]=gsl_vector_get(x,i);
            parser.DefineVar(parNames[i].ascii(), &parameters[i]);
        }

		QMapIterator<QString, double> i(fitter->constants());
 		while (i.hasNext()){
     		i.next();
			parser.DefineConst(i.key().ascii(), i.value());
 		}

        parser.SetExpr(function);
        for (int j = 0; j < n; j++) {
            xvar = X[j];
			double s = 1.0/sqrt(sigma[j]);
            try {
				double t = (parser.EvalRemoveSingularity(&xvar) - Y[j])/s;
				val += t*t;
			} catch (MyParser::Pole) {
				return GSL_POSINF; //weird, I know. blame gsl.
			}
        }
        delete[] parameters;
    } catch (mu::ParserError &e) {
        QMessageBox::critical(0,"QtiPlot - Input function error",QString::fromStdString(e.GetMsg()));
        return GSL_EINVAL;
    }
    return val;
}

bool FunctionCurve::loadData(int points, bool xLog10Scale)
{
    if (!points)
        points = dataSize();

	double *X = (double *)malloc(points*sizeof(double));
	if (!X){
		QMessageBox::critical(0, QObject::tr("QtiPlot - Memory Allocation Error"),
		QObject::tr("Not enough memory, operation aborted!"));
		return false;
	}
	double *Y = (double *)malloc(points*sizeof(double));
	if (!Y){
		QMessageBox::critical(0, QObject::tr("QtiPlot - Memory Allocation Error"),
		QObject::tr("Not enough memory, operation aborted!"));
		free(X);
		return false;
	}

	double step = (d_to - d_from)/(double)(points - 1.0);
	if (d_function_type == Normal){
		MyParser parser;
		double x = d_from;
		try {
			parser.DefineVar(d_variable.ascii(), &x);
			QMapIterator<QString, double> i(d_constants);
			while (i.hasNext()){
				i.next();
				parser.DefineConst(i.key().ascii(), i.value());
			}
			parser.SetExpr(d_formulas[0].ascii());

			int lastButOne = points - 1;
			try {
				double xl = x, xr;
				double y = parser.EvalRemoveSingularity(&x, false);
				bool wellDefinedFunction = true;
				if (!gsl_finite(y)){// try to find a first well defined point (might help for some not really bad functions)
					wellDefinedFunction = false;
					for (int i = 0; i < lastButOne; i++){
						xl = x;
						x += step;
						xr = x;
						y = parser.Eval();
						if (gsl_finite(y)){
							wellDefinedFunction = true;
							int iter = 0;
							double x0 = x, y0 = y;
							while(fabs(xr - xl)/step > 1e-15 && iter < points){
								x = 0.5*(xl + xr);
								y = parser.Eval();
								if (gsl_finite(y)){
									xr = x;
									x0 = x;
									y0 = y;
								} else
									xl = x;
								iter++;
							}
							d_from = x0;
							X[0] = x0;
							Y[0] = y0;
							step = (d_to - d_from)/(double)(lastButOne);
							break;
						}
					}
					if (!wellDefinedFunction){
						QMessageBox::critical(0, QObject::tr("QtiPlot"),
						QObject::tr("The function %1 is not defined in the specified interval!").arg(d_formulas[0]));
						free(X); free(Y);
						return false;
					}
				} else {
					X[0] = d_from;
					Y[0] = y;
				}
			} catch (MyParser::Pole) {}

			ScaleEngine *sc_engine = 0;
			if (plot())
				sc_engine = (ScaleEngine *)plot()->axisScaleEngine(xAxis());

			if (xLog10Scale || (d_from > 0 && d_to > 0 && sc_engine &&
				sc_engine->type() == ScaleTransformation::Log10)){
				step = log10(d_to/d_from)/(double)(points - 1);
				for (int i = 1; i < lastButOne; i++ ){
					x = d_from*pow(10, i*step);
					X[i] = x;
					try {
						Y[i] = parser.EvalRemoveSingularity(&x, false);
					} catch (MyParser::Pole){}
				}
			} else {
				for (int i = 1; i < lastButOne; i++ ){
					x += step;
					X[i] = x;
					try {
						Y[i] = parser.EvalRemoveSingularity(&x, false);
					} catch (MyParser::Pole){}
				}
//.........这里部分代码省略.........