--- /dev/null
+
+ASSIMILATION_STUDY(Study_name='ADAO skeleton case',
+ Study_repertory='@prefix@/share/salome/adao_examples/daSkeletons/External_data_definition_by_scripts',
+ Debug=0,
+ Algorithm='3DVAR',
+ Background=_F(INPUT_TYPE='Vector',
+ data=_F(FROM='Script',
+ SCRIPT_FILE='Script_Background_xb.py',),),
+ BackgroundError=_F(INPUT_TYPE='Matrix',
+ data=_F(FROM='Script',
+ SCRIPT_FILE='Script_BackgroundError_B.py',),),
+ Observation=_F(INPUT_TYPE='Vector',
+ data=_F(FROM='Script',
+ SCRIPT_FILE='Script_Observation_yo.py',),),
+ ObservationError=_F(INPUT_TYPE='Matrix',
+ data=_F(FROM='Script',
+ SCRIPT_FILE='Script_ObservationError_R.py',),),
+ ObservationOperator=_F(INPUT_TYPE='Function',
+ data=_F(FROM='FunctionDict',
+ FUNCTIONDICT_FILE='Script_ObservationOperator_H.py',),),
+ AlgorithmParameters=_F(INPUT_TYPE='Dict',
+ data=_F(FROM='Script',
+ SCRIPT_FILE='Script_AlgorithmParameters.py',),),
+ UserPostAnalysis=_F(FROM='Script',
+ SCRIPT_FILE='Script_UserPostAnalysis.py',),);
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2011 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+# Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
+
+__doc__ = """
+ ADAO skeleton case, for wide script usage in case definition
+ ------------------------------------------------------------
+
+ External definition of physical data and simple error covariance matrices
+ """
+__author__ = "Jean-Philippe ARGAUD"
+#
+# ==============================================================================
+#
+import numpy
+#
+def True_state():
+ """
+ Arbitrary values and names, as a tuple of two series of same length
+ """
+ return (numpy.array([1, 2, 3]), ['Para1', 'Para2', 'Para3'])
+#
+def Simple_Matrix( size, diagonal=None ):
+ """
+ Diagonal matrix, with either 1 or a given vector on the diagonal
+ """
+ if diagonal is not None:
+ S = numpy.diag( diagonal )
+ else:
+ S = numpy.matrix(numpy.identity(int(size)))
+ return S
+#
+# ==============================================================================
+if __name__ == "__main__":
+
+ print
+ print "AUTODIAGNOSTIC"
+ print "=============="
+
+ print
+ print "True_state = ", True_state()
+ print
+ print "B or R =\n",Simple_Matrix(3)
+ print
+ print "B or R =\n",Simple_Matrix(4, diagonal=numpy.arange(4,dtype=float))
+ print
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2011 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+# Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
+
+__doc__ = """
+ ADAO skeleton case, for wide script usage in case definition
+ ------------------------------------------------------------
+
+ External definition of the physical simulation operator, its tangent and
+ its adjoint form. In case adjoint is not known, a finite difference version
+ is given by default in this skeleton.
+ """
+__author__ = "Jean-Philippe ARGAUD"
+#
+# ==============================================================================
+#
+import os, numpy, time
+#
+def FunctionH( XX ):
+ """ Direct non-linear simulation operator """
+ #
+ # NEED TO BE COMPLETED
+ # NEED TO BE COMPLETED
+ # NEED TO BE COMPLETED
+ #
+ # --------------------------------------> EXAMPLE TO BE REMOVED
+ if type(XX) is type(numpy.matrix([])): # EXAMPLE TO BE REMOVED
+ HX = XX.A1.tolist() # EXAMPLE TO BE REMOVED
+ elif type(XX) is type(numpy.array([])): # EXAMPLE TO BE REMOVED
+ HX = numpy.matrix(XX).A1.tolist() # EXAMPLE TO BE REMOVED
+ else: # EXAMPLE TO BE REMOVED
+ HX = XX # EXAMPLE TO BE REMOVED
+ # --------------------------------------> EXAMPLE TO BE REMOVED
+ #
+ return numpy.array( HX )
+
+# ==============================================================================
+def TangentH( X, increment = 0.01, centeredDF = False ):
+ """
+ Calcul de l'opérateur tangent comme la Jacobienne par différences finies,
+ c'est-à-dire le gradient de H en X. On utilise des différences finies
+ directionnelles autour du point X.
+
+ Différences finies centrées :
+ 1/ Pour chaque composante i de X, on ajoute et on enlève la perturbation
+ dX[i] à la composante X[i], pour composer X_plus_dXi et X_moins_dXi, et
+ on calcule les réponses HX_plus_dXi = H( X_plus_dXi ) et HX_moins_dXi =
+ H( X_moins_dXi )
+ 2/ On effectue les différences (HX_plus_dXi-HX_moins_dXi) et on divise par
+ le pas 2*dXi
+ 3/ Chaque résultat, par composante, devient une colonne de la Jacobienne
+
+ Différences finies non centrées :
+ 1/ Pour chaque composante i de X, on ajoute la perturbation dX[i] à la
+ composante X[i] pour composer X_plus_dXi, et on calcule la réponse
+ HX_plus_dXi = H( X_plus_dXi )
+ 2/ On calcule la valeur centrale HX = H(X)
+ 3/ On effectue les différences (HX_plus_dXi-HX) et on divise par
+ le pas dXi
+ 4/ Chaque résultat, par composante, devient une colonne de la Jacobienne
+
+ """
+ print
+ print " == Calcul de la Jacobienne avec un incrément de %s*X"%increment
+ #
+ dX = increment * X.A1
+ #
+ if centeredDF:
+ #
+ # Boucle de calcul des colonnes de la Jacobienne
+ # ----------------------------------------------
+ Jacobienne = []
+ for i in range( len(dX) ):
+ X_plus_dXi = X.A1
+ X_plus_dXi[i] = X[i] + dX[i]
+ X_moins_dXi = X.A1
+ X_moins_dXi[i] = X[i] - dX[i]
+ #
+ HX_plus_dXi = FunctionH( X_plus_dXi )
+ HX_moins_dXi = FunctionH( X_moins_dXi )
+ #
+ HX_Diff = ( HX_plus_dXi - HX_moins_dXi ) / (2.*dX[i])
+ #
+ Jacobienne.append( HX_Diff )
+ #
+ else:
+ #
+ # Boucle de calcul des colonnes de la Jacobienne
+ # ----------------------------------------------
+ HX_plus_dX = []
+ for i in range( len(dX) ):
+ X_plus_dXi = X.A1
+ X_plus_dXi[i] = X[i] + dX[i]
+ #
+ HX_plus_dXi = FunctionH( X_plus_dXi )
+ #
+ HX_plus_dX.append( HX_plus_dXi )
+ #
+ # Calcul de la valeur centrale
+ # ----------------------------
+ HX = FunctionH( X )
+ #
+ # Calcul effectif de la Jacobienne par différences finies
+ # -------------------------------------------------------
+ Jacobienne = []
+ for i in range( len(dX) ):
+ Jacobienne.append( ( HX_plus_dX[i] - HX ) / dX[i] )
+ #
+ Jacobienne = numpy.matrix( Jacobienne )
+ print
+ print " == Fin du calcul de la Jacobienne"
+ #
+ return Jacobienne
+
+# ==============================================================================
+def AdjointH( (X, Y) ):
+ """
+ Calcul de l'adjoint à l'aide de la Jacobienne.
+ """
+ Jacobienne = TangentH( X, centeredDF = False )
+ #
+ # Calcul de la valeur de l'adjoint en X appliqué à Y
+ # --------------------------------------------------
+ Y = numpy.asmatrix(Y).flatten().T
+ HtY = numpy.dot(Jacobienne, Y)
+ #
+ return HtY.A1
+
+# ==============================================================================
+if __name__ == "__main__":
+
+ print
+ print "AUTODIAGNOSTIC"
+ print "=============="
+
+ from Physical_data_and_covariance_matrices import True_state
+ X0, noms = True_state()
+
+ FX = FunctionH( X0 )
+ print "FX =", FX
+ print
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2011 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+# Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
+
+__doc__ = """
+ ADAO skeleton case, for wide script usage in case definition
+ ------------------------------------------------------------
+
+ Script defining the AlgorithmParameters
+ """
+__author__ = "Jean-Philippe ARGAUD"
+#
+# ==============================================================================
+#
+# Creating the required ADAO variable
+# -----------------------------------
+AlgorithmParameters = {
+ "Minimizer" : "TNC", # Possible : "LBFGSB", "TNC", "CG", "BFGS"
+ "MaximumNumberOfSteps" : 15, # Number of iterative steps
+ "Bounds" : [
+ [ None, None ], # Bound on the first parameter
+ [ 0., 4. ], # Bound on the second parameter
+ [ 0., None ], # Bound on the third parameter
+ ],
+}
+#
+# ==============================================================================
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2011 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+# Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
+
+__doc__ = """
+ ADAO skeleton case, for wide script usage in case definition
+ ------------------------------------------------------------
+
+ Script defining the BackgroundError
+ """
+__author__ = "Jean-Philippe ARGAUD"
+#
+# ==============================================================================
+#
+import sys, os ; sys.path.insert(0,os.path.dirname(os.path.abspath(__file__)))
+from Physical_data_and_covariance_matrices import True_state, Simple_Matrix
+#
+xt, names = True_state()
+#
+B = 0.1 * Simple_Matrix( size = len(xt) )
+#
+# Creating the required ADAO variable
+# -----------------------------------
+BackgroundError = B
+#
+# ==============================================================================
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2011 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+# Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
+
+__doc__ = """
+ ADAO skeleton case, for wide script usage in case definition
+ ------------------------------------------------------------
+
+ Script defining the Background
+ """
+__author__ = "Jean-Philippe ARGAUD"
+#
+# ==============================================================================
+#
+import sys, os ; sys.path.insert(0,os.path.dirname(os.path.abspath(__file__)))
+from Physical_data_and_covariance_matrices import True_state
+import numpy
+#
+xt, names = True_state()
+#
+Standard_deviation = 0.2*xt # 20% for each variable
+#
+xb = xt + abs(numpy.random.normal(0.,Standard_deviation,size=(len(xt),)))
+#
+# Creating the required ADAO variable
+# -----------------------------------
+Background = list(xb)
+#
+# ==============================================================================
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2011 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+# Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
+
+__doc__ = """
+ ADAO skeleton case, for wide script usage in case definition
+ ------------------------------------------------------------
+
+ Script defining the ObservationError
+ """
+__author__ = "Jean-Philippe ARGAUD"
+#
+# ==============================================================================
+#
+import sys, os ; sys.path.insert(0,os.path.dirname(os.path.abspath(__file__)))
+from Physical_data_and_covariance_matrices import True_state, Simple_Matrix
+from Physical_simulation_functions import FunctionH
+#
+xt, names = True_state()
+#
+yo = FunctionH( xt )
+#
+R = 0.0001 * Simple_Matrix( size = len(yo) )
+#
+# Creating the required ADAO variable
+# -----------------------------------
+ObservationError = R
+#
+# ==============================================================================
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2011 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+# Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
+
+__doc__ = """
+ ADAO skeleton case, for wide script usage in case definition
+ ------------------------------------------------------------
+
+ Script defining the ObservationOperator
+ """
+__author__ = "Jean-Philippe ARGAUD"
+#
+# ==============================================================================
+#
+import sys, os ; sys.path.insert(0,os.path.dirname(os.path.abspath(__file__)))
+import Physical_simulation_functions
+import numpy, logging
+#
+# -----------------------------------------------------------------------
+# SALOME input data and parameters: all information are the required input
+# variable "computation", containing for example:
+# {'inputValues': [[[[0.0, 0.0, 0.0]]]],
+# 'inputVarList': ['adao_default'],
+# 'outputVarList': ['adao_default'],
+# 'specificParameters': [{'name': 'method', 'value': 'Direct'}]}
+# -----------------------------------------------------------------------
+#
+# Recovering the type of computation: "Direct", "Tangent" or "Adjoint"
+# --------------------------------------------------------------------
+method = ""
+for param in computation["specificParameters"]:
+ if param["name"] == "method":
+ method = param["value"]
+logging.info("ComputationFunctionNode: Found method is \'%s\'"%method)
+#
+# Loading the H operator functions from external definitions
+# ----------------------------------------------------------
+logging.info("ComputationFunctionNode: Loading operator functions")
+FunctionH = Physical_simulation_functions.FunctionH
+AdjointH = Physical_simulation_functions.AdjointH
+#
+# Executing the possible computations
+# -----------------------------------
+if method == "Direct":
+ logging.info("ComputationFunctionNode: Direct computation")
+ Xcurrent = computation["inputValues"][0][0][0]
+ data = FunctionH(numpy.matrix( Xcurrent ).T)
+#
+if method == "Tangent":
+ logging.info("ComputationFunctionNode: Tangent computation")
+ Xcurrent = computation["inputValues"][0][0][0]
+ data = FunctionH(numpy.matrix( Xcurrent ).T)
+#
+if method == "Adjoint":
+ logging.info("ComputationFunctionNode: Adjoint computation")
+ Xcurrent = computation["inputValues"][0][0][0]
+ Ycurrent = computation["inputValues"][0][0][1]
+ data = AdjointH((numpy.matrix( Xcurrent ).T, numpy.matrix( Ycurrent ).T))
+#
+# Formatting the output
+# ---------------------
+logging.info("ComputationFunctionNode: Formatting the output")
+it = data.flat
+outputValues = [[[[]]]]
+for val in it:
+ outputValues[0][0][0].append(val)
+#
+# Creating the required ADAO variable
+# -----------------------------------
+result = {}
+result["outputValues"] = outputValues
+result["specificOutputInfos"] = []
+result["returnCode"] = 0
+result["errorMessage"] = ""
+#
+# ==============================================================================
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2011 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+# Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
+
+__doc__ = """
+ ADAO skeleton case, for wide script usage in case definition
+ ------------------------------------------------------------
+
+ Script defining the Observation
+ """
+__author__ = "Jean-Philippe ARGAUD"
+#
+# ==============================================================================
+#
+import sys, os ; sys.path.insert(0,os.path.dirname(os.path.abspath(__file__)))
+from Physical_data_and_covariance_matrices import True_state
+from Physical_simulation_functions import FunctionH
+#
+xt, noms = True_state()
+#
+yo = FunctionH( xt )
+#
+# Creating the required ADAO variable
+# -----------------------------------
+Observation = list(yo)
+#
+# ==============================================================================
--- /dev/null
+#-*-coding:iso-8859-1-*-
+#
+# Copyright (C) 2008-2011 EDF R&D
+#
+# This library is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License as published by the Free Software Foundation; either
+# version 2.1 of the License.
+#
+# This library is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with this library; if not, write to the Free Software
+# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+#
+# See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+#
+# Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
+
+__doc__ = """
+ ADAO skeleton case, for wide script usage in case definition
+ ------------------------------------------------------------
+
+ Script defining the UserPostAnalysis
+ """
+__author__ = "Jean-Philippe ARGAUD"
+#
+# ==============================================================================
+#
+import sys, os ; sys.path.insert(0,os.path.dirname(os.path.abspath(__file__)))
+from Physical_data_and_covariance_matrices import True_state
+import numpy
+#
+xt, names = True_state()
+xa = ADD.get("Analysis").valueserie(-1)
+x_series = ADD.get("CurrentState").valueserie()
+J = ADD.get("CostFunctionJ").valueserie()
+#
+# Verifying the results by printing
+# ---------------------------------
+print
+print "xt = %s"%xt
+print "xa = %s"%numpy.array(xa)
+print
+for i in range( len(x_series) ):
+ print "Step %2i : J = %.5e et X = %s"%(i, J[i], x_series[i])
+print
+#
+# ==============================================================================
