10691143420

Committed 03 Sep 2024 09:48PM UTC coverage: 76.813% (+0.8%) from 76.01%

Build # 10691143420

Build Type

Pull #35

github

Committed by

web-flow

Commit Message

Merge 1ebfe99fb into 81a2299af

Pull Request Pull Request #35: Bump actions/download-artifact from 2 to 4.1.7 in /.github/workflows

Run Details

1928 of 2510 relevant lines covered (76.81%)

10.59 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

79.67

/src/welltestpy/estimate/steady_lib.py

"""welltestpy subpackage providing base classes for steady state estimations."""
import os
import time as timemodule
from copy import deepcopy as dcopy

import numpy as np
import spotpy

from ..data import testslib
from ..process import processlib
from ..tools import plotter
from . import spotpylib

__all__ = [
    "SteadyPumping",
]


class SteadyPumping:
    """Class to estimate steady Type-Curve parameters.

    Parameters
    ----------
    name : :class:`str`
        Name of the Estimation.
    campaign : :class:`welltestpy.data.Campaign`
        The pumping test campaign which should be used to estimate the
        parameters
    type_curve : :any:`callable`
        The given type-curve. Output will be reshaped to flat array.
    val_ranges : :class:`dict`
        Dictionary containing the fit-ranges for each value in the type-curve.
        Names should be as in the type-curve signature.
        Ranges should be a tuple containing min and max value.
    make_steady : :class:`bool`, optional
        State if the tests should be converted to steady observations.
        See: :any:`PumpingTest.make_steady`.
        Default: True
    val_fix : :class:`dict` or :any:`None`
        Dictionary containing fixed values for the type-curve.
        Names should be as in the type-curve signature.
        Default: None
    val_fit_type : :class:`dict` or :any:`None`
        Dictionary containing fitting transformation type for each value.
        Names should be as in the type-curve signature.
        val_fit_type can be "lin", "log", "exp", "sqrt", "quad", "inv"
        or a tuple of two callable functions where the
        first is the transformation and the second is its inverse.
        "log" is for example equivalent to ``(np.log, np.exp)``.
        By default, values will be fitted linear.
        Default: None
    val_fit_name : :class:`dict` or :any:`None`
        Display name of the fitting transformation.
        Will be the val_fit_type string if it is a predefined one,
        or ``f`` if it is a given callable as default for each value.
        Default: None
    val_plot_names : :class:`dict` or :any:`None`
        Dictionary containing keyword names in the type-curve for each value.

            {value-name: string for plot legend}

        This is useful to get better plots.
        By default, parameter names will be value names.
        Default: None
    testinclude : :class:`dict`, optional
        Dictionary of which tests should be included. If ``None`` is given,
        all available tests are included.
        Default: ``None``
    generate : :class:`bool`, optional
        State if time stepping, processed observation data and estimation
        setup should be generated with default values.
        Default: ``False``
    """

    def __init__(
        self,
        name,
        campaign,
        type_curve,
        val_ranges,
        make_steady=True,
        val_fix=None,
        val_fit_type=None,
        val_fit_name=None,
        val_plot_names=None,
        testinclude=None,
        generate=False,
    ):
        val_fix = {} if val_fix is None else val_fix
        self.setup_kw = {
            "type_curve": type_curve,
            "val_ranges": val_ranges,
            "val_fix": val_fix,
            "val_fit_type": val_fit_type,
            "val_fit_name": val_fit_name,
            "val_plot_names": val_plot_names,
        }
        """:class:`dict`: TypeCurve Spotpy Setup definition"""
        self.name = name
        """:class:`str`: Name of the Estimation"""
        self.campaign_raw = dcopy(campaign)
        """:class:`welltestpy.data.Campaign`:\
        Copy of the original input campaign"""
        self.campaign = dcopy(campaign)
        """:class:`welltestpy.data.Campaign`:\
        Copy of the input campaign to be modified"""

        self.prate = None
        """:class:`float`: Pumpingrate at the pumping well"""

        self.rad = None
        """:class:`numpy.ndarray`: array of the radii from the wells"""
        self.data = None
        """:class:`numpy.ndarray`: observation data"""
        self.radnames = None
        """:class:`numpy.ndarray`: names of the radii well combination"""
        self.r_ref = None
        """:class:`float`: reference radius of the biggest distance"""
        self.h_ref = None
        """:class:`float`: reference head at the biggest distance"""

        self.estimated_para = {}
        """:class:`dict`: estimated parameters by name"""
        self.result = None
        """:class:`list`: result of the spotpy estimation"""
        self.sens = None
        """:class:`dict`: result of the spotpy sensitivity analysis"""
        self.testinclude = {}
        """:class:`dict`: dictionary of which tests should be included"""

        if testinclude is None:
            tests = list(self.campaign.tests.keys())
            self.testinclude = {}
            for test in tests:
                self.testinclude[test] = self.campaign.tests[
                    test
                ].observationwells
        elif not isinstance(testinclude, dict):
            self.testinclude = {}
            for test in testinclude:
                self.testinclude[test] = self.campaign.tests[
                    test
                ].observationwells
        else:
            self.testinclude = testinclude

        for test in self.testinclude:
            if not isinstance(self.campaign.tests[test], testslib.PumpingTest):
                raise ValueError(test + " is not a pumping test.")
            if make_steady is not False:
                if make_steady is True:
                    make_steady = "latest"
                self.campaign.tests[test].make_steady(make_steady)
            if not self.campaign.tests[test].constant_rate:
                raise ValueError(test + " is not a constant rate test.")
            if (
                not self.campaign.tests[test].state(
                    wells=self.testinclude[test]
                )
                == "steady"
            ):
                raise ValueError(test + ": selection is not steady.")

        rwell_list = []
        rinf_list = []
        for test in self.testinclude:
            pwell = self.campaign.tests[test].pumpingwell
            rwell_list.append(self.campaign.wells[pwell].radius)
            rinf_list.append(self.campaign.tests[test].radius)
        self.rwell = min(rwell_list)
        """:class:`float`: radius of the pumping wells"""
        self.rinf = max(rinf_list)
        """:class:`float`: radius of the furthest wells"""

        if generate:
            self.setpumprate()
            self.gen_data()
            self.gen_setup()

    def setpumprate(self, prate=-1.0):
        """Set a uniform pumping rate at all pumpingwells wells.

        We assume linear scaling by the pumpingrate.

        Parameters
        ----------
        prate : :class:`float`, optional
            Pumping rate. Default: ``-1.0``
        """
        for test in self.testinclude:
            processlib.normpumptest(
                self.campaign.tests[test], pumpingrate=prate
            )
        self.prate = prate

    def gen_data(self):
        """Generate the observed drawdown.

        It will also generate an array containing all radii of all well
        combinations.
        """
        rad = np.array([])
        data = np.array([])

        radnames = []

        for test in self.testinclude:
            pwell = self.campaign.wells[self.campaign.tests[test].pumpingwell]
            for obs in self.testinclude[test]:
                temphead = self.campaign.tests[test].observations[obs]()
                data = np.hstack((data, temphead))

                owell = self.campaign.wells[obs]
                if pwell == owell:
                    temprad = pwell.radius
                else:
                    temprad = pwell - owell
                rad = np.hstack((rad, temprad))

                tempname = (self.campaign.tests[test].pumpingwell, obs)
                radnames.append(tempname)

        # sort everything by the radii
        idx = rad.argsort()
        radnames = np.array(radnames)
        self.rad = rad[idx]
        self.data = data[idx]
        self.radnames = radnames[idx]
        self.r_ref = self.rad[-1]
        self.h_ref = self.data[-1]

    def gen_setup(
        self,
        prate_kw="rate",
        rad_kw="rad",
        r_ref_kw="r_ref",
        h_ref_kw="h_ref",
        dummy=False,
    ):
        """Generate the Spotpy Setup.

        Parameters
        ----------
        prate_kw : :class:`str`, optional
            Keyword name for the pumping rate in the used type curve.
            Default: "rate"
        rad_kw : :class:`str`, optional
            Keyword name for the radius in the used type curve.
            Default: "rad"
        r_ref_kw : :class:`str`, optional
            Keyword name for the reference radius in the used type curve.
            Default: "r_ref"
        h_ref_kw : :class:`str`, optional
            Keyword name for the reference head in the used type curve.
            Default: "h_ref"
        dummy : :class:`bool`, optional
            Add a dummy parameter to the model. This could be used to equalize
            sensitivity analysis.
            Default: False
        """
        self.extra_kw_names = {"rad": rad_kw}
        setup_kw = dcopy(self.setup_kw)  # create a copy here
        setup_kw["val_fix"].setdefault(prate_kw, self.prate)
        setup_kw["val_fix"].setdefault(rad_kw, self.rad)
        setup_kw["val_fix"].setdefault(r_ref_kw, self.r_ref)
        setup_kw["val_fix"].setdefault(h_ref_kw, self.h_ref)
        setup_kw.setdefault("data", self.data)
        setup_kw["dummy"] = dummy
        self.setup = spotpylib.TypeCurve(**setup_kw)

    def run(
        self,
        rep=5000,
        parallel="seq",
        run=True,
        folder=None,
        dbname=None,
        traceplotname=None,
        fittingplotname=None,
        interactplotname=None,
        estname=None,
        plot_style="WTP",
    ):
        """Run the estimation.

        Parameters
        ----------
        rep : :class:`int`, optional
            The number of repetitions within the SCEua algorithm in spotpy.
            Default: ``5000``
        parallel : :class:`str`, optional
            State if the estimation should be run in parallel or not. Options:

                    * ``"seq"``: sequential on one CPU
                    * ``"mpi"``: use the mpi4py package

            Default: ``"seq"``
        run : :class:`bool`, optional
            State if the estimation should be executed. Otherwise all plots
            will be done with the previous results.
            Default: ``True``
        folder : :class:`str`, optional
            Path to the output folder. If ``None`` the CWD is used.
            Default: ``None``
        dbname : :class:`str`, optional
            File-name of the database of the spotpy estimation.
            If ``None``, it will be the current time +
            ``"_db"``.
            Default: ``None``
        traceplotname : :class:`str`, optional
            File-name of the parameter trace plot of the spotpy estimation.
            If ``None``, it will be the current time +
            ``"_paratrace.pdf"``.
            Default: ``None``
        fittingplotname : :class:`str`, optional
            File-name of the fitting plot of the estimation.
            If ``None``, it will be the current time +
            ``"_fit.pdf"``.
            Default: ``None``
        interactplotname : :class:`str`, optional
            File-name of the parameter interaction plot
            of the spotpy estimation.
            If ``None``, it will be the current time +
            ``"_parainteract.pdf"``.
            Default: ``None``
        estname : :class:`str`, optional
            File-name of the results of the spotpy estimation.
            If ``None``, it will be the current time +
            ``"_estimate"``.
            Default: ``None``
        plot_style : str, optional
            Plot style. The default is "WTP".
        """
        if self.setup.dummy:
            raise ValueError(
                "Estimate: for parameter estimation"
                " you can't use a dummy parameter."
            )
        act_time = timemodule.strftime("%Y-%m-%d_%H-%M-%S")

        # generate the filenames
        if folder is None:
            folder = os.path.join(os.getcwd(), self.name)
        folder = os.path.abspath(folder)
        if not os.path.exists(folder):
            os.makedirs(folder)

        if dbname is None:
            dbname = os.path.join(folder, act_time + "_db")
        elif not os.path.isabs(dbname):
            dbname = os.path.join(folder, dbname)
        if traceplotname is None:
            traceplotname = os.path.join(folder, act_time + "_paratrace.pdf")
        elif not os.path.isabs(traceplotname):
            traceplotname = os.path.join(folder, traceplotname)
        if fittingplotname is None:
            fittingplotname = os.path.join(folder, act_time + "_fit.pdf")
        elif not os.path.isabs(fittingplotname):
            fittingplotname = os.path.join(folder, fittingplotname)
        if interactplotname is None:
            interactplotname = os.path.join(folder, act_time + "_interact.pdf")
        elif not os.path.isabs(interactplotname):
            interactplotname = os.path.join(folder, interactplotname)
        if estname is None:
            paraname = os.path.join(folder, act_time + "_estimate.txt")
        elif not os.path.isabs(estname):
            paraname = os.path.join(folder, estname)

        # generate the parameter-names for plotting
        paranames = dcopy(self.setup.para_names)
        paralabels = []
        for name in paranames:
            p_label = self.setup.val_plot_names[name]
            fit_n = self.setup.val_fit_name[name]
            paralabels.append(f"{fit_n}({p_label})" if fit_n else p_label)

        if parallel == "mpi":
            # send the dbname of rank0
            from mpi4py import MPI

            comm = MPI.COMM_WORLD
            rank = comm.Get_rank()
            size = comm.Get_size()
            if rank == 0:
                print(rank, "send dbname:", dbname)
                for i in range(1, size):
                    comm.send(dbname, dest=i, tag=0)
            else:
                dbname = comm.recv(source=0, tag=0)
                print(rank, "got dbname:", dbname)
        else:
            rank = 0

        # initialize the sampler
        sampler = spotpy.algorithms.sceua(
            self.setup,
            dbname=dbname,
            dbformat="csv",
            parallel=parallel,
            save_sim=True,
            db_precision=np.float64,
        )
        # start the estimation with the sce-ua algorithm
        if run:
            sampler.sample(rep, ngs=10, kstop=100, pcento=1e-4, peps=1e-3)

        if rank == 0:
            if run:
                self.result = sampler.getdata()
            else:
                self.result = np.genfromtxt(
                    dbname + ".csv", delimiter=",", names=True
                )
            para_opt = spotpy.analyser.get_best_parameterset(
                self.result, maximize=False
            )
            void_names = para_opt.dtype.names
            para = []
            header = []
            for name in void_names:
                para.append(para_opt[0][name])
                fit_n = self.setup.val_fit_name[name[3:]]
                header.append(f"{fit_n}-{name[3:]}" if fit_n else name[3:])
                self.estimated_para[name[3:]] = para[-1]
            np.savetxt(paraname, para, header=" ".join(header))
            # plot the estimation-results
            plotter.plotparatrace(
                result=self.result,
                parameternames=paranames,
                parameterlabels=paralabels,
                stdvalues=self.estimated_para,
                plotname=traceplotname,
                style=plot_style,
            )
            plotter.plotfit_steady(
                setup=self.setup,
                data=self.data,
                para=self.estimated_para,
                rad=self.rad,
                radnames=self.radnames,
                extra=self.extra_kw_names,
                plotname=fittingplotname,
                style=plot_style,
            )
            plotter.plotparainteract(
                self.result, paralabels, interactplotname, style=plot_style
            )

    def sensitivity(
        self,
        rep=None,
        parallel="seq",
        folder=None,
        dbname=None,
        plotname=None,
        traceplotname=None,
        sensname=None,
        plot_style="WTP",
    ):
        """Run the sensitivity analysis.

        Parameters
        ----------
        rep : :class:`int`, optional
            The number of repetitions within the FAST algorithm in spotpy.
            Default: estimated
        parallel : :class:`str`, optional
            State if the estimation should be run in parallel or not. Options:

                    * ``"seq"``: sequential on one CPU
                    * ``"mpi"``: use the mpi4py package

            Default: ``"seq"``
        folder : :class:`str`, optional
            Path to the output folder. If ``None`` the CWD is used.
            Default: ``None``
        dbname : :class:`str`, optional
            File-name of the database of the spotpy estimation.
            If ``None``, it will be the current time +
            ``"_sensitivity_db"``.
            Default: ``None``
        plotname : :class:`str`, optional
            File-name of the result plot of the sensitivity analysis.
            If ``None``, it will be the current time +
            ``"_sensitivity.pdf"``.
            Default: ``None``
        traceplotname : :class:`str`, optional
            File-name of the parameter trace plot of the spotpy sensitivity
            analysis.
            If ``None``, it will be the current time +
            ``"_senstrace.pdf"``.
            Default: ``None``
        sensname : :class:`str`, optional
            File-name of the results of the FAST estimation.
            If ``None``, it will be the current time +
            ``"_estimate"``.
            Default: ``None``
        plot_style : str, optional
            Plot style. The default is "WTP".
        """
        if len(self.setup.para_names) == 1 and not self.setup.dummy:
            raise ValueError(
                "Sensitivity: for estimation with only one parameter"
                " you have to use a dummy parameter."
            )
        if rep is None:
            rep = spotpylib.fast_rep(
                len(self.setup.para_names) + int(self.setup.dummy)
            )

        act_time = timemodule.strftime("%Y-%m-%d_%H-%M-%S")
        # generate the filenames
        if folder is None:
            folder = os.path.join(os.getcwd(), self.name)
        folder = os.path.abspath(folder)
        if not os.path.exists(folder):
            os.makedirs(folder)

        if dbname is None:
            dbname = os.path.join(folder, act_time + "_sensitivity_db")
        elif not os.path.isabs(dbname):
            dbname = os.path.join(folder, dbname)
        if plotname is None:
            plotname = os.path.join(folder, act_time + "_sensitivity.pdf")
        elif not os.path.isabs(plotname):
            plotname = os.path.join(folder, plotname)
        if traceplotname is None:
            traceplotname = os.path.join(folder, act_time + "_senstrace.pdf")
        elif not os.path.isabs(traceplotname):
            traceplotname = os.path.join(folder, traceplotname)
        if sensname is None:
            sensname = os.path.join(folder, act_time + "_FAST_estimate.txt")
        elif not os.path.isabs(sensname):
            sensname = os.path.join(folder, sensname)

        sens_base, sens_ext = os.path.splitext(sensname)
        sensname1 = sens_base + "_S1" + sens_ext

        # generate the parameter-names for plotting
        paranames = dcopy(self.setup.para_names)
        paralabels = []
        for name in paranames:
            p_label = self.setup.val_plot_names[name]
            fit_n = self.setup.val_fit_name[name]
            paralabels.append(f"{fit_n}({p_label})" if fit_n else p_label)

        if self.setup.dummy:
            paranames.append("dummy")
            paralabels.append("dummy")

        if parallel == "mpi":
            # send the dbname of rank0
            from mpi4py import MPI

            comm = MPI.COMM_WORLD
            rank = comm.Get_rank()
            size = comm.Get_size()
            if rank == 0:
                print(rank, "send dbname:", dbname)
                for i in range(1, size):
                    comm.send(dbname, dest=i, tag=0)
            else:
                dbname = comm.recv(source=0, tag=0)
                print(rank, "got dbname:", dbname)
        else:
            rank = 0

        # initialize the sampler
        sampler = spotpy.algorithms.fast(
            self.setup,
            dbname=dbname,
            dbformat="csv",
            parallel=parallel,
            save_sim=True,
            db_precision=np.float64,
        )
        sampler.sample(rep)

        if rank == 0:
            data = sampler.getdata()
            parmin = sampler.parameter()["minbound"]
            parmax = sampler.parameter()["maxbound"]
            bounds = list(zip(parmin, parmax))
            sens_est = sampler.analyze(
                bounds, np.nan_to_num(data["like1"]), len(paranames), paranames
            )
            self.sens = {}
            for sen_typ in sens_est:
                self.sens[sen_typ] = {
                    par: sen for par, sen in zip(paranames, sens_est[sen_typ])
                }
            header = " ".join(paranames)
            np.savetxt(sensname, sens_est["ST"], header=header)
            np.savetxt(sensname1, sens_est["S1"], header=header)
            plotter.plotsensitivity(
                paralabels, sens_est, plotname, style=plot_style
            )
            plotter.plotparatrace(
                data,
                parameternames=paranames,
                parameterlabels=paralabels,
                stdvalues=None,
                plotname=traceplotname,
                style=plot_style,
            )

1	"""welltestpy subpackage providing base classes for steady state estimations."""
2	import os	14✔
3	import time as timemodule	14✔
4	from copy import deepcopy as dcopy	14✔
5
6	import numpy as np	14✔
7	import spotpy	14✔
8
9	from ..data import testslib	14✔
10	from ..process import processlib	14✔
11	from ..tools import plotter	14✔
12	from . import spotpylib	14✔
13
14	__all__ = [	14✔
15	"SteadyPumping",
16	]
17
18
19	class SteadyPumping:	14✔
20	"""Class to estimate steady Type-Curve parameters.
21
22	Parameters
23	----------
24	name : :class:`str`
25	Name of the Estimation.
26	campaign : :class:`welltestpy.data.Campaign`
27	The pumping test campaign which should be used to estimate the
28	parameters
29	type_curve : :any:`callable`
30	The given type-curve. Output will be reshaped to flat array.
31	val_ranges : :class:`dict`
32	Dictionary containing the fit-ranges for each value in the type-curve.
33	Names should be as in the type-curve signature.
34	Ranges should be a tuple containing min and max value.
35	make_steady : :class:`bool`, optional
36	State if the tests should be converted to steady observations.
37	See: :any:`PumpingTest.make_steady`.
38	Default: True
39	val_fix : :class:`dict` or :any:`None`
40	Dictionary containing fixed values for the type-curve.
41	Names should be as in the type-curve signature.
42	Default: None
43	val_fit_type : :class:`dict` or :any:`None`
44	Dictionary containing fitting transformation type for each value.
45	Names should be as in the type-curve signature.
46	val_fit_type can be "lin", "log", "exp", "sqrt", "quad", "inv"
47	or a tuple of two callable functions where the
48	first is the transformation and the second is its inverse.
49	"log" is for example equivalent to ``(np.log, np.exp)``.
50	By default, values will be fitted linear.
51	Default: None
52	val_fit_name : :class:`dict` or :any:`None`
53	Display name of the fitting transformation.
54	Will be the val_fit_type string if it is a predefined one,
55	or ``f`` if it is a given callable as default for each value.
56	Default: None
57	val_plot_names : :class:`dict` or :any:`None`
58	Dictionary containing keyword names in the type-curve for each value.
59
60	{value-name: string for plot legend}
61
62	This is useful to get better plots.
63	By default, parameter names will be value names.
64	Default: None
65	testinclude : :class:`dict`, optional
66	Dictionary of which tests should be included. If ``None`` is given,
67	all available tests are included.
68	Default: ``None``
69	generate : :class:`bool`, optional
70	State if time stepping, processed observation data and estimation
71	setup should be generated with default values.
72	Default: ``False``
73	"""
74
75	def __init__(	14✔
76	self,
77	name,
78	campaign,
79	type_curve,
80	val_ranges,
81	make_steady=True,
82	val_fix=None,
83	val_fit_type=None,
84	val_fit_name=None,
85	val_plot_names=None,
86	testinclude=None,
87	generate=False,
88	):
89	val_fix = {} if val_fix is None else val_fix	14✔
90	self.setup_kw = {	14✔
91	"type_curve": type_curve,
92	"val_ranges": val_ranges,
93	"val_fix": val_fix,
94	"val_fit_type": val_fit_type,
95	"val_fit_name": val_fit_name,
96	"val_plot_names": val_plot_names,
97	}
98	""":class:`dict`: TypeCurve Spotpy Setup definition"""	6✔
99	self.name = name	14✔
100	""":class:`str`: Name of the Estimation"""	6✔
101	self.campaign_raw = dcopy(campaign)	14✔
102	""":class:`welltestpy.data.Campaign`:\	6✔
103	Copy of the original input campaign"""
104	self.campaign = dcopy(campaign)	14✔
105	""":class:`welltestpy.data.Campaign`:\	6✔
106	Copy of the input campaign to be modified"""
107
108	self.prate = None	14✔
109	""":class:`float`: Pumpingrate at the pumping well"""	6✔
110
111	self.rad = None	14✔
112	""":class:`numpy.ndarray`: array of the radii from the wells"""	6✔
113	self.data = None	14✔
114	""":class:`numpy.ndarray`: observation data"""	6✔
115	self.radnames = None	14✔
116	""":class:`numpy.ndarray`: names of the radii well combination"""	6✔
117	self.r_ref = None	14✔
118	""":class:`float`: reference radius of the biggest distance"""	6✔
119	self.h_ref = None	14✔
120	""":class:`float`: reference head at the biggest distance"""	6✔
121
122	self.estimated_para = {}	14✔
123	""":class:`dict`: estimated parameters by name"""	6✔
124	self.result = None	14✔
125	""":class:`list`: result of the spotpy estimation"""	6✔
126	self.sens = None	14✔
127	""":class:`dict`: result of the spotpy sensitivity analysis"""	6✔
128	self.testinclude = {}	14✔
129	""":class:`dict`: dictionary of which tests should be included"""	6✔
130
131	if testinclude is None:	14✔
132	tests = list(self.campaign.tests.keys())	14✔
133	self.testinclude = {}	14✔
134	for test in tests:	14✔
135	self.testinclude[test] = self.campaign.tests[	14✔
136	test
137	].observationwells
138	elif not isinstance(testinclude, dict):	×
139	self.testinclude = {}	×
140	for test in testinclude:	×
141	self.testinclude[test] = self.campaign.tests[	×
142	test
143	].observationwells
144	else:
145	self.testinclude = testinclude	×
146
147	for test in self.testinclude:	14✔
148	if not isinstance(self.campaign.tests[test], testslib.PumpingTest):	14✔
149	raise ValueError(test + " is not a pumping test.")	×
150	if make_steady is not False:	14✔
151	if make_steady is True:	14✔
152	make_steady = "latest"	14✔
153	self.campaign.tests[test].make_steady(make_steady)	14✔
154	if not self.campaign.tests[test].constant_rate:	14✔
155	raise ValueError(test + " is not a constant rate test.")	×
156	if (	14✔
157	not self.campaign.tests[test].state(
158	wells=self.testinclude[test]
159	)
160	== "steady"
161	):
162	raise ValueError(test + ": selection is not steady.")	×
163
164	rwell_list = []	14✔
165	rinf_list = []	14✔
166	for test in self.testinclude:	14✔
167	pwell = self.campaign.tests[test].pumpingwell	14✔
168	rwell_list.append(self.campaign.wells[pwell].radius)	14✔
169	rinf_list.append(self.campaign.tests[test].radius)	14✔
170	self.rwell = min(rwell_list)	14✔
171	""":class:`float`: radius of the pumping wells"""	6✔
172	self.rinf = max(rinf_list)	14✔
173	""":class:`float`: radius of the furthest wells"""	6✔
174
175	if generate:	14✔
176	self.setpumprate()	14✔
177	self.gen_data()	14✔
178	self.gen_setup()	14✔
179
180	def setpumprate(self, prate=-1.0):	14✔
181	"""Set a uniform pumping rate at all pumpingwells wells.
182
183	We assume linear scaling by the pumpingrate.
184
185	Parameters
186	----------
187	prate : :class:`float`, optional
188	Pumping rate. Default: ``-1.0``
189	"""
190	for test in self.testinclude:	14✔
191	processlib.normpumptest(	14✔
192	self.campaign.tests[test], pumpingrate=prate
193	)
194	self.prate = prate	14✔
195
196	def gen_data(self):	14✔
197	"""Generate the observed drawdown.
198
199	It will also generate an array containing all radii of all well
200	combinations.
201	"""
202	rad = np.array([])	14✔
203	data = np.array([])	14✔
204
205	radnames = []	14✔
206
207	for test in self.testinclude:	14✔
208	pwell = self.campaign.wells[self.campaign.tests[test].pumpingwell]	14✔
209	for obs in self.testinclude[test]:	14✔
210	temphead = self.campaign.tests[test].observations[obs]()	14✔
211	data = np.hstack((data, temphead))	14✔
212
213	owell = self.campaign.wells[obs]	14✔
214	if pwell == owell:	14✔
215	temprad = pwell.radius	14✔
216	else:
217	temprad = pwell - owell	14✔
218	rad = np.hstack((rad, temprad))	14✔
219
220	tempname = (self.campaign.tests[test].pumpingwell, obs)	14✔
221	radnames.append(tempname)	14✔
222
223	# sort everything by the radii
224	idx = rad.argsort()	14✔
225	radnames = np.array(radnames)	14✔
226	self.rad = rad[idx]	14✔
227	self.data = data[idx]	14✔
228	self.radnames = radnames[idx]	14✔
229	self.r_ref = self.rad[-1]	14✔
230	self.h_ref = self.data[-1]	14✔
231
232	def gen_setup(	14✔
233	self,
234	prate_kw="rate",
235	rad_kw="rad",
236	r_ref_kw="r_ref",
237	h_ref_kw="h_ref",
238	dummy=False,
239	):
240	"""Generate the Spotpy Setup.
241
242	Parameters
243	----------
244	prate_kw : :class:`str`, optional
245	Keyword name for the pumping rate in the used type curve.
246	Default: "rate"
247	rad_kw : :class:`str`, optional
248	Keyword name for the radius in the used type curve.
249	Default: "rad"
250	r_ref_kw : :class:`str`, optional
251	Keyword name for the reference radius in the used type curve.
252	Default: "r_ref"
253	h_ref_kw : :class:`str`, optional
254	Keyword name for the reference head in the used type curve.
255	Default: "h_ref"
256	dummy : :class:`bool`, optional
257	Add a dummy parameter to the model. This could be used to equalize
258	sensitivity analysis.
259	Default: False
260	"""
261	self.extra_kw_names = {"rad": rad_kw}	14✔
262	setup_kw = dcopy(self.setup_kw) # create a copy here	14✔
263	setup_kw["val_fix"].setdefault(prate_kw, self.prate)	14✔
264	setup_kw["val_fix"].setdefault(rad_kw, self.rad)	14✔
265	setup_kw["val_fix"].setdefault(r_ref_kw, self.r_ref)	14✔
266	setup_kw["val_fix"].setdefault(h_ref_kw, self.h_ref)	14✔
267	setup_kw.setdefault("data", self.data)	14✔
268	setup_kw["dummy"] = dummy	14✔
269	self.setup = spotpylib.TypeCurve(**setup_kw)	14✔
270
271	def run(	14✔
272	self,
273	rep=5000,
274	parallel="seq",
275	run=True,
276	folder=None,
277	dbname=None,
278	traceplotname=None,
279	fittingplotname=None,
280	interactplotname=None,
281	estname=None,
282	plot_style="WTP",
283	):
284	"""Run the estimation.
285
286	Parameters
287	----------
288	rep : :class:`int`, optional
289	The number of repetitions within the SCEua algorithm in spotpy.
290	Default: ``5000``
291	parallel : :class:`str`, optional
292	State if the estimation should be run in parallel or not. Options:
293
294	* ``"seq"``: sequential on one CPU
295	* ``"mpi"``: use the mpi4py package
296
297	Default: ``"seq"``
298	run : :class:`bool`, optional
299	State if the estimation should be executed. Otherwise all plots
300	will be done with the previous results.
301	Default: ``True``
302	folder : :class:`str`, optional
303	Path to the output folder. If ``None`` the CWD is used.
304	Default: ``None``
305	dbname : :class:`str`, optional
306	File-name of the database of the spotpy estimation.
307	If ``None``, it will be the current time +
308	``"_db"``.
309	Default: ``None``
310	traceplotname : :class:`str`, optional
311	File-name of the parameter trace plot of the spotpy estimation.
312	If ``None``, it will be the current time +
313	``"_paratrace.pdf"``.
314	Default: ``None``
315	fittingplotname : :class:`str`, optional
316	File-name of the fitting plot of the estimation.
317	If ``None``, it will be the current time +
318	``"_fit.pdf"``.
319	Default: ``None``
320	interactplotname : :class:`str`, optional
321	File-name of the parameter interaction plot
322	of the spotpy estimation.
323	If ``None``, it will be the current time +
324	``"_parainteract.pdf"``.
325	Default: ``None``
326	estname : :class:`str`, optional
327	File-name of the results of the spotpy estimation.
328	If ``None``, it will be the current time +
329	``"_estimate"``.
330	Default: ``None``
331	plot_style : str, optional
332	Plot style. The default is "WTP".
333	"""
334	if self.setup.dummy:	14✔
335	raise ValueError(	×
336	"Estimate: for parameter estimation"
337	" you can't use a dummy parameter."
338	)
339	act_time = timemodule.strftime("%Y-%m-%d_%H-%M-%S")	14✔
340
341	# generate the filenames
342	if folder is None:	14✔
343	folder = os.path.join(os.getcwd(), self.name)	14✔
344	folder = os.path.abspath(folder)	14✔
345	if not os.path.exists(folder):	14✔
346	os.makedirs(folder)	14✔
347
348	if dbname is None:	14✔
349	dbname = os.path.join(folder, act_time + "_db")	14✔
350	elif not os.path.isabs(dbname):	×
351	dbname = os.path.join(folder, dbname)	×
352	if traceplotname is None:	14✔
353	traceplotname = os.path.join(folder, act_time + "_paratrace.pdf")	14✔
354	elif not os.path.isabs(traceplotname):	×
355	traceplotname = os.path.join(folder, traceplotname)	×
356	if fittingplotname is None:	14✔
357	fittingplotname = os.path.join(folder, act_time + "_fit.pdf")	14✔
358	elif not os.path.isabs(fittingplotname):	×
359	fittingplotname = os.path.join(folder, fittingplotname)	×
360	if interactplotname is None:	14✔
361	interactplotname = os.path.join(folder, act_time + "_interact.pdf")	14✔
362	elif not os.path.isabs(interactplotname):	×
363	interactplotname = os.path.join(folder, interactplotname)	×
364	if estname is None:	14✔
365	paraname = os.path.join(folder, act_time + "_estimate.txt")	14✔
366	elif not os.path.isabs(estname):	×
367	paraname = os.path.join(folder, estname)	×
368
369	# generate the parameter-names for plotting
370	paranames = dcopy(self.setup.para_names)	14✔
371	paralabels = []	14✔
372	for name in paranames:	14✔
373	p_label = self.setup.val_plot_names[name]	14✔
374	fit_n = self.setup.val_fit_name[name]	14✔
375	paralabels.append(f"{fit_n}({p_label})" if fit_n else p_label)	14✔
376
377	if parallel == "mpi":	14✔
378	# send the dbname of rank0
379	from mpi4py import MPI	×
380
381	comm = MPI.COMM_WORLD	×
382	rank = comm.Get_rank()	×
383	size = comm.Get_size()	×
384	if rank == 0:	×
385	print(rank, "send dbname:", dbname)	×
386	for i in range(1, size):	×
387	comm.send(dbname, dest=i, tag=0)	×
388	else:
389	dbname = comm.recv(source=0, tag=0)	×
390	print(rank, "got dbname:", dbname)	×
391	else:
392	rank = 0	14✔
393
394	# initialize the sampler
395	sampler = spotpy.algorithms.sceua(	14✔
396	self.setup,
397	dbname=dbname,
398	dbformat="csv",
399	parallel=parallel,
400	save_sim=True,
401	db_precision=np.float64,
402	)
403	# start the estimation with the sce-ua algorithm
404	if run:	14✔
405	sampler.sample(rep, ngs=10, kstop=100, pcento=1e-4, peps=1e-3)	14✔
406
407	if rank == 0:	14✔
408	if run:	14✔
409	self.result = sampler.getdata()	14✔
410	else:
411	self.result = np.genfromtxt(	×
412	dbname + ".csv", delimiter=",", names=True
413	)
414	para_opt = spotpy.analyser.get_best_parameterset(	14✔
415	self.result, maximize=False
416	)
417	void_names = para_opt.dtype.names	14✔
418	para = []	14✔
419	header = []	14✔
420	for name in void_names:	14✔
421	para.append(para_opt[0][name])	14✔
422	fit_n = self.setup.val_fit_name[name[3:]]	14✔
423	header.append(f"{fit_n}-{name[3:]}" if fit_n else name[3:])	14✔
424	self.estimated_para[name[3:]] = para[-1]	14✔
425	np.savetxt(paraname, para, header=" ".join(header))	14✔
426	# plot the estimation-results
427	plotter.plotparatrace(	14✔
428	result=self.result,
429	parameternames=paranames,
430	parameterlabels=paralabels,
431	stdvalues=self.estimated_para,
432	plotname=traceplotname,
433	style=plot_style,
434	)
435	plotter.plotfit_steady(	14✔
436	setup=self.setup,
437	data=self.data,
438	para=self.estimated_para,
439	rad=self.rad,
440	radnames=self.radnames,
441	extra=self.extra_kw_names,
442	plotname=fittingplotname,
443	style=plot_style,
444	)
445	plotter.plotparainteract(	14✔
446	self.result, paralabels, interactplotname, style=plot_style
447	)
448
449	def sensitivity(	14✔
450	self,
451	rep=None,
452	parallel="seq",
453	folder=None,
454	dbname=None,
455	plotname=None,
456	traceplotname=None,
457	sensname=None,
458	plot_style="WTP",
459	):
460	"""Run the sensitivity analysis.
461
462	Parameters
463	----------
464	rep : :class:`int`, optional
465	The number of repetitions within the FAST algorithm in spotpy.
466	Default: estimated
467	parallel : :class:`str`, optional
468	State if the estimation should be run in parallel or not. Options:
469
470	* ``"seq"``: sequential on one CPU
471	* ``"mpi"``: use the mpi4py package
472
473	Default: ``"seq"``
474	folder : :class:`str`, optional
475	Path to the output folder. If ``None`` the CWD is used.
476	Default: ``None``
477	dbname : :class:`str`, optional
478	File-name of the database of the spotpy estimation.
479	If ``None``, it will be the current time +
480	``"_sensitivity_db"``.
481	Default: ``None``
482	plotname : :class:`str`, optional
483	File-name of the result plot of the sensitivity analysis.
484	If ``None``, it will be the current time +
485	``"_sensitivity.pdf"``.
486	Default: ``None``
487	traceplotname : :class:`str`, optional
488	File-name of the parameter trace plot of the spotpy sensitivity
489	analysis.
490	If ``None``, it will be the current time +
491	``"_senstrace.pdf"``.
492	Default: ``None``
493	sensname : :class:`str`, optional
494	File-name of the results of the FAST estimation.
495	If ``None``, it will be the current time +
496	``"_estimate"``.
497	Default: ``None``
498	plot_style : str, optional
499	Plot style. The default is "WTP".
500	"""
501	if len(self.setup.para_names) == 1 and not self.setup.dummy:	14✔
502	raise ValueError(	×
503	"Sensitivity: for estimation with only one parameter"
504	" you have to use a dummy parameter."
505	)
506	if rep is None:	14✔
507	rep = spotpylib.fast_rep(	14✔
508	len(self.setup.para_names) + int(self.setup.dummy)
509	)
510
511	act_time = timemodule.strftime("%Y-%m-%d_%H-%M-%S")	14✔
512	# generate the filenames
513	if folder is None:	14✔
514	folder = os.path.join(os.getcwd(), self.name)	14✔
515	folder = os.path.abspath(folder)	14✔
516	if not os.path.exists(folder):	14✔
517	os.makedirs(folder)	×
518
519	if dbname is None:	14✔
520	dbname = os.path.join(folder, act_time + "_sensitivity_db")	14✔
521	elif not os.path.isabs(dbname):	×
522	dbname = os.path.join(folder, dbname)	×
523	if plotname is None:	14✔
524	plotname = os.path.join(folder, act_time + "_sensitivity.pdf")	14✔
525	elif not os.path.isabs(plotname):	×
526	plotname = os.path.join(folder, plotname)	×
527	if traceplotname is None:	14✔
528	traceplotname = os.path.join(folder, act_time + "_senstrace.pdf")	14✔
529	elif not os.path.isabs(traceplotname):	×
530	traceplotname = os.path.join(folder, traceplotname)	×
531	if sensname is None:	14✔
532	sensname = os.path.join(folder, act_time + "_FAST_estimate.txt")	14✔
533	elif not os.path.isabs(sensname):	×
534	sensname = os.path.join(folder, sensname)	×
535
536	sens_base, sens_ext = os.path.splitext(sensname)	14✔
537	sensname1 = sens_base + "_S1" + sens_ext	14✔
538
539	# generate the parameter-names for plotting
540	paranames = dcopy(self.setup.para_names)	14✔
541	paralabels = []	14✔
542	for name in paranames:	14✔
543	p_label = self.setup.val_plot_names[name]	14✔
544	fit_n = self.setup.val_fit_name[name]	14✔
545	paralabels.append(f"{fit_n}({p_label})" if fit_n else p_label)	14✔
546
547	if self.setup.dummy:	14✔
548	paranames.append("dummy")	14✔
549	paralabels.append("dummy")	14✔
550
551	if parallel == "mpi":	14✔
552	# send the dbname of rank0
553	from mpi4py import MPI	×
554
555	comm = MPI.COMM_WORLD	×
556	rank = comm.Get_rank()	×
557	size = comm.Get_size()	×
558	if rank == 0:	×
559	print(rank, "send dbname:", dbname)	×
560	for i in range(1, size):	×
561	comm.send(dbname, dest=i, tag=0)	×
562	else:
563	dbname = comm.recv(source=0, tag=0)	×
564	print(rank, "got dbname:", dbname)	×
565	else:
566	rank = 0	14✔
567
568	# initialize the sampler
569	sampler = spotpy.algorithms.fast(	14✔
570	self.setup,
571	dbname=dbname,
572	dbformat="csv",
573	parallel=parallel,
574	save_sim=True,
575	db_precision=np.float64,
576	)
577	sampler.sample(rep)	14✔
578
579	if rank == 0:	14✔
580	data = sampler.getdata()	14✔
581	parmin = sampler.parameter()["minbound"]	14✔
582	parmax = sampler.parameter()["maxbound"]	14✔
583	bounds = list(zip(parmin, parmax))	14✔
584	sens_est = sampler.analyze(	14✔
585	bounds, np.nan_to_num(data["like1"]), len(paranames), paranames
586	)
587	self.sens = {}	14✔
588	for sen_typ in sens_est:	14✔
589	self.sens[sen_typ] = {	14✔
590	par: sen for par, sen in zip(paranames, sens_est[sen_typ])
591	}
592	header = " ".join(paranames)	14✔
593	np.savetxt(sensname, sens_est["ST"], header=header)	14✔
594	np.savetxt(sensname1, sens_est["S1"], header=header)	14✔
595	plotter.plotsensitivity(	14✔
596	paralabels, sens_est, plotname, style=plot_style
597	)
598	plotter.plotparatrace(	14✔
599	data,
600	parameternames=paranames,
601	parameterlabels=paralabels,
602	stdvalues=None,
603	plotname=traceplotname,
604	style=plot_style,
605	)

GeoStat-Framework / welltestpy / 10691143420

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous