class TanimotoScoringFunction(ScoringFunctionBasedOnRdkitMol):
    """
    Scoring function that looks at the fingerprint similarity against a target molecule.
    """

    def __init__(self, target, fp_type, score_modifier: ScoreModifier = None) -> None:
        """
        Args:
            target: target molecule
            fp_type: fingerprint type
            score_modifier: score modifier
        """
        super().__init__(score_modifier=score_modifier)

        self.target = target
        self.fp_type = fp_type
        target_mol = smiles_to_rdkit_mol(target)
        if target_mol is None:
            raise RuntimeError(f'The similarity target {target} is not a valid molecule.')

        self.ref_fp = get_fingerprint(target_mol, self.fp_type)

    def score_mol(self, mol: Chem.Mol) -> float:
        fp = get_fingerprint(mol, self.fp_type)
        return TanimotoSimilarity(fp, self.ref_fp) class RdkitScoringFunction(ScoringFunctionBasedOnRdkitMol):
    """
    Scoring function wrapping RDKit descriptors.
    """

    def __init__(self, descriptor: Callable[[Chem.Mol], float], score_modifier: ScoreModifier = None) -> None:
        """
        Args:
            descriptor: molecular descriptors, such as the ones in descriptors.py
            score_modifier: score modifier
        """
        super().__init__(score_modifier=score_modifier)
        self.descriptor = descriptor

    def score_mol(self, mol: Chem.Mol) -> float:
        return self.descriptor(mol) class CNS_MPO_ScoringFunction(ScoringFunctionBasedOnRdkitMol):
    """
    CNS MPO scoring function
    """

    def __init__(self, max_logP=5.0, maxMW=360, min_tpsa=40, max_tpsa=90, max_hbd=0) -> None:
        super().__init__()

        self.logP_gauss = MinGaussianModifier(max_logP, 1)
        self.molW_gauss = MinGaussianModifier(maxMW, 60)
        self.tpsa_maxgauss = MaxGaussianModifier(min_tpsa, 20)
        self.tpsa_mingauss = MinGaussianModifier(max_tpsa, 30)
        self.hbd_gauss = MinGaussianModifier(max_hbd, 2.0)

    def score_mol(self, mol: Chem.Mol) -> float:
        mw = mol_weight(mol)
        lp = logP(mol)
        hbd = num_H_donors(mol)
        mol_tpsa = tpsa(mol)

        o1 = self.tpsa_mingauss(mol_tpsa)
        o2 = self.tpsa_maxgauss(mol_tpsa)
        o3 = self.hbd_gauss(hbd)
        o4 = self.logP_gauss(lp)
        o5 = self.molW_gauss(mw)

        return 0.2 * (o1 + o2 + o3 + o4 + o5) class IsomerScoringFunction(MoleculewiseScoringFunction):
    """
    Scoring function for closeness to a molecular formula.

    The score penalizes deviations from the required number of atoms for each element type, and for the total
    number of atoms.

    F.i., if the target formula is C2H4, the scoring function is the average of three contributions:
    - number of C atoms with a Gaussian modifier with mu=2, sigma=1
    - number of H atoms with a Gaussian modifier with mu=4, sigma=1
    - total number of atoms with a Gaussian modifier with mu=6, sigma=2
    """

    def __init__(self, molecular_formula: str, mean_function='geometric') -> None:
        """
        Args:
            molecular_formula: target molecular formula
            mean_function: which function to use for averaging: 'arithmetic' or 'geometric'
        """
        super().__init__()

        self.mean_function = self.determine_mean_function(mean_function)
        self.scoring_functions = self.determine_scoring_functions(molecular_formula)

    @staticmethod
    def determine_mean_function(mean_function: str) -> Callable[[List[float]], float]:
        if mean_function == 'arithmetic':
            return arithmetic_mean
        if mean_function == 'geometric':
            return geometric_mean
        raise ValueError(f'Invalid mean function: "{mean_function}"')

    @staticmethod
    def determine_scoring_functions(molecular_formula: str) -> List[RdkitScoringFunction]:
        element_occurrences = parse_molecular_formula(molecular_formula)

        total_number_atoms = sum(element_tuple[1] for element_tuple in element_occurrences)

        # scoring functions for each element
        functions = [RdkitScoringFunction(descriptor=AtomCounter(element),
                                          score_modifier=GaussianModifier(mu=n_atoms, sigma=1.0))
                     for element, n_atoms in element_occurrences]

        # scoring functions for the total number of atoms
        functions.append(RdkitScoringFunction(descriptor=num_atoms,
                                              score_modifier=GaussianModifier(mu=total_number_atoms, sigma=2.0)))

        return functions

    def raw_score(self, smiles: str) -> float:
        # return the average of all scoring functions
        scores = [f.score(smiles) for f in self.scoring_functions]
        if self.corrupt_score in scores:
            return self.corrupt_score
        return self.mean_function(scores) class SMARTSScoringFunction(ScoringFunctionBasedOnRdkitMol):
    """
    Tests for SMARTS which should be or should not be present in the compound.


    """

    def __init__(self, target: str, inverse=False) -> None:
        """

        :param target: The SMARTS string to match.
        :param inverse: Specifies whether the SMARTS is desired (False) or an antipattern, which we don't want to see
                        in the molecules (inverse=False)
        """
        super().__init__()
        self.inverse = inverse
        self.smarts = target
        self.target = Chem.MolFromSmarts(target)

        assert target is not None

    def score_mol(self, mol: Chem.Mol) -> float:

        matches = mol.GetSubstructMatches(self.target)

        if len(matches) > 0:
            if self.inverse:
                return 0.0
            else:
                return 1.0
        else:
            if self.inverse:
                return 1.0
            else:
                return 0.0 class DistributionLearningBenchmark:
    """
    Base class for assessing how well a model is able to generate molecules matching a molecule distribution.

    Derived class should implement the assess_molecules function
    """

    def __init__(self, name: str, number_samples: int) -> None:
        self.name = name
        self.number_samples = number_samples

    @abstractmethod
    def assess_model(self, model: DistributionMatchingGenerator) -> DistributionLearningBenchmarkResult:
        """
        Assess a distribution-matching generator model.

        Args:
            model: model to assess
        """ class NoveltyBenchmark(DistributionLearningBenchmark):
    def __init__(self, number_samples: int, training_set: Iterable[str]) -> None:
        """
        Args:
            number_samples: number of samples to generate from the model
            training_set: molecules from the training set
        """
        super().__init__(name='Novelty', number_samples=number_samples)
        self.training_set_molecules = set(canonicalize_list(training_set, include_stereocenters=False))

    def assess_model(self, model: DistributionMatchingGenerator) -> DistributionLearningBenchmarkResult:
        """
        Assess a distribution-matching generator model.

        Args:
            model: model to assess
        """
        start_time = time.time()
        molecules = sample_unique_molecules(model=model, number_molecules=self.number_samples, max_tries=2)
        end_time = time.time()

        if len(molecules) != self.number_samples:
            logger.warning('The model could not generate enough unique molecules. The score will be penalized.')

        # canonicalize_list in order to remove stereo information (also removes duplicates and invalid molecules, but there shouldn't be any)
        unique_molecules = set(canonicalize_list(molecules, include_stereocenters=False))

        novel_molecules = unique_molecules.difference(self.training_set_molecules)

        novel_ratio = len(novel_molecules) / self.number_samples

        metadata = {
            'number_samples': self.number_samples,
            'number_novel': len(novel_molecules)
        }

        return DistributionLearningBenchmarkResult(benchmark_name=self.name,
                                                   score=novel_ratio,
                                                   sampling_time=end_time - start_time,
                                                   metadata=metadata) class KLDivBenchmark(DistributionLearningBenchmark):
    """
    Computes the KL divergence between a number of samples and the training set for physchem descriptors
    """

    def __init__(self, number_samples: int, training_set: List[str]) -> None:
        """
        Args:
            number_samples: number of samples to generate from the model
            training_set: molecules from the training set
        """
        super().__init__(name='KL divergence', number_samples=number_samples)
        self.training_set_molecules = canonicalize_list(get_random_subset(training_set, self.number_samples, seed=42),
                                                        include_stereocenters=False)
        self.pc_descriptor_subset = [
            'BertzCT',
            'MolLogP',
            'MolWt',
            'TPSA',
            'NumHAcceptors',
            'NumHDonors',
            'NumRotatableBonds',
            'NumAliphaticRings',
            'NumAromaticRings'
        ]

    def assess_model(self, model: DistributionMatchingGenerator) -> DistributionLearningBenchmarkResult:
        """
        Assess a distribution-matching generator model.

        Args:
            model: model to assess
        """
        start_time = time.time()
        molecules = sample_unique_molecules(model=model, number_molecules=self.number_samples, max_tries=2)
        end_time = time.time()

        if len(molecules) != self.number_samples:
            logger.warning('The model could not generate enough unique molecules. The score will be penalized.')

        # canonicalize_list in order to remove stereo information (also removes duplicates and invalid molecules, but there shouldn't be any)
        unique_molecules = set(canonicalize_list(molecules, include_stereocenters=False))

        # first we calculate the descriptors, which are np.arrays of size n_samples x n_descriptors
        d_sampled = calculate_pc_descriptors(unique_molecules, self.pc_descriptor_subset)
        d_chembl = calculate_pc_descriptors(self.training_set_molecules, self.pc_descriptor_subset)

        kldivs = {}

        # now we calculate the kl divergence for the float valued descriptors ...
        for i in range(4):
            kldiv = continuous_kldiv(X_baseline=d_chembl[:, i], X_sampled=d_sampled[:, i])
            kldivs[self.pc_descriptor_subset[i]] = kldiv

        # ... and for the int valued ones.
        for i in range(4, 9):
            kldiv = discrete_kldiv(X_baseline=d_chembl[:, i], X_sampled=d_sampled[:, i])
            kldivs[self.pc_descriptor_subset[i]] = kldiv

        # pairwise similarity

        chembl_sim = calculate_internal_pairwise_similarities(self.training_set_molecules)
        chembl_sim = chembl_sim.max(axis=1)

        sampled_sim = calculate_internal_pairwise_similarities(unique_molecules)
        sampled_sim = sampled_sim.max(axis=1)

        kldiv_int_int = continuous_kldiv(X_baseline=chembl_sim, X_sampled=sampled_sim)
        kldivs['internal_similarity'] = kldiv_int_int

        # for some reason, this runs into problems when both sets are identical.
        # cross_set_sim = calculate_pairwise_similarities(self.training_set_molecules, unique_molecules)
        # cross_set_sim = cross_set_sim.max(axis=1)
        #
        # kldiv_ext = discrete_kldiv(chembl_sim, cross_set_sim)
        # kldivs['external_similarity'] = kldiv_ext
        # kldiv_sum += kldiv_ext

        metadata = {
            'number_samples': self.number_samples,
            'kl_divs': kldivs
        }

        # Each KL divergence value is transformed to be in [0, 1].
        # Then their average delivers the final score.
        partial_scores = [np.exp(-score) for score in kldivs.values()]
        score = sum(partial_scores) / len(partial_scores)

        return DistributionLearningBenchmarkResult(benchmark_name=self.name,
                                                   score=score,
                                                   sampling_time=end_time - start_time,
                                                   metadata=metadata) class FrechetBenchmark(DistributionLearningBenchmark):
    """
    Calculates the Fréchet ChemNet Distance.

    See http://dx.doi.org/10.1021/acs.jcim.8b00234 for the publication.
    """

    def __init__(self, training_set: List[str],
                 chemnet_model_filename='ChemNet_v0.13_pretrained.h5',
                 sample_size=10000) -> None:
        """
        Args:
            training_set: molecules from the training set
            chemnet_model_filename: name of the file for trained ChemNet model.
                Must be present in the 'fcd' package, since it will be loaded directly from there.
            sample_size: how many molecules to generate the distribution statistics from (both reference data and model)
        """
        self.chemnet_model_filename = chemnet_model_filename
        self.sample_size = sample_size
        super().__init__(name='Frechet ChemNet Distance', number_samples=self.sample_size)

        self.reference_molecules = get_random_subset(training_set, self.sample_size, seed=42)

    def assess_model(self, model: DistributionMatchingGenerator) -> DistributionLearningBenchmarkResult:
        chemnet = self._load_chemnet()

        start_time = time.time()
        generated_molecules = sample_valid_molecules(model=model, number_molecules=self.number_samples)
        end_time = time.time()

        if len(generated_molecules) != self.number_samples:
            logger.warning('The model could not generate enough valid molecules.')

        mu_ref, cov_ref = self._calculate_distribution_statistics(chemnet, self.reference_molecules)
        mu, cov = self._calculate_distribution_statistics(chemnet, generated_molecules)

        FCD = fcd.calculate_frechet_distance(mu1=mu_ref, mu2=mu,
                                             sigma1=cov_ref, sigma2=cov)
        score = np.exp(-0.2 * FCD)

        metadata = {
            'number_reference_molecules': len(self.reference_molecules),
            'number_generated_molecules': len(generated_molecules),
            'FCD': FCD
        }

        return DistributionLearningBenchmarkResult(benchmark_name=self.name,
                                                   score=score,
                                                   sampling_time=end_time - start_time,
                                                   metadata=metadata)

    def _load_chemnet(self):
        """
        Load the ChemNet model from the file specified in the init function.

        This file lives inside a package but to use it, it must always be an actual file.
        The safest way to proceed is therefore:
        1. read the file with pkgutil
        2. save it to a temporary file
        3. load the model from the temporary file
        """
        model_bytes = pkgutil.get_data('fcd', self.chemnet_model_filename)
        assert model_bytes is not None

        tmpdir = tempfile.gettempdir()
        model_path = os.path.join(tmpdir, self.chemnet_model_filename)

        with open(model_path, 'wb') as f:
            f.write(model_bytes)

        logger.info(f'Saved ChemNet model to \'{model_path}\'')

        return fcd.load_ref_model(model_path)

    def _calculate_distribution_statistics(self, model, molecules: List[str]):
        sample_std = fcd.canonical_smiles(molecules)
        gen_mol_act = fcd.get_predictions(model, sample_std)

        mu = np.mean(gen_mol_act, axis=0)
        cov = np.cov(gen_mol_act.T)
        return mu, cov class GoalDirectedBenchmark:
    """
    This class assesses how well a model is able to generate molecules satisfying a given objective.
    """

    def __init__(self, name: str, objective: ScoringFunction,
                 contribution_specification: ScoreContributionSpecification,
                 starting_population: Optional[List[str]] = None) -> None:
        """
        Args:
            name: Benchmark name
            objective: Objective for the goal-directed optimization
            contribution_specification: Specifies how to calculate the global benchmark score
        """
        self.name = name
        self.objective = objective
        self.wrapped_objective = ScoringFunctionWrapper(scoring_function=objective)
        self.contribution_specification = contribution_specification
        self.starting_population = starting_population

    def assess_model(self, model: GoalDirectedGenerator) -> GoalDirectedBenchmarkResult:
        """
        Assess the given model by asking it to generate molecules optimizing a scoring function.
        The number of molecules to generate is determined automatically from the score contribution specification.

        Args:
            model: model to assess
        """
        number_molecules_to_generate = max(self.contribution_specification.top_counts)
        start_time = time.time()
        molecules = model.generate_optimized_molecules(scoring_function=self.wrapped_objective,
                                                       number_molecules=number_molecules_to_generate,
                                                       starting_population=self.starting_population
                                                       )
        end_time = time.time()

        canonicalized_molecules = canonicalize_list(molecules, include_stereocenters=False)
        unique_molecules = remove_duplicates(canonicalized_molecules)
        scores = self.objective.score_list(unique_molecules)

        if len(unique_molecules) != number_molecules_to_generate:
            number_missing = number_molecules_to_generate - len(unique_molecules)
            logger.warning(f'An incorrect number of distinct molecules was generated: '
                           f'{len(unique_molecules)} instead of {number_molecules_to_generate}. '
                           f'Padding scores with {number_missing} zeros...')
            scores.extend([0.0] * number_missing)

        global_score, top_x_dict = compute_global_score(self.contribution_specification, scores)

        scored_molecules = zip(unique_molecules, scores)
        sorted_scored_molecules = sorted(scored_molecules, key=lambda x: (x[1], x[0]), reverse=True)

        internal_similarities = calculate_internal_pairwise_similarities(unique_molecules)

        # accumulate internal_similarities in metadata
        int_simi_histogram = np.histogram(internal_similarities, bins=10, range=(0, 1), density=True)

        metadata: Dict[str, Any] = {}
        metadata.update(top_x_dict)
        metadata['internal_similarity_max'] = internal_similarities.max()
        metadata['internal_similarity_mean'] = internal_similarities.mean()
        metadata["internal_similarity_histogram_density"] = int_simi_histogram[0].tolist(),
        metadata["internal_similarity_histogram_bins"] = int_simi_histogram[1].tolist(),

        return GoalDirectedBenchmarkResult(benchmark_name=self.name,
                                           score=global_score,
                                           optimized_molecules=sorted_scored_molecules,
                                           execution_time=end_time - start_time,
                                           number_scoring_function_calls=self.wrapped_objective.evaluations,
                                           metadata=metadata) def uniform_specification(*top_counts: int) -> ScoreContributionSpecification:
    """
    Creates an instance of ScoreContributionSpecification where all the top-x contributions have equal weight

    Args:
        top_counts: list of values, where each value x will correspond to the top-x contribution
    """
    contributions = [(x, 1.0) for x in top_counts]
    return ScoreContributionSpecification(contributions=contributions) class ScoreModifier:
class ChainedModifier(ScoreModifier):
class LinearModifier(ScoreModifier):
class SquaredModifier(ScoreModifier):
class AbsoluteScoreModifier(ScoreModifier):
class GaussianModifier(ScoreModifier):
class MinMaxGaussianModifier(ScoreModifier):
class ClippedScoreModifier(ScoreModifier):
class SmoothClippedScoreModifier(ScoreModifier):
class ThresholdedLinearModifier(ScoreModifier):
    def __call__(self, x):
    def __init__(self, modifiers: List[ScoreModifier]) -> None:
    def __call__(self, x):
    def __init__(self, slope=1.0):
    def __call__(self, x):
    def __init__(self, target_value: float, coefficient=1.0) -> None:
    def __call__(self, x):
    def __init__(self, target_value: float) -> None:
    def __call__(self, x):
    def __init__(self, mu: float, sigma: float) -> None:
    def __call__(self, x):
    def __init__(self, mu: float, sigma: float, minimize=False) -> None:
    def __call__(self, x):
    def __init__(self, upper_x: float, lower_x=0.0, high_score=1.0, low_score=0.0) -> None:
    def __call__(self, x):
    def __init__(self, upper_x: float, lower_x=0.0, high_score=1.0, low_score=0.0) -> None:
    def __call__(self, x):
    def __init__(self, threshold: float) -> None:
    def __call__(self, x): class ArithmeticMeanScoringFunction(BatchScoringFunction):
    """
    Scoring function that combines multiple scoring functions linearly.
    """

    def __init__(self, scoring_functions: List[ScoringFunction], weights=None) -> None:
        """
        Args:
            scoring_functions: scoring functions to combine
            weights: weight for the corresponding scoring functions. If None, all will have the same weight.
        """
        super().__init__()

        self.scoring_functions = scoring_functions
        number_scoring_functions = len(scoring_functions)

        self.weights = np.ones(number_scoring_functions) if weights is None else weights
        assert number_scoring_functions == len(self.weights)

    def raw_score_list(self, smiles_list: List[str]) -> List[float]:
        scores = []

        for function, weight in zip(self.scoring_functions, self.weights):
            res = function.score_list(smiles_list)
            scores.append(weight * np.array(res))

        scores = np.array(scores).sum(axis=0) / np.sum(self.weights)

        return list(scores) class GeometricMeanScoringFunction(MoleculewiseScoringFunction):
    """
    Scoring function that combines multiple scoring functions multiplicatively.
    """

    def __init__(self, scoring_functions: List[ScoringFunction]) -> None:
        """
        Args:
            scoring_functions: scoring functions to combine
        """
        super().__init__()

        self.scoring_functions = scoring_functions

    def raw_score(self, smiles: str) -> float:
        partial_scores = [f.score(smiles) for f in self.scoring_functions]
        if self.corrupt_score in partial_scores:
            return self.corrupt_score

        return geometric_mean(partial_scores) class ScoringFunction:
    """
    Base class for an objective function.

    In general, do not inherit directly from this class. Prefer `MoleculewiseScoringFunction` or `BatchScoringFunction`.
    """

    def __init__(self, score_modifier: ScoreModifier = None) -> None:
        """
        Args:
            score_modifier: Modifier to apply to the score. If None, will be LinearModifier()
        """
        self.score_modifier = score_modifier
        self.corrupt_score = -1.0

    @property
    def score_modifier(self):
        return self._score_modifier

    @score_modifier.setter
    def score_modifier(self, modifier: Optional[ScoreModifier]):
        self._score_modifier = LinearModifier() if modifier is None else modifier

    def modify_score(self, raw_score: float) -> float:
        return self._score_modifier(raw_score)

    @abstractmethod
    def score(self, smiles: str) -> float:
        """
        Score a single molecule as smiles
        """
        raise NotImplementedError

    @abstractmethod
    def score_list(self, smiles_list: List[str]) -> List[float]:
        """
        Score a list of smiles.

        Args:
            smiles_list: list of smiles [smiles1, smiles2,...]

        Returns: a list of scores

        the order of the input smiles is matched in the output.

        """
        raise NotImplementedError def num_rotatable_bonds(mol: Mol) -> int:
    return Descriptors.NumRotatableBonds(mol) def num_aromatic_rings(mol: Mol) -> int:
    return rdMolDescriptors.CalcNumAromaticRings(mol) def logP(mol: Mol) -> float:
    return Descriptors.MolLogP(mol) def qed(mol: Mol) -> float:
    return Descriptors.qed(mol) def tpsa(mol: Mol) -> float:
    return Descriptors.TPSA(mol) def bertz(mol: Mol) -> float:
    return Descriptors.BertzCT(mol) def mol_weight(mol: Mol) -> float:
    return Descriptors.MolWt(mol) class AtomCounter:

    def __init__(self, element: str) -> None:
        """
        Args:
            element: element to count within a molecule
        """
        self.element = element

    def __call__(self, mol: Mol) -> int:
        """
        Count the number of atoms of a given type.

        Args:
            mol: molecule

        Returns:
            The number of atoms of the given type.
        """
        # if the molecule contains H atoms, they may be implicit, so add them
        if self.element == 'H':
            mol = Chem.AddHs(mol)

        return sum(1 for a in mol.GetAtoms() if a.GetSymbol() == self.element) def num_rings(mol: Mol) -> int:
    return rdMolDescriptors.CalcNumRings(mol)
from rdkit import Chem
from guacamol.common_scoring_functions import TanimotoScoringFunction, RdkitScoringFunction, CNS_MPO_ScoringFunction, \
    IsomerScoringFunction, SMARTSScoringFunction
from guacamol.distribution_learning_benchmark import DistributionLearningBenchmark, NoveltyBenchmark, KLDivBenchmark
from guacamol.frechet_benchmark import FrechetBenchmark
from guacamol.goal_directed_benchmark import GoalDirectedBenchmark
from guacamol.goal_directed_score_contributions import uniform_specification
from guacamol.score_modifier import MinGaussianModifier, MaxGaussianModifier, ClippedScoreModifier, GaussianModifier
from guacamol.scoring_function import ArithmeticMeanScoringFunction, GeometricMeanScoringFunction, ScoringFunction
from guacamol.utils.descriptors import num_rotatable_bonds, num_aromatic_rings, logP, qed, tpsa, bertz, mol_weight, \
    AtomCounter, num_rings


def tpsa_benchmark(target: float) -> GoalDirectedBenchmark:
    benchmark_name = f'TPSA (target: {target})'
    objective = RdkitScoringFunction(descriptor=tpsa,
                                     score_modifier=GaussianModifier(mu=target, sigma=20.0))

    specification = uniform_specification(1, 10, 100)

    return GoalDirectedBenchmark(name=benchmark_name,
                                 objective=objective,
                                 contribution_specification=specification)


def cns_mpo(max_logP=5.0) -> GoalDirectedBenchmark:
    specification = uniform_specification(1, 10, 100)
    return GoalDirectedBenchmark(name='CNS MPO', objective=CNS_MPO_ScoringFunction(max_logP=max_logP),
                                 contribution_specification=specification)


def qed_benchmark() -> GoalDirectedBenchmark:
    specification = uniform_specification(1, 10, 100)
    return GoalDirectedBenchmark(name='QED',
                                 objective=RdkitScoringFunction(descriptor=qed),
                                 contribution_specification=specification)


def median_camphor_menthol(mean_cls=GeometricMeanScoringFunction) -> GoalDirectedBenchmark:
    t_camphor = TanimotoScoringFunction('CC1(C)C2CCC1(C)C(=O)C2', fp_type='ECFP4')
    t_menthol = TanimotoScoringFunction('CC(C)C1CCC(C)CC1O', fp_type='ECFP4')
    median = mean_cls([t_menthol, t_camphor])

    specification = uniform_specification(1, 10, 100)

    return GoalDirectedBenchmark(name='Median molecules 1',
                                 objective=median,
                                 contribution_specification=specification)


def novelty_benchmark(training_set_file: str, number_samples: int) -> DistributionLearningBenchmark:
    smiles_list = [s.strip() for s in open(training_set_file).readlines()]
    return NoveltyBenchmark(number_samples=number_samples, training_set=smiles_list)


def kldiv_benchmark(training_set_file: str, number_samples: int) -> DistributionLearningBenchmark:
    smiles_list = [s.strip() for s in open(training_set_file).readlines()]
    return KLDivBenchmark(number_samples=number_samples, training_set=smiles_list)


def frechet_benchmark(training_set_file: str, number_samples: int) -> DistributionLearningBenchmark:
    smiles_list = [s.strip() for s in open(training_set_file).readlines()]
    return FrechetBenchmark(training_set=smiles_list, sample_size=number_samples)


def perindopril_rings() -> GoalDirectedBenchmark:
    # perindopril with two aromatic rings
    perindopril = TanimotoScoringFunction('O=C(OCC)C(NC(C(=O)N1C(C(=O)O)CC2CCCCC12)C)CCC',
                                          fp_type='ECFP4')
    arom_rings = RdkitScoringFunction(descriptor=num_aromatic_rings,
                                      score_modifier=GaussianModifier(mu=2, sigma=0.5))

    specification = uniform_specification(1, 10, 100)

    return GoalDirectedBenchmark(name='Perindopril MPO',
                                 objective=GeometricMeanScoringFunction([perindopril, arom_rings]),
                                 contribution_specification=specification)


def amlodipine_rings() -> GoalDirectedBenchmark:
    # amlodipine with 3 rings
    amlodipine = TanimotoScoringFunction(r'Clc1ccccc1C2C(=C(/N/C(=C2/C(=O)OCC)COCCN)C)\C(=O)OC', fp_type='ECFP4')
    rings = RdkitScoringFunction(descriptor=num_rings,
                                 score_modifier=GaussianModifier(mu=3, sigma=0.5))

    specification = uniform_specification(1, 10, 100)

    return GoalDirectedBenchmark(name='Amlodipine MPO',
                                 objective=GeometricMeanScoringFunction([amlodipine, rings]),
                                 contribution_specification=specification)


def sitagliptin_replacement() -> GoalDirectedBenchmark:
    # Find a molecule dissimilar to sitagliptin, but with the same properties
    smiles = 'Fc1cc(c(F)cc1F)CC(N)CC(=O)N3Cc2nnc(n2CC3)C(F)(F)F'
    sitagliptin = Chem.MolFromSmiles(smiles)
    target_logp = logP(sitagliptin)
    target_tpsa = tpsa(sitagliptin)

    similarity = TanimotoScoringFunction(smiles, fp_type='ECFP4',
                                         score_modifier=GaussianModifier(mu=0, sigma=0.1))
    lp = RdkitScoringFunction(descriptor=logP,
                              score_modifier=GaussianModifier(mu=target_logp, sigma=0.2))
    tp = RdkitScoringFunction(descriptor=tpsa,
                              score_modifier=GaussianModifier(mu=target_tpsa, sigma=5))
    isomers = IsomerScoringFunction('C16H15F6N5O')

    specification = uniform_specification(1, 10, 100)

    return GoalDirectedBenchmark(name='Sitagliptin MPO',
                                 objective=GeometricMeanScoringFunction([similarity, lp, tp, isomers]),
                                 contribution_specification=specification)


def zaleplon_with_other_formula() -> GoalDirectedBenchmark:
    # zaleplon_with_other_formula with other formula
    zaleplon = TanimotoScoringFunction('O=C(C)N(CC)C1=CC=CC(C2=CC=NC3=C(C=NN23)C#N)=C1',
                                       fp_type='ECFP4')
    formula = IsomerScoringFunction('C19H17N3O2')

    specification = uniform_specification(1, 10, 100)

    return GoalDirectedBenchmark(name='Zaleplon MPO',
                                 objective=GeometricMeanScoringFunction([zaleplon, formula]),
                                 contribution_specification=specification)


def smarts_with_other_target(smarts: str, other_molecule: str) -> ScoringFunction:
    smarts_scoring_function = SMARTSScoringFunction(target=smarts)
    other_mol = Chem.MolFromSmiles(other_molecule)
    target_logp = logP(other_mol)

    target_tpsa = tpsa(other_mol)