Implementation of ContinuousPeriodic1DFunction: Common platform

platforms/common/src/ExpressionUtilities.cpp

@@ -70,7 +70,7 @@ void ExpressionUtilities::processExpression(stringstream& out, const ExpressionT
     string name = prefix+context.intToString(temps.size());
     bool hasRecordedNode = false;
     bool isVecType = (tempType[tempType.size()-1] == '3');
-    
+
     out << tempType << " " << name << " = ";
     switch (node.getOperation().getId()) {
         case Operation::CONSTANT:
@@ -215,7 +215,7 @@ void ExpressionUtilities::processExpression(stringstream& out, const ExpressionT
             }
             else {
                 // This is a tabulated function.
-                
+
                 int i;
                 for (i = 0; i < (int) functionNames.size() && functionNames[i].first != node.getOperation().getName(); i++)
                     ;
@@ -255,6 +255,24 @@ void ExpressionUtilities::processExpression(stringstream& out, const ExpressionT
                         }
                         out << "}\n";
                     }
+                    else if (dynamic_cast<const ContinuousPeriodic1DFunction*>(functions[i]) != NULL) {
+                        out << "real x = " << getTempName(node.getChildren()[0], temps) << suffix << ";\n";
+                        out << "if (x >= " << paramsFloat[0] << " && x <= " << paramsFloat[1] << ") {\n";
+                        out << "x = (x - " << paramsFloat[0] << ")*" << paramsFloat[2] << ";\n";
+                        out << "int index = (int) (floor(x));\n";
+                        out << "index = min(index, (int) " << paramsInt[3] << ");\n";
+                        out << "float4 coeff = " << functionNames[i].second << "[index];\n";
+                        out << "real b = x-index;\n";
+                        out << "real a = 1.0f-b;\n";
+                        for (int j = 0; j < nodes.size(); j++) {
+                            const vector<int>& derivOrder = dynamic_cast<const Operation::Custom*>(&nodes[j]->getOperation())->getDerivOrder();
+                            if (derivOrder[0] == 0)
+                                out << nodeNames[j] << suffix << " = a*coeff.x+b*coeff.y+((a*a*a-a)*coeff.z+(b*b*b-b)*coeff.w)/(" << paramsFloat[2] << "*" << paramsFloat[2] << ");\n";
+                            else
+                                out << nodeNames[j] << suffix << " = (coeff.y-coeff.x)*" << paramsFloat[2] << "+((1.0f-3.0f*a*a)*coeff.z+(3.0f*b*b-1.0f)*coeff.w)/" << paramsFloat[2] << ";\n";
+                        }
+                        out << "}\n";
+                    }
                     else if (dynamic_cast<const Continuous2DFunction*>(functions[i]) != NULL) {
                         out << "real x = " << getTempName(node.getChildren()[0], temps) << suffix << ";\n";
                         out << "real y = " << getTempName(node.getChildren()[1], temps) << suffix << ";\n";
@@ -446,7 +464,7 @@ void ExpressionUtilities::processExpression(stringstream& out, const ExpressionT
             out << "-" << getTempName(node.getChildren()[0], temps);
             break;
         case Operation::SQRT:
-            callFunction(out, "sqrtf", "sqrt", getTempName(node.getChildren()[0], temps), tempType); 
+            callFunction(out, "sqrtf", "sqrt", getTempName(node.getChildren()[0], temps), tempType);
             break;
         case Operation::EXP:
             callFunction(out, "expf", "exp", getTempName(node.getChildren()[0], temps), tempType);
@@ -675,19 +693,19 @@ void ExpressionUtilities::findRelatedCustomFunctions(const ExpressionTreeNode& n
             vector<const Lepton::ExpressionTreeNode*>& nodes) {
     if (searchNode.getOperation().getId() == Operation::CUSTOM && node.getOperation().getName() == searchNode.getOperation().getName()) {
         // Make sure the arguments are identical.
-        
+
         for (int i = 0; i < (int) node.getChildren().size(); i++)
             if (node.getChildren()[i] != searchNode.getChildren()[i])
                 return;
-        
+
         // See if we already have an identical node.
-        
+
         for (int i = 0; i < (int) nodes.size(); i++)
             if (*nodes[i] == searchNode)
                 return;
-        
+
         // Add the node.
-        
+
         nodes.push_back(&searchNode);
     }
     else
@@ -735,6 +753,28 @@ vector<float> ExpressionUtilities::computeFunctionCoefficients(const TabulatedFu
         width = 4;
         return f;
     }
+    if (dynamic_cast<const ContinuousPeriodic1DFunction*>(&function) != NULL) {
+        // Compute the spline coefficients.
+
+        const ContinuousPeriodic1DFunction& fn = dynamic_cast<const ContinuousPeriodic1DFunction&>(function);
+        vector<double> values;
+        double min, max;
+        fn.getFunctionParameters(values, min, max);
+        int numValues = values.size();
+        vector<double> x(numValues), derivs;
+        for (int i = 0; i < numValues; i++)
+            x[i] = min+i*(max-min)/(numValues-1);
+        SplineFitter::createPeriodicSpline(x, values, derivs);
+        vector<float> f(4*(numValues-1));
+        for (int i = 0; i < (int) values.size()-1; i++) {
+            f[4*i] = (float) values[i];
+            f[4*i+1] = (float) values[i+1];
+            f[4*i+2] = (float) (derivs[i]/6.0);
+            f[4*i+3] = (float) (derivs[i+1]/6.0);
+        }
+        width = 4;
+        return f;
+    }
     if (dynamic_cast<const Continuous2DFunction*>(&function) != NULL) {
         // Compute the spline coefficients.
 
@@ -785,7 +825,7 @@ vector<float> ExpressionUtilities::computeFunctionCoefficients(const TabulatedFu
     }
     if (dynamic_cast<const Discrete1DFunction*>(&function) != NULL) {
         // Record the tabulated values.
-        
+
         const Discrete1DFunction& fn = dynamic_cast<const Discrete1DFunction&>(function);
         vector<double> values;
         fn.getFunctionParameters(values);
@@ -798,7 +838,7 @@ vector<float> ExpressionUtilities::computeFunctionCoefficients(const TabulatedFu
     }
     if (dynamic_cast<const Discrete2DFunction*>(&function) != NULL) {
         // Record the tabulated values.
-        
+
         const Discrete2DFunction& fn = dynamic_cast<const Discrete2DFunction&>(function);
         int xsize, ysize;
         vector<double> values;
@@ -812,7 +852,7 @@ vector<float> ExpressionUtilities::computeFunctionCoefficients(const TabulatedFu
     }
     if (dynamic_cast<const Discrete3DFunction*>(&function) != NULL) {
         // Record the tabulated values.
-        
+
         const Discrete3DFunction& fn = dynamic_cast<const Discrete3DFunction&>(function);
         int xsize, ysize, zsize;
         vector<double> values;
@@ -840,6 +880,16 @@ vector<vector<double> > ExpressionUtilities::computeFunctionParameters(const vec
             params[i].push_back((values.size()-1)/(max-min));
             params[i].push_back(values.size()-2);
         }
+        else if (dynamic_cast<const ContinuousPeriodic1DFunction*>(functions[i]) != NULL) {
+            const ContinuousPeriodic1DFunction& fn = dynamic_cast<const ContinuousPeriodic1DFunction&>(*functions[i]);
+            vector<double> values;
+            double min, max;
+            fn.getFunctionParameters(values, min, max);
+            params[i].push_back(min);
+            params[i].push_back(max);
+            params[i].push_back((values.size()-1)/(max-min));
+            params[i].push_back(values.size()-2);
+        }
         else if (dynamic_cast<const Continuous2DFunction*>(functions[i]) != NULL) {
             const Continuous2DFunction& fn = dynamic_cast<const Continuous2DFunction&>(*functions[i]);
             vector<double> values;
@@ -906,6 +956,8 @@ vector<vector<double> > ExpressionUtilities::computeFunctionParameters(const vec
 Lepton::CustomFunction* ExpressionUtilities::getFunctionPlaceholder(const TabulatedFunction& function) {
     if (dynamic_cast<const Continuous1DFunction*>(&function) != NULL)
         return &fp1;
+    if (dynamic_cast<const ContinuousPeriodic1DFunction*>(&function) != NULL)
+        return &fp1;
     if (dynamic_cast<const Continuous2DFunction*>(&function) != NULL)
         return &fp2;
     if (dynamic_cast<const Continuous3DFunction*>(&function) != NULL)