text:
import copy
import json
import uuid
import multiprocessing as mp
from datetime import datetime
from pathlib import Path

import numpy as np

from pymoo.optimize import minimize
from pymoo.util.display import MultiObjectiveDisplay
from pymoo.util.termination.default import MultiObjectiveDefaultTermination

from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split

from competitor.models.loader import load_env, load_model

from backport.problem_factory import ProblemFactory
from backport.utils import RepresentationTranslator
from backport.wrapper.actions_wrapper import ActionWrapper
from backport.wrapper.datasets_wrapper import DatasetWrapper
from sequential.adult_dependency_graph import get_adult_dependency_graph
from feature_cost_model.feature_dependency_graph import FeatureDependencyGraph

from cscf.CSCF import CSCF

from scipy.special import softmax

from sacred import Experiment

experiment_name = "Comparison to competitor"
ex = Experiment(
    experiment_name,
    interactive=False,
)

def write_instances_to_file(
    instances, dataset_name, target_class, identifier, prefix=""
):
    Path(f"experiment_data/test_instances/").mkdir(parents=True, exist_ok=True)
    now = datetime.now()
    file_path = (
        f"experiment_data/test_instances/{prefix}_{dataset_name}_n={len(instances)}.csv"
    )

    instances = np.array(instances)

    # check if already exists and output a warning if so in case the content is different than
    # what we want to write
    if Path(file_path).is_file():
        existing_instances = np.loadtxt(file_path, delimiter=",")
        assert np.array_equal(
            instances, existing_instances
        ), "Initial instances are not the same. Probably the random seed is not set or wrong"
        if not np.array_equal(instances, existing_instances):
            print(f"Existing array at {file_path} is not the same as the current one!")
            # add exp_id
            file_path = f"experiment_data/test_instances/{prefix}_UNEQUAL_{identifier}_{now.strftime('%d-%m-%Y')}_{dataset_name}_n={len(instances)}.csv"
    else:
        np.savetxt(file_path, instances, delimiter=",")
        print(f"Saved {len(instances)} to {file_path}")

#############
#  CONFIG   #
#############

@ex.config
def general_cfg():
    # General
    experiment_description = "Evaluation against competitor"
    notes = ""
    seed = 1
    # ! If the competitor model is used, we cannot use parallelization
    # ! as it seems to not work with tensorflow. I.e., only with RFC!
    parallelized = False

    # Specific
    n_experiment_instances = 100
    # This causes usage of a lot of memory
    save_optimal_population_trace = True

## EA mains
@ex.named_config
def comp_adult_problem_main_cfg():
    """
    Adult dataset experiments, for the original and main
    EA variant that I want to present
    """
    dataset_name = "adult"
    target_class = 1

    chosen_objectives = [
        "feature_tweaking_frequencies",
        "summed_costs",
        "gowers_distance",
    ]

@ex.named_config
def comp_german_problem_main_cfg():
    """
    German dataset experiments, for the original and main
    EA variant that I want to present
    """
    dataset_name = "german"
    target_class = 0

    chosen_objectives = [
        "feature_tweaking_frequencies",
        "summed_costs",
        "gowers_distance",
    ]

## EA alternatives
@ex.named_config
def comp_adult_problem_alternative_cfg():
    """
    Adult dataset experiments, for the alternative
    EA variant which is more similar to the competitor
    """
    dataset_name = "adult"
    target_class = 1

    chosen_objectives = [
        "feature_tweaking_frequencies",
        "summed_costs_discounted",
        "gowers_distance",
    ]

@ex.config
def optimization_cfg():
    # For the EA
    total_pop = 500
    # 1, because we use NDS, so it's changed to the number of ND solutions anyways
    n_elites_frac = 0.2
    # n_elites = int(n_elites_frac * total_pop)
    n_elites = 1  # doesn't make a difference since we use NDS and not crowding distance
    # offsprings_frac = 0.7
    offsprings_frac = 0.8
    n_offspring = int(offsprings_frac * total_pop)
    # mutants_frac = 0.1
    mutants_frac = 0.2
    n_mutants = int(mutants_frac * total_pop)
    bias = 0.7
    eliminate_duplicates = True

    # For the minimization
    n_generations = 150

@ex.config
def model_cfg():
    """
    Legacy now replaced
    """
    # comp-NN or rfc
    # model_type = "rfc"
    model_type = "comp-NN"

    # RFC
    n_estimators = 100

#############
#  SETUPS   #
#############

@ex.capture
def setup_experiment_instances(
    dataset, model, target_class, n_experiment_instances, seed
):

    # original represenation of the instances
    data_in_x = dataset.get_optimizer_data()
    _, labels = dataset.get_classifier_data()
    # ! Use the transformed, exact, instances to make the comparison fair
    # ! Now both methods use the exact same starting conditions
    data_in_z = dataset.encode_features(data_in_x)

    full_idx = np.arange(len(labels))

    # Create conditions for the instances we want to pick from, i.e. candidates
    is_opposite_originally = labels != target_class
    is_opposite_in_blackbox_z = (
        np.argmax(model.predict(data_in_z), axis=1) != target_class
    )

    # * Actually not needed anymore, but kept since it doesn't change anything
    is_opposite_in_blackbox_x = (
        np.argmax(model.predict(dataset.encode_features(data_in_x)), axis=1)
        != target_class
    )

    # Only use target class ones
    full_idx = full_idx[
        (is_opposite_originally)
        & (is_opposite_in_blackbox_z)
        & (is_opposite_in_blackbox_x)
    ]
    # Random order
    random_order = np.random.permutation(len(full_idx))[:n_experiment_instances]

    final_chosen_idx = full_idx[random_order]

    initial_instances_in_x = data_in_x[final_chosen_idx]
    initial_instances_in_z = data_in_z[final_chosen_idx]
    print(final_chosen_idx)
    print(
        np.argmax(
            model.predict(dataset.encode_features(initial_instances_in_x)), axis=1
        )
    )
    assert (
        np.argmax(
            model.predict(dataset.encode_features(initial_instances_in_x)), axis=1
        )
        != target_class
    ).all()
    print(np.argmax(model.predict(initial_instances_in_z), axis=1))
    assert (
        np.argmax(model.predict(initial_instances_in_z), axis=1) != target_class
    ).all()
    print(labels[final_chosen_idx])
    assert len(initial_instances_in_x) == n_experiment_instances
    assert len(initial_instances_in_z) == n_experiment_instances
    assert len(final_chosen_idx) == n_experiment_instances

    return initial_instances_in_x, initial_instances_in_z, final_chosen_idx

@ex.capture
def setup_boundaries_and_allowed_values(dataset, dataset_name):
    # if dataset_name == "adult":
    #     bounds_and_values = {
    #         # 0: np.unique(dataset.xs[:, 2]).tolist(),  # Edu unique values
    #         0: [-max(abs(dataset.xs[:, 2])), max(abs(dataset.xs[:, 2]))],  # Edu
    #         1: [-max(abs(dataset.xs[:, 9])), max(abs(dataset.xs[:, 9]))],  # Cap loss
    #         2: [-max(abs(dataset.xs[:, 10])), max(abs(dataset.xs[:, 10]))],  # Work hrs
    #         3: np.unique(dataset.xs[:, 4]).tolist(),  # Enlist
    #         4: [-max(abs(dataset.xs[:, 8])), max(abs(dataset.xs[:, 8]))],  # Cap gain
    #         5: [-max(abs(dataset.xs[:, 0])), max(abs(dataset.xs[:, 0]))],  # Years
    #     }
    # elif dataset_name == "german":
    #     bounds_and_values = {
    #         0: [-max(abs(dataset.xs[:, 12])), max(abs(dataset.xs[:, 12]))],  # years
    #         1: np.unique(dataset.xs[:, 19]).tolist(),
    #         2: [-max(abs(dataset.xs[:, 4])), max(abs(dataset.xs[:, 4]))],
    #         3: [-max(abs(dataset.xs[:, 1])), max(abs(dataset.xs[:, 1]))],
    #         4: [-max(abs(dataset.xs[:, 4])), max(abs(dataset.xs[:, 4]))],
    #         5: np.unique(dataset.xs[:, 9]).tolist(),
    #         6: np.unique(dataset.xs[:, 16]).tolist(),
    #     }

    # * As z-values
    if dataset_name == "adult":
        bounds_and_values = {
            # 0: np.unique(dataset.xs[:, 2]).tolist(),  # Edu unique values
            0: [-max(abs(dataset.zs[:, 8])), max(abs(dataset.zs[:, 8]))],  # Edu
            1: [-max(abs(dataset.zs[:, 44])), max(abs(dataset.zs[:, 44]))],  # Cap loss
            2: [-max(abs(dataset.zs[:, 45])), max(abs(dataset.zs[:, 45]))],  # Work hrs
            3: np.unique(dataset.xs[:, 4]).tolist(),  # Enlist
            4: [-max(abs(dataset.zs[:, 43])), max(abs(dataset.zs[:, 43]))],  # Cap gain
            5: [-max(abs(dataset.zs[:, 0])), max(abs(dataset.zs[:, 0]))],  # Years
        }
    elif dataset_name == "german":
        bounds_and_values = {
            0: [-max(abs(dataset.zs[:, 46])), max(abs(dataset.zs[:, 46]))],  # years
            1: np.unique(dataset.xs[:, 19]).tolist(),
            2: [-max(abs(dataset.zs[:, 21])), max(abs(dataset.zs[:, 21]))],
            3: [-max(abs(dataset.zs[:, 4])), max(abs(dataset.zs[:, 4]))],
            4: [-max(abs(dataset.zs[:, 21])), max(abs(dataset.zs[:, 21]))],
            5: np.unique(dataset.xs[:, 9]).tolist(),
            6: np.unique(dataset.xs[:, 16]).tolist(),
        }
    return bounds_and_values

@ex.capture
def setup_dependency_graph(features, dataset_name):
    if dataset_name == "adult":
        dependency_graph = get_adult_dependency_graph(features)
    elif dataset_name == "german":
        dependency_graph = FeatureDependencyGraph(None)
    return dependency_graph

@ex.capture
def setup_actions(x0, legacy_actions, legacy_features, features, dataset_name):
    translator = RepresentationTranslator(legacy_features)

    if dataset_name == "german":
        (
            waitYears,
            naturalize,
            chCreditAm,
            chLoanPeriod,
            adjLoanPeriod,
            getGuarantor,
            getUnskilledJob,
        ) = legacy_actions

        actions = {
            0: {
                "action_id": 0,
                "legacy_action": waitYears,
                "translator": translator,
                "feature_idx": features["age_in_years"],
                "feature_idx_in_z": 46,
                "initial_instance": x0,
                "is_categorical": False,
                "feat_name": "age_in_years",
                "key": 0,
            },
            1: {
                "action_id": 1,
                "legacy_action": naturalize,
                "translator": translator,
                "feature_idx": features["foreign_worker"],
                "feature_idx_in_z": 19,
                "initial_instance": x0,
                "is_categorical": True,
                "feat_name": "foreign_worker",
                "key": 1,
            },
            2: {
                "action_id": 2,
                "legacy_action": chCreditAm,
                "translator": translator,
                "feature_idx": features["credit_amount"],
                "feature_idx_in_z": 21,
                "initial_instance": x0,
                "is_categorical": False,
                "feat_name": "credit_amount",
                "key": 2,
            },
            3: {
                "action_id": 3,
                "legacy_action": chLoanPeriod,
                "translator": translator,
                "feature_idx": features["loan_duration"],
                "feature_idx_in_z": 4,
                "initial_instance": x0,
                "is_categorical": False,
                "feat_name": "loan_duration",
                "key": 3,
            },
            4: {
                "action_id": 4,
                "legacy_action": adjLoanPeriod,
                "translator": translator,
                "feature_idx": features["credit_amount"],
                "feature_idx_in_z": 21,
                "initial_instance": x0,
                "is_categorical": False,
                "feat_name": "credit_amount",
                "key": 4,
            },
            5: {
                "action_id": 5,
                "legacy_action": getGuarantor,
                "translator": translator,
                "feature_idx": features["other_debtors_guarantors"],
                "feature_idx_in_z": 9,
                "initial_instance": x0,
                "is_categorical": True,
                "feat_name": "other_debtors_guarantors",
                "key": 5,
            },
            6: {
                "action_id": 6,
                "legacy_action": getUnskilledJob,
                "translator": translator,
                "feature_idx": features["job"],
                "feature_idx_in_z": 16,
                "initial_instance": x0,
                "is_categorical": True,
                "feat_name": "job",
                "key": 6,
            },
        }
    elif dataset_name == "adult":
        addEdu, chWorkHrs, chCapLoss, enlist, incCapGain, waitYears = legacy_actions

        actions = {
            0: {
                "action_id": 0,
                "legacy_action": addEdu,
                "translator": translator,
                "feature_idx": features["Education num"],
                "feature_idx_in_z": 8,
                "initial_instance": x0,
                "is_categorical": False,
                "feat_name": "Education num",
                "key": 0,
            },
            1: {
                "action_id": 1,
                "legacy_action": chCapLoss,
                "translator": translator,
                "feature_idx": features["Capital Loss"],
                "feature_idx_in_z": 44,
                "initial_instance": x0,
                "is_categorical": False,
                "feat_name": "Capital Loss",
                "key": 1,
            },
            2: {
                "action_id": 2,
                "legacy_action": chWorkHrs,
                "translator": translator,
                "feature_idx": features["Hours/Week"],
                "feature_idx_in_z": 45,
                "initial_instance": x0,
                "is_categorical": False,
                "feat_name": "Hours/Week",
                "key": 2,
            },
            3: {
                "action_id": 3,
                "legacy_action": enlist,
                "translator": translator,
                "feature_idx": features["Occupation"],
                "feature_idx_in_z": 4,
                "initial_instance": x0,
                "is_categorical": True,
                "feat_name": "Occupation",
                "key": 3,
            },
            4: {
                "action_id": 4,
                "legacy_action": incCapGain,
                "translator": translator,
                "feature_idx": features["Capital Gain"],
                "feature_idx_in_z": 43,
                "initial_instance": x0,
                "is_categorical": False,
                "feat_name": "Capital Gain",
                "key": 4,
            },
            5: {
                "action_id": 5,
                "legacy_action": waitYears,
                "translator": translator,
                "feature_idx": features["Age"],
                "feature_idx_in_z": 0,
                "initial_instance": x0,
                "is_categorical": False,
                "feat_name": "Age",
                "key": 5,
            },
        }
    return actions

@ex.capture
def setup_problem(
    dataset,
    actions,
    x0,
    blackbox,
    blackbox_prob,
    bounds_and_values,
    features,
    G,
    target_class,
    chosen_objectives,
    seed,
):

    problem = ProblemFactory(
        copy.copy(dataset),
        copy.copy(actions),
        copy.copy(x0),
        target_class=target_class,
        blackbox_classifier=blackbox,
        blackbox_probabilitiy=blackbox_prob,
        bounds_and_values=copy.copy(bounds_and_values),
        feature=features,
        categorical_features=dataset.cat_features,
        dependency_graph=G,
        chosen_objectives=chosen_objectives,
        legacy_mode=True,  # important!
    )
    return problem

@ex.capture
def setup_dataset(dataset_name, seed):
    if dataset_name == "german":
        legacy_dataset, legacy_actions, legacy_features, _ = load_env(
            dataset_name, "data.npy"
        )
    elif dataset_name == "adult":
        legacy_dataset, legacy_actions, legacy_features, _ = load_env(
            dataset_name, "data.npy"
        )

    dataset = DatasetWrapper(dataset_name, legacy_dataset, legacy_features)

    return dataset, legacy_actions, legacy_features

@ex.capture
def setup_model(dataset, dataset_name, model_type, n_estimators, seed):
    X_train, y_train = dataset.get_classifier_data()
    X_test, y_test = dataset.get_test_data()

    if model_type == "rfc":
        from sklearn.ensemble import RandomForestClassifier

        model = RandomForestClassifier(n_estimators=n_estimators, random_state=seed)
        # Train it
        model.fit(X_train, y_train)
    elif model_type == "comp-NN":
        model = load_model(dataset_name, "model.h5")
    return model

@ex.capture
def setup_optimizer(n_elites, n_offspring, n_mutants, bias, eliminate_duplicates, seed):
    algorithm = CSCF(
        n_elites=n_elites,
        n_offsprings=n_offspring,
        n_mutants=n_mutants,
        bias=bias,
        eliminate_duplicates=eliminate_duplicates,
    )
    return algorithm

@ex.capture
def get_blackboxes(model, dataset, model_type):
    if model_type == "rfc":

        def blackbox(X):
            _X = X.copy()
            if _X.ndim == 1:
                _X = _X.reshape(1, -1)
            return model.predict(dataset.encode_features(_X))

        def blackbox_prob(X):
            _X = X.copy()
            if _X.ndim == 1:
                _X = _X.reshape(1, -1)
            return model.predict_proba(dataset.encode_features(_X))

    elif model_type == "comp-NN":

        def blackbox(X):
            _X = X.copy()
            if _X.ndim == 1:
                _X = _X.reshape(1, -1)
            return np.argmax(model.predict(dataset.encode_features(_X)), axis=1)

        def blackbox_prob(X):
            _X = X.copy()
            if _X.ndim == 1:
                _X = _X.reshape(1, -1)
            # * Use softmax to get the actual probabilities
            return softmax(model.predict(dataset.encode_features(_X)), axis=1)

    return blackbox, blackbox_prob

def run_optimization_problem(
    i,
    x0,
    legacy_actions,
    legacy_features,
    dataset,
    bounds_and_values,
    target_class,
    model_type,
    n_generations,
    seed,
    experiment_id,
    save_optimal_population_trace,
    features,
    model=None,
):
    if model is None:
        model = setup_model(dataset)
    blackbox, blackbox_prob = get_blackboxes(model, dataset)
    actions = setup_actions(x0, legacy_actions, legacy_features, features)
    dependency_graph = setup_dependency_graph(features)
    problem = setup_problem(
        dataset,
        actions,
        x0,
        blackbox,
        blackbox_prob,
        bounds_and_values,
        features,
        dependency_graph,
    )
    algorithm = setup_optimizer()
    # termination = MultiObjectiveDefaultTermination(
    #     x_tol=1e-8,
    #     cv_tol=1e-6,
    #     f_tol=0.0025,
    #     nth_gen=5,
    #     n_last=30,
    #     n_max_gen=n_generations,
    #     n_max_evals=100000,
    # )
    termination = ("n_gen", n_generations)

    print(f"Running experiment {i}")
    res = minimize(
        problem,
        algorithm,
        termination,
        seed=seed,
        display=MultiObjectiveDisplay(),
        save_history=save_optimal_population_trace,
        verbose=False,
    )

    now = datetime.now()
    # Create dir if not exists
    folder_name = (
        f"{now.strftime('%d-%m-%Y')} {experiment_id} {dataset.name}/Final Results"
    )
    Path(f"results/{folder_name}/").mkdir(parents=True, exist_ok=True)
    try:
        file_path = f"results/{folder_name}/{experiment_id}_{now.strftime('%d-%m-%Y')}_{dataset.name}_res_opt_{i}.json"
        info = {
            "experiment": experiment_name,
            "experiment_id": experiment_id,
            "date": now.strftime("%d/%m/%Y at %H:%M:%S"),
            "target_class": target_class,
            "model_type": model_type,
            "x0_id": i,
            "x0": x0.tolist(),
            "iterations": n_generations,
            "seed": seed,
            "dataset": dataset.name,
            "success": True,
            "runtime": res.exec_time,
            "n_solutions": len(res.X),
            "X": res.X.tolist(),
            "F": res.F.tolist(),
            "G": res.G.tolist(),
            "pheno": res.opt.get("pheno").tolist(),
            "tweaked_res": res.opt.get("tweaked_instances").tolist(),
            "summed_costs": res.opt.get("summed_costs").tolist(),
            "summed_costs_discounted": res.opt.get("summed_costs_discounted").tolist(),
        }
    except Exception as e:
        file_path = f"results/{folder_name}/{experiment_id}_{now.strftime('%d-%m-%Y')}_{dataset.name}_res_opt_{i}_fail.json"
        info = {
            "experiment": experiment_name,
            "experiment_id": experiment_id,
            "date": now.strftime("%d/%m/%Y at %H:%M:%S"),
            "target_class": target_class,
            "model_type": model_type,
            "x0_id": i,
            "x0": x0.tolist(),
            "iterations": n_generations,
            "seed": seed,
            "dataset": dataset.name,
            "success": False,
            "exception": str(e),
            "runtime": res.exec_time,
            "n_solutions": 0,
            "X": [],
            "F": [],
            "G": [],
            "pheno": [],
            "tweaked_res": [],
            "summed_costs": [],
            "summed_costs_discounted": [],
        }
        print(f"Run {i}:", e)

    with open(file_path, "w+") as f:
        json.dump(info, f, indent=4)

    # Also save optimal population info per generation
    if save_optimal_population_trace:
        # Create dir if not exists
        folder_name = (
            f"{now.strftime('%d-%m-%Y')} {experiment_id} {dataset.name}/History Results"
        )
        Path(f"results/{folder_name}/").mkdir(parents=True, exist_ok=True)
        try:
            # i denotes the generation here
            all_X = {i: h.opt.get("X").tolist() for i, h in enumerate(res.history)}
            all_F = {i: h.opt.get("F").tolist() for i, h in enumerate(res.history)}
            all_G = {i: h.opt.get("G").tolist() for i, h in enumerate(res.history)}
            all_pheno = {
                i: h.opt.get("pheno").tolist() for i, h in enumerate(res.history)
            }
            all_tweaked = {
                i: h.opt.get("tweaked_instances").tolist()
                for i, h in enumerate(res.history)
            }
            all_summed_costs = {
                i: h.opt.get("summed_costs").tolist() for i, h in enumerate(res.history)
            }
            all_summed_costs_discounted = {
                i: h.opt.get("summed_costs_discounted").tolist()
                for i, h in enumerate(res.history)
            }

            file_path = f"results/{folder_name}/{experiment_id}_{now.strftime('%d-%m-%Y')}_{dataset.name}_history_opt_{i}.json"
            info = {
                "experiment": experiment_name,
                "experiment_id": experiment_id,
                "date": now.strftime("%d/%m/%Y at %H:%M:%S"),
                "target_class": target_class,
                "model_type": model_type,
                "x0_id": i,
                "x0": x0.tolist(),
                "iterations": n_generations,
                "seed": seed,
                "dataset": dataset.name,
                "success": True,
                "runtime": res.exec_time,
                "n_solutions": len(res.X),
                "X": all_X,
                "F": all_F,
                "G": all_G,
                "pheno": all_pheno,
                "tweaked_res": all_tweaked,
                "summed_costs": all_summed_costs,
                "summed_costs_discounted": all_summed_costs_discounted,
            }
        except Exception as e:
            file_path = f"results/{folder_name}/{experiment_id}_{now.strftime('%d-%m-%Y')}_{dataset.name}_history_opt_{i}_fail.json"
            info = {
                "experiment": experiment_name,
                "experiment_id": experiment_id,
                "date": now.strftime("%d/%m/%Y at %H:%M:%S"),
                "target_class": target_class,
                "model_type": model_type,
                "x0_id": i,
                "x0": x0.tolist(),
                "iterations": n_generations,
                "seed": seed,
                "dataset": dataset.name,
                "success": False,
                "exception": str(e),
                "runtime": res.exec_time,
                "n_solutions": 0,
                "X": [],
                "F": [],
                "G": [],
                "pheno": [],
                "tweaked_res": [],
                "summed_costs": [],
                "summed_costs_discounted": [],
            }
            print(f"Run {i}:", e)

        with open(file_path, "w+") as f:
            json.dump(info, f, indent=4)

    return res.exec_time

@ex.automain
def run(_run):
    np.random.seed(_run.config["seed"])

    # Dump config
    identifier = str(uuid.uuid1())[:7]
    ex.add_config({"experiment_id": identifier})
    now = datetime.now()

    dataset, legacy_actions, legacy_features = setup_dataset()
    features = {b: int(a) for a, b in enumerate(dataset.columns)}
    assert np.all([x in dataset.columns for x in features.keys()])
    model = setup_model(dataset)
    blackbox, blackbox_prob = get_blackboxes(model, dataset)
    bounds_and_values = setup_boundaries_and_allowed_values(dataset)

    (
        initial_instances_orig,
        initial_instances_proc,
        instances_random_idx,
    ) = setup_experiment_instances(dataset, model)
    # write the initial instances to a file for later usage
    write_instances_to_file(
        initial_instances_orig,
        _run.config["dataset_name"],
        _run.config["target_class"],
        identifier,
        prefix="original",
    )
    write_instances_to_file(
        initial_instances_proc,
        _run.config["dataset_name"],
        _run.config["target_class"],
        identifier,
        prefix="processed",
    )

    # Write config
    Path(f"results/").mkdir(parents=True, exist_ok=True)
    d_name = _run.config["dataset_name"]
    file_path = f"results/{identifier}_config_{now.strftime('%d-%m-%Y at %H-%M')}_{d_name}_{experiment_name}.json"
    additional_info = {
        "experiment_id": identifier,
        "initial_random_idx": instances_random_idx.tolist(),
    }
    config_dict = {**additional_info, **_run.config.copy()}
    with open(file_path, "w+") as f:
        json.dump(config_dict, f, indent=4)

    if _run.config["parallelized"]:
        # Parellelize
        nprocs = mp.cpu_count()
        print(f"Number of CPU cores: {nprocs}")
        with mp.Pool(processes=nprocs) as pool:
            experiments = [
                (
                    i,
                    x0,
                    legacy_actions,
                    legacy_features,
                    dataset,
                    bounds_and_values,
                    _run.config["target_class"],
                    _run.config["model_type"],
                    _run.config["n_generations"],
                    _run.config["seed"],
                    identifier,
                    _run.config["save_optimal_population_trace"],
                    features,
                    model,
                )
                for i, x0 in enumerate(initial_instances_orig)
            ]
            runtimes = pool.starmap(run_optimization_problem, experiments)
    else:
        # Sequentially
        runtimes = []
        for i, x0 in enumerate(initial_instances_orig):
            print(f"Running experiment for instance {i}")
            time = run_optimization_problem(
                i,
                x0,
                legacy_actions,
                legacy_features,
                dataset,
                bounds_and_values,
                _run.config["target_class"],
                _run.config["model_type"],
                _run.config["n_generations"],
                _run.config["seed"],
                identifier,
                _run.config["save_optimal_population_trace"],
                features,
                model,
            )
            runtimes.append(time)
            print(f"Experiment finished after {time}")
        print(runtimes)
        print(np.mean(runtimes))
        print("Finished")

max_stars_repo_name: ppnaumann/CSCF

max_issues_repo_name: ppnaumann/CSCF

max_forks_repo_name: ppnaumann/CSCF


text:
import matplotlib.pyplot as plt
import numpy as np

m = np.arange(1.0,11.0,1)           #m modes
R = np.arange(30.0,35.0,0.01)       #RoC mirror

# Creating 2-D grid of features

r1 = 0.998                          #R of mirror
r2 = 0.990
L = 9.1                             #length of cav
g = 1-(L/R)                         #g
F = (np.pi*((r1*r2)**0.25))/(1-(r1*r2)**0.5)    #Finesse

print("F",F)
#print("G",g)
[X, Y] = np.meshgrid(m, R)

fig, ax = plt.subplots(1,1)

#Z = 1/(1+((2*F/np.pi)**2)*np.sin(X*np.arccos((1-L/Y)**0.5))**2)**0.5
#Z=  1/(1+((2*F/np.pi)**2)*np.sin(X*np.arccos((1-L/Y)**0.5))**2)**0.5

#relative to fun HOM

Z= 1/((1+((2*F/np.pi)**2)*(np.sin(X*np.arccos((1-L/Y)**0.5)))**2)**0.5)

# plots filled contour plot
cf = ax.contourf(X, Y, Z)

ax.set_title('RoC vs TEM')
ax.set_xlabel('TEM (m)')
ax.set_ylabel('RoC ')
fig.colorbar(cf, ax=ax)
plt.figure(figsize=(8,8))
plt.show()

ext: py

max_stars_repo_path: Fabry Perot Cavity/m+n_analytical.py

max_issues_repo_path: Fabry Perot Cavity/m+n_analytical.py

max_forks_repo_path: Fabry Perot Cavity/m+n_analytical.py

reason: import numpy


text:
# AUTOGENERATED! DO NOT EDIT! File to edit: bcml_model.ipynb (unless otherwise specified).

__all__ = ['bcml_model', 'refresh_model']

# Cell
import numpy as np
import sklearn.metrics as skm
import scipy.interpolate
import scipy.optimize
import pandas as pd
import joblib

# Cell
class bcml_model:
    """
    Represents a machine learning (ML) binary classification (BC) model.
    """

    def __init__(self, model):
        self.model = model

    def fit(self, preds, labels):
        r"""
        Fits `model` to data.

        `preds` should be ? $\times$ $m$ for model with $m$ features. `labels` should
        be ? $\times$ 1 and take values in $\{0,1\}$.
        """
        self.model.fit(preds, labels)

    def predict_proba(self, preds):
        r"""
        Predicts signal probability for each element of a dataset ($? \times m$ `numpy` array).

        Returns `numpy` array of with values in $[0,1]$ giving predicted signal probabilities for
        each data point.
        """
        try:
            return self.model.predict_proba(preds)[:,1]
        except:
            print("self.model doesn't have a predict_proba function")

    def predict(self, preds, threshold=None):
        r"""
        Predicts signal ($1$) or background ($0$) for each element of a dataset ($? \times m$ `numpy` array).

        Returns `numpy` array of length with values in $\{0,1\}$ giving predicted classifications.

        Uses the `predict` method built into `scikit-learn` models.
        """
        if threshold is not None:
            probs = self.model.predict_proba(preds)[:,1]
            return np.where(probs > threshold, np.ones_like(probs), np.zeros_like(probs))
        else:
            try:
                return self.model.predict(preds)
            except:
                print("self.model doesn't have a predict function")

    def predict_hist(self, preds, labels, num_bins=100, sepbg=False, sig_norm=1, bg_norm=1, dataframe=False):
        r"""
        Constructs a histogram of predicted signal probabilities for signal and background constituents of
        a dataset ($? \times M$ `numpy` array).

        If `sepbg` is `False` (the default), labels are assumed to take values in $\{0,1\}$. Backgrounds are treated in combination:
        a list of $3$ $?_i \times m$ `numpy` arrays are returned, containing bin edges (partitioning $[0,1]$),
        signal bin contents, and background bin contents.

        If `sepbg` is `True`, labels are assumed to take values in $\{-n,\dots,-1,1\}$ (if there are $n$ backgrounds) while
        `bg_norm` should be a list of length $n$. Backgrounds are then differentiated: a list of $2 + n$ `numpy` arrays of shape
        $?_i \times m$ are returned, containing bin edges (partitioning $[0,1]$), signal bin contents, and
        $n$ background bin contents.
        """
        predictions = self.predict_proba(preds)
        labels = np.array(labels)
        sig_bins, bin_edges = np.histogram(predictions[labels==1], bins=num_bins, density=True)
        sig_bins *= sig_norm
        if sepbg:
            bg_norms = bg_norm
            bg_binss = [
                bg_norm * np.histogram(predictions[labels==-(i+1)], bins=num_bins, density=True)[0]
                for i, bg_norm in enumerate(bg_norms)]
            if dataframe:
                return pd.DataFrame(data=[bin_edges, sig_bins] + bg_binss, columns=['Bin Edges', 'Signal'] + ['Background {}'.format(i) for i in range(1, self.num_bgs+1)])
            else:
                return [bin_edges, sig_bins] + bg_binss
        else:
            bg_bins = np.histogram(predictions[labels!=1], bins=num_bins, density=True)[0]
            if dataframe:
                return pd.DataFrame(data=[bin_edges, sig_bins, bg_bins], columns=['Bin Edges', 'Signal', 'Background'])
            else:
                return [bin_edges, sig_bins, bg_bins]

    def feature_importance(self):
        """
        Returns the importance of the $M$ features used to train `self.model`.
        """
        try:
            return self.model.feature_importances_
        except:
            try:
                return self.model[-1].feature_importances_
            except:
                print("It looks like self.model doesn't have an attribute 'feature_importances_'")

    def sorted_feature_importance(self, features):
        r"""
        Returns list of features sorted by importance.

        Given arguments `features` and `importances`, lists of length $M$, returns list of size $M \times 2$ where
        the first column gives features and the second their associated importances, sorted by importance.
        """
        importances = self.feature_importance()
        ranked_indices = np.argsort(-np.abs(importances))
        return [[features[i], importances[i]] for i in ranked_indices]

    def accuracy(self, preds, labels, threshold=None):
        r"""
        Computes model accuracy on a dataset ($? x m$ predictors, length $?$ labels).

        Returns value in $[0,1]$ giving model accuracy on the provided predictors and labels.
        """
        predictions = self.predict(preds=preds, threshold=threshold)
        return len(preds) - np.sum(np.abs(predictions - labels)) / len(preds)

    def conf_matrix(self, labels, predictions=None, preds=None):
        r"""
        Computes the confusion matrix of the trained model on a dataset ($? x M$ predictors, length $?$ labels).

        Returns $2 \times 2$ confusion matrix using `sklearn.metrics.confusion_matrix`.

        If predictors `preds` aren't provided, `self.test_preds` is used.
        If `labels` aren't provided, `self.test_labels` is used.
        """
        if predictions is not None:
            return skm.confusion_matrix(labels, predictions, labels=[0,1])
        elif preds is not None:
            return skm.confusion_matrix(labels, self.predict(preds), labels=[0,1])
        else:
            raise ValueError('Either predictions or preds must be passed.')

    def tpr_cm(self, conf_matrix):
        """
        Computes the true positive rate (tpr; correctly identified signal/total signal)
        of a trained model given a confusion matrix.

        Returns value in $[0,1]$.
        """
        return conf_matrix[1,1]/np.sum(conf_matrix[1])

    def fpr_cm(self, conf_matrix):
        """
        Computes the false positive rate (fpr; misidentified background/total background)
        of a trained model given a confusion matrix.

        Returns value in $[0,1]$.
        """
        return conf_matrix[0,1]/np.sum(conf_matrix[0])

    def tpr(self, labels, predictions=None, preds=None):
        r"""
        Computes the true positive rate (tpr; correctly identified signal/total signal)
        of a trained model given predictions and labels (both `numpy` array of length $?$ with values in $\{0,1\}$)

        Returns value in $[0,1]$.
        """
        return self.tpr_cm(self.conf_matrix(labels, predictions=predictions, preds=preds))

    def fpr(self, labels, predictions=None, preds=None):
        r"""
        Computes the false positive rate (fpr; misidentified background/total background)
        of a trained model given predictions and labels (both `numpy` array of length $?$ with values in $\{0,1\}$)

        Returns value in $[0,1]$.
        """
        return self.fpr_cm(self.conf_matrix(labels, predictions=predictions, preds=preds))

    def significance(self, signal, background, tpr, fpr, sepbg=False):
        r"""
        Computes signal significance of a trained model given signal and background yield.

        Returns a positive real number computed by
        $$\frac{S \cdot TPR}{\sqrt{S \cdot TPR + B \cdot FPR}}$$
        which corresponds to signal significance after selecting only datapoints the model identifies as signal.

        If `sepbg` is `False`, `background` should be a single real number and is multiplied by `fpr`. If `sepbg` is `True`,
        `background` should be a list of length `self.num_bgs` where the $i$th element contains background yield of the $i$th
        background type. `fpr`, if passed, is then also a list of length `self.num_bgs` giving false positive rates for each
        of the background types.
        """
        if sepbg:
            backgrounds = background
            fprs = fpr
            fprXbackground = np.sum(np.multiply(fprs, backgrounds), axis=-1)
            return (signal * tpr) / np.sqrt(signal * tpr + fprXbackground + 1e-10)
        else:
            return (signal * tpr) / np.sqrt(signal * tpr + background * fpr + 1e-10)

    def newvar2thresh(self, newvar):
        r"""
        Helper method for `bcml.max_allowable_threshold()`, `bcml.get_tprs_fprs()`, and `bcml.best_threshold()`,
        performing change of variables from `newvar` to `threshold`

        In particular, threshold $= 1 - 10^{\text{newvar}}$
        """
        return 1 - np.power(10, newvar)

    def thresh2newvar(self, thresh):
        r"""
        Helper method for `bcml.max_allowable_threshold()`, `bcml.get_tprs_fprs()`, and `bcml.best_threshold()`,
        performing change of variables from `threshold` to `newvar`

        In particular, newvar $= \log_{10}(1 - \text{threhold})$
        """
        return np.log10(1 - thresh)

    def max_allowable_threshold(self, preds, labels, sigYield):
        """
        Returns the highest threshold such that only labelling elements of `self.test_pred` with predicted
        probabilities higher than that threshold as signal still yields 25 signal.

        To achieve a discovery potential of $5\sigma$, even in the best case scenario ($TPR = 1, FPR = 0$) we still
        require $5^2 = 25$ signal events, hence we cannot chose a threshold so high that we do not keep at least
        25 signal events.
        """

        sig_indx = np.where(np.array(labels)==1)[0]
        preds = np.array(preds)[sig_indx]
        probs = self.predict_proba(preds)
        num_predicts = np.array(
            [np.sum(np.where(probs > self.newvar2thresh(newvar),np.ones_like(probs), np.zeros_like(probs)))
             for newvar in newvars])
        num_sig_yields = (num_predicts / len(preds)) * sigYield
        return [newvars, num_sig_yields]
#         f = scipy.interpolate.interp1d(num_sig_yield, newvars, kind='cubic')
#         return self.newvar2thresh(f(25))

    def get_tprs_fprs(self, preds, labels, sepbg=False):
        """
        Produces (true positive rate, false positive rate) pairs for various thresholds
        for the trained model on data sets.

        If `sepbg` is `True`, labels should take values in $\{-n,\dots,-1,1\}$. Background is combined and a list of length $4$
        is returned containing a list of $L$ sampled newvars (a convenient change of variable to approach arbitrarily close to 1:
        related to thresholds by  `bcml_model.newvar2thresh()`), an $L$-list of tprs associated to those thresholds, an $L$-list of fprs
        related to those thresholds, and an $L$-list of length $?$ `numpy` arrays giving the predicted signal probabilities
        for the given data set.

        If `sepbg` is `Frue`, labels should take values in $\{0,1\}$. Background is split and a list of length $4$ `self.num_bgs`
        is returned containing a  list of $L$ sampled newvars, an $L$-list of tprs associated to those thresholds, an $L$-list of lists of length
        $n$ (number of backgrounds) containing fprs for each background type for each threshold, and an $L$-list of length $?$
        """
        # setting up variables
        min_newvar, max_newvar = [-10, 0]
        newvars = np.concatenate((np.linspace(min_newvar, -2, 10, endpoint=False), np.linspace(-2, max_newvar, 15, endpoint=False)))

        # computing tprs, fprs
        if sepbg:
            num_bgs = len(np.unique(labels)) - 1
            labelTypes = [1] + [-(i+1) for i in range(num_bgs)]
            labelsIndices = [np.where(np.array(labels)==i)[0] for i in labelTypes]
#             print(np.array(labelsIndices).shape)
            predss = [preds[indices] for indices in labelsIndices]
#             print(predss[0][:10])
            probss = [self.predict_proba(preds) for preds in predss]
#             print(probsss[0][:10])
            predictionsss = np.array(
                [[np.where(probs > self.newvar2thresh(newvar),
                          np.ones_like(probs), np.zeros_like(probs)) for probs in probss] for newvar in newvars])
            sig_conf_matrices = [
                self.conf_matrix(labels=np.ones_like(predictionss[0]), predictions=predictionss[0]) for predictionss in predictionsss]
            bg_conf_matricess = [
                [self.conf_matrix(labels=np.zeros_like(predictions), predictions=predictions) for i, predictions in enumerate(predictionss[1:])] for predictionss in predictionsss]
            tprs = np.array([self.tpr_cm(conf_matrix) for conf_matrix in sig_conf_matrices])
            fprss = np.array([[self.fpr_cm(conf_matrix) for conf_matrix in conf_matrices] for conf_matrices in bg_conf_matricess])
            sums = tprs + np.sum(fprss, axis=1)
            cutoff = len(sums) - np.argmax(np.flip(sums)==0) + 1 if 0 in sums else 0
            return [newvars[cutoff:], tprs[cutoff:], fprss[cutoff:], probss]
        else:
            probs = self.predict_proba(preds)
            predictionss = np.array(
                [np.where(probs > self.newvar2thresh(newvar),
                          np.ones_like(probs), np.zeros_like(probs)) for newvar in newvars])
            conf_matrices = [self.conf_matrix(labels=labels, predictions=predictions) for predictions in predictionss]
            tprs = np.array([self.tpr_cm(conf_matrix) for conf_matrix in conf_matrices])
            fprs = np.array([self.fpr_cm(conf_matrix) for conf_matrix in conf_matrices])
            sums = tprs + fprs
            cutoff = len(sums) - np.argmax(np.flip(sums)==0) + 1 if 0 in sums else 0
            return [newvars[cutoff:], tprs[cutoff:], fprs[cutoff:], probs]

    def best_threshold(self, signal, background, preds, labels, sepbg=False):
        """
        Optimizes the threshold on a given data set ($? x M$ predictors, length $?$ labels).
        """
        newvars, tprs, fprs, probs = self.get_tprs_fprs(preds, labels, sepbg)
        significances = -self.significance(signal, background, tprs, fprs, sepbg=sepbg)

        # interpolating significance as a function of threshold, then maximizing
        max_sig = np.amin(significances)
        significances_list = list(significances)
        i = significances_list.index(max_sig)
        min_i, max_i = [max(0,i-4),min(len(significances),i+5)]
        f = scipy.interpolate.interp1d(newvars[min_i:max_i], significances[min_i:max_i], kind='cubic')
        res = scipy.optimize.minimize(f, [0.5 * (newvars[min_i] + newvars[max_i-1])], bounds=[(newvars[min_i] + 1e-1, newvars[max_i-1] - 1e-1)])

        # computing significance, tpr, fpr for optimized threshold
        best_threshold = self.newvar2thresh(res.x[0])
        if sepbg:
            probss = probs
            fprss = fprs
            best_predictss = [np.where(probs > best_threshold, np.ones_like(probs), np.zeros_like(probs)) for probs in probss]
            sig_conf_matrix = self.conf_matrix(labels=np.ones_like(best_predictss[0]), predictions=best_predictss[0])
            bg_conf_matrices = [self.conf_matrix(labels=np.zeros_like(best_predicts), predictions=best_predicts) for i, best_predicts in enumerate(best_predictss[1:])]
            tpr = self.tpr_cm(sig_conf_matrix)
            fprs = [self.fpr_cm(conf_matrix) for conf_matrix in bg_conf_matrices]
            best_sig = self.significance(signal, background, tpr, fprs, sepbg=sepbg)
            return [best_threshold, best_sig, tpr, fprs, tprs, fprss]
        else:
            best_predicts = np.where(probs > best_threshold, np.ones_like(probs), np.zeros_like(probs))
            conf_matrix = self.conf_matrix(labels=labels, predictions=best_predicts)
            tpr = self.tpr_cm(conf_matrix)
            fpr = self.fpr_cm(conf_matrix)
            best_sig = self.significance(signal, background, tpr, fpr, sepbg=sepbg)
            return [best_threshold, best_sig, tpr, fpr, tprs, fprs]

    def req_sig_cs(self, lumi, bg_cs, tpr, fpr, sig=5, sepbg=False):
        """
        Given a luminosity (in fb$^{-1}$), a background cross section (in pb), a true positive rate, a false positive rate,
        and a signal significance, computes the signal cross section required for the signal significance to be achieved.

        If `sepbg` is False, background is combined and a single FPR is used; if `sepbg` is True, it is assumed that
        `bg_cs`, `fpr` are each lists of length $n$ (number of backgrounds) and their vector dot product is used for background yield.
        """
        conv = 10**15 / 10**12
        if sepbg:
            bg = np.sum(np.multiply(bg_cs, fpr))
            coef = [-tpr**2 * lumi * conv**2, sig**2 * tpr * conv, sig**2 * bg * conv]
        else:
            coef = [-tpr**2 * lumi * conv**2, sig**2 * tpr * conv, sig**2 * fpr * bg_cs * conv]
        return np.amax(np.roots(coef))

    def save_model(self, filename):
        """
        Saves the model to `filename.joblib`
        """
        joblib.dump(self.model, filename + '.joblib')

# Cell
def refresh_model(model):
        """
        If this class gets updated, run this function on your already trained model to have it reflect the updated
        class without retraining being necessary.
        """
        return bcml_model(model.model)

reason: import numpy,import scipy


text:
import matplotlib.pyplot as plt
import argparse
import pandas as pd
import numpy as np
import os
from tqdm import tqdm
from pint import UnitRegistry
ur = UnitRegistry()

parser = argparse.ArgumentParser(description='Compute the smoothing tree of a halo.')
parser.add_argument('--in', dest='infile', type=str, help='Path to the output file of compute_extrema (in hdf format)', required=True)
parser.add_argument('--out', '-o', type=str, help='Prefix of the outputs')
parser.add_argument('--infofile', type=str, help='Path to the information file (the one containing units of RAMSES, …)')

args = parser.parse_args()

# read the info file
infos = dict()
infos['headers'] = pd.read_csv(args.infofile, sep=' *= *', nrows=19, names=['key', 'value'], index_col='key').T
infos['domain'] = pd.read_csv(args.infofile, delim_whitespace=True, skiprows=20)

# read the center
from scipy.io import FortranFile
ff = FortranFile('data/halo_536-centers.bin')
ff.read_ints() # don't care
outputs = ff.read_ints()
centers = ff.read_reals().reshape(len(outputs), 3)
mins    = ff.read_reals().reshape(len(outputs), 3)
span    = ff.read_reals().reshape(len(outputs), 3)
maxs    = ff.read_reals().reshape(len(outputs), 3)

# create the output dir if required
# if not os.path.isdir(args.out):
#     os.mkdir(args.out)

HDF = pd.HDFStore(args.infile)
# read the data
df    = HDF['extremas']
dens  = HDF['dens']
edges = HDF['edges']

def to(v, unit):
    return infos['headers']['unit_l'].value * np.array(v) * (1*ur.cm).to(unit)

nbin = len(edges) - 1

# estimate positions
def realPos(S, edges):
    ''' Convert the sery into real positions in cm'''
    m, M = edges.min(), edges.max()
    return to(((S-1) / (nbin - 1) * M + m), 'Mpc')

# compute 'curvature' (product of eigen values of the hessian)
df['curvature'] = np.abs(df.evx * df.evy * df.evz)
df['rx'] = realPos(df.x, edges.x)
df['ry'] = realPos(df.y, edges.y)
df['rz'] = realPos(df.z, edges.z)

# some values have x,y,z above maximum, remove them
posMask = ((df.x >= 0) & (df.x < nbin) &
           (df.y >= 0) & (df.y < nbin) &
           (df.z >= 0) & (df.z < nbin))

# plots configuration
import matplotlib
from mpl_toolkits.mplot3d import Axes3D
matplotlib.style.use(['ggplot'])
import matplotlib.pyplot as plt

fig = plt.figure()
############## plot sigma vs x
# plt.scatter('x', 'sigma', data=df[df.type == 4], c='r', marker='o', label='Maximum')
# plt.scatter('x', 'sigma', data=df[df.type == 3], c='b', marker='x', alpha=0.5, label='Point selle, 2 maximums')
# plt.xlim(0, len(edges))

# plt.legend(loc='best')
# plt.xlabel('x')
# plt.ylabel('$\sigma$')

############## plot sigma vs dist to center
# plt.scatter('dist', 'sigma', data=df[df.type == 4], c='r', marker='o', label='Maximum')
# plt.scatter('dist', 'sigma', data=df[df.type == 3], c='b', marker='.', alpha=0.5, label='Point selle, 2 maximums')

############## plot sigma vs x,y
ax = fig.add_subplot(111, projection='3d')

t = 4
# map evz onto 20-120
size_evz = 50*(df[df.type==t].evz.max()/df[df.type==t].evz) + 10
m = posMask & (df.type == t) & (df.sigma >= 5) & (df.i < nbin) & (df.j < nbin) & (df.k < nbin)
x, y = df[m].rx, df[m].ry
ax.scatter3D(x, y, df[m].sigma,
             s=size_evz,
             c='r', marker='o', label='Maximum')

t = 3
m = posMask & (df.type == t) & (df.sigma >= 5) & (df.i < nbin) & (df.j < nbin) & (df.k < nbin)
x, y = df[m].rx, df[m].ry
ax.scatter3D(x, y, df[m].sigma,
             c='b', marker='.', alpha=0.5, label='Point selle, 2 maximums')

t = 2
m = posMask & (df.type == t) & (df.sigma >= 5) & (df.i < nbin) & (df.j < nbin) & (df.k < nbin)
x, y = df[m].rx, df[m].ry
ax.scatter3D(x, y, df[m].sigma,
             c='g', marker='.', alpha=0.1, label='Point selle, 1 maximums')

# t = 1
# m = posMask & (df.type == t) & (df.sigma >= 5) & (df.i < nbin) & (df.j < nbin) & (df.k < nbin)
# x, y = df[m].rx, df[m].ry
# ax.scatter3D(x, y, df[m].sigma,
#              c='grey', marker='+', alpha=1, label='Minimum')

# ax.set_xlim(to(edges.x.min(), 'Mpc').m, to(edges.x.max(), 'Mpc').m)
# ax.set_ylim(to(edges.y.min(), 'Mpc').m, to(edges.y.max(), 'Mpc').m)
ax.set_xlabel('x [Mpc/ h]')
ax.set_ylabel('y [Mpc/ h]')
ax.set_zlabel('$\sigma$')

cMpc = to(centers, 'Mpc').m
ax.plot([cMpc[-1, 0], cMpc[0, 0]], [cMpc[-1, 1], cMpc[0, 1]], [df.sigma.min(), df.sigma.max()],
        label='Halo', c='black')

ax.legend()
plt.tight_layout()
# for azim in tqdm(np.linspace(0, 360)):
#     ax.view_init(elev=0, azim=azim)
#     ax.legend()

#     # plt.savefig('/tmp/plots/%s.pdf' % azim)
#     plt.savefig('/tmp/plots/%s.png' % azim, dpi=180)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

sc = ax.scatter3D(df[m].rx, df[m].ry, df[m].rz, c=df[m].sigma, s=5)
cb = fig.colorbar(sc)

ax.set_xlabel('x [Mpc/h]')
ax.set_ylabel('y [Mpc/h]')
ax.set_zlabel('z [Mpc/h]')

ax.set_xlim()
cb.set_label('$\sigma$')
fig.tight_layout()

plt.show()

reason: import numpy,from scipy


text:
"""
This script plots the average and standard deviation of the norm of the
difference between the true and estimated parameter values.

There are separate plots for the difference in delta, the offset parameter,
and sigma, the standard deviation of the generating process.

There are two sets of plots: one comparing the results across several methods
and the other comparing HVAE with and without tempering.
"""
import os
import pickle
import numpy as np
import matplotlib.pyplot as plt

from conf import params
from plot_utils import parse_pickle_save, plot_avg_diff

# Plot specifications
colours = ['b', 'g', 'r', 'c', 'k']
markers = ['d', '+', 's', '*', 'x']
fig_size = (6, 4.5)

# Method names for the tests and plot legends
#   'all' corresponds to the test of all methods
#   'temp' corresponding to the tempering test
methods_all = ['HVAE_1', 'HVAE_2', 'NF_1', 'NF_2', 'VB']
methods_temp = ['HVAE_1', 'HVAE_2', 'HVAE_notemp_1', 'HVAE_notemp_2']
plotnames_all = [r'HVAE, $K={0}$'.format(params['HVAE_K_1']),
                 r'HVAE, $K={0}$'.format(params['HVAE_K_2']),
                 r'NF, $K={0}$'.format(params['NF_K_1']),
                 r'NF, $K={0}$'.format(params['NF_K_2']),
                 'VB']
plotnames_temp = [r'HVAE, $K={0}$'.format(params['HVAE_K_1']),
                  r'HVAE, $K={0}$'.format(params['HVAE_K_2']),
                  r'HVAE, $K={0}$, $\beta_0=1$'.format(params['HVAE_K_1']),
                  r'HVAE, $K={0}$, $\beta_0=1$'.format(params['HVAE_K_2'])]

# Extract the parameter difference information from the pickle files
diffs_all = parse_pickle_save(methods_all)
diffs_temp = parse_pickle_save(methods_temp)

# Make the plots testing all methods with error bars
if not os.path.exists('plots'):
    os.makedirs('plots')

plot_avg_diff(colours, markers, 1, fig_size, plotnames_all, diffs_all,
    r'Error in estimate of $\theta$')
plt.savefig('plots/all_methods.png')

# Make the plot comparing tempering vs. no tempering
# \beta_0 = 1 corresponds to no tempering
plot_avg_diff(colours, markers, 1, fig_size, plotnames_temp, diffs_temp,
    r'Error in estimate of $\theta$')
plt.savefig('plots/temp_vs_no_temp.png')

max_stars_repo_name: anthonycaterini/hvae-nips

max_issues_repo_name: anthonycaterini/hvae-nips

max_forks_repo_name: anthonycaterini/hvae-nips


text:
from flask import Flask, jsonify, request, send_from_directory
import omdb;
from pymongo import MongoClient
import math;
import datetime
import hashlib
from flask_uploads import (UploadSet, configure_uploads, IMAGES,
                              UploadNotAllowed)
from sklearn.cluster import KMeans
import numpy as np
import sys
import pyfpgrowth
# App

UPLOADED_PHOTOS_DEST = '/tmp/photolog'
ALLOWED_EXTENSIONS = set(['txt', 'pdf', 'png', 'jpg', 'jpeg', 'gif'])

app = Flask(__name__, static_url_path='')

# Global Vars
omdb.api._client.params_map['apikey'] = 'apikey'
omdb.set_default('apikey', '26029d08')
photos = UploadSet('photos', IMAGES,default_dest=lambda app:app.instance_path)
configure_uploads(app, photos)
count = 0
user_count = 0
batch_size = 12
client = MongoClient('mongodb://localhost:27017/')
mydb = client['movie_database']
user_rating_matrix = []
movie_index_dict = {}
user_index_dict = {}
transactions = []
genres = []
lastCall = 0;
recommendation = {
    "predicted_numpy_matrix" : np.array(user_rating_matrix)
}

def init_user_rating_matrix():
    user_rating_cursor = mydb.user_ratings.find()

    for user_rating_item in user_rating_cursor:
        insert_into_matrix(user_rating_item.get("user_id"), user_rating_item.get("movie_id"), user_rating_item.get("rating"))
    print "[info] Loading initial user matrix complete"

def convert_to_numpy(list_matrix):
    max_len = len(max(list_matrix,key=len))
    for list_item in list_matrix:
        current_len = len(list_item)
        for i in xrange(max_len - current_len):
            list_item.append(0)
    return np.array(list_matrix)

def key_from_value(value, dict):
    for key in dict.keys():
        if(dict.get(key) == value):
            return key

def matrix_factorization(R, P, Q, K, steps=5000, alpha=0.0002, beta=0.02):
    Q = Q.T
    for step in xrange(steps):
        for i in xrange(len(R)):
            for j in xrange(len(R[i])):
                if R[i][j] > 0:
                    eij = R[i][j] - np.dot(P[i,:],Q[:,j])
                    for k in xrange(K):
                        P[i][k] = P[i][k] + alpha * (2 * eij * Q[k][j] - beta * P[i][k])
                        Q[k][j] = Q[k][j] + alpha * (2 * eij * P[i][k] - beta * Q[k][j])
        eR = np.dot(P,Q)
        e = 0
        for i in xrange(len(R)):
            for j in xrange(len(R[i])):
                if R[i][j] > 0:
                    e = e + pow(R[i][j] - np.dot(P[i,:],Q[:,j]), 2)
                    for k in xrange(K):
                        e = e + (beta/2) * (pow(P[i][k],2) + pow(Q[k][j],2))
        if e < 0.001:
            break
    return P, Q.T
def predict_numpy_matrix():
    print "[info] Predictng numpy matrix"

    R = convert_to_numpy(user_rating_matrix)
    N = len(R)
    M = len(R[0])
    K = 5
    P = np.random.rand(N,K)
    Q = np.random.rand(M,K)
    nP, nQ = matrix_factorization(R, P, Q, K)
    nR = np.dot(nP, nQ.T)
    nR = np.around(nR, decimals = 3)
    recommendation["predicted_numpy_matrix"] = nR
    print "[Done] Predictng numpy matrix"

def insert_into_matrix(user_id, movie_id, rating):
    matrix_len = len(user_rating_matrix)

    user_index = get_item_index(user_id, user_index_dict)
    movie_index = get_item_index(movie_id, movie_index_dict)
    if user_index >= matrix_len:
        itr_count = user_index - matrix_len + 1
        for i in xrange(itr_count):
            user_rating_matrix.append([])
    tempArray = user_rating_matrix[user_index]
    tempArray_len = len(tempArray)
    if movie_index < tempArray_len:
        tempArray[movie_index] = rating
    else:
        itr_count = movie_index - tempArray_len
        for i in xrange(itr_count):
            tempArray.append(0)
        tempArray.append(rating)
    user_rating_matrix[user_index] = tempArray

def get_item_index(item_id, dict):
    if item_id in dict.keys():
        return dict[item_id]
    else:
        current_length = len(dict)
        dict[item_id] = current_length
        return current_length
# Routes
@app.route('/_uploads/photos/<path:path>')
def send_js(path):
    return send_from_directory('instance', path)

@app.route('/')
def index():
    return "Hello, World!"

@app.route('/login', methods=['POST'])
def login():
    print "I do came here.."
    content = request.get_json(force=True)
    print request
    print content

    username = content.get("username", None)
    password = content.get("password", None)
    password = computeMD5hash(password)
    print username

    user = mydb.users.find_one({"_id" : username, "password" : password});
    if user :
        valid = True;
        return jsonify({"valid" : valid, "mod" : user.get("moderator")})
    else:
        return jsonify({"valid" : False})

@app.route('/rateMovie', methods=["POST"])
def rateMovie():
    content = request.get_json(force=True)
    print content
    username = content.get("username", None)
    movie_id = content.get("movie_id", None)
    valid_user = mydb.users.find_one({"_id" : username});
    valid_movie = mydb.movies.find_one({"_id": movie_id})

    if not valid_user or not valid_movie:
        return jsonify({"message" : "Invalid User or Movie"})

    rating = content.get("rating", 0)
    if username and movie_id and rating :
        rating_item = mydb.user_ratings.find_one({"user_id" : username, "movie_id" : movie_id});
        if rating_item :
            mydb.user_ratings.update_one({"_id" : rating_item.get("_id")}, {"$set": {"rating" : rating, "time_stamp" : datetime.datetime.now().isoformat()}})
            createActivity(valid_user.get("profile_url"),username, "rated", valid_movie.get("title"), rating)

        else:
            mydb.user_ratings.insert_one({"user_id" : username, "movie_id" : movie_id,"rating": rating, "time_stamp" : datetime.datetime.now().isoformat()})
            createActivity(valid_user.get("profile_url"),username, "rated", valid_movie.get("title"), rating)
        insert_into_matrix(username, movie_id, rating)
        convert_to_numpy(user_rating_matrix)
        return jsonify({"message" : "success"})
    else:
        return jsonify({"message" : "failure"})

@app.route('/register', methods=["POST"])
def register():
    content = request.get_json(force=True)
    print content
    username = content.get("username", None)
    password = content.get("password", None)
    moderator = content.get("moderator", False)
    print moderator
    url =  content.get("profile_url", None)
    if request.method == 'POST' and 'photo' in request.files:
        filename = photos.save(request.files['photo'])
        url = photos.url(filename)

    password = computeMD5hash(password)
    if moderator :
        moderator = True
    user = mydb.users.find_one({"_id" : username});
    if user:
        return jsonify({"message" : "User Already Exists!"})
    else :
        valid = True
        mydb.users.insert_one({
            "_id" : username,
            "password" : password,
            "moderator" : moderator,
            "time_stamp" : datetime.datetime.now().isoformat(),
            "profile_url" : url
        })
    return jsonify({"message" : "success"})

@app.route('/save_movie/<string:id>')
def saveMovie(id):
    result = mydb.movies.find({"_id" : id})
    if count_iterable(result) < 1:
        movieInfo = omdb.imdbid(id)
        movieInfo["_id"] = movieInfo["imdb_id"]
        movieInfo["time_stamp"] = datetime.datetime.now().isoformat();
        mydb.movies.insert_one(movieInfo)
        print movieInfo.title +" Inserted!"
        calculate_movies_count();
    return jsonify({"message" : "success"})

@app.route('/recommend/<string:id>')
def recommend(id):
    global lastCall
    lastCall += 1
    i = get_item_index(id, user_index_dict)
    if i >= len(user_rating_matrix):
        return jsonify({"movies":[], "message" : "Insufficient movies"})
    print str(i) + " <<>> " + str(len(recommendation["predicted_numpy_matrix"]) )

    print str(lastCall) + " <<LastCall>> "
    if(lastCall % 10 == 0 or i >= len(recommendation["predicted_numpy_matrix"])):
        print "Re Calculating"
        predict_numpy_matrix()
    else:
        print "Re Using"
    nR = recommendation["predicted_numpy_matrix"]

    a = nR[i].tolist()
    max_i =  np.argsort(a)[::-1]
    result = []
    max_i = max_i.tolist()
    total_movies = 0
    for index in max_i:
        if total_movies > 5:
            break
        if user_rating_matrix[i][index] == 0:

            if a[index] > 2.5:
                total_movies = total_movies + 1
                movie_id = key_from_value(index, movie_index_dict)
                movie = mydb.movies.find_one({"_id" : movie_id})
                if a[index] > 5:
                    movie["predicted_rating"] = 5
                else:
                    movie["predicted_rating"] = a[index]
                result.append(movie)

    return jsonify({"movies" : result})
@app.route('/delete_movie/<string:id>')
def deleteMovie(id):
    result = mydb.movies.find({"_id" : id})
    if count_iterable(result) > 0:
        movieInfo = omdb.imdbid(id)
        movieInfo["_id"] = movieInfo["imdb_id"]
        mydb.movies.delete_many({"_id" : id})
        print movieInfo.title +" Deleted!"
        calculate_movies_count();
    return jsonify({"message" : "success"})

@app.route('/movies/<int:page_no>')
def get_movies(page_no ):
    if get_movies_count() < 1:
        calculate_movies_count();

    if request.args :
        username = request.args.get("username", None)
    print username
    ratedMovies = []
    if username:
        movie_iter = mydb.user_ratings.find({"user_id" : username})
        for movie in movie_iter:
            print movie
            ratedMovies.append(movie)

    movies = mydb.movies.find().sort("time_stamp", -1).skip((page_no - 1) * batch_size).limit(batch_size)
    result = []
    for movie in movies:
        rating_item = isInArray(ratedMovies, movie, "movie_id", "_id")
        if rating_item:
            movie["user_rating"]=rating_item.get("rating")
        result.append(movie)
    return jsonify({"count": get_movies_count(), "movies": result})

@app.route('/users/<int:page_no>')
def get_users(page_no ):
    if get_users_count() < 1:
        calculate_users_count();

    users = mydb.users.find().sort("time_stamp", -1).skip((page_no - 1) * batch_size).limit(batch_size)
    result = []
    for user in users:
        result.append(user)
    return jsonify({"count": get_users_count(), "users": result})

@app.route('/fpRecommender/<string:id>')
def getFPRecommendations(id):
    enrichedRecList = []
    minsup = 3
    if request.args and request.args.get("minsup"):
        minsup = request.args.get("minsup")

    recommendations = fpRecommender(id, int(minsup))
    movies_sample = []
    for recommendation in recommendations:
        enrichedRec = []
        for movie in recommendation.get("reason", []):
            enrichedRec.append(mydb.movies.find_one({"_id" : movie})["title"])
        if(len(enrichedRec) > 1):
            enrichedRecList.append({
                "reason" : enrichedRec,
                "movies" : recommendation.get("movies", [])
            })
        movies_sample = concatWithoutDuplicates(movies_sample, recommendation.get("movies", []))
    print movies_sample
    enrichedMoviesSample = []
    for movie in movies_sample:
        enrichedMoviesSample.append(mydb.movies.find_one({"_id" : movie}))
    labels =  formClusters(n_cl = 2, movies = enrichedMoviesSample).labels_
    return jsonify({"recommendation" : enrichedRecList, "movie_sample" : enrichedMoviesSample, "labels" : labels.tolist()})

@app.route('/similar/<string:id>')
def getSimilar(id):
    movies_c = mydb.movies.find();
    movie_sample = []
    index = 0;
    i = 0;
    for movie in movies_c:
        if movie["_id"] == id:
            index = i
            print "selected index " + str(index)
        i += 1
        movie_sample.append(movie)

    print i;
    n_c = i/5;
    print n_c
    labels = formClusters(n_cl = n_c, movies = movie_sample).labels_
    print labels
    class_m = labels[index]
    result = []
    i = 0
    for l in labels:
        if l == class_m and index != i:
            result.append(movie_sample[i])
        i += 1

    return jsonify(result)

def concatWithoutDuplicates(list1, list2):
    return list1 + list(set(list2) - set(list1))

def updateTransactions(minsup):
    transactions = []
    users = mydb.users.find();
    for user in users:
        movie_id_list = [];
        ratings = mydb.user_ratings.find({"user_id" : user.get("_id"), "rating" : {"$gt" : 2.9}})
        for rating in ratings:
            movie_id_list.append(rating.get("movie_id"))
        if(len(movie_id_list) > 0):
            transactions.append(movie_id_list)
    size = 0;
    itemsets = pyfpgrowth.find_frequent_patterns(transactions, minsup)
    patterns = []
    for itemset in itemsets:
        if(len(itemset) > 1):
            patterns.append(itemset)

    return patterns

def intersect(a, b):
    return list(set(a) & set(b))

def difference(a,b):
    return list(set(a).symmetric_difference(set(b)))

def fpRecommender(user_id, minsup):
    ratings = mydb.user_ratings.find({"user_id" : user_id, "rating" : {"$gt" : 2.9}})
    user_movies = []
    recommended_movies = []
    for rating in ratings:
        user_movies.append(rating.get("movie_id"))
    patterns = updateTransactions(minsup)
    for pattern in patterns:
        intersection = intersect(user_movies, pattern)
        if(len(intersection) > 0 and len(intersection) < len(pattern)):

            print float(len(intersection))/len(pattern)
            if(float(len(intersection))/len(pattern) > 0.4):
                diff = difference(pattern, intersection)
                recommendation = {
                    "reason" : intersection,
                    "movies" : diff
                }
                recommended_movies.append(recommendation)
    return recommended_movies

@app.route('/movie/<string:id>')
def get_movie(id):
    movie = mydb.movies.find_one({"_id" : id})
    return jsonify(movie)
@app.route('/search_movies')
def search_movies():
    search_string = ""
    response = []
    print "hello"
    if request.args and request.args.get("s"):
        search_string = request.args.get("s")
        auto_store = request.args.get("auto_store")

        result = omdb.search_movie(search_string)

        for movie in result:
            try:
                print movie.title
                movie["_id"] = movie["imdb_id"]
                if movie["poster"].endswith(".jpg"):
                    response.append(movie)
                    alreadyExists = mydb.movies.find({"_id" : movie["imdb_id"]})
                    alreadyExists = count_iterable(alreadyExists) > 0
                    if(not alreadyExists and auto_store == "1"):
                        print "going to save a movie " + movie["_id"]
                        movieInfo = omdb.imdbid(movie["_id"])
                        movieInfo["_id"] = movieInfo["imdb_id"]
                        movieInfo["time_stamp"] = datetime.datetime.now().isoformat();
                        mydb.movies.insert_one(movieInfo)
                        print movieInfo.title +" Inserted!"
                        calculate_movies_count();
                    if alreadyExists:
                        movie.stored = "true";
                    else:
                        movie.stored = "false";

            except Exception as e:
                print "Exception has occured"
                print str(e)
    return jsonify(response)

# Helpers

def count_iterable(i):
    return sum(1 for e in i)

def get_users_count():
    return user_count;

def get_movies_count():
    return count;

def calculate_movies_count():
    global count;
    count = math.ceil(mydb.movies.count() / float(batch_size));
    print count;

def calculate_users_count():
    global user_count;
    user_count = math.ceil(mydb.users.count() / float(batch_size));
    print user_count;

def computeMD5hash(string):
    m = hashlib.md5()
    m.update(string.encode('utf-8'))
    return m.hexdigest()

@app.route('/upload', methods=['GET','POST'])
def upload():
    if request.method == 'POST' and 'photo' in request.files:
        filename = photos.save(request.files['photo'])
        return jsonify({"message" : "success", "url" : photos.url(filename)})
    else:
        return jsonify({"message" : "invalid"})
@app.route('/activities', methods=['GET'])
def get_activites():
    username = None
    if request.args:
        username = request.args("username", None)
    result = []
    if username:
        activity_iter = mydb.activity.find({"user_id" : username}).sort("time_stamp",-1).limit(20);
        for activity in activity_iter:
            result.append(activity)
    else:
        activity_iter = mydb.activity.find().sort("time_stamp",-1).limit(20);
        for activity in activity_iter:
            activity["_id"] = "[HIDDEN]"
            result.append(activity)
    return jsonify({"activities" : result})

# Main
def isInArray(source, target, id1, id2):
    for sourceItem in source:
        if(sourceItem.get(id1) == target.get(id2)):
            return sourceItem
    else:
        return False

def createActivity(userImage, username, verb,moviename, rating):
    mydb.activity.insert_one({
        "profile_url" : userImage,
        "username" : username,
        "verb" : verb,
        "moviename" : moviename,
        "rating" : rating,
        "time_stamp" : datetime.datetime.now().isoformat()
    })
def transformMovie(movie):
    movie_genre = movie["genre"].split(", ")
    #Sci-Fi, Crime, Mystery, Thriller, Action, Comedy , Adventure, Drama, Horror, Family Short Animation rating, year, runtime
    movie_t = []
    genres = ["Sci-Fi", "Crime", "Thriller", "Action", "Comedy","Adventure","Drama","Horror","Family","Short","Animation"]
    for genre in genres:
        if(genre in movie_genre):
            movie_t.append(1)
        else:
            movie_t.append(0)

    if movie["imdb_rating"] != "N/A":
        movie_t.append(float(movie["imdb_rating"])/ 10)
    else:
        movie_t.append(0)
    movie_t.append(int(movie["year"]) / 2017)
    if movie["runtime"] != "N/A":
        movie_t.append(int(movie["runtime"].split(" ")[0])/ 160)
    else:
        movie_t.append(0)

    return movie_t

def formClusters(n_cl = 2,movies = []):
    movies_a = []
    for movie in movies:
        movies_a.append(transformMovie(movie))

    x = np.array(movies_a)

    kmeans = KMeans(n_clusters=n_cl, random_state=0).fit(x)
    return kmeans

if __name__ == '__main__':
    init_user_rating_matrix()
    #updateTransactions()
    app.run(host='0.0.0.0', debug = True)

def isInArray(source, target, id1, id2):
    print target
    for sourceItem in source:
        print sourceItem
        if(sourceItem.get(id1) == target.get(id2)):
            return sourceItem
    else:
        return False

max_stars_repo_name: tanmayv/movie-treat-node

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: tanmayv/movie-treat-node

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: tanmayv/movie-treat-node

max_forks_repo_licenses: MIT Unlicense


text:
#################################################################################
# The Institute for the Design of Advanced Energy Systems Integrated Platform
# Framework (IDAES IP) was produced under the DOE Institute for the
# Design of Advanced Energy Systems (IDAES), and is copyright (c) 2018-2021
# by the software owners: The Regents of the University of California, through
# Lawrence Berkeley National Laboratory,  National Technology & Engineering
# Solutions of Sandia, LLC, Carnegie Mellon University, West Virginia University
# Research Corporation, et al.  All rights reserved.
#
# Please see the files COPYRIGHT.md and LICENSE.md for full copyright and
# license information.
#################################################################################

__author__ = "John Eslick"

import pytest
from pyomo.environ import ConcreteModel, value
import idaes.generic_models.properties.swco2 as swco2
import idaes

class TestSWCO2(object):

    @pytest.fixture(scope="class")
    def model2(self):
        model = ConcreteModel()
        model.prop_param = swco2.SWCO2ParameterBlock()
        model.prop_in = swco2.SWCO2StateBlock(
            default={"parameters": model.prop_param}
        )
        return model

    @pytest.mark.unit
    def test_transport(self, model2):
        """Transport property tests
        """

        data = (
           #(T,      P,        mu,          tc,         tol_mu,  tol_tc)
            (200,    0.1e6,    10.06e-6,    9.63e-3,    0.02,    0.02),
            (240,    0.1e6,    12.07e-6,    12.23e-3,   0.02,    0.02),
            (240,    10.0e6,   188.91e-6,   159.24e-3,  0.02,    0.02),
            (240,    25e6,     214.16e-6,   171.59e-3,  0.02,    0.02),
            (280,    0.1e6,    14.05e-6,    15.19e-3,   0.02,    0.02),
            (280,    2.5e6,    14.51e-6,    17.46e-3,   0.02,    0.04),
            (280,    30.0e6,   134.98e-6,   136.63e-3,  0.02,    0.02),
            (280,    50.0e6,   160.73e-6,   152.67e-3,  0.02,    0.02),
            (306,    7.5e6,    24.06e-6,    23.93e-3,   0.02,    0.40), # supercritcal near critical
            (400,    5.0e6,    20.37e-6,    27.46e-3,   0.02,    0.02),
            (400,    50.0e6,   65.67e-6,    86.48e-3,   0.02,    0.02),
            (400,    100e6,    101.41e-6,   121.41e-3,  0.02,    0.02),
            (600,    50e6,     42.20e-6,    64.11e-3,   0.02,    0.02),
            (1000,   5e6,      41.42e-6,    71.18e-3,   0.02,    0.02),
            (1000,   50e6,     46.16e-6,    80.33e-3,   0.02,    0.02),
            (1000,   100e6,    54.79e-6,    92.16e-3,   0.02,    0.02),
        )

        for d in data:
            tol_mu = d[4]
            tol_tc = d[5]
            P = d[1]
            T = d[0]
            mu_data = d[2]
            tc_data = d[3]
            model2.prop_in.temperature.set_value(T)
            model2.prop_in.pressure = P
            Tsat = value(model2.prop_in.temperature_sat)
            if P >= 7.377e6 and T >= 304.1282:
                # super critical, which we clasify as liquid
                ph = "Liq"
            if P <= 0.5179e6:
                # below triple point
                ph = "Vap"
            elif T > Tsat + 0.5:
                ph = "Vap"
            elif T < Tsat - 0.5:
                ph = "Liq"
            else:
                # if too close to sat, don't want to get into a situation where
                # I don't know if it's liquid or vapor, so using extreme
                # caution
                # If we're woried about it, we can add tests on the sat curve,
                # and test both liquid and vapor
                continue

            rho = value(model2.prop_in.dens_mass_phase[ph])
            mu = value(model2.prop_in.visc_d_phase[ph])
            tc = value(model2.prop_in.therm_cond_phase[ph])
            print("T = {}, P = {}, mu = {}, tc = {}, rho = {}, phase = {}".
                  format(T, P, mu, tc, rho, ph)
                  )
            print("Data: mu = {}, tc = {}".format(mu_data, tc_data))
            assert tc == pytest.approx(tc_data, rel=tol_tc)
            assert mu == pytest.approx(mu_data, rel=tol_mu)

max_stars_repo_name: carldlaird/idaes-pse

max_issues_repo_name: carldlaird/idaes-pse

max_forks_repo_name: carldlaird/idaes-pse

reason: from pyomo


text:
from typing import List, Optional, Tuple

import numpy as np

from ..data import (
    filter_agents_by_labels,
    filter_tl_faces_by_frames,
    get_agents_slice_from_frames,
    get_tl_faces_slice_from_frames,
)
from ..data.filter import filter_agents_by_frames, filter_agents_by_track_id
from ..geometry import angular_distance, compute_agent_pose, rotation33_as_yaw, transform_points
from ..kinematic import Perturbation
from ..rasterization import EGO_EXTENT_HEIGHT, EGO_EXTENT_LENGTH, EGO_EXTENT_WIDTH, Rasterizer, RenderContext
from .slicing import get_future_slice, get_history_slice

def get_agent_context(
    state_index: int,
    frames: np.ndarray,
    agents: np.ndarray,
    tl_faces: np.ndarray,
    history_num_frames: int,
    future_num_frames: int,
) -> Tuple[np.ndarray, np.ndarray, List[np.ndarray], List[np.ndarray], List[np.ndarray], List[np.ndarray]]:
    """
    Slice zarr or numpy arrays to get the context around the agent onf interest (both in space and time)

    Args:
        state_index (int): frame index inside the scene
        frames (np.ndarray): frames from the scene
        agents (np.ndarray): agents from the scene
        tl_faces (np.ndarray): tl_faces from the scene
        history_num_frames (int): how many frames in the past to slice
        future_num_frames (int): how many frames in the future to slice

    Returns:
        Tuple[np.ndarray, np.ndarray, List[np.ndarray], List[np.ndarray], List[np.ndarray], List[np.ndarray]]
    """

    #  the history slice is ordered starting from the latest frame and goes backward in time., ex. slice(100, 91, -2)
    history_slice = get_history_slice(state_index, history_num_frames, 1, include_current_state=True)
    future_slice = get_future_slice(state_index, future_num_frames, 1)

    history_frames = frames[history_slice].copy()  # copy() required if the object is a np.ndarray
    future_frames = frames[future_slice].copy()

    sorted_frames = np.concatenate((history_frames[::-1], future_frames))  # from past to future

    # get agents (past and future)
    agent_slice = get_agents_slice_from_frames(sorted_frames[0], sorted_frames[-1])
    agents = agents[agent_slice].copy()
    # sync interval with the agents array
    history_frames["agent_index_interval"] -= agent_slice.start
    future_frames["agent_index_interval"] -= agent_slice.start
    history_agents = filter_agents_by_frames(history_frames, agents)
    future_agents = filter_agents_by_frames(future_frames, agents)

    # get traffic lights (past and future)
    tl_slice = get_tl_faces_slice_from_frames(sorted_frames[0], sorted_frames[-1])
    tl_faces = tl_faces[tl_slice].copy()
    # sync interval with the traffic light faces array
    history_frames["traffic_light_faces_index_interval"] -= tl_slice.start
    future_frames["traffic_light_faces_index_interval"] -= tl_slice.start
    history_tl_faces = filter_tl_faces_by_frames(history_frames, tl_faces)
    future_tl_faces = filter_tl_faces_by_frames(future_frames, tl_faces)

    return history_frames, future_frames, history_agents, future_agents, history_tl_faces, future_tl_faces

def compute_agent_velocity(
    history_positions_m: np.ndarray, future_positions_m: np.ndarray, step_time: float
) -> Tuple[np.ndarray, np.ndarray]:
    """
    Compute estimated velocities by finite differentiation on future positions as(pos(T+t) - pos(T))/t.
    This simple approach gives less than 0.5% velocity difference
    compared to (pos(T+t) - pos(T-t))/2t on v1.1/sample.zarr.tar.

    Args:
        history_positions_m (np.ndarray): history XY positions in meters
        future_positions_m (np.ndarray): future XY positions in meters
        step_time (np.ndarray): length of a step in second

    Returns:
        Tuple[np.ndarray, np.ndarray]: history and future XY speeds

    """
    assert step_time > np.finfo(float).eps, f"step_time must be greater then eps, got {step_time}"

    # [future_num_frames, 2]
    future_positions_diff_m = np.concatenate((future_positions_m[:1], np.diff(future_positions_m, axis=0)))
    # [future_num_frames, 2]
    future_vels_mps = np.float32(future_positions_diff_m / step_time)

    # current position is included in history positions
    # [history_num_frames, 2]
    history_positions_diff_m = np.diff(history_positions_m, axis=0)
    # [history_num_frames, 2]
    history_vels_mps = np.float32(history_positions_diff_m / step_time)

    return history_vels_mps, future_vels_mps

def get_relative_poses(
    num_frames: int,
    frames: np.ndarray,
    selected_track_id: Optional[int],
    agents: List[np.ndarray],
    agent_from_world: np.ndarray,
    current_agent_yaw: float,
) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
    """
    Internal function that creates the targets and availability masks for deep prediction-type models.
    The futures/history offset (in meters) are computed. When no info is available (e.g. agent not in frame)
    a 0 is set in the availability array (1 otherwise).

    Note: return dtype is float32, even if the provided args are float64. Still, the transformation
    between space is performed in float64 to ensure precision

    Args:
        num_frames (int): number of offset we want in the future/history
        frames (np.ndarray): available frames. This may be less than num_frames
        selected_track_id (Optional[int]): agent track_id or AV (None)
        agents (List[np.ndarray]): list of agents arrays (same len of frames)
        agent_from_world (np.ndarray): local from world matrix
        current_agent_yaw (float): angle of the agent at timestep 0

    Returns:
        Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: position offsets, angle offsets, extent, availabilities

    """
    # How much the coordinates differ from the current state in meters.
    positions_m = np.zeros((num_frames, 2), dtype=agent_from_world.dtype)
    yaws_rad = np.zeros((num_frames, 1), dtype=np.float32)
    extents_m = np.zeros((num_frames, 2), dtype=np.float32)
    availabilities = np.zeros((num_frames,), dtype=np.float32)

    for i, (frame, frame_agents) in enumerate(zip(frames, agents)):
        if selected_track_id is None:
            agent_centroid_m = frame["ego_translation"][:2]
            agent_yaw_rad = rotation33_as_yaw(frame["ego_rotation"])
            agent_extent = (EGO_EXTENT_LENGTH, EGO_EXTENT_WIDTH)
        else:
            # it's not guaranteed the target will be in every frame
            try:
                agent = filter_agents_by_track_id(frame_agents, selected_track_id)[0]
                agent_centroid_m = agent["centroid"]
                agent_yaw_rad = agent["yaw"]
                agent_extent = agent["extent"][:2]
            except IndexError:
                availabilities[i] = 0.0  # keep track of invalid futures/history
                continue

        positions_m[i] = agent_centroid_m
        yaws_rad[i] = agent_yaw_rad
        extents_m[i] = agent_extent
        availabilities[i] = 1.0

    # batch transform to speed up
    positions_m = transform_points(positions_m, agent_from_world) * availabilities[:, np.newaxis]
    yaws_rad = angular_distance(yaws_rad, current_agent_yaw) * availabilities[:, np.newaxis]
    return positions_m.astype(np.float32), yaws_rad, extents_m, availabilities

def generate_agent_sample(
    state_index: int,
    frames: np.ndarray,
    agents: np.ndarray,
    tl_faces: np.ndarray,
    selected_track_id: Optional[int],
    render_context: RenderContext,
    history_num_frames: int,
    future_num_frames: int,
    step_time: float,
    filter_agents_threshold: float,
    rasterizer: Optional[Rasterizer] = None,
    perturbation: Optional[Perturbation] = None,
) -> dict:
    """Generates the inputs and targets to train a deep prediction model. A deep prediction model takes as input
    the state of the world (here: an image we will call the "raster"), and outputs where that agent will be some
    seconds into the future.

    This function has a lot of arguments and is intended for internal use, you should try to use higher level classes
    and partials that use this function.

    Args:
        state_index (int): The anchor frame index, i.e. the "current" timestep in the scene
        frames (np.ndarray): The scene frames array, can be numpy array or a zarr array
        agents (np.ndarray): The full agents array, can be numpy array or a zarr array
        tl_faces (np.ndarray): The full traffic light faces array, can be numpy array or a zarr array
        selected_track_id (Optional[int]): Either None for AV, or the ID of an agent that you want to
        predict the future of. This agent is centered in the raster and the returned targets are derived from
        their future states.
        render_context (RenderContext): The context for rasterisation
        history_num_frames (int): Amount of history frames to draw into the rasters
        future_num_frames (int): Amount of history frames to draw into the rasters
        step_time (float): seconds between consecutive steps
        filter_agents_threshold (float): Value between 0 and 1 to use as cutoff value for agent filtering
        based on their probability of being a relevant agent
        rasterizer (Optional[Rasterizer]): Rasterizer of some sort that draws a map image
        perturbation (Optional[Perturbation]): Object that perturbs the input and targets, used
        to train models that can recover from slight divergence from training set data

    Raises:
        IndexError: An IndexError is returned if the specified ``selected_track_id`` is not present in the scene
        or was filtered by applying the ``filter_agent_threshold`` probability filtering.

    Returns:
        dict: a dict object with the raster array, the future offset coordinates (meters),
        the future yaw angular offset, the future_availability as a binary mask
    """
    (
        history_frames,
        future_frames,
        history_agents,
        future_agents,
        history_tl_faces,
        future_tl_faces,
    ) = get_agent_context(state_index, frames, agents, tl_faces, history_num_frames, future_num_frames,)

    if perturbation is not None:
        history_frames, future_frames = perturbation.perturb(
            history_frames=history_frames, future_frames=future_frames
        )

    # State you want to predict the future of.
    cur_frame = history_frames[0]
    cur_agents = history_agents[0]

    if selected_track_id is None:
        agent_centroid_m = cur_frame["ego_translation"][:2]
        agent_yaw_rad = rotation33_as_yaw(cur_frame["ego_rotation"])
        agent_extent_m = np.asarray((EGO_EXTENT_LENGTH, EGO_EXTENT_WIDTH, EGO_EXTENT_HEIGHT))
        selected_agent = None
    else:
        # this will raise IndexError if the agent is not in the frame or under agent-threshold
        # this is a strict error, we cannot recover from this situation
        filtered_agents = filter_agents_by_labels(cur_agents, filter_agents_threshold)
        agent = filter_agents_by_track_id(filtered_agents, selected_track_id)[0]
        agent_centroid_m = agent["centroid"]
        agent_yaw_rad = float(agent["yaw"])
        agent_extent_m = agent["extent"]
        selected_agent = agent

    input_im = (
        None
        if not rasterizer
        else rasterizer.rasterize(history_frames, history_agents, history_tl_faces, selected_agent)
    )

    world_from_agent = compute_agent_pose(agent_centroid_m, agent_yaw_rad)
    agent_from_world = np.linalg.inv(world_from_agent)
    raster_from_world = render_context.raster_from_world(agent_centroid_m, agent_yaw_rad)

    future_positions_m, future_yaws_rad, future_extents, future_availabilities = get_relative_poses(
        future_num_frames, future_frames, selected_track_id, future_agents, agent_from_world, agent_yaw_rad,
    )
    # history_num_frames + 1 because it also includes the current frame
    history_positions_m, history_yaws_rad, history_extents, history_availabilities = get_relative_poses(
        history_num_frames + 1, history_frames, selected_track_id, history_agents, agent_from_world, agent_yaw_rad,
    )

    history_vels_mps, future_vels_mps = compute_agent_velocity(history_positions_m, future_positions_m, step_time)

    return {
        "frame_index": state_index,
        "image": input_im,
        "target_positions": future_positions_m,
        "target_yaws": future_yaws_rad,
        "target_velocities": future_vels_mps,
        "target_availabilities": future_availabilities,
        "history_positions": history_positions_m,
        "history_yaws": history_yaws_rad,
        "history_velocities": history_vels_mps,
        "history_availabilities": history_availabilities,
        "world_to_image": raster_from_world,  # TODO deprecate
        "raster_from_agent": raster_from_world @ world_from_agent,
        "raster_from_world": raster_from_world,
        "agent_from_world": agent_from_world,
        "world_from_agent": world_from_agent,
        "centroid": agent_centroid_m,
        "yaw": agent_yaw_rad,
        "speed": np.linalg.norm(future_vels_mps[0]),
        "extent": agent_extent_m,
        "history_extents": history_extents,
        "future_extents": future_extents,
    }


text:
import bisect
from copy import deepcopy
from fractions import Fraction
from functools import reduce
import heapq as hq
import io
from itertools import combinations, permutations
import math
from math import factorial
import re
import sys
sys.setrecursionlimit(10000)

#from numba import njit
import numpy as np

_INPUT_1 = """\
12
"""
_INPUT_2 = """\
100128
"""
_INPUT_3 = """\
1 2 3
"""

def solve():
    N = int(input())
    SMALL = 1
    LARGE = 1000000000
    while True:
        i = (SMALL+LARGE) // 2
        v = ((1+i)*i)//2
        if v == N:
            day = i-
        elif v > N:
            LARGE = i
        elif v < N:

if __file__ != './Main.py':
    if '_INPUT_1' in globals():
        sys.stdin = io.StringIO(_INPUT_1)
        solve()
    if '_INPUT_2' in globals():
        sys.stdin = io.StringIO(_INPUT_2)
        solve()
    if '_INPUT_3' in globals():
        sys.stdin = io.StringIO(_INPUT_3)
        solve()
else:
    solve()

reason: import numpy,from numba


text:
import numpy as np
from gym.spaces import Box

from metaworld.envs import reward_utils
from metaworld.envs.asset_path_utils import full_v2_path_for
from metaworld.envs.mujoco.sawyer_xyz.sawyer_xyz_env import SawyerXYZEnv, _assert_task_is_set
from scipy.spatial.transform import Rotation

class SawyerPegInsertionSideEnvV2(SawyerXYZEnv):
    TARGET_RADIUS = 0.07
    """
    Motivation for V2:
        V1 was difficult to solve because the observation didn't say where
        to insert the peg (the hole's location). Furthermore, the hole object
        could be initialized in such a way that it severely restrained the
        sawyer's movement.
    Changelog from V1 to V2:
        - (8/21/20) Updated to Byron's XML
        - (7/7/20) Removed 1 element vector. Replaced with 3 element position
            of the hole (for consistency with other environments)
        - (6/16/20) Added a 1 element vector to the observation. This vector
            points from the end effector to the hole in the Y direction.
            i.e. (self._target_pos - pos_hand)[1]
        - (6/16/20) Used existing goal_low and goal_high values to constrain
            the hole's position, as opposed to hand_low and hand_high
    """
    def __init__(self):
        hand_init_pos = (0, 0.6, 0.2)

        hand_low = (-0.5, 0.40, 0.05)
        hand_high = (0.5, 1, 0.5)
        obj_low = (.0, 0.5, 0.02)
        obj_high = (.2, 0.7, 0.02)
        goal_low = (-0.35, 0.4, -0.001)
        goal_high = (-0.25, 0.7, 0.001)

        super().__init__(
            self.model_name,
            hand_low=hand_low,
            hand_high=hand_high,
        )

        self.init_config = {
            'obj_init_pos': np.array([0, 0.6, 0.02]),
            'hand_init_pos': np.array([0, .6, .2]),
        }

        self.goal = np.array([-0.3, 0.6, 0.0])

        self.obj_init_pos = self.init_config['obj_init_pos']
        self.hand_init_pos = self.init_config['hand_init_pos']

        self.hand_init_pos = np.array(hand_init_pos)

        self._random_reset_space = Box(
            np.hstack((obj_low, goal_low)),
            np.hstack((obj_high, goal_high)),
        )
        self.goal_space = Box(
            np.array(goal_low) + np.array([.03, .0, .13]),
            np.array(goal_high) + np.array([.03, .0, .13])
        )

    @property
    def model_name(self):
        return full_v2_path_for('sawyer_xyz/sawyer_peg_insertion_side.xml')

    @_assert_task_is_set
    def evaluate_state(self, obs, action):
        obj = obs[4:7]

        reward, tcp_to_obj, tcp_open, obj_to_target, grasp_reward, in_place_reward, collision_box_front, ip_orig= (
            self.compute_reward(action, obs))
        grasp_success = float(tcp_to_obj < 0.02 and (tcp_open > 0) and (obj[2] - 0.01 > self.obj_init_pos[2]))
        success = float(obj_to_target <= 0.07)
        near_object = float(tcp_to_obj <= 0.03)

        info = {
            'success': success,
            'near_object': near_object,
            'grasp_success': grasp_success,
            'grasp_reward': grasp_reward,
            'in_place_reward': in_place_reward,
            'obj_to_target': obj_to_target,
            'unscaled_reward': reward,
        }

        return reward, info

    def _get_pos_objects(self):
        return self._get_site_pos('pegGrasp')

    def _get_quat_objects(self):
        return Rotation.from_matrix(self.data.get_site_xmat('pegGrasp')).as_quat()

    def reset_model(self):
        self._reset_hand()

        pos_peg = self.obj_init_pos
        pos_box = self.goal
        if self.random_init:
            pos_peg, pos_box = np.split(self._get_state_rand_vec(), 2)
            while np.linalg.norm(pos_peg[:2] - pos_box[:2]) < 0.1:
                pos_peg, pos_box = np.split(self._get_state_rand_vec(), 2)

        self.obj_init_pos = pos_peg
        self.peg_head_pos_init = self._get_site_pos('pegHead')
        self._set_obj_xyz(self.obj_init_pos)

        self.sim.model.body_pos[self.model.body_name2id('box')] = pos_box
        self._target_pos = pos_box + np.array([.03, .0, .13])

        return self._get_obs()

    def compute_reward(self, action, obs):
        tcp = self.tcp_center
        obj = obs[4:7]
        obj_head = self._get_site_pos('pegHead')
        tcp_opened = obs[3]
        target = self._target_pos
        tcp_to_obj = np.linalg.norm(obj - tcp)
        scale = np.array([1., 2., 2.])
        #  force agent to pick up object then insert
        obj_to_target = np.linalg.norm((obj_head - target) * scale)

        in_place_margin = np.linalg.norm((self.peg_head_pos_init - target) * scale)
        in_place = reward_utils.tolerance(obj_to_target,
                                    bounds=(0, self.TARGET_RADIUS),
                                    margin=in_place_margin,
                                    sigmoid='long_tail',)
        ip_orig = in_place
        brc_col_box_1 = self._get_site_pos('bottom_right_corner_collision_box_1')
        tlc_col_box_1 = self._get_site_pos('top_left_corner_collision_box_1')

        brc_col_box_2 = self._get_site_pos('bottom_right_corner_collision_box_2')
        tlc_col_box_2 = self._get_site_pos('top_left_corner_collision_box_2')
        collision_box_bottom_1 = reward_utils.rect_prism_tolerance(curr=obj_head,
                                                                   one=tlc_col_box_1,
                                                                   zero=brc_col_box_1)
        collision_box_bottom_2 = reward_utils.rect_prism_tolerance(curr=obj_head,
                                                                   one=tlc_col_box_2,
                                                                   zero=brc_col_box_2)
        collision_boxes = reward_utils.hamacher_product(collision_box_bottom_2,
                                                        collision_box_bottom_1)
        in_place = reward_utils.hamacher_product(in_place,
                                                 collision_boxes)

        pad_success_margin = 0.03
        object_reach_radius=0.01
        x_z_margin = 0.005
        obj_radius = 0.0075

        object_grasped = self._gripper_caging_reward(action,
                                                     obj,
                                                     object_reach_radius=object_reach_radius,
                                                     obj_radius=obj_radius,
                                                     pad_success_thresh=pad_success_margin,
                                                     xz_thresh=x_z_margin,
                                                     high_density=True)
        if tcp_to_obj < 0.08 and (tcp_opened > 0) and (obj[2] - 0.01 > self.obj_init_pos[2]):
            object_grasped = 1.
        in_place_and_object_grasped = reward_utils.hamacher_product(object_grasped,
                                                                    in_place)
        reward = in_place_and_object_grasped

        if tcp_to_obj < 0.08 and (tcp_opened > 0) and (obj[2] - 0.01 > self.obj_init_pos[2]):
            reward += 1. + 5 * in_place

        if obj_to_target <= 0.07:
            reward = 10.

        return [reward, tcp_to_obj, tcp_opened, obj_to_target, object_grasped, in_place, collision_boxes, ip_orig]

max_stars_repo_name: yiwc/robotics-world

max_issues_repo_name: yiwc/robotics-world

max_forks_repo_name: yiwc/robotics-world