add back missing parametrizations

Huangzizhou · Oct 10, 2024 · b33ee89 · b33ee89
1 parent f0369b1
commit b33ee89
Show file tree

Hide file tree

Showing 4 changed files with 171 additions and 1 deletion.
diff --git a/json-specs/opt-input-spec.json b/json-specs/opt-input-spec.json
@@ -335,9 +335,31 @@
             "sdf-to-mesh",
             "periodic-mesh-tile",
             "mesh-affine",
-            "matrix-product"
+            "matrix-product",
+            "box-constraint-reparametrization"
         ]
     },
+    {
+        "pointer": "/variable_to_simulation/*/composition/*",
+        "type_name": "box-constraint-reparametrization",
+        "type": "object",
+        "required": [
+            "bounds"
+        ],
+        "doc": "Reparametrize variable to fit into box constraints."
+    },
+    {
+        "pointer": "/variable_to_simulation/*/composition/*/bounds",
+        "type": "list"
+    },
+    {
+        "pointer": "/variable_to_simulation/*/composition/*/bounds/*",
+        "type": "list"
+    },
+    {
+        "pointer": "/variable_to_simulation/*/composition/*/bounds/*/*",
+        "type": "float"
+    },
     {
         "pointer": "/variable_to_simulation/*/composition/*",
         "type_name": "matrix-product",
@@ -641,6 +663,28 @@
         "default": 0,
         "type": "float"
     },
+    {
+        "pointer": "/variable_to_simulation/*/composition/*",
+        "type_name": "periodic-mesh-tile",
+        "type": "object",
+        "required": [
+            "dimensions"
+        ],
+        "doc": "Append a list of constants at the end of the input vector"
+    },
+    {
+        "pointer": "/variable_to_simulation/*/composition/*",
+        "type_name": "mesh-affine",
+        "type": "object",
+        "required": [
+            "type",
+            "dimension"
+        ],
+        "optional": [
+            "transformation"
+        ],
+        "doc": "Append a list of constants at the end of the input vector"
+    },
     {
         "pointer": "/variable_to_simulation/*/composition/*/dimensions",
         "type": "list"

diff --git a/src/polyfem/solver/Optimizations.cpp b/src/polyfem/solver/Optimizations.cpp
@@ -484,6 +484,14 @@ namespace polyfem::solver
 		{
 			map = std::make_shared<InsertConstantMap>(args["size"], args["value"], args["start"]);
 		}
+		else if (type == "matrix-product")
+		{
+			map = std::make_shared<MatrixProduct>(args);
+		}
+		else if (type == "box-constraint-reparametrization")
+		{
+			map = std::make_shared<BoxConstraintReparametrization>(args);
+		}
 		else if (type == "linear-filter")
 		{
 			map = std::make_shared<LinearFilter>(*(states[args["state"]]->mesh), args["radius"]);

diff --git a/src/polyfem/solver/forms/parametrization/Parametrizations.cpp b/src/polyfem/solver/forms/parametrization/Parametrizations.cpp
@@ -3,6 +3,7 @@
 #include <polyfem/basis/ElementBases.hpp>
 #include <polyfem/utils/Logger.hpp>
 #include <polyfem/Common.hpp>
+#include <polyfem/utils/JSONUtils.hpp>
 #include <polyfem/utils/ElasticityUtils.hpp>
 
 namespace polyfem::solver
@@ -493,4 +494,92 @@ namespace polyfem::solver
 		return hess.transpose() * grad;
 	}
 
+	MatrixProduct::MatrixProduct(const json &args)
+	{
+		int rows = args["matrix"].size();
+		assert (args["matrix"][0].is_array());
+		int cols = args["matrix"][0].size();
+
+		mat.setZero(rows, cols);
+		Eigen::VectorXd tmp;
+		for (int i = 0; i < rows; i++)
+		{
+			tmp = args["matrix"][i];
+			assert(tmp.size() == cols);
+			mat.row(i) = tmp;
+		}
+
+		left_mult = args["left_multiply"];
+	}
+	int MatrixProduct::size(const int x_size) const
+	{
+		if (left_mult)
+			return (x_size * mat.rows()) / mat.cols();
+		else
+			return (x_size * mat.cols()) / mat.rows();
+	}
+	Eigen::VectorXd MatrixProduct::eval(const Eigen::VectorXd &x) const
+	{
+		if (left_mult)
+		{
+			assert(x.size() % mat.cols() == 0);
+			Eigen::MatrixXd input = utils::unflatten(x, x.size() / mat.cols());
+			return utils::flatten(mat * input);
+		}
+		else
+		{
+			assert(x.size() % mat.rows() == 0);
+			Eigen::MatrixXd input = utils::unflatten(x, mat.rows());
+			return utils::flatten(input * mat);
+		}
+	}
+	Eigen::VectorXd MatrixProduct::apply_jacobian(const Eigen::VectorXd &grad, const Eigen::VectorXd &x) const
+	{
+		if (left_mult)
+		{
+			Eigen::MatrixXd grad_mat = utils::unflatten(grad, x.size() / mat.cols());
+			return utils::flatten(mat.transpose() * grad_mat);
+		}
+		else
+		{
+			Eigen::MatrixXd grad_mat = utils::unflatten(grad, mat.cols());
+			return utils::flatten(grad_mat * mat.transpose());
+		}
+	}
+
+	BoxConstraintReparametrization::BoxConstraintReparametrization(const json &args)
+	{
+		bounds = args["bounds"];
+		if (bounds.rows() != 2)
+			log_and_throw_adjoint_error("Invalid bounds size!");
+	}
+	int BoxConstraintReparametrization::size(const int x_size) const
+	{
+		return x_size;
+	}
+	Eigen::VectorXd BoxConstraintReparametrization::eval(const Eigen::VectorXd &x) const
+	{
+		Eigen::VectorXd y = x;
+		if (bounds.cols() != x.size())
+			log_and_throw_adjoint_error("Inconsistent input size, {} vs {}!", x.size(), bounds.cols());
+		for (int i = 0; i < x.size(); i++)
+			y(i) = logistic(x(i), bounds(0, i), bounds(1, i));
+		return y;
+	}
+	Eigen::VectorXd BoxConstraintReparametrization::apply_jacobian(const Eigen::VectorXd &grad, const Eigen::VectorXd &x) const
+	{
+		Eigen::VectorXd gradv = grad;
+		for (int i = 0; i < x.size(); i++)
+			gradv(i) *= logistic_grad(x(i), bounds(0, i), bounds(1, i));
+		return gradv;
+	}
+	double BoxConstraintReparametrization::logistic(double x, double lower, double higher)
+	{
+		return lower + (higher - lower) / (1. + exp(-x));
+	}
+	double BoxConstraintReparametrization::logistic_grad(double x, double lower, double higher)
+	{
+		const double e = exp(-x);
+		return (higher - lower) * e / (1 + e) / (1 + e);
+	}
 } // namespace polyfem::solver
diff --git a/src/polyfem/solver/forms/parametrization/Parametrizations.hpp b/src/polyfem/solver/forms/parametrization/Parametrizations.hpp
@@ -182,4 +182,33 @@ namespace polyfem::solver
 		const double start_val_;
 		const double dt_;
 	};
+
+	class MatrixProduct : public Parametrization
+	{
+	public:
+		MatrixProduct(const json &args);
+
+		int size(const int x_size) const override;
+		Eigen::VectorXd eval(const Eigen::VectorXd &x) const override;
+		Eigen::VectorXd apply_jacobian(const Eigen::VectorXd &grad, const Eigen::VectorXd &x) const override;
+
+	private:
+		Eigen::MatrixXd mat;
+		bool left_mult;
+	};
+
+	class BoxConstraintReparametrization : public Parametrization
+	{
+	public:
+		BoxConstraintReparametrization(const json &args);
+
+		int size(const int x_size) const override;
+		Eigen::VectorXd eval(const Eigen::VectorXd &x) const override;
+		Eigen::VectorXd apply_jacobian(const Eigen::VectorXd &grad, const Eigen::VectorXd &x) const override;
+	private:
+		static double logistic(double x, double lower, double higher);
+		static double logistic_grad(double x, double lower, double higher);
+
+		Eigen::MatrixXd bounds;
+	};
 } // namespace polyfem::solver