1 // @HEADER 2 // ************************************************************************ 3 // 4 // Rapid Optimization Library (ROL) Package 5 // Copyright (2014) Sandia Corporation 6 // 7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive 8 // license for use of this work by or on behalf of the U.S. Government. 9 // 10 // Redistribution and use in source and binary forms, with or without 11 // modification, are permitted provided that the following conditions are 12 // met: 13 // 14 // 1. Redistributions of source code must retain the above copyright 15 // notice, this list of conditions and the following disclaimer. 16 // 17 // 2. Redistributions in binary form must reproduce the above copyright 18 // notice, this list of conditions and the following disclaimer in the 19 // documentation and/or other materials provided with the distribution. 20 // 21 // 3. Neither the name of the Corporation nor the names of the 22 // contributors may be used to endorse or promote products derived from 23 // this software without specific prior written permission. 24 // 25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY 26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE 29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 36 // 37 // Questions? Contact lead developers: 38 // Drew Kouri (dpkouri@sandia.gov) and 39 // Denis Ridzal (dridzal@sandia.gov) 40 // 41 // ************************************************************************ 42 // @HEADER 43 44 #ifndef ROL_PD_RANDVARFUNCTIONAL_HPP 45 #define ROL_PD_RANDVARFUNCTIONAL_HPP 46 47 #include "ROL_RandVarFunctional.hpp" 48 49 namespace ROL { 50 51 template<class Real> 52 class PD_RandVarFunctional : public RandVarFunctional<Real> { 53 typedef typename std::vector<Real>::size_type uint; 54 private: 55 Real pen_; 56 int update_; 57 bool setData_; 58 59 Ptr<SampledScalar<Real>> values_; 60 Ptr<SampledScalar<Real>> multipliers_; 61 Ptr<SampledScalar<Real>> multipliers_new_; 62 63 protected: 64 // Set value data at current parameter setValue(const Real val,const std::vector<Real> & pt)65 void setValue(const Real val, const std::vector<Real> &pt) { 66 values_->set(val, pt); 67 } 68 69 // Get multiplier at current parameter getMultiplier(Real & lam,const std::vector<Real> & pt) const70 void getMultiplier(Real &lam, const std::vector<Real> &pt) const { 71 multipliers_->get(lam, pt); 72 } 73 setMultiplier(Real & lam,const std::vector<Real> & pt)74 void setMultiplier(Real &lam, const std::vector<Real> &pt) { 75 multipliers_new_->set(lam, pt); 76 } 77 78 // Get penalty parameter getPenaltyParameter(void) const79 Real getPenaltyParameter(void) const { 80 return pen_; 81 } 82 83 // Smooth plus function approximation ppf(const Real x,const Real t,const Real r,const int deriv=0) const84 Real ppf(const Real x, const Real t, const Real r, const int deriv = 0) const { 85 const Real zero(0), half(0.5), one(1), arg(r*x+t); 86 Real val(0); 87 if ( arg < zero ) { 88 val = (deriv==0 ? -half*t*t/r : zero); 89 } 90 else if ( zero <= arg && arg <= one ) { 91 val = (deriv==0 ? half*r*x*x+t*x 92 : (deriv==1 ? arg : r)); 93 } 94 else { 95 val = (deriv==0 ? (arg-half*(t*t+one))/r 96 : (deriv==1 ? one : zero)); 97 } 98 return val; 99 } 100 101 public: PD_RandVarFunctional(void)102 PD_RandVarFunctional(void) 103 : RandVarFunctional<Real>(), pen_(1.0), update_(0), setData_(true) { 104 values_ = makePtr<SampledScalar<Real>>(); 105 multipliers_ = makePtr<SampledScalar<Real>>(); 106 multipliers_new_ = makePtr<SampledScalar<Real>>(); 107 } 108 setData(SampleGenerator<Real> & sampler,const Real pen,const Real lam=0.0)109 void setData(SampleGenerator<Real> &sampler, const Real pen, const Real lam = 0.0) { 110 if (setData_) { 111 pen_ = pen; 112 for (int i = sampler.start(); i < sampler.numMySamples(); ++i) { 113 multipliers_->set(lam, sampler.getMyPoint(i)); 114 } 115 setData_ = false; 116 } 117 } 118 computeDual(SampleGenerator<Real> & sampler)119 virtual Real computeDual(SampleGenerator<Real> &sampler) { 120 const Real zero(0), one(1); 121 Real val(0), lold(0), lnew(0), mdiff(0), gdiff(0); 122 for (int i = sampler.start(); i < sampler.numMySamples(); ++i) { 123 values_->get(val, sampler.getMyPoint(i)); 124 multipliers_->get(lold, sampler.getMyPoint(i)); 125 if (update_ == 0) { 126 //lnew = ppf(val, lold, pen_, 1); 127 lnew = std::min(one, std::max(zero, pen_*val+lold)); 128 } 129 else { 130 lnew = (val < zero ? zero : one); 131 } 132 mdiff += sampler.getMyWeight(i) * std::pow(lnew-lold,2); 133 multipliers_new_->set(lnew, sampler.getMyPoint(i)); 134 } 135 sampler.sumAll(&mdiff,&gdiff,1); 136 gdiff = std::sqrt(gdiff); 137 return gdiff; 138 } 139 updateDual(SampleGenerator<Real> & sampler)140 void updateDual(SampleGenerator<Real> &sampler) { 141 Real lam(0); 142 for (int i = sampler.start(); i < sampler.numMySamples(); ++i) { 143 multipliers_new_->get(lam, sampler.getMyPoint(i)); 144 multipliers_->set(lam, sampler.getMyPoint(i)); 145 } 146 } 147 updatePenalty(const Real pen)148 void updatePenalty(const Real pen) { 149 pen_ = pen; 150 } 151 setStorage(const Ptr<SampledScalar<Real>> & value_storage,const Ptr<SampledVector<Real>> & gradient_storage)152 virtual void setStorage(const Ptr<SampledScalar<Real>> &value_storage, 153 const Ptr<SampledVector<Real>> &gradient_storage) { 154 RandVarFunctional<Real>::setStorage(value_storage,gradient_storage); 155 } 156 setHessVecStorage(const Ptr<SampledScalar<Real>> & gradvec_storage,const Ptr<SampledVector<Real>> & hessvec_storage)157 virtual void setHessVecStorage(const Ptr<SampledScalar<Real>> &gradvec_storage, 158 const Ptr<SampledVector<Real>> &hessvec_storage) { 159 RandVarFunctional<Real>::setHessVecStorage(gradvec_storage,hessvec_storage); 160 } 161 initialize(const Vector<Real> & x)162 virtual void initialize(const Vector<Real> &x) { 163 RandVarFunctional<Real>::initialize(x); 164 } 165 }; 166 167 } 168 169 #endif 170