1# _______________________________________________________________________ 2# 3# PECOS: Parallel Environment for Creation Of Stochastics 4# Copyright (c) 2011, Sandia National Laboratories. 5# This software is distributed under the GNU Lesser General Public License. 6# For more information, see the README file in the top Pecos directory. 7# _______________________________________________________________________ 8 9import numpy 10from PyDakota.approximation import Function, PyFunction 11class GenzFunction(Function): 12 def __init__( self, func_type, num_vars, c=None, w=None, name=None ): 13 """ 14 having c and w as arguments are required to enable pickling 15 name allows name to be consistent when pickling 16 because name set here is appended to when set_coefficients is called 17 and that function is not called when pickling 18 """ 19 Function.__init__( self ) # must call c++ object wrapper init function 20 21 self.func_type = func_type 22 self.num_vars = num_vars 23 self.num_qoi = 1 24 if name is None: 25 self.name = 'genz-' + self.func_type + '-%d' %self.num_vars 26 else: 27 self.name=name 28 self.c=c 29 self.w=w 30 31 self.evaluator=PyFunction(self.value_) 32 33 def value(self, x, opts=dict()): 34 eval_type=opts.get('eval_type','value') 35 if ( ('grad' in eval_type) and 36 (self.func_type == "discontinuous" or self.func_type == "continuous") ): 37 msg = "gradients cannot be computed for %s Genz function"%self.func_type 38 raise Exception, msg 39 assert x.min()>=0 and x.max()<=1. 40 vals = self.evaluator.value(x) 41 if eval_type=='value': 42 return vals[:,:1] 43 if eval_type=='value-grad': 44 return vals 45 if eval_type=='grad': 46 return vals[:,1:] 47 48 def value_( self, x ): 49 assert x.ndim == 1 50 if (self.func_type == "discontinuous" or self.func_type == "continuous"): 51 result = numpy.empty( ( 1 ) , numpy.double ) 52 else: 53 result = numpy.empty( ( self.num_vars+1 ) , numpy.double ); 54 self.num_vars = self.c.shape[0] 55 if ( self.func_type == "oscillatory" ): 56 result[0] = 2.0 * numpy.pi * self.w[0] 57 for d in range( self.num_vars ): 58 result[0] += self.c[d] * x[d]; 59 for d in range( self.num_vars ): 60 result[1+d] = -self.c[d] * numpy.sin( result[0] ); 61 result[0] = numpy.cos( result[0] ); 62 elif ( self.func_type == "product-peak" ): 63 result[0] = 1.0 64 for d in range( self.num_vars ): 65 result[0] *= ( 1.0 / ( self.c[d] * self.c[d] ) + 66 ( x[d] - self.w[d] ) * ( x[d] - self.w[d] ) ) 67 for d in range( self.num_vars ): 68 result[1+d] = 2. * ( x[d] - self.w[d] ) / result[0] 69 result[0] = 1.0 / result[0] 70 elif ( self.func_type == "corner-peak" ): 71 result[0] = 1.0; 72 for d in range( self.num_vars ): 73 result[0] += self.c[d] * x[d]; 74 for d in range( self.num_vars ): 75 result[1+d] = -self.c[d] * (self.num_vars+1) / \ 76 result[0]**( self.num_vars+2 ) 77 result[0] = 1.0 / result[0]**( self.num_vars+1 ) 78 elif ( self.func_type == "gaussian-peak" ): 79 result[0] = 0.0; 80 for d in range( self.num_vars ): 81 result[0] += self.c[d] * self.c[d] * ( x[d] - self.w[d] ) * \ 82 ( x[d] - self.w[d] ) 83 for d in range( self.num_vars ): 84 result[1+d] = 2. * self.c[d]**2 * (self.w[d]-x[d]) * \ 85 numpy.exp( -result[0] ) 86 result[0] = numpy.exp( -result[0] ); 87 elif ( self.func_type == "continuous" ): 88 result[0] = 0.0; 89 for d in range( self.num_vars ): 90 result[0] += self.c[d] * numpy.abs( x[d] - self.w[d] ); 91 result[0] = numpy.exp( -result[0] ); 92 result[1:] = 0. 93 elif ( self.func_type == "discontinuous" ): 94 result[0] = 0.0; 95 if ( self.num_vars == 1 ): 96 if ( x[0] <= self.w[0] ): 97 for d in range( self.num_vars ): 98 result[0] += self.c[d] * x[d]; 99 result[0] = numpy.exp( result[0] ); 100 else: 101 if ( ( x[0] <= self.w[0] ) and ( x[1] <= self.w[1] ) ): 102 for d in range( self.num_vars ): 103 result[0] += self.c[d] * x[d]; 104 result[0] = numpy.exp( result[0] ); 105 elif (self.func_type=='corner-peak-second-order'): 106 result[0]=0. 107 for d in xrange(self.num_vars-1): 108 result[0] += (1.+self.c[d]*x[d] +self.c[d+1]*x[d+1])**(-3) 109 else: 110 msg = "ensure func_num in [\"oscillatory\",\"product-peak\"," 111 msg += "\"corner-peak\",\"gaussian-peak\"," 112 msg += "\"continuous\",\"discontinuous\"]" 113 raise Exception, msg 114 115 return result 116 117 def set_coefficients( self, c_factor, coef_type, w_factor = 0., seed=0 ): 118 119 self.c = numpy.empty( ( self.num_vars ), numpy.double ); 120 self.w = numpy.empty( ( self.num_vars ), numpy.double ); 121 self.name += '-' + coef_type 122 123 if ( coef_type == "no-decay" ): 124 csum = 0.0 125 for d in range( self.num_vars ): 126 self.w[d] = w_factor 127 self.c[d] = ( d + 0.5 ) / self.num_vars 128 csum += self.c[d] 129 130 self.c *= ( c_factor / csum ) 131 elif ( coef_type == "quadratic-decay" ): 132 csum = 0.0 133 for d in range( self.num_vars ): 134 self.w[d] = w_factor 135 self.c[d] = 1.0 / ( d + 1. )**2 136 csum += self.c[d] 137 for d in range( self.num_vars ): 138 self.c[d] *= ( c_factor / csum ) 139 elif ( coef_type == "quartic-decay" ): 140 csum = 0.0 141 for d in range( self.num_vars ): 142 self.w[d] = w_factor 143 self.c[d] = 1.0 / ( d + 1. )**4 144 csum += self.c[d] 145 for d in range( self.num_vars ): 146 self.c[d] *= ( c_factor / csum ) 147 elif ( coef_type == "squared-exponential-decay" ): 148 csum = 0.0 149 for d in range( self.num_vars ): 150 self.w[d] = w_factor 151 self.c[d]=numpy.exp((d+1)**2*numpy.log(1.e-15)/self.num_vars**2) 152 csum += self.c[d] 153 for d in range( self.num_vars ): 154 self.c[d] *= ( c_factor / csum ) 155 elif ( coef_type == "exponential-decay" ): 156 csum = 0.0 157 for d in range( self.num_vars ): 158 self.w[d] = w_factor 159 #self.c[d] = numpy.exp( (d+1)*numpy.log( 1.e-15 )/self.num_vars ) 160 self.c[d] = numpy.exp( (d+1)*numpy.log( 1.e-8 )/self.num_vars ) 161 csum += self.c[d] 162 for d in range( self.num_vars ): 163 self.c[d] *= ( c_factor / csum ) 164 elif (coef_type=='random'): 165 numpy.random.seed(seed) 166 self.c = numpy.random.uniform(0.,1.,(self.num_vars)) 167 csum = numpy.sum(self.c) 168 self.c *= c_factor/csum 169 #self.w = numpy.random.uniform( 0.,1.,(self.num_vars)) 170 self.w = 0. 171 else: 172 msg = "Ensure coef_type in [\"no-decay\",\"quadratic-decay\"" 173 msg += ",\"quartic-decay\", \"squared-exponential-decay,\"" 174 msg += ",\"exponential-decay\"]" 175 raise Exception, msg 176 177 def variance( self ): 178 mean = self.integrate() 179 return self.recursive_uncentered_moment(self.num_vars,0.0,self.num_vars+1)-mean**2 180 181 def recursive_uncentered_moment( self, d, integral_sum, r ): 182 if ( self.func_type == "corner-peak" ): 183 if ( d > 0 ): 184 return ( self.recursive_uncentered_moment( 185 d-1, integral_sum, r ) - 186 self.recursive_uncentered_moment( 187 d-1, integral_sum + self.c[d-1], r ) ) / \ 188 ( (d+r) * self.c[d-1] ) 189 else: 190 return 1.0 / ( 1.0 + integral_sum )**(1+r); 191 else: 192 return 0. 193 194 def integrate( self ): 195 return self.recursive_integrate( self.num_vars, 0.0 ) 196 197 def recursive_integrate( self, d, integral_sum ): 198 if ( self.func_type == "oscillatory" ): 199 if ( d > 0 ): 200 return ( self.recursive_integrate( d-1, 201 integral_sum+self.c[d-1] ) 202 - self.recursive_integrate( d-1, 203 integral_sum ) ) /\ 204 self.c[d-1]; 205 else: 206 case = self.num_vars % 4 207 if (case == 0 ): 208 return numpy.cos( 2.0 * numpy.pi * self.w[0] + integral_sum) 209 if (case == 1 ): 210 return numpy.sin( 2.0 * numpy.pi * self.w[0] + integral_sum) 211 if (case == 2 ): 212 return -numpy.cos( 2.0 * numpy.pi * self.w[0] + integral_sum) 213 if (case == 3 ): 214 return -numpy.sin( 2.0 * numpy.pi * self.w[0] + integral_sum) 215 216 elif ( self.func_type == "product-peak" ): 217 prod = 1.0 218 for i in xrange( self.num_vars ): 219 prod = prod * self.c[i] * ( numpy.arctan( self.c[i] * 220 (1.0 - self.w[i] ) ) + 221 numpy.arctan( self.c[i] * self.w[i])) 222 return prod 223 elif ( self.func_type == "corner-peak" ): 224 if ( d > 0 ): 225 return ( self.recursive_integrate( d-1, integral_sum ) - 226 self.recursive_integrate( d-1, integral_sum + 227 self.c[d-1] ) ) / \ 228 ( d * self.c[d-1] ) 229 else: 230 return 1.0 / ( 1.0 + integral_sum ); 231 232 elif ( self.func_type == "gaussian-peak" ): 233 msg = "GenzModel::recursive_integrate() integration of " 234 msg += "gaussian_peak function has not been implmemented."; 235 raise Exception, msg 236 237 elif ( self.func_type == "continuous" ): 238 prod = 1.0; 239 for i in xrange( self.num_vars ): 240 prod /= ( self.c[i] * ( 2.0 - numpy.exp( -self.c[i]*self.w[i])- 241 numpy.exp( self.c[i]*(self.w[i]-1.0)))); 242 return prod; 243 elif ( self.func_type == "discontinuous" ): 244 prod = 1.0; 245 if ( self.num_vars < 2 ): 246 for i in xrange( self.num_vars ): 247 prod *= ( numpy.exp( self.c[i] * self.w[i] )-1.0 )/self.c[i]; 248 else: 249 for i in xrange( 2 ): 250 prod *= ( numpy.exp( self.c[i] * self.w[i] ) - 1.0)/self.c[i] 251 for i in xrange( 2, self.num_vars ): 252 prod *= ( numpy.exp( self.c[i] ) - 1.0 ) / self.c[i]; 253 return prod; 254 elif (self.func_type=='corner-peak-second-order'): 255 self.func_type='corner-peak' 256 c = self.c.copy() 257 num_vars = self.num_vars 258 result = 0 259 for d in xrange(num_vars-1): 260 self.num_vars=2; self.c = c[d:d+2] 261 result += self.recursive_integrate(2,0.) 262 self.func_type='corner-peak-second-order' 263 self.c = c.copy() 264 self.num_vars = num_vars 265 return result 266 else: 267 msg = "GenzModel::recursive_integrate() incorrect f_type_ "; 268 msg += "was provided"; 269 raise Exception, msg 270 271 def __reduce__(self): 272 return (type(self),(self.func_type, self.num_vars, self.c, self.w, self.name)) 273 274# add unit test like to test pickling 275#python -c "from PyDakota.models.genz import GenzFunction; g = GenzFunction('oscillatory',2); import pickle; pickle.dump(g,open('function.pkl','wb')); g1 = pickle.load(open('function.pkl','rb'))" 276