From 320b961ab1fa544129afaf4e3e78a614240110ef Mon Sep 17 00:00:00 2001 From: Dingcheng Date: Thu, 3 Aug 2023 14:36:58 -0500 Subject: [PATCH 1/9] changed function names costGrad to grad and costHessian to hessian --- soupy/modeling/controlCostFunctional.py | 22 +++++++++---------- soupy/modeling/controlCostHessian.py | 2 +- soupy/modeling/meanVarRiskMeasure.py | 4 ++-- soupy/modeling/meanVarRiskMeasureSAA.py | 8 +++---- soupy/modeling/riskMeasure.py | 6 ++--- soupy/modeling/superquantileRiskMeasureSAA.py | 4 ++-- soupy/optimization/bfgs.py | 2 +- soupy/optimization/inexactNewtonCG.py | 8 +++---- soupy/optimization/sgd.py | 4 ++-- soupy/optimization/steepestDescent.py | 4 ++-- soupy/test/ptest_meanVarRiskMeasureSAA_MPI.py | 6 ++--- soupy/test/ptest_superquantileSAA_MPI.py | 6 ++--- soupy/test/test_bfgs.py | 4 ++-- soupy/test/test_controlCostFunctional.py | 6 ++--- soupy/test/test_inexactNewtonCG.py | 4 ++-- soupy/test/test_meanVarRiskMeasure.py | 2 +- soupy/test/test_meanVarRiskMeasureSAA.py | 6 ++--- .../test_meanVarRiskMeasureSAACostFunction.py | 8 +++---- soupy/test/test_riskMeasureCostFunctional.py | 2 +- soupy/test/test_sgd.py | 4 ++-- soupy/test/test_steepestDescent.py | 4 ++-- soupy/test/test_superquantileSAA.py | 6 ++--- soupy/utils/scipyCostWrapper.py | 14 ++++++------ soupy/utils/stochasticCostFiniteDifference.py | 4 ++-- 24 files changed, 70 insertions(+), 70 deletions(-) diff --git a/soupy/modeling/controlCostFunctional.py b/soupy/modeling/controlCostFunctional.py index 52292dc..1cb3c0b 100644 --- a/soupy/modeling/controlCostFunctional.py +++ b/soupy/modeling/controlCostFunctional.py @@ -28,21 +28,21 @@ def cost(self, z, order=0): Given the control variable z evaluate the cost functional. Order specifies \ the order of derivatives required after the computation of the cost """ - raise NotImplementedError("Child class should implement method costValue") + raise NotImplementedError("Child class should implement method cost") - def costGrad(self, g): + def grad(self, g): """ Evaluate the gradient of the cost functional. Assumes :code:`cost` is called \ with order >=1 """ - raise NotImplementedError("Child class should implement method costGrad") + raise NotImplementedError("Child class should implement method grad") - def costHessian(self, zhat, Hzhat): + def hessian(self, zhat, Hzhat): """ Evaluate the Hessian of the cost functional acting in direction :code:`zhat`. \ Assumes :code:`cost` is called with order >=2 """ - raise NotImplementedError("Child class should implement method costHessian") + raise NotImplementedError("Child class should implement method hessian") class DeterministicControlCostFunctional(ControlCostFunctional): @@ -167,7 +167,7 @@ def cost(self, z, order=0, **kwargs): penalization = self.penalization.cost(z) return objective + penalization - def costGrad(self, g): + def grad(self, g): """ Computes the gradient of the cost functional @@ -186,7 +186,7 @@ def costGrad(self, g): gradnorm = np.sqrt(g.inner(g)) return gradnorm - def costHessian(self, zhat, Hzhat): + def hessian(self, zhat, Hzhat): """ Apply the the reduced Hessian to the vector :math:`zhat` @@ -268,7 +268,7 @@ def cost(self, z, order=0, **kwargs): return cost_risk+cost_penalization - def costGrad(self, g): + def grad(self, g): """ Computes the gradient of the cost functional @@ -282,7 +282,7 @@ def costGrad(self, g): g.zero() # Risk measure gradient - self.risk_measure.costGrad(g) + self.risk_measure.grad(g) if self.penalization is not None: self.penalization.grad(self.z, self.z_help) @@ -291,7 +291,7 @@ def costGrad(self, g): gradnorm = np.sqrt(g.inner(g)) return gradnorm - def costHessian(self, zhat, Hzhat): + def hessian(self, zhat, Hzhat): """ Apply the the reduced Hessian to the vector :math:`zhat` @@ -303,7 +303,7 @@ def costHessian(self, zhat, Hzhat): .. note:: Assumes :code:`self.cost` has been called with :code:`order >= 2` """ Hzhat.zero() - self.risk_measure.costHessian(zhat, Hzhat) + self.risk_measure.hessian(zhat, Hzhat) if self.penalization is not None: self.penalization.hessian(self.z, zhat, self.z_help) diff --git a/soupy/modeling/controlCostHessian.py b/soupy/modeling/controlCostHessian.py index b303f3c..362cffc 100644 --- a/soupy/modeling/controlCostHessian.py +++ b/soupy/modeling/controlCostHessian.py @@ -63,7 +63,7 @@ def mult(self, zhat, Hzhat): :type Hzhat: :py:class:`dolfin.Vector` or similar """ - self.cost_functional.costHessian(zhat, Hzhat) + self.cost_functional.hessian(zhat, Hzhat) self._ncalls += 1 diff --git a/soupy/modeling/meanVarRiskMeasure.py b/soupy/modeling/meanVarRiskMeasure.py index 4f9a205..d2afb89 100644 --- a/soupy/modeling/meanVarRiskMeasure.py +++ b/soupy/modeling/meanVarRiskMeasure.py @@ -139,7 +139,7 @@ def cost(self): """ return self.q_bar + self.beta * np.std(self.q_samples)**2 - def costGrad(self, g): + def grad(self, g): """ Evaluates the gradient of the risk measure once components have been computed @@ -153,5 +153,5 @@ def costGrad(self, g): g.axpy(2*self.beta, self.qg_bar) g.axpy(-2*self.beta*self.q_bar, self.g_bar) - def costHessian(self, zhat, Hzhat): + def hessian(self, zhat, Hzhat): logging.warning("Hessian not implemented for MeanVarRiskMeasure") diff --git a/soupy/modeling/meanVarRiskMeasureSAA.py b/soupy/modeling/meanVarRiskMeasureSAA.py index 39db992..4ef0b57 100644 --- a/soupy/modeling/meanVarRiskMeasureSAA.py +++ b/soupy/modeling/meanVarRiskMeasureSAA.py @@ -235,7 +235,7 @@ def cost(self): """ return self.q_bar + self.beta * (self.q2_bar - self.q_bar**2) - def costGrad(self, g): + def grad(self, g): """ Evaluates the gradient of the risk measure once components have been computed @@ -250,7 +250,7 @@ def costGrad(self, g): g.axpy(-2*self.beta*self.q_bar, self.g_bar) - def costHessian(self, zhat, Hzhat): + def hessian(self, zhat, Hzhat): """ Apply the hessian of the risk measure once components have been computed \ in the direction :code:`zhat` @@ -462,7 +462,7 @@ def cost(self): """ return self.q_bar + self.beta * np.std(self.q_samples)**2 - def costGrad(self, g): + def grad(self, g): """ Compute the gradient of the risk measure once components have been computed @@ -476,7 +476,7 @@ def costGrad(self, g): g.axpy(2*self.beta, self.qg_bar) g.axpy(-2*self.beta*self.q_bar, self.g_bar) - def costHessian(self, zhat, Hzhat): + def hessian(self, zhat, Hzhat): """ Apply the hessian of the risk measure once components have been computed \ in the direction :code:`zhat` diff --git a/soupy/modeling/riskMeasure.py b/soupy/modeling/riskMeasure.py index 777ddb5..458b9a1 100644 --- a/soupy/modeling/riskMeasure.py +++ b/soupy/modeling/riskMeasure.py @@ -44,7 +44,7 @@ def cost(self): """ raise NotImplementedError("Child class should implement method costValue") - def costGrad(self): + def grad(self): """ Evaluates the gradient of the risk measure once components have been computed @@ -53,9 +53,9 @@ def costGrad(self): .. note:: Assumes :code:`self.computeComponents` has been called with :code:`order >= 1` """ - raise NotImplementedError("Child class should implement method costGrad") + raise NotImplementedError("Child class should implement method grad") - def costHessian(self, zhat, Hzhat): + def hessian(self, zhat, Hzhat): """ Apply the hessian of the risk measure once components have been computed \ in the direction :code:`zhat` diff --git a/soupy/modeling/superquantileRiskMeasureSAA.py b/soupy/modeling/superquantileRiskMeasureSAA.py index 1a190e0..5b48397 100644 --- a/soupy/modeling/superquantileRiskMeasureSAA.py +++ b/soupy/modeling/superquantileRiskMeasureSAA.py @@ -263,7 +263,7 @@ def cost(self): t = self.zt.get_scalar() return t + 1/(1-self.beta) * self.s_bar - def costGrad(self, gt): + def grad(self, gt): """ Evaluates the gradient of the risk measure once components have been computed @@ -279,7 +279,7 @@ def costGrad(self, gt): gt.set_local(dz_np) - def costHessian(self, zt_hat, Hzt_hat): + def hessian(self, zt_hat, Hzt_hat): """ Apply the hessian of the risk measure once components have been computed \ in the direction :code:`zhat` diff --git a/soupy/optimization/bfgs.py b/soupy/optimization/bfgs.py index 602490b..f8224dc 100644 --- a/soupy/optimization/bfgs.py +++ b/soupy/optimization/bfgs.py @@ -293,7 +293,7 @@ def solve(self, z, H0inv=RescaledIdentity(), box_bounds=None, constraint_project self.BFGSop.H0inv.setPoint(z) g_old = g.copy() - gradnorm = self.cost_functional.costGrad(g) + gradnorm = self.cost_functional.grad(g) if gradnorm is None: gradnorm = np.sqrt(g.inner(g)) gradnorms.append(gradnorm) diff --git a/soupy/optimization/inexactNewtonCG.py b/soupy/optimization/inexactNewtonCG.py index 329cf80..f40541d 100644 --- a/soupy/optimization/inexactNewtonCG.py +++ b/soupy/optimization/inexactNewtonCG.py @@ -83,8 +83,8 @@ class InexactNewtonCG: More specifically the model object should implement following methods: - :code:`cost(z)` -> evaluate the cost functional - - :code:`costGrad(g)` -> evaluate the gradient and store to :code:`g` - - :code:`costHessian(zhat, Hzhat)` -> apply the cost Hessian + - :code:`grad(g)` -> evaluate the gradient and store to :code:`g` + - :code:`hessian(zhat, Hzhat)` -> apply the cost Hessian Type :code:`help(Model)` for additional information """ @@ -174,7 +174,7 @@ def _solve_ls(self, z): while (self.it < max_iter) and (self.converged == False): # Compute the gradient - gradnorm = self.cost_functional.costGrad(mg) + gradnorm = self.cost_functional.grad(mg) if self.it == 0: gradnorm_ini = gradnorm @@ -292,7 +292,7 @@ def _solve_tr(self,z): cost_old = self.cost_functional.cost(z, order=2) while (self.it < max_iter) and (self.converged == False): - gradnorm = self.cost_functional.costGrad(mg) + gradnorm = self.cost_functional.grad(mg) if self.it == 0: gradnorm_ini = gradnorm diff --git a/soupy/optimization/sgd.py b/soupy/optimization/sgd.py index df4d8fb..9583e81 100644 --- a/soupy/optimization/sgd.py +++ b/soupy/optimization/sgd.py @@ -48,7 +48,7 @@ class SGD: More specifically the cost functional object should implement following methods: - :code:`generate_vector()` -> generate the object containing state, parameter, adjoint, and control. - :code:`cost(z, order, rng)` -> evaluate the cost functional which allows a given :code:`rng` - - :code:`costGrad(g)` -> evaluate the gradient of the cost functional + - :code:`grad(g)` -> evaluate the gradient of the cost functional """ termination_reasons = [ "Maximum number of Iteration reached", #0 @@ -136,7 +136,7 @@ def solve(self, z, box_bounds=None, constraint_projection=None): # Run optimization while (self.it < max_iter) and (self.converged == False): - gradnorm = self.cost_functional.costGrad(g) + gradnorm = self.cost_functional.grad(g) if gradnorm is None: gradnorm = np.sqrt(g.inner(g)) diff --git a/soupy/optimization/steepestDescent.py b/soupy/optimization/steepestDescent.py index ce842a4..58e852f 100644 --- a/soupy/optimization/steepestDescent.py +++ b/soupy/optimization/steepestDescent.py @@ -46,7 +46,7 @@ class SteepestDescent: More specifically the cost functional object should implement following methods: - :code:`generate_vector()` -> generate the object containing state, parameter, adjoint, and control. - :code:`cost(z)` -> evaluate the cost functional - - :code:`costGrad(g)` -> evaluate the gradient of the cost functional + - :code:`grad(g)` -> evaluate the gradient of the cost functional """ termination_reasons = [ "Maximum number of Iteration reached", #0 @@ -139,7 +139,7 @@ def solve(self, z, box_bounds=None, constraint_projection=None): # Run optimization while (self.it < max_iter) and (self.converged == False): - gradnorm = self.cost_functional.costGrad(g) + gradnorm = self.cost_functional.grad(g) if gradnorm is None: gradnorm = np.sqrt(g.inner(g)) diff --git a/soupy/test/ptest_meanVarRiskMeasureSAA_MPI.py b/soupy/test/ptest_meanVarRiskMeasureSAA_MPI.py index 5ff8453..28a1233 100644 --- a/soupy/test/ptest_meanVarRiskMeasureSAA_MPI.py +++ b/soupy/test/ptest_meanVarRiskMeasureSAA_MPI.py @@ -195,14 +195,14 @@ def finiteDifferenceCheck(self, sample_size, is_fwd_linear=True): risk.computeComponents(z0, order=1) c0 = risk.cost() - risk.costGrad(g0) - risk.costHessian(dz, Hdz) + risk.grad(g0) + risk.hessian(dz, Hdz) rng = hp.Random() risk.computeComponents(z1, order=1) c1 = risk.cost() - risk.costGrad(g1) + risk.grad(g1) # Test gradients dcdz_fd = (c1 - c0)/self.delta diff --git a/soupy/test/ptest_superquantileSAA_MPI.py b/soupy/test/ptest_superquantileSAA_MPI.py index e1c743f..44660c9 100644 --- a/soupy/test/ptest_superquantileSAA_MPI.py +++ b/soupy/test/ptest_superquantileSAA_MPI.py @@ -205,12 +205,12 @@ def finiteDifferenceCheck(self, sample_size, smoothplus_type, epsilon, is_fwd_li risk.computeComponents(zt0, order=1) c0 = risk.cost() - risk.costGrad(gt0) - risk.costHessian(dzt, Hdzt) + risk.grad(gt0) + risk.hessian(dzt, Hdzt) risk.computeComponents(zt1, order=1) c1 = risk.cost() - risk.costGrad(gt1) + risk.grad(gt1) dcdz_fd = (c1 - c0)/self.delta dcdz_ad = gt0.inner(dzt) diff --git a/soupy/test/test_bfgs.py b/soupy/test/test_bfgs.py index 1a089e4..a081dbd 100644 --- a/soupy/test/test_bfgs.py +++ b/soupy/test/test_bfgs.py @@ -41,12 +41,12 @@ def cost(self, z, order=0, sample_size=1, rng=None): self.z = z return z.inner(z) - def costGrad(self, g): + def grad(self, g): g.zero() g.axpy(2.0, self.z) return g.inner(g) - def costHessian(self, zhat, Hzhat): + def hessian(self, zhat, Hzhat): Hzhat.zero() Hzhat.axpy(1.0, zhat) diff --git a/soupy/test/test_controlCostFunctional.py b/soupy/test/test_controlCostFunctional.py index 724eae4..d490633 100644 --- a/soupy/test/test_controlCostFunctional.py +++ b/soupy/test/test_controlCostFunctional.py @@ -116,15 +116,15 @@ def finiteDifferenceCheck(self, qoi_varf, is_fwd_linear=True): c0 = cost.cost(z0, order=2) g0 = model.generate_vector(CONTROL) - cost.costGrad(g0) - cost.costHessian(dz, Hdz) + cost.grad(g0) + cost.hessian(dz, Hdz) z1.axpy(1.0, z0) z1.axpy(self.delta, dz) c1 = cost.cost(z1, order=1) g1 = model.generate_vector(CONTROL) - cost.costGrad(g1) + cost.grad(g1) dcdz_fd = (c1 - c0)/self.delta dcdz_ad = g0.inner(dz) diff --git a/soupy/test/test_inexactNewtonCG.py b/soupy/test/test_inexactNewtonCG.py index 5595dc7..0a483b1 100644 --- a/soupy/test/test_inexactNewtonCG.py +++ b/soupy/test/test_inexactNewtonCG.py @@ -41,12 +41,12 @@ def cost(self, z, order=0, sample_size=1, rng=None): self.z = z return z.inner(z) - def costGrad(self, g): + def grad(self, g): g.zero() g.axpy(2.0, self.z) return g.inner(g) - def costHessian(self, zhat, Hzhat): + def hessian(self, zhat, Hzhat): Hzhat.zero() Hzhat.axpy(1.0, zhat) diff --git a/soupy/test/test_meanVarRiskMeasure.py b/soupy/test/test_meanVarRiskMeasure.py index 65bf7a4..7434bfe 100644 --- a/soupy/test/test_meanVarRiskMeasure.py +++ b/soupy/test/test_meanVarRiskMeasure.py @@ -114,7 +114,7 @@ def l2norm(u,m,z): risk.computeComponents(z0, order=1, sample_size=sample_size, rng=rng) c0 = risk.cost() - risk.costGrad(g0) + risk.grad(g0) rng = hp.Random() risk.computeComponents(z1, order=0, sample_size=sample_size, rng=rng) diff --git a/soupy/test/test_meanVarRiskMeasureSAA.py b/soupy/test/test_meanVarRiskMeasureSAA.py index ca043e5..370c1b6 100644 --- a/soupy/test/test_meanVarRiskMeasureSAA.py +++ b/soupy/test/test_meanVarRiskMeasureSAA.py @@ -190,14 +190,14 @@ def finiteDifferenceCheck(self, sample_size, is_fwd_linear=True): risk.computeComponents(z0, order=1) c0 = risk.cost() - risk.costGrad(g0) - risk.costHessian(dz, Hdz) + risk.grad(g0) + risk.hessian(dz, Hdz) rng = hp.Random() risk.computeComponents(z1, order=1, sample_size=sample_size) c1 = risk.cost() - risk.costGrad(g1) + risk.grad(g1) # Test gradients dcdz_fd = (c1 - c0)/self.delta diff --git a/soupy/test/test_meanVarRiskMeasureSAACostFunction.py b/soupy/test/test_meanVarRiskMeasureSAACostFunction.py index 2c70aa6..e266b35 100644 --- a/soupy/test/test_meanVarRiskMeasureSAACostFunction.py +++ b/soupy/test/test_meanVarRiskMeasureSAACostFunction.py @@ -150,12 +150,12 @@ def finiteDifferenceCheck(self, sample_size, is_fwd_linear=True): z1.axpy(self.delta, dz) c0 = costFun.cost(z0, order=2) - costFun.costGrad(g0) - costFun.costHessian(dz, Hdz) + costFun.grad(g0) + costFun.hessian(dz, Hdz) c1 = costFun.cost(z1, order=1) - costFun.costGrad(g1) - costFun.costGrad(g1) + costFun.grad(g1) + costFun.grad(g1) dcdz_fd = (c1 - c0)/self.delta dcdz_ad = g0.inner(dz) diff --git a/soupy/test/test_riskMeasureCostFunctional.py b/soupy/test/test_riskMeasureCostFunctional.py index 9a6d059..ef61c87 100644 --- a/soupy/test/test_riskMeasureCostFunctional.py +++ b/soupy/test/test_riskMeasureCostFunctional.py @@ -148,7 +148,7 @@ def l2norm(u,m,z): rng = hp.Random() c0 = costFun.cost(z0, order=1, sample_size=sample_size, rng=rng) - costFun.costGrad(g0) + costFun.grad(g0) rng = hp.Random() c1 = costFun.cost(z1, order=0, sample_size=sample_size, rng=rng) diff --git a/soupy/test/test_sgd.py b/soupy/test/test_sgd.py index ff6327a..e037615 100644 --- a/soupy/test/test_sgd.py +++ b/soupy/test/test_sgd.py @@ -40,12 +40,12 @@ def cost(self, z, order=0, sample_size=1, rng=None): self.z = z return z.inner(z) - def costGrad(self, g): + def grad(self, g): g.zero() g.axpy(2.0, self.z) return g.inner(g) - def costHessian(self, zhat, Hzhat): + def hessian(self, zhat, Hzhat): Hzhat.zero() Hzhat.axpy(1.0, zhat) diff --git a/soupy/test/test_steepestDescent.py b/soupy/test/test_steepestDescent.py index 1261e26..8ec079a 100644 --- a/soupy/test/test_steepestDescent.py +++ b/soupy/test/test_steepestDescent.py @@ -41,12 +41,12 @@ def cost(self, z, order=0, sample_size=1, rng=None): self.z = z return z.inner(z) - def costGrad(self, g): + def grad(self, g): g.zero() g.axpy(2.0, self.z) return g.inner(g) - def costHessian(self, zhat, Hzhat): + def hessian(self, zhat, Hzhat): Hzhat.zero() Hzhat.axpy(1.0, zhat) diff --git a/soupy/test/test_superquantileSAA.py b/soupy/test/test_superquantileSAA.py index dc6f2e7..a62282a 100644 --- a/soupy/test/test_superquantileSAA.py +++ b/soupy/test/test_superquantileSAA.py @@ -177,13 +177,13 @@ def finiteDifferenceCheck(self, sample_size, smoothplus_type, epsilon, is_fwd_li risk.computeComponents(zt0, order=1) c0 = risk.cost() - risk.costGrad(gt0) - risk.costHessian(dzt, Hdzt) + risk.grad(gt0) + risk.hessian(dzt, Hdzt) rng = hp.Random() risk.computeComponents(zt1, order=1) c1 = risk.cost() - risk.costGrad(gt1) + risk.grad(gt1) dcdz_fd = (c1 - c0)/self.delta dcdz_ad = gt0.inner(dzt) diff --git a/soupy/utils/scipyCostWrapper.py b/soupy/utils/scipyCostWrapper.py index 1071827..82b3604 100644 --- a/soupy/utils/scipyCostWrapper.py +++ b/soupy/utils/scipyCostWrapper.py @@ -16,18 +16,18 @@ class ScipyCostWrapper: """ - Class to interface the controlCostFunctional with a \ + Class to interface the :py:class:`soupy.ControlCostFunctional` with a \ scipy optimizer. Converts inputs to functions taking \ and returning :code:`numpy` arrays """ - def __init__(self, controlCostFunctional, verbose=False): + def __init__(self, cost_functional, verbose=False): """ Constructor - :param controlCostFunctional: The cost functional to wrap - :type controlCostFunctional: :py:class:`ControlCostFunctional` + :param cost_functional: The cost functional to wrap + :type cost_functional: :py:class:`ControlCostFunctional` """ - self.cost_functional = controlCostFunctional + self.cost_functional = cost_functional self.z_help = self.cost_functional.generate_vector(CONTROL) self.z_hat_help = self.cost_functional.generate_vector(CONTROL) self.z_out_help = self.cost_functional.generate_vector(CONTROL) @@ -70,7 +70,7 @@ def grad(self, z_np): self.z_help.set_local(z_np) self.z_help.apply("") self.cost_functional.cost(self.z_help, order=1) - self.cost_functional.costGrad(self.z_out_help) + self.cost_functional.grad(self.z_out_help) if self.verbose and MPI.COMM_WORLD.Get_rank() == 0: print("Gradient evaluation") self.n_grad += 1 @@ -85,7 +85,7 @@ def hessian(self, z_np, zhat_np): self.z_hat_help.apply("") self.cost_functional.cost(self.z_help, order=2) - self.cost_functional.costHessian(self.z_hat_help, self.z_out_help) + self.cost_functional.hessian(self.z_hat_help, self.z_out_help) if self.verbose and MPI.COMM_WORLD.Get_rank() == 0: print("Hessian action evaluation") self.n_hess += 1 diff --git a/soupy/utils/stochasticCostFiniteDifference.py b/soupy/utils/stochasticCostFiniteDifference.py index 322be6a..a503ef9 100644 --- a/soupy/utils/stochasticCostFiniteDifference.py +++ b/soupy/utils/stochasticCostFiniteDifference.py @@ -38,7 +38,7 @@ def stochasticCostFiniteDifference(pde_cost, z, dz, delta=1e-3, sample_size=1): rng = hp.Random() c0 = pde_cost.cost(z0, rng=rng, order=1, sample_size=sample_size) - pde_cost.costGrad(gz) + pde_cost.grad(gz) rng = hp.Random() delta = 1e-3 @@ -70,7 +70,7 @@ def SAACostFiniteDifference(pde_cost, z, dz, delta=1e-3): rng = hp.Random() c0 = pde_cost.cost(z0, order=1) - pde_cost.costGrad(gz) + pde_cost.grad(gz) rng = hp.Random() delta = 1e-3 From 0023f044acd0ecd392761b9f52c0bbec492babe3 Mon Sep 17 00:00:00 2001 From: Dingcheng Date: Thu, 3 Aug 2023 15:10:29 -0500 Subject: [PATCH 2/9] deprecated meanVarRiskMeasureSAA and renamed other risk measures --- soupy/modeling/__init__.py | 6 +- soupy/modeling/meanVarRiskMeasureSAA.py | 215 +----------------- ...ure.py => meanVarRiskMeasureStochastic.py} | 8 +- soupy/modeling/superquantileRiskMeasureSAA.py | 2 +- ..._MPI.py => ptest_meanVarRiskMeasureSAA.py} | 8 +- soupy/test/ptest_scipyCostWrapper.py | 4 +- ...leSAA_MPI.py => ptest_superquantileSAA.py} | 90 ++++---- ...tFunction.py => test_SAACostFunctional.py} | 0 ...y => test_meanVarRiskMeasureStochastic.py} | 8 +- ...al.py => test_stochasticCostFunctional.py} | 10 +- soupy/test/test_superquantileSAA.py | 4 +- 11 files changed, 69 insertions(+), 286 deletions(-) rename soupy/modeling/{meanVarRiskMeasure.py => meanVarRiskMeasureStochastic.py} (97%) rename soupy/test/{ptest_meanVarRiskMeasureSAA_MPI.py => ptest_meanVarRiskMeasureSAA.py} (96%) rename soupy/test/{ptest_superquantileSAA_MPI.py => ptest_superquantileSAA.py} (84%) rename soupy/test/{test_meanVarRiskMeasureSAACostFunction.py => test_SAACostFunctional.py} (100%) rename soupy/test/{test_meanVarRiskMeasure.py => test_meanVarRiskMeasureStochastic.py} (94%) rename soupy/test/{test_riskMeasureCostFunctional.py => test_stochasticCostFunctional.py} (94%) diff --git a/soupy/modeling/__init__.py b/soupy/modeling/__init__.py index e073a20..5b7153f 100644 --- a/soupy/modeling/__init__.py +++ b/soupy/modeling/__init__.py @@ -25,9 +25,9 @@ from .controlQoI import ControlQoI, L2MisfitVarfHandler, VariationalControlQoI, L2MisfitControlQoI -from .meanVarRiskMeasure import meanVarRiskMeasureSettings, MeanVarRiskMeasure +from .meanVarRiskMeasureStochastic import meanVarRiskMeasureStochasticSettings, MeanVarRiskMeasureStochastic -from .meanVarRiskMeasureSAA import meanVarRiskMeasureSAASettings, MeanVarRiskMeasureSAA, MeanVarRiskMeasureSAA_MPI +from .meanVarRiskMeasureSAA import meanVarRiskMeasureSAASettings, MeanVarRiskMeasureSAA from .penalization import Penalization, L2Penalization, WeightedL2Penalization, MultiPenalization @@ -35,7 +35,7 @@ from .smoothPlusApproximation import SmoothPlusApproximationQuartic, SmoothPlusApproximationSoftplus -from .superquantileRiskMeasureSAA import SuperquantileRiskMeasureSAA_MPI, superquantileRiskMeasureSAASettings, sampleSuperquantile, sampleSuperquantileByMinimization +from .superquantileRiskMeasureSAA import SuperquantileRiskMeasureSAA, superquantileRiskMeasureSAASettings, sampleSuperquantile, sampleSuperquantileByMinimization from .variables import STATE, PARAMETER, ADJOINT, CONTROL diff --git a/soupy/modeling/meanVarRiskMeasureSAA.py b/soupy/modeling/meanVarRiskMeasureSAA.py index 4ef0b57..5d713a9 100644 --- a/soupy/modeling/meanVarRiskMeasureSAA.py +++ b/soupy/modeling/meanVarRiskMeasureSAA.py @@ -34,9 +34,7 @@ def meanVarRiskMeasureSAASettings(data = {}): return ParameterList(data) - - -class MeanVarRiskMeasureSAA_MPI(RiskMeasure): +class MeanVarRiskMeasureSAA(RiskMeasure): """ Mean + variance risk measure using sample average approximation @@ -305,214 +303,3 @@ def gatherSamples(self): return q_all - - - - -class MeanVarRiskMeasureSAA(RiskMeasure): - """ - The mean + variance risk measure using sample average approximation - - .. math:: \\rho[Q](z) = \mathbb{E}_m[Q(m,z)] + \\beta \mathbb{V}_m[Q(m,z)] - """ - - def __init__(self, control_model, prior, settings = meanVarRiskMeasureSAASettings()): - """ - Constructor: - - :param control_model: control model containing the :code:`soupy.PDEVariationalControlProblem` - and :code:`soupy.ControlQoI` - :type control_model: :py:class:`soupy.ControlModel` - :param prior: The prior distribution for the random parameter - :type prior: :py:class:`hippylib.Prior` - :param settings: additional settings - :type settings: :py:class:`hippylib.ParameterList` - """ - self.model = control_model - self.prior = prior - self.settings = settings - self.sample_size = self.settings['sample_size'] - self.beta = settings['beta'] - - # Aggregate components for computing cost, grad, hess - self.g = self.model.generate_vector(CONTROL) - self.q_samples = np.zeros(self.sample_size) - - self.q_bar = 0 - self.g_bar = self.model.generate_vector(CONTROL) - self.qg_bar = self.model.generate_vector(CONTROL) - - # For sampling - self.noise = dl.Vector() - self.prior.init_vector(self.noise, "noise") - rng = Random(seed=self.settings['seed']) - - # Generate samples for m - self.x_mc = [] - self.g_mc = [] - self.z = self.model.generate_vector(CONTROL) - logging.info("Initial sampling of stochastic parameter") - for i in range(self.sample_size): - ui = self.model.generate_vector(STATE) - mi = self.model.generate_vector(PARAMETER) - pi = self.model.generate_vector(ADJOINT) - x = [ui, mi, pi, self.z] - rng.normal(1.0, self.noise) - self.prior.sample(self.noise, mi) - self.x_mc.append(x) - - gi = self.model.generate_vector(CONTROL) - self.g_mc.append(gi) - - # Helper variable - self.Hz_helper = self.model.generate_vector(CONTROL) - self.diff_helper= self.model.generate_vector(CONTROL) - - - self.control_model_hessian = None - self.Hi_zhat = self.model.generate_vector(CONTROL) - - - self.has_adjoint_solve = False - self.has_forward_solve = False - - def generate_vector(self, component = "ALL"): - return self.model.generate_vector(component) - - def computeComponents(self, z, order=0, **kwargs): - """ - - Computes the components for the evaluation of the risk measure - - :param z: the control variable - :type z: :py:class:`dolfin.Vector` - :param order: Order of the derivatives needed. - :code:`0` for cost, :code:`1` for gradient, :code:`2` for Hessian - :type order: int - :param kwargs: dummy keyword arguments for compatibility - """ - - # Check if a new control variable is used - new_forward_solve = False - self.diff_helper.zero() - self.diff_helper.axpy(1.0, self.z) - self.diff_helper.axpy(-1.0, z) - diff_norm = np.sqrt(self.diff_helper.inner(self.diff_helper)) - - # Check if new forward solve is needed - if diff_norm > dl.DOLFIN_EPS or not self.has_forward_solve: - # Update control variable (changes all samples) - # Ask that new forward and adjoint solves are computed - self.z.zero() - self.z.axpy(1.0, z) - new_forward_solve = True - logging.info("Using new forward solve") - - # Check if a new adjoint solve is needed - if new_forward_solve: - # Make sure previous adjoint solve is no longer valid - self.has_adjoint_solve = False - new_adjoint_solve = True - self.g_bar.zero() - self.qg_bar.zero() - - elif self.has_adjoint_solve: - # If it's not a new forward solve, check if we already have an adjoint - new_adjoint_solve = False - else: - # If we don't already have an adjoint, compute a new one - new_adjoint_solve = True - # In this case, self.has_adjoint_solve is already False - self.g_bar.zero() - self.qg_bar.zero() - - - if new_forward_solve or new_adjoint_solve: - # Now actually compute the solves - for i in range(self.sample_size): - x = self.x_mc[i] - # Solve state - if new_forward_solve: - self.model.solveFwd(x[STATE], x) - self.q_samples[i] = self.model.cost(x) - - if order >= 1 and new_adjoint_solve: - self.model.solveAdj(x[ADJOINT], x) - self.model.evalGradientControl(x, self.g_mc[i]) - self.g_bar.axpy(1/self.sample_size, self.g_mc[i]) - self.qg_bar.axpy(self.q_samples[i]/self.sample_size, self.g_mc[i]) - - self.q_bar = np.mean(self.q_samples) - - if order >= 1: - # We have computed a new adjoint solve - # Don't need to alter if it's any other case - self.has_adjoint_solve = True - - # Will always be true after solving - self.has_forward_solve = True - - def cost(self): - """ - Computes the value of the risk measure once components have been computed - - :return: Value of the cost functional - - .. note:: Assumes :code:`computeComponents` has been called with :code:`order>=0` - """ - return self.q_bar + self.beta * np.std(self.q_samples)**2 - - def grad(self, g): - """ - Compute the gradient of the risk measure once components have been computed - - :param g: (Dual of) the gradient of the risk measure to store the result in - :type g: :py:class:`dolfin.Vector` - - .. note:: Assumes :code:`self.computeComponents` has been called with :code:`order >= 1` - """ - g.zero() - g.axpy(1.0, self.g_bar) - g.axpy(2*self.beta, self.qg_bar) - g.axpy(-2*self.beta*self.q_bar, self.g_bar) - - def hessian(self, zhat, Hzhat): - """ - Apply the hessian of the risk measure once components have been computed \ - in the direction :code:`zhat` - - :param zhat: Direction for application of Hessian action of the risk measure - :type zhat: :py:class:`dolfin.Vector` - :param Hzhat: (Dual of) Result of the Hessian action of the risk measure - to store the result in - :type Hzhat: :py:class:`dolfin.Vector` - - .. note:: Assumes :code:`self.computeComponents` has been called with :code:`order >= 1` - """ - if self.control_model_hessian is None: - self.control_model_hessian = ControlModelHessian(self.model) - - Hzhat.zero() - for i in range(self.sample_size): - xi = self.x_mc[i] - gi = self.g_mc[i] - qi = self.q_samples[i] - - self.model.setLinearizationPoint(xi) - self.control_model_hessian.mult(zhat, self.Hi_zhat) - - hessian_scale_factor = 1 + 2*self.beta*qi - 2*self.beta*self.q_bar - - # Apply :math:`H_i z` the qoi hessian at sample - Hzhat.axpy(hessian_scale_factor/self.sample_size, self.Hi_zhat) - - # Apply :math:`g_i g_i^T`for qoi gradients at sample - gradient_scale_factor = gi.inner(zhat) * 2*self.beta - Hzhat.axpy(gradient_scale_factor/self.sample_size, gi) - - # Apply :math:`\bar{g} \bar{g}^T` for mean qoi gradients - mean_gradient_scale_factor = self.g_bar.inner(zhat) * 2*self.beta - Hzhat.axpy(-mean_gradient_scale_factor, self.g_bar) - - - diff --git a/soupy/modeling/meanVarRiskMeasure.py b/soupy/modeling/meanVarRiskMeasureStochastic.py similarity index 97% rename from soupy/modeling/meanVarRiskMeasure.py rename to soupy/modeling/meanVarRiskMeasureStochastic.py index d2afb89..ff629cf 100644 --- a/soupy/modeling/meanVarRiskMeasure.py +++ b/soupy/modeling/meanVarRiskMeasureStochastic.py @@ -23,12 +23,12 @@ from .variables import STATE, PARAMETER, ADJOINT, CONTROL -def meanVarRiskMeasureSettings(data = {}): +def meanVarRiskMeasureStochasticSettings(data = {}): data['beta'] = [0,'Weighting factor for variance'] return ParameterList(data) -class MeanVarRiskMeasure(RiskMeasure): +class MeanVarRiskMeasureStochastic(RiskMeasure): """ Class for memory efficient evaluation of the Mean + Variance risk measure @@ -37,7 +37,7 @@ class MeanVarRiskMeasure(RiskMeasure): that allows for stochastic approximations and SGD """ - def __init__(self, control_model, prior, settings = meanVarRiskMeasureSettings()): + def __init__(self, control_model, prior, settings = meanVarRiskMeasureStochasticSettings()): """ Constructor: @@ -154,4 +154,4 @@ def grad(self, g): g.axpy(-2*self.beta*self.q_bar, self.g_bar) def hessian(self, zhat, Hzhat): - logging.warning("Hessian not implemented for MeanVarRiskMeasure") + logging.warning("Hessian not implemented for MeanVarRiskMeasureStochastic") diff --git a/soupy/modeling/superquantileRiskMeasureSAA.py b/soupy/modeling/superquantileRiskMeasureSAA.py index 5b48397..46e5a5f 100644 --- a/soupy/modeling/superquantileRiskMeasureSAA.py +++ b/soupy/modeling/superquantileRiskMeasureSAA.py @@ -60,7 +60,7 @@ def superquantileRiskMeasureSAASettings(data = {}): return ParameterList(data) -class SuperquantileRiskMeasureSAA_MPI(RiskMeasure): +class SuperquantileRiskMeasureSAA(RiskMeasure): """ Risk measure for the sample average approximation of the superquantile risk measure (CVaR) with sample parallelism using MPI diff --git a/soupy/test/ptest_meanVarRiskMeasureSAA_MPI.py b/soupy/test/ptest_meanVarRiskMeasureSAA.py similarity index 96% rename from soupy/test/ptest_meanVarRiskMeasureSAA_MPI.py rename to soupy/test/ptest_meanVarRiskMeasureSAA.py index 28a1233..69daf7f 100644 --- a/soupy/test/ptest_meanVarRiskMeasureSAA_MPI.py +++ b/soupy/test/ptest_meanVarRiskMeasureSAA.py @@ -28,7 +28,7 @@ sys.path.append('../../') from soupy import VariationalControlQoI, ControlModel, \ meanVarRiskMeasureSAASettings, \ - MeanVarRiskMeasureSAA_MPI,\ + MeanVarRiskMeasureSAA,\ PDEVariationalControlProblem, \ STATE, PARAMETER, CONTROL @@ -43,7 +43,7 @@ def qoi_for_testing(u,m,z): return u**2*dl.dx + dl.exp(m) * dl.inner(dl.grad(u), dl.grad(u))*dl.ds -class TestMeanVarRiskMeasureSAA_MPI(unittest.TestCase): +class TestMeanVarRiskMeasureSAA(unittest.TestCase): def setUp(self): self.reltol = 1e-3 self.fdtol = 1e-2 @@ -78,7 +78,7 @@ def l2norm(u,m,z): qoi = VariationalControlQoI(self.mesh, self.Vh, l2norm) model = ControlModel(pde, qoi) - risk = MeanVarRiskMeasureSAA_MPI(model, prior, comm_sampler=self.comm_sampler) + risk = MeanVarRiskMeasureSAA(model, prior, comm_sampler=self.comm_sampler) z = model.generate_vector(CONTROL) c_val = risk.cost() print("Before computing: ", c_val) @@ -104,7 +104,7 @@ def _setup_mean_var_risk_measure(self, model, prior, sample_size, beta): rm_settings = meanVarRiskMeasureSAASettings() rm_settings['sample_size'] = sample_size rm_settings['beta'] = beta - risk = MeanVarRiskMeasureSAA_MPI(model, prior, rm_settings, comm_sampler=self.comm_sampler) + risk = MeanVarRiskMeasureSAA(model, prior, rm_settings, comm_sampler=self.comm_sampler) return risk def testSavedSolution(self): diff --git a/soupy/test/ptest_scipyCostWrapper.py b/soupy/test/ptest_scipyCostWrapper.py index 3bdcc96..fbe12d6 100644 --- a/soupy/test/ptest_scipyCostWrapper.py +++ b/soupy/test/ptest_scipyCostWrapper.py @@ -30,7 +30,7 @@ sys.path.append('../../') from soupy import ControlCostFunctional, PDEVariationalControlProblem, \ VariationalControlQoI, L2Penalization, \ - ControlModel, MeanVarRiskMeasureSAA_MPI, meanVarRiskMeasureSAASettings, \ + ControlModel, MeanVarRiskMeasureSAA, meanVarRiskMeasureSAASettings, \ RiskMeasureControlCostFunctional, ScipyCostWrapper, \ STATE, PARAMETER, ADJOINT, CONTROL, \ MultipleSerialPDEsCollective @@ -101,7 +101,7 @@ def _setup_mean_var_risk_measure(self, model, prior, sample_size, beta): rm_settings = meanVarRiskMeasureSAASettings() rm_settings['sample_size'] = sample_size rm_settings['beta'] = beta - risk = MeanVarRiskMeasureSAA_MPI(model, prior, rm_settings, comm_sampler=self.comm_sampler) + risk = MeanVarRiskMeasureSAA(model, prior, rm_settings, comm_sampler=self.comm_sampler) return risk diff --git a/soupy/test/ptest_superquantileSAA_MPI.py b/soupy/test/ptest_superquantileSAA.py similarity index 84% rename from soupy/test/ptest_superquantileSAA_MPI.py rename to soupy/test/ptest_superquantileSAA.py index 44660c9..7d3d409 100644 --- a/soupy/test/ptest_superquantileSAA_MPI.py +++ b/soupy/test/ptest_superquantileSAA.py @@ -31,7 +31,7 @@ sys.path.append('../../') from soupy import VariationalControlQoI, ControlModel, \ - superquantileRiskMeasureSAASettings, SuperquantileRiskMeasureSAA_MPI, \ + superquantileRiskMeasureSAASettings, SuperquantileRiskMeasureSAA, \ STATE, PARAMETER, ADJOINT, CONTROL, MultipleSerialPDEsCollective from setupPoissonControlProblem import poisson_control_settings, setupPoissonPDEProblem @@ -65,9 +65,9 @@ def qoi_for_testing(u,m,z): return u**2*dl.dx + dl.exp(m) * dl.inner(dl.grad(u), dl.grad(u))*dl.ds -class TestSuperquantileSAA_MPI(unittest.TestCase): +class TestSuperquantileSAA(unittest.TestCase): """ - Testing :code:`soupy.SuperquantileRiskMeasureSAA_MPI` with MPI \ + Testing :code:`soupy.SuperquantileRiskMeasureSAA` with MPI \ capabilities enabled """ def setUp(self): @@ -115,7 +115,7 @@ def _setup_superquantile_risk_measure(self, model, prior, sample_size, beta, rm_settings['beta'] = beta rm_settings['smoothplus_type'] rm_settings['epsilon'] - risk = SuperquantileRiskMeasureSAA_MPI(model, prior, rm_settings, comm_sampler=self.comm_sampler) + risk = SuperquantileRiskMeasureSAA(model, prior, rm_settings, comm_sampler=self.comm_sampler) return risk def testSavedSolution(self): @@ -265,7 +265,7 @@ def testFiniteDifference(self): def computeSuperquantileValue(self, sample_size, beta): """ - Evaluate superquantile using :code:`soupy.SuperquantileRiskMeasureSAA_MPI.superquantile` \ + Evaluate superquantile using :code:`soupy.SuperquantileRiskMeasureSAA.superquantile` \ for samples drawn from a normal distribution """ @@ -285,7 +285,7 @@ def l2norm(u,m,z): rm_settings = superquantileRiskMeasureSAASettings() rm_settings['sample_size'] = sample_size rm_settings['beta'] = beta - risk = SuperquantileRiskMeasureSAA_MPI(model, prior, rm_settings, comm_sampler=self.comm_sampler) + risk = SuperquantileRiskMeasureSAA(model, prior, rm_settings, comm_sampler=self.comm_sampler) np.random.seed(1) # Sample all of the samples @@ -343,47 +343,43 @@ def testSuperquantileValue(self): self.assertTrue(np.abs(sq_normal - sq) < tol) -# def testCostValue(self): -# beta = 0.95 -# sample_size = 100 -# smoothplus_types = ["softplus", "quartic"] -# epsilons = [1e-3, 1e-4] -# settings = poisson_control_settings() -# settings['N_WELLS_PER_SIDE'] = self.n_wells_per_side -# -# tol = 5e-2 -# def l2norm(u,m,z): -# return u**2*dl.dx -# -# for epsilon, smoothplus_type in zip(epsilons, smoothplus_types): -# rm_settings = superquantileRiskMeasureSAASettings() -# rm_settings['sample_size'] = sample_size -# rm_settings['beta'] = beta -# rm_settings['smoothplus_type'] = smoothplus_type -# rm_settings['epsilon'] = epsilon -# -# pde, prior, control_dist = setupPoissonPDEProblem(self.Vh, settings) -# qoi = VariationalControlQoI(self.mesh, self.Vh, l2norm) -# model = ControlModel(pde, qoi) -# risk = SuperquantileRiskMeasureSAA_MPI(model, prior, rm_settings) -# -# zt = risk.generate_vector(CONTROL) -# z = zt.get_vector() -# control_dist.sample(z) -# # z.set_local(np.random.randn(len(z.get_local()))) -# -# risk.computeComponents(zt, order=0) -# quantile = np.percentile(risk.q_samples, beta * 100) -# zt.set_scalar(quantile) -# -# risk.computeComponents(zt, order=0) -# sq_by_cost = risk.cost() -# sq_by_samples = risk.superquantile() -# -# print("Superquantile from cost: ", sq_by_cost) -# print("Superquantile from samples: ", sq_by_samples) -# err = np.abs(sq_by_cost - sq_by_samples) -# self.assertTrue(err/np.abs(sq_by_samples) < tol) + def testCostValue(self): + BETA = 0.95 + SAMPLE_SIZE = 100 + IS_FWD_LINEAR = False + + pde, prior, control_dist = self._setup_pde_and_distributions(IS_FWD_LINEAR) + qoi, model = self._setup_control_model(pde, l2_norm) + + + smoothplus_types = ["softplus", "quartic"] + epsilons = [1e-3, 1e-4] + settings = poisson_control_settings() + settings['N_WELLS_PER_SIDE'] = self.n_wells_per_side + + tol = 5e-2 + def l2norm(u,m,z): + return u**2*dl.dx + + for epsilon, smoothplus_type in zip(epsilons, smoothplus_types): + risk = self._setup_superquantile_risk_measure(model, prior, + SAMPLE_SIZE, BETA, smoothplus_type=smoothplus_type, epsilon=epsilon) + + zt = risk.generate_vector(CONTROL) + sample_zt_and_bcast(control_dist, zt, self.collective_sampler) + + risk.computeComponents(zt, order=0) + quantile = np.percentile(risk.gatherSamples(), BETA * 100) + zt.set_scalar(quantile) + + risk.computeComponents(zt, order=0) + sq_by_cost = risk.cost() + sq_by_samples = risk.superquantile() + + print("Superquantile from cost: ", sq_by_cost) + print("Superquantile from samples: ", sq_by_samples) + err = np.abs(sq_by_cost - sq_by_samples) + self.assertTrue(err/np.abs(sq_by_samples) < tol) if __name__ == "__main__": diff --git a/soupy/test/test_meanVarRiskMeasureSAACostFunction.py b/soupy/test/test_SAACostFunctional.py similarity index 100% rename from soupy/test/test_meanVarRiskMeasureSAACostFunction.py rename to soupy/test/test_SAACostFunctional.py diff --git a/soupy/test/test_meanVarRiskMeasure.py b/soupy/test/test_meanVarRiskMeasureStochastic.py similarity index 94% rename from soupy/test/test_meanVarRiskMeasure.py rename to soupy/test/test_meanVarRiskMeasureStochastic.py index 7434bfe..3e74c21 100644 --- a/soupy/test/test_meanVarRiskMeasure.py +++ b/soupy/test/test_meanVarRiskMeasureStochastic.py @@ -27,7 +27,7 @@ sys.path.append('../../') from soupy import VariationalControlQoI, ControlModel, \ - meanVarRiskMeasureSettings, MeanVarRiskMeasure,\ + meanVarRiskMeasureStochasticSettings, MeanVarRiskMeasureStochastic,\ STATE, PARAMETER, ADJOINT, CONTROL from setupPoissonControlProblem import poisson_control_settings, setupPoissonPDEProblem @@ -66,7 +66,7 @@ def l2norm(u,m,z): qoi = VariationalControlQoI(self.mesh, self.Vh, l2norm) model = ControlModel(pde, qoi) - risk = MeanVarRiskMeasure(model, prior) + risk = MeanVarRiskMeasureStochastic(model, prior) z = model.generate_vector(CONTROL) c_val = risk.cost() print("Before computing: ", c_val) @@ -94,10 +94,10 @@ def l2norm(u,m,z): qoi = VariationalControlQoI(self.mesh, self.Vh, l2norm) model = ControlModel(pde, qoi) - rm_settings = meanVarRiskMeasureSettings() + rm_settings = meanVarRiskMeasureStochasticSettings() rm_settings['beta'] = 0.5 - risk = MeanVarRiskMeasure(model, prior, rm_settings) + risk = MeanVarRiskMeasureStochastic(model, prior, rm_settings) z0 = model.generate_vector(CONTROL) dz = model.generate_vector(CONTROL) diff --git a/soupy/test/test_riskMeasureCostFunctional.py b/soupy/test/test_stochasticCostFunctional.py similarity index 94% rename from soupy/test/test_riskMeasureCostFunctional.py rename to soupy/test/test_stochasticCostFunctional.py index ef61c87..7451d73 100644 --- a/soupy/test/test_riskMeasureCostFunctional.py +++ b/soupy/test/test_stochasticCostFunctional.py @@ -29,7 +29,7 @@ sys.path.append('../../') from soupy import ControlCostFunctional, PDEVariationalControlProblem, \ VariationalControlQoI, L2Penalization, \ - ControlModel, MeanVarRiskMeasure, meanVarRiskMeasureSettings, \ + ControlModel, MeanVarRiskMeasureStochastic, meanVarRiskMeasureStochasticSettings, \ RiskMeasureControlCostFunctional, \ STATE, PARAMETER, ADJOINT, CONTROL @@ -67,7 +67,7 @@ def l2norm(u,m,z): qoi = VariationalControlQoI(self.mesh, self.Vh, l2norm) model = ControlModel(pde, qoi) - risk = MeanVarRiskMeasure(model, prior) + risk = MeanVarRiskMeasureStochastic(model, prior) alpha = 2.0 if use_penalization: penalty = L2Penalization(self.Vh, alpha) @@ -84,7 +84,7 @@ def l2norm(u,m,z): rng = hp.Random() c_val = costFun.cost(z, order=0, sample_size=sample_size, rng=rng) - if isinstance(risk, MeanVarRiskMeasure): + if isinstance(risk, MeanVarRiskMeasureStochastic): print("Stochastic approximation risk measure") print("Sample size = %d" %(len(risk.q_samples))) self.assertEqual(len(risk.q_samples), sample_size) @@ -127,10 +127,10 @@ def l2norm(u,m,z): qoi = VariationalControlQoI(self.mesh, self.Vh, l2norm) model = ControlModel(pde, qoi) - rm_settings = meanVarRiskMeasureSettings() + rm_settings = meanVarRiskMeasureStochasticSettings() rm_settings['beta'] = 0.5 - risk = MeanVarRiskMeasure(model, prior, rm_settings) + risk = MeanVarRiskMeasureStochastic(model, prior, rm_settings) alpha = 2.0 penalty = L2Penalization(self.Vh, alpha) costFun = RiskMeasureControlCostFunctional(risk, penalty) diff --git a/soupy/test/test_superquantileSAA.py b/soupy/test/test_superquantileSAA.py index a62282a..99342be 100644 --- a/soupy/test/test_superquantileSAA.py +++ b/soupy/test/test_superquantileSAA.py @@ -29,7 +29,7 @@ sys.path.append('../../') from soupy import VariationalControlQoI, ControlModel, \ - superquantileRiskMeasureSAASettings, SuperquantileRiskMeasureSAA_MPI, \ + superquantileRiskMeasureSAASettings, SuperquantileRiskMeasureSAA, \ STATE, PARAMETER, ADJOINT, CONTROL from setupPoissonControlProblem import poisson_control_settings, setupPoissonPDEProblem @@ -84,7 +84,7 @@ def _setup_superquantile_risk_measure(self, model, prior, sample_size, beta, rm_settings['beta'] = beta rm_settings['smoothplus_type'] rm_settings['epsilon'] - risk = SuperquantileRiskMeasureSAA_MPI(model, prior, rm_settings) + risk = SuperquantileRiskMeasureSAA(model, prior, rm_settings) return risk def testSavedSolution(self): From acc2b53919d8e1cf21dc82b34ac9671a46164de2 Mon Sep 17 00:00:00 2001 From: Dingcheng Date: Thu, 3 Aug 2023 15:16:18 -0500 Subject: [PATCH 3/9] updated naming --- .github/workflows/run_tests.sh | 4 ++-- doc/soupy.modeling.rst | 6 +++--- examples/poisson/compare_results.py | 4 ++-- examples/poisson/driver_poisson_mean_mpi.py | 2 +- examples/semilinear_elliptic/driver_semilinear_cvar.py | 2 +- soupy/modeling/__init__.py | 2 +- soupy/modeling/meanVarRiskMeasureSAA.py | 2 +- soupy/modeling/superquantileRiskMeasureSAA.py | 10 +++++----- soupy/test/ptest_meanVarRiskMeasureSAA.py | 2 +- soupy/test/ptest_superquantileSAA.py | 4 ++-- 10 files changed, 19 insertions(+), 19 deletions(-) diff --git a/.github/workflows/run_tests.sh b/.github/workflows/run_tests.sh index 2513a29..e217304 100644 --- a/.github/workflows/run_tests.sh +++ b/.github/workflows/run_tests.sh @@ -10,6 +10,6 @@ cd soupy/test $PYTHON -m unittest discover -v $PYTHON -m unittest discover -v -p 'ptest_*' -mpirun -n 4 $PYTHON ptest_meanVarRiskMeasureSAA_MPI.py +mpirun -n 4 $PYTHON ptest_meanVarRiskMeasureSAA.py mpirun -n 4 $PYTHON ptest_scipyCostWrapper.py -mpirun -n 4 $PYTHON ptest_superquantileSAA_MPI.py +mpirun -n 4 $PYTHON ptest_superquantileSAA.py diff --git a/doc/soupy.modeling.rst b/doc/soupy.modeling.rst index 60f9032..8bbb3e1 100644 --- a/doc/soupy.modeling.rst +++ b/doc/soupy.modeling.rst @@ -62,10 +62,10 @@ soupy.modeling.riskMeasure :show-inheritance: -soupy.modeling.meanVarRiskMeasure ---------------------------------------- +soupy.modeling.meanVarRiskMeasureStochastic +-------------------------------------------- -.. automodule:: soupy.modeling.meanVarRiskMeasure +.. automodule:: soupy.modeling.meanVarRiskMeasureStochastic :members: :undoc-members: :show-inheritance: diff --git a/examples/poisson/compare_results.py b/examples/poisson/compare_results.py index 51a0efd..e82eef1 100644 --- a/examples/poisson/compare_results.py +++ b/examples/poisson/compare_results.py @@ -34,7 +34,7 @@ def main(): N_ELEMENTS_Y = 20 TOL = 1e-12 - RESULTS_DIRECTORY_MPI = "results_mpi" + RESULTS_DIRECTORY = "results_mpi" RESULTS_DIRECTORY_SERIAL = "results_serial" mesh = dl.UnitSquareMesh(N_ELEMENTS_X, N_ELEMENTS_Y) @@ -46,7 +46,7 @@ def main(): load_file.read(z_opt_no_mpi, "control") z_opt_with_mpi = dl.Function(Vh_STATE) - with dl.HDF5File(mesh.mpi_comm(), "%s/z_opt.h5" %(RESULTS_DIRECTORY_MPI), "r") as load_file: + with dl.HDF5File(mesh.mpi_comm(), "%s/z_opt.h5" %(RESULTS_DIRECTORY), "r") as load_file: load_file.read(z_opt_with_mpi, "control") z_diff = z_opt_no_mpi.vector().copy() diff --git a/examples/poisson/driver_poisson_mean_mpi.py b/examples/poisson/driver_poisson_mean_mpi.py index f6e7981..2c2e4cd 100644 --- a/examples/poisson/driver_poisson_mean_mpi.py +++ b/examples/poisson/driver_poisson_mean_mpi.py @@ -93,7 +93,7 @@ def boundary(x, on_boundary): risk_settings = soupy.meanVarRiskMeasureSAASettings() risk_settings["beta"] = VARIANCE_WEIGHT risk_settings["sample_size"] = SAMPLE_SIZE -risk_measure = soupy.MeanVarRiskMeasureSAA_MPI(control_model, prior, risk_settings, comm_sampler=comm_sampler) +risk_measure = soupy.MeanVarRiskMeasureSAA(control_model, prior, risk_settings, comm_sampler=comm_sampler) # 7. Define the penalization term penalty = soupy.L2Penalization(Vh, PENALTY_WEIGHT) diff --git a/examples/semilinear_elliptic/driver_semilinear_cvar.py b/examples/semilinear_elliptic/driver_semilinear_cvar.py index cf12817..a699b02 100644 --- a/examples/semilinear_elliptic/driver_semilinear_cvar.py +++ b/examples/semilinear_elliptic/driver_semilinear_cvar.py @@ -139,7 +139,7 @@ def print_on_root(print_str, mpi_comm=MPI.COMM_WORLD): rm_param["sample_size"] = args.N_sample rm_param["seed"] = seed rm_param["epsilon"] = 1e-4 - pde_rm = soupy.SuperquantileRiskMeasureSAA_MPI(control_model, prior, settings=rm_param, comm_sampler=comm_sampler) + pde_rm = soupy.SuperquantileRiskMeasureSAA(control_model, prior, settings=rm_param, comm_sampler=comm_sampler) pde_cost = soupy.RiskMeasureControlCostFunctional(pde_rm, None) # Scipy cost diff --git a/soupy/modeling/__init__.py b/soupy/modeling/__init__.py index 5b7153f..8060a14 100644 --- a/soupy/modeling/__init__.py +++ b/soupy/modeling/__init__.py @@ -35,7 +35,7 @@ from .smoothPlusApproximation import SmoothPlusApproximationQuartic, SmoothPlusApproximationSoftplus -from .superquantileRiskMeasureSAA import SuperquantileRiskMeasureSAA, superquantileRiskMeasureSAASettings, sampleSuperquantile, sampleSuperquantileByMinimization +from .superquantileRiskMeasureSAA import SuperquantileRiskMeasureSAA, superquantileRiskMeasureSAASettings, sample_superquantile, sample_superquantile_by_minimization from .variables import STATE, PARAMETER, ADJOINT, CONTROL diff --git a/soupy/modeling/meanVarRiskMeasureSAA.py b/soupy/modeling/meanVarRiskMeasureSAA.py index 5d713a9..5148272 100644 --- a/soupy/modeling/meanVarRiskMeasureSAA.py +++ b/soupy/modeling/meanVarRiskMeasureSAA.py @@ -291,7 +291,7 @@ def hessian(self, zhat, Hzhat): Hzhat.axpy(-mean_gradient_scale_factor, self.g_bar) - def gatherSamples(self): + def gather_samples(self): """ Gather the QoI samples from all processes diff --git a/soupy/modeling/superquantileRiskMeasureSAA.py b/soupy/modeling/superquantileRiskMeasureSAA.py index 46e5a5f..712d677 100644 --- a/soupy/modeling/superquantileRiskMeasureSAA.py +++ b/soupy/modeling/superquantileRiskMeasureSAA.py @@ -30,7 +30,7 @@ -def sampleSuperquantile(samples, beta): +def sample_superquantile(samples, beta): """ Evaluate superquantile from samples """ @@ -38,7 +38,7 @@ def sampleSuperquantile(samples, beta): return quantile + np.mean(np.maximum(samples - quantile, 0))/(1-beta) -def sampleSuperquantileByMinimization(samples, beta, epsilon=1e-2): +def sample_superquantile_by_minimization(samples, beta, epsilon=1e-2): """ Evaluate superquantile from samples by minimization """ @@ -336,7 +336,7 @@ def hessian(self, zt_hat, Hzt_hat): Hzt_hat.set_scalar(Ht_hat_all) - def gatherSamples(self): + def gather_samples(self): """ Gather the QoI samples from all processes @@ -356,8 +356,8 @@ def superquantile(self): .. note:: Assumes :code:`computeComponents` has been called with :code:`order>=0` """ - q_all = self.gatherSamples() - value = sampleSuperquantile(q_all, self.beta) + q_all = self.gather_samples() + value = sample_superquantile(q_all, self.beta) return value diff --git a/soupy/test/ptest_meanVarRiskMeasureSAA.py b/soupy/test/ptest_meanVarRiskMeasureSAA.py index 69daf7f..ce5fc3f 100644 --- a/soupy/test/ptest_meanVarRiskMeasureSAA.py +++ b/soupy/test/ptest_meanVarRiskMeasureSAA.py @@ -248,7 +248,7 @@ def checkGatherSamples(self, sample_size): rank_sampler = self.comm_sampler.Get_rank() risk.q_samples = np.zeros(risk.q_samples.shape) + rank_sampler - q_all = risk.gatherSamples() + q_all = risk.gather_samples() number_equal_to_rank = np.sum(np.isclose(q_all, rank_sampler)) print("Testing samples are gathering correctly") diff --git a/soupy/test/ptest_superquantileSAA.py b/soupy/test/ptest_superquantileSAA.py index 7d3d409..fce940d 100644 --- a/soupy/test/ptest_superquantileSAA.py +++ b/soupy/test/ptest_superquantileSAA.py @@ -313,7 +313,7 @@ def checkGatherSamples(self, sample_size): rank_sampler = self.comm_sampler.Get_rank() risk.q_samples = np.zeros(risk.q_samples.shape) + rank_sampler - q_all = risk.gatherSamples() + q_all = risk.gather_samples() number_equal_to_rank = np.sum(np.isclose(q_all, rank_sampler)) print("Testing samples are gathering correctly") @@ -369,7 +369,7 @@ def l2norm(u,m,z): sample_zt_and_bcast(control_dist, zt, self.collective_sampler) risk.computeComponents(zt, order=0) - quantile = np.percentile(risk.gatherSamples(), BETA * 100) + quantile = np.percentile(risk.gather_samples(), BETA * 100) zt.set_scalar(quantile) risk.computeComponents(zt, order=0) From cd2fb709ca6d5b775ba5a96a9ec6f23236dae0f2 Mon Sep 17 00:00:00 2001 From: Dingcheng Date: Thu, 3 Aug 2023 15:23:32 -0500 Subject: [PATCH 4/9] removed duplicate tests for SAA risk measures --- soupy/test/test_meanVarRiskMeasureSAA.py | 239 ------------------ soupy/test/test_superquantileSAA.py | 296 ----------------------- 2 files changed, 535 deletions(-) delete mode 100644 soupy/test/test_meanVarRiskMeasureSAA.py delete mode 100644 soupy/test/test_superquantileSAA.py diff --git a/soupy/test/test_meanVarRiskMeasureSAA.py b/soupy/test/test_meanVarRiskMeasureSAA.py deleted file mode 100644 index 370c1b6..0000000 --- a/soupy/test/test_meanVarRiskMeasureSAA.py +++ /dev/null @@ -1,239 +0,0 @@ -# Copyright (c) 2023, The University of Texas at Austin -# & Georgia Institute of Technology -# -# All Rights reserved. -# See file COPYRIGHT for details. -# -# This file is part of the SOUPy package. For more information see -# https://github.com/hippylib/soupy/ -# -# SOUPy is free software; you can redistribute it and/or modify it under the -# terms of the GNU General Public License (as published by the Free -# Software Foundation) version 3.0 dated June 2007. - -import unittest -import dolfin as dl -import numpy as np -import matplotlib.pyplot as plt - -import logging -logging.getLogger('FFC').setLevel(logging.WARNING) -logging.getLogger('UFL').setLevel(logging.WARNING) -dl.set_log_active(False) - -import os, sys -sys.path.append(os.environ.get('HIPPYLIB_PATH')) -import hippylib as hp - -sys.path.append('../../') -from soupy import VariationalControlQoI, ControlModel, \ - meanVarRiskMeasureSAASettings, MeanVarRiskMeasureSAA,\ - PDEVariationalControlProblem, \ - STATE, PARAMETER, CONTROL - -from setupPoissonControlProblem import poisson_control_settings, setupPoissonPDEProblem - -def l2_norm(u,m,z): - return u**2*dl.dx - -def qoi_for_testing(u,m,z): - return u**2*dl.dx + dl.exp(m) * dl.inner(dl.grad(u), dl.grad(u))*dl.ds - - -class TestMeanVarRiskMeasureSAA(unittest.TestCase): - def setUp(self): - self.reltol = 1e-3 - self.fdtol = 1e-1 - self.delta = 1e-3 - self.n_wells_per_side = 3 - self.nx = 20 - self.ny = 20 - self.n_control = self.n_wells_per_side**2 - - # Make spaces - self.mesh = dl.UnitSquareMesh(self.nx, self.ny) - Vh_STATE = dl.FunctionSpace(self.mesh, "CG", 1) - Vh_PARAMETER = dl.FunctionSpace(self.mesh, "CG", 1) - Vh_CONTROL = dl.VectorFunctionSpace(self.mesh, "R", degree=0, dim=self.n_control) - self.Vh = [Vh_STATE, Vh_PARAMETER, Vh_STATE, Vh_CONTROL] - - def _setup_pde_and_distributions(self, is_fwd_linear): - settings = poisson_control_settings() - settings['nx'] = self.nx - settings['ny'] = self.ny - settings['N_WELLS_PER_SIDE'] = self.n_wells_per_side - settings['LINEAR'] = is_fwd_linear - pde, prior, control_dist = setupPoissonPDEProblem(self.Vh, settings) - return pde, prior, control_dist - - def _setup_control_model(self, pde, qoi_varf): - qoi = VariationalControlQoI(self.mesh, self.Vh, qoi_varf) - model = ControlModel(pde, qoi) - return qoi, model - - def _setup_mean_var_risk_measure(self, model, prior, sample_size, beta): - rm_settings = meanVarRiskMeasureSAASettings() - rm_settings['sample_size'] = sample_size - rm_settings['beta'] = beta - risk = MeanVarRiskMeasureSAA(model, prior, rm_settings) - return risk - - def testCostValue(self): - settings = poisson_control_settings() - settings['nx'] = self.nx - settings['ny'] = self.ny - settings['N_WELLS_PER_SIDE'] = self.n_wells_per_side - settings['LINEAR'] = True - - pde, prior, _ = setupPoissonPDEProblem(self.Vh, settings) - def l2norm(u,m,z): - return u**2*dl.dx + (m - dl.Constant(1.0))**2*dl.dx - - qoi = VariationalControlQoI(self.mesh, self.Vh, l2norm) - model = ControlModel(pde, qoi) - - risk = MeanVarRiskMeasureSAA(model, prior) - z = model.generate_vector(CONTROL) - c_val = risk.cost() - print("Before computing: ", c_val) - risk.computeComponents(z) - c_val = risk.cost() - print("After computing: ", c_val) - - def testSavedSolution(self): - IS_FWD_LINEAR = False - SAMPLE_SIZE = 5 - BETA = 0.5 - - pde, prior, control_dist = self._setup_pde_and_distributions(IS_FWD_LINEAR) - qoi, model = self._setup_control_model(pde, l2_norm) - risk = self._setup_mean_var_risk_measure(model, prior, SAMPLE_SIZE, BETA) - - z0 = model.generate_vector(CONTROL) - z1 = model.generate_vector(CONTROL) - - # np.random.seed(1) - control_dist.sample(z0) - control_dist.sample(z1) - - print("Test if correct solution and adjoint solves are being stored") - - # initial cost - risk.computeComponents(z0, order=0) - self.assertFalse(risk.has_adjoint_solve) - c0 = risk.cost() - - # same place, check same - risk.computeComponents(z0, order=0) - self.assertFalse(risk.has_adjoint_solve) - self.assertAlmostEqual(c0, risk.cost()) - - # now with gradient - risk.computeComponents(z0, order=1) - self.assertAlmostEqual(c0, risk.cost()) - self.assertTrue(risk.has_adjoint_solve) - - # new cost, no gradient - risk.computeComponents(z1, order=0) - c1 = risk.cost() - self.assertFalse(risk.has_adjoint_solve) - - # back to old cost - risk.computeComponents(z0, order=0) - self.assertAlmostEqual(c0, risk.cost()) - self.assertFalse(risk.has_adjoint_solve) - - # new cost, with gradient - risk.computeComponents(z1, order=1) - self.assertTrue(risk.has_adjoint_solve) - self.assertAlmostEqual(c1, risk.cost()) - - # new cost, no gradient - risk.computeComponents(z1, order=0) - self.assertTrue(risk.has_adjoint_solve) - self.assertAlmostEqual(c1, risk.cost()) - - # old cost, no gradient - risk.computeComponents(z0, order=1) - self.assertAlmostEqual(c0, risk.cost()) - self.assertTrue(risk.has_adjoint_solve) - - - - - def finiteDifferenceCheck(self, sample_size, is_fwd_linear=True): - print("-" * 80) - print("Finite difference check with %d samples" %(sample_size)) - if is_fwd_linear: - print("Linear problem") - else: - print("Nonlinear problem") - - BETA = 0.5 - pde, prior, control_dist = self._setup_pde_and_distributions(is_fwd_linear) - qoi, model = self._setup_control_model(pde, qoi_for_testing) - risk = self._setup_mean_var_risk_measure(model, prior, sample_size, BETA) - - z0 = model.generate_vector(CONTROL) - dz = model.generate_vector(CONTROL) - z1 = model.generate_vector(CONTROL) - g0 = model.generate_vector(CONTROL) - g1 = model.generate_vector(CONTROL) - Hdz = model.generate_vector(CONTROL) - - # np.random.seed(1) - control_dist.sample(z0) - control_dist.sample(dz) - z1.axpy(1.0, z0) - z1.axpy(self.delta, dz) - - risk.computeComponents(z0, order=1) - - c0 = risk.cost() - risk.grad(g0) - risk.hessian(dz, Hdz) - - rng = hp.Random() - risk.computeComponents(z1, order=1, sample_size=sample_size) - - c1 = risk.cost() - risk.grad(g1) - - # Test gradients - dcdz_fd = (c1 - c0)/self.delta - dcdz_ad = g0.inner(dz) - print("Initial cost: ", c0) - print("New cost: ", c1) - print("Finite difference derivative: %g" %(dcdz_fd)) - print("Adjoint derivative: %g" %(dcdz_ad)) - self.assertTrue(abs((dcdz_fd - dcdz_ad)/dcdz_ad) < self.fdtol) - - # Test hessians - Hdz_fd = (g1.get_local() - g0.get_local())/self.delta - Hdz_ad = Hdz.get_local() - print("Finite difference Hessian action \n", Hdz_fd) - print("Adjoint Hessian action \n", Hdz_ad) - err_hess = np.linalg.norm(Hdz_fd - Hdz_ad) - print("Norm error: %g" %(err_hess)) - self.assertTrue(err_hess/np.linalg.norm(Hdz_ad) < self.fdtol) - - def testFiniteDifferenceLinearProblem(self): - is_fwd_linear = True - n_samples = 1 - self.finiteDifferenceCheck(n_samples, is_fwd_linear) - - n_samples = 10 - self.finiteDifferenceCheck(n_samples, is_fwd_linear) - - def testFiniteDifferenceNonlinearProblem(self): - is_fwd_linear = False - n_samples = 1 - self.finiteDifferenceCheck(n_samples, is_fwd_linear) - - n_samples = 10 - self.finiteDifferenceCheck(n_samples, is_fwd_linear) - - -if __name__ == "__main__": - unittest.main() - diff --git a/soupy/test/test_superquantileSAA.py b/soupy/test/test_superquantileSAA.py deleted file mode 100644 index 99342be..0000000 --- a/soupy/test/test_superquantileSAA.py +++ /dev/null @@ -1,296 +0,0 @@ -# Copyright (c) 2023, The University of Texas at Austin -# & Georgia Institute of Technology -# -# All Rights reserved. -# See file COPYRIGHT for details. -# -# This file is part of the SOUPy package. For more information see -# https://github.com/hippylib/soupy/ -# -# SOUPy is free software; you can redistribute it and/or modify it under the -# terms of the GNU General Public License (as published by the Free -# Software Foundation) version 3.0 dated June 2007. - -import unittest -import dolfin as dl -import numpy as np -import matplotlib.pyplot as plt -import scipy.stats -import math - -import logging -logging.getLogger('FFC').setLevel(logging.WARNING) -logging.getLogger('UFL').setLevel(logging.WARNING) -dl.set_log_active(False) - -import os, sys -sys.path.append(os.environ.get('HIPPYLIB_PATH')) -import hippylib as hp - -sys.path.append('../../') -from soupy import VariationalControlQoI, ControlModel, \ - superquantileRiskMeasureSAASettings, SuperquantileRiskMeasureSAA, \ - STATE, PARAMETER, ADJOINT, CONTROL - -from setupPoissonControlProblem import poisson_control_settings, setupPoissonPDEProblem - - -def standardNormalSuperquantile(beta): - quantile = scipy.stats.norm.ppf(beta) - return np.exp(-quantile**2/2)/(1-beta)/np.sqrt(2*math.pi) - -def l2_norm(u,m,z): - return u**2*dl.dx - -def qoi_for_testing(u,m,z): - return u**2*dl.dx + dl.exp(m) * dl.inner(dl.grad(u), dl.grad(u))*dl.ds - -class TestSuperquantileSAA(unittest.TestCase): - def setUp(self): - self.reltol = 1e-3 - self.fdtol = 1e-2 - self.delta = 1e-3 - self.n_wells_per_side = 3 - self.nx = 20 - self.ny = 20 - self.n_control = self.n_wells_per_side**2 - - # Make spaces - self.mesh = dl.UnitSquareMesh(self.nx, self.ny) - Vh_STATE = dl.FunctionSpace(self.mesh, "CG", 1) - Vh_PARAMETER = dl.FunctionSpace(self.mesh, "CG", 1) - Vh_CONTROL = dl.VectorFunctionSpace(self.mesh, "R", degree=0, dim=self.n_control) - self.Vh = [Vh_STATE, Vh_PARAMETER, Vh_STATE, Vh_CONTROL] - - def _setup_pde_and_distributions(self, is_fwd_linear): - settings = poisson_control_settings() - settings['nx'] = self.nx - settings['ny'] = self.ny - settings['N_WELLS_PER_SIDE'] = self.n_wells_per_side - settings['LINEAR'] = is_fwd_linear - pde, prior, control_dist = setupPoissonPDEProblem(self.Vh, settings) - return pde, prior, control_dist - - def _setup_control_model(self, pde, qoi_varf): - qoi = VariationalControlQoI(self.mesh, self.Vh, qoi_varf) - model = ControlModel(pde, qoi) - return qoi, model - - - def _setup_superquantile_risk_measure(self, model, prior, sample_size, beta, - smoothplus_type='quartic', epsilon=1e-4): - rm_settings = superquantileRiskMeasureSAASettings() - rm_settings['sample_size'] = sample_size - rm_settings['beta'] = beta - rm_settings['smoothplus_type'] - rm_settings['epsilon'] - risk = SuperquantileRiskMeasureSAA(model, prior, rm_settings) - return risk - - def testSavedSolution(self): - """ - Testing correct solutions are saved internally - """ - IS_FWD_LINEAR = False - SAMPLE_SIZE = 100 - BETA = 0.5 - - pde, prior, control_dist = self._setup_pde_and_distributions(IS_FWD_LINEAR) - qoi, model = self._setup_control_model(pde, l2_norm) - risk = self._setup_superquantile_risk_measure(model, prior, SAMPLE_SIZE, BETA) - - zt0 = risk.generate_vector(CONTROL) - zt1 = risk.generate_vector(CONTROL) - - # np.random.seed(1) - control_dist.sample(zt0.get_vector()) - control_dist.sample(zt1.get_vector()) - zt0.set_scalar(np.random.randn()) - zt1.set_scalar(np.random.randn()) - - print("Test if correct solution and adjoint solves are being stored") - - # initial cost - risk.computeComponents(zt0, order=0) - self.assertFalse(risk.has_adjoint_solve) - c0 = risk.cost() - - # same place, check same - risk.computeComponents(zt0, order=0) - self.assertFalse(risk.has_adjoint_solve) - self.assertAlmostEqual(c0, risk.cost()) - - # now with gradient - risk.computeComponents(zt0, order=1) - self.assertAlmostEqual(c0, risk.cost()) - self.assertTrue(risk.has_adjoint_solve) - - # new cost, no gradient - risk.computeComponents(zt1, order=0) - c1 = risk.cost() - self.assertFalse(risk.has_adjoint_solve) - - # back to old cost - risk.computeComponents(zt0, order=0) - self.assertAlmostEqual(c0, risk.cost()) - self.assertFalse(risk.has_adjoint_solve) - - # new cost, with gradient - risk.computeComponents(zt1, order=1) - self.assertTrue(risk.has_adjoint_solve) - self.assertAlmostEqual(c1, risk.cost()) - - # new cost, no gradient - risk.computeComponents(zt1, order=0) - self.assertTrue(risk.has_adjoint_solve) - self.assertAlmostEqual(c1, risk.cost()) - - # old cost, no gradient - risk.computeComponents(zt0, order=1) - self.assertAlmostEqual(c0, risk.cost()) - self.assertTrue(risk.has_adjoint_solve) - - - def finiteDifferenceCheck(self, sample_size, smoothplus_type, epsilon, is_fwd_linear=True): - BETA = 0.5 - - pde, prior, control_dist = self._setup_pde_and_distributions(is_fwd_linear) - qoi, model = self._setup_control_model(pde, qoi_for_testing) - risk = self._setup_superquantile_risk_measure(model, prior, sample_size, BETA, - smoothplus_type=smoothplus_type, epsilon=epsilon) - - zt0 = risk.generate_vector(CONTROL) - dzt = risk.generate_vector(CONTROL) - zt1 = risk.generate_vector(CONTROL) - gt0 = risk.generate_vector(CONTROL) - gt1 = risk.generate_vector(CONTROL) - Hdzt = risk.generate_vector(CONTROL) - - control_dist.sample(zt0.get_vector()) - zt0.set_scalar(np.random.randn()) - - control_dist.sample(dzt.get_vector()) - dzt.set_scalar(np.random.randn()) - - zt1.axpy(1.0, zt0) - zt1.axpy(self.delta, dzt) - - risk.computeComponents(zt0, order=1) - c0 = risk.cost() - risk.grad(gt0) - risk.hessian(dzt, Hdzt) - - rng = hp.Random() - risk.computeComponents(zt1, order=1) - c1 = risk.cost() - risk.grad(gt1) - - dcdz_fd = (c1 - c0)/self.delta - dcdz_ad = gt0.inner(dzt) - - print(40*"-") - if is_fwd_linear: - print("Linear, %d Samples, Smooth plus type: %s" %(sample_size, smoothplus_type)) - else: - print("Nonlinear, %d Samples, Smooth plus type: %s" %(sample_size, smoothplus_type)) - - print("Q: ", risk.q_samples) - print("Initial cost: ", c0) - print("New cost: ", c1) - print("Finite difference derivative: %g" %(dcdz_fd)) - print("Adjoint derivative: %g" %(dcdz_ad)) - self.assertTrue(abs((dcdz_fd - dcdz_ad)/dcdz_ad) < self.fdtol) - - Hdz_fd = (gt1.get_local() - gt0.get_local())/self.delta - Hdz_ad = Hdzt.get_local() - print("Finite difference Hessian action \n", Hdz_fd) - print("Adjoint Hessian action \n", Hdz_ad) - err_hess = np.linalg.norm(Hdz_fd - Hdz_ad) - print("Norm error: %g" %(err_hess)) - - if np.allclose(Hdz_ad, 0): - self.assertTrue(err_hess < self.fdtol) - else: - norm_hess = np.linalg.norm(Hdz_ad) - self.assertTrue(err_hess/norm_hess < self.fdtol) - - print(40*"-") - - - def testFiniteDifference(self): - sample_sizes = [1, 10] - linearities = [True, False] - smoothpluses = ["softplus", "quartic"] - epsilons = [1e-1, 1e-4] - - for sample_size in sample_sizes: - for linearity in linearities: - for epsilon, smoothplus in zip(epsilons, smoothpluses): - self.finiteDifferenceCheck(sample_size, smoothplus, epsilon, linearity) - - def computeSuperquantileValue(self, sample_size, beta): - IS_FWD_LINEAR = False - pde, prior, control_dist = self._setup_pde_and_distributions(IS_FWD_LINEAR) - qoi, model = self._setup_control_model(pde, l2_norm) - risk = self._setup_superquantile_risk_measure(model, prior, sample_size, beta) - np.random.seed(1) - risk.q_samples = np.random.randn(len(risk.q_samples)) - sq = risk.superquantile() - return sq - - - def testSuperquantileValue(self): - sample_size = 10000 - beta = 0.2 - sq_normal = standardNormalSuperquantile(beta) - sq = self.computeSuperquantileValue(sample_size, beta) - print("Computed superquantile: ", sq) - print("Analytic superquantile: ", sq_normal) - tol = 1e-2 - self.assertTrue(np.abs(sq_normal - sq) < tol) - - - def testCostValue(self): - BETA = 0.95 - SAMPLE_SIZE = 100 - IS_FWD_LINEAR = False - - pde, prior, control_dist = self._setup_pde_and_distributions(IS_FWD_LINEAR) - qoi, model = self._setup_control_model(pde, l2_norm) - - - smoothplus_types = ["softplus", "quartic"] - epsilons = [1e-3, 1e-4] - settings = poisson_control_settings() - settings['N_WELLS_PER_SIDE'] = self.n_wells_per_side - - tol = 5e-2 - def l2norm(u,m,z): - return u**2*dl.dx - - for epsilon, smoothplus_type in zip(epsilons, smoothplus_types): - risk = self._setup_superquantile_risk_measure(model, prior, - SAMPLE_SIZE, BETA, smoothplus_type=smoothplus_type, epsilon=epsilon) - - zt = risk.generate_vector(CONTROL) - z = zt.get_vector() - control_dist.sample(z) - # z.set_local(np.random.randn(len(z.get_local()))) - - risk.computeComponents(zt, order=0) - quantile = np.percentile(risk.q_samples, BETA * 100) - zt.set_scalar(quantile) - - risk.computeComponents(zt, order=0) - sq_by_cost = risk.cost() - sq_by_samples = risk.superquantile() - - print("Superquantile from cost: ", sq_by_cost) - print("Superquantile from samples: ", sq_by_samples) - err = np.abs(sq_by_cost - sq_by_samples) - self.assertTrue(err/np.abs(sq_by_samples) < tol) - - -if __name__ == "__main__": - unittest.main() - From 5c81eca3887ad8dc81c1b88d58e12595f830f541 Mon Sep 17 00:00:00 2001 From: Dingcheng Date: Thu, 3 Aug 2023 15:30:09 -0500 Subject: [PATCH 5/9] removed serial example from poisson examples --- examples/poisson/compare_results.py | 61 ------- ...son_mean_mpi.py => driver_poisson_mean.py} | 4 +- .../poisson/driver_poisson_mean_serial.py | 157 ------------------ 3 files changed, 2 insertions(+), 220 deletions(-) delete mode 100644 examples/poisson/compare_results.py rename examples/poisson/{driver_poisson_mean_mpi.py => driver_poisson_mean.py} (98%) delete mode 100644 examples/poisson/driver_poisson_mean_serial.py diff --git a/examples/poisson/compare_results.py b/examples/poisson/compare_results.py deleted file mode 100644 index e82eef1..0000000 --- a/examples/poisson/compare_results.py +++ /dev/null @@ -1,61 +0,0 @@ -# Copyright (c) 2023, The University of Texas at Austin -# & Georgia Institute of Technology -# -# All Rights reserved. -# See file COPYRIGHT for details. -# -# This file is part of the SOUPy package. For more information see -# https://github.com/hippylib/soupy/ -# -# SOUPy is free software; you can redistribute it and/or modify it under the -# terms of the GNU General Public License (as published by the Free -# Software Foundation) version 3.0 dated June 1991. - - -import sys -import os - -import dolfin as dl -import numpy as np -import matplotlib.pyplot as plt - -sys.path.append(os.environ.get("HIPPYLIB_PATH")) -sys.path.append("../../") - -import hippylib as hp -import soupy - - -def main(): - """ - Compares the solutions with/without mpi - """ - N_ELEMENTS_X = 20 - N_ELEMENTS_Y = 20 - TOL = 1e-12 - - RESULTS_DIRECTORY = "results_mpi" - RESULTS_DIRECTORY_SERIAL = "results_serial" - - mesh = dl.UnitSquareMesh(N_ELEMENTS_X, N_ELEMENTS_Y) - Vh_STATE = dl.FunctionSpace(mesh, "CG", 1) - - - z_opt_no_mpi = dl.Function(Vh_STATE) - with dl.HDF5File(mesh.mpi_comm(), "%s/z_opt.h5" %(RESULTS_DIRECTORY_SERIAL), "r") as load_file: - load_file.read(z_opt_no_mpi, "control") - - z_opt_with_mpi = dl.Function(Vh_STATE) - with dl.HDF5File(mesh.mpi_comm(), "%s/z_opt.h5" %(RESULTS_DIRECTORY), "r") as load_file: - load_file.read(z_opt_with_mpi, "control") - - z_diff = z_opt_no_mpi.vector().copy() - z_diff.axpy(-1.0, z_opt_with_mpi.vector()) - - z_error = z_diff.inner(z_diff) - print("Error between MPI and no MPI runs: %.3e" %(z_error)) - assert z_error < TOL - - -if __name__ == "__main__": - main() diff --git a/examples/poisson/driver_poisson_mean_mpi.py b/examples/poisson/driver_poisson_mean.py similarity index 98% rename from examples/poisson/driver_poisson_mean_mpi.py rename to examples/poisson/driver_poisson_mean.py index 2c2e4cd..bd63dd0 100644 --- a/examples/poisson/driver_poisson_mean_mpi.py +++ b/examples/poisson/driver_poisson_mean.py @@ -45,7 +45,7 @@ SAMPLE_SIZE = 4 PENALTY_WEIGHT = 1e-3 -RESULTS_DIRECTORY = "results_mpi" +RESULTS_DIRECTORY = "results" # 0. MPI Communicators comm_mesh = MPI.COMM_SELF @@ -103,7 +103,7 @@ def boundary(x, on_boundary): # 9. Define the optimizer -optimizer = soupy.BFGS(cost_functional) +optimizer = soupy.InexactNewtonCG(cost_functional) # 10. Provide initial guess and solve if comm_sampler.rank == 0 : diff --git a/examples/poisson/driver_poisson_mean_serial.py b/examples/poisson/driver_poisson_mean_serial.py deleted file mode 100644 index 42d4764..0000000 --- a/examples/poisson/driver_poisson_mean_serial.py +++ /dev/null @@ -1,157 +0,0 @@ -# Copyright (c) 2023, The University of Texas at Austin -# & Georgia Institute of Technology -# -# All Rights reserved. -# See file COPYRIGHT for details. -# -# This file is part of the SOUPy package. For more information see -# https://github.com/hippylib/soupy/ -# -# SOUPy is free software; you can redistribute it and/or modify it under the -# terms of the GNU General Public License (as published by the Free -# Software Foundation) version 3.0 dated June 1991. - -""" -Minimization of the mean risk measure using the SAA approximation for the -linear poisson problem with a log-normal conductivity field. - -This driver is purely serial. For a parallel implementation, see -:code:`driver_poisson_mean_mpi.py` -""" - -import sys -import os - -import dolfin as dl -import numpy as np -import matplotlib.pyplot as plt - -sys.path.append(os.environ.get("HIPPYLIB_PATH")) -sys.path.append("../../") -import hippylib as hp -import soupy -dl.set_log_active(False) - -N_ELEMENTS_X = 20 -N_ELEMENTS_Y = 20 -IS_FWD_LINEAR = True - -PRIOR_GAMMA = 10.0 -PRIOR_DELTA = 50.0 -PRIOR_MEAN = -2.0 - -VARIANCE_WEIGHT = 0.0 -SAMPLE_SIZE = 4 -PENALTY_WEIGHT = 1e-3 - -RESULTS_DIRECTORY = "results_serial" - -# 1. Setup function spaces -mesh = dl.UnitSquareMesh(N_ELEMENTS_X, N_ELEMENTS_Y) -Vh_STATE = dl.FunctionSpace(mesh, "CG", 1) -Vh_PARAMETER = dl.FunctionSpace(mesh, "CG", 1) -Vh_CONTROL = dl.FunctionSpace(mesh, "CG", 1) -Vh = [Vh_STATE, Vh_PARAMETER, Vh_STATE, Vh_CONTROL] - -# 2. Define PDE Problem -# This is a linear poisson equation where the control variable -# is the distributed source term. Dirichlet Boundary conditions -# are on the left and bottom boundaries - -def residual(u,m,p,z): - return dl.exp(m)*dl.inner(dl.grad(u), dl.grad(p))*dl.dx - z * p *dl.dx - -def boundary(x, on_boundary): - return on_boundary and (dl.near(x[0], 0) or dl.near(x[1], 0)) - -boundary_value = dl.Expression("x[1]", degree=1) - -bc = dl.DirichletBC(Vh_STATE, boundary_value, boundary) -bc0 = dl.DirichletBC(Vh_STATE, dl.Constant(0.0), boundary) -pde = soupy.PDEVariationalControlProblem(Vh, residual, [bc], [bc0], is_fwd_linear=IS_FWD_LINEAR) - -# 3. Define a Gaussian prior for the random parameter -mean_vector = dl.interpolate(dl.Constant(PRIOR_MEAN), Vh_PARAMETER).vector() -prior = hp.BiLaplacianPrior(Vh_PARAMETER, PRIOR_GAMMA, PRIOR_DELTA, mean=mean_vector, robin_bc=True) - -# 4. Define the QoI -u_target = dl.Expression("x[1] + sin(k*x[0]) * sin(k*x[1])", k=1.5*np.pi, degree=2) -u_target_function = dl.interpolate(u_target, Vh_STATE) -u_target_vector = u_target_function.vector() -qoi = soupy.L2MisfitControlQoI(mesh, Vh, u_target_vector) - -# 5. Define the ControlModel -control_model = soupy.ControlModel(pde, qoi) - -# 6. Choose the risk measure -risk_settings = soupy.meanVarRiskMeasureSAASettings() -risk_settings["beta"] = VARIANCE_WEIGHT -risk_settings["sample_size"] = SAMPLE_SIZE -risk_measure = soupy.MeanVarRiskMeasureSAA(control_model, prior, risk_settings) - -# 7. Define the penalization term -penalty = soupy.L2Penalization(Vh, PENALTY_WEIGHT) - -# 8. Assemble the cost functional -cost_functional = soupy.RiskMeasureControlCostFunctional(risk_measure, penalty) - -# 9. Define the optimizer -optimizer = soupy.BFGS(cost_functional) - -# 10. Provide initial guess and solve -print("Starting optimization") -z = cost_functional.generate_vector(soupy.CONTROL) -optimizer.solve(z) -print("Done") - -# 11. Post processing -risk_measure.computeComponents(z, order=0) -estimate_risk = risk_measure.cost() - -# Saving output to disk -z_fun = hp.vector2Function(z, Vh[soupy.CONTROL]) -with dl.HDF5File(mesh.mpi_comm(), "%s/z_opt.h5" %(RESULTS_DIRECTORY), "w") as save_file: - save_file.write(z_fun, "control") - -x = cost_functional.generate_vector() -x[soupy.CONTROL].axpy(1.0, z) - -# Sample from the prior -noise = dl.Vector() -prior.init_vector(noise, "noise") -rng = hp.Random(seed=0) -rng.normal(1.0, noise) -prior.sample(noise, x[soupy.PARAMETER]) - -# solve the forward problem -control_model.solveFwd(x[soupy.STATE], x) - -# save some figures -os.makedirs(RESULTS_DIRECTORY, exist_ok=True) - -plt.figure() -hp.nb.plot(hp.vector2Function(x[soupy.CONTROL], Vh[soupy.CONTROL])) -plt.title("Optimal control") -plt.savefig("%s/optimal_control.png" %(RESULTS_DIRECTORY)) -plt.close() - -plt.figure() -hp.nb.plot(hp.vector2Function(x[soupy.PARAMETER], Vh[soupy.PARAMETER])) -plt.title("Sample parameter at optimal") -plt.savefig("%s/sample_parameter.png" %(RESULTS_DIRECTORY)) -plt.close() - -plt.figure() -hp.nb.plot(hp.vector2Function(x[soupy.STATE], Vh[soupy.STATE])) -plt.title("Sample state at optimal") -plt.savefig("%s/sample_state.png" %(RESULTS_DIRECTORY)) -plt.close() - -plt.figure() -hp.nb.plot(u_target_function) -plt.title("Target state") -plt.savefig("%s/target_state.png" %(RESULTS_DIRECTORY)) -plt.close() - - - From 60dfb24820cda054c26d603619937fcaa8f17eff Mon Sep 17 00:00:00 2001 From: Dingcheng Date: Thu, 3 Aug 2023 16:05:50 -0500 Subject: [PATCH 6/9] updated examples --- ...rministic.py => driver_hyperelasticity.py} | 78 +++++++++++++++---- .../hyperelasticityControlPDE.py | 2 +- .../setupHyperelasticityProblem.py | 4 +- examples/poisson/driver_poisson_mean.py | 8 +- .../driver_semilinear_cvar.py | 7 +- .../semilinearEllipticControlPDE.py | 2 +- .../semilinearEllipticOUU.py | 2 +- 7 files changed, 80 insertions(+), 23 deletions(-) rename examples/hyperelasticity/{driver_hyperelasticity_deterministic.py => driver_hyperelasticity.py} (53%) diff --git a/examples/hyperelasticity/driver_hyperelasticity_deterministic.py b/examples/hyperelasticity/driver_hyperelasticity.py similarity index 53% rename from examples/hyperelasticity/driver_hyperelasticity_deterministic.py rename to examples/hyperelasticity/driver_hyperelasticity.py index ae59613..8ee6d40 100644 --- a/examples/hyperelasticity/driver_hyperelasticity_deterministic.py +++ b/examples/hyperelasticity/driver_hyperelasticity.py @@ -9,20 +9,43 @@ # # SOUPy is free software; you can redistribute it and/or modify it under the # terms of the GNU General Public License (as published by the Free -# Software Foundation) version 3.0 dated June 1991. +# Software Foundation) version 3.0 dated June 2007. """ -This example implements a deterministic optimization problem for material -design in hyperelasticity. The objective is to minimize the compliance of -the structure subject. The design variable takes values in [0, 1], which -selects between two different material properties. +This example implements an optimization problem for material +design of a hyperelastic beam. The objective is to minimize the compliance of +the structure subject to uncertain loading. + +The design variable takes values in [0, 1], which selects between +two different material properties. A nominal external load +is prescribed to be at the center of the beam, but is scaled by +a multiplicative Gaussian random field. + +- Set risk measure flag :code:`-r` to :code:`deterministic` + for Deterministic optimization using the mean parameter value + +- Set risk measure flag :code:`-r` to :code:`mean_var` for optimization of + the mean + variance risk measure + +- Flag :code:`beta` controls the variance weighting + +- Flag :code:`penalization` controls the l2 penalization on the design variable See :code:`setupHyperelasticityProblem.py` for problem settings, including mesh, geometry, and solver properties. The example uses a custom forward solver. See :code:`hyperelasticityControlPDE.py` for the PDE definition and solver + +This example also shows how to use the :code:`ScipyCostWrapper` to convert +a :code:`ControlCostFunctional` to be compatible with :code:`scipy.optimize` + +This driver supports MPI to parallelize the sampling of the parameter field. +e.g. + +mpirun -n 4 python driver_hyperelasticity.py + """ @@ -53,64 +76,91 @@ if __name__ == "__main__": parser = argparse.ArgumentParser() - parser.add_argument('-n', '--maxiter', type=int, default=100, help="Maximum number of SD iterations") + parser.add_argument('-r', '--risk_measure', type=str, default="deterministic", choices=["deterministic", "mean_var"], help="Risk measure type") + parser.add_argument('-n', '--sample_size', type=int, default=32, help="Sample size for risk measure computation") + parser.add_argument('-b', '--beta', type=float, default=1.0, help="Variance weight for risk measure") parser.add_argument('-q', '--qoi_type', type=str, default="stiffness", choices=["all", "stiffness", "point"]) parser.add_argument('-p', '--penalization', type=float, default=1e-1, help="Scaling of penalization") + + parser.add_argument('--maxiter', type=int, default=100, help="Maximum number of SD iterations") parser.add_argument('--show', default=False, action="store_true", help="Show figure") args = parser.parse_args() + save_dir = "results_%s" %(args.risk_measure) + os.makedirs(save_dir, exist_ok=True) + # Create mesh and setup model components - comm_mesh = MPI.COMM_WORLD + comm_mesh = MPI.COMM_SELF + comm_sampler = MPI.COMM_WORLD settings = hyperelasticity_problem_settings() settings["qoi_type"] = args.qoi_type mesh, Vh, hyperelasticity_varf, control_model, prior = setup_hyperelasticity_problem(settings, comm_mesh) l2_penalty = soupy.L2Penalization(Vh, args.penalization) - pde_cost = soupy.DeterministicControlCostFunctional(control_model, prior, l2_penalty) + if args.risk_measure == "deterministic": + pde_cost = soupy.DeterministicControlCostFunctional(control_model, prior, l2_penalty) + + else: + # Use the mean variance risk measure to assemble cost + risk_settings = soupy.meanVarRiskMeasureSAASettings() + risk_settings["beta"] = args.beta + risk_settings["sample_size"] = args.sample_size + risk_measure = soupy.MeanVarRiskMeasureSAA(control_model, prior, risk_settings, comm_sampler=comm_sampler) + pde_cost = soupy.RiskMeasureControlCostFunctional(risk_measure, l2_penalty) + + # ------------------ Using the scipy cost interface ------------------ # scipy_cost = soupy.ScipyCostWrapper(pde_cost, verbose=True) box_bounds = scipy.optimize.Bounds(lb=0.0, ub=1.0) + # Solve the PDE for the beam using entirely soft material x = control_model.generate_vector() x[soupy.PARAMETER].axpy(1.0, prior.mean) control_model.solveFwd(x[soupy.STATE], x) disp_fun_init = hp.vector2Function(x[soupy.STATE], Vh[soupy.STATE]) + # Solve optimal design problem using 0.5 as initial guess z0_np = x[soupy.CONTROL].get_local() + 0.5 options = {'gtol' : 1e-12, 'disp' : True, 'maxiter' : args.maxiter} results = scipy.optimize.minimize(scipy_cost.function(), z0_np, method="L-BFGS-B", jac=scipy_cost.jac(), bounds=box_bounds, options=options) + + # Optimal design z_opt = results['x'] + # --------------------------------------------------------------------- # + x[soupy.CONTROL].set_local(z_opt) x[soupy.CONTROL].apply("") control_model.solveFwd(x[soupy.STATE], x) z_fun = dl.Function(Vh[soupy.CONTROL], x[soupy.CONTROL]) - - with dl.HDF5File(mesh.mpi_comm(), "z_opt.h5", "w") as save_file: - save_file.write(z_fun, "control") + # ---------- postprocessing the plotting ----------- # + plt.figure() ax = dl.plot(disp_fun_init, mode="displacement") plt.colorbar(ax) plt.title("Displacement with soft material") - plt.savefig("u_soft.png") + plt.savefig("%s/u_soft.png" %(save_dir)) plt.figure() ax = dl.plot(hp.vector2Function(x[soupy.CONTROL], Vh[soupy.CONTROL])) plt.colorbar(ax) plt.title("Optimal design") - plt.savefig("z_opt.png") + plt.savefig("%s/z_opt.png" %(save_dir)) plt.figure() ax = dl.plot(hp.vector2Function(x[soupy.STATE], Vh[soupy.STATE]), mode="displacement") plt.colorbar(ax) plt.title("Optimal displacement") - plt.savefig("u_opt.png") + plt.savefig("%s/u_opt.png" %(save_dir)) + + if comm_sampler.Get_rank() == 0: + np.save('%s/z_opt.npy' %(save_dir), z_opt) if args.show: plt.show() diff --git a/examples/hyperelasticity/hyperelasticityControlPDE.py b/examples/hyperelasticity/hyperelasticityControlPDE.py index 0355284..a709026 100644 --- a/examples/hyperelasticity/hyperelasticityControlPDE.py +++ b/examples/hyperelasticity/hyperelasticityControlPDE.py @@ -9,7 +9,7 @@ # # SOUPy is free software; you can redistribute it and/or modify it under the # terms of the GNU General Public License (as published by the Free -# Software Foundation) version 3.0 dated June 1991. +# Software Foundation) version 3.0 dated June 2007. import pickle import sys, os diff --git a/examples/hyperelasticity/setupHyperelasticityProblem.py b/examples/hyperelasticity/setupHyperelasticityProblem.py index 262525f..85020e9 100644 --- a/examples/hyperelasticity/setupHyperelasticityProblem.py +++ b/examples/hyperelasticity/setupHyperelasticityProblem.py @@ -9,7 +9,7 @@ # # SOUPy is free software; you can redistribute it and/or modify it under the # terms of the GNU General Public License (as published by the Free -# Software Foundation) version 3.0 dated June 1991. +# Software Foundation) version 3.0 dated June 2007. import pickle import sys, os @@ -67,7 +67,7 @@ def hyperelasticity_problem_settings(): settings["geometry"] = {"lx" : 2.0, "ly" : 0.5, "lz" : 0.25, "dim" : 2} settings["mesh"] = {"nx" : 96, "ny" : 24, "nz" : 12} settings["solver"] = {"backtrack" : True, "load_steps" : [0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 0.9375, 1.0]} - settings["uncertainty"] = {"gamma" : 0.1, "delta" : 1, "robin_bc" : True} + settings["uncertainty"] = {"gamma" : 0.2, "delta" : 1.0, "robin_bc" : True} return settings diff --git a/examples/poisson/driver_poisson_mean.py b/examples/poisson/driver_poisson_mean.py index bd63dd0..1ad0020 100644 --- a/examples/poisson/driver_poisson_mean.py +++ b/examples/poisson/driver_poisson_mean.py @@ -9,14 +9,16 @@ # # SOUPy is free software; you can redistribute it and/or modify it under the # terms of the GNU General Public License (as published by the Free -# Software Foundation) version 3.0 dated June 1991. +# Software Foundation) version 3.0 dated June 2007. """ Minimization of the mean risk measure using the SAA approximation for the linear poisson problem with a log-normal conductivity field. -This driver uses MPI to parallelize the sampling of the parameter field. -For a purely serial implementation, see :code:`driver_poisson_mean_serial.py` +This driver supports MPI to parallelize the sampling of the parameter field. +e.g. + +mpirun -n 4 python driver_poisson_mean.py """ import sys diff --git a/examples/semilinear_elliptic/driver_semilinear_cvar.py b/examples/semilinear_elliptic/driver_semilinear_cvar.py index a699b02..2f0fec3 100644 --- a/examples/semilinear_elliptic/driver_semilinear_cvar.py +++ b/examples/semilinear_elliptic/driver_semilinear_cvar.py @@ -9,7 +9,7 @@ # # SOUPy is free software; you can redistribute it and/or modify it under the # terms of the GNU General Public License (as published by the Free -# Software Foundation) version 3.0 dated June 1991. +# Software Foundation) version 3.0 dated June 2007. """ Minimization of the CVaR risk measure using the SAA approximation for a @@ -29,6 +29,11 @@ This example also shows how to use the :code:`ScipyCostWrapper` to convert a :code:`ControlCostFunctional` to be compatible with :code:`scipy.optimize` + +This driver supports MPI to parallelize the sampling of the parameter field. +e.g. + +mpirun -n 4 python driver_semilinear_cvar.py """ import random diff --git a/examples/semilinear_elliptic/semilinearEllipticControlPDE.py b/examples/semilinear_elliptic/semilinearEllipticControlPDE.py index 8c01984..b335ed0 100644 --- a/examples/semilinear_elliptic/semilinearEllipticControlPDE.py +++ b/examples/semilinear_elliptic/semilinearEllipticControlPDE.py @@ -9,7 +9,7 @@ # # SOUPy is free software; you can redistribute it and/or modify it under the # terms of the GNU General Public License (as published by the Free -# Software Foundation) version 3.0 dated June 1991. +# Software Foundation) version 3.0 dated June 2007. import sys, os import math diff --git a/examples/semilinear_elliptic/semilinearEllipticOUU.py b/examples/semilinear_elliptic/semilinearEllipticOUU.py index d8b9e80..1ce3526 100644 --- a/examples/semilinear_elliptic/semilinearEllipticOUU.py +++ b/examples/semilinear_elliptic/semilinearEllipticOUU.py @@ -9,7 +9,7 @@ # # SOUPy is free software; you can redistribute it and/or modify it under the # terms of the GNU General Public License (as published by the Free -# Software Foundation) version 3.0 dated June 1991. +# Software Foundation) version 3.0 dated June 2007. import time From 7958cd9cd5541382ab16dc209c17de85f9392ce1 Mon Sep 17 00:00:00 2001 From: Dingcheng Date: Thu, 3 Aug 2023 16:06:33 -0500 Subject: [PATCH 7/9] updated running script for examples ci --- .github/workflows/run_examples.sh | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/.github/workflows/run_examples.sh b/.github/workflows/run_examples.sh index efff82e..e3aa810 100644 --- a/.github/workflows/run_examples.sh +++ b/.github/workflows/run_examples.sh @@ -5,8 +5,7 @@ export HIPPYLIB_PATH=$(pwd)/hippylib # Run poisson example cd examples/poisson -python3 driver_poisson_mean_serial.py +python3 driver_poisson_mean.py -mpirun -n 2 python3 driver_poisson_mean_mpi.py +mpirun -n 2 python3 driver_poisson.py -python3 compare_results.py From d23cff8280b6f9d3f7a2eb169d2b5ebd9e6980af Mon Sep 17 00:00:00 2001 From: Dingcheng Date: Thu, 3 Aug 2023 16:20:37 -0500 Subject: [PATCH 8/9] naming changes --- soupy/solver/newtonBacktrackSolver.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/soupy/solver/newtonBacktrackSolver.py b/soupy/solver/newtonBacktrackSolver.py index 35551c7..08dad30 100644 --- a/soupy/solver/newtonBacktrackSolver.py +++ b/soupy/solver/newtonBacktrackSolver.py @@ -22,7 +22,7 @@ import hippylib as hp -def applyBC(bcs, u): +def apply_BC(bcs, u): """ Applies a single or a list of :code:`dolfin.DirichletBC` to the vector :code:`u` @@ -35,7 +35,7 @@ def applyBC(bcs, u): else: bcs.apply(u) -def homogenizeBC(bcs): +def homogenize_BC(bcs): """ Converts a single or a list of :code:`dolfin.DirichletBC` to their homogenized (zeroed) forms @@ -131,14 +131,14 @@ def solve(self, residual_form, u, bcs=None, J_form = None, energy_form=None, sol my_rank = MPI.COMM_WORLD.Get_rank() # Apply BC to initial guess - applyBC(bcs, u.vector()) + apply_BC(bcs, u.vector()) # Homogenize the dirichlet BCs if they exist - bcs0 = homogenizeBC(bcs) + bcs0 = homogenize_BC(bcs) # Assemble initial residual residual = dl.assemble(residual_form) - # applyBC(bcs0, residual) + # apply_BC(bcs0, residual) r_norm = residual.norm("l2") tol = max(r_norm*rtol, atol) if energy_form is not None: @@ -204,7 +204,7 @@ def solve(self, residual_form, u, bcs=None, J_form = None, energy_form=None, sol else: # Backtrack on residual norm residual_new = dl.assemble(residual_form) - applyBC(bcs0, residual_new) + apply_BC(bcs0, residual_new) r_new_norm = residual_new.norm('l2') if r_new_norm < r_norm: backtrack_converged = True @@ -226,7 +226,7 @@ def solve(self, residual_form, u, bcs=None, J_form = None, energy_form=None, sol # Compute residual norm residual = dl.assemble(residual_form) - applyBC(bcs0, residual) + apply_BC(bcs0, residual) r_norm = residual.norm("l2") # Print new residual norms From c2d1ff18b5b9b5a83b5a8fc2b00766786b59097e Mon Sep 17 00:00:00 2001 From: Dingcheng Date: Thu, 3 Aug 2023 16:26:30 -0500 Subject: [PATCH 9/9] fixed typo in ci run example driver --- .github/workflows/run_examples.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/run_examples.sh b/.github/workflows/run_examples.sh index e3aa810..0781a27 100644 --- a/.github/workflows/run_examples.sh +++ b/.github/workflows/run_examples.sh @@ -7,5 +7,5 @@ cd examples/poisson python3 driver_poisson_mean.py -mpirun -n 2 python3 driver_poisson.py +mpirun -n 2 python3 driver_poisson_mean.py