]> SALOME platform Git repositories - modules/adao.git/blob - src/daComposant/daNumerics/pst4mod/modelica_calibration/case.py
Salome HOME
Documentation and models update, pst4mod
[modules/adao.git] / src / daComposant / daNumerics / pst4mod / modelica_calibration / case.py
1 # -*- coding: utf-8 -*-
2 #
3 # Copyright (C) 2016-2024 EDF R&D
4 #
5 # This library is free software; you can redistribute it and/or
6 # modify it under the terms of the GNU Lesser General Public
7 # License as published by the Free Software Foundation; either
8 # version 2.1 of the License.
9 #
10 # This library is distributed in the hope that it will be useful,
11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
13 # Lesser General Public License for more details.
14 #
15 # You should have received a copy of the GNU Lesser General Public
16 # License along with this library; if not, write to the Free Software
17 # Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
18 #
19 # Author: Jean-Philippe Argaud, jean-philippe.argaud@edf.fr, EDF R&D
20 # Author: Luis Corona Mesa-Molles, luis.corona-mesa-molles@edf.fr, EDF R&D
21
22 __all__ = ["Calibration", "name", "version", "year"]
23
24 # ==============================================================================
25
26 name     = "Modelica and Dymola Calibration Tools"
27 version  = "1.0.9.7.0" # "x.x"+"adao version"
28 year     = "2021"
29
30 # ==============================================================================
31 # Default configuration
32 # ---------------------
33 import os, sys, io, shutil, numpy, scipy, tempfile, time, logging, pandas, subprocess, copy, csv
34 from datetime import timedelta
35 from datetime import datetime
36
37 try:
38     import adao
39     from adao import adaoBuilder
40     from daCore.Interfaces import ImportFromFile, ImportScalarLinesFromFile
41 except ImportError:
42     raise ImportError("ADAO module not found, please install it first.")
43
44 try:
45     from buildingspy.io.outputfile import Reader
46 except ImportError:
47     raise ImportError("buildingspy module not found, please install it first.")
48
49 try:
50     from pst4mod.modelica_libraries.automatic_simulation import Automatic_Simulation
51     from pst4mod.modelica_libraries.around_simulation import Around_Simulation
52     from pst4mod.modelica_libraries import write_in_dsin
53     from pst4mod.modelica_libraries.change_working_point import Dict_Var_To_Fix
54     from pst4mod.modelica_libraries.change_working_point import Working_Point_Modification
55 except ImportError:
56     raise ImportError("modelica_libraries module not found, please install it first.")
57
58 try:
59     from fmpy import simulate_fmu
60     from fmpy.fmi1 import FMICallException
61 # except ImportError:
62 #     raise ImportError("fmpy library not found, please install it first")
63 except:
64     pass
65
66 _configuration = {
67     "DirectoryModels"   : "Models",
68     "DirectoryMeasures" : "Measures",
69     "DirectoryMethods"  : "Methods",
70     "DirectoryResults"  : "Results",
71     "Launcher"          : "configuration.py",
72     }
73 __paths = []
74 __paths.append( os.path.abspath(os.path.join(os.path.dirname(__file__), '..','..','buildingspy-1.5.0')) )
75 __paths.append( os.path.abspath(os.path.join(os.path.dirname(__file__), '..','..','Packages')) )
76 for __path in __paths:
77     if os.path.exists(__path): sys.path.append(__path)
78
79 # ==============================================================================
80 class Calibration(object):
81     "ADAO case based parameters inverse estimation tools"
82     def __init__(self, Name = "Calibration case", SaveStdoutIn = None, Verbose = False):
83         """
84         Name           : name as a string
85         SaveAllStdoutIn: filename to be used for stdout stream
86         Verbose        : verbose output
87         """
88         self.__name     = str(Name)
89         self.__measures = {}
90         self.__model    = {}
91         self.__link     = {}
92         self.__method   = {}
93         self.__result   = {}
94         self._setConfigurationDefaults()
95         self.__verbose  = bool(Verbose)
96         self.__stdoutid = sys.stdout
97         if SaveStdoutIn is not None:
98             sys.stdout = open(SaveStdoutIn,"w")
99         if self.__verbose:
100             __msg = "[VERBOSE] %s"%self.__name
101             print("")
102             print("  %s"%__msg)
103             print("  %s"%("-"*len(__msg),))
104             print("")
105             VersionsLogicielles()
106
107     def _setConfigurationDefaults(self, Configuration = None):
108         """
109         Impose case directory and files structure configuration defaults based
110         on argument dictionary
111         """
112         if Configuration is None or type(Configuration) is not dict:
113             Configuration = _configuration
114         for __k,__v in Configuration.items():
115             setattr(self, __k, __v)
116
117     def _setModelTmpDir(self, __model_name="dsin.txt", __model_format="DYMOSIM", __model_dest = "dsin.txt", __before_tmp = None, __suffix=None):
118         if __model_name is None:
119             raise ValueError('Model file or directory has to be set and can not be None')
120         elif os.path.isfile(__model_name):
121             __model_nam = os.path.basename(__model_name)
122             __model_dir = os.path.abspath(os.path.dirname(__model_name))
123             __model_dst = __model_dest
124         elif os.path.isdir(__model_name):
125             __model_nam = "dsin.txt"
126             __model_dir = os.path.abspath(__model_name)
127             __model_dst = __model_dest
128         else:
129             raise ValueError('Model file or directory not found using %s'%str(__model_name))
130         #
131         if __before_tmp is not None: # Mettre "../.." si nécessaire
132             __mtmp = os.path.join(__model_dir, __before_tmp, "tmp")
133         else:
134             __mtmp = os.path.join(__model_dir, "tmp")
135         if not os.path.exists(__mtmp):
136             os.mkdir(__mtmp)
137         __prefix = time.strftime('%Y%m%d_%Hh%Mm%Ss_tmp_',time.localtime())
138         if __suffix is not None:
139             __prefix = __prefix + __suffix + "_"
140         __ltmp = tempfile.mkdtemp( prefix=__prefix, dir=__mtmp )
141         #
142         for sim, pmf, dst in (
143             (__model_nam,__model_format,__model_dst),
144             ("dymosim.exe","DYMOSIM","dymosim.exe"),
145             ("pythonsim.exe","PYSIM","pythonsim.exe"),
146             ):
147             # Recherche un fichier de données ou un simulateur dans cet ordre
148
149             if os.path.isfile(os.path.join(__model_dir, sim)) and (pmf.upper() in [__model_format, "GUESS"]):
150                 shutil.copy(
151                     os.path.join(__model_dir, sim),
152                     os.path.join(__ltmp, dst)
153                     )
154
155                 # _secure_copy_textfile(
156                 #         os.path.join(__model_dir, sim),
157                 #         os.path.join(__ltmp, dst)
158                 #         )
159
160                 break
161         return __ltmp
162
163     def _setModelTmpDir_simple(self, __model_name="dsin.txt", __model_format="DYMOSIM", __model_dest = "dsin.txt"):
164         if __model_name is None:
165             raise ValueError('Model file or directory has to be set and can not be None')
166         elif os.path.isfile(__model_name):
167             if __model_format == "OPENMODELICA": #same structure as in  the directory case
168                 __model_dir = os.path.abspath(os.path.dirname(__model_name))
169
170                 __model_nam_1 = os.path.basename(__model_name)
171                 __model_nam_2 = __model_nam_1[:-9]+'.exe'
172                 __model_nam = [__model_nam_1, __model_nam_2] #the first found model is kept
173
174                 __model_nam_3 = __model_nam_1[:-9]+'_info.json' #check if this file exists as well
175                 if os.path.exists(os.path.join(__model_dir,__model_nam_3)):
176                     __model_nam.append(__model_nam_3)                      #get the three files necessar for simulation
177                 __model_dst = __model_nam #the same file name is kept
178
179             else: #cas classique
180                 __model_nam = os.path.basename(__model_name)
181                 __model_dir = os.path.abspath(os.path.dirname(__model_name))
182                 __model_dst = os.path.basename(__model_dest)
183         elif os.path.isdir(__model_name):
184             if __model_format == "DYMOSIM":
185                 __model_nam = "dsin.txt"
186                 __model_dir = os.path.abspath(__model_name)
187                 __model_dst = os.path.basename(__model_dest)
188
189             elif __model_format == "FMI" or __model_format == "FMU":
190                 __model_dir = os.path.abspath(__model_name)
191                 __model_nam = [files for files in os.listdir(__model_dir) if files[-4:] =='.fmu'][0] #the first found model is kept
192                 __model_dst = __model_nam #the same file name is kept
193
194             elif __model_format == "OPENMODELICA":
195                 __model_dir = os.path.abspath(__model_name)
196                 __model_nam_1 = [files for files in os.listdir(__model_dir) if files[-4:] =='.xml'][0] #the first found model is kept, .xml extension is considered
197                 __model_nam_2 = __model_nam_1[:-9]+'.exe'
198                 __model_nam = [__model_nam_1, __model_nam_2] #the first found model is kept
199
200                 __model_nam_3 = __model_nam_1[:-9]+'_info.json' #check if this file exists as well
201                 if os.path.exists(os.path.join(__model_dir,__model_nam_3)):
202                     __model_nam.append(__model_nam_3)                      #get the three files necessar for simulation
203                 __model_dst = __model_nam #the same file name is kept
204
205             else :
206                 raise ValueError('If the model is not provided in DYMOSIM, OpenModelica or FMI format, the name of the model must refer to a file name (in a relative way with respect to the configuration file)')
207         else:
208             raise ValueError('Model file or directory not found using %s'%str(__model_name))
209         #
210         __mtmp = os.path.join(__model_dir, "tmp")
211         if not os.path.exists(__mtmp):
212             os.mkdir(__mtmp)
213         __ltmp = __mtmp
214         #
215         if type(__model_nam) == list: #it means that it is an OpenModelica model
216             for sim_element in __model_nam:
217                 shutil.copy(
218                     os.path.join(__model_dir, sim_element),
219                     os.path.join(__ltmp, sim_element) #the same name is kept for the files
220                     )
221         else:
222             for sim, pmf, dst in (
223                 (__model_nam,__model_format,__model_dst),
224                 ("dsin.txt","DYMOSIM","dsin.txt"),
225                 ("pythonsim.exe","PYSIM","pythonsim.exe"),
226                 ):
227                 # Recherche un fichier de données ou un simulateur dans cet ordre
228                 if os.path.isfile(os.path.join(__model_dir, sim)) and (pmf.upper() in [__model_format, "GUESS"]):
229
230                     shutil.copy(
231                         os.path.join(__model_dir, sim),
232                         os.path.join(__ltmp, dst)
233                         )
234                     # _secure_copy_textfile(
235                     #         os.path.join(__model_dir, sim),
236                     #         os.path.join(__ltmp, dst)
237                     #         )
238
239                     break
240         return __ltmp
241
242     def _get_Name_ModelTmpDir_simple(self, __model_name="dsin.txt"):
243         if __model_name is None:
244             raise ValueError('Model file or directory has to be set and can not be None')
245         elif os.path.isfile(__model_name):
246             __model_dir = os.path.abspath(os.path.dirname(__model_name))
247         elif os.path.isdir(__model_name):
248             __model_dir = os.path.abspath(__model_name)
249         else:
250             raise ValueError('Model file or directory not found using %s'%str(__model_name))
251         #
252         __mtmp = os.path.join(__model_dir, "tmp")
253         __ltmp = __mtmp
254         return __ltmp
255
256     def _setModelTmpDir_REF_CALCULATION(self, __model_name="dsin.txt", __model_format="DYMOSIM", __suffix=None, __model_dest = "dsin.txt"):
257         if __model_name is None:
258             raise ValueError('Model file or directory has to be set and can not be None')
259         elif os.path.isfile(__model_name):
260             __model_dir = os.path.abspath(os.path.dirname(__model_name))
261             __model_nam = os.path.basename(__model_name)
262         elif os.path.isdir(__model_name):
263             __model_dir = os.path.abspath(__model_name)
264             __model_nam = "dsin.txt"
265         else:
266             raise ValueError('Model file or directory not found using %s'%str(__model_name))
267         if __suffix is None:
268             raise ValueError('Observation files must be named')
269         #
270         __mtmp = os.path.join(__model_dir, "tmp")
271         if not os.path.exists(__mtmp):
272             os.mkdir(__mtmp)
273         for sim, pmf in ((__model_nam,__model_format), ("dymosim.exe","DYMOSIM"), ("pythonsim.exe","PYSIM")):
274             # Recherche un simulateur dans cet ordre
275             if os.path.isfile(os.path.join(__model_dir, sim)) and (pmf.upper() in [__model_format, "GUESS"]):
276                 __pref = time.strftime('%Y%m%d_%Hh%Mm%Ss_REF_CALCULATION_',time.localtime())
277                 if __suffix is not None:
278                     __pref = __pref+__suffix+"_"
279                 __ltmp = tempfile.mkdtemp( prefix=__pref, dir=__mtmp )
280                 shutil.copy(
281                     os.path.join(__model_dir, sim),
282                     os.path.join(__ltmp, __model_dest)
283                     )
284                 # _secure_copy_textfile(
285                 #         os.path.join(__model_dir, sim),
286                 #         os.path.join(__ltmp, __model_dest)
287                 #         )
288
289                 __ok = True
290                 break
291             else:
292                 __ok = False
293         if not __ok:
294             raise ValueError("Model simulator not found using \"%s\" path and \"%s\" format"%(__model_name,__model_format))
295         return __ltmp
296
297     def describeMeasures(self, SourceName, Variables, Format="Guess"):
298         "To store description of measures"
299         self.__measures["SourceName"] = str(SourceName) # File or data name
300         self.__measures["Variables"]  = tuple(Variables) # Name of variables
301         self.__measures["Format"]     = str(Format) # File or data format
302         if self.__verbose: print("  [VERBOSE] Measures described")
303         return self.__measures
304
305     def describeMultiMeasures(self, SourceNames, Variables, Format="Guess"):
306         "To store description of measures"
307         self.__measures["SourceName"] = tuple(SourceNames) # List of file or data name
308         self.__measures["Variables"]  = tuple(Variables) # Name of variables
309         self.__measures["Format"]     = str(Format) # File or data format
310         if self.__verbose: print("  [VERBOSE] Measures described")
311         return self.__measures
312
313     def describeModel(self, SourceName, VariablesToCalibrate, OutputVariables, Format="Guess"):
314         "To store description of model"
315         self.__model["SourceName"]           = str(SourceName) # Model command file
316         self.__model["VariablesToCalibrate"] = tuple(VariablesToCalibrate) # Name of variables
317         self.__model["OutputVariables"]      = tuple(OutputVariables) # Name of variables
318         self.__model["Format"]               = str(Format) # File or model format
319         if os.path.isfile(SourceName):
320             self.__model["Verbose_printing_name"] = os.path.basename(SourceName)
321         elif os.path.isdir(SourceName):
322             if self.__model["Format"] == "DYMOSIM":
323                 self.__model["Verbose_printing_name"] = "dsin.txt + dymosim.exe in " + os.path.basename(SourceName)
324             elif self.__model["Format"] == "FMI" or self.__model["Format"] == "FMU":
325                 self.__model["Verbose_printing_name"] = [files for files in os.listdir(os.path.abspath(SourceName)) if files[-4:] =='.fmu'][0]
326             elif self.__model["Format"] == "OPENMODELICA":
327                 xml_file_name = [files for files in os.listdir(os.path.abspath(SourceName)) if files[-4:] =='.xml'][0]
328                 exe_file_name = xml_file_name[:-9] + '.exe'
329                 json_file_name = xml_file_name[:-9] + '_info.json'
330                 if os.path.exists(os.path.join(os.path.abspath(SourceName),json_file_name)):
331                     self.__model["Verbose_printing_name"] = xml_file_name + ' ' + exe_file_name + ' ' + json_file_name
332                 else:
333                     self.__model["Verbose_printing_name"] = xml_file_name + ' ' + exe_file_name + ' '
334             else:
335                 print("  [VERBOSE] Model described: ", "/!\\ Unknown ModelFormat /!\\ ")
336                 exit()
337         if self.__verbose: print("  [VERBOSE] Model described: ", self.__model["Verbose_printing_name"] )
338         return self.__model
339
340     def describeLink(self, MeasuresNames, ModelName, LinkName, LinkVariable, Format="Guess"):
341         "To store description of link between measures and model input data"
342         self.__link["SourceName"]    = str(LinkName) # File or link name
343         self.__link["MeasuresNames"] = tuple(MeasuresNames) # List of file or name
344         self.__link["ModelName"]     = tuple(ModelName) # File or model name
345         self.__link["Variable"]      = str(LinkVariable) # Name of variable names column
346         self.__link["Format"]        = str(Format) # File or model format
347         if self.__verbose: print("  [VERBOSE] Link between model and multiple measures described")
348         return self.__link
349
350     def describeMethod(self, SourceName, MethodFormat="PY", BackgroundName=None, BackgroundFormat="Guess"):
351         "To store description of calibration method"
352         self.__method["SourceName"]       = str(SourceName) # File or method name
353         self.__method["Format"]           = str(MethodFormat) # File or model format
354         self.__method["BackgroundName"]   = str(BackgroundName) # File or background name
355         self.__method["BackgroundFormat"] = str(BackgroundFormat) # File or model format
356         if self.__verbose: print("  [VERBOSE] Optimization method settings described")
357         return self.__method
358
359     def describeResult(self, ResultName, Level=None, Format="TXT"):
360         "To store description of results"
361         self.__result["SourceName"] = str(ResultName) # File or method name
362         self.__result["Level"]      = str(Level) # Level = "Final", "IntermediaryFinal"
363         self.__result["Format"]     = str(Format) # File or model format
364         if self.__verbose: print("  [VERBOSE] Results settings described")
365         return self.__result
366
367     def calibrate(self,
368             ModelName            = None, # Model command file
369             ModelFormat          = "Guess",
370             DataName             = None, # Measures
371             DataFormat           = "Guess",
372             BackgroundName       = None, # Background
373             BackgroundFormat     = "Guess",
374             MethodName           = None, # Assimilation
375             MethodFormat         = "Guess",
376             VariablesToCalibrate = None,
377             OutputVariables      = None,
378             ResultName           = "results.txt", # Results
379             ResultFormat         = "Guess",
380             ResultLevel          = None,
381             Verbose              = False,
382             ):
383         """
384         General method to set and realize a calibration (optimal estimation
385         task) of model parameters using measures data.
386         """
387         if Verbose:
388             self.__verbose = True
389             print("  [VERBOSE] Verbose mode activated")
390         else:
391             self.__verbose = False
392         if DataName is not None:
393             # On force l'utilisateur à nommer dans son fichier de mesures
394             # les variables avec le même nom que des sorties dans le modèle
395             self.describeMeasures(DataName, OutputVariables, DataFormat)
396         if ModelName is not None:
397             self.describeModel(ModelName, VariablesToCalibrate, OutputVariables, ModelFormat)
398         if MethodName is not None:
399             self.describeMethod(MethodName, MethodFormat, BackgroundName, BackgroundFormat)
400         if ResultName is not None:
401             self.describeResult(ResultName, ResultLevel, ResultFormat)
402         #
403         ONames, Observations = _readMeasures(
404             self.__measures["SourceName"],
405             self.__measures["Variables"],
406             self.__measures["Format"],
407             )
408         Algo, Params, CovB, CovR = _readMethod(
409             self.__method["SourceName"],
410             self.__method["Format"],
411             )
412         __bgfile = None
413         if self.__method["BackgroundFormat"] == "DSIN":  __bgfile = self.__model["SourceName"]
414         if self.__method["BackgroundFormat"] == "ADAO":  __bgfile = self.__method["SourceName"]
415         if self.__method["BackgroundFormat"] == "USER":  __bgfile = self.__method["BackgroundName"]
416         if self.__method["BackgroundFormat"] == "Guess":
417             if self.__method["BackgroundName"] is not None: __bgfile = self.__method["BackgroundName"]
418             if self.__method["SourceName"] is not None:     __bgfile = self.__method["SourceName"]
419             if self.__model["SourceName"] is not None:      __bgfile = self.__model["SourceName"]
420         BNames, Background, Bounds = _readBackground(
421             __bgfile,
422             self.__model["VariablesToCalibrate"],
423             self.__method["BackgroundFormat"],
424             )
425         if "Bounds" not in Params: # On force la priorité aux bornes utilisateur
426             Params["Bounds"] = Bounds
427         if self.__verbose:
428             print("  [VERBOSE] Measures read")
429             print("  [VERBOSE] Optimization information read")
430             if "MaximumNumberOfSteps" in Params:
431                 print("  [VERBOSE] Maximum possible number of iteration:",Params['MaximumNumberOfSteps'])
432             print("  [VERBOSE] Background read:",Background)
433             print("  [VERBOSE] Bounds read:",Params['Bounds'])
434         #
435         if "DifferentialIncrement" not in Params:
436             Params["DifferentialIncrement"] = 0.001
437         if "StoreSupplementaryCalculations" not in Params:
438             Params["StoreSupplementaryCalculations"] = ["SimulatedObservationAtOptimum",]
439         if "StoreSupplementaryCalculations" in Params and \
440             "SimulatedObservationAtOptimum" not in Params["StoreSupplementaryCalculations"]:
441             Params["StoreSupplementaryCalculations"].append("SimulatedObservationAtOptimum")
442         if self.__verbose and "StoreSupplementaryCalculations" in Params and \
443             "CurrentOptimum" not in Params["StoreSupplementaryCalculations"]:
444             Params["StoreSupplementaryCalculations"].append("CurrentOptimum")
445         if self.__verbose and "StoreSupplementaryCalculations" in Params and \
446             "CostFunctionJAtCurrentOptimum" not in Params["StoreSupplementaryCalculations"]:
447             Params["StoreSupplementaryCalculations"].append("CostFunctionJAtCurrentOptimum")
448         #
449         def exedymosim( x_values_matrix ):
450             "Appel du modèle et restitution des résultats"
451             from buildingspy.io.outputfile import Reader
452             from pst4mod.modelica_libraries.automatic_simulation import Automatic_Simulation
453             #
454             # x_values = list(numpy.ravel(x_values_matrix) * Background)
455             x_values = list(numpy.ravel(x_values_matrix))
456             dict_inputs=dict(zip(VariablesToCalibrate,x_values))
457             # if Verbose: print("  Simulation for %s"%numpy.ravel(x_values))
458             #
459             simudir = self._setModelTmpDir(ModelName, ModelFormat)
460             auto_simul = Automatic_Simulation(
461                 simu_dir=simudir,
462                 dict_inputs=dict_inputs,
463                 logfile=True,
464                 timeout=60,
465                 without_modelicares=True) #linux=true is removed
466             auto_simul.single_simulation()
467             reader = Reader(os.path.join(simudir,'dsres.mat'),'dymola')
468             y_values = [reader.values(y_name)[1][-1] for y_name in OutputVariables]
469             y_values_matrix = numpy.asarray(y_values)
470             if not Verbose:
471                 shutil.rmtree(simudir, ignore_errors=True)
472             return y_values_matrix
473         #
474         def exepython( x_values_matrix ):
475             "Appel du modèle et restitution des résultats"
476             # x_values = list(numpy.ravel(x_values_matrix))
477             # dict_inputs=dict(zip(VariablesToCalibrate,x_values))
478             #
479             simudir = self._setModelTmpDir(ModelName, ModelFormat)
480             auto_simul = Python_Simulation(
481                 simu_dir=simudir,
482                 val_inputs=x_values_matrix,
483                 logfile=True,
484                 timeout=60,
485                 verbose=self.__verbose)
486             y_values_matrix = auto_simul.single_simulation()
487             if not Verbose:
488                 shutil.rmtree(simudir, ignore_errors=True)
489             return y_values_matrix
490         #
491         if self.__model["Format"].upper() in ["DYMOSIM", "GUESS"]:
492             simulateur = exedymosim
493         elif self.__model["Format"].upper() == "PYSIM":
494             simulateur = exepython
495         else:
496             raise ValueError("Unknown model format \"%s\""%self.__model["Format"].upper())
497         #
498         __adaocase = adaoBuilder.New()
499         __adaocase.set( 'AlgorithmParameters',  Algorithm=Algo, Parameters=Params )
500         __adaocase.set( 'Background',           Vector=Background)
501         __adaocase.set( 'Observation',          Vector=Observations )
502         if type(CovB) is float:
503             __adaocase.set( 'BackgroundError',  ScalarSparseMatrix=CovB )
504         else:
505             __adaocase.set( 'BackgroundError',  Matrix=CovB )
506         if type(CovR) is float:
507             __adaocase.set( 'ObservationError', ScalarSparseMatrix=CovR )
508         else:
509             __adaocase.set( 'ObservationError', Matrix=CovR )
510         __adaocase.set( 'ObservationOperator',  OneFunction=simulateur, Parameters = {"DifferentialIncrement":Params["DifferentialIncrement"]} )
511         if self.__verbose:
512             __adaocase.set( 'Observer', Variable="CostFunctionJAtCurrentOptimum", String="print('\\n  -----------------------------------\\n  [VERBOSE] Current iteration: %i'%(len(var),))" )
513             __adaocase.set( 'Observer', Variable="CostFunctionJAtCurrentOptimum", Template="ValuePrinter", Info="\n  [VERBOSE] Current best cost value:" )
514             __adaocase.set( 'Observer', Variable="CurrentOptimum",                Template="ValuePrinter", Info="\n  [VERBOSE] Current optimal state..:" )
515         __adaocase.execute()
516         #
517         __resultats = dict(
518             InitialParameters = numpy.asarray(Background),
519             OptimalParameters = __adaocase.get("Analysis")[-1], # * Background,
520             OptimalSimulation = __adaocase.get("SimulatedObservationAtOptimum")[-1],
521             Measures          = Observations,
522             NamesOfParameters = BNames,
523             NamesOfMeasures   = ONames,
524             )
525         if self.__verbose: print("")
526         del __adaocase
527         #
528         _saveResults(
529             __resultats,
530             self.__result["SourceName"],
531             self.__result["Level"],
532             self.__result["Format"],
533             self.__verbose,
534             )
535         #
536         return __resultats
537
538
539
540     def calibrateMultiObs(self,
541             ModelName            = None, # Model command file
542             ModelFormat          = "Guess",
543             DataNames            = None, # Multiple measures
544             DataFormat           = "Guess",
545             LinkName             = None, # CL with data names
546             LinkVariable         = "Variable", # Name of variable names column
547             LinkFormat           = "Guess",
548             BackgroundName       = None, # Background
549             BackgroundFormat     = "Guess",
550             MethodName           = None, # Assimilation
551             MethodFormat         = "Guess",
552             VariablesToCalibrate = None,
553             OutputVariables      = None,
554             ResultName           = "results.txt", # Results
555             ResultFormat         = "Guess",
556             ResultLevel          = None,
557             Verbose              = False,
558             IntermediateInformation = False,
559             ComplexModel         = False,
560             FMUInput             = None,
561             NeedInitVariables         = False,
562             KeepCalculationFolders = False,
563             Linux = False,
564             CreateOptimaldsin = False,
565             InitialSimulation = False,
566             VerifyFunction = False,
567             GlobalSensitivityCheck = False,
568             ParamSensitivityCheck = False,
569             CreateCSVoutputResults = False,
570             AdaptObservationError = False,
571             ResultsSummary = False,
572             AdvancedDebugModel = False,
573             TimeoutModelExecution = 300,
574             ModelStabilityCheck = False,
575             ModelStabilityCheckingPoints = 10
576             ):
577         """
578         General method to set and realize a calibration (optimal estimation
579         task) of model parameters using multiple measures data, with an
580         explicit link between measures and data for the model.
581         """
582         if Verbose:
583             self.__verbose = True
584             print("  [VERBOSE] Verbose mode activated")
585         if IntermediateInformation:
586             if self.__verbose:
587                 self.__IntermediateInformation = True
588                 print("  [VERBOSE] IntermediateInformation provided")
589             else:
590                 self.__IntermediateInformation = False
591                 print("  [VERBOSE] IntermediateInformation requires Verbose mode activated")
592         else:
593             self.__IntermediateInformation = False
594
595         if TimeoutModelExecution < 0:
596             raise ValueError("TimeoutModelExecution must be positive")
597
598         if DataNames is not None:
599             # On force l'utilisateur à nommer dans son fichier de mesures
600             # les variables avec le même nom que des sorties dans le modèle
601             self.describeMultiMeasures(DataNames, OutputVariables, DataFormat)
602         #
603         if self.__verbose:
604             print("  [VERBOSE] Timeoout for model execution is:", TimeoutModelExecution, "seconds")
605         #
606         if ModelFormat is not None: #included sot that the model format is allways in upper format
607             ModelFormat = ModelFormat.upper()
608             if ModelFormat == "DYMOLA": #to handle situations in which name Dymola is indicated
609                 ModelFormat = "DYMOSIM"
610         if ModelName is not None:
611             self.describeModel(ModelName, VariablesToCalibrate, OutputVariables, ModelFormat)
612         if LinkName is not None:
613             # On force l'utilisateur à utiliser les fichiers de mesures requis
614             # et non pas tous ceux présents dans le Link
615             self.describeLink(DataNames, ModelName, LinkName, LinkVariable, LinkFormat)
616         if MethodName is not None:
617             self.describeMethod(MethodName, MethodFormat, BackgroundName, BackgroundFormat)
618         if ResultName is not None:
619             self.describeResult(ResultName, ResultLevel, ResultFormat)
620         #
621         #handling of the new option for complex models while keeping the previous key word in the code
622         if ComplexModel:
623             print(" The option ComplexModel  is deprecated and has to be replaced by NeedInitVariables. Please update your code.")
624             NeedInitVariables = ComplexModel
625         ComplexModel = NeedInitVariables
626        #
627
628         ONames, Observations = _readMultiMeasures(
629             self.__measures["SourceName"],
630             self.__measures["Variables"],
631             self.__measures["Format"],
632             )
633
634
635         Algo, Params, CovB, CovR = _readMethod(
636             self.__method["SourceName"],
637             self.__method["Format"],
638             )
639         __bgfile = None
640         if self.__method["BackgroundFormat"] == "DSIN":  __bgfile = self.__model["SourceName"]
641         if self.__method["BackgroundFormat"] == "ADAO":  __bgfile = self.__method["SourceName"]
642         if self.__method["BackgroundFormat"] == "USER":  __bgfile = self.__method["BackgroundName"]
643         if self.__method["BackgroundFormat"] == "Guess":
644             if self.__method["BackgroundName"] is not None: __bgfile = self.__method["BackgroundName"]
645             if self.__method["SourceName"] is not None:     __bgfile = self.__method["SourceName"]
646             if self.__model["SourceName"] is not None:      __bgfile = self.__model["SourceName"]
647         BNames, Background, Bounds = _readBackground(
648             __bgfile,
649             self.__model["VariablesToCalibrate"],
650             self.__method["BackgroundFormat"],
651             )
652         if "Bounds" not in Params: # On force la priorité aux bornes utilisateur
653             Params["Bounds"] = Bounds
654         LNames, LColumns, LVariablesToChange = _readLink(
655             self.__link["SourceName"],
656             self.__link["MeasuresNames"],
657             self.__link["Variable"],
658             self.__link["Format"],
659             )
660         if self.__verbose:
661             print("  [VERBOSE] Measures read")
662             print("  [VERBOSE] Links read")
663             print("  [VERBOSE] Optimization information read")
664             print("  [VERBOSE] Background read:",Background)
665             print("  [VERBOSE] Bounds read:",Params['Bounds'])
666         #
667         if "DifferentialIncrement" not in Params:
668             Params["DifferentialIncrement"] = 0.001
669         if "EnableMultiProcessing" not in Params:
670             Params["EnableMultiProcessing"] = 0
671         if "NumberOfProcesses" not in Params:
672             Params["NumberOfProcesses"] = 0
673         if "StoreSupplementaryCalculations" not in Params:
674             Params["StoreSupplementaryCalculations"] = ["SimulatedObservationAtOptimum",]
675         if "StoreSupplementaryCalculations" in Params and \
676             "SimulatedObservationAtOptimum" not in Params["StoreSupplementaryCalculations"]:
677             Params["StoreSupplementaryCalculations"].append("SimulatedObservationAtOptimum")
678         if self.__verbose and "StoreSupplementaryCalculations" in Params and \
679             "CurrentOptimum" not in Params["StoreSupplementaryCalculations"]:
680             Params["StoreSupplementaryCalculations"].append("CurrentOptimum")
681         if self.__verbose and "StoreSupplementaryCalculations" in Params and \
682             "CostFunctionJAtCurrentOptimum" not in Params["StoreSupplementaryCalculations"]:
683             Params["StoreSupplementaryCalculations"].append("CostFunctionJAtCurrentOptimum")
684
685         dict_dsin_path_ref ={} #creation of the dict_dsin_ref (for storing the path of the dsfinal file for each dataset
686         #
687
688         VariablesToCalibrate = list(BNames)
689
690         def exedymosimMultiobs_REF_CALCULATION( x_values_matrix, dict_dsin_paths = dict_dsin_path_ref): #2ème argument à ne pas modifier
691             "Appel du modèle et restitution des résultats"
692
693             x_values = list(numpy.ravel(x_values_matrix))
694             x_inputs = dict(zip(VariablesToCalibrate,x_values))
695             dict_inputs_for_CWP = {}
696
697 #            multi_y_values_matrix = []
698             for etat in range(len(LNames)):
699                 simudir = self._setModelTmpDir_REF_CALCULATION(ModelName, ModelFormat, LNames[etat])
700
701                 try: #to handle situations in which there is only one CL (boundary conditions) file
702                     var_int = numpy.ravel(LColumns[:,etat])
703                 except :
704                     var_int = numpy.ravel(LColumns)
705                 dict_inputs = dict(zip(LVariablesToChange, var_int))
706
707                 dict_inputs.update( x_inputs )
708
709                 auto_simul = Automatic_Simulation(
710                     simu_dir=simudir,
711                     dymosim_path = os.path.join(simudir, os.pardir, os.pardir) ,
712                     dsres_path = os.path.join(simudir, str("dsres" +LNames[etat]+ ".mat")),
713                     dict_inputs=dict_inputs,
714                     logfile=True,
715                     timeout=TimeoutModelExecution,
716                     without_modelicares=True,
717                     linux = Linux)
718                 auto_simul.single_simulation()
719
720
721
722                 if AdvancedDebugModel:
723                     dslog_file = os.path.join(simudir, 'dslog.txt')
724                     debug_model_file = os.path.join(simudir, os.pardir, 'log_debug_model.txt')
725                     try:
726                         shutil.copy(dslog_file,debug_model_file)
727                     except:
728                         pass
729
730                 if auto_simul.success_code == 2 :
731                     raise ValueError("The simulation falis after initialization, please check the log file and your model in order to be sure that the whole simulation can be performed (specially for dynamic models)")
732
733                 if auto_simul.success_code == 3 :
734                     raise ValueError("The simulation did not reach the final time, probably due to a lack of time, please increase the time for model execution by modifying the TimeoutModelExecution value in configuration.py")
735
736                 if auto_simul.success_code == 0:
737
738                     path_around_simu_no_CWP = os.path.join(simudir, os.pardir, os.pardir)
739                     around_simu_no_CWP = Around_Simulation(dymola_version='2012',
740                                       curdir=path_around_simu_no_CWP,
741                                       source_list_iter_var = 'from_dymtranslog',
742                                       copy_from_dym_trans_log= os.path.join(path_around_simu_no_CWP,'ini.txt'),
743                                       mat = os.path.join(simudir, str("dsres" +LNames[etat]+ ".mat")),
744                                       iter_var_values_options={'source':'from_mat', 'moment':'initial'},
745                                       verbose = False)
746                     around_simu_no_CWP.set_list_iter_var()
747                     around_simu_no_CWP.set_dict_iter_var()
748
749                     writer_no_CWP = write_in_dsin.Write_in_dsin(dict_inputs =around_simu_no_CWP.dict_iter_var,filedir=simudir,dsin_name='dsin.txt',old_file_name='old_dsin.txt',new_file_name=str("dsin_updated_" +LNames[etat]+ ".txt"))
750                     writer_no_CWP.write_in_dsin
751                     writer_no_CWP.write_in_dsin()
752
753                     dict_dsin_paths[LNames[etat]] = os.path.join(simudir,str("dsin_updated_" +LNames[etat]+ ".txt"))
754
755                 else:
756                     from pst4mod.modelica_libraries.change_working_point import Dict_Var_To_Fix
757                     from pst4mod.modelica_libraries.change_working_point import Working_Point_Modification
758
759                     path_around_simu = os.path.join(simudir,os.path.pardir,os.path.pardir)
760
761                     if not(os.path.exists(os.path.join(path_around_simu, str("ref" + ".mat")))):
762                                         auto_simul_ref = Automatic_Simulation(
763                                                                         simu_dir=path_around_simu,
764                                                                         dymosim_path = os.path.join(simudir, os.pardir, os.pardir) ,
765                                                                         dict_inputs={},
766                                                                         logfile=True,
767                                                                         dsres_path = os.path.join(path_around_simu, str("ref" + ".mat")),
768                                                                         timeout=TimeoutModelExecution,
769                                                                         without_modelicares=True,
770                                                                         linux = Linux)
771                                         auto_simul_ref.single_simulation()
772
773                                         if not(auto_simul_ref.success_code == 0):
774                                             raise ValueError("A working dsin.txt file must be provided, provide a dsin.txt that makes it possible for the initial simulation to converge")
775
776                                         temporary_files_to_remove = [os.path.join(path_around_simu,'dsfinal.txt'),
777                                                                      os.path.join(path_around_simu,'dsin_old.txt'),
778                                                                      os.path.join(path_around_simu,'dslog.txt'),
779                                                                      os.path.join(path_around_simu,'status'),
780                                                                      os.path.join(path_around_simu,'success')
781                                                                      ]
782                                         for path_to_remove in temporary_files_to_remove:
783                                             if os.path.exists(path_to_remove):
784                                                 os.remove(path_to_remove)
785
786                     around_simu = Around_Simulation(dymola_version='2012',
787                                       curdir=path_around_simu,
788                                       source_list_iter_var = 'from_dymtranslog',
789                                       copy_from_dym_trans_log= os.path.join(path_around_simu,'ini.txt'),
790                                       mat = os.path.join(path_around_simu, str("ref" + ".mat")),
791                                       iter_var_values_options={'source':'from_mat'},
792                                       verbose = False)
793                     around_simu.set_list_iter_var()
794                     around_simu.set_dict_iter_var()
795
796                     reader_CWP = Reader(os.path.join(path_around_simu, str("ref" + ".mat")),'dymola')
797                     x_values_intemediary = [reader_CWP.values(x_name)[1][-1] for x_name in VariablesToCalibrate]
798                     x_values_from_ref_calculation = numpy.asarray(x_values_intemediary)
799                     x_inputs_from_ref_calculation = dict(zip(VariablesToCalibrate,x_values_from_ref_calculation))
800
801                     cl_values_intemediary = [reader_CWP.values(cl_name)[1][-1] for cl_name in LVariablesToChange]
802                     cl_values_from_ref_calculation = numpy.asarray(cl_values_intemediary)
803                     cl_inputs_from_ref_calculation = dict(zip(LVariablesToChange,cl_values_from_ref_calculation))
804
805
806                     for var_calib in x_inputs_from_ref_calculation.keys():
807                         dict_inputs_for_CWP[var_calib] = (x_inputs_from_ref_calculation[var_calib], x_inputs[var_calib])
808
809                     try: #to handle situations in which there is only one CL (boundary conditions) file
810                         var_int = numpy.ravel(LColumns[:,etat])
811                     except :
812                         var_int = numpy.ravel(LColumns)
813
814                     dict_inputs = dict(zip(LVariablesToChange, var_int))
815
816                     for var_cl in cl_inputs_from_ref_calculation.keys():
817                         dict_inputs_for_CWP[var_cl] = (cl_inputs_from_ref_calculation[var_cl], dict_inputs[var_cl])
818
819                     dict_var_to_fix2 = Dict_Var_To_Fix(option='automatic',
820                                           dict_auto_var_to_fix=dict_inputs_for_CWP)
821
822                     dict_var_to_fix2.set_dict_var_to_fix()
823
824                     if Verbose:
825                         LOG_FILENAME = os.path.join(simudir,'change_working_point.log')
826                         for handler in logging.root.handlers[:]:
827                             logging.root.removeHandler(handler)
828                         logging.basicConfig(filename=LOG_FILENAME,level=logging.DEBUG,filemode='w')
829
830                     working_point_modif = Working_Point_Modification(main_dir = simudir,
831                                  simu_material_dir = simudir,
832                                  dymosim_path = os.path.join(simudir, os.pardir, os.pardir),
833                                  simu_dir = 'SIMUS',
834                                  store_res_dir = 'RES',
835                                  dict_var_to_fix = dict_var_to_fix2.dict_var_to_fix,
836                                  around_simulation0 = around_simu,
837                                  var_to_follow_path= os.path.join(simudir,'var_to_follow.csv'),
838                                  gen_scripts_ini= False,
839                                  nit_max= 1000000000000000,
840                                  min_step_val = 0.00000000000005,
841                                  timeout = TimeoutModelExecution,
842                                  linux = Linux)
843
844                     working_point_modif.create_working_directory()
845                     working_point_modif.working_point_modification(skip_reference_simulation=True)
846
847                     if AdvancedDebugModel:
848                         dslog_file = os.path.join(simudir, 'SIMUS', 'dslog.txt')
849                         debug_model_file = os.path.join(simudir, os.pardir, 'log_debug_model.txt')
850                         try:
851                             shutil.copy(dslog_file,debug_model_file)
852                         except:
853                             pass
854
855                     if not(os.path.exists(os.path.join(simudir, 'RES','1.0.mat'))):
856                         raise ValueError(str("Simulation with Background values does not converge automatically, try to give a dsin.txt file that makes it possible to run a simulation with boundary conditions of " + LNames[etat]))
857
858                     shutil.copy(
859                         os.path.join(simudir, 'RES','1.0.mat'),
860                         os.path.join(simudir, str("dsres" +LNames[etat]+ ".mat"))
861                         )
862
863                     path_around_simu_after_CWP = os.path.join(simudir,os.path.pardir,os.path.pardir)
864                     around_simu_after_CWP = Around_Simulation(dymola_version='2012',
865                                       curdir=path_around_simu_after_CWP,
866                                       source_list_iter_var = 'from_dymtranslog',
867                                       copy_from_dym_trans_log= os.path.join(path_around_simu_after_CWP,'ini.txt'),
868                                       mat = os.path.join(simudir, 'RES','1.0.mat'),  #Le fichier .mat qui vient d'être créé et pour lequel on atteint les valeurs cibles de la variable
869                                       iter_var_values_options={'source':'from_mat', 'moment':'initial'},
870                                       verbose = False)
871                     around_simu_after_CWP.set_list_iter_var()
872                     around_simu_after_CWP.set_dict_iter_var()
873
874                     writer_after_CWP = write_in_dsin.Write_in_dsin(dict_inputs =around_simu_after_CWP.dict_iter_var,filedir=os.path.join(simudir, 'SIMUS'),dsin_name='dsin.txt',old_file_name='old_dsin.txt',new_file_name='dsin_updated.txt')
875                     writer_after_CWP.write_in_dsin
876                     writer_after_CWP.write_in_dsin()
877
878                     shutil.copy(
879                         os.path.join(simudir, 'SIMUS','dsin_updated.txt'),
880                         os.path.join(simudir, str("dsin_updated_" +LNames[etat]+ ".txt"))
881                         )
882                     dict_dsin_paths[LNames[etat]] = os.path.join(simudir,str("dsin_updated_" +LNames[etat]+ ".txt"))
883
884                 os.rename(
885                 os.path.join(os.path.dirname(dict_dsin_paths[LNames[etat]]), str("dsres" +LNames[etat]+ ".mat")),
886                 os.path.join(os.path.dirname(dict_dsin_paths[LNames[etat]]), str("REF_for_CWP" + ".mat"))
887                         )
888
889             return
890
891
892
893         def preparation_exedymosimMultiobs_simple():
894             from pst4mod.modelica_libraries import write_in_dsin
895
896             for etat in range(len(LNames)):
897                 simudir = self._setModelTmpDir_simple(ModelName, ModelFormat, str("dsin_" +LNames[etat]+ ".txt"))
898                 #
899                 try: #to handle situations in which there is only one CL (boundary conditions) file
900                     var_int = numpy.ravel(LColumns[:,etat])
901                 except :
902                     var_int = numpy.ravel(LColumns)
903                 dict_inputs = dict(zip(LVariablesToChange, var_int))
904                 writer_no_CWP = write_in_dsin.Write_in_dsin(dict_inputs = dict_inputs, filedir=simudir, dsin_name=str("dsin_" +LNames[etat]+ ".txt"), old_file_name='old_dsin.txt',new_file_name = str("dsin_" +LNames[etat]+ ".txt"))
905                 writer_no_CWP.write_in_dsin
906                 writer_no_CWP.write_in_dsin()
907             return
908
909         def preparation_exefmu_om_Multiobs():
910             self._setModelTmpDir_simple(ModelName, ModelFormat, ModelName)
911
912         def exepythonMultiobs( x_values_matrix ):
913             "Appel du modèle et restitution des résultats"
914             #
915             x_values = list(numpy.ravel(x_values_matrix))
916             x_inputs = dict(zip(VariablesToCalibrate,x_values))
917             #
918             multi_y_values_matrix = []
919             for etat in range(len(LNames)):
920                 simudir = self._setModelTmpDir(ModelName, ModelFormat)
921                 #
922                 dict_inputs = dict(zip(LVariablesToChange, numpy.ravel(LColumns[:,etat])))
923                 dict_inputs.update( x_inputs )
924                 #
925                 auto_simul = Python_Simulation(
926                     simu_dir=simudir,
927                     val_inputs=x_values_matrix,
928                     logfile=True,
929                     timeout=TimeoutModelExecution,
930                     verbose=self.__verbose)
931                 y_values_matrix = auto_simul.single_simulation()
932                 multi_y_values_matrix.append( y_values_matrix )
933             if not Verbose:
934                 shutil.rmtree(simudir, ignore_errors=True)
935             y_values_matrix = numpy.ravel(numpy.array(multi_y_values_matrix))
936             return y_values_matrix
937         #
938
939
940         #Test for output variables - not repetition  - This test must be before the redefinition of OutputVariables
941         OutputVariables_set = set(OutputVariables)
942         if len(OutputVariables) != len(OutputVariables_set):
943             diff_number_values = len(OutputVariables) - len(OutputVariables_set)
944             raise ValueError("There are " + str(diff_number_values) + " repeated output variables in the definition of OutputVariables, please remove repeated names.")
945
946         #Handle several times same obs -START
947         OutputVariables_new_tmp = []
948
949         whole_obs_in_fileobs = readObsnamesfile(self.__measures["SourceName"])
950
951         for outputname in OutputVariables:
952             if outputname not in whole_obs_in_fileobs:
953                 print("  [VERBOSE] /!\\ ", outputname," is not defined in the measurements file")
954         for obsname in whole_obs_in_fileobs:
955             if obsname in OutputVariables:
956                 OutputVariables_new_tmp.append(obsname)
957
958         OutputVariables_new = []
959         for obs_name in OutputVariables: #to get the same order as for measures
960             count_obs_name = OutputVariables_new_tmp.count(obs_name)
961             for i in range(count_obs_name):
962                 OutputVariables_new.append(obs_name)
963         OutputVariables = OutputVariables_new
964
965         ONames = OutputVariables
966         #Handle several times same obs - END
967
968
969         ODates,OHours = readMultiDatesHours(self.__measures["SourceName"])
970
971         List_Multideltatime =[]
972
973         for etat in range(len(ODates)):
974             time_initial = datetime.strptime(ODates[etat][0]+ ' ' + OHours[etat][0], '%d/%m/%Y %H:%M:%S')
975             list_deltatime = []
976
977             if len(ODates[etat])>1:
978                 list_deltatime.append(0)
979                 for index_tmp in range(len(ODates[etat])-1): #the first date is considered as reference
980                     time_index_tmp = datetime.strptime(ODates[etat][index_tmp+1] + ' ' + OHours[etat][index_tmp+1], '%d/%m/%Y %H:%M:%S')
981                     list_deltatime.append((time_index_tmp-time_initial).total_seconds())
982
983                     if (time_index_tmp-time_initial).total_seconds() < 0:
984                          raise ValueError('The initial date is not chronological for state %s'%str(etat))
985
986             List_Multideltatime.append(list_deltatime)
987
988
989         #Bricolage Debut
990         Observations_new =[]
991         number_observations_by_cl = [] #to be used in the test function, number of measurements by boundary condition (i.e. number of lines in the measurements file)
992
993         for etat in range(len(ODates)): #correspond au nombre de conditions aux limites
994             Obs_etat = Observations[etat].tolist()
995             list_obs_etat = []
996             for obs in range(len(list(ONames))):
997                 if len(ODates[etat])==1: #cas classique, statique
998                     list_obs_etat = Obs_etat
999                     if not(isinstance(list_obs_etat, list)):
1000                         list_obs_etat = [list_obs_etat]
1001
1002                 else:
1003                     for time_etat in range(len(Obs_etat)):
1004                         if not(isinstance(Obs_etat[time_etat], list)): #une seule grandeur observée
1005                             list_obs_etat.append(Obs_etat[time_etat])
1006                         else:
1007                             list_obs_etat.append(Obs_etat[time_etat][obs])
1008
1009             #EXTRA for verify tests - START
1010             if len(ODates[etat])==1:
1011                 number_observations_by_cl.append(1) #in this case only one line in the measurements file, this is just to be used as input of the test functions
1012             else:
1013                 number_observations_by_cl.append(len(Obs_etat)) #we get the number of lines of the measurements file, this is just to be used as input of the test functions
1014             #EXTRA for verify tests - END
1015
1016             Observations_new = Observations_new + list_obs_etat
1017         Observations = Observations_new
1018         #Bricolage Fin
1019
1020
1021
1022
1023         def prev_adao_version_TOP_LEVEL_exedymosimMultiobs( x_values_matrix):
1024             y_values_matrix = TOP_LEVEL_exedymosimMultiobs( x_values_matrix=x_values_matrix, VariablesToCalibrate=VariablesToCalibrate, OutputVariables=OutputVariables, LNames = LNames, ModelFormat = ModelFormat, KeepCalculationFolders = KeepCalculationFolders, Verbose = Verbose, dict_dsin_paths = dict_dsin_path_ref, Linux = Linux, List_Multideltatime = List_Multideltatime, AdvancedDebugModel = AdvancedDebugModel, TimeoutModelExecution = TimeoutModelExecution)
1025             return y_values_matrix
1026
1027         def prev_adao_version_TOPLEVEL_exedymosimMultiobs_simple(x_values_matrix):
1028             y_values_matrix = TOPLEVEL_exedymosimMultiobs_simple(x_values_matrix = x_values_matrix, VariablesToCalibrate=VariablesToCalibrate, OutputVariables=OutputVariables, LNames = LNames, ref_simudir = self._get_Name_ModelTmpDir_simple(ModelName), KeepCalculationFolders = KeepCalculationFolders, Linux = Linux, List_Multideltatime = List_Multideltatime, AdvancedDebugModel = AdvancedDebugModel, TimeoutModelExecution = TimeoutModelExecution)
1029             return y_values_matrix
1030
1031         if self.__model["Format"].upper() in ["DYMOSIM", "GUESS"]:
1032
1033             if adao.version[:5] < '9.7.0' and adao.version[5]=='.': #to handle differences in the way directoperator is mangaged after modifications brought to direct operator
1034                 if ComplexModel:
1035                     simulateur = prev_adao_version_TOP_LEVEL_exedymosimMultiobs
1036                     exedymosimMultiobs_REF_CALCULATION(Background)
1037                 else:
1038                     simulateur = prev_adao_version_TOPLEVEL_exedymosimMultiobs_simple
1039                     preparation_exedymosimMultiobs_simple()
1040
1041             else:
1042                 if ComplexModel:
1043                     simulateur = TOP_LEVEL_exedymosimMultiobs #to handle parallel computing
1044                     exedymosimMultiobs_REF_CALCULATION(Background)
1045                 else:
1046                     simulateur = TOPLEVEL_exedymosimMultiobs_simple #to handle parallel computing
1047                     preparation_exedymosimMultiobs_simple()
1048
1049         elif self.__model["Format"].upper() == "PYSIM":
1050             simulateur = exepythonMultiobs
1051         elif self.__model["Format"].upper() == "FMI" or self.__model["Format"].upper() == "FMU":
1052             simulateur = TOP_LEVEL_exefmuMultiobs
1053             preparation_exefmu_om_Multiobs()
1054         elif self.__model["Format"].upper() == "OPENMODELICA":
1055             preparation_exefmu_om_Multiobs()
1056             simulateur = TOP_LEVEL_exeOpenModelicaMultiobs
1057
1058         else:
1059             raise ValueError("Unknown model format \"%s\""%self.__model["Format"].upper())
1060
1061         def verify_function(x_values, model_format):
1062
1063
1064             simulation_results_all = []
1065             list_simulation_time = []
1066             for i in list(range(3)): #we consider only three iterations
1067                 start_time_verify = time.time()
1068                 if model_format in ["DYMOSIM"]:
1069                     if ComplexModel:
1070                         simulation_results = prev_adao_version_TOP_LEVEL_exedymosimMultiobs(x_values)
1071                     else:
1072                         simulation_results = prev_adao_version_TOPLEVEL_exedymosimMultiobs_simple(x_values)
1073                 elif model_format in ["FMI","FMU"]:
1074                     simulation_results = TOP_LEVEL_exefmuMultiobs(x_values, VariablesToCalibrate=VariablesToCalibrate, OutputVariables=OutputVariables,  LNames = LNames, LColumns = LColumns, LVariablesToChange=LVariablesToChange, ref_simudir = self._get_Name_ModelTmpDir_simple(ModelName),  ModelName = ModelName, List_Multideltatime = List_Multideltatime, AdvancedDebugModel = AdvancedDebugModel, TimeoutModelExecution = TimeoutModelExecution, FMUInput = FMUInput)
1075                 elif model_format in ["OPENMODELICA"]:
1076                     simulation_results = TOP_LEVEL_exeOpenModelicaMultiobs(x_values, KeepCalculationFolders = KeepCalculationFolders, VariablesToCalibrate=VariablesToCalibrate, OutputVariables=OutputVariables,  LNames = LNames, LColumns = LColumns, LVariablesToChange=LVariablesToChange, ref_simudir = self._get_Name_ModelTmpDir_simple(ModelName),  ModelName = ModelName, List_Multideltatime = List_Multideltatime, Linux = Linux, AdvancedDebugModel = AdvancedDebugModel, TimeoutModelExecution = TimeoutModelExecution)
1077                 else:
1078                     raise NotImplementedError("Not yet implemented for current model format: ", model_format)
1079                 end_time_verify = time.time()
1080
1081                 list_simulation_time.append(round(end_time_verify - start_time_verify,3))
1082                 simulation_results_all.append(simulation_results)
1083
1084
1085             elapsed_time = round(sum(list_simulation_time)/len(list_simulation_time),3)
1086                 #not necessary, already in good format
1087             # numpy.set_printoptions(precision=2,linewidth=5000, threshold = 10000)
1088             # reshaped_simulation_results_all = numpy.array(simulation_results_all).reshape((len(OutputVariables),-1)).transpose()
1089
1090             # numpy.set_printoptions(precision=2,linewidth=5000)
1091
1092             list_check_output = []
1093             list_false_output = []
1094             simulation_results_all_one_boundary_condition =[]
1095
1096             for simulation_output in simulation_results_all: #iteration over the three simulations performed
1097                 # THIS DOES NOT WORK IF OBSERVATIONS FILES ARE OF DIFFERENT SIZE# simultation_ouput_reshaped = simulation_output.reshape((len(LNames),-1))[0]    #consideration of only one set of boundary conditions (otherwise outputs could be modified by boundary condition modifications)
1098                 total_number_measurements_first_simulation = number_observations_by_cl[0]*len(OutputVariables)
1099                 first_simulation_results = simulation_output[:total_number_measurements_first_simulation]
1100                 simulation_results_all_one_boundary_condition.append(first_simulation_results)
1101
1102             for index_line_time_simulation in range(len(simulation_results_all_one_boundary_condition[0].reshape(len(OutputVariables),-1).transpose())): #iteration over each time results individually
1103                 for i in range(len(OutputVariables)):
1104                     if (simulation_results_all_one_boundary_condition[0].reshape(len(OutputVariables),-1).transpose()[index_line_time_simulation][i] == simulation_results_all_one_boundary_condition[1].reshape(len(OutputVariables),-1).transpose()[index_line_time_simulation][i]) & (simulation_results_all_one_boundary_condition[0].reshape(len(OutputVariables),-1).transpose()[index_line_time_simulation][i] == simulation_results_all_one_boundary_condition[2].reshape(len(OutputVariables),-1).transpose()[index_line_time_simulation][i]):
1105                         list_check_output.append("True")
1106
1107                     else:
1108                         list_check_output.append("False")
1109                         if OutputVariables[i] in list_false_output: #problematic variable already included
1110                             pass
1111                         else:
1112                                 list_false_output.append(OutputVariables[i])
1113
1114             if "False" in list_check_output:
1115                 verify_function_result = "False"
1116             else:
1117                 verify_function_result = "True"
1118             return verify_function_result, list_false_output, elapsed_time
1119
1120         def verify_function_sensitivity(x_values, increment, model_format):
1121
1122             if model_format in ["DYMOSIM"]:
1123                 if ComplexModel:
1124                     simulation_results = prev_adao_version_TOP_LEVEL_exedymosimMultiobs(x_values)
1125                 else:
1126                     simulation_results = prev_adao_version_TOPLEVEL_exedymosimMultiobs_simple(x_values)
1127             elif model_format in ["FMI","FMU"]:
1128                 simulation_results = TOP_LEVEL_exefmuMultiobs(x_values, VariablesToCalibrate=VariablesToCalibrate, OutputVariables=OutputVariables,  LNames = LNames, LColumns = LColumns, LVariablesToChange=LVariablesToChange, ref_simudir = self._get_Name_ModelTmpDir_simple(ModelName),  ModelName = ModelName, List_Multideltatime = List_Multideltatime, AdvancedDebugModel = AdvancedDebugModel, TimeoutModelExecution = TimeoutModelExecution, FMUInput = FMUInput)
1129             elif model_format in ["OPENMODELICA"]:
1130                 simulation_results = TOP_LEVEL_exeOpenModelicaMultiobs(x_values, KeepCalculationFolders = KeepCalculationFolders, VariablesToCalibrate=VariablesToCalibrate, OutputVariables=OutputVariables,  LNames = LNames, LColumns = LColumns, LVariablesToChange=LVariablesToChange, ref_simudir = self._get_Name_ModelTmpDir_simple(ModelName),  ModelName = ModelName, List_Multideltatime = List_Multideltatime, Linux = Linux, AdvancedDebugModel = AdvancedDebugModel, TimeoutModelExecution = TimeoutModelExecution)
1131             else:
1132                 raise NotImplementedError("Not yet implemented for current model format: ", model_format)
1133
1134
1135             x_values_modif = [x*(1+increment) for x in x_values]
1136
1137             if model_format in ["DYMOSIM"]:
1138                 if ComplexModel:
1139                     simulation_results_modif = prev_adao_version_TOP_LEVEL_exedymosimMultiobs(x_values_modif)
1140                 else:
1141                     simulation_results_modif = prev_adao_version_TOPLEVEL_exedymosimMultiobs_simple(x_values_modif)
1142             elif model_format in ["FMI","FMU"]:
1143                 simulation_results_modif = TOP_LEVEL_exefmuMultiobs(x_values_modif, VariablesToCalibrate=VariablesToCalibrate, OutputVariables=OutputVariables,  LNames = LNames, LColumns = LColumns, LVariablesToChange=LVariablesToChange, ref_simudir = self._get_Name_ModelTmpDir_simple(ModelName),  ModelName = ModelName, List_Multideltatime = List_Multideltatime, AdvancedDebugModel = AdvancedDebugModel, TimeoutModelExecution = TimeoutModelExecution, FMUInput = FMUInput)
1144             elif model_format in ["OPENMODELICA"]:
1145                 simulation_results_modif = TOP_LEVEL_exeOpenModelicaMultiobs(x_values_modif, KeepCalculationFolders = KeepCalculationFolders, VariablesToCalibrate=VariablesToCalibrate, OutputVariables=OutputVariables,  LNames = LNames, LColumns = LColumns, LVariablesToChange=LVariablesToChange, ref_simudir = self._get_Name_ModelTmpDir_simple(ModelName),  ModelName = ModelName, List_Multideltatime = List_Multideltatime, Linux = Linux, AdvancedDebugModel = AdvancedDebugModel, TimeoutModelExecution = TimeoutModelExecution)
1146             else:
1147                 raise NotImplementedError("Not yet implemented for current model format: ", model_format)
1148
1149             list_check_output = ["False" for i in OutputVariables] #everything initialized as not good and then if something changes (at least one point), then it is OK
1150             list_not_modified_variables =[]
1151
1152             total_number_measurements_first_simulation = number_observations_by_cl[0]*len(OutputVariables)
1153
1154             # THIS DOES NOT WORK IF OBSERVATIONS FILES ARE OF DIFFERENT SIZE# simulation_results_one_boundary_condition = simulation_results.reshape((len(LNames),-1))[0]
1155             # THIS DOES NOT WORK IF OBSERVATIONS FILES ARE OF DIFFERENT SIZE# simulation_results_modif_one_boundary_condition = simulation_results_modif.reshape((len(LNames),-1))[0]
1156             simulation_results_one_boundary_condition = simulation_results[:total_number_measurements_first_simulation]
1157             simulation_results_modif_one_boundary_condition = simulation_results_modif[:total_number_measurements_first_simulation]
1158
1159             for index_line_time_simulation in range(len(simulation_results_one_boundary_condition.reshape(len(OutputVariables),-1).transpose())):
1160
1161                 for i in range(len(OutputVariables)):
1162                     if simulation_results_one_boundary_condition.reshape(len(OutputVariables),-1).transpose()[index_line_time_simulation][i] != simulation_results_modif_one_boundary_condition.reshape(len(OutputVariables),-1).transpose()[index_line_time_simulation][i] :
1163                         list_check_output[i] = "True"
1164
1165
1166             for i in range(len(OutputVariables)):
1167                 if list_check_output[i] == "False": #it means the variable has not changed
1168                     list_not_modified_variables.append(OutputVariables[i])
1169
1170             if "False" in list_check_output:
1171                 verify_function_sensitivity_result = "False"
1172             else:
1173                 verify_function_sensitivity_result = "True"
1174
1175             return verify_function_sensitivity_result, list_not_modified_variables
1176
1177         def verify_function_sensitivity_parameter(x_values, increment, model_format):
1178
1179             if model_format in ["DYMOSIM"]:
1180                 if ComplexModel:
1181                     simulation_results_ref = prev_adao_version_TOP_LEVEL_exedymosimMultiobs(x_values)
1182                 else:
1183                     simulation_results_ref = prev_adao_version_TOPLEVEL_exedymosimMultiobs_simple(x_values)
1184             elif model_format in ["FMI","FMU"]:
1185                 simulation_results_ref = TOP_LEVEL_exefmuMultiobs(x_values, VariablesToCalibrate=VariablesToCalibrate, OutputVariables=OutputVariables,  LNames = LNames, LColumns = LColumns, LVariablesToChange=LVariablesToChange, ref_simudir = self._get_Name_ModelTmpDir_simple(ModelName),  ModelName = ModelName, List_Multideltatime = List_Multideltatime, AdvancedDebugModel = AdvancedDebugModel, TimeoutModelExecution = TimeoutModelExecution, FMUInput = FMUInput)
1186             elif model_format in ["OPENMODELICA"]:
1187                 simulation_results_ref = TOP_LEVEL_exeOpenModelicaMultiobs(x_values, KeepCalculationFolders = KeepCalculationFolders, VariablesToCalibrate=VariablesToCalibrate, OutputVariables=OutputVariables,  LNames = LNames, LColumns = LColumns, LVariablesToChange=LVariablesToChange, ref_simudir = self._get_Name_ModelTmpDir_simple(ModelName),  ModelName = ModelName, List_Multideltatime = List_Multideltatime, Linux = Linux, AdvancedDebugModel = AdvancedDebugModel, TimeoutModelExecution = TimeoutModelExecution)
1188             else:
1189                 raise NotImplementedError("Not yet implemented for current model format: ", model_format)
1190
1191             list_param_with_no_impact = []
1192             Check_params_impact = "True"
1193
1194             x_values_modif = [x*(1+increment) for x in x_values]
1195
1196
1197             for param_index in range(len(x_values)):   #the goal is to modify each param individually and check if it is possible to observe variations
1198                 x_values_param = list(x_values) #avoid memory surprises
1199                 x_values_param[param_index] = x_values_modif[param_index]
1200                 if model_format in ["DYMOSIM"]:
1201                     if ComplexModel:
1202                         simulation_results_param = prev_adao_version_TOP_LEVEL_exedymosimMultiobs(x_values_param)
1203                     else:
1204                         simulation_results_param = prev_adao_version_TOPLEVEL_exedymosimMultiobs_simple(x_values_param)
1205                 elif model_format in ["FMI","FMU"]:
1206                     simulation_results_param = TOP_LEVEL_exefmuMultiobs(x_values_param, VariablesToCalibrate=VariablesToCalibrate, OutputVariables=OutputVariables,  LNames = LNames, LColumns = LColumns, LVariablesToChange=LVariablesToChange, ref_simudir = self._get_Name_ModelTmpDir_simple(ModelName),  ModelName = ModelName, List_Multideltatime = List_Multideltatime, AdvancedDebugModel = AdvancedDebugModel, TimeoutModelExecution = TimeoutModelExecution, FMUInput = FMUInput)
1207                 elif model_format in ["OPENMODELICA"]:
1208                     simulation_results_param = TOP_LEVEL_exeOpenModelicaMultiobs(x_values_param, KeepCalculationFolders = KeepCalculationFolders, VariablesToCalibrate=VariablesToCalibrate, OutputVariables=OutputVariables,  LNames = LNames, LColumns = LColumns, LVariablesToChange=LVariablesToChange, ref_simudir = self._get_Name_ModelTmpDir_simple(ModelName),  ModelName = ModelName, List_Multideltatime = List_Multideltatime, Linux = Linux, AdvancedDebugModel = AdvancedDebugModel, TimeoutModelExecution = TimeoutModelExecution)
1209                 else:
1210                     raise NotImplementedError("Not yet implemented for current model format: ", model_format)
1211
1212                 if numpy.array_equal(simulation_results_param,simulation_results_ref): #there is not impact on the whole simulation
1213                     list_param_with_no_impact.append(VariablesToCalibrate[param_index])
1214                     Check_params_impact = "False"
1215
1216             return Check_params_impact, list_param_with_no_impact
1217
1218
1219         def check_model_stability(x_values_names, background_values, bounds, number_of_tests, model_format):
1220
1221             dict_check_model_stability = {}
1222             for index_x in range(len(x_values_names)):
1223                 x_bg_tested = x_values_names[index_x]
1224
1225                 print("  [VERBOSE] Model stability check is being performed for the following variable to be calibrated: ", x_bg_tested)
1226
1227                 if bounds[index_x][0] == None or bounds[index_x][0] == "None":
1228                     bounds_inf_x_bg_tested = background_values[index_x] - abs(background_values[index_x]*100)
1229                 else:
1230                     bounds_inf_x_bg_tested = bounds[index_x][0]
1231
1232                 if bounds[index_x][1] == None or bounds[index_x][1] == "None":
1233                     bounds_sup_x_bg_tested = background_values[index_x] + abs(background_values[index_x]*100)
1234                 else:
1235                     bounds_sup_x_bg_tested = bounds[index_x][1]
1236
1237                 #avoid problems with values to be tested hen the bounds are strictly equal to 0
1238                 if bounds_inf_x_bg_tested == 0: #which means that bounds_sup_x_bg_tested > 0
1239                     bounds_inf_x_bg_tested = min(1e-9,bounds_sup_x_bg_tested*1e-9)
1240                 if bounds_sup_x_bg_tested == 0: #which means that bounds_inf_x_bg_tested < 0
1241                     bounds_sup_x_bg_tested = max(-1e-9,bounds_inf_x_bg_tested*1e-9)
1242
1243                 list_x_bg_tested = []
1244
1245                 if bounds_inf_x_bg_tested < 0:
1246                     if bounds_sup_x_bg_tested<0:
1247                         list_x_bg_tested = list(numpy.geomspace(bounds_inf_x_bg_tested,bounds_sup_x_bg_tested,int(number_of_tests)))
1248                     else:
1249                         list_x_bg_tested = list(numpy.geomspace(bounds_inf_x_bg_tested,bounds_inf_x_bg_tested*1e-4,int(number_of_tests/2)))
1250                         list_x_bg_tested = list_x_bg_tested + list((numpy.geomspace(bounds_sup_x_bg_tested*1e-4,bounds_sup_x_bg_tested,int(number_of_tests/2))))
1251                 else: #both bounds are >0
1252                     list_x_bg_tested = list(numpy.geomspace(bounds_inf_x_bg_tested,bounds_sup_x_bg_tested,int(number_of_tests)))
1253
1254                 #for a linear partition, we consider a log repartition better
1255                 # for index_partition in range(number_of_tests):
1256                 #     list_x_bg_tested.append(bounds_inf_x_bg_tested + index_partition*(bounds_sup_x_bg_tested-bounds_inf_x_bg_tested)/(number_of_tests-1))
1257
1258                 list_failed_model_evaluation = []
1259                 for x_value_bg_tested in list_x_bg_tested:
1260
1261                     x_whole_values_with_bg = copy.deepcopy(background_values) #avoid modifiction of backgroundvalues
1262                     x_whole_values_with_bg[index_x] = x_value_bg_tested
1263                     try:
1264                         if model_format in ["DYMOSIM"]:
1265                             if ComplexModel:
1266                                 simulation_results = prev_adao_version_TOP_LEVEL_exedymosimMultiobs(x_whole_values_with_bg)
1267                             else:
1268                                 simulation_results = prev_adao_version_TOPLEVEL_exedymosimMultiobs_simple(x_whole_values_with_bg)
1269                         elif model_format in ["FMI","FMU"]:
1270                             simulation_results = TOP_LEVEL_exefmuMultiobs(x_whole_values_with_bg, VariablesToCalibrate=VariablesToCalibrate, OutputVariables=OutputVariables,  LNames = LNames, LColumns = LColumns, LVariablesToChange=LVariablesToChange, ref_simudir = self._get_Name_ModelTmpDir_simple(ModelName),  ModelName = ModelName, List_Multideltatime = List_Multideltatime, AdvancedDebugModel = AdvancedDebugModel, TimeoutModelExecution = TimeoutModelExecution, FMUInput = FMUInput)
1271                         elif model_format in ["OPENMODELICA"]:
1272                             simulation_results = TOP_LEVEL_exeOpenModelicaMultiobs(x_whole_values_with_bg, KeepCalculationFolders = KeepCalculationFolders, VariablesToCalibrate=VariablesToCalibrate, OutputVariables=OutputVariables,  LNames = LNames, LColumns = LColumns, LVariablesToChange=LVariablesToChange, ref_simudir = self._get_Name_ModelTmpDir_simple(ModelName),  ModelName = ModelName, List_Multideltatime = List_Multideltatime, Linux = Linux, AdvancedDebugModel = AdvancedDebugModel, TimeoutModelExecution = TimeoutModelExecution)
1273                         else:
1274                             raise NotImplementedError("Not yet implemented for current model format: ", model_format)
1275                     except:
1276                         list_failed_model_evaluation.append(x_value_bg_tested)
1277
1278                 dict_check_model_stability[x_bg_tested] = list_failed_model_evaluation
1279             return dict_check_model_stability
1280
1281         if VerifyFunction:
1282             check_function, list_not_same_model_outputs, elapsed_time_for_simu = verify_function(x_values = Background, model_format = self.__model["Format"].upper())
1283             if check_function == "True":
1284                 print("  [VERBOSE] Results of function checking: ", "OK for all model outputs, ","simulation time (including all boundary conditions) is ",elapsed_time_for_simu, " seconds" )
1285             else:
1286                 print("  [VERBOSE] Results of function checking: ", "NOT OK for all model outputs")
1287                 print("  [VERBOSE] The following model outputs are not the same when repeating a simulation: ", list_not_same_model_outputs )
1288
1289         if GlobalSensitivityCheck:
1290             check_sensitivity, list_same_model_outputs = verify_function_sensitivity(x_values = Background, increment = Params["DifferentialIncrement"], model_format = self.__model["Format"].upper())
1291             if VerifyFunction == False:
1292                 print("  [VERBOSE] /!\\ Activate VerifyFunction option and check that the result of the check is OK to get reliable results from GlobalSensitivityCheck option")
1293             if check_sensitivity == "True":
1294                 print("  [VERBOSE] Results of function global sensitivity analysis: ", "All model outputs vary when background values are modified")
1295             else:
1296                 print("  [VERBOSE] Results of function global sensitivity analysis: ", "Some model outputs DO NOT vary when background values are modified")
1297                 print("  [VERBOSE] The following model outputs DO NOT vary when background values are modified: ", list_same_model_outputs )
1298
1299         if ParamSensitivityCheck:
1300             check_function, list_params_without_impact = verify_function_sensitivity_parameter(x_values = Background, increment = Params["DifferentialIncrement"], model_format = self.__model["Format"].upper())
1301             if VerifyFunction == False:
1302                 print("  [VERBOSE] /!\\ Activate VerifyFunction option and check that the result of the check is OK to get reliable results from ParamSensitivityCheck option")
1303             if check_function == "True":
1304                 print("  [VERBOSE] Results of parameter sensitivity checking: ", "OK, all parameters have an impact on model outputs")
1305             else:
1306                 print("  [VERBOSE] Results of parameter sensitivity checking: ", "NOT OK, some parameters do not have an impact on model outputs")
1307                 print("  [VERBOSE] The following parameters do not have an impact on model outputs: ", list_params_without_impact )
1308
1309         if ModelStabilityCheck:
1310             print("  [VERBOSE] Model stability check is being performed, this check might take several minutes depending on the model and number of variables to calibrate (you can reduce the value of ModelStabilityCheckingPoints option to reduce checking time, by default its value is 10)")
1311             check_stability = check_model_stability(x_values_names = VariablesToCalibrate, background_values = Background, bounds = Params['Bounds'], number_of_tests = ModelStabilityCheckingPoints, model_format = self.__model["Format"].upper())
1312             if check_stability == dict(zip(VariablesToCalibrate,[list() for i in VariablesToCalibrate])):
1313                 print("  [VERBOSE] Results of Model stability check: ", "The model simulates correctly when modifying individually the values of the variables to calibrate within the range in which their optimal value should be found")
1314             else:
1315                 print("  [VERBOSE] Results of Model stability check: ", "The model FAILS to converge with the following values of variables to calibrate", check_stability)
1316                 print("  [VERBOSE] The data assimilation procedure might not succeed if the model fails to converge, please revise the model and/or the range in which the optimal value of the variables to calibrate should be found")
1317
1318         #
1319
1320         __adaocase = adaoBuilder.New()
1321         __adaocase.set( 'AlgorithmParameters',  Algorithm=Algo, Parameters=Params )
1322         __adaocase.set( 'Background',           Vector=Background)
1323         __adaocase.set( 'Observation',          Vector=Observations )
1324         if type(CovB) is float:
1325             __adaocase.set( 'BackgroundError',  ScalarSparseMatrix=CovB )
1326         else:
1327             __adaocase.set( 'BackgroundError',  Matrix=CovB )
1328
1329         if AdaptObservationError:
1330             Square_Observations = [x*x for x  in Observations]
1331
1332             # if (type(CovR) is float and CovR!= 1.0 ):# by default we still consider the square, better not to modify it
1333             if type(CovR) is float:
1334                 ObsError = CovR*numpy.diag(Square_Observations)
1335                 __adaocase.set( 'ObservationError', Matrix = ObsError )
1336             else:
1337                 ObsError = 1e-15*numpy.diag(Square_Observations) #arbitary low value to eventually handle pressures in Pa
1338                 __adaocase.set( 'ObservationError', Matrix = ObsError )
1339             print("  [VERBOSE] AdaptObservationError was set to True, the following ObservationError matrix has been considered: ")
1340             print(ObsError)
1341         else: #classical situation
1342             if type(CovR) is float:
1343                 __adaocase.set( 'ObservationError', ScalarSparseMatrix=CovR )
1344             else:
1345                 __adaocase.set( 'ObservationError', Matrix=CovR )
1346
1347         if Params["EnableMultiProcessing"]==1 :
1348             ParallelComputationGradient = True
1349         else:
1350             ParallelComputationGradient = False
1351
1352
1353         if self.__model["Format"].upper() in ["DYMOSIM"]:
1354             if ParallelComputationGradient:
1355
1356                 if adao.version[:5] < '9.7.0' and adao.version[5]=='.':
1357                     raise ValueError("ADAO version must be at least 9.7.0 to support parallel computation of the gradient, current version is %s"%str(adao.version[:5]))
1358
1359                 if ComplexModel:
1360                       __adaocase.set( 'ObservationOperator',  OneFunction=simulateur,
1361                                Parameters = {"DifferentialIncrement":Params["DifferentialIncrement"], "EnableMultiProcessing":Params["EnableMultiProcessing"], "NumberOfProcesses" : Params["NumberOfProcesses"]},
1362                                ExtraArguments = {'VariablesToCalibrate':VariablesToCalibrate, 'OutputVariables' : OutputVariables, 'LNames': LNames, 'ModelFormat': ModelFormat, 'KeepCalculationFolders': KeepCalculationFolders, 'Verbose' : Verbose, 'dict_dsin_paths' : dict_dsin_path_ref, 'Linux': Linux, 'List_Multideltatime' : List_Multideltatime, 'AdvancedDebugModel':AdvancedDebugModel, 'TimeoutModelExecution':TimeoutModelExecution   },
1363                                )
1364                 else:
1365                     __adaocase.set( 'ObservationOperator',  OneFunction=simulateur,
1366                                Parameters = {"DifferentialIncrement":Params["DifferentialIncrement"], "EnableMultiProcessing":Params["EnableMultiProcessing"], "NumberOfProcesses" : Params["NumberOfProcesses"]},
1367                                ExtraArguments = {'VariablesToCalibrate':VariablesToCalibrate, 'OutputVariables' : OutputVariables, 'LNames': LNames, 'ref_simudir': self._get_Name_ModelTmpDir_simple(ModelName), 'KeepCalculationFolders': KeepCalculationFolders, 'Linux': Linux, 'List_Multideltatime' : List_Multideltatime, 'AdvancedDebugModel':AdvancedDebugModel, 'TimeoutModelExecution':TimeoutModelExecution  },
1368                                )
1369
1370             else:
1371                 if adao.version[:5] < '9.7.0' and adao.version[5]=='.':
1372                     __adaocase.set( 'ObservationOperator',  OneFunction=simulateur, Parameters = {"DifferentialIncrement":Params["DifferentialIncrement"]})
1373
1374                 else:
1375                     if ComplexModel:
1376                           __adaocase.set( 'ObservationOperator',  OneFunction=simulateur,
1377                                    Parameters = {"DifferentialIncrement":Params["DifferentialIncrement"]},
1378                                    ExtraArguments = {'VariablesToCalibrate':VariablesToCalibrate, 'OutputVariables' : OutputVariables, 'LNames': LNames, 'ModelFormat': ModelFormat, 'KeepCalculationFolders': KeepCalculationFolders, 'Verbose' : Verbose, 'dict_dsin_paths' : dict_dsin_path_ref, 'Linux': Linux, 'List_Multideltatime' : List_Multideltatime, 'AdvancedDebugModel':AdvancedDebugModel, 'TimeoutModelExecution':TimeoutModelExecution  },
1379                                    )
1380                     else:
1381                         __adaocase.set( 'ObservationOperator',  OneFunction=simulateur,
1382                                    Parameters = {"DifferentialIncrement":Params["DifferentialIncrement"]},
1383                                    ExtraArguments = {'VariablesToCalibrate':VariablesToCalibrate, 'OutputVariables' : OutputVariables, 'LNames': LNames, 'ref_simudir': self._get_Name_ModelTmpDir_simple(ModelName), 'KeepCalculationFolders': KeepCalculationFolders, 'Linux': Linux, 'List_Multideltatime' : List_Multideltatime, 'AdvancedDebugModel':AdvancedDebugModel, 'TimeoutModelExecution':TimeoutModelExecution  })
1384
1385         elif self.__model["Format"].upper() in ["OPENMODELICA"]: #OPENMODELICA (different from FMI since there are different keywords for the function: Linux and KeepCalculationFolders)
1386
1387             if ComplexModel:
1388                 print("ComplexModel option is only valid for DYMOSIM format (it has no effect for the current calculation)")
1389
1390             if ParallelComputationGradient:
1391                 if adao.version[:5] < '9.7.0' and adao.version[5]=='.':
1392                     raise ValueError("ADAO version must be at least 9.7.0 to support parallel computation of the gradient, current version is %s"%str(adao.version[:5]))
1393                 else:
1394                     __adaocase.set( 'ObservationOperator',  OneFunction=simulateur,
1395                                Parameters = {"DifferentialIncrement":Params["DifferentialIncrement"]},
1396                                ExtraArguments = {'VariablesToCalibrate':VariablesToCalibrate, 'OutputVariables' : OutputVariables, 'LNames': LNames, 'LColumns':LColumns, 'LVariablesToChange':LVariablesToChange, 'ref_simudir': self._get_Name_ModelTmpDir_simple(ModelName), 'ModelName': ModelName, 'List_Multideltatime' : List_Multideltatime, 'Linux': Linux, 'KeepCalculationFolders': KeepCalculationFolders, 'AdvancedDebugModel':AdvancedDebugModel, 'TimeoutModelExecution':TimeoutModelExecution})
1397
1398             else:
1399                 if adao.version[:5] < '9.7.0' and adao.version[5]=='.':
1400                     raise ValueError("ADAO version must be at least 9.7.0 to support other formats than DYMOSIM")
1401                 else:
1402                     __adaocase.set( 'ObservationOperator',  OneFunction=simulateur,
1403                                Parameters = {"DifferentialIncrement":Params["DifferentialIncrement"]},
1404                                ExtraArguments = {'VariablesToCalibrate':VariablesToCalibrate, 'OutputVariables' : OutputVariables, 'LNames': LNames, 'LColumns':LColumns, 'LVariablesToChange':LVariablesToChange, 'ref_simudir': self._get_Name_ModelTmpDir_simple(ModelName), 'ModelName': ModelName, 'List_Multideltatime' : List_Multideltatime, 'Linux': Linux, 'KeepCalculationFolders': KeepCalculationFolders, 'AdvancedDebugModel':AdvancedDebugModel, 'TimeoutModelExecution':TimeoutModelExecution})
1405
1406         else: #FMI or GUESS format
1407
1408             if ComplexModel:
1409                 print("ComplexModel option is only valid for DYMOSIM format (it has no effect for the current calculation)")
1410
1411             if ParallelComputationGradient:
1412                 if adao.version[:5] < '9.7.0' and adao.version[5]=='.':
1413                     raise ValueError("ADAO version must be at least 9.7.0 to support parallel computation of the gradient, current version is %s"%str(adao.version[:5]))
1414                 else:
1415                     __adaocase.set( 'ObservationOperator',  OneFunction=simulateur,
1416                                Parameters = {"DifferentialIncrement":Params["DifferentialIncrement"]},
1417                                ExtraArguments = {'VariablesToCalibrate':VariablesToCalibrate, 'OutputVariables' : OutputVariables, 'LNames': LNames, 'LColumns':LColumns, 'LVariablesToChange':LVariablesToChange, 'ref_simudir': self._get_Name_ModelTmpDir_simple(ModelName), 'ModelName': ModelName, 'List_Multideltatime' : List_Multideltatime, 'AdvancedDebugModel':AdvancedDebugModel, 'TimeoutModelExecution':TimeoutModelExecution, 'FMUInput':FMUInput})
1418
1419             else:
1420                 if adao.version[:5] < '9.7.0' and adao.version[5]=='.':
1421                     raise ValueError("ADAO version must be at least 9.7.0 to support other formats than DYMOSIM")
1422                 else:
1423                     __adaocase.set( 'ObservationOperator',  OneFunction=simulateur,
1424                                Parameters = {"DifferentialIncrement":Params["DifferentialIncrement"]},
1425                                ExtraArguments = {'VariablesToCalibrate':VariablesToCalibrate, 'OutputVariables' : OutputVariables, 'LNames': LNames, 'LColumns':LColumns, 'LVariablesToChange':LVariablesToChange, 'ref_simudir': self._get_Name_ModelTmpDir_simple(ModelName), 'ModelName': ModelName, 'List_Multideltatime' : List_Multideltatime, 'AdvancedDebugModel':AdvancedDebugModel, 'TimeoutModelExecution':TimeoutModelExecution, 'FMUInput':FMUInput})
1426
1427
1428         #
1429         if self.__verbose:
1430             __adaocase.set( 'Observer', Variable="CostFunctionJAtCurrentOptimum", String="print('\\n  -----------------------------------\\n  [VERBOSE] Current iteration number: %i'%len(var))" )
1431             __adaocase.set( 'Observer', Variable="CostFunctionJAtCurrentOptimum",  Template="ValuePrinter", Info="\n  [VERBOSE] Current optimal cost value.:" )
1432             __adaocase.set( 'Observer', Variable="CurrentOptimum",                 Template="ValuePrinter", Info="\n  [VERBOSE] Current optimal state......:" )        #
1433
1434         if self.__IntermediateInformation:
1435             __adaocase.set( 'Observer', Variable="CostFunctionJ",  Template="ValuePrinter", Info="\n  [VERBOSE] Current cost value.:" )
1436             __adaocase.set( 'Observer', Variable="CurrentState",  Template="ValuePrinter", Info="\n  [VERBOSE] Current state......:" )
1437
1438         __adaocase.execute()
1439         #
1440
1441         if InitialSimulation: #the prev_adao_version_TOP_LEVEL_XXX is already defined with proper arguments (boundary conditions etc.)
1442
1443             if self.__model["Format"].upper() in ["DYMOSIM"]:
1444                 if ComplexModel:
1445                     initialsimulation_results = prev_adao_version_TOP_LEVEL_exedymosimMultiobs(Background)
1446                 else:
1447                     initialsimulation_results = prev_adao_version_TOPLEVEL_exedymosimMultiobs_simple(Background)
1448             elif self.__model["Format"].upper() in ["FMI", "FMU"]:
1449                 initialsimulation_results = TOP_LEVEL_exefmuMultiobs(Background, VariablesToCalibrate=VariablesToCalibrate, OutputVariables=OutputVariables,  LNames = LNames, LColumns = LColumns, LVariablesToChange=LVariablesToChange, ref_simudir = self._get_Name_ModelTmpDir_simple(ModelName),  ModelName = ModelName, List_Multideltatime = List_Multideltatime, AdvancedDebugModel = AdvancedDebugModel, TimeoutModelExecution = TimeoutModelExecution, FMUInput = FMUInput)
1450
1451             elif self.__model["Format"].upper() in ["OPENMODELICA"]:
1452                 initialsimulation_results = TOP_LEVEL_exeOpenModelicaMultiobs(Background, KeepCalculationFolders = KeepCalculationFolders, VariablesToCalibrate=VariablesToCalibrate, OutputVariables=OutputVariables,  LNames = LNames, LColumns = LColumns, LVariablesToChange=LVariablesToChange, ref_simudir = self._get_Name_ModelTmpDir_simple(ModelName),  ModelName = ModelName, List_Multideltatime = List_Multideltatime, Linux = Linux, AdvancedDebugModel = AdvancedDebugModel,  TimeoutModelExecution = TimeoutModelExecution )
1453
1454             __resultats = dict(
1455                 InitialParameters = numpy.asarray(Background),
1456                 OptimalParameters = __adaocase.get("Analysis")[-1], # * Background,
1457                 InitialSimulation = numpy.asarray(initialsimulation_results), #the simulation results are given for backgroung values, can be used to evaluate the overall performance of the calibration procedure
1458                 OptimalSimulation = __adaocase.get("SimulatedObservationAtOptimum")[-1],
1459                 Measures          = Observations,
1460                 NamesOfParameters = BNames,
1461                 NamesOfMeasures   = ONames,
1462                 )
1463         else:
1464             __resultats = dict(
1465                 InitialParameters = numpy.asarray(Background),
1466                 OptimalParameters = __adaocase.get("Analysis")[-1], # * Background,
1467                 OptimalSimulation = __adaocase.get("SimulatedObservationAtOptimum")[-1],
1468                 Measures          = Observations,
1469                 NamesOfParameters = BNames,
1470                 NamesOfMeasures   = ONames,
1471                 )
1472
1473         if "APosterioriCovariance" in Params["StoreSupplementaryCalculations"]:
1474             __resultats["APosterioriCovariance"] = __adaocase.get("APosterioriCovariance")[-1]
1475
1476         if ResultsSummary:
1477             if InitialSimulation:
1478                 measures_input = __resultats ["Measures"]
1479                 initial_simu = __resultats ["InitialSimulation"]
1480                 optimal_simu = __resultats ["OptimalSimulation"]
1481
1482                 diff_measures_initial_rmse = numpy.array([(measures_input[i]-initial_simu[i])**2 for i in range(len(measures_input))])
1483                 diff_measures_initial_reldiff = numpy.array([abs((measures_input[i]-initial_simu[i])/measures_input[i]) for i in range(len(measures_input)) if measures_input[i] !=0])
1484
1485                 diff_measures_optimal_rmse = numpy.array([(measures_input[i]-optimal_simu[i])**2 for i in range(len(measures_input))])
1486                 diff_measures_optimal_reldiff = numpy.array([abs((measures_input[i]-optimal_simu[i])/measures_input[i]) for i in range(len(measures_input)) if measures_input[i] !=0])
1487
1488                 diff_measures_initial_rmse_res = str(numpy.format_float_scientific(numpy.sqrt(diff_measures_initial_rmse.mean()), precision = 3))
1489                 diff_measures_initial_reldiff_res = str(numpy.round(diff_measures_initial_reldiff.mean()*100,decimals =2))
1490
1491                 diff_measures_optimal_rmse_res = str(numpy.format_float_scientific(numpy.sqrt(diff_measures_optimal_rmse.mean()), precision = 3))
1492                 diff_measures_optimal_reldiff_res = str(numpy.round(diff_measures_optimal_reldiff.mean()*100,decimals =2))
1493
1494                 __resultats["ResultsSummary_RMSE"] = ["Average_RMSE_INITIAL = " + diff_measures_initial_rmse_res, "Average_RMSE_OPTIMAL = " + diff_measures_optimal_rmse_res ]
1495                 __resultats["ResultsSummary_RelativeDifference"] = ["Average_RelativeDifference_INITIAL = " +  diff_measures_initial_reldiff_res + "%", "Average_RelativeDifference_OPTIMAL = " + diff_measures_optimal_reldiff_res + "%" ]
1496             else:
1497                 measures_input = __resultats ["Measures"]
1498                 optimal_simu = __resultats ["OptimalSimulation"]
1499
1500                 diff_measures_optimal_rmse = numpy.array([(measures_input[i]-optimal_simu[i])**2 for i in range(len(measures_input))])
1501                 diff_measures_optimal_reldiff = numpy.array([abs((measures_input[i]-optimal_simu[i])/measures_input[i]) for i in range(len(measures_input)) if measures_input[i] !=0])
1502
1503
1504                 diff_measures_optimal_rmse_res = str(numpy.format_float_scientific(numpy.sqrt(diff_measures_optimal_rmse.mean()), precision = 3))
1505                 diff_measures_optimal_reldiff_res = str(numpy.round(diff_measures_optimal_reldiff.mean()*100,decimals =2))
1506
1507                 __resultats["ResultsSummary_RMSE"] = ["Average_RMSE_OPTIMAL = " + diff_measures_optimal_rmse_res ]
1508                 __resultats["ResultsSummary_RelativeDifference"] = ["Average_RelativeDifference_OPTIMAL = " + diff_measures_optimal_reldiff_res + "%" ]
1509
1510         #
1511         _saveResults(
1512             __resultats,
1513             self.__result["SourceName"],
1514             self.__result["Level"],
1515             self.__result["Format"],
1516             self.__verbose,
1517             )
1518         #
1519
1520         if CreateOptimaldsin:
1521              if self.__model["Format"].upper() in ["DYMOSIM"]:
1522                  for etat in range(len(LNames)):
1523                      dict_optimal_params = dict(zip(VariablesToCalibrate, __adaocase.get("Analysis")[-1])) #dans la boucle, sinon il se vide à chaque fois que l'on fait itère
1524                      if not(ComplexModel):
1525                          dir_for_dsin_opti = os.path.abspath(os.path.join(self._get_Name_ModelTmpDir_simple(ModelName),os.pardir))
1526                          simu_dir_opti = os.path.abspath(self._get_Name_ModelTmpDir_simple(ModelName))
1527
1528                          writer_opti = write_in_dsin.Write_in_dsin(dict_inputs =dict_optimal_params,filedir=simu_dir_opti,dsin_name=str("dsin_" +LNames[etat]+ ".txt"),old_file_name='dsin.txt',new_file_name= str("dsin_" +LNames[etat]+ "_optimal.txt"))
1529                          writer_opti.write_in_dsin
1530
1531                          writer_opti.write_in_dsin()
1532
1533                          try:
1534                              os.remove(os.path.join(dir_for_dsin_opti,str("dsin_" +LNames[etat]+ "_optimal.txt")))
1535                          except:
1536                              pass
1537
1538                          shutil.copyfile(
1539                                     os.path.join(simu_dir_opti, str("dsin_" +LNames[etat]+ "_optimal.txt")),
1540                                     os.path.join(dir_for_dsin_opti, str("dsin_" +LNames[etat]+ "_optimal.txt"))
1541                                     )
1542
1543                          auto_simul_opti = Automatic_Simulation(
1544                                                      simu_dir=simu_dir_opti,
1545                                                      dsin_name= str("dsin_" +LNames[etat]+ "_optimal.txt"),
1546                                                      dymosim_path = os.path.join(simu_dir_opti, os.pardir) ,
1547                                                      dsres_path = str("dsres_" +LNames[etat]+ "_opti.mat"),
1548                                                      dict_inputs= {},
1549                                                      logfile=False,
1550                                                      timeout=TimeoutModelExecution,
1551                                                      without_modelicares=True,
1552                                                      linux=Linux
1553                                                      )
1554                          auto_simul_opti.single_simulation()
1555
1556                          if not(auto_simul_opti.success_code == 0):
1557                              print("The dsin.txt with the optimal values does not converge directly, please update the value of iteration variables using Change_working_Point script (consider both calibrated paramaters and boundary conditions in the procedure)")
1558                      else:
1559                          try: #to handle situations in which there is only one CL (boundary conditions) file
1560                              var_int = numpy.ravel(LColumns[:,etat])
1561                          except :
1562                              var_int = numpy.ravel(LColumns)
1563                          dict_inputs_boundary_conditions_optimal_params = dict(zip(LVariablesToChange, var_int))
1564                          dict_inputs_boundary_conditions_optimal_params.update(dict_optimal_params)
1565
1566                          try:
1567                              os.remove(os.path.join(dir_for_dsin_opti,str("dsin_" +LNames[etat]+ "_optimal.txt")))
1568                          except:
1569                              pass
1570                          dir_for_dsin_opti = os.path.abspath(os.path.join(self._get_Name_ModelTmpDir_simple(ModelName),os.pardir))
1571                          shutil.copyfile(
1572                                     os.path.join(dir_for_dsin_opti, str("dsin.txt")),
1573                                     os.path.join(dir_for_dsin_opti, str("dsin_" +LNames[etat]+ "_optimal.txt"))
1574                                     )
1575                          writer_opti = write_in_dsin.Write_in_dsin(dict_inputs =dict_inputs_boundary_conditions_optimal_params,filedir=dir_for_dsin_opti,dsin_name=str("dsin_" +LNames[etat]+ "_optimal.txt"),old_file_name=str("dsin_" +LNames[etat]+ "_old.txt"),new_file_name= str("dsin_" +LNames[etat]+ "_optimal.txt"))
1576                          writer_opti.write_in_dsin
1577                          writer_opti.write_in_dsin()
1578                          os.remove(os.path.join(dir_for_dsin_opti,str("dsin_" +LNames[etat]+ "_old.txt")))
1579                          print("The dsin.txt with the optimal values has not been tested, updating the values of iteration variables might be necessary: this can be done by using Change_working_Point script (consider both calibrated paramaters and boundary conditions in the procedure)")
1580
1581              else: #FMI, OM or GUESS format
1582                  print("CreateOptimaldsin option is only valid for DYMOSIM format (it has no effect for the current calculation)")
1583
1584         if self.__verbose: print("")
1585         del __adaocase
1586
1587         if not(KeepCalculationFolders):
1588             shutil.rmtree(self._get_Name_ModelTmpDir_simple(ModelName), ignore_errors=True)
1589
1590         if CreateCSVoutputResults:
1591             print("Creation of CSV output results")
1592
1593
1594             optimal_results = __resultats["OptimalSimulation"]
1595             for obs_file_index in range(len(self.__measures["SourceName"])):
1596                 obs_filename = self.__measures["SourceName"][obs_file_index] #for the name of the outputresult
1597
1598                 time_measures_number = len(OHours[obs_file_index]) #get the number of lines in the obs file
1599                 total_measures_number_per_cl = time_measures_number*len(__resultats["NamesOfMeasures"]) #size of the list to be removed from the whole results file
1600                 optimal_results_cl = optimal_results[:total_measures_number_per_cl] # output results for the current cl
1601                 optimal_results = optimal_results[total_measures_number_per_cl:] # update of the whole results file (we remove te results already treated)
1602
1603                 optimal_results_cl_to_csv = optimal_results_cl.reshape((len(__resultats["NamesOfMeasures"]),-1)).transpose() #reshape to get one line per timedate
1604                 df_optimal_results_cl_to_csv = pandas.DataFrame(optimal_results_cl_to_csv, columns = __resultats["NamesOfMeasures"]) #creation of the df with the names of measures
1605                 df_optimal_results_cl_to_csv.insert(0,"Hour", OHours[obs_file_index]) #add Hour column
1606                 df_optimal_results_cl_to_csv.insert(0,"Date", ODates[obs_file_index]) # add Date column
1607
1608                 optimal_file_results_name_cl_tmp = "Optimal_simulation_tmp.csv"
1609                 df_optimal_results_cl_to_csv.to_csv(optimal_file_results_name_cl_tmp, sep = ";", index = False)
1610
1611                 optimal_file_results_name_cl = "Optimal_simulation_"+obs_filename
1612                 with open(optimal_file_results_name_cl_tmp, 'r') as infile:
1613                     readie=csv.reader(infile, delimiter=';')
1614                     with open(optimal_file_results_name_cl, 'wt', newline='') as output:
1615                         outwriter=csv.writer(output, delimiter=';')
1616                         outwriter.writerow(["#"])
1617                         for row in readie:
1618                              outwriter.writerow(row)
1619                 try:
1620                     os.remove(optimal_file_results_name_cl_tmp)
1621                 except:
1622                     pass
1623
1624             if InitialSimulation: #If initial simulation is required as well, the Initial results files are given as well
1625                 initial_results = __resultats["InitialSimulation"]
1626                 for obs_file_index in range(len(self.__measures["SourceName"])):
1627                     obs_filename = self.__measures["SourceName"][obs_file_index] #for the name of the outputresult
1628
1629                     time_measures_number = len(OHours[obs_file_index]) #get the number of lines in the obs file
1630                     total_measures_number_per_cl = time_measures_number*len(__resultats["NamesOfMeasures"]) #size of the list to be removed from the whole results file
1631                     initial_results_cl = initial_results[:total_measures_number_per_cl] # output results for the current cl
1632                     initial_results = initial_results[total_measures_number_per_cl:] # update of the whole results file (we remove te results already treated)
1633
1634                     initial_results_cl_to_csv = initial_results_cl.reshape((len(__resultats["NamesOfMeasures"]),-1)).transpose() #reshape to get one line per timedate
1635                     df_initial_results_cl_to_csv = pandas.DataFrame(initial_results_cl_to_csv, columns = __resultats["NamesOfMeasures"]) #creation of the df with the names of measures
1636                     df_initial_results_cl_to_csv.insert(0,"Hour", OHours[obs_file_index]) #add Hour column
1637                     df_initial_results_cl_to_csv.insert(0,"Date", ODates[obs_file_index]) # add Date column
1638
1639                     initial_file_results_name_cl_tmp = "Initial_simulation_tmp.csv"
1640                     df_initial_results_cl_to_csv.to_csv(initial_file_results_name_cl_tmp, sep = ";", index = False)
1641
1642                     initial_file_results_name_cl = "Initial_simulation_"+obs_filename
1643                     with open(initial_file_results_name_cl_tmp, 'r') as infile:
1644                         readie=csv.reader(infile, delimiter=';')
1645                         with open(initial_file_results_name_cl, 'wt', newline='') as output:
1646                             outwriter=csv.writer(output, delimiter=';')
1647                             outwriter.writerow(["#"])
1648                             for row in readie:
1649                                  outwriter.writerow(row)
1650                     try:
1651                         os.remove(initial_file_results_name_cl_tmp)
1652                     except:
1653                         pass
1654
1655         return __resultats
1656
1657
1658         def verify_gradient(self):
1659             raise NotImplementedError("Not yet implemented")
1660
1661 # ==============================================================================
1662 def _readLink(__filename=None, __colnames=None, __indexname="Variable", __format="Guess"):
1663     """
1664     Read file of link between model and measures
1665
1666     Arguments:
1667         - File name
1668         - Names of the columns to read, known by their variable header
1669         - Name of the column containing the variables to adapt
1670         - Format of the data (CSV, TSV... or can be guessed)
1671     """
1672     #
1673     if __colnames is None:
1674         raise
1675     if __indexname is None: __indexname="Variable"
1676
1677     try: #Solution provisoire pou gérer les cas sans CL précisées par l'utilisateur
1678         with ImportFromFile(__filename, __colnames, __indexname, __format, False) as reading:
1679             colnames, columns, indexname, variablestochange = reading.getvalue()
1680     except:
1681         colnames =__colnames
1682         columns =[]
1683         variablestochange =[]
1684         print("/!\\ Boundary conditions indicated in ", __filename, " are not considered for calculation. Indicate at least two boundary conditions so that they are considered. /!\\ ")
1685     #
1686     variablestochange = tuple([v.strip() for v in variablestochange])
1687     return (colnames, columns, variablestochange)
1688
1689 # ==============================================================================
1690 def _readMultiMeasures(__filenames=None, __colnames=None, __format="Guess"):
1691     """
1692     Read files of measures, using only some named variables by columns
1693
1694     Arguments:
1695         - File name
1696         - Names of the columns to read, known by their variable header
1697         - Format of the data (CSV, TSV... or can be guessed)
1698     """
1699     #
1700     MultiMeasures = []
1701     for __filename in __filenames:
1702         colnames, columns = _readMeasures(__filename, __colnames, __format)
1703         MultiMeasures.append( columns )
1704     #
1705     return (colnames, MultiMeasures)
1706
1707 # ==============================================================================
1708 def _readMeasures(__filename=None, __colnames=None, __format="Guess"):
1709     """
1710     Read files of measures, using only some named variables by columns
1711
1712     Arguments:
1713         - File name
1714         - Names of the columns to read, known by their variable header
1715         - Format of the data (CSV, TSV... or can be guessed)
1716     """
1717     #
1718     with ImportFromFile(__filename, __colnames, None, __format) as reading:
1719         colnames, columns, indexname, index = reading.getvalue()
1720     #
1721     return (colnames, columns)
1722
1723 # ==============================================================================
1724 def _readBackground(__filename="dsin.txt", __varnames=None, __backgroundformat="Guess"):
1725     """
1726     Read background in model definition. The priority is "USER", then "ADAO",
1727     then "DSIN", and the default one correspond to "DSIN".
1728
1729     Arguments:
1730         - File name
1731         - Names of the variables to be found in the model data
1732         - Format of the data (can be guessed)
1733     """
1734     __format = None
1735     if __backgroundformat.upper() in ["GUESS", "DSIN"]:
1736         __format = "DSIN"
1737         if __filename is not None and os.path.isfile(__filename):
1738             __model_dir = os.path.dirname(__filename)
1739             __model_fil = os.path.basename(__filename)
1740         elif  __filename is not None and os.path.isdir(__filename) and os.path.isfile(os.path.join(__filename, "dsin.txt")):
1741             __model_dir = os.path.abspath(__filename)
1742             __model_fil = "dsin.txt"
1743         else:
1744             raise ValueError('No such file or directory: %s'%str(__filename))
1745         __bgfile = os.path.abspath(os.path.join(__model_dir, __model_fil))
1746         if (len(__bgfile)<7) and (__bgfile[-8:].lower() != "dsin.txt" ):
1747             raise ValueError("Model definition file \"dsin.txt\" can not be found.")
1748         if __varnames is None:
1749             raise ValueError("Model variables to be read has to be given")
1750     if __backgroundformat.upper() == "ADAO":
1751         __format = "ADAO"
1752         if __filename is None or not os.path.isfile(__filename):
1753             raise ValueError('No such file or directory: %s'%str(__filename))
1754         __bgfile = os.path.abspath(__filename)
1755         if not( len(__bgfile)>4 ):
1756             raise ValueError("User file name \"%s\" is too short and seems not to be valid."%__bgfile)
1757     if __backgroundformat.upper() == "USER":
1758         __format = "USER"
1759         if __filename is None or not os.path.isfile(__filename):
1760             raise ValueError('No such file or directory: %s'%str(__filename))
1761         __bgfile = os.path.abspath(__filename)
1762         if not( len(__bgfile)>5 ):
1763             raise ValueError("User file name \"%s\" is too short and seems not to be valid."%__bgfile)
1764     if __format not in ["DSIN", "ADAO", "USER"]:
1765         raise ValueError("Background initial values asked format \"%s\" is not valid."%__format)
1766     #
1767     #---------------------------------------------
1768     if __format == "DSIN":
1769         if __varnames is None:
1770             raise ValueError("Names of variables to read has to be given")
1771         #
1772         __names, __background, __bounds = [], [], []
1773         with open(__bgfile, 'r') as fid:
1774             __content = fid.readlines()
1775             for v in tuple(__varnames):
1776                 if _verify_existence_of_variable_in_dsin(v, __content):
1777                     __value = _read_existing_variable_in_dsin(v, __content)
1778                     __names.append(v)
1779                     __background.append(__value[0])
1780                     __bounds.append(__value[1:3])
1781                     # print "%s = "%v, __value
1782         __names = tuple(__names)
1783         __background = numpy.ravel(__background)
1784         __bounds = tuple(__bounds)
1785     #---------------------------------------------
1786     elif __format.upper() == "ADAO":
1787         with open(__bgfile, 'r') as fid:
1788             exec(fid.read())
1789         #
1790         __background = locals().get("Background", None)
1791         if __background is None:
1792             raise ValueError("Background is not defined")
1793         else:
1794             __background = numpy.ravel(__background)
1795         __Params     = locals().get("Parameters", {})
1796         if "Bounds" not in __Params:
1797             raise ValueError("Bounds can not be find in file \"%s\""%str(__filename))
1798         else:
1799             __bounds = tuple(__Params["Bounds"])
1800             if len(__bounds) != len(__background):
1801                 raise ValueError("Initial background length does not match bounds number")
1802         #
1803         if __varnames is None:
1804             __names = ()
1805         else:
1806             __names = tuple(__varnames)
1807             if len(__names) != len(__background):
1808                 raise ValueError("Initial background length does not match names number")
1809     #---------------------------------------------
1810     elif __format.upper() == "USER":
1811         __names, __background, __bounds = ImportScalarLinesFromFile(__bgfile).getvalue(__varnames)
1812     #---------------------------------------------
1813     for i,b in enumerate(__bounds) :
1814         if b[0] is not None and b[1] is not None and not (b[0] < b[1]) :
1815             raise ValueError(f'Inconsistent boundaries values for "{__names[i]}": {b[0]} not < {b[1]}')
1816
1817     return (__names, __background, __bounds)
1818
1819 def _verify_existence_of_variable_in_dsin(__variable, __content, __number = 2):
1820     "Verify if the variable exist in the model file content"
1821     if "".join(__content).count(__variable) >= __number:
1822         return True
1823     else:
1824         return False
1825
1826 def _read_existing_variable_in_dsin(__variable, __content):
1827     "Return the value of the real variable"
1828     __previousline, __currentline = "", ""
1829     internalOrder = 0
1830     initialValuePart = False
1831     for line in __content:
1832         if not initialValuePart and "initialValue" not in line: continue
1833         if not initialValuePart and "initialValue" in line: initialValuePart = True
1834         __previousline = __currentline
1835         __currentline  = line
1836         if __variable in __currentline and "#" in __currentline and "#" not in __previousline:
1837             # Informations ecrites sur deux lignes successives
1838             vini, value, vmin, vmax = __previousline.split()
1839             # vcat, vtype, diez, vnam = __currentline.split()
1840             value = float(value)
1841             if float(vmin) >= float(vmax):
1842                 vmin, vmax = None, None
1843             else:
1844                 vmin, vmax = float(vmin), float(vmax)
1845             return (value, vmin, vmax)
1846         elif __variable in __currentline and "#" in __currentline and "#" in __previousline:
1847             # Informations ecrites sur une seule ligne
1848             vini, value, vmin, vmax, reste = __previousline.split(maxsplit=5)
1849             value = float(value)
1850             if float(vmin) >= float(vmax):
1851                 vmin, vmax = None, None
1852             else:
1853                 vmin, vmax = float(vmin), float(vmax)
1854             return (value, vmin, vmax)
1855     #
1856     raise ValueError("Value of variable %s not found in the file content"%__variable)
1857
1858 # ==============================================================================
1859 def _readMethod(__filename="parameters.py", __format="ADAO"):
1860     """
1861     Read data assimilation method parameters
1862
1863     Arguments:
1864         - File name
1865         - Format of the data (can be guessed)
1866     """
1867     if not os.path.isfile(__filename):
1868         raise ValueError('No such file or directory: %s'%str(__filename))
1869     #
1870     #---------------------------------------------
1871     if __format.upper() in ["GUESS", "ADAO"]:
1872         with open(__filename, 'r') as fid:
1873             exec(fid.read())
1874         __Algo   = locals().get("Algorithm", "3DVAR")
1875         __Params = locals().get("Parameters", {})
1876         __CovB   = locals().get("BackgroundError", 1.)
1877         __CovR   = locals().get("ObservationError", 1.)
1878         # __Xb     = locals().get("Background", __background)
1879         # if not isinstance(__Params, dict): __Params = {"Bounds":__bounds}
1880         # if "Bounds" not in __Params:       __Params["Bounds"] = __bounds
1881     #---------------------------------------------
1882     else:
1883         __Algo, __Params, __CovB, __CovR = "3DVAR", {}, 1., 1.
1884     #---------------------------------------------
1885     # if __Xb is None:
1886     #     raise ValueError("Background is not defined")
1887     # else:
1888     #     __Xb = numpy.ravel(__Xb)
1889     #
1890     return (__Algo, __Params, __CovB, __CovR)
1891
1892 # ==============================================================================
1893 def _saveResults(__resultats, __filename=None, __level=None, __format="Guess", __verbose=False):
1894     """
1895     Save results relatively to the level required
1896
1897     Arguments!
1898         - File name
1899         - Level of saving : final output only or intermediary with final output
1900         - Format of the data
1901     """
1902     if __filename is None:
1903         raise ValueError('A file to save results has to be named, please specify one')
1904     if __format.upper() == "GUESS":
1905         if __filename.split(".")[-1].lower() == "txt":
1906             __format = "TXT"
1907         elif __filename.split(".")[-1].lower() == "py":
1908             __format = "PY"
1909         else:
1910             raise ValueError("Can not guess the file format of \"%s\", please specify the good one"%__filename)
1911     else:
1912         __format     = str(__format).upper()
1913     if __format not in ["TXT", "PY"]:
1914         raise ValueError("Result file format \"%s\" is not valid."%__format)
1915     if __format == "PY" and os.path.splitext(__filename)[1] != ".py":
1916         __filename = os.path.splitext(__filename)[0] + ".py"
1917     #
1918     #---------------------------------------------
1919     if __format.upper() == "TXT":
1920         output = ["# Final results",]
1921         keys = list(__resultats.keys())
1922         keys.sort()
1923         for k in keys:
1924             output.append("%18s = %s"%(k,tuple(__resultats[k])))
1925         output.append("")
1926         with open(__filename, 'w') as fid:
1927             fid.write( "\n".join(output) )
1928             fid.flush()
1929     #---------------------------------------------
1930     elif __format.upper() == "PY":
1931         output = ["# Final results",]
1932         keys = list(__resultats.keys())
1933         keys.sort()
1934         for k in keys:
1935             output.append("%s = %s"%(k,tuple(__resultats[k])))
1936         output.append("")
1937         with open(__filename, 'w') as fid:
1938             fid.write( "\n".join(output) )
1939             fid.flush()
1940     #---------------------------------------------
1941     elif __format.upper() == "CSV":
1942         raise NotImplementedError("Not yet implemented")
1943     #---------------------------------------------
1944     elif __format.upper() == "DICT":
1945         raise NotImplementedError("Not yet implemented")
1946     #---------------------------------------------
1947     if __verbose: # Format TXT
1948         output = ["# Final results",]
1949         keys = list(__resultats.keys())
1950         keys.sort()
1951         for k in keys:
1952             output.append("%18s = %s"%(k,tuple(__resultats[k])))
1953         output.append("")
1954         print( "\n".join(output) )
1955     #---------------------------------------------
1956     return True
1957
1958 def _secure_copy_textfile(__original_file, __destination_file):
1959     "Copy the file guranteeing that it is copied"
1960
1961
1962     shutil.copy(
1963                 __original_file,
1964                 __destination_file
1965                 )
1966     shutil.copystat( # commande pour assurer que le fichier soit bien copié
1967         __original_file,
1968         __destination_file
1969         )
1970     os.path.getsize(__destination_file)# commande bis pour assurer que le fichier soit bien copié
1971
1972     def file_len(fname):
1973         count = 0
1974         with open(fname) as f:
1975             for line in f:
1976                 count += 1
1977         return count
1978
1979     while file_len(__original_file)!= file_len(__destination_file):
1980         shutil.copy(
1981             __original_file,
1982             __destination_file
1983             )
1984
1985 # ==============================================================================
1986 class Python_Simulation(object):
1987     def __init__(self, simu_dir=".", val_inputs=(), logfile=True, timeout=60, verbose=True):
1988         self.__simu_dir = os.path.abspath(simu_dir)
1989         self.__inputs = numpy.ravel(val_inputs)
1990         if verbose:
1991             print()
1992             print("  [VERBOSE] Input values %s"%self.__inputs)
1993
1994     def single_simulation(self, __inputs = None):
1995         with open(os.path.join(self.__simu_dir,"pythonsim.exe"), 'r') as fid:
1996             exec(fid.read())
1997         __directoperator = locals().get("DirectOperator", None)
1998         if __inputs is None: __inputs = self.__inputs
1999         return __directoperator(__inputs)
2000
2001 # ==============================================================================
2002 class VersionsLogicielles (object):
2003     "Modules version information"
2004     def __init__(self):
2005         print("  [VERBOSE] System configuration:")
2006         print("      - Python...:",sys.version.split()[0])
2007         print("      - Numpy....:",numpy.version.version)
2008         print("      - Scipy....:",scipy.version.version)
2009         print("      - ADAO.....:",adao.version)
2010         print("")
2011
2012 def TOP_LEVEL_exefmuMultiobs_ref( x_values_matrix , VariablesToCalibrate=None, OutputVariables=None, ref_simudir = None,  ModelName = None ):
2013     "Appel du modèle en format FMU et restitution des résultats"
2014     x_values = list(numpy.ravel(x_values_matrix))
2015     x_inputs=dict(zip(VariablesToCalibrate,x_values))
2016
2017     fmuname = os.path.basename(ModelName) #in case ModelName is given as a path (works also if it is a file name)
2018     fmu = os.path.join(ref_simudir,fmuname)
2019
2020     reader = simulate_fmu(fmu, output = OutputVariables, start_values = x_inputs)
2021     y_values = [reader[y_name][-1] for y_name in OutputVariables]
2022     y_values_matrix = numpy.asarray(y_values)
2023
2024     return y_values_matrix
2025
2026 def TOP_LEVEL_exefmuMultiobs( x_values_matrix , VariablesToCalibrate=None, OutputVariables=None,  LNames = None, LColumns = None, LVariablesToChange=None, ref_simudir = None,  ModelName = None, List_Multideltatime = None, AdvancedDebugModel = None, TimeoutModelExecution = None, FMUInput = None):
2027     "Appel du modèle en format FMU et restitution des résultats"
2028
2029     if VariablesToCalibrate is None:
2030         raise ValueError("VariablesToCalibrate")
2031
2032     if OutputVariables is None:
2033         raise ValueError("OutputVariables")
2034
2035     if LNames is None:
2036         raise ValueError("LNames")
2037
2038     if LColumns is None:
2039         raise ValueError("LColumns")
2040
2041     if LVariablesToChange is None:
2042         raise ValueError("LVariablesToChange")
2043
2044     if ref_simudir is None:
2045         raise ValueError("ref_simudir")
2046
2047     if List_Multideltatime is None:
2048         raise ValueError("Problem defining simulation output results")
2049
2050     if AdvancedDebugModel is None:
2051         raise ValueError("AdvancedDebugModel")
2052
2053     if TimeoutModelExecution is None:
2054         raise ValueError("TimeoutModelExecution")
2055
2056     x_values = list(numpy.ravel(x_values_matrix))
2057     x_inputs=dict(zip(VariablesToCalibrate,x_values))
2058
2059     if os.path.isfile(ModelName):
2060         fmuname = os.path.basename(ModelName)
2061     elif os.path.isdir(ModelName):
2062         fmuname = [files for files in os.listdir(os.path.abspath(ModelName)) if files[-4:] =='.fmu'][0]
2063
2064     fmu = os.path.join(ref_simudir,fmuname)
2065
2066     multi_y_values_matrix = []
2067
2068     for etat in range(len(LNames)):
2069         y_values =[]
2070
2071         try: #to handle situations in which there is only one CL (boundary conditions) file
2072             var_int = numpy.ravel(LColumns[:,etat])
2073         except :
2074             var_int = numpy.ravel(LColumns)
2075
2076         dict_inputs = dict(zip(LVariablesToChange, var_int))
2077         #
2078         dict_inputs.update(x_inputs)
2079
2080         if FMUInput:
2081             fmu_inputs = FMUInput[LNames[etat]]
2082             timestep = fmu_inputs[1][0] - fmu_inputs[0][0]  #Assuming constant timestep
2083             if AdvancedDebugModel: print(f'The timestep for {LNames[etat]} is {timestep} seconds')
2084         else:
2085             fmu_inputs = None
2086             timestep = None
2087
2088         try:
2089             stoptime_fmu = List_Multideltatime[etat][-1]
2090
2091         except IndexError:
2092             stoptime_fmu = None
2093
2094
2095         if AdvancedDebugModel == True:
2096             old_stdout = sys.stdout
2097             new_stdout = io.StringIO()
2098             sys.stdout = new_stdout
2099
2100             start_time_simulation_fmi = time.time()#timeout manangement since fmpy does not raise an error for this, it just ends the simulations and continues
2101
2102             try :
2103                 reader = simulate_fmu(fmu, output = OutputVariables, start_values = dict_inputs, debug_logging = True, timeout = TimeoutModelExecution, input = fmu_inputs, stop_time= stoptime_fmu, output_interval= timestep)
2104             except FMICallException as e:
2105                 output = new_stdout.getvalue()
2106                 sys.stdout = old_stdout
2107                 print('[ERROR] Failed simulation with the following input:\n',dict_inputs)
2108                 raise e
2109
2110             elapsed_time_simulation_fmi =  time.time() - start_time_simulation_fmi
2111             if elapsed_time_simulation_fmi > TimeoutModelExecution:
2112                 raise TimeoutError("Timeout for simulation reached, please increase it in order to be able to simulate your model and/or check if your model is correct (use TimeoutModelExecution option in configuration.py file)")
2113
2114             output = new_stdout.getvalue()
2115             sys.stdout = old_stdout
2116
2117             dir_run=ref_simudir
2118             log_file=os.path.join(dir_run,
2119                                   'log_debug_model.txt')
2120             try: #try toremove the previous file
2121                 os.remove(log_file)
2122             except:
2123                 pass
2124
2125             f=open(log_file,'a')
2126             for line in output:
2127                 f.write(line)
2128
2129         else:
2130             start_time_simulation_fmi = time.time()
2131
2132             try :
2133                 reader = simulate_fmu(fmu, output = OutputVariables, start_values = dict_inputs, timeout = TimeoutModelExecution, input = fmu_inputs, stop_time= stoptime_fmu, output_interval= timestep)
2134             except FMICallException as e:
2135                 print('[ERROR] Failed simulation with the following input:\n',dict_inputs)
2136                 raise e
2137
2138             elapsed_time_simulation_fmi =  time.time() - start_time_simulation_fmi
2139             if elapsed_time_simulation_fmi > TimeoutModelExecution:
2140                 raise TimeoutError("Timeout for simulation reached, please increase it in order to be able to simulate your model and/or check if your model is correct (use TimeoutModelExecution option in configuration.py file)")
2141
2142         y_whole = [reader[y_name] for y_name in OutputVariables]
2143
2144         for y_ind in range(len(OutputVariables)):
2145             y_ind_whole_values = y_whole[y_ind]
2146             y_ind_whole_time = reader['time'] #standard output of fmu
2147
2148             if len(List_Multideltatime[etat])>1:
2149                 index_y_values = [find_nearest_index(y_ind_whole_time,time)[0] for time in List_Multideltatime[etat]]
2150             else:
2151                 index_y_values = [-1] #we only take the last point if one measure point
2152             y_ind_values = [y_ind_whole_values[i] for i in index_y_values]
2153             y_values = y_values + y_ind_values
2154
2155         y_values_matrix = y_values #pbs in results management
2156         multi_y_values_matrix = multi_y_values_matrix + y_values_matrix
2157
2158     y_values_matrix_def = numpy.ravel(numpy.array(multi_y_values_matrix))
2159     return y_values_matrix_def
2160
2161 def TOP_LEVEL_exeOpenModelicaMultiobs( x_values_matrix , KeepCalculationFolders = None, VariablesToCalibrate=None, OutputVariables=None,  LNames = None, LColumns = None, LVariablesToChange=None, ref_simudir = None,  ModelName = None, List_Multideltatime = None, Linux = None, AdvancedDebugModel = None, TimeoutModelExecution = None):
2162     "Appel du modèle en format OpenModelica et restitution des résultats"
2163
2164     if VariablesToCalibrate is None:
2165         raise ValueError("VariablesToCalibrate")
2166
2167     if OutputVariables is None:
2168         raise ValueError("OutputVariables")
2169
2170     if LNames is None:
2171         raise ValueError("LNames")
2172
2173     if LColumns is None:
2174         raise ValueError("LColumns")
2175
2176     if LVariablesToChange is None:
2177         raise ValueError("LVariablesToChange")
2178
2179     if ref_simudir is None:
2180         raise ValueError("ref_simudir")
2181
2182     if KeepCalculationFolders is None:
2183         raise ValueError("KeepCalculationFolders")
2184
2185     if Linux is None:
2186         raise ValueError("Linux")
2187
2188     if List_Multideltatime is None:
2189         raise ValueError("Problem defining simulation output results")
2190
2191     if AdvancedDebugModel is None:
2192         raise ValueError("AdvancedDebugModel")
2193
2194     if TimeoutModelExecution is None:
2195         raise ValueError("TimeoutModelExecution")
2196
2197     x_values = list(numpy.ravel(x_values_matrix))
2198     x_inputs=dict(zip(VariablesToCalibrate,x_values))
2199
2200     if os.path.isfile(ModelName):
2201         om_xml = os.path.basename(ModelName)
2202     elif os.path.isdir(ModelName):
2203         om_xml = [files for files in os.listdir(os.path.abspath(ModelName)) if files[-4:] =='.xml'][0]
2204
2205     om_xml = os.path.abspath(os.path.join(ref_simudir,om_xml))
2206
2207     multi_y_values_matrix = []
2208
2209     for etat in range(len(LNames)):
2210         y_values =[]
2211
2212         try: #to handle situations in which there is only one CL (boundary conditions) file
2213             var_int = numpy.ravel(LColumns[:,etat])
2214         except :
2215             var_int = numpy.ravel(LColumns)
2216
2217         dict_inputs = dict(zip(LVariablesToChange, var_int))
2218         #
2219
2220         dict_inputs.update(x_inputs)
2221
2222
2223         results_dir = tempfile.mkdtemp( dir=os.path.dirname(om_xml) ) #to handle parallel computation
2224
2225         results_file_name = os.path.join(results_dir,os.path.basename(om_xml)[:-9] + '_' + LNames[etat].replace('.','_') + ".mat")
2226
2227         OM_execution = run_OM_model(om_xml, dict_inputs = dict_inputs, result_file_name = results_file_name , Linux = Linux, AdvancedDebugModel = AdvancedDebugModel, TimeoutModelExecution = TimeoutModelExecution)
2228
2229         try:
2230             reader = Reader(results_file_name,'dymola') #dymola even if it is OpenModelica
2231         except:
2232             raise ValueError("Simulation cannot be performed: reduce the number of parameters to calibrate and/or the range in which their optimal value should be found (or modify and simplify the model to make it easier to simulate)" )
2233
2234         y_whole = [reader.values(y_name) for y_name in OutputVariables]
2235
2236         for y_ind in range(len(OutputVariables)):
2237             y_ind_whole_values = y_whole[y_ind][1]
2238             y_ind_whole_time = y_whole[y_ind][0]
2239
2240             if len(List_Multideltatime[etat])>1:
2241                 index_y_values = [find_nearest_index(y_ind_whole_time,time)[0] for time in List_Multideltatime[etat]]
2242             else:
2243                 index_y_values = [-1] #we only take the last point if one measure point
2244
2245             y_ind_values = [y_ind_whole_values[i] for i in index_y_values]
2246             y_values = y_values + y_ind_values
2247 #        y_values_matrix = numpy.asarray(y_values)
2248         y_values_matrix = y_values #pbs in results management
2249         multi_y_values_matrix = multi_y_values_matrix + y_values_matrix
2250
2251         if not KeepCalculationFolders:
2252             shutil.rmtree(results_dir, ignore_errors=True)
2253
2254     y_values_matrix_def = numpy.ravel(numpy.array(multi_y_values_matrix))
2255
2256     return y_values_matrix_def
2257
2258 def TOPLEVEL_exedymosimMultiobs_simple( x_values_matrix , VariablesToCalibrate=None, OutputVariables=None, LNames = None, ref_simudir = None, KeepCalculationFolders = None, Linux = None, List_Multideltatime = None, AdvancedDebugModel = None, TimeoutModelExecution = None):
2259     "Appel du modèle en format DYMOSIM et restitution des résultats"
2260
2261
2262     if VariablesToCalibrate is None:
2263         raise ValueError("VariablesToCalibrate")
2264
2265     if OutputVariables is None:
2266         raise ValueError("OutputVariables")
2267
2268     if LNames is None:
2269         raise ValueError("LNames")
2270
2271     if ref_simudir is None:
2272         raise ValueError("ref_simudir")
2273
2274     if KeepCalculationFolders is None:
2275         raise ValueError("KeepCalculationFolders")
2276
2277     if Linux is None:
2278         raise ValueError("Linux")
2279
2280     if List_Multideltatime is None:
2281         raise ValueError("Problem defining simulation output results")
2282
2283     if AdvancedDebugModel is None:
2284         raise ValueError("AdvancedDebugModel")
2285
2286     if TimeoutModelExecution is None:
2287         raise ValueError("TimeoutModelExecution")
2288     #
2289     # x_values = list(numpy.ravel(x_values_matrix) * Background)
2290     x_values = list(numpy.ravel(x_values_matrix))
2291     x_inputs=dict(zip(VariablesToCalibrate,x_values))
2292     # if Verbose: print("  Simulation for %s"%numpy.ravel(x_values))
2293     #
2294
2295     multi_y_values_matrix = []
2296
2297     for etat in range(len(LNames)):
2298         y_values =[]
2299
2300         simudir = tempfile.mkdtemp( dir=ref_simudir )
2301         #
2302         dict_inputs ={}
2303         dict_inputs.update( x_inputs )
2304         #
2305         shutil.copy(
2306             os.path.join(ref_simudir, str("dsin_" +LNames[etat]+ ".txt")),
2307             os.path.join(simudir, str("dsin_" +LNames[etat]+ ".txt"))
2308             )
2309         _secure_copy_textfile(
2310                         os.path.join(ref_simudir, str("dsin_" +LNames[etat]+ ".txt")),
2311                         os.path.join(simudir, str("dsin_" +LNames[etat]+ ".txt"))
2312                         )
2313
2314         auto_simul = Automatic_Simulation(
2315                                         simu_dir = simudir,
2316                                         dymosim_path = os.path.join(ref_simudir,os.pardir) ,
2317                                         dsin_name = str("dsin_" +LNames[etat]+ ".txt"),
2318                                         dsres_path = os.path.join(simudir, str("dsres" +LNames[etat]+ ".mat")),
2319                                         dict_inputs = dict_inputs,
2320                                         logfile = True,
2321                                         timeout = TimeoutModelExecution,
2322                                         without_modelicares = True,
2323                                         linux=Linux,
2324                                         )
2325         auto_simul.single_simulation()
2326
2327         if AdvancedDebugModel:
2328             dslog_file = os.path.join(simudir, 'dslog.txt')
2329             debug_model_file = os.path.join(ref_simudir,'log_debug_model.txt')
2330             try:
2331                 shutil.copy(dslog_file,debug_model_file)
2332             except:
2333                 pass
2334
2335         if auto_simul.success_code == 2 :
2336             raise ValueError("The simulation falis after initialization, please check the log file and your model in order to be sure that the whole simulation can be performed (specially for dynamic models)")
2337
2338         if auto_simul.success_code == 3 :
2339             raise ValueError("The simulation did not reach the final time, probably due to a lack of time, please increase the time for model execution by modifying the TimeoutModelExecution value in configuration.py")
2340
2341         if auto_simul.success_code == 0:
2342             pass #means that everything OK, we keep going
2343         else  :
2344             raise ValueError("The NeedInitVariables iteration should be set to True, convergence issues found for %s" % str(LNames[etat])) #means this value is 1--> not correct initialization
2345
2346         reader = Reader(os.path.join(simudir, str("dsres" +LNames[etat]+ ".mat")),'dymola')
2347
2348         y_whole = [reader.values(y_name) for y_name in OutputVariables]
2349
2350         for y_ind in range(len(OutputVariables)):
2351             y_ind_whole_values = y_whole[y_ind][1]
2352             y_ind_whole_time = y_whole[y_ind][0]
2353
2354             if len(List_Multideltatime[etat])>1:
2355                 index_y_values = [find_nearest_index(y_ind_whole_time,time)[0] for time in List_Multideltatime[etat]]
2356             else:
2357                 index_y_values = [-1] #we only take the last point if one measure point
2358
2359             y_ind_values = [y_ind_whole_values[i] for i in index_y_values]
2360             y_values = y_values + y_ind_values
2361 #        y_values_matrix = numpy.asarray(y_values)
2362         y_values_matrix = y_values #pbs in results management
2363         multi_y_values_matrix = multi_y_values_matrix + y_values_matrix
2364
2365         if not KeepCalculationFolders:
2366             shutil.rmtree(simudir, ignore_errors=True)
2367
2368     y_values_matrix_def = numpy.ravel(numpy.array(multi_y_values_matrix))
2369     return y_values_matrix_def
2370
2371 def TOP_LEVEL_exedymosimMultiobs( x_values_matrix , VariablesToCalibrate=None, OutputVariables=None, LNames = None, ModelFormat = None, KeepCalculationFolders = None, Verbose = None, dict_dsin_paths = None, Linux = None, List_Multideltatime = None, AdvancedDebugModel = None, TimeoutModelExecution = None): #2ème argument à ne pas modifier
2372     "Appel du modèle en format DYMOSIM et restitution des résultats"
2373
2374     if VariablesToCalibrate is None:
2375         raise ValueError("VariablesToCalibrate")
2376
2377     if OutputVariables is None:
2378         raise ValueError("OutputVariables")
2379
2380     if LNames is None:
2381         raise ValueError("LNames")
2382
2383     if ModelFormat is None:
2384         raise ValueError("ModelFormat")
2385
2386     if KeepCalculationFolders is None:
2387         raise ValueError("KeepCalculationFolders")
2388
2389     if Verbose is None:
2390         raise ValueError("Verbose")
2391
2392     if dict_dsin_paths is None:
2393         raise ValueError("dict_dsin_paths")
2394
2395     if Linux is None:
2396         raise ValueError("Linux")
2397
2398     if List_Multideltatime is None:
2399         raise ValueError("Problem defining simulation output results")
2400
2401     if AdvancedDebugModel is None:
2402         raise ValueError("AdvancedDebugModel")
2403
2404     if TimeoutModelExecution is None:
2405         raise ValueError("TimeoutModelExecution")
2406 #
2407     # x_values = list(numpy.ravel(x_values_matrix) * Background)
2408     x_values = list(numpy.ravel(x_values_matrix))
2409     x_inputs = dict(zip(VariablesToCalibrate,x_values))
2410     x_inputs_for_CWP = {}
2411     # if Verbose: print("  Simulation for %s"%numpy.ravel(x_values))
2412
2413     multi_y_values_matrix = []
2414
2415
2416
2417     for etat in range(len(LNames)):
2418         y_values =[]
2419
2420         ref_simudir = Calibration._setModelTmpDir( None,dict_dsin_paths[LNames[etat]], ModelFormat,   'dsin.txt', os.path.join(os.path.pardir,os.path.pardir),  LNames[etat])
2421         simudir = ref_simudir
2422         dict_inputs={}
2423         dict_inputs.update(x_inputs)
2424
2425
2426         auto_simul = Automatic_Simulation(
2427             simu_dir=simudir,
2428             dymosim_path = os.path.join(simudir, os.pardir, os.pardir) ,
2429             dsres_path = os.path.join(simudir, str("dsres" +LNames[etat]+ ".mat")),
2430             dict_inputs=dict_inputs,
2431             logfile=True,
2432             timeout=TimeoutModelExecution,
2433             without_modelicares=True,
2434             linux=Linux
2435             )
2436         auto_simul.single_simulation()
2437
2438         if AdvancedDebugModel:
2439             dslog_file = os.path.join(simudir, 'dslog.txt')
2440             debug_model_file = os.path.join(simudir, os.pardir, 'log_debug_model.txt')
2441             try:
2442                 shutil.copy(dslog_file,debug_model_file)
2443             except:
2444                 pass
2445
2446         if auto_simul.success_code == 2 :
2447             raise ValueError("The simulation falis after initialization, please check the log file and your model in order to be sure that the whole simulation can be performed (specially for dynamic models)")
2448
2449         if auto_simul.success_code == 3 :
2450             raise ValueError("The simulation did not reach the final time, probably due to a lack of time, please increase the time for model execution by modifying the TimeoutModelExecution value in configuration.py")
2451
2452
2453         if auto_simul.success_code == 0:
2454             reader = Reader(os.path.join(simudir, str("dsres" +LNames[etat]+ ".mat")),'dymola')
2455
2456         else:
2457
2458             path_around_simu = os.path.join(simudir,os.path.pardir,os.path.pardir)
2459             around_simu = Around_Simulation(dymola_version='2012',
2460                               curdir=path_around_simu,
2461                               source_list_iter_var = 'from_dymtranslog',
2462                               copy_from_dym_trans_log= os.path.join(path_around_simu,'ini.txt'),
2463                               mat = os.path.join(os.path.dirname(dict_dsin_paths[LNames[etat]]), str("REF_for_CWP" + ".mat")),
2464                               iter_var_values_options={'source':'from_mat','moment':'initial'},
2465                               verbose = False)
2466             around_simu.set_list_iter_var()
2467             around_simu.set_dict_iter_var()
2468
2469             reader_CWP = Reader(os.path.join(os.path.dirname(dict_dsin_paths[LNames[etat]]), str("REF_for_CWP" + ".mat")),'dymola')
2470             x_values_intemediary = [reader_CWP.values(x_name)[1][-1] for x_name in VariablesToCalibrate]
2471             x_values_from_last_calculation = numpy.asarray(x_values_intemediary)
2472             x_inputs_from_last_calculation = dict(zip(VariablesToCalibrate,x_values_from_last_calculation))
2473
2474
2475             for var_calib in x_inputs_from_last_calculation.keys():
2476                 x_inputs_for_CWP[var_calib] = (x_inputs_from_last_calculation[var_calib], x_inputs[var_calib])
2477
2478             dict_var_to_fix2 = Dict_Var_To_Fix(option='automatic',
2479                                   dict_auto_var_to_fix=x_inputs_for_CWP)
2480
2481             dict_var_to_fix2.set_dict_var_to_fix()
2482
2483             if Verbose:
2484                 LOG_FILENAME = os.path.join(simudir,'change_working_point.log')
2485                 for handler in logging.root.handlers[:]:
2486                     logging.root.removeHandler(handler)
2487                 logging.basicConfig(filename=LOG_FILENAME,level=logging.DEBUG,filemode='w')
2488
2489             working_point_modif = Working_Point_Modification(main_dir = simudir,
2490                          simu_material_dir = simudir,
2491                          dymosim_path = os.path.join(simudir, os.pardir, os.pardir),
2492                          simu_dir = 'SIMUS',
2493                          store_res_dir = 'RES',
2494                          dict_var_to_fix = dict_var_to_fix2.dict_var_to_fix,
2495                          around_simulation0 = around_simu,
2496                          var_to_follow_path= os.path.join(simudir,'var_to_follow.csv'),
2497                          gen_scripts_ini= False,
2498                          nit_max= 1000000000000000,
2499                          min_step_val = 0.00000000000005,
2500                          timeout = TimeoutModelExecution,
2501                          linux=Linux)
2502
2503             working_point_modif.create_working_directory()
2504             working_point_modif.working_point_modification(skip_reference_simulation=True)
2505
2506             if AdvancedDebugModel:
2507                 dslog_file = os.path.join(simudir,'SIMUS', 'dslog.txt')
2508                 debug_model_file = os.path.join(simudir, os.pardir, 'log_debug_model.txt')
2509                 try:
2510                     shutil.copy(dslog_file,debug_model_file)
2511                 except:
2512                     pass
2513             try:
2514                 reader = Reader(os.path.join(simudir, 'RES','1.0.mat'),'dymola')
2515             except:
2516                 raise ValueError("Simulation cannot be performed: reduce the number of parameters to calibrate and/or the range in which their optimal value should be found (or modify and simplify the model to make it easier to simulate)" )
2517
2518
2519         y_whole = [reader.values(y_name) for y_name in OutputVariables]
2520
2521         for y_ind in range(len(OutputVariables)):
2522             y_ind_whole_values = y_whole[y_ind][1]
2523             y_ind_whole_time = y_whole[y_ind][0]
2524             if len(List_Multideltatime[etat])>1:
2525                 index_y_values = [find_nearest_index(y_ind_whole_time,time)[0] for time in List_Multideltatime[etat]]
2526             else:
2527                 index_y_values = [-1] #we only take the last point if one measure point
2528
2529             y_ind_values = [y_ind_whole_values[i] for i in index_y_values]
2530             y_values = y_values + y_ind_values
2531 #        y_values_matrix = numpy.asarray(y_values)
2532         y_values_matrix = y_values #pbs in results management
2533         multi_y_values_matrix = multi_y_values_matrix + y_values_matrix
2534
2535         if not KeepCalculationFolders:
2536             shutil.rmtree(simudir, ignore_errors=True)
2537     y_values_matrix_def = numpy.ravel(numpy.array(multi_y_values_matrix))
2538     return y_values_matrix_def
2539
2540 def run_OM_model(xml_file, #req arg, file path to _init.xml file
2541                  result_file_name = None, #file path for storing results
2542                  dict_inputs = None, #optionnal argument : for overriding paramters
2543                  Linux = False,
2544                  AdvancedDebugModel = None,
2545                  TimeoutModelExecution = None):
2546
2547     Command=[] #initialize command with list of argument
2548
2549     xml_file = os.path.abspath(xml_file)
2550
2551     #set main command to be runned
2552     if Linux:
2553         OM_binary = xml_file.replace('_init.xml','') # not tested
2554     else: #WIndows
2555         OM_binary = xml_file.replace('_init.xml','.exe')
2556     Command.append(OM_binary)
2557
2558     #Get base path for binaries and input
2559
2560
2561     inputPath='-inputPath='+os.path.dirname(xml_file)
2562     Command.append(inputPath)
2563
2564     #Generate override command
2565     if dict_inputs !=None:
2566         Override='-override='
2567         for elt in dict_inputs.items():
2568             Override=Override+elt[0]+'='+str(elt[1])+','
2569         Override=Override.rstrip(',') #remove last ',' if any
2570         Command.append(Override)
2571
2572     #Generate name of result file
2573     if result_file_name != None:
2574         result='-r='+result_file_name
2575         Command.append(result)
2576     else:
2577         result='-r='+OM_binary+'.mat' #by default result is stored next to _init and bin file
2578         Command.append(result)
2579
2580     result='-w' #by default result is stored next to _init and bin file
2581     Command.append(result)
2582
2583     result='-lv='+'LOG_STATS,LOG_NLS_V' #by default result is stored next to _init and bin file
2584     Command.append(result)
2585
2586     inputPath='-outputPath='+os.path.dirname(xml_file)
2587     Command.append(inputPath)
2588
2589    #launch calculation
2590     proc = subprocess.Popen(Command, #Command in the form of a text list
2591                             stderr=subprocess.STDOUT,
2592                             stdout=subprocess.PIPE,
2593                             universal_newlines=True)
2594     # for line in proc.stdout.readlines():
2595     #     print(line)
2596
2597     if AdvancedDebugModel == True:
2598         dir_run=os.path.join(os.path.dirname(result_file_name),os.pardir)
2599         log_file=os.path.join(dir_run,
2600                               'log_debug_model.txt')
2601
2602         try:
2603             os.remove(log_file)
2604         except:
2605             pass
2606         try:
2607             f=open(log_file,'a')
2608             for line in proc.stdout.readlines():
2609                 f.write(line)
2610             f.close()
2611             proc.communicate(timeout = TimeoutModelExecution)
2612         except subprocess.TimeoutExpired:
2613             raise ValueError("Timeout for simulation reached, please increase it in order to be able to simulate your model and/or check if your model is correct (use TimeoutModelExecution option in configuration.py file)")
2614
2615     else:
2616         try:
2617             proc.communicate(timeout = TimeoutModelExecution)
2618         except:
2619             raise ValueError("Timeout for simulation reached, please increase it in order to be able to simulate your model and/or check if your model is correct (use TimeoutModelExecution option in configuration.py file)")
2620
2621     # proc.communicate()
2622
2623     # r = proc.poll()
2624     # while r !=0:
2625     #     r = proc.poll()
2626
2627     # if proc.returncode !=0:
2628     #     raise ValueError("Simulation not ended")
2629
2630
2631     # sys.stdout.flush()
2632     # sys.stderr.flush()
2633
2634     # dir_run=os.path.dirname(result_file_name)
2635     # log_file=os.path.join(dir_run,
2636     #                       'log.txt')
2637     # try:
2638     #     proc.wait(timeout=50)
2639
2640     #     f=open(log_file,'a')
2641     #     for line in proc.stdout.readlines():
2642     #         f.write(line)
2643
2644     #     f.close()
2645     #     return_code=proc.returncode
2646
2647     # except subprocess.TimeoutExpired:
2648
2649     #     f=open(log_file,'a')
2650     #     for line in proc.stdout.readlines():
2651     #         f.write(line)
2652
2653     #     f.write('Simulation stoppe because of TimeOut')
2654     #     f.close()
2655
2656     #     return_code=2
2657
2658     return #return_code
2659
2660 def readObsnamesfile(__filenames=None):
2661     "Provisionnel pour lire les noms des observations d'un fichier"
2662     for __filename in __filenames:
2663         try:
2664             df = pandas.read_csv(__filename, sep = ";", header =0) #so that repeated names are not modified
2665             obsnames_infile = df.loc[0, :].values.tolist()[2:]
2666         except:
2667             df_tmp = pandas.read_csv(__filename, sep = ";", header =1)
2668             with open(__filename, 'r') as infile:
2669                 readie=csv.reader(infile, delimiter=';')
2670                 with open('tmp_obs_file.csv', 'wt', newline='') as output:
2671                     outwriter=csv.writer(output, delimiter=';')
2672                     list_row = ['#'] + df_tmp.columns
2673                     outwriter.writerow(list_row)
2674                     for row in readie:
2675                          outwriter.writerow(row)
2676
2677             df = pandas.read_csv('tmp_obs_file.csv', sep = ";", header =0) #so that repeated names are not modified
2678             df = df.iloc[1: , :]
2679             obsnames_infile = df.iloc[0, :].values.tolist()[2:]
2680
2681             try:
2682                 os.remove('tmp_obs_file.csv')
2683             except:
2684                 pass
2685     return obsnames_infile
2686
2687 def readDatesHours(__filename=None):
2688     "Provisionnel pour lire les heures et les dates d'une simulation dynamique"
2689     df = pandas.read_csv(__filename, sep = ";", header =1)
2690     dates = numpy.array(df['Date'])
2691     hours = numpy.array(df['Hour'])
2692     return (dates,hours)
2693
2694 def readMultiDatesHours(__filenames=None):
2695     "Provisionnel pour lire les heures et les dates d'une simulation dynamique dans le cas multi-observations"
2696     Multidates =[]
2697     Multihours = []
2698     for __filename in __filenames:
2699         df = pandas.read_csv(__filename, sep = ";", header =1)
2700         dates = numpy.array(df[df.columns[0]])
2701         hours = numpy.array(df[df.columns[1]])
2702         Multidates.append( dates )
2703         Multihours.append( hours )
2704     return (Multidates,Multihours)
2705
2706 def find_nearest_index(array, value):
2707     "Trouver le point de la simulation qui se rapproche le plus aux observations"
2708     array = numpy.asarray(array)
2709     idx = (numpy.abs(array - value)).argmin()
2710     return [idx, array[idx]]
2711 # ==============================================================================
2712 if __name__ == "__main__":
2713     print('\n  AUTODIAGNOSTIC\n  ==============\n')
2714     print("  Module name..............: %s"%name)
2715     print("  Module version...........: %s"%version)
2716     print("  Module year..............: %s"%year)
2717     print("")
2718     print("  Configuration attributes.:")
2719     for __k, __v in _configuration.items():
2720         print("%26s : %s"%(__k,__v))
2721     print("")