script:
    "../scripts/00.calculate_theta_watterson.py"

script:
    "../scripts/01.create_abc_simulation_randints.py"

script:
    "../scripts/01.create_abc_pod_randints.py"

script:
    "../scripts/01.simulate_treeseqs.py"

script:
    "../scripts/01.confirm_params_for_loci_being_same_and_save_params.py"

script:
    "../scripts/01.simulate_treeseq_pods.py"

script:
    "../scripts/01.generate_discretizing_breakpoints_for_sumstats.py"

script:
    "../scripts/01.calculate_sumstats.py"

script:
    "../scripts/01.calculate_masked_sumstats.py"

script:
    "../scripts/01.aggregate_sumstats.py"

script:
    "../scripts/01.aggregate_sumstats.py"

script:
    "../scripts/01.create_abc_simulation_randints.py"

script:
    "../scripts/01.alternative_model_simulate_treeseqs.py"

script:
    "../scripts/01.calculate_sumstats.py"

script:
    "../scripts/01.calculate_masked_sumstats.py"

script:
    "../scripts/01.confirm_params_for_loci_being_same_and_save_params.py"

script:
    "../scripts/01.aggregate_sumstats.py"

script:
    "../scripts/01.aggregate_sumstats.py"

script:
    "../scripts/01.calculate_podstats.py"

script:
    "../scripts/01.aggregate_podstats.py"

script:
    "../scripts/01.calculate_masked_podstats.py"

script:
    "../scripts/01.aggregate_podstats.py"

script:
    "../scripts/02.model_choice.R"

script:
    "../scripts/02.model_choice.R"

script:
    "../scripts/02.aggregate_model_choice.R"

script:
    "../scripts/02.aggregate_model_choice.R"

script:
    "../scripts/03.parameter_estimation.R"

script:
    "../scripts/03.parameter_estimation.R"

script:
    "../scripts/03.aggregate_parameter_estimation.R"

script:
    "../scripts/03.aggregate_parameter_estimation.R"

script:
    "../scripts/01.create_abc_simulation_randints.py"

script:
    "../scripts/04.simulate_sixparmodel_treeseqs.py"

script:
    "../scripts/01.confirm_params_for_loci_being_same_and_save_params.py"

script:
    "../scripts/01.calculate_sumstats.py"

script:
    "../scripts/01.calculate_masked_sumstats.py"

script:
    "../scripts/01.aggregate_sumstats.py"

script:
    "../scripts/01.aggregate_sumstats.py"

script:
    "../scripts/04.calculate_athal_observations.py"

script:
    "../scripts/04.calculate_athal_observations.py"

script:
    "../scripts/04.calculate_athal_all_region_observations.py"

script:
    "../scripts/04.calculate_athal_all_region_observations.py"

script:
    "../scripts/04.aggregate_athal_sumstats.py"

script:
    "../scripts/04.aggregate_athal_sumstats.py"

script:
    "../scripts/04.aggregate_athal_sumstats.py"

script:
    "../scripts/04.aggregate_athal_sumstats.py"

script:
    "../scripts/04.parameter_estimation.R"

script:
    "../scripts/04.parameter_estimation.R"

script:
    "../scripts/04.parameter_estimation.R"

script:
    "../scripts/04.parameter_estimation.R"

script:
    "../scripts/04.aggregate_athal_parameter_estimation.R"

script:
    "../scripts/04.aggregate_athal_parameter_estimation.R"

run:
    sys.exit(
        "#" * 600 + "inside aggregate_athal_parameter_estimation_masked\n" + ""
    )

run:
    sys.exit(
        "#" * 600 + "inside aggregate_athal_parameter_estimation_masked\n" + ""
    )

script:
    "../scripts/05.visualize_podstats.R"

script:
    "../scripts/05.visualize_podstats.R"

script:
    "../scripts/05.visualize_model_choice.R"

script:
    "../scripts/05.visualize_model_choice.R"

script:
    "../scripts/05.visualize_parameter_estimation.R"

script:
    "../scripts/05.visualize_parameter_estimation.R"

import numpy as np
import pandas as pd
import re
import datetime

# log, create log file
with open(snakemake.log.log1, "w", encoding="utf-8") as logfile:
    print(datetime.datetime.now(), end="\t", file=logfile)
    print(
        "start aggregating params and sumstats",
        file=logfile,
    )


# obtain simid and locid from filenames
nloci = len(snakemake.params.locid_wc)
nsimid = len(snakemake.params.simid_wc)

# create data list to read and store by locus
read_data_list = [[[] for _ in range(nloci)] for _ in range(nsimid)]
for npyid, npy in enumerate(snakemake.input.npys):
    split = re.split(r"_|\.|\/", npy)
    simid = int(split[np.where(np.array(split) == "sim")[0].max() + 1])
    locid = int(split[np.where(np.array(split) == "locus")[0].max() + 1])
    read_data_list[np.where(simid == np.array(snakemake.params.simid_wc))[0][0]][
        locid
    ] = np.load(npy)
read_data = np.array(read_data_list, dtype=float)


# take mean value over the independent loci
read_data = read_data.mean(axis=1)


# reshape the data into a 2d-array, because each input file provides a cluster
# batch of independent simulations
read_data = read_data.reshape(
    (read_data.shape[0] * read_data.shape[1], read_data.shape[2])
)

with open(snakemake.log.log1, "a", encoding="utf-8") as logfile:
    print(datetime.datetime.now(), end="\t", file=logfile)
    print("worked out sumstats", file=logfile)

# read parameters
read_param_list = [[] for _ in range(nsimid)]
for npy in snakemake.input.params:
    split = re.split(r"_|\.|\/", npy)
    simid = int(split[np.where(np.array(split) == "sim")[0].max() + 1])
    read_param_list[
        np.where(simid == np.array(snakemake.params.simid_wc))[0][0]
    ] = np.load(npy)
read_param = np.array(read_param_list)
num_provided_params = read_param.shape[2]
read_param = read_param.reshape(
    (read_param.shape[0] * read_param.shape[1], read_param.shape[2])
)

with open(snakemake.log.log1, "a", encoding="utf-8") as logfile:
    print(datetime.datetime.now(), end="\t", file=logfile)
    print("worked out params", file=logfile)

# concatenate parameters and sumstatas
result_table = np.concatenate((read_param, read_data), axis=1)

with open(snakemake.log.log1, "a", encoding="utf-8") as logfile:
    print(datetime.datetime.now(), end="\t", file=logfile)
    print("concatenated out params and sumstats", file=logfile)


# read sumstat names
sumstat_names_list = []
for name_file in snakemake.input.sumstat_names:
    sumstat_names_list.append(np.load(name_file))
    if len(sumstat_names_list) >= 2:
        assert all(
            sumstat_names_list[-2] == sumstat_names_list[-1]
        ), "sumstats are not congruent"
sumstat_names = sumstat_names_list[0]

# add param names
column_names = np.concatenate(
    (
        np.array([f"param_{parid}" for parid in range(num_provided_params)]),
        sumstat_names,
    ),
    axis=0,
)

result_table = pd.DataFrame(data=result_table, columns=column_names)
result_table.to_feather(snakemake.output.sumstats, compression="lz4")

with open(snakemake.log.log1, "a", encoding="utf-8") as logfile:
    print(datetime.datetime.now(), end="\t", file=logfile)
    print(f"saved dataframe to {snakemake.output.sumstats}", file=logfile)

import sys
import datetime
import pickle
import numpy as np
import tskit
import msprime
import pyfuncs  # from file


# log, create log file
with open(snakemake.log.log1, "w", encoding="utf-8") as logfile:
    print(datetime.datetime.now(), end="\t", file=logfile)
    print("start simulating tree sequences", file=logfile)


def main(simid):
    """Heart of this script

    The main function will create a single tree sequence. The main function has
    to be repeatedly executed until the list of treesequences has been created.
    """
    # provide the correct seed to create the random number generator
    seed = np.load(snakemake.input.npy)[simid][int(float(snakemake.wildcards.locid))]
    seed_params = np.load(snakemake.input.npy)[simid][
        0
    ]  # the parameters for the independent loci must be the same
    rng = np.random.default_rng(seed)
    rng_params = np.random.default_rng(seed_params)

    # draw the parameters of the model
    params = [
        pyfuncs.draw_parameter_from_prior(prior_definition, rng_params)
        for prior_definition in snakemake.params.model["priors"]
    ]
    del (
        seed_params,
        rng_params,
    )  # only the parameter drawing relies on the same seed, the simulations must be independent

    # simulate
    ts = pyfuncs.simulate_treesequence_under_alternative_model(
        params, snakemake.params, rng, snakemake.log.log1
    )

    return ts, params


# define the loop for the clustering of consequential simulations
start = int(float(snakemake.wildcards.simid))
end = min(
    start + int(float(snakemake.config["ABC"]["jobluster_simulations"])),
    int(float(snakemake.config["ABC"]["alternative_model"]["nsim"])),
)

assert start < end, "simid (Simulation IDs) maldefined, you will not simulate anything"

# run the simulations
treesequence_list = []
parameters_list = []
for simid in range(start, end):
    # log
    with open(snakemake.log.log1, "a", encoding="utf-8") as logfile:
        print(datetime.datetime.now(), end="\t", file=logfile)
        print(f"creating treesequence no {simid}", file=logfile)

    treeseq, params = main(simid)
    treesequence_list.append(treeseq)
    parameters_list.append(params)


# save treeseq list
with open(snakemake.output.tsl, "wb") as tsl_file:
    pickle.dump(treesequence_list, tsl_file)


# save parameters list
np.save(snakemake.output.params, np.array(parameters_list))

import sys
import datetime
import warnings
import itertools
import pickle
import re
import math
import numpy as np
import tskit
import pyfuncs  # from file


# log, create log file
with open(snakemake.log.log1, "w", encoding="utf-8") as logfile:
    print(datetime.datetime.now(), end="\t", file=logfile)
    print("start calculating summarizing statistics", file=logfile)


# Loading list of tree sequences
with open(snakemake.input.tsl, "rb") as tsl_file:
    tsl = np.array(pickle.load(tsl_file))


# sample only one haplotype per individual
rng = np.random.default_rng(snakemake.params.seed)
tsl_haploid = np.empty(tsl.shape, dtype=tskit.TreeSequence)
for treeid, treeseq in np.ndenumerate(tsl):
    sample_set = [
        this_individual.nodes[rng.integers(low=0, high=2)]
        for this_individual in treeseq.individuals()
    ]
    tsl_haploid[treeid] = treeseq.simplify(samples=sample_set)
del tsl
tsl = tsl_haploid
del tsl_haploid


# read and prepare mask files
mask = []
for mask_file in snakemake.input.mask:
    # the id of the chromosome from input file
    split = re.split(r"_|\.|\/", mask_file)
    chromid = int(split[np.where(np.array(split) == "locus")[0].max() + 1])

    # load mask file
    this_mask = np.loadtxt(mask_file, dtype=int)
    region_start, region_end = snakemake.config["ABC"]["athaliana"]["observations"][
        "treeseq_1001"
    ]["chosen_region"][chromid]
    region_mask = this_mask[
        (region_start <= this_mask[:, 0]) & (region_end > this_mask[:, 1])
    ]
    this_mask = region_mask - region_start

    # filter mask for disjoint intervals
    this_mask = pyfuncs.filter_mask_for_disjoint_intervals(
        this_mask, log=snakemake.log.log1
    )
    mask.append(this_mask)
    del region_mask, this_mask


# log
with open(snakemake.log.log1, "a", encoding="utf-8") as logfile:
    print(datetime.datetime.now(), end="\t", file=logfile)
    print("loaded all masks for pods", file=logfile)


# log
with open(snakemake.log.log1, "a", encoding="utf-8") as logfile:
    print(datetime.datetime.now(), end="\t", file=logfile)
    print("read and prepared the mask file", file=logfile)


tsl_masked = np.empty(tsl.shape, dtype=tskit.TreeSequence)
for treeid, treeseq in np.ndenumerate(tsl):
    chromid = treeid[1]
    tsl_masked[treeid] = treeseq.delete_intervals(
        mask[chromid], simplify=True, record_provenance=True
    )

del tsl
tsl = tsl_masked
del tsl_masked


# log
with open(snakemake.log.log1, "a", encoding="utf-8") as logfile:
    print(datetime.datetime.now(), end="\t", file=logfile)
    print("read and masked the tree sequences", file=logfile)


# read the breakpoint files
for breaks_filename in snakemake.input.breakpoints:
    if "SFS" in breaks_filename:
        with warnings.catch_warnings():  # will prevent warning if file is empty
            warnings.simplefilter("ignore")
            breaks_sfs = np.loadtxt(breaks_filename, dtype=float)
    elif "LD" in breaks_filename:
        with warnings.catch_warnings():
            warnings.simplefilter("ignore")
            breaks_ld = np.loadtxt(breaks_filename, dtype=float)
    elif "TM_WIN" in breaks_filename:
        with warnings.catch_warnings():
            warnings.simplefilter("ignore")
            breaks_tm_win = np.loadtxt(breaks_filename, dtype=float)

# log
with open(snakemake.log.log1, "a", encoding="utf-8") as logfile:
    print(datetime.datetime.now(), end="\t", file=logfile)
    print("read the breaks for binning", file=logfile)


# which sumstats to calculate
listed_sumstats = set(
    itertools.chain.from_iterable(
        [sumstat_set.split("/") for sumstat_set in snakemake.params.sumstats]
    )
)


# log
with open(snakemake.log.log1, "a", encoding="utf-8") as logfile:
    print(datetime.datetime.now(), end="\t", file=logfile)
    print(
        "calculating summarizing statistics",
        file=logfile,
    )

# calculate the SFS
if "SFS" in listed_sumstats:
    breaks = breaks_sfs

    dataframe_sumstats_sfs = np.empty(tsl.shape, dtype=tskit.TreeSequence)
    for treeid, treeseq in np.ndenumerate(tsl):
        dataframe_sumstats_sfs[treeid] = treeseq.allele_frequency_spectrum(
            sample_sets=None,
            windows=None,
            mode="site",
            span_normalise=False,
            polarised=snakemake.config["ABC"]["sumstats_specs"]["SFS"]["polarised"],
        )

    # check if the sfs shall be discretized as well; this usually is not
    # necessary as the sfs is a discrete statistic already
    if len(breaks):
        # for each sfs sum up the values in between the bin_edges (=breaks)
        assert (
            len(breaks) >= 2
        ), "expect closed interval for bin_edges; at least 2 breaks must be provided"
        new_sfs_dataframe = []
        for sfs in dataframe_sumstats_sfs:
            sfs_vals = []
            for breakid in range(1, len(breaks)):
                low = int(breaks[breakid - 1])
                high = breaks[breakid]
                if np.isinf(high):
                    high = len(sfs) + 1
                else:
                    high = math.ceil(high)
                sfs_vals.append(sum(sfs[low:high]))
            new_sfs_dataframe.append(sfs_vals)

        # reassign to the dataframe sumstat
        dataframe_sumstats_sfs = np.array(new_sfs_dataframe)

    # get average of sfs over the different loci
    dataframe_sumstats_sfs = dataframe_sumstats_sfs.mean(axis=1)

    dataframe_sumstats_sfs = np.concatenate(dataframe_sumstats_sfs, axis=0).reshape(
        (len(dataframe_sumstats_sfs), dataframe_sumstats_sfs[0].shape[0])
    )
else:
    dataframe_sumstats_sfs = np.empty(shape=(0, 0))


# log
with open(snakemake.log.log1, "a", encoding="utf-8") as logfile:
    print(datetime.datetime.now(), end="\t", file=logfile)
    print(
        "calculated SFS statistics",
        file=logfile,
    )


# calculate the LD
if "LD" in listed_sumstats:
    specs = snakemake.config["ABC"]["sumstats_specs"]["LD"]
    breaks = breaks_ld

    dataframe_sumstats_ld = np.empty(tsl.shape, dtype=tskit.TreeSequence)
    for treeid, treeseq in np.ndenumerate(tsl):
        dataframe_sumstats_ld[treeid] = pyfuncs.calculate_ld(
            treeseq, specs, breaks, rng, snakemake.log.log1
        )

    # get average of ld over the different loci
    dataframe_sumstats_ld = dataframe_sumstats_ld.mean(axis=1)

    dataframe_sumstats_ld = np.concatenate(dataframe_sumstats_ld, axis=0).reshape(
        (len(dataframe_sumstats_ld), dataframe_sumstats_ld[0].shape[0])
    )
else:
    dataframe_sumstats_ld = np.empty(shape=(0, 0))


# log
with open(snakemake.log.log1, "a", encoding="utf-8") as logfile:
    print(datetime.datetime.now(), end="\t", file=logfile)
    print(
        "calculated LD statistics",
        file=logfile,
    )


# calculate the TM_WIN
if "TM_WIN" in listed_sumstats:
    specs = snakemake.config["ABC"]["sumstats_specs"]["TM_WIN"]
    breaks = breaks_tm_win

    dataframe_sumstats_tm_win = np.empty(tsl.shape, dtype=tskit.TreeSequence)
    for treeid, treeseq in np.ndenumerate(tsl):
        dataframe_sumstats_tm_win[treeid] = pyfuncs.calculate_tm_win(
            treeseq, specs, breaks_tm_win
        )

    # get average of tm_win over the different loci
    dataframe_sumstats_tm_win = dataframe_sumstats_tm_win.mean(axis=1)
    dataframe_sumstats_tm_win = np.concatenate(
        dataframe_sumstats_tm_win, axis=0
    ).reshape((len(dataframe_sumstats_tm_win), dataframe_sumstats_tm_win[0].shape[0]))
else:
    dataframe_sumstats_tm_win = np.empty(shape=(0, 0))


# log
with open(snakemake.log.log1, "a", encoding="utf-8") as logfile:
    print(datetime.datetime.now(), end="\t", file=logfile)
    print(
        "calculated TM_WIN statistics",
        file=logfile,
    )


# fuse the summarizing stats list into a 2d-np.array, each row containing the
# summarizing stats of a single simulated tree sequence
dataframe_sumstats = np.concatenate(
    [
        npy
        for npy in (
            dataframe_sumstats_sfs,
            dataframe_sumstats_ld,
            dataframe_sumstats_tm_win,
        )
        if npy.size != 0
    ],
    axis=1,
)


# log
with open(snakemake.log.log1, "a", encoding="utf-8") as logfile:
    print(datetime.datetime.now(), end="\t", file=logfile)
    print("calculated summarizing statisitcs", file=logfile)


# save treeseq list
np.save(snakemake.output.npy, dataframe_sumstats)


# log
with open(snakemake.log.log1, "a", encoding="utf-8") as logfile:
    print(datetime.datetime.now(), end="\t", file=logfile)
    print(f"saved sumstats to {snakemake.output.npy}", file=logfile)


# save the names of the sumstats
sumstat_type = (
    ["sfs"] * dataframe_sumstats_sfs.shape[1]
    + ["ld"] * dataframe_sumstats_ld.shape[1]
    + ["tm_win"] * dataframe_sumstats_tm_win.shape[1]
)
sumstat_type_id = [
    str(element)
    for element in list(range(dataframe_sumstats_sfs.shape[1]))
    + list(range(dataframe_sumstats_ld.shape[1]))
    + list(range(dataframe_sumstats_tm_win.shape[1]))
]

sumstat_names = np.array(list(map("_".join, zip(sumstat_type, sumstat_type_id))))


assert (
    len(sumstat_names) == dataframe_sumstats.shape[1]
), "name vector for summary statistics has the wrong dimension"


np.save(snakemake.output.sumstat_count, sumstat_names)

import sys
import datetime
import warnings
import itertools
import pickle
import math
import numpy as np
import tskit
import pyfuncs  # from file


# log, create log file
with open(snakemake.log.log1, "w", encoding="utf-8") as logfile:
    print(datetime.datetime.now(), end="\t", file=logfile)
    print("start calculating summarizing statisticss", file=logfile)


# Loading list of tree sequences
with open(snakemake.input.tsl, "rb") as tsl_file:
    tsl = pickle.load(tsl_file)

# sample only one haplotype per individual
rng = np.random.default_rng(snakemake.params.seed)
tsl_haploid = []
for treeseq in tsl:
    sample_set = [
        this_individual.nodes[rng.integers(low=0, high=2)]
        for this_individual in treeseq.individuals()
    ]
    tsl_haploid.append(treeseq.simplify(samples=sample_set))
del tsl
tsl = tsl_haploid
del tsl_haploid


# read and prepare mask files
mask = np.loadtxt(snakemake.input.mask, dtype=int)
region_start, region_end = snakemake.config["ABC"]["athaliana"]["observations"][
    "treeseq_1001"
]["chosen_region"][int(float(snakemake.wildcards.locid))]
region_mask = mask[(region_start <= mask[:, 0]) & (region_end > mask[:, 1])]
mask = region_mask - region_start
del region_mask

# filter mask for disjoint intervals
mask = pyfuncs.filter_mask_for_disjoint_intervals(mask, log=snakemake.log.log1)


# log
with open(snakemake.log.log1, "a", encoding="utf-8") as logfile:
    print(datetime.datetime.now(), end="\t", file=logfile)
    print(
        "proportion of exons to mask in chromosome(1-indexed) "
        + f"{int(float(snakemake.wildcards.locid)) + 1}: "
        + f"{(mask[:, 1] - mask[:, 0]).sum()/(region_end-region_start)}",
        file=logfile,
    )


# log
with open(snakemake.log.log1, "a", encoding="utf-8") as logfile:
    print(datetime.datetime.now(), end="\t", file=logfile)
    print("read and prepared the mask file", file=logfile)


# mask for regions provided
tsl_masked = []
for treeseq in tsl:
    tsl_masked.append(
        treeseq.delete_intervals(mask, simplify=True, record_provenance=True)
    )
tsl = tsl_masked
del tsl_masked


# log
with open(snakemake.log.log1, "a", encoding="utf-8") as logfile:
    print(datetime.datetime.now(), end="\t", file=logfile)
    print("read and masked the tree sequences", file=logfile)


# read the breakpoint files
for breaks_filename in snakemake.input.breakpoints:
    if "SFS" in breaks_filename:
        with warnings.catch_warnings():  # will prevent warning if file is empty
            warnings.simplefilter("ignore")
            breaks_sfs = np.loadtxt(breaks_filename, dtype=float)
    elif "LD" in breaks_filename:
        with warnings.catch_warnings():
            warnings.simplefilter("ignore")
            breaks_ld = np.loadtxt(breaks_filename, dtype=float)
    elif "TM_WIN" in breaks_filename:
        with warnings.catch_warnings():
            warnings.simplefilter("ignore")
            breaks_tm_win = np.loadtxt(breaks_filename, dtype=float)

# log
with open(snakemake.log.log1, "a", encoding="utf-8") as logfile:
    print(datetime.datetime.now(), end="\t", file=logfile)
    print("read the breaks for binning", file=logfile)


# which sumstats to calculate
listed_sumstats = set(
    itertools.chain.from_iterable(
        [sumstat_set.split("/") for sumstat_set in snakemake.params.sumstats]
    )
)


# log
with open(snakemake.log.log1, "a", encoding="utf-8") as logfile:
    print(datetime.datetime.now(), end="\t", file=logfile)
    print(
        "calculating summarizing statistics",
        file=logfile,
    )

# calculate the SFS
if "SFS" in listed_sumstats:
    breaks = breaks_sfs
    dataframe_sumstats_sfs = np.array(
        [
            treeseq.allele_frequency_spectrum(
                sample_sets=None,
                windows=None,
                mode="site",
                span_normalise=False,
                polarised=snakemake.config["ABC"]["sumstats_specs"]["SFS"]["polarised"],
            )
            for treeseq in tsl
        ]
    )

    # check if the sfs shall be discretized as well; this usually is not
    # necessary as the sfs is a discrete statistic already
    if len(breaks):
        # for each sfs sum up the values in between the bin_edges (=breaks)
        assert (
            len(breaks) >= 2
        ), "expect closed interval for bin_edges; at least 2 breaks must be provided"
        new_sfs_dataframe = []
        for sfs in dataframe_sumstats_sfs:
            sfs_vals = []
            for breakid in range(1, len(breaks)):
                low = int(breaks[breakid - 1])
                high = breaks[breakid]
                if np.isinf(high):
                    high = len(sfs) + 1
                else:
                    high = math.ceil(high)
                sfs_vals.append(sum(sfs[low:high]))
            new_sfs_dataframe.append(sfs_vals)

        # reassign to the dataframe sumstat
        dataframe_sumstats_sfs = np.array(new_sfs_dataframe)
else:
    dataframe_sumstats_sfs = np.empty(shape=(0, 0))


# log
with open(snakemake.log.log1, "a", encoding="utf-8") as logfile:
    print(datetime.datetime.now(), end="\t", file=logfile)
    print(
        "calculated SFS statistics",
        file=logfile,
    )


# calculate the LD
if "LD" in listed_sumstats:
    specs = snakemake.config["ABC"]["sumstats_specs"]["LD"]
    breaks = breaks_ld

    dataframe_sumstats_ld = np.empty(np.array(tsl).shape, dtype=tskit.TreeSequence)
    for treeid, treeseq in np.ndenumerate(tsl):
        dataframe_sumstats_ld[treeid] = pyfuncs.calculate_ld(
            treeseq, specs, breaks, rng, snakemake.log.log1
        )

    dataframe_sumstats_ld = np.concatenate(dataframe_sumstats_ld, axis=0).reshape(
        (len(dataframe_sumstats_ld), dataframe_sumstats_ld[0].shape[0])
    )
else:
    dataframe_sumstats_ld = np.empty(shape=(0, 0))


# log
with open(snakemake.log.log1, "a", encoding="utf-8") as logfile:
    print(datetime.datetime.now(), end="\t", file=logfile)
    print(
        "calculated LD statistics",
        file=logfile,
    )


# calculate the TM_WIN
if "TM_WIN" in listed_sumstats:
    specs = snakemake.config["ABC"]["sumstats_specs"]["TM_WIN"]
    breaks = breaks_tm_win
    dataframe_sumstats_tm_win = np.array(
        [pyfuncs.calculate_tm_win(treeseq, specs, breaks_tm_win) for treeseq in tsl]
    )
else:
    dataframe_sumstats_tm_win = np.empty(shape=(0, 0))


# log
with open(snakemake.log.log1, "a", encoding="utf-8") as logfile:
    print(datetime.datetime.now(), end="\t", file=logfile)
    print(
        "calculated TM_WIN statistics",
        file=logfile,
    )


# fuse the summarizing stats list into a 2d-np.array, each row containing the
# summarizing stats of a single simulated tree sequence
dataframe_sumstats = np.concatenate(
    [
        npy
        for npy in (
            dataframe_sumstats_sfs,
            dataframe_sumstats_ld,
            dataframe_sumstats_tm_win,
        )
        if npy.size != 0
    ],
    axis=1,
)


# log
with open(snakemake.log.log1, "a", encoding="utf-8") as logfile:
    print(datetime.datetime.now(), end="\t", file=logfile)
    print("calculated summarizing statisitcs", file=logfile)


# save treeseq list
np.save(snakemake.output.npy, dataframe_sumstats)


# log
with open(snakemake.log.log1, "a", encoding="utf-8") as logfile:
    print(datetime.datetime.now(), end="\t", file=logfile)
    print(f"saved sumstats to {snakemake.output.npy}", file=logfile)


# save the names of the sumstats
sumstat_type = (
    ["sfs"] * dataframe_sumstats_sfs.shape[1]
    + ["ld"] * dataframe_sumstats_ld.shape[1]
    + ["tm_win"] * dataframe_sumstats_tm_win.shape[1]
)
sumstat_type_id = [
    str(element)
    for element in list(range(dataframe_sumstats_sfs.shape[1]))
    + list(range(dataframe_sumstats_ld.shape[1]))
    + list(range(dataframe_sumstats_tm_win.shape[1]))
]

sumstat_names = np.array(list(map("_".join, zip(sumstat_type, sumstat_type_id))))


assert (
    len(sumstat_names) == dataframe_sumstats.shape[1]
), "name vector for summary statistics has the wrong dimension"


np.save(snakemake.output.sumstat_count, sumstat_names)