Added option to specify a filter function for the target objective.

FEBioOpt/FEObjectiveFunction.cpp

@@ -109,7 +109,9 @@ double FEObjectiveFunction::RegressionCoefficient(const std::vector<double>& y0,
 	xyb /= ndata;
 	x2b /= ndata; y2b /= ndata;
 
-	double rsq = pow(xyb - xb * yb, 2) / (x2b - xb * xb) / (y2b - yb * yb);
+	double D = (x2b - xb * xb) * (y2b - yb * yb);
+
+	double rsq = (D != 0 ? pow(xyb - xb * yb, 2) / D : 0);
 	return rsq;
 }
 
@@ -207,42 +209,74 @@ void FEDataFitObjective::EvaluateFunctions(vector<double>& f)
 }
 
 //=============================================================================
-FEMinimizeObjective::FEMinimizeObjective(FEModel* fem) : FEObjectiveFunction(fem)
+
+bool FEMinimizeObjective::ParamFunction::Init()
 {
+	if (!Function::Init()) return false;
+
+	FEParamValue val = m_fem->GetParameterValue(ParamString(m_name.c_str()));
+	if (val.isValid() == false) return false;
+	if (val.type() != FE_PARAM_DOUBLE) return false;
+	m_var = (double*)val.data_ptr();
+	if (m_var == nullptr) return false;
+	return true;
 }
 
-bool FEMinimizeObjective::AddFunction(const char* szname, double targetValue)
+FEMinimizeObjective::FilterAvgFunction::FilterAvgFunction(FEModel* fem, FELogElemData* pd, FEElementSet* elemSet, double trg) : Function(fem)
 {
-	// make sure we have a model
-	FEModel* fem = GetFEModel();
-	if (fem == 0) return false;
+	m_pd = pd;
+	m_elemSet = elemSet;
+	m_y0 = trg;
+}
+
+bool FEMinimizeObjective::FilterAvgFunction::Init()
+{
+	if (!Function::Init()) return false;
+	if (m_pd == nullptr) return false;
+	if (m_elemSet == nullptr) return false;
+	return true;
+}
+
+double FEMinimizeObjective::FilterAvgFunction::Value() const
+{
+	int NE = m_elemSet->Elements();
+	double sum = 0.0;
+	for (int i = 0; i < NE; ++i)
+	{
+		sum += m_pd->value(m_elemSet->Element(i));
+	}
+	sum /= (double)NE;
 
-	// create a new function
-	Function fnc;
-	fnc.name = szname;
-	fnc.y0 = targetValue;
+	return sum;
+}
 
-	m_Func.push_back(fnc);
+FEMinimizeObjective::FEMinimizeObjective(FEModel* fem) : FEObjectiveFunction(fem)
+{
+}
 
-	return true;
+FEMinimizeObjective::~FEMinimizeObjective()
+{
+	for (Function* f : m_Func) delete f;
+	m_Func.clear();
+}
+
+void FEMinimizeObjective::AddFunction(FEMinimizeObjective::Function* func)
+{
+	m_Func.push_back(func);
 }
 
 bool FEMinimizeObjective::Init()
 {
 	if (FEObjectiveFunction::Init() == false) return false;
 
 	FEModel* fem = GetFEModel();
-	if (fem == 0) return false;
+	if (fem == nullptr) return false;
 
 	int N = (int) m_Func.size();
 	for (int i=0; i<N; ++i)
 	{
-		Function& Fi = m_Func[i];
-		FEParamValue val = fem->GetParameterValue(ParamString(Fi.name.c_str()));
-		if (val.isValid() == false) return false;
-		if (val.type() != FE_PARAM_DOUBLE) return false;
-		Fi.var = (double*) val.data_ptr();
-		if (Fi.var == 0) return false;
+		Function& Fi = *m_Func[i];
+		if (Fi.Init() == false) return false;
 	}
 
 	return true;
@@ -253,16 +287,14 @@ int FEMinimizeObjective::Measurements()
 	return (int) m_Func.size();
 }
 
-
 void FEMinimizeObjective::EvaluateFunctions(vector<double>& f)
 {
 	int N = (int)m_Func.size();
 	f.resize(N);
 	for (int i=0; i<N; ++i)
 	{
-		Function& Fi = m_Func[i];
-		if (Fi.var) f[i] = *Fi.var;
-		else f[i] = 0;
+		Function& Fi = *m_Func[i];
+		f[i] = Fi.Value();
 	}
 }
 
@@ -272,7 +304,7 @@ void FEMinimizeObjective::GetMeasurements(vector<double>& y)
 	y.resize(N);
 	for (int i=0; i<N; ++i)
 	{
-		y[i] = m_Func[i].y0;
+		y[i] = m_Func[i]->Target();
 	}
 }
 

FEBioOpt/FEObjectiveFunction.h

@@ -148,24 +148,58 @@ class FEMinimizeObjective : public FEObjectiveFunction
 	class Function
 	{
 	public:
-		string		name;
-		double*		var;
+		Function(FEModel* fem) : m_fem(fem) {}
+		virtual ~Function() {}
 
-		double	y0;		// target value (i.e. "measurment")
+		virtual bool Init() { return (m_fem != nullptr); }
+
+		virtual double Target() const = 0;
+		virtual double Value() const = 0;
+
+	protected:
+		FEModel* m_fem;
+	};
+
+public:
+	class ParamFunction : public Function
+	{
+	public:
+		string	m_name;
+		double* m_var = nullptr;
+
+		double	m_y0 = 0;		// target value (i.e. "measurment")
+
+		ParamFunction(FEModel* fem, const string& name, double trg) : Function(fem), m_name(name), m_y0(trg) {}
+
+	public:
+		bool Init() override;
+		double Target() const override { return m_y0; }
+		double Value() const override { return *m_var; }
+	};
+
+	class FilterAvgFunction : public Function
+	{
+	public:
+		FELogElemData* m_pd = nullptr;
+		FEElementSet* m_elemSet = nullptr;
+		double	m_y0 = 0;
+
+		FilterAvgFunction(FEModel* fem, FELogElemData* pd, FEElementSet* elemSet, double trg);
 
 	public:
-		Function() : var(0), y0(0.0) {}
-		Function(const Function& f) { name = f.name; var = f.var; y0 = f.y0; }
-		void operator = (const Function& f) { name = f.name; var = f.var; y0 = f.y0; }
+		bool Init() override;
+		double Target() const override { return m_y0; }
+		double Value() const override;
 	};
 
 public:
 	FEMinimizeObjective(FEModel* fem);
+	~FEMinimizeObjective();
 
 	// one-time initialization
 	bool Init() override;
 
-	bool AddFunction(const char* szname, double targetValue);
+	void AddFunction(Function* func);
 
 public:
 	// return number of measurements
@@ -178,7 +212,7 @@ class FEMinimizeObjective : public FEObjectiveFunction
 	void GetMeasurements(std::vector<double>& y) override;
 
 private:
-	std::vector<Function>	m_Func;
+	std::vector<Function*>	m_Func;
 };
 
 //=============================================================================