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.61

/src/welltestpy/estimate/transient_lib.py

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

import anaflow as ana
import numpy as np
import spotpy

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

__all__ = [
    "TransientPumping",
]


class TransientPumping:
    """Class to estimate transient 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.
    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,
        val_fix=None,
        val_fit_type=None,
        val_fit_name=None,
        val_plot_names=None,
        testinclude=None,
        generate=False,
    ):
        val_fix = val_fix or {}
        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.time = None
        """:class:`numpy.ndarray`: time points of the observation"""
        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.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 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]
                )
                == "transient"
            ):
                raise ValueError(test + ": selection is not transient.")

        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.settime()
            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 settime(self, time=None, tmin=10.0, tmax=np.inf, typ="quad", steps=10):
        """Set uniform time points for the observations.

        Parameters
        ----------
        time : :class:`numpy.ndarray`, optional
            Array of specified time points. If ``None`` is given, they will
            be determined by the observation data.
            Default: ``None``
        tmin : :class:`float`, optional
            Minimal time value. It will set a minimal value of 10s.
            Default: ``10``
        tmax : :class:`float`, optional
            Maximal time value.
            Default: ``inf``
        typ : :class:`str` or :class:`float`, optional
            Typ of the time selection. You can select from:

                * ``"exp"``: for exponential behavior
                * ``"log"``: for logarithmic behavior
                * ``"geo"``: for geometric behavior
                * ``"lin"``: for linear behavior
                * ``"quad"``: for quadratic behavior
                * ``"cub"``: for cubic behavior
                * :class:`float`: here you can specifi any exponent
                  ("quad" would be equivalent to 2)

            Default: "quad"

        steps : :class:`int`, optional
            Number of generated time steps. Default: 10
        """
        if time is None:
            for test in self.testinclude:
                for obs in self.testinclude[test]:
                    _, temptime = self.campaign.tests[test].observations[obs]()
                    tmin = max(tmin, temptime.min())
                    tmax = min(tmax, temptime.max())
                    tmin = tmax if tmin > tmax else tmin
            time = ana.specialrange(tmin, tmax, steps, typ)

        for test in self.testinclude:
            for obs in self.testinclude[test]:
                processlib.filterdrawdown(
                    self.campaign.tests[test].observations[obs], tout=time
                )

        self.time = time

    def gen_data(self):
        """Generate the observed drawdown at given time points.

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

        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]()
                temphead = np.array(temphead).reshape(-1)[np.newaxis].T

                if data is None:
                    data = dcopy(temphead)
                else:
                    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]

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

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