Coverage for pySDC/projects/soft_failure/generate_statistics.py: 75%
199 statements
« prev ^ index » next coverage.py v7.5.0, created at 2024-04-29 09:02 +0000
1from __future__ import division
3import dill
4import numpy as np
6from pySDC.helpers.stats_helper import get_sorted
8from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
9from pySDC.implementations.problem_classes.GeneralizedFisher_1D_FD_implicit import generalized_fisher
10from pySDC.implementations.problem_classes.HeatEquation_ND_FD import heatNd_unforced
11from pySDC.implementations.problem_classes.Van_der_Pol_implicit import vanderpol
12from pySDC.projects.soft_failure.FaultHooks import fault_hook
13from pySDC.projects.soft_failure.implicit_sweeper_faults import implicit_sweeper_faults
14from pySDC.projects.soft_failure.visualization_helper import (
15 show_residual_across_simulation,
16 show_min_max_residual_across_simulation,
17 show_iter_hist,
18)
21def diffusion_setup():
22 """
23 Setup routine for diffusion test
24 """
25 # initialize level parameters
26 level_params = dict()
27 level_params['restol'] = 1e-10
28 level_params['dt'] = 0.25
29 level_params['nsweeps'] = 1
31 # initialize sweeper parameters
32 sweeper_params = dict()
33 sweeper_params['quad_type'] = 'RADAU-RIGHT'
34 sweeper_params['num_nodes'] = 3
35 sweeper_params['QI'] = 'LU' # For the IMEX sweeper, the LU-trick can be activated for the implicit part
36 sweeper_params['initial_guess'] = 'spread'
37 sweeper_params['detector_threshold'] = 1e-10
39 # initialize problem parameters
40 problem_params = dict()
41 problem_params['nu'] = 0.1 # diffusion coefficient
42 problem_params['freq'] = 4 # frequency for the test value
43 problem_params['nvars'] = 127 # number of degrees of freedom for each level
44 problem_params['bc'] = 'dirichlet-zero' # boundary conditions
46 # initialize step parameters
47 step_params = dict()
48 step_params['maxiter'] = 50
50 # initialize controller parameters
51 controller_params = dict()
52 controller_params['logger_level'] = 30
54 # fill description dictionary for easy step instantiation
55 description = dict()
56 description['problem_class'] = heatNd_unforced # pass problem class
57 description['problem_params'] = problem_params # pass problem parameters
58 description['sweeper_class'] = implicit_sweeper_faults # pass sweeper
59 description['sweeper_params'] = sweeper_params # pass sweeper parameters
60 description['level_params'] = level_params # pass level parameters
61 description['step_params'] = step_params # pass step parameters
63 return description, controller_params
66def reaction_setup():
67 """
68 Setup routine for diffusion-reaction test with Newton solver
69 """
70 # initialize level parameters
71 level_params = dict()
72 level_params['restol'] = 1e-10
73 level_params['dt'] = 0.25
74 level_params['nsweeps'] = 1
76 # initialize sweeper parameters
77 sweeper_params = dict()
78 sweeper_params['quad_type'] = 'RADAU-RIGHT'
79 sweeper_params['num_nodes'] = 3
80 sweeper_params['QI'] = 'LU' # For the IMEX sweeper, the LU-trick can be activated for the implicit part
81 sweeper_params['initial_guess'] = 'spread'
82 sweeper_params['detector_threshold'] = 1e-10
84 # initialize problem parameters
85 problem_params = dict()
86 problem_params['nu'] = 1.0
87 problem_params['lambda0'] = 2.0
88 problem_params['newton_maxiter'] = 20
89 problem_params['newton_tol'] = 1e-10
90 problem_params['stop_at_nan'] = False
91 problem_params['interval'] = (-5, 5)
92 problem_params['nvars'] = 127
94 # initialize step parameters
95 step_params = dict()
96 step_params['maxiter'] = 50
98 # initialize controller parameters
99 controller_params = dict()
100 controller_params['logger_level'] = 30
102 # fill description dictionary for easy step instantiation
103 description = dict()
104 description['problem_class'] = generalized_fisher # pass problem class
105 description['problem_params'] = problem_params # pass problem parameters
106 description['sweeper_class'] = implicit_sweeper_faults # pass sweeper
107 description['sweeper_params'] = sweeper_params # pass sweeper parameters
108 description['level_params'] = level_params # pass level parameters
109 description['step_params'] = step_params # pass step parameters
111 return description, controller_params
114def vanderpol_setup():
115 """
116 Van der Pol's oscillator
117 """
118 # initialize level parameters
119 level_params = dict()
120 level_params['restol'] = 1e-10
121 level_params['dt'] = 0.25
122 level_params['nsweeps'] = 1
124 # initialize sweeper parameters
125 sweeper_params = dict()
126 sweeper_params['quad_type'] = 'RADAU-RIGHT'
127 sweeper_params['num_nodes'] = 3
128 sweeper_params['QI'] = 'LU'
129 sweeper_params['detector_threshold'] = 1e-10
131 # initialize problem parameters
132 problem_params = dict()
133 problem_params['newton_tol'] = 1e-10
134 problem_params['newton_maxiter'] = 50
135 problem_params['stop_at_nan'] = False
136 problem_params['mu'] = 18
137 problem_params['u0'] = (1.0, 0.0)
139 # initialize step parameters
140 step_params = dict()
141 step_params['maxiter'] = 50
143 # initialize controller parameters
144 controller_params = dict()
145 controller_params['logger_level'] = 30
147 # Fill description dictionary for easy hierarchy creation
148 description = dict()
149 description['problem_class'] = vanderpol
150 description['problem_params'] = problem_params
151 description['sweeper_class'] = implicit_sweeper_faults
152 description['sweeper_params'] = sweeper_params
153 description['level_params'] = level_params
154 description['step_params'] = step_params
156 return description, controller_params
159def run_clean_simulations(type=None):
160 """
161 A simple code to run fault-free simulations
163 Args:
164 type (str): setup type
165 f: file handler
166 """
168 if type == 'diffusion':
169 description, controller_params = diffusion_setup()
170 elif type == 'reaction':
171 description, controller_params = reaction_setup()
172 elif type == 'vanderpol':
173 description, controller_params = vanderpol_setup()
174 else:
175 raise ValueError('No valid setup type provided, aborting..')
177 # set time parameters
178 t0 = 0.0
179 Tend = description['level_params']['dt']
181 # instantiate controller
182 controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description)
184 # get initial values on finest level
185 P = controller.MS[0].levels[0].prob
186 uinit = P.u_exact(t0)
188 # this is where the iteration is happening
189 uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
191 # filter statistics by type (number of iterations)
192 iter_counts = get_sorted(stats, type='niter', sortby='time')
194 # print('This clean run took %s iterations!' % iter_counts[0][1])
196 return iter_counts[0][1]
199def run_faulty_simulations(type=None, niters=None, cwd=''):
200 """
201 A simple program to run faulty simulations
203 Args:
204 type (str): setup type
205 niters (int): number of iterations in clean run
206 f: file handler
207 cwd (str): current workind directory
208 """
210 if type == 'diffusion':
211 description, controller_params = diffusion_setup()
212 elif type == 'reaction':
213 description, controller_params = reaction_setup()
214 elif type == 'vanderpol':
215 description, controller_params = vanderpol_setup()
216 else:
217 raise ValueError('No valid setup type provided, aborting..')
219 # set time parameters
220 t0 = 0.0
221 Tend = description['level_params']['dt']
223 filehandle_injections = open(cwd + 'data/dump_injections_' + type + '.txt', 'w')
225 controller_params['hook_class'] = fault_hook
226 description['sweeper_params']['allow_fault_correction'] = True
227 description['sweeper_params']['dump_injections_filehandle'] = filehandle_injections
228 description['sweeper_params']['niters'] = niters
230 # instantiate controller
231 controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description)
233 # get initial values on finest level
234 P = controller.MS[0].levels[0].prob
235 uinit = P.u_exact(t0)
237 # number of runs
238 nruns = 500
239 results = []
240 for _ in range(nruns):
241 # this is where the iteration is happening
242 uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
244 results.append(stats)
246 filehandle_injections.close()
248 dill.dump(results, open(cwd + "data/results_" + type + ".pkl", "wb"))
251def process_statistics(type=None, cwd=''):
252 results = dill.load(open(cwd + "data/results_" + type + ".pkl", "rb"))
254 # get minimal length of residual vector
255 minlen = 1000
256 nruns = 0
257 for stats in results:
258 residuals = get_sorted(stats, type='residual_post_iteration', sortby='iter')
259 minlen = min(minlen, len(residuals))
260 nruns += 1
262 # initialize minimal residual vector
263 minres = np.zeros(minlen)
264 minres[:] = 1000
265 # initialize maximal residual vector
266 maxres = np.zeros(minlen)
267 # initialize mean residual vector
268 meanres = np.zeros(minlen)
269 # initialize median residual vector
270 medianres = np.zeros(minlen)
271 # initialize helper list
272 median_list = [[] for _ in range(minlen)]
274 for stats in results:
275 # Some black magic to extract fault stats out of monstrous stats object
276 # fault_stats = get_sorted(stats, type='fault_stats', sortby='type')[0][1]
277 # Some black magic to extract residuals dependent on iteration out of monstrous stats object
278 residuals_iter = get_sorted(stats, type='residual_post_iteration', sortby='iter')
279 # extract residuals ouf of residuals_iter
280 residuals = np.array([item[1] for item in residuals_iter])
282 # calculate minimal, maximal, mean residual vectors
283 for i in range(minlen):
284 if np.isnan(residuals[i]) or np.isinf(residuals[i]):
285 residuals[i] = 1000
286 minres[i] = min(minres[i], residuals[i])
287 maxres[i] = max(maxres[i], residuals[i])
288 meanres[i] += residuals[i]
289 median_list[i].append(residuals[i])
291 # Example output of what we now can do
292 # print(fault_stats.nfaults_injected_u, fault_stats.nfaults_injected_f, fault_stats.nfaults_detected,
293 # fault_stats.nfalse_positives, fault_stats.nfalse_positives_in_correction,
294 # fault_stats.nfaults_missed, fault_stats.nclean_steps)
296 # print()
298 # call helper routine to produce residual plot
299 # fname = 'residuals.png'
300 fname = cwd + 'data/' + type + '_' + str(nruns) + '_' + 'runs' + '_' + 'residuals.png'
301 show_residual_across_simulation(stats=stats, fname=fname)
302 meanres /= nruns
303 # print(minres)
304 # print(maxres)
305 # print(meanres)
306 # calculate median residual vector
307 for i in range(minlen):
308 medianres[i] = np.median(median_list[i])
309 # print(median_list)
310 # print(medianres)
311 # call helper routine to produce residual plot of minres, maxres, meanres and medianres
312 # fname = 'min_max_residuals.png'
313 fname = cwd + 'data/' + type + '_' + str(nruns) + '_' + 'runs' + '_' + 'min_max_residuals.png'
314 show_min_max_residual_across_simulation(
315 fname=fname, minres=minres, maxres=maxres, meanres=meanres, medianres=medianres, maxiter=minlen
316 )
318 # calculate maximum number of iterations per test run
319 maxiter = []
320 for stats in results:
321 residuals = get_sorted(stats, type='residual_post_iteration', sortby='iter')
322 maxiters = max(np.array([item[0] for item in residuals]))
323 maxiter.append(maxiters)
324 # print(maxiter)
325 # call helper routine to produce histogram of maxiter
326 # fname = 'iter_hist.png'
327 fname = cwd + 'data/' + type + '_' + str(nruns) + '_' + 'runs' + '_' + 'iter_hist.png'
328 show_iter_hist(fname=fname, maxiter=maxiter, nruns=nruns)
330 # initialize sum of nfaults_detected
331 nfd = 0
332 # initialize sum of nfalse_positives
333 nfp = 0
334 # initialize sum of nfaults_missed
335 nfm = 0
336 # initialize sum of nfalse_positives_in_correction
337 nfpc = 0
338 # calculate sum of nfaults_detected, sum of nfalse_positives, sum of nfaults_missed
339 for stats in results:
340 # Some black magic to extract fault stats out of monstrous stats object
341 fault_stats = get_sorted(stats, type='fault_stats', sortby='type')[0][1]
342 nfd += fault_stats.nfaults_detected
343 nfp += fault_stats.nfalse_positives
344 nfm += fault_stats.nfaults_missed
345 nfpc += fault_stats.nfalse_positives_in_correction
347 g = open(cwd + 'data/' + type + '_' + str(nruns) + '_' + 'runs' + '_' + 'Statistics.txt', 'w')
348 out = 'Type: ' + type + ' ' + str(nruns) + ' runs'
349 g.write(out + '\n')
350 # detector metrics (Sloan, Kumar, Bronevetsky 2012)
351 # nfaults_detected
352 out = 'true positives: ' + str(nfd)
353 g.write(out + '\n')
354 # nfaults_positives
355 out = 'false positives: ' + str(nfp)
356 g.write(out + '\n')
357 # nfaults_missed
358 out = 'false negatives: ' + str(nfm)
359 g.write(out + '\n')
360 # nfalse_positives_in_correction
361 out = 'false positives in correction: ' + str(nfpc)
362 g.write(out + '\n')
363 # F-Score
364 f_score = 2 * nfd / (2 * nfd + nfp + nfm)
365 out = 'F-Score: ' + str(f_score)
366 g.write(out + '\n')
367 # false positive rate (FPR)
368 fpr = nfp / (nfd + nfp)
369 out = 'False positive rate: ' + str(fpr)
370 g.write(out + '\n')
371 # true positive rate (TPR)
372 tpr = nfd / (nfd + nfp)
373 out = 'True positive rate: ' + str(tpr)
374 g.write(out + '\n')
375 g.close()
378def main():
379 # type = 'diffusion'
380 # niters = run_clean_simulations(type=type)
381 # run_faulty_simulations(type=type, niters=niters)
382 # process_statistics(type=type)
384 # type = 'reaction'
385 # niters = run_clean_simulations(type=type)
386 # run_faulty_simulations(type=type, niters=niters)
387 # process_statistics(type=type)
389 type = 'vanderpol'
390 niters = run_clean_simulations(type=type)
391 run_faulty_simulations(type=type, niters=niters)
392 process_statistics(type=type)
395if __name__ == "__main__":
396 main()