How to use the pastas.StressModel2 function in pastas

def test_rfunc(rfunc_name):
    if rfunc_name not in []:
        obs = read_csv("tests/data/obs.csv", index_col=0, parse_dates=True,
                       squeeze=True)
        rain = read_csv("tests/data/rain.csv", index_col=0, parse_dates=True,
                        squeeze=True)
        evap = read_csv("tests/data/evap.csv", index_col=0, parse_dates=True,
                        squeeze=True)
        # Create the time series model
        ml = ps.Model(obs, name="Test_Model")

        ## Create stress
        rfunc = getattr(ps.rfunc, rfunc_name)
        sm = ps.StressModel2(stress=[rain, evap], rfunc=rfunc, name='test_sm')
        ml.add_stressmodel(sm)

        # Solve the time series model
        ml.solve()

import pastas as ps
from pastas.stressmodels import WellModel

fname = 'data/MenyanthesTest.men'
meny = ps.read.MenyData(fname)

# Create the time series model
H = meny.H['Obsevation well']
ml = ps.Model(H['values'])

# Add precipitation
IN = meny.IN['Precipitation']['values']
IN.index = IN.index.round("D")
IN2 = meny.IN['Evaporation']['values']
IN2.index = IN2.index.round("D")
sm = ps.StressModel2([IN, IN2], ps.Gamma, 'Recharge')
ml.add_stressmodel(sm)

stresses = [meny.IN['Extraction 1']["values"],
            meny.IN['Extraction 2']["values"],
            meny.IN['Extraction 3']["values"]]

# Get distances from metadata
xo = meny.H["Obsevation well"]['xcoord']
yo = meny.H["Obsevation well"]['ycoord']
distances = []
for extr in ['Extraction 1', 'Extraction 2', 'Extraction 3']:
    xw = meny.IN[extr]["xcoord"]
    yw = meny.IN[extr]["ycoord"]
    distances.append(np.sqrt((xo-xw)**2 + (yo-yw)**2))

w = WellModel(stresses, ps.Hantush, distances=distances, name="Wells")

import pastas as ps

# Create a simple model taken from example.py
obs = ps.read_dino('data/B58C0698001_1.csv')
rain = ps.read_knmi('data/neerslaggeg_HEIBLOEM-L_967-2.txt', variables='RD')
evap = ps.read_knmi('data/etmgeg_380.txt', variables='EV24')

# Create a Pastas Project
mls = ps.Project(name="test_project")

mls.add_series(obs, "GWL", kind="oseries", metadata=dict())
mls.add_series(rain, name="Prec", kind="prec", metadata=dict())
mls.add_series(evap, name="Evap", kind="evap", metadata=dict())

ml = mls.add_model(oseries="GWL")
sm = ps.StressModel2([mls.stresses.loc["Prec", "series"],
                      mls.stresses.loc["Evap", "series"]],
                     ps.Exponential, name='recharge')
ml.add_stressmodel(sm)
n = ps.NoiseModel()
ml.add_noisemodel(n)
ml.solve(freq="D", warmup=1000, report=False)

mls.to_file("test_project.pas")

fname = 'data/MenyanthesTest.men'
meny = ps.read.MenyData(fname)

# Create the time series model
H = meny.H['Obsevation well']
ml = ps.Model(H['values'])

# Add precipitation
IN = meny.IN['Precipitation']['values']
IN.index = IN.index.round("D")
IN.name = 'Precipitation'
IN2 = meny.IN['Evaporation']['values']
IN2.index = IN2.index.round("D")
IN2.name = 'Evaporation'
sm = ps.StressModel2([IN, IN2], ps.Gamma, 'Recharge')
ml.add_stressmodel(sm)

#Add well extraction 1
# IN = meny.IN['Extraction 1']
# well = ps.TimeSeries(IN["values"], freq_original="M", settings="well")
# # extraction amount counts for the previous month
# sm = ps.StressModel(well, ps.Hantush, 'Extraction_1', up=False)

# Add well extraction 2
IN = meny.IN['Extraction 2']
well = ps.TimeSeries(IN["values"], freq_original="M", settings="well")
# extraction amount counts for the previous month
sm1 = ps.StressModel(well, ps.Hantush, 'Extraction_2', up=False)

# Add well extraction 3
IN = meny.IN['Extraction 3']

fname = 'data/B32D0136001_1.csv'
obs = ps.read_dino(fname)

# Create the time series model
ml = ps.Model(obs)

# read climate data
fname = 'data/KNMI_Bilt.txt'
RH = ps.read_knmi(fname, variables='RH')
EV24 = ps.read_knmi(fname, variables='EV24')
#rech = RH.series - EV24.series

# Create stress
#sm = ps.Recharge(RH, EV24, ps.Gamma, ps.Linear, name='recharge')
#sm = Recharge(RH, EV24, Gamma, Combination, name='recharge')
sm = ps.StressModel2([RH, EV24], ps.Gamma, name='recharge')
#sm = ps.StressModel(RH, ps.Gamma, name='precip')
#sm1 = ps.StressModel(EV24, ps.Gamma, name='evap')
ml.add_stressmodel(sm)
#ml.add_tseries(sm1)

# Add noise model
n = ps.NoiseModel()
ml.add_noisemodel(n)

# Solve
ml.solve(freq="7D")
ml.plot()