fix the stochastic collocation with gamma distribution

doc/tests/analytic_tests.tex

@@ -303,6 +303,7 @@
 \input{projectile.tex}
 \input{attenuate.tex}
 \input{tensor_poly.tex}
+\input{gamma_scgpc.tex}
 \input{sobol_sens.tex}
 \input{rims.tex}
 \input{parabolas.tex}

doc/tests/gamma_scgpc.tex

@@ -0,0 +1,58 @@
+\section{Stochastic Collocation with Gamma Distribution}
+Associated external model: \texttt{poly\_scgpc\_gamma.py}
+
+Recall that the \textit{Gamma} distribution has the probability density function
+\begin{align}
+  f(x) = \frac{x^{\alpha-1}e^{-x/\beta}}{\beta^\alpha\Gamma\left(\alpha\right)}, \alpha > 0, \beta > 0
+\end{align}
+
+The following two polynomials are used to compute the analytic statistical moments:
+\textit{Gamma} distribution:
+\begin{align}
+  u_1(x, y) = x + y \\
+  u_2(x, y) = x^2 + y^2
+\end{align}
+where $x$ and $y$ are two mutually independent \textit{Gamma} variates, i.e.
+\begin{align}
+  x \thicksim \Gamma\left(\alpha_{1},\beta_{1}\right) \notag \\
+  y \thicksim \Gamma\left(\alpha_{2},\beta_{2}\right) \notag
+\end{align}
+
+\subsection{Mean and Variance}
+The first two statistical moments of $u_1(x,y)$ and $u_2{x,y}$ are:
+\begin{align}
+  \expv{u_1(x,y)} &= \int_{0}^\infty dxdyP\left(x,y\right)u_1\left(x,y\right), \notag \\
+  &= \int_{0}^\infty dxdy \Gamma\left(\alpha_{1},\beta_{1}\right) \Gamma\left(\alpha_{2},\beta_{2}\right) u_1\left(x,y\right), \notag \\
+  &= \frac{\alpha_{1}}{\beta_{1}} + \frac{\alpha_{2}}{\beta_{2}}
+\end{align}
+\begin{align}
+  \expv{u_2(x,y)} &= \int_{0}^\infty dxdyP\left(x,y\right)u_2\left(x,y\right), \notag \\
+  &= \int_{0}^\infty dxdy \Gamma\left(\alpha_{1},\beta_{1}\right) \Gamma\left(\alpha_{2},\beta_{2}\right) u_2\left(x,y\right), \notag \\
+  &= \frac{\left(\alpha_{1} + 1\right)\alpha_1}{\beta_{1}^2} + \frac{\left(\alpha_{2}+1\right)\alpha_2}{\beta_{2}^2}
+\end{align}
+
+\begin{align}
+  \text{var}[u_1(x,y)] &= \int_{0}^\infty dxdyP\left(x,y\right)\left[u_1\left(x,y\right) - \expv{u_1(x,y)}\right]^2, \notag \\
+  &= \frac{\alpha_{1}}{\beta_{1}^2} + \frac{\alpha_{2}}{\beta_{2}^2}
+\end{align}
+\begin{align}
+  \text{var}[u_2(x,y)] &= \int_{0}^\infty dxdyP\left(x,y\right)\left[u_2\left(x,y\right) - \expv{u_2(x,y)}\right]^2, \notag \\
+  &= \frac{\left(4\alpha_{1} + 6.0\right)\left(\alpha_{1} + 1\right)\alpha_1}{\beta_{1}^4} + \frac{\left(4\alpha_{2} + 6.0\right)\left(\alpha_{2}+1\right)\alpha_2}{\beta_{2}^4}
+\end{align}
+
+\subsection{numeric values}
+Some numeric values for the mean and variance are listed below for given distributions:
+\begin{align}
+  x \thicksim \Gamma\left(11, 5\right) \notag \\
+  y \thicksim \Gamma\left(2, 0.8\right) \notag
+\end{align}
+\begin{table}[h!]
+  \centering
+  \begin{tabular}{c|c|c}
+    $Function$ & mean & variance \\ \hline
+    $u_1$ & 4.7    & 3.565 \\
+    $u_2$ & 14.655 & 215.638125 \\
+  \end{tabular}
+\end{table}
+
+

framework/Distributions.py

@@ -877,7 +877,6 @@ def stdProbabilityNorm(self):
     @ In, None, None
     @ Out float, norm
     """
-    #return self.beta**self.alpha/factorial(self.alpha-1.)
     return 1./factorial(self.alpha-1)
 
 DistributionsCollection.addSub(Gamma.getInputSpecification())

framework/OrthoPolynomials.py

@@ -377,7 +377,7 @@ def norm(self,order):
       @ In, order, int, polynomial order to get norm of
       @ Out, norm, float, value of poly norm
     """
-    return np.sqrt(gamma(1.0+order)/gamma(1.0+order+self.params[0]))
+    return np.sqrt(gamma(1.0+order)*gamma(1.0+self.params[0])/gamma(1.0+order+self.params[0]))
 
 
 

tests/framework/AnalyticModels/poly_scgpc_gamma.py

@@ -0,0 +1,30 @@
+# Copyright 2017 Battelle Energy Alliance, LLC
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+# http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#***************************************
+#* Simple analytic test ExternalModule *
+#***************************************
+# This model tests the stochastic collocation with gamma distribution
+
+import numpy as np
+
+def eval(inp,exp):
+  return sum(n**exp for n in inp)
+
+def run(self,Input):
+  self.ans = eval((self.x1,self.x2),1.0)
+  self.ans2 = eval((self.x1,self.x2),2.0)
+
+#
+#  This model has analytic mean and variance documented in raven/docs/tests
+#

tests/framework/Samplers/SparseGrid/gold/scgpc/dumpGamma.csv

@@ -0,0 +1,10 @@
+x1,x2,ans,ans2
+1.42393039323,1.16977778441,2.59370817764,3.39595782966
+1.42393039323,4.13175911167,5.5556895049,19.0990111216
+1.42393039323,9.69846310393,11.1223934972,96.0877643431
+2.46509275174,1.16977778441,3.63487053615,7.44506233958
+2.46509275174,4.13175911167,6.59685186341,23.1481156315
+2.46509275174,9.69846310393,12.1635558557,100.136868853
+3.91097685502,1.16977778441,5.08075463943,16.6641200254
+3.91097685502,4.13175911167,8.04273596669,32.3671733174
+3.91097685502,9.69846310393,13.609439959,109.355926539

tests/framework/Samplers/SparseGrid/gold/scgpc/dumpromGamma.csv

@@ -0,0 +1,10 @@
+x1,x2,ans,ans2
+1.42393039323,1.16977778441,2.59370817764,3.39595782966
+1.42393039323,4.13175911167,5.5556895049,19.0990111216
+1.42393039323,9.69846310393,11.1223934972,96.0877643431
+2.46509275174,1.16977778441,3.63487053615,7.44506233958
+2.46509275174,4.13175911167,6.59685186341,23.1481156315
+2.46509275174,9.69846310393,12.1635558557,100.136868853
+3.91097685502,1.16977778441,5.08075463943,16.6641200254
+3.91097685502,4.13175911167,8.04273596669,32.3671733174
+3.91097685502,9.69846310393,13.609439959,109.355926539

tests/framework/Samplers/SparseGrid/gold/scgpc/dumpromGammaXml.xml

@@ -0,0 +1,71 @@
+<ROM type="Static">
+  <type>GAUSSgpcROM(ans-ans2)</type>
+
+  <ans>
+    <mean>4.7</mean>
+    <expectedValue>4.7</expectedValue>
+    <variance>3.565</variance>
+    <samples>9</samples>
+    <polyCoeffs>
+      <inputVariables>x1,x2</inputVariables>
+      <_0_0_>4.7</_0_0_>
+      <_0_1_>-1.76776695297</_0_1_>
+      <_0_2_>2.35794922682e-15</_0_2_>
+      <_1_0_>-0.663324958071</_1_0_>
+      <_1_1_>6.4161883723e-17</_1_1_>
+      <_1_2_>-3.52890360476e-16</_1_2_>
+      <_2_0_>2.41409087508e-15</_2_0_>
+      <_2_1_>-7.37861662814e-16</_2_1_>
+      <_2_2_>6.4161883723e-17</_2_2_>
+    </polyCoeffs>
+    <partialVariance>
+      <x2>3.125</x2>
+      <x1>0.44</x1>
+      <x1.x2>6.77204934613e-31</x1.x2>
+    </partialVariance>
+    <sobolIndices>
+      <x2>0.876577840112</x2>
+      <x1>0.123422159888</x1>
+      <x1.x2>1.89959308447e-31</x1.x2>
+    </sobolIndices>
+    <sobolTotalIndices>
+      <x2>0.876577840112</x2>
+      <x1>0.123422159888</x1>
+      <x1.x2>1.89959308447e-31</x1.x2>
+    </sobolTotalIndices>
+  </ans>
+
+  <ans2>
+    <mean>14.655</mean>
+    <expectedValue>14.655</expectedValue>
+    <variance>215.638125</variance>
+    <samples>9</samples>
+    <polyCoeffs>
+      <inputVariables>x1,x2</inputVariables>
+      <_0_0_>14.655</_0_0_>
+      <_0_1_>-13.2582521472</_0_1_>
+      <_0_2_>5.41265877365</_0_2_>
+      <_1_0_>-3.18395979874</_1_0_>
+      <_1_1_>1.02659013957e-15</_1_1_>
+      <_1_2_>-1.53988520935e-15</_1_2_>
+      <_2_0_>0.649923072371</_2_0_>
+      <_2_1_>-5.90289330251e-15</_2_1_>
+      <_2_2_>1.79653274424e-15</_2_2_>
+    </polyCoeffs>
+    <partialVariance>
+      <x2>205.078125</x2>
+      <x1>10.56</x1>
+      <x1.x2>4.14968130146e-29</x1.x2>
+    </partialVariance>
+    <sobolIndices>
+      <x2>0.951029067796</x2>
+      <x1>0.0489709322041</x1>
+      <x1.x2>1.92437274321e-31</x1.x2>
+    </sobolIndices>
+    <sobolTotalIndices>
+      <x2>0.951029067796</x2>
+      <x1>0.0489709322041</x1>
+      <x1.x2>1.92437274321e-31</x1.x2>
+    </sobolTotalIndices>
+  </ans2>
+</ROM>

tests/framework/Samplers/SparseGrid/scgpc/polynomial.py

@@ -15,7 +15,6 @@
 
 def eval(inp,exp):
   return sum(n**exp for n in inp)
-  #return np.exp(-sum(inp)/len(inp))
 
 def run(self,Input):
   self.ans = eval((self.x1,self.x2),1.0)

tests/framework/Samplers/SparseGrid/test_scgpc_gamma_large_scale.xml

@@ -0,0 +1,127 @@
+<?xml version="1.0" ?>
+<Simulation printTimeStamps="false">
+  <RunInfo>
+    <WorkingDir>scgpc</WorkingDir>
+    <Sequence>makeG,trainG,testG,printG,printRomG</Sequence>
+    <batchSize>1</batchSize>
+  </RunInfo>
+  <TestInfo>
+    <name>framework/Samplers/SparseGrid.gamma_large_scale</name>
+    <author>wangc</author>
+    <created>2018-07-26</created>
+    <classesTested>Samplers.SparseGridCollocation</classesTested>
+    <description>
+        This tests using SparseGridCollocation with gamma-distributed inputs with large scale parameters
+        to test the fix for issue #693
+    </description>
+    <analytic>
+        This test is analytic in mean and variance documented in the analytic tests.
+        Two models are used:
+        ans = x1 + x2 with analytic mean = 4.7, variance = 3.565
+        ans = x1**2 + x2**2 with analytic mean = 14.655, variance = 215.638125
+    </analytic>
+  </TestInfo>
+
+  <Distributions>
+    <Gamma name="GammaDist1">
+      <low>0</low>
+      <alpha>11</alpha>
+      <beta>5</beta>
+    </Gamma>
+    <Gamma name="GammaDist2">
+      <low>0</low>
+      <alpha>2</alpha>
+      <beta>0.8</beta>
+    </Gamma>
+  </Distributions>
+
+  <Samplers>
+    <SparseGridCollocation name="GSG" parallel="1">
+      <variable name="x1">
+        <distribution>GammaDist1</distribution>
+      </variable>
+      <variable name="x2">
+        <distribution>GammaDist2</distribution>
+      </variable>
+      <ROM class="Models" type="ROM">GROM</ROM>
+    </SparseGridCollocation>
+  </Samplers>
+
+  <Models>
+    <Dummy name="MyDummy" subType=""/>
+    <ExternalModel ModuleToLoad="../../../AnalyticModels/poly_scgpc_gamma" name="polynomial" subType="">
+      <variables>x1,x2,ans,ans2</variables>
+    </ExternalModel>
+    <ROM name="GROM" subType="GaussPolynomialRom">
+      <Target>ans,ans2</Target>
+      <Features>x1,x2</Features>
+      <IndexSet>TensorProduct</IndexSet>
+      <PolynomialOrder>2</PolynomialOrder>
+      <Interpolation poly="Laguerre" quad="Laguerre">x1</Interpolation>
+      <Interpolation poly="Laguerre" quad="Laguerre" weight="1">x2</Interpolation>
+    </ROM>
+  </Models>
+
+  <Steps>
+    <MultiRun name="makeG" pauseAtEnd="false">
+      <Input class="DataObjects" type="PointSet">dummyIN</Input>
+      <Model class="Models" type="ExternalModel">polynomial</Model>
+      <Sampler class="Samplers" type="SparseGridCollocation">GSG</Sampler>
+      <Output class="DataObjects" type="PointSet">Gsolns</Output>
+    </MultiRun>
+    <MultiRun name="testG" pauseAtEnd="false">
+      <Input class="DataObjects" type="PointSet">dummyIN</Input>
+      <Model class="Models" type="ROM">GROM</Model>
+      <Sampler class="Samplers" type="SparseGridCollocation">GSG</Sampler>
+      <Output class="DataObjects" type="PointSet">Gtests</Output>
+    </MultiRun>
+    <IOStep name="printRomG">
+      <Input class="DataObjects" type="PointSet">Gtests</Input>
+      <Input class="Models" type="ROM">GROM</Input>
+      <Output class="OutStreams" type="Print">dumpromGamma</Output>
+      <Output class="OutStreams" type="Print">dumpromGammaXml</Output>
+    </IOStep>
+    <IOStep name="printG">
+      <Input class="DataObjects" type="PointSet">Gsolns</Input>
+      <Output class="OutStreams" type="Print">dumpGamma</Output>
+    </IOStep>
+    <RomTrainer name="trainG">
+      <Input class="DataObjects" type="PointSet">Gsolns</Input>
+      <Output class="Models" type="ROM">GROM</Output>
+    </RomTrainer>
+  </Steps>
+
+  <DataObjects>
+    <PointSet name="dummyIN">
+      <Input>x1,x2</Input>
+      <Output>OutputPlaceHolder</Output>
+    </PointSet>
+    <PointSet name="Gtests">
+      <Input>x1,x2</Input>
+      <Output>ans,ans2</Output>
+    </PointSet>
+    <PointSet name="Gsolns">
+      <Input>x1,x2</Input>
+      <Output>ans,ans2</Output>
+    </PointSet>
+  </DataObjects>
+
+  <OutStreams>
+    <Print name="dumpromGamma">
+      <type>csv</type>
+      <source>Gtests</source>
+      <what>input,output</what>
+    </Print>
+    <Print name="dumpGamma">
+      <type>csv</type>
+      <source>Gsolns</source>
+      <what>input,output</what>
+    </Print>
+    <Print name="dumpromGammaXml">
+      <type>xml</type>
+      <source>GROM</source>
+      <what>all</what>
+    </Print>
+  </OutStreams>
+
+</Simulation>

tests/framework/Samplers/SparseGrid/tests

@@ -63,6 +63,12 @@
   UnorderedXml = 'scgpc/Gdumprom.xml'
   csv = 'scgpc/Gdump.csv'
  [../]
+ [./gamma_large_scale]
+  type = 'RavenFramework'
+  input = 'test_scgpc_gamma_large_scale.xml'
+  UnorderedXml = 'scgpc/dumpromGammaXml.xml'
+  csv = 'scgpc/dumpGamma.csv scgpc/dumpromGamma.csv'
+ [../]
  [./beta]
   type = 'RavenFramework'
   input = 'test_scgpc_beta.xml'