real scores pending

.github/workflows/verify.yaml

@@ -0,0 +1,20 @@
+name: Verify
+
+on:
+  push:
+    branches: [ "*" ]
+
+jobs:
+  build-backend:
+    runs-on: ubuntu-latest
+    steps:
+      - name: Check out the repository
+        uses: actions/checkout@v5
+      - name: Set up Python 3.10
+        uses: actions/setup-python@v6
+        with:
+          python-version: '3.10'
+      - name: Install dependencies
+        run: |
+          python -m pip install --upgrade pip
+          pip install -r requirements.txt

.gitignore

@@ -1,3 +1,4 @@
 **/.idea
 **/.venv
-**/venv
+**/venv
+**/local-run-data

app.py

@@ -1,21 +1,39 @@
 import datetime
-import random
+import os.path
+import sys
 import uuid
 from os import PathLike
 
 import gradio as gr
 import pandas as pd
+import torch
 
 from config import APP_CONFIG
 from data_repository import REPOSITORY_INSTANCE, ModelScoringResult
 from designs_submission_validations import validate_github_link, validate_user_designs
+from domain_constants import SCORE_NAMES_MAP, USER_GEN_DESIGNS_COLUMNS
+
+sys.path.append(os.path.join(os.path.dirname(__file__), "bike_bench_internal/src/"))
+from bikebench.benchmarking import benchmarking_utils
+
 
 def compute_scores(user_gen_designs: pd.DataFrame) -> ModelScoringResult:
+    user_gen_designs = pd.DataFrame(user_gen_designs, columns=USER_GEN_DESIGNS_COLUMNS)
+    designs_length = len(user_gen_designs)
+    if designs_length < 10_000:
+        raise Exception(f"Too few designs to evaluate. Expected > 10,000, got {designs_length}")
+    data_tens = torch.tensor(user_gen_designs.values, dtype=torch.float32)
+    main_scores, detailed_scores, all_evaluation_scores = benchmarking_utils.evaluate(data_tens, device="cpu",
+                                                                                      evaluate_as_aggregate=False)
     return ModelScoringResult(
-        score=random.randint(50, 5000),
-        scoring_time=datetime.datetime.now(),
+        uuid=str(uuid.uuid4()),
         submission_time=datetime.datetime.now(),
-        uuid=str(uuid.uuid4())
+        design_quality=main_scores[SCORE_NAMES_MAP["design_quality"]],
+        diversity_dpp=main_scores[SCORE_NAMES_MAP["diversity_dpp"]],
+        mean_novelty=main_scores[SCORE_NAMES_MAP["mean_novelty"]],
+        sim_to_data_mmd=main_scores[SCORE_NAMES_MAP["sim_to_data_mmd"]],
+        mean_violations=main_scores[SCORE_NAMES_MAP["mean_violations"]],
+        binary_validity=main_scores[SCORE_NAMES_MAP["binary_validity"]],
     )
 
 
@@ -29,21 +47,28 @@ def process_generated_designs(github_link: str, file: PathLike[str]):
     return f"File uploaded successfully, uuid {scores.uuid}"
 
 
-with gr.Blocks() as gradio_app:
-    with gr.Tab("Bike Bench Leaderboard"):
-        gr.Markdown("Hello beautiful people!")
-        gr.Dataframe(REPOSITORY_INSTANCE.get_data_to_display, label="Scores of Previous Files")
-
-    with gr.Tab("Upload File"):
-        gr.Interface(
-            fn=process_generated_designs,
-            inputs=[
-                gr.Textbox(label="Github Link"),
-                gr.File(label="Upload a file"),
-            ],
-            outputs="text",
-            title="Bike Bench Leaderboard",
-            description="Upload a file to see the result."
-        )
-
-gradio_app.launch(debug=(not APP_CONFIG.production))
+def build_approval_app():
+    pass
+
+
+def build_app():
+    with gr.Blocks() as gradio_app:
+        with gr.Tab("Bike Bench Leaderboard"):
+            gr.Markdown("Hello beautiful people!")
+            gr.Dataframe(REPOSITORY_INSTANCE.get_data_to_display, label="Scores of Previous Files")
+
+        with gr.Tab("Upload File"):
+            gr.Interface(
+                fn=process_generated_designs,
+                inputs=[
+                    gr.Textbox(label="Github Link"),
+                    gr.File(label="Upload a file"),
+                ],
+                outputs="text",
+                title="Bike Bench Leaderboard",
+                description="Upload a file to see the result."
+            )
+    return gradio_app
+
+
+build_app().launch(debug=(not APP_CONFIG.production))

bike_bench_internal/.gitignore

@@ -0,0 +1,9 @@
+**/__pycache__/
+__pycache__/
+src/resources/datasets/Generative_Modeling_Datasets/
+src/resources/datasets/Predictive_Modeling_Datasets/
+src/resources/datasets/Original_BIKED_Data/
+src/resources/datasets/Real_Extended_Data/
+src/resources/datasets/Synthetic_Extended_Data/
+**/BikeCAD_17.1_configuration/
+src/bikebench/introductory_notebooks/frames

bike_bench_internal/README.md

@@ -0,0 +1 @@
+# Bike-Bench-Internal

bike_bench_internal/benchmark_models/benchmarking_utils.py

@@ -0,0 +1,201 @@
+import os
+import torch
+import pandas as pd
+from bikebench.design_evaluation.scoring import construct_scorer, MainScores, DetailedScores
+from bikebench.design_evaluation.design_evaluation import get_standard_evaluations
+from bikebench.conditioning import conditioning
+from tqdm import trange, tqdm
+
+
+
+def get_condition_by_idx(idx=0):
+    rider_condition = conditioning.sample_riders(10, split="test")
+    use_case_condition = conditioning.sample_use_case(10, split="test")
+    image_embeddings = conditioning.sample_image_embedding(10, split="test")
+    condition = {"Rider": rider_condition[idx], "Use Case": use_case_condition[idx], "Embedding": image_embeddings[idx]}
+    return condition
+
+def get_conditions_10k():
+    rider_condition = conditioning.sample_riders(10000, split="test")
+    use_case_condition = conditioning.sample_use_case(10000, split="test")
+    image_embeddings = conditioning.sample_image_embedding(10000, split="test")
+    conditions = {"Rider": rider_condition, "Use Case": use_case_condition, "Embedding": image_embeddings}
+    return conditions
+
+def evaluate_uncond(result_tens, name, cond_idx, data_columns, device, save=True):
+
+    condition = get_condition_by_idx(cond_idx)
+
+    result_dir = os.path.join("results", "unconditional", f"cond_{cond_idx}", name)
+    os.makedirs(result_dir, exist_ok=True)
+    
+    main_scorer = construct_scorer(MainScores, get_standard_evaluations(device), data_columns)
+    detailed_scorer = construct_scorer(DetailedScores, get_standard_evaluations(device), data_columns)
+
+    main_scores = main_scorer(result_tens, condition)
+    
+    detailed_scores = detailed_scorer(result_tens, condition)
+    
+    if save:
+        result_tens = result_tens.cpu()
+        torch.save(result_tens, os.path.join(result_dir, "result_tens.pt"))
+        main_scores.to_csv(os.path.join(result_dir, "main_scores.csv"), index_label=False, header=False)
+        detailed_scores.to_csv(os.path.join(result_dir, "detailed_scores.csv"), index_label=False, header=False)
+    return main_scores, detailed_scores
+
+def evaluate_cond(result_tens, name, data_columns, device, save=True):
+    condition = get_conditions_10k()
+
+    condition = {"Rider": condition["Rider"], "Use Case": condition["Use Case"], "Embedding": condition["Embedding"]}
+
+    result_dir = os.path.join("results", "conditional", name)
+    os.makedirs(result_dir, exist_ok=True)
+
+    main_scorer = construct_scorer(MainScores, get_standard_evaluations(device), data_columns, device)
+    detailed_scorer = construct_scorer(DetailedScores, get_standard_evaluations(device), data_columns, device)
+
+    main_scores = main_scorer(result_tens, condition)
+    detailed_scores = detailed_scorer(result_tens, condition)
+
+    if save:
+        result_tens = result_tens.cpu()
+        torch.save(result_tens, os.path.join(result_dir, "result_tens.pt"))
+        main_scores.to_csv(os.path.join(result_dir, "main_scores.csv"), index_label=False, header=False)
+        detailed_scores.to_csv(os.path.join(result_dir, "detailed_scores.csv"), index_label=False, header=False)
+
+    return main_scores, detailed_scores
+
+
+def create_score_report_conditional():
+    """
+    Looks through the results folder and creates a score report for each conditional result.
+    """
+    all_scores = []
+    result_dir = os.path.join("results", "conditional")
+    for name in os.listdir(result_dir):
+        if os.path.isdir(os.path.join(result_dir, name)):
+            main_scores = pd.read_csv(os.path.join(result_dir, name, "main_scores.csv"), header=None)
+            main_scores.columns = ["Metric", "Score"]
+            main_scores["Model"] = name
+            all_scores.append(main_scores)
+    all_scores = pd.concat(all_scores, axis=0)
+    #make metric names the three columns, make models the rows
+    all_scores = all_scores.pivot(index="Model", columns="Metric", values="Score")
+    #drop the index name and the column name
+    all_scores.columns.name = None
+    all_scores.index.name = None
+    
+    return all_scores
+
+def create_score_report_unconditional():
+    """
+    Looks through the results folder and creates a score report for each unconditional result.
+    """
+    all_scores = []
+    result_dir = os.path.join("results", "unconditional")
+    for i in range(10):
+        c_dir = os.path.join(result_dir, f"cond_{i}")
+        for name in os.listdir(c_dir):
+            dirname = os.path.join(c_dir, name)
+            if os.path.isdir(dirname):
+                main_scores = pd.read_csv(os.path.join(dirname, "main_scores.csv"), header=None)
+                main_scores.columns = ["Metric", "Score"]
+                main_scores["Model"] = name
+                main_scores["Condition"] = i
+                all_scores.append(main_scores)
+    all_scores = pd.concat(all_scores, axis=0)
+    #average over condition 
+    all_scores = all_scores.groupby(["Model", "Metric"]).mean().reset_index()
+    #make metric names the three columns, make models the rows
+    all_scores = all_scores.pivot(index="Model", columns="Metric", values="Score")
+    #drop the index name and the column name
+    all_scores.columns.name = None
+    all_scores.index.name = None
+    return all_scores
+
+
+def rescore_unconditional(data_columns, device, cond_idxs = None, model_names = None, results_root="results/unconditional"):
+    """
+    Recompute main and detailed scores for all unconditional results.
+    Overwrites only the CSV score files, leaves result_tens.pt untouched.
+    """
+
+    evals = get_standard_evaluations(device)
+    main_scorer     = construct_scorer(MainScores,    evals, data_columns, device)
+    detailed_scorer = construct_scorer(DetailedScores, evals, data_columns, device)
+    device = torch.device(device)
+    if cond_idxs is None:
+        cond_idxs = range(10)
+    for cond_idx in tqdm(cond_idxs):
+        cond_dir = os.path.join(results_root, f"cond_{cond_idx}")
+        if not os.path.isdir(cond_dir):
+            continue
+        # fetch the one shared condition for this index
+        condition = get_condition_by_idx(cond_idx)
+
+        if model_names is not None:
+            models = model_names
+        else:
+            models = os.listdir(cond_dir)
+
+        for model_name in models:
+            model_dir = os.path.join(cond_dir, model_name)
+            tensor_path = os.path.join(model_dir, "result_tens.pt")
+            if not os.path.isdir(model_dir) or not os.path.isfile(tensor_path):
+                continue
+
+            # load results
+            result_tens = torch.load(tensor_path, map_location=device)
+
+            # rescore
+            main_scores     = main_scorer(result_tens, condition)
+            detailed_scores = detailed_scorer(result_tens, condition)
+
+            # overwrite only the CSVs
+            main_scores.to_csv(
+                os.path.join(model_dir, "main_scores.csv"), header=False
+            )
+            detailed_scores.to_csv(
+                os.path.join(model_dir, "detailed_scores.csv"), header=False
+            )
+
+
+def rescore_conditional(data_columns, device, model_names, results_root="results/conditional"):
+    """
+    Recompute main and detailed scores for all conditional results.
+    Overwrites only the CSV score files, leaves result_tens.pt untouched.
+    """
+    device = torch.device(device)
+    # fetch the full 10k‐point condition set once
+    condition = get_conditions_10k()
+
+    # build scorers
+    evals = get_standard_evaluations(device)
+    main_scorer     = construct_scorer(MainScores, evals, data_columns, device)
+    detailed_scorer = construct_scorer(DetailedScores, evals, data_columns, device)
+
+
+    if model_names is not None:
+        models = model_names
+    else:
+        models = os.listdir(results_root)
+    for model_name in models:
+        model_dir  = os.path.join(results_root, model_name)
+        tensor_path = os.path.join(model_dir, "result_tens.pt")
+        if not os.path.isdir(model_dir) or not os.path.isfile(tensor_path):
+            continue
+
+        # load results
+        result_tens = torch.load(tensor_path, map_location=device)
+
+        # rescore
+        main_scores     = main_scorer(result_tens, condition)
+        detailed_scores = detailed_scorer(result_tens, condition)
+
+        # overwrite only the CSVs
+        main_scores.to_csv(
+            os.path.join(model_dir, "main_scores.csv"), header=False
+        )
+        detailed_scores.to_csv(
+            os.path.join(model_dir, "detailed_scores.csv"), header=False
+        )

bike_bench_internal/benchmark_models/generative_modeling_utils.py

@@ -0,0 +1,641 @@
+from tqdm import tqdm, trange
+import torch
+import torch.nn as nn
+from torch.utils.data import DataLoader, Dataset, TensorDataset
+import numpy as np
+from torch.autograd import grad
+from diffusers import DDPMScheduler
+from torch.nn import MSELoss
+
+
+from bikebench.conditioning import conditioning
+from bikebench.design_evaluation.design_evaluation import *
+from bikebench.conditioning import conditioning
+from bikebench.design_evaluation import scoring
+
+class TorchScaler:
+    def __init__(self, data):
+        self.data = data
+        self.mean = torch.mean(data, dim=0)
+        self.std = torch.std(data, dim=0)
+
+    def scale(self, x):
+        return (x - self.mean) / self.std
+
+    def unscale(self, x):
+        return x * self.std + self.mean
+
+def sample_continuous(num_samples, split="test", randomize = False):
+    emb = conditioning.sample_image_embedding(num_samples, split, randomize)
+    rider = conditioning.sample_riders(num_samples, split, randomize)
+    use_case = conditioning.sample_use_case(num_samples, split, randomize)
+    all = torch.cat((emb, rider, use_case), dim=1)
+    return all
+
+def sample_continuous_text(text_strings, device, split="test", randomize = False):
+    embedding_model = clip_embedding_calculator.ClipEmbeddingCalculator(
+            device=device, batch_size=64
+        )
+    text_embedding = embedding_model.embed_texts(text_strings)
+    num_samples = text_embedding.shape[0]
+    rider = conditioning.sample_riders(num_samples, split, randomize)
+    use_case = conditioning.sample_use_case(num_samples, split, randomize)
+    all = torch.cat((text_embedding, rider, use_case), dim=1)
+    return all
+
+def sample_standard(num_samples, split="test", randomize = False):
+    emb = conditioning.sample_image_embedding(num_samples, split, randomize)
+    rider = conditioning.sample_riders(num_samples, split, randomize)
+    use_case = conditioning.sample_use_case(num_samples, split, randomize)
+    condition = {"Rider": rider, "Use Case": use_case, "Embedding": emb}
+    return condition
+
+def parse_continuous_condition(condition):
+    image_embeddings = condition[:, :512]
+    use_case_condition = condition[:, -3:]
+    rider_condition = condition[:, 512:-3]
+    condition = {"Rider": rider_condition, "Use Case": use_case_condition, "Embedding": image_embeddings}
+    return condition
+
+def piecewise_constraint_score(constraint_scores, constraint_falloff = 10):
+    constraint_scores_safexp = torch.clamp(constraint_scores, max=0.0)
+    piece1 = torch.exp(constraint_scores_safexp * constraint_falloff)/constraint_falloff
+    piece2 = constraint_scores + 1/constraint_falloff
+    mask = constraint_scores < 0.0
+    mask = mask.float()
+    result = piece1 * mask + piece2 * (1 - mask)
+    return result
+
+def get_composite_score_fn(scaler, columns, evaluations, constrant_vs_objective_weight = 10.0, constraint_falloff=10.0, device="cpu"):
+    evaluator, requirement_names, requirement_types = construct_tensor_evaluator(evaluations, columns)
+
+    isobjective = torch.tensor(requirement_types) == 1
+
+    weights = scoring.get_ref_point(evaluator, requirement_names, requirement_names, reduction="meanabs", device=device)
+    weights = torch.tensor(weights, dtype=torch.float32, device=device)
+
+    assert weights.min() > 0, "Ref point should be greater than 0"
+
+    def composite_score_fn(x, continuous_condition, evaluator = evaluator, scaler = scaler):
+        #print if there are any NaN values in x
+        if torch.isnan(x).any():
+            print("NaN values in x")
+        x = scaler.unscale(x)
+        condition = parse_continuous_condition(continuous_condition)
+        eval_scores = evaluator(x, condition)
+        scaled_scores = eval_scores / weights
+        objective_scores = scaled_scores[:, isobjective]
+        constraint_scores_raw = scaled_scores[:, ~isobjective]
+        constraint_scores = piecewise_constraint_score(constraint_scores_raw, constraint_falloff)
+        total_scores = torch.sum(objective_scores, dim=1) + torch.sum(constraint_scores, dim=1) * constrant_vs_objective_weight
+        composite_scores = total_scores / constrant_vs_objective_weight
+
+        quality_scores = 1/composite_scores
+        # print("Quality scores: ", quality_scores)
+        if torch.isnan(quality_scores).any():
+            print("NaN values in quality scores")
+
+        mean_comp_scores = torch.mean(composite_scores)
+        constraint_satisfaction_rate = torch.mean(torch.all(constraint_scores_raw <= 0, dim=1).float())
+        report = {"CSR": constraint_satisfaction_rate, "MCS": mean_comp_scores}
+        return quality_scores, report
+
+    return composite_score_fn
+
+def get_uneven_batch_sizes(total_data_points, batch_size):
+    """
+    Given the total number of data points and a target batch size, 
+    returns a list of batch sizes that sum up to the total number of data points.
+    The batch sizes are distributed as evenly as possible but may be uneven.
+
+    :param total_data_points: Total number of data points (int).
+    :param batch_size: Target batch size (int).
+    :return: A list of batch sizes (list of int).
+    """
+    # Calculate the number of batches needed
+    num_batches = total_data_points // batch_size
+    remainder = total_data_points % batch_size
+
+    # Initialize the batch sizes
+    batch_sizes = [batch_size] * num_batches
+
+    # Distribute the remainder across the batches
+    for i in range(remainder):
+        batch_sizes[i%num_batches] += 1
+
+    return batch_sizes
+
+
+def get_diversity_loss_fn(scaler:TorchScaler, columns, evaluations, diversity_weight=0.1, score_weight=0.1, constraint_vs_objective_weight=10.0, constraint_falloff=10.0, dpp_batch=16, device="cpu"):
+    composite_score_fn = get_composite_score_fn(scaler, columns, evaluations, constrant_vs_objective_weight = constraint_vs_objective_weight, constraint_falloff = constraint_falloff, device=device)
+
+    def diversity_loss_fn(x, condition, diversity_weight=diversity_weight, score_weight=score_weight):
+        if diversity_weight == 0:
+            return torch.tensor(0.0), {"DIV-OFF": 0.0}
+        scores, report= composite_score_fn(x, condition)
+
+        # Initialize the total loss
+        total_loss = 0.0
+
+        # Get uneven batch sizes based on the total number of data points
+        batch_sizes = get_uneven_batch_sizes(x.size(0), dpp_batch)
+
+        # Split the data into uneven batches
+        start_idx = 0
+        for batch_size in batch_sizes:
+            # Get the current batch
+            end_idx = start_idx + batch_size
+            x_batch = x[start_idx:end_idx]
+            scores_batch = scores[start_idx:end_idx]
+            # Compute pairwise squared Euclidean distances for the batch
+            r = torch.sum(x_batch ** 2, dim=1, keepdim=True)
+            D = r - 2 * torch.matmul(x_batch, x_batch.T) + r.T
+            D_norm = D / x_batch.size(1)  # Normalize by the number of features
+            # Compute the similarity matrix using RBF for the batch
+            S = torch.exp(-0.5 * D_norm ** 2) / 2
+
+            # Compute the quality matrix for the batch
+            Q = torch.matmul(scores_batch, scores_batch.T)
+            Q = torch.pow(Q, score_weight)
+            L = S * Q
+
+            L = (L + L.T) / 2.0
+
+            L_stable = L + 1e-6 * torch.eye(L.size(0), device=L.device)  
+
+            # Compute the eigenvalues of the similarity matrix for the batch
+            try:
+                eig_val = torch.linalg.eigvalsh(L_stable)
+            except:
+                print(f"Eigenvalue computation failed for batch with size {batch_size}")
+                eig_val = torch.ones(x_batch.size(0), device=x.device)
+            if torch.isnan(eig_val).any():
+                print("NaNs detected in eig_val")
+            # if (eig_val <= 0).any():
+            #     print("Nonpositive eigenvalues:", eig_val)
+            # Compute the loss for the batch as the negative mean log of the eigenvalues
+            if torch.isinf(torch.log(eig_val)).any():
+                print("Log produced inf! Min/max eig_val:", eig_val.min().item(), eig_val.max().item())
+
+            batch_loss = -torch.mean(torch.log(torch.clamp(eig_val, min=1e-6, max=1e6)))
+
+            total_loss += batch_loss
+
+            # Update the start index for the next batch
+            start_idx = end_idx
+        # Compute the final loss by averaging across batches
+        loss = total_loss / len(batch_sizes)* diversity_weight
+        return loss, report
+    return diversity_loss_fn
+
+class Down_Model(nn.Module):
+    def __init__(self, in_dim, out_dim, hidden_dim=400, num_hidden_layers=2):
+        super(Down_Model, self).__init__()
+        
+        self.layers = nn.ModuleList([nn.Linear(in_dim, hidden_dim), nn.LeakyReLU()])
+        
+        for _ in range(num_hidden_layers - 1):
+            self.layers.append(nn.Linear(hidden_dim, hidden_dim))
+            self.layers.append(nn.LeakyReLU())
+        
+        self.layers.append(nn.Linear(hidden_dim, out_dim))
+
+    def forward(self, inputs):
+        x = inputs
+        for layer in self.layers:
+            x = layer(x)
+        return x
+
+class Up_Model(nn.Module):
+    def __init__(self, in_dim, out_dim, hidden_dim=400, num_hidden_layers=2):
+        super(Up_Model, self).__init__()
+        
+        self.layers = nn.ModuleList([nn.Linear(in_dim, hidden_dim), nn.LeakyReLU()])
+        
+        for _ in range(num_hidden_layers - 1):
+            self.layers.append(nn.Linear(hidden_dim, hidden_dim))
+            self.layers.append(nn.LeakyReLU())
+        
+        self.layers.append(nn.Linear(hidden_dim, out_dim))
+
+    def forward(self, inputs):
+        x = inputs
+        for layer in self.layers:
+            x = layer(x)
+        return x
+
+
+
+def GAN_step(D, G, D_opt, G_opt, data_batch, cond_batch, noise_batch, batch_size, device, auxiliary_loss_fn):
+    criterion = nn.BCEWithLogitsLoss()
+    D.zero_grad()
+    real_label = torch.full((batch_size,), 1, dtype=torch.float, device=device)
+    fake_label = torch.full((batch_size,), 0, dtype=torch.float, device=device)
+
+    data_and_condition = torch.cat([data_batch, cond_batch], dim=1)
+    noise_and_condition = torch.cat([noise_batch, cond_batch], dim=1)
+
+    output = D(data_and_condition).view(-1)
+    L_D_real = criterion(output, real_label)
+
+    fake_data = G(noise_and_condition)
+    fake_data_and_condition = torch.cat([fake_data, cond_batch], dim=1)
+    output = D(fake_data_and_condition.detach()).view(-1)
+    L_D_fake = criterion(output, fake_label)
+
+    L_D_tot = L_D_real + L_D_fake
+    L_D_tot.backward()
+    D_opt.step()
+
+    G.zero_grad()
+    fake_data = G(noise_and_condition)
+    fake_data_and_condition = torch.cat([fake_data, cond_batch], dim=1)
+    output = D(fake_data_and_condition).view(-1)
+    L_G = criterion(output, real_label)
+
+    if auxiliary_loss_fn is not None:
+        L_aux, rep= auxiliary_loss_fn(fake_data, cond_batch)
+        L_G_tot = L_G + L_aux
+
+        report = {"L_D_real": L_D_real.item(), "L_D_fake": L_D_fake.item(), "L_G": L_G.item(), "L_aux": L_aux.item()}
+        report.update(rep)
+    else:
+        L_G_tot = L_G
+        report = {"L_D_real": L_D_real.item(), "L_D_fake": L_D_fake.item(), "L_G": L_G.item()}
+
+    L_G_tot.backward()
+
+    torch.nn.utils.clip_grad_norm_(G.parameters(), max_norm=20)
+
+    G_opt.step()
+
+    return report
+
+
+def VAE_step(D, G, D_opt, G_opt, data_batch, cond_batch, noise_batch, batch_size, device, auxiliary_loss_fn):
+    
+    D.zero_grad()
+    G.zero_grad()
+    
+    alpha = 0.2
+
+    data_and_condition = torch.cat([data_batch, cond_batch], dim=1)
+
+    encoded = D(data_and_condition)
+    latent_dim = encoded.shape[1] // 2
+    mu = encoded[:, :latent_dim]
+    logvar = encoded[:, latent_dim:]
+    
+    std = torch.exp(0.5 * logvar)
+    eps = torch.randn_like(std)
+    z = mu + eps * std  # z = mu + sigma * epsilon
+    
+    # Forward pass through decoder (G)
+    noise_and_condition = torch.cat([z, cond_batch], dim=1)
+
+    reconstructed = G(noise_and_condition)
+    
+    # Compute losses
+    L_R = nn.MSELoss()(reconstructed, data_batch)
+    L_KL = -0.5 * torch.sum(1 + logvar - mu.pow(2) - logvar.exp()) / data_batch.size(0)
+
+    if auxiliary_loss_fn is not None:
+        L_aux, rep = auxiliary_loss_fn(reconstructed, cond_batch)
+
+        L_tot = alpha * L_KL + L_R + L_aux
+
+        report = {"L_KL": L_KL.item(), "L_R": L_R.item(), "L_tot": L_tot.item(), "L_aux": L_aux.item()}
+        report.update(rep)
+    else:
+        L_tot = alpha * L_KL + L_R
+
+        report = {"L_KL": L_KL.item(), "L_R": L_R.item(), "L_tot": L_tot.item()}
+
+    L_tot.backward()
+
+    # total_norm_D = torch.sqrt(sum(p.grad.norm()**2 for p in D.parameters() if p.grad is not None))
+    # total_norm_G = torch.sqrt(sum(p.grad.norm()**2 for p in G.parameters() if p.grad is not None))
+    # print(f"Gradient norm for D: {total_norm_D.item():.4f}")
+    # print(f"Gradient norm for G: {total_norm_G.item():.4f}")
+
+    torch.nn.utils.clip_grad_norm_(D.parameters(), max_norm=20)
+    torch.nn.utils.clip_grad_norm_(G.parameters(), max_norm=20)
+
+    D_opt.step()
+    G_opt.step()
+    
+    return report
+
+
+def DDPM_step_wrapper(scheduler: DDPMScheduler):
+    def DDPM_step(D, G, D_opt, G_opt, data_batch, cond_batch, noise_batch, batch_size, device, auxiliary_loss_fn):
+        # sample random t
+        t = torch.randint(0, scheduler.config.num_train_timesteps, (data_batch.size(0),), device=device).long()
+
+        # compute x_t using q_sample
+        noise = torch.randn_like(data_batch, device=device)
+        x_t = scheduler.add_noise(data_batch, noise, t)
+
+        # embed timestep and concat with cond
+        t_embedded = t.unsqueeze(-1).float() / scheduler.config.num_train_timesteps
+        x_input = torch.cat([x_t, t_embedded], dim=-1)
+
+        # predict noise
+        noise_pred = D(x_input)
+
+        # MSE loss with optional weighting (scheduler does not expose beta_t directly)
+        mse = MSELoss(reduction="none")(noise_pred, noise)
+        base_loss = mse.mean()
+
+        # reconstruct x0 and compute auxiliary loss
+        alpha_cumprod = scheduler.alphas_cumprod.to(device)
+        sqrt_alpha_cumprod_t = alpha_cumprod[t].sqrt().unsqueeze(-1)
+        sqrt_one_minus_alpha_cumprod_t = (1 - alpha_cumprod[t]).sqrt().unsqueeze(-1)
+        x0_pred = (x_t - sqrt_one_minus_alpha_cumprod_t * noise_pred) / sqrt_alpha_cumprod_t
+
+        if auxiliary_loss_fn is not None:
+            thresh = int(0.9 * scheduler.config.num_train_timesteps)
+            valid = (t < thresh)
+
+            total_loss = base_loss
+            report = {
+                "loss": base_loss.item(),
+            }
+        else:
+            total_loss = base_loss
+            report = {"loss": base_loss.item()}
+
+        D.zero_grad()
+        total_loss.backward()
+        D_opt.step()
+
+        return report
+
+    return DDPM_step
+
+
+def DDPM_step_cond_wrapper(scheduler: DDPMScheduler):
+    def DDPM_step_cond(D, G, D_opt, G_opt, data_batch, cond_batch, noise_batch, batch_size, device, auxiliary_loss_fn):
+        # sample random t
+        t = torch.randint(0, scheduler.config.num_train_timesteps, (data_batch.size(0),), device=device).long()
+
+        # compute x_t using q_sample
+        noise = torch.randn_like(data_batch, device=device)
+        x_t = scheduler.add_noise(data_batch, noise, t)
+
+        # embed timestep and concat with cond
+        t_embedded = t.unsqueeze(-1).float() / scheduler.config.num_train_timesteps
+        x_input = torch.cat([x_t, cond_batch, t_embedded], dim=-1)
+
+        # predict noise
+        noise_pred = D(x_input)
+
+        # MSE loss with optional weighting (scheduler does not expose beta_t directly)
+        mse = MSELoss(reduction="none")(noise_pred, noise)
+        base_loss = mse.mean()
+
+        # reconstruct x0 and compute auxiliary loss
+        alpha_cumprod = scheduler.alphas_cumprod.to(device)
+        sqrt_alpha_cumprod_t = alpha_cumprod[t].sqrt().unsqueeze(-1)
+        sqrt_one_minus_alpha_cumprod_t = (1 - alpha_cumprod[t]).sqrt().unsqueeze(-1)
+        x0_pred = (x_t - sqrt_one_minus_alpha_cumprod_t * noise_pred) / sqrt_alpha_cumprod_t
+
+        if auxiliary_loss_fn is not None:
+            thresh = int(0.5 * scheduler.config.num_train_timesteps)
+            valid = (t < thresh)
+
+            if valid.any():
+                x0_sub, cond_sub = x0_pred[valid], cond_batch[valid]
+                L_aux, rep = auxiliary_loss_fn(x0_sub, cond_sub)
+            else:
+                L_aux = torch.tensor(0.0, device=device)
+                rep = {}
+
+            total_loss = base_loss + L_aux
+            report = {
+                "loss": base_loss.item(),
+                "L_aux": L_aux.item(),
+                **rep
+            }
+        else:
+            total_loss = base_loss
+            report = {"loss": base_loss.item()}
+
+        D.zero_grad()
+        total_loss.backward()
+        D_opt.step()
+
+        return report
+
+    return DDPM_step_cond
+
+
+class ReusableDataLoader:
+    def __init__(self, dataset, batch_size, shuffle=True):
+        self.dataset = dataset
+        self.batch_size = batch_size
+        self.shuffle = shuffle
+        self.indices = list(range(len(self.dataset)))
+        self.previous_indices = []
+
+    def _shuffle_indices(self):
+        self.indices = torch.randperm(len(self.dataset)).tolist()
+
+    def get_batch(self):
+        queued = self.previous_indices
+        while len(queued) < self.batch_size:
+            if self.shuffle:
+                self._shuffle_indices()
+            queued.extend(self.indices)  # Add individual elements to queued list
+        
+        self.previous_indices = queued[self.batch_size:]  # Store remaining indices for the next batch
+        batch_indices = queued[:self.batch_size]  # Get the batch of the correct size
+        return torch.stack([self.dataset[i][0] for i in batch_indices])
+
+def train(D, G, D_opt, G_opt, loader, num_steps, batch_size, noise_dim, train_step_fn, device, auxiliary_loss_fn, condition_sampler):
+    # Loss function
+    D.train()
+    G.train()
+    steps_range = trange(num_steps, position=0, leave=True)
+
+    for step in steps_range:
+        data_batch = loader.get_batch().to(device)
+        noise_batch = torch.randn(batch_size, noise_dim).to(device)
+        cond_batch = condition_sampler(batch_size).to(device)
+        effective_auxiliary_loss_fn = auxiliary_loss_fn
+        report = train_step_fn(D, G, D_opt, G_opt, data_batch, cond_batch, noise_batch, batch_size, device, effective_auxiliary_loss_fn)
+        postfix = {key: "{:.4f}".format(value) for key, value in report.items()}
+        steps_range.set_postfix(postfix)
+    return D, G
+
+
+def get_DDPM_generate_cond(scheduler: DDPMScheduler, data_dim, batch_size=64):
+    def DDPM_generate_cond(D, G, cond_batch, latent_dim, device, batch_size=batch_size):
+        with torch.no_grad():
+            results = []
+            numgen = cond_batch.shape[0]
+            for start_idx in range(0, numgen, batch_size):
+                end_idx = min(start_idx + batch_size, numgen)
+                current_batch_size = end_idx - start_idx
+
+                x = torch.randn(current_batch_size, data_dim).to(device)
+
+                for t in reversed(range(scheduler.config.num_train_timesteps)):
+                    t_tensor = torch.full((current_batch_size,), t, dtype=torch.long, device=device)
+                    t_embedded = t_tensor.unsqueeze(-1).float() / scheduler.config.num_train_timesteps
+                    x_input = torch.cat([x, cond_batch, t_embedded], dim=-1)
+
+                    with torch.no_grad():
+                        noise_pred = D(x_input)
+
+                    x = scheduler.step(model_output=noise_pred, timestep=t, sample=x).prev_sample
+                results.append(x)
+
+            return torch.cat(results, dim=0)
+    return DDPM_generate_cond
+
+def get_DDPM_generate_guided(scheduler: DDPMScheduler, data_dim, auxiliary_loss_fn, batch_size=64):
+    def DDPM_generate_guided(D, G, cond_batch, latent_dim, device, auxiliary_loss_fn=auxiliary_loss_fn, batch_size=batch_size):
+        results = []
+        numgen = cond_batch.shape[0]
+
+        num_guided_timesteps = int(0.5 * scheduler.config.num_train_timesteps)
+
+        for start_idx in range(0, numgen, batch_size):
+            end_idx = min(start_idx + batch_size, numgen)
+            current_batch_size = end_idx - start_idx
+
+            # Start with pure noise
+            x = torch.randn(current_batch_size, data_dim, device=device)
+
+            for t in tqdm(reversed(range(scheduler.config.num_train_timesteps))):
+                t_tensor = torch.full((current_batch_size,), t, dtype=torch.long, device=device)
+                t_embedded = t_tensor.unsqueeze(-1).float() / scheduler.config.num_train_timesteps
+                x_input = torch.cat([x, t_embedded], dim=-1)
+
+                # Model prediction
+                noise_pred = D(x_input)
+
+                # Apply auxiliary loss only if t < threshold (later timesteps)
+                threshold = int(0.5 * scheduler.config.num_train_timesteps)
+                if t < threshold:
+                    alpha_cumprod = scheduler.alphas_cumprod.to(device)
+                    sqrt_alpha_cumprod_t = alpha_cumprod[t].sqrt().unsqueeze(-1)
+                    sqrt_one_minus_alpha_cumprod_t = (1 - alpha_cumprod[t]).sqrt().unsqueeze(-1)
+                    x0_pred = (x - sqrt_one_minus_alpha_cumprod_t * noise_pred) / sqrt_alpha_cumprod_t
+
+                    aux_loss, _ = auxiliary_loss_fn(x0_pred, cond_batch)
+                    aux_loss.backward(retain_graph=True)
+                    grad = x.grad / num_guided_timesteps
+
+                    # Update x based on this local gradient only
+                    x = (x - grad).detach().requires_grad_(True)
+                    x.retain_grad()
+
+                # Always apply scheduler step
+                x = scheduler.step(model_output=noise_pred, timestep=t, sample=x).prev_sample
+                x.retain_grad()
+
+            results.append(x)
+
+        return torch.cat(results, dim=0)
+
+    return DDPM_generate_guided
+
+
+
+
+
+
+def VAE_generate(D, G, cond_batch, latent_dim, device):
+    with torch.no_grad():
+        numgen = cond_batch.shape[0]
+        z = torch.randn(numgen, latent_dim).to(device)
+        z_and_condition = torch.cat([z, cond_batch], dim=1)
+        generated_data = G(z_and_condition)
+        return generated_data
+
+# def VAE_generate_cond(D, G, cond_batch, latent_dim, device):
+#     numgen = cond_batch.shape[0]
+#     z = torch.randn(numgen, latent_dim).to(device)
+#     labels = torch.ones(numgen, 1).to(device)
+#     z = torch.cat([z, labels], dim=1)
+#     generated_data = G(z)
+#     return generated_data
+    
+def GAN_generate(D, G, cond_batch, noise_dim, device):
+    with torch.no_grad():
+        numgen = cond_batch.shape[0]
+        noise = torch.randn(numgen, noise_dim).to(device)
+        noise_and_condition = torch.cat([noise, cond_batch], dim=1)
+        labels = torch.ones(numgen, 1).to(device)
+        generated_data = G(noise_and_condition)
+        return generated_data
+
+def train_model(data, model_type, train_params, auxiliary_loss_fn, condition_sampler, device):
+    batch_size, disc_lr, gen_lr, noise_dim, num_epochs, n_hidden, layer_size= train_params
+
+    
+
+    data = torch.tensor(data).float()
+    sample_condition = sample_continuous(data.shape[0], randomize=True).to(device)
+
+    loader = ReusableDataLoader(TensorDataset(data), batch_size)
+
+    data_dim = data.shape[1]
+    cond_dim = sample_condition.shape[1]
+
+    if model_type in ["GAN"]:
+        train_step = GAN_step
+        generate_fn = GAN_generate
+        D_in = data_dim + cond_dim
+        D_out = 1
+        G_in = noise_dim +cond_dim
+        G_out = data_dim
+    elif model_type in ["VAE"]:
+        train_step = VAE_step
+        generate_fn = VAE_generate
+        D_in = data_dim + cond_dim
+        D_out = 2*noise_dim
+        G_in = noise_dim + cond_dim
+        G_out = data_dim
+    elif model_type in ["DDPM_guided"]:
+        scheduler = DDPMScheduler(num_train_timesteps=100)
+        train_step = DDPM_step_wrapper(scheduler)
+        generate_fn = get_DDPM_generate_guided(scheduler, data_dim, auxiliary_loss_fn, batch_size=batch_size)
+        D_in = data_dim + 1
+        D_out = data_dim
+        G_in = 1 #unused
+        G_out = 1 #unused
+    elif model_type in ["DDPM_conditional"]:
+        scheduler = DDPMScheduler(num_train_timesteps=1000)
+        train_step = DDPM_step_cond_wrapper(scheduler)
+        generate_fn = get_DDPM_generate_cond(scheduler, data_dim, batch_size=batch_size)
+        D_in = data_dim + cond_dim + 1
+        D_out = data_dim
+        G_in = 1 #unused
+        G_out = 1 #unused
+    # else:
+    #     raise ValueError("Invalid model_type")
+
+
+    D = Down_Model(D_in, D_out, layer_size, n_hidden)
+    G = Up_Model(G_in, G_out, layer_size, n_hidden)
+
+    
+    D.to(device)
+    G.to(device)
+    D_opt = torch.optim.Adam(D.parameters(), lr=disc_lr, betas=(0.5,0.999))
+    G_opt = torch.optim.Adam(G.parameters(), lr=gen_lr, betas=(0.5,0.999))
+
+    
+
+    if num_epochs>0:
+        num_steps = num_epochs*len(data)//batch_size
+    else:
+        num_steps = -num_epochs #hacky way to specify fixed number of steps rather than epochs
+    
+
+    train(D, G, D_opt, G_opt, loader, num_steps, batch_size, noise_dim, train_step, device, auxiliary_loss_fn, condition_sampler)
+
+    return D, G, generate_fn

bike_bench_internal/benchmark_models/libmoon/example.py

@@ -0,0 +1,4 @@
+def add_one(number):
+    return number + 1
+
+

bike_bench_internal/benchmark_models/libmoon/problem/mop.py

@@ -0,0 +1,93 @@
+import numpy as np
+import torch
+
+class mop():
+    def __init__(self,
+                 n_var: int,
+                 n_obj: int,
+                 lbound: np.ndarray,
+                 ubound: np.ndarray,
+                 n_cons: int=0,
+                 ) -> None:
+
+        self.n_var = n_var
+        self.n_obj = n_obj
+        self.n_cons = n_cons
+
+        self.lbound=lbound
+        self.ubound=ubound
+
+
+    @property
+    def get_number_variable(self) -> int:
+        return self.n_var
+
+    @property
+    def get_number_objective(self) -> int:
+        return self.n_obj
+
+    @property
+    def get_lower_bound(self) -> np.ndarray:
+        return self.lbound
+
+    @property
+    def get_upper_bound(self) -> np.ndarray:
+        return self.ubound
+
+    @property
+    def has_constraint(self) -> bool:
+        return self.n_cons > 0
+
+    def evaluate(self, x):
+        raise NotImplementedError("Subclasses should implement this method.")
+
+    def __call__(self, x):
+        return self.evaluate(x)
+
+    def evaluate(self, x: any) -> any:
+        """
+            Evaluate the objectives for x
+            Parameters
+            ----------
+            x : any
+                Tensor or ndarray
+            Returns
+            -------
+            any
+                Tensor or ndarray correspondingly
+            Raises
+            ------
+            ValueError
+                wrong type of x
+        """
+
+        if type(x) == torch.Tensor:
+            return self._evaluate_torch(torch.atleast_2d(x))
+        elif isinstance(x, np.ndarray):
+            return self._evaluate_numpy(np.atleast_2d(x))
+        else:
+            raise ValueError("Input has to be in the form of Tensor or ndarray!")
+
+    def get_pf(self, n_points: int=100) -> np.ndarray:
+        """
+        Get Pareto front
+        Parameters
+        ----------
+        num_points : int, optional
+            _description_, by default 100
+        Returns
+        -------
+        np.ndarray
+            _description_
+        """
+        # TODO
+        # if method=='uniform':
+        if hasattr(self, "_get_pf"): return self._get_pf(n_points)
+        else: raise NotImplementedError("Subclasses should implement this method.")
+
+
+
+class mop_noCons(mop):
+
+    def __init__(self, n_var: int, n_obj: int, lbound: np.ndarray, ubound: np.ndarray, n_cons: int = 0) -> None:
+        super().__init__(n_var, n_obj, lbound, ubound, n_cons)

bike_bench_internal/benchmark_models/libmoon/problem/mtl/loaders/__init__.py

@@ -0,0 +1 @@
+from .multimnist_loader import MultiMNISTData

bike_bench_internal/benchmark_models/libmoon/problem/mtl/loaders/adult_loader.py

@@ -0,0 +1,124 @@
+import os
+import pandas as pd
+import torch
+
+from sklearn.model_selection import train_test_split
+from sklearn.preprocessing import StandardScaler
+from torch.utils import data
+from libmoon.util_global.constant import root_name
+
+
+
+
+def load_dataset(path, s_label):
+    # print()
+    data = pd.read_csv(path)
+    # Preprocessing taken from https://www.kaggle.com/islomjon/income-prediction-with-ensembles-of-decision-trees
+    # replace missing values with majority class
+    data['workclass'] = data['workclass'].replace('?','Private')
+    data['occupation'] = data['occupation'].replace('?','Prof-specialty')
+    data['native-country'] = data['native-country'].replace('?','United-States')
+
+    # education category
+    data.education = data.education.replace(['Preschool','1st-4th','5th-6th','7th-8th','9th','10th','11th','12th'],'left')
+    data.education = data.education.replace('HS-grad','school')
+    data.education = data.education.replace(['Assoc-voc','Assoc-acdm','Prof-school','Some-college'],'higher')
+    data.education = data.education.replace('Bachelors','undergrad')
+    data.education = data.education.replace('Masters','grad')
+    data.education = data.education.replace('Doctorate','doc')
+
+    # marital status
+    data['marital-status'] = data['marital-status'].replace(['Married-civ-spouse','Married-AF-spouse'],'married')
+    data['marital-status'] = data['marital-status'].replace(['Never-married','Divorced','Separated','Widowed', 'Married-spouse-absent'], 'not-married')
+
+    # income
+    data.income = data.income.replace('<=50K', 0)
+    data.income = data.income.replace('>50K', 1)
+
+    # sex
+    data.gender = data.gender.replace('Male', 0)
+    data.gender = data.gender.replace('Female', 1)
+    
+    # mtldata.race = mtldata.race.replace('White', 0)
+    # mtldata.race = mtldata.race.replace('Black', 1)
+    # mtldata.race = mtldata.race.astype(int)
+
+    # encode categorical values
+    data1 = data.copy()
+    data1 = pd.get_dummies(data1)
+    data1 = data1.drop(['income', s_label], axis=1)
+    # data1 = data1.drop(['income', s_label], axis=1)
+
+    X = StandardScaler().fit(data1).transform(data1)
+    y = data['income'].values
+    s1 = data[s_label].values
+    # s2 = mtldata['race'].values
+
+    return X, y, s1
+
+
+
+
+class ADULT(data.Dataset):
+
+
+    def __init__(self, split="train", sensible_attribute="gender"):
+        assert split in ["train", "val", "test"]
+
+        # folder_name = os.path.dirname( os.path.dirname(__file__) )
+
+
+
+
+        path = os.path.join(root_name, 'mtldata', "adult.csv")
+
+        x, y, s1 = load_dataset(path, sensible_attribute)
+
+
+        x = torch.from_numpy(x).float()
+        y = torch.from_numpy(y).long()
+        s1 = torch.from_numpy(s1).long()
+        # s2 = torch.from_numpy(s2).long()
+
+        # train/val/test split: 70/10/20 %
+        x_train, x_test, y_train, y_test, s1_train, s1_test  = train_test_split(x, y, s1,test_size=.2, random_state=1)
+        x_train, x_val, y_train, y_val, s1_train, s1_val= train_test_split(x_train, y_train, s1_train, test_size=.125, random_state=1)
+
+        if split == 'train':
+            self.x = x_train
+            self.y = y_train
+            self.s1 = s1_train
+            # self.s2 = s2_train
+            
+        elif split == 'val':
+            self.x = x_val
+            self.y = y_val
+            self.s1 = s1_val
+            # self.s2 = s2_val
+        elif split == 'test':
+            self.x = x_test
+            self.y = y_test
+            self.s1 = s1_test
+            # self.s2 = s2_test
+        
+        print("loaded {} instances for split {}. y positives={}, {} positives={}".format(
+            len(self.y), split, sum(self.y), sensible_attribute, sum(self.s1)))
+
+    def __len__(self):
+        """__len__"""
+        return len(self.x)
+
+    def __getitem__(self, index):
+        return dict(data=self.x[index], labels=self.y[index], sensible_attribute=self.s1[index])
+
+    def task_names(self):
+        return None
+
+
+    
+if __name__ == "__main__":
+    dataset = ADULT(split="train")
+    trainloader = data.DataLoader(dataset, batch_size=256, num_workers=0)
+
+    for i, data in enumerate(trainloader):
+        break

bike_bench_internal/benchmark_models/libmoon/problem/mtl/loaders/compas_loader.py

@@ -0,0 +1,101 @@
+import torch
+import os
+import pandas as pd
+
+from datetime import datetime
+from sklearn.model_selection import train_test_split
+from sklearn.preprocessing import StandardScaler
+
+
+def load_dataset(path, s_label):
+
+    # following the preprocessing on 
+    # https://github.com/propublica/compas-analysis/blob/master/Compas%20Analysis.ipynb
+    raw_data = pd.read_csv(path)
+
+    data = raw_data[(
+            (raw_data['days_b_screening_arrest'] <= 30) &
+            (raw_data['days_b_screening_arrest'] >= -30) &
+            (raw_data['is_recid'] != -1) &
+            (raw_data['c_charge_degree'] != 'O') &
+            (raw_data['score_text'] != 'N/A')
+    )]
+
+    # Only some columns are relevant
+    data = data[['age', 'c_charge_degree', 'race', 'age_cat', 'score_text', 'sex', 'priors_count',
+                 'days_b_screening_arrest', 'decile_score', 'is_recid', 'two_year_recid', 'c_jail_in', 'c_jail_out']]
+
+    # convert c_jail_in and c_jail_out to time_in_jail, mesured in hours
+    def date_from_str(s):
+        return datetime.strptime(s, '%Y-%m-%d %H:%M:%S')
+
+    data['c_jail_in'] = data['c_jail_in'].apply(date_from_str)
+    data['c_jail_out'] = data['c_jail_out'].apply(date_from_str)
+
+    data['length_of_stay'] = data['c_jail_out'] - data['c_jail_in']
+
+    # data['length_of_stay'] = data['length_of_stay'].astype('timedelta64[h]')  # modified by xz, 11.6
+    data['length_of_stay'] = data['length_of_stay'].dt.days
+
+
+    data = data.drop(['c_jail_in', 'c_jail_out'], axis=1)
+
+    # encode sex
+    data['sex'] = data['sex'].replace('Male', 0)
+    data['sex'] = data['sex'].replace('Female', 1)
+
+    # one-hot encode categorical variables
+    data1 = data.copy()
+    data1 = data1.drop(['two_year_recid', 'sex'], axis=1)
+    data1 = pd.get_dummies(data1)
+
+    x = StandardScaler().fit(data1).transform(data1)
+    y = data['two_year_recid'].values
+    s = data['sex'].values
+
+    return x, y, s
+
+
+class Compas(torch.utils.data.Dataset):
+
+    def __init__(self, split, sensible_attribute='sex'):
+        assert split in ['train', 'val', 'test']
+
+        folder_name = os.path.dirname(os.path.dirname(__file__))
+        path = os.path.join(folder_name, 'mtldata', "compas.csv")
+
+        x, y, s = load_dataset(path, sensible_attribute)
+
+        x = torch.from_numpy(x).float()
+        y = torch.from_numpy(y).long()
+        s = torch.from_numpy(s).long()
+
+        # train/val/test split: 70/10/20 %
+        x_train, x_test, y_train, y_test, s_train, s_test = train_test_split(x, y, s, test_size=.2, random_state=1)
+        x_train, x_val, y_train, y_val, s_train, s_val = train_test_split(x_train, y_train, s_train, test_size=.125, random_state=1)
+
+        if split == 'train':
+            self.x = x_train
+            self.y = y_train
+            self.s = s_train
+        elif split == 'val':
+            self.x = x_val
+            self.y = y_val
+            self.s = s_val
+        elif split == 'test':
+            self.x = x_test
+            self.y = y_test
+            self.s = s_test
+        
+        print("loaded {} instances for split {}. y positives={}, {} positives={}".format(
+            len(self.y), split, sum(self.y), sensible_attribute, sum(self.s)))
+        
+    
+    def __len__(self):
+        return len(self.y)
+    
+    def __getitem__(self, index):
+        return dict(data=self.x[index], labels=self.y[index], sensible_attribute=self.s[index])
+    
+    def task_names(self):
+        return None

bike_bench_internal/benchmark_models/libmoon/problem/mtl/loaders/credit_loader.py

@@ -0,0 +1,80 @@
+import torch
+import os
+import pandas as pd
+
+from sklearn.model_selection import train_test_split
+from sklearn.preprocessing import StandardScaler
+
+
+
+def load_dataset(path, s_label):
+    data = pd.read_csv( path )
+
+    # convert categorical columns
+    to_categorical = [
+        'EDUCATION', 'MARRIAGE', 'PAY_0', 'PAY_2', 
+        'PAY_3', 'PAY_4', 'PAY_5', 'PAY_6'
+    ]
+    for column in to_categorical:
+        data[column] = data[column].astype('category')
+
+    data.SEX = data.SEX.replace(1, 0)   # male
+    data.SEX = data.SEX.replace(2, 1)   # female
+
+    # Scale and split
+    data1 = data.copy()
+    data1 = data1.drop(['default.payment.next.month', s_label], axis=1)
+    data1 = pd.get_dummies(data1)
+
+    x = StandardScaler().fit(data1).transform(data1)
+    y = data['default.payment.next.month'].values
+    s = data[s_label].values
+
+    return x, y, s
+
+
+class Credit(torch.utils.data.Dataset):
+
+    def __init__(self, split, sensible_attribute='SEX'):
+        assert split in ['train', 'val', 'test']
+
+        folder_name = os.path.dirname(os.path.dirname(__file__))
+        path = os.path.join(folder_name, 'mtldata', "credit.csv")
+
+        x, y, s = load_dataset(path, sensible_attribute)
+
+        x = torch.from_numpy(x).float()
+        y = torch.from_numpy(y).long()
+        s = torch.from_numpy(s).long()
+
+        # train/val/test split: 70/10/20 %
+        x_train, x_test, y_train, y_test, s_train, s_test = train_test_split(x, y, s, test_size=.2, random_state=1)
+        x_train, x_val, y_train, y_val, s_train, s_val = train_test_split(x_train, y_train, s_train, test_size=.125, random_state=1)
+
+        if split == 'train':
+            self.x = x_train
+            self.y = y_train
+            self.s = s_train
+        elif split == 'val':
+            self.x = x_val
+            self.y = y_val
+            self.s = s_val
+        elif split == 'test':
+            self.x = x_test
+            self.y = y_test
+            self.s = s_test
+        
+        print("loaded {} instances for split {}. y positives={}, {} positives={}".format(
+            len(self.y), split, sum(self.y), sensible_attribute, sum(self.s)))
+
+    
+    def __len__(self):
+        return len(self.y)
+    
+    def __getitem__(self, index):
+        return dict(data=self.x[index], labels=self.y[index], sensible_attribute=self.s[index])
+    
+    def task_names(self):
+        return None
+
+

bike_bench_internal/benchmark_models/libmoon/problem/mtl/loaders/multimnist_loader.py

@@ -0,0 +1,101 @@
+import torch
+import pickle
+from sklearn.model_selection import train_test_split
+# taken from https://github.com/Xi-L/ParetoMTL and adapted
+import os
+
+from libmoon.util_global.constant import root_name
+
+
+class MultiMNISTData(torch.utils.data.Dataset):
+    """
+    The datasets from ParetoMTL
+    """
+    def __init__(self, dataset, split, root='data/multi', **kwargs):
+        assert dataset in ['mnist', 'fashion', 'fmnist']
+        assert split in ['train', 'val', 'test']
+
+        # equal size of val and test split
+        train_split = .9
+
+        if dataset == 'mnist':
+            self.path = os.path.join(root_name, 'problem', 'mtl', 'data', 'multimnist', 'mnist.pickle')
+
+        elif dataset == 'fashion':
+            self.path = os.path.join(root_name, 'problem', 'mtl', 'data', 'multimnist', 'fashion.pickle')
+
+        elif dataset == 'fmnist':
+            self.path = os.path.join(root_name, 'problem', 'mtl', 'data', 'multimnist', 'fmnist.pickle')
+
+
+        self.val_size = .1
+        with open(self.path, 'rb') as f:
+            trainX, trainLabel, testX, testLabel = pickle.load(f)
+
+        n_train = len(trainX)
+        if self.val_size > 0:
+            trainX, valX, trainLabel, valLabel = train_test_split(
+                trainX, trainLabel, test_size=self.val_size, random_state=42
+            )
+            n_train = len(trainX)
+            n_val = len(valX)
+
+        trainX = torch.from_numpy(trainX.reshape(n_train, 1, 36, 36)).float()
+        trainLabel = torch.from_numpy(trainLabel).long()
+        testX = torch.from_numpy(testX.reshape(20000, 1, 36, 36)).float()
+        testLabel = torch.from_numpy(testLabel).long()
+
+        if self.val_size > 0:
+            valX = torch.from_numpy(valX.reshape(n_val, 1, 36, 36)).float()
+            valLabel = torch.from_numpy(valLabel).long()
+
+        if split in ['train', 'val']:
+            n = int(len(trainX) * train_split)
+            if split == 'val':
+                self.X = valX
+                self.y = valLabel
+            elif split == 'train':
+                self.X = trainX
+                self.y = trainLabel
+        elif split == 'test':
+            self.X = testX
+            self.y = testLabel
+
+    def __getitem__(self, index):
+        return dict(data=self.X[index], labels_l=self.y[index, 0], labels_r=self.y[index, 1])
+
+    def __len__(self):
+        return len(self.X)
+
+    def task_names(self):
+        return ['l', 'r']
+
+
+
+
+
+
+
+if __name__ == '__main__':
+    import matplotlib.pyplot as plt
+    dst = MultiMNISTData(dataset='mnist', split='val')
+    loader = torch.utils.data.DataLoader(dst, batch_size=10, shuffle=True, num_workers=0)
+
+    for dat in loader:
+        ims = dat['data'].view(10, 36, 36).numpy()
+        labs_l = dat['labels_l']
+        labs_r = dat['labels_r']
+
+        f, axarr = plt.subplots(2, 5)
+        for j in range(5):
+            for i in range(2):
+                axarr[i][j].imshow(ims[j * 2 + i, :, :], cmap='gray')
+                axarr[i][j].set_title('{}_{}'.format(labs_l[j * 2 + i], labs_r[j * 2 + i]))
+        plt.show()
+        a = input()
+
+        if a == 'ex':
+            break
+        else:
+            plt.close()
+

bike_bench_internal/benchmark_models/libmoon/problem/mtl/mnist.py

@@ -0,0 +1,189 @@
+import matplotlib.pyplot as plt
+from libmoon.util_global.constant import root_name
+from libmoon.problem.mtl.loaders.multimnist_loader import MultiMNISTData
+import torch
+
+from libmoon.problem.mtl.objectives import CrossEntropyLoss
+from libmoon.problem.mtl.model.simple import MultiLeNet
+from libmoon.util_global.weight_factor.funs import uniform_pref
+from libmoon.util_global.constant import FONT_SIZE
+
+loss_1 = CrossEntropyLoss(label_name='labels_l', logits_name='logits_l')
+loss_2 = CrossEntropyLoss(label_name='labels_r', logits_name='logits_r')
+
+from tqdm import tqdm
+import numpy as np
+from numpy import array
+import os
+from solver.gradient import get_core_solver
+from solver.gradient.utils.util import get_grads_from_model, numel_params
+from libmoon.util_global.constant import is_pref_based
+import itertools
+
+class MultiMnistProblem:
+
+    # How to train at the same time.
+    def __init__(self, args, prefs):
+        self.dataset = MultiMNISTData('mnist', 'train')
+        self.args = args
+        self.loader = torch.utils.data.DataLoader(self.dataset, batch_size=self.args.batch_size, shuffle=True,
+                                                  num_workers=0)
+        self.dataset_test = MultiMNISTData('mnist', 'test')
+        self.loader_test = torch.utils.data.DataLoader(self.dataset_test, batch_size=args.batch_size, shuffle=True,
+                                                       num_workers=0)
+
+        self.lr = args.lr
+        self.prefs = prefs
+        self.n_prob = len(prefs)
+
+        self.model_arr = [MultiLeNet([1, 36, 36]) for _ in range(self.n_prob)]
+
+        num_params = numel_params(self.model_arr[0])
+        print('num_params: ', num_params)
+
+        for model in self.model_arr:
+            model.to(args.device)
+
+        self.is_pref_flag = is_pref_based(args.mtd)
+
+        if self.is_pref_flag:
+            self.core_solver_arr = [get_core_solver(args, pref) for pref in prefs]
+            self.optimizer_arr = [torch.optim.Adam(self.model_arr[idx].parameters(), lr=self.lr) for idx in
+                                  range(self.n_prob)]
+        else:
+            self.set_core_solver = get_core_solver(args)
+
+            params = [model.parameters() for model in self.model_arr]
+            self.set_optimizer = torch.optim.Adam(itertools.chain(*params), lr=0.01)
+
+
+    def optimize(self):
+        loss_all = []
+        for _ in tqdm(range(self.args.num_epoch)):
+            if self.is_pref_flag:
+                loss_hostory = [[] for i in range(self.n_prob)]
+            else:
+                loss_hostory = []
+
+            for data in self.loader:
+                data_ = {k: v.to(self.args.device) for k, v in data.items()}
+
+                # pref based mtd
+                if self.is_pref_flag:
+                    for pref_idx, (pref, model, optimizer) in enumerate(
+                            zip(self.prefs, self.model_arr, self.optimizer_arr)):
+                        logits_dict = self.model_arr[pref_idx](data_)
+                        logits_dict['labels_l'] = data_['labels_l']
+                        logits_dict['labels_r'] = data_['labels_r']
+                        l1 = loss_1(**logits_dict)
+                        l2 = loss_2(**logits_dict)
+
+                        l_contains_grad = [l1, l2]
+                        G = get_grads_from_model(l_contains_grad, model)
+
+                        l1_np = np.array(l1.cpu().detach().numpy(), copy=True)
+                        l2_np = np.array(l2.cpu().detach().numpy(), copy=True)
+                        losses = array([l1_np, l2_np])
+                        alpha = self.core_solver_arr[pref_idx].get_alpha(G = G, losses=losses)
+                        self.optimizer_arr[pref_idx].zero_grad()
+                        (alpha[0] * l1 + alpha[1] * l2).backward()
+                        self.optimizer_arr[pref_idx].step()
+                        l1_np = np.array(l1.cpu().detach().numpy(), copy=True)
+                        l2_np = np.array(l2.cpu().detach().numpy(), copy=True)
+                        loss_hostory[pref_idx].append([l1_np, l2_np])
+                else:
+                    # set based method is more complicated.
+                    losses = [0,] * self.n_prob
+                    losses_ts = [0] * self.n_prob
+
+                    for model_idx, model in enumerate(self.model_arr):
+                        logits_dict = self.model_arr[model_idx](data_)
+                        logits_dict['labels_l'] = data_['labels_l']
+                        logits_dict['labels_r'] = data_['labels_r']
+                        l1 = loss_1(**logits_dict)
+                        l2 = loss_2(**logits_dict)
+
+                        losses_ts[model_idx] = torch.stack([l1, l2])
+
+                        l1_np, l2_np = np.array(l1.cpu().detach().numpy(), copy=True), np.array(l2.cpu().detach().numpy(), copy=True)
+                        losses[model_idx] = [l1_np, l2_np]
+
+                    losses_ts = torch.stack(losses_ts)
+                    losses = np.array(losses)
+                    alpha = self.set_core_solver.get_alpha(losses).to(self.args.device)
+                    self.set_optimizer.zero_grad()
+                    torch.sum(alpha * losses_ts).backward()
+                    self.set_optimizer.step()
+                    loss_hostory.append(losses)
+            loss_hostory = np.array(loss_hostory)
+            if args.is_pref_based:
+                loss_history_mean = np.mean(loss_hostory, axis=1)
+            else:
+                loss_history_mean = np.mean(loss_hostory, axis=0)
+            loss_all.append(loss_history_mean)
+        return loss_all
+
+
+
+if __name__ == '__main__':
+    import argparse
+    parser = argparse.ArgumentParser()
+
+    parser.add_argument('--problem', default='mnist', type=str)  # For attribute in args, we all call problem.
+    parser.add_argument('--split', default='train', type=str)
+    parser.add_argument('--batch_size', default=512, type=int)
+    parser.add_argument('--shuffle', default=True, type=bool)
+    parser.add_argument('--lr', default=1e-2, type=float)
+    parser.add_argument('--num_epoch', default=10, type=int)
+    parser.add_argument('--use-cuda', default=True, type=bool)
+
+    parser.add_argument('--mtd', default='hvgrad', type=str)
+    parser.add_argument('--agg-mtd', default='ls', type=str)   # This att is only valid when args.mtd=agg.
+    parser.add_argument('--n-obj', default=2, type=int)   # This att is only valid when args.mtd=agg.
+
+    args = parser.parse_args()
+    args.is_pref_based = is_pref_based(args.mtd)
+    if torch.cuda.is_available() and args.use_cuda:
+        args.device = torch.device("cuda")  # Use the GPU
+        print('cuda is available')
+    else:
+        args.device = torch.device("cpu")  # Use the CPU
+        print('cuda is not available')
+
+    prefs = uniform_pref(n_partition=10, n_obj=2, clip_eps=0.1)
+    args.n_prob = len(prefs)
+
+    problem = MultiMnistProblem(args, prefs)
+    # args.n_obj = problem.n_obj
+
+    loss_history = problem.optimize()
+    loss_history = np.array(loss_history)
+
+    final_solution = loss_history[-1,:,:]
+    # plt.scatter(final_solution[:,0], final_solution[:,1], label='final solution')
+    for idx in range(loss_history.shape[1]):
+        plt.plot(loss_history[:,idx,0], loss_history[:,idx,1], 'o-', label='pref {}'.format(idx))
+
+    plt.plot(final_solution[:,0], final_solution[:,1], color='k', linewidth=3)
+
+    plt.legend(fontsize=FONT_SIZE)
+    # draw pref
+    solution_norm = np.linalg.norm(final_solution, axis=1, keepdims=True)
+    prefs_norm = prefs / np.linalg.norm(prefs, axis=1, keepdims=True) * solution_norm
+
+    if args.is_pref_based:
+        for pref in prefs_norm:
+            plt.plot([0, pref[0]], [0, pref[1]], color='k')
+
+
+    plt.xlabel('$L_1$', fontsize=FONT_SIZE)
+    plt.ylabel('$L_2$', fontsize=FONT_SIZE)
+
+
+    folder_name = os.path.join( root_name, 'output', args.problem, args.mtd)
+    os.makedirs(folder_name, exist_ok=True)
+    fig_name = os.path.join(folder_name, 'final_solution.svg')
+    plt.savefig(fig_name)
+    print('saved in ', fig_name)
+
+    plt.show()

bike_bench_internal/benchmark_models/libmoon/problem/mtl/model/simple.py

@@ -0,0 +1,77 @@
+import torch.nn as nn
+
+'''
+    MTL problems need . 
+'''
+
+
+
+class MultiLeNet(nn.Module):
+
+    def __init__(self, dim, **kwargs):
+
+        '''
+            :param dim: a 3d-array. [chanel, height, width]
+            :param kwargs:
+        '''
+        super().__init__()
+        self.shared = nn.Sequential(
+            nn.Conv2d(dim[0], 10, kernel_size=5),
+            nn.MaxPool2d(kernel_size=2),
+            nn.ReLU(),
+            nn.Conv2d(10, 20, kernel_size=5),
+            nn.MaxPool2d(kernel_size=2),
+            nn.ReLU(),
+            nn.Flatten(),
+            nn.Linear(720, 50),
+            nn.ReLU(),
+        )
+        self.private_left = nn.Linear(50, 10)
+        self.private_right = nn.Linear(50, 10)
+
+
+    def forward(self, batch):
+        x = batch['data']
+        x = self.shared(x)
+        return dict(logits_l=self.private_left(x), logits_r=self.private_right(x))
+
+    def private_params(self):
+        return ['private_left.weight', 'private_left.bias', 'private_right.weight', 'private_right.bias']
+
+
+
+
+
+class FullyConnected(nn.Module):
+    def __init__(self, dim, **kwargs):
+        super().__init__()
+        self.f = nn.Sequential(
+            nn.Linear(dim[0], 60),
+            nn.ReLU(),
+            nn.Linear(60, 25),
+            nn.ReLU(),
+            nn.Linear(25, 1),
+        )
+
+    def forward(self, batch):
+        x = batch['data']
+        return dict(logits=self.f(x))
+
+
+
+
+
+if __name__ == '__main__':
+    from libmoon.util_global.constant import root_name
+    import os
+    import pickle
+
+
+
+
+    pickle_name = os.path.join(root_name, 'problem', 'mtl', 'data', 'multimnist', 'mnist.pickle')
+    with open(pickle_name, 'rb') as f:
+        data = pickle.load(f)
+
+    model = MultiLeNet([3, 32, 32])
+    print('hello world')

bike_bench_internal/benchmark_models/libmoon/problem/mtl/objectives.py

@@ -0,0 +1,179 @@
+import torch
+
+def from_name(names, task_names):
+    objectives = {
+        'CrossEntropyLoss': CrossEntropyLoss,
+        'BinaryCrossEntropyLoss': BinaryCrossEntropyLoss,
+        'L1Regularization': L1Regularization,
+        'L2Regularization': L2Regularization,
+        'ddp': DDPHyperbolicTangentRelaxation,
+        'deo': DEOHyperbolicTangentRelaxation,
+    }
+
+    if task_names is not None:
+        return [objectives[n]("labels_{}".format(t), "logits_{}".format(t)) for n, t in zip(names, task_names)]
+    else:
+        return [ objectives[n]() for n in names ]
+
+
+class CrossEntropyLoss(torch.nn.CrossEntropyLoss):
+    def __init__(self, label_name='labels', logits_name='logits'):
+        super().__init__(reduction='mean')
+        self.label_name = label_name
+        self.logits_name = logits_name
+
+    def __call__(self, **kwargs):
+        logits = kwargs[self.logits_name]
+        labels = kwargs[self.label_name]
+        return super().__call__(logits, labels)
+
+
+
+class BinaryCrossEntropyLoss(torch.nn.BCEWithLogitsLoss):
+
+    def __init__(self, label_name='labels', logits_name='logits', pos_weight=None):
+        super().__init__(reduction='mean', pos_weight=torch.Tensor([pos_weight]).cuda() if pos_weight else None)
+        self.label_name = label_name
+        self.logits_name = logits_name
+
+    def __call__(self, **kwargs):
+        logits = kwargs[self.logits_name]
+        labels = kwargs[self.label_name]
+        if logits.ndim == 2:
+            logits = torch.squeeze(logits)
+        if labels.dtype != torch.float:
+            labels = labels.float()
+        return super().__call__(logits, labels)
+
+
+
+
+
+class MSELoss(torch.nn.MSELoss):
+
+    def __init__(self, label_name='labels'):
+        super().__init__()
+        self.label_name = label_name
+
+    def __call__(self, **kwargs):
+        logits = kwargs['logits']
+        labels = kwargs[self.label_name]
+        if logits.ndim == 2:
+            logits = torch.squeeze(logits)
+        return super().__call__(logits, labels)
+
+
+class L1Regularization():
+
+    def __call__(self, **kwargs):
+        model = kwargs['model']
+        return torch.linalg.norm(torch.cat([p.view(-1) for p in model.parameters()]), ord=1)
+
+
+class L2Regularization():
+
+    def __call__(self, **kwargs):
+        model = kwargs['model']
+        return torch.linalg.norm(torch.cat([p.view(-1) for p in model.parameters()]), ord=2)
+
+
+class DDPHyperbolicTangentRelaxation():
+
+    def __init__(self, label_name='labels', logits_name='logits', s_name='sensible_attribute', c=1):
+        self.label_name = label_name
+        self.logits_name = logits_name
+        self.s_name = s_name
+        self.c = c
+
+    def __call__(self, **kwargs):
+        logits = kwargs[self.logits_name]
+        labels = kwargs[self.label_name]
+        sensible_attribute = kwargs[self.s_name]
+
+        n = logits.shape[0]
+        logits = torch.sigmoid(logits)
+        s_negative = logits[sensible_attribute.bool()]
+        s_positive = logits[~sensible_attribute.bool()]
+
+        return 1 / n * torch.abs(torch.sum(torch.tanh(self.c * torch.relu(s_positive))) - torch.sum(
+            torch.tanh(self.c * torch.relu(s_negative))))
+
+
+class DEOHyperbolicTangentRelaxation():
+
+    def __init__(self, label_name='labels', logits_name='logits', s_name='sensible_attribute', c=1):
+        self.label_name = label_name
+        self.logits_name = logits_name
+        self.s_name = s_name
+        self.c = c
+
+    def __call__(self, **kwargs):
+        logits = kwargs[self.logits_name]
+        labels = kwargs[self.label_name]
+        sensible_attribute = kwargs[self.s_name]
+
+        n = logits.shape[0]
+        logits = torch.sigmoid(logits)
+        s_negative = logits[(sensible_attribute.bool()) & (labels == 1)]
+        s_positive = logits[(~sensible_attribute.bool()) & (labels == 1)]
+
+        return 1 / n * torch.abs(torch.sum(torch.tanh(self.c * torch.relu(s_positive))) - torch.sum(
+            torch.tanh(self.c * torch.relu(s_negative))))
+
+
+"""
+Popular problem proposed by
+
+    Carlos Manuel Mira da Fonseca. Multiobjective genetic algorithms with 
+    application to controlengineering problems.PhD thesis, University of Sheffield, 1995.
+
+with a concave pareto front.
+
+$ \mathcal{L}_1(\theta) = 1 - \exp{ - || \theta - 1 / \sqrt{d} || $
+$ \mathcal{L}_1(\theta) = 1 - \exp{ - || \theta + 1 / \sqrt{d} || $
+
+with $\theta \in R^d$ and $ d = 100$
+"""
+
+
+#
+# class Fonseca1():
+#
+#     def f1(theta):
+#         d = len(theta)
+#         sum1 = autograd.numpy.sum([(theta[i] - 1.0 / autograd.numpy.sqrt(d)) ** 2 for i in range(d)])
+#         f1 = 1 - autograd.numpy.exp(- sum1)
+#         return f1
+#
+#     f1_dx = autograd.grad(f1)
+#
+#     def __call__(self, **kwargs):
+#         return f1(kwargs['parameters'])
+#
+#     def gradient(self, **kwargs):
+#         return f1_dx(kwargs['parameters'])
+#
+#
+# class Fonseca2():
+#
+#     def f2(theta):
+#         d = len(theta)
+#         sum1 = autograd.numpy.sum([(theta[i] + 1.0 / autograd.numpy.sqrt(d)) ** 2 for i in range(d)])
+#         f1 = 1 - autograd.numpy.exp(- sum1)
+#         return f1
+#
+#     f2_dx = autograd.grad(f2)
+#
+#     def __call__(self, **kwargs):
+#         return f2(kwargs['parameters'])
+#
+#     def gradient(self, **kwargs):
+#         return f2_dx(kwargs['parameters'])
+
+
+
+
+
+#
+# if __name__ == '__main__':
+#     problem = Fonseca1()

bike_bench_internal/benchmark_models/libmoon/problem/synthetic/__init__.py

@@ -0,0 +1,4 @@
+from .vlmop import VLMOP1, VLMOP2
+from .zdt import ZDT1, ZDT2, ZDT3, ZDT4, ZDT6
+from .maf import MAF1
+

bike_bench_internal/benchmark_models/libmoon/problem/synthetic/dtlz.py

@@ -0,0 +1,150 @@
+
+import numpy as np
+import torch
+from ..mop import mop
+
+
+
+
+class DTLZ1(mop):
+
+    def __init__(self, n_var=30, n_obj=3, lbound=np.zeros(30),
+                 ubound=np.ones(30)):
+
+        super().__init__(n_var=n_var,
+                         n_obj=n_obj,
+                         lbound=lbound,
+                         ubound=ubound)
+        self.problem_name= 'DTLZ1'
+
+    def _evaluate_torch(self, x: torch.Tensor):
+        x1 = x[:,0]
+        x2 = x[:,1]
+        xm = x[:, 2:]
+
+        g = 100 * (self.n_var - 2 + torch.sum(torch.pow(xm - 0.5, 2) -
+                                              torch.cos(20 * np.pi * (xm - 0.5)), dim=1))
+
+        f1 = 0.5 * x1 * x2 * (1+g)
+        f2 = 0.5 * x1 * (1 - x2) * (1+g)
+        f3 = 0.5 * (1 - x1) * (1+g)
+        return torch.stack((f1, f2, f3), dim=1)
+
+    def _evaluate_numpy(self, x: np.ndarray):
+        x1 = x[:, 0]
+        x2 = x[:, 1]
+        xm = x[:, 2:]
+
+        g = 100 * (self.n_var - 2 + np.sum(np.power(xm - 0.5, 2) -
+                                              np.cos(20 * np.pi * (xm - 0.5)), axis=1))
+
+        f1 = 0.5 * x1 * x2 * (1+g)
+        f2 = 0.5 * x1 * (1 - x2) * (1+g)
+        f3 = 0.5 * (1 - x1) * (1+g)
+        return np.stack((f1, f2, f3), axis=1)
+
+
+
+class DTLZ2(mop):
+    def __init__(self, n_var=30, n_obj=3, lbound=np.zeros(30),
+                 ubound=np.ones(30)):
+        super().__init__(n_var=n_var,
+                         n_obj=n_obj,
+                         lbound=lbound,
+                         ubound=ubound, )
+        self.problem_name = 'DTLZ2'
+
+    def _evaluate_torch(self, x):
+        xm = x[:, 2:]
+        g = torch.sum(torch.pow(xm - 0.5, 2), dim=1)
+        f1 = torch.cos(x[:, 0] * np.pi / 2) * torch.cos(x[:, 1] * np.pi / 2) * (1 + g)
+        f2 = torch.cos(x[:, 0] * np.pi / 2) * torch.sin(x[:, 1] * np.pi / 2) * (1 + g)
+        f3 = torch.sin(x[:, 0] * np.pi / 2) * (1 + g)
+        return torch.stack((f1, f2, f3), dim=1)
+
+    def _evaluate_numpy(self, x):
+        xm = x[:, 2:]
+        g = np.sum(np.power(xm - 0.5, 2), axis=1)
+        f1 = np.cos(x[:, 0] * np.pi / 2) * np.cos(x[:, 1] * np.pi / 2) * (1 + g)
+        f2 = np.cos(x[:, 0] * np.pi / 2) * np.sin(x[:, 1] * np.pi / 2) * (1 + g)
+        f3 = np.sin(x[:, 0] * np.pi / 2) * (1 + g)
+        return np.stack((f1, f2, f3), axis=1)
+
+
+class DTLZ3(mop):
+    def __init__(self, n_var=30, n_obj=3, lbound=np.zeros(30),
+                 ubound=np.ones(30)):
+        super().__init__(n_var=n_var,
+                         n_obj=n_obj,
+                         lbound=lbound,
+                         ubound=ubound, )
+        self.problem_name = 'DTLZ3'
+
+
+    def _evaluate_torch(self, x):
+        xm = x[:, 2:]
+        g = 100 * (self.n_var - 2 + torch.sum(torch.pow(xm - 0.5, 2) -
+                                              torch.cos(20 * np.pi * (xm - 0.5)), dim=1))
+        f1 = torch.cos(x[:, 0] * np.pi / 2) * torch.cos(x[:, 1] * np.pi / 2) * (1 + g)
+        f2 = torch.cos(x[:, 0] * np.pi / 2) * torch.sin(x[:, 1] * np.pi / 2) * (1 + g)
+        f3 = torch.sin(x[:, 0] * np.pi / 2) * (1 + g)
+        return torch.stack((f1, f2, f3), dim=1)
+
+    def _evaluate_numpy(self, x):
+        xm = x[:, 2:]
+
+        g = 100 * (self.n_var - 2 + np.sum(np.power(xm - 0.5, 2) -
+                                              np.cos(20 * np.pi * (xm - 0.5)), axis=1))
+
+        f1 = np.cos(x[:, 0] * np.pi / 2) * np.cos(x[:, 1] * np.pi / 2) * (1 + g)
+        f2 = np.cos(x[:, 0] * np.pi / 2) * np.sin(x[:, 1] * np.pi / 2) * (1 + g)
+        f3 = np.sin(x[:, 0] * np.pi / 2) * (1 + g)
+        return np.stack((f1, f2, f3), axis=1)
+
+
+class DTLZ4(mop):
+    def __init__(self, n_var=30, n_obj=3, lbound=np.zeros(30),
+                 ubound=np.ones(30)):
+        super().__init__(n_var=n_var,
+                         n_obj=n_obj,
+                         lbound=lbound,
+                         ubound=ubound, )
+        self.problem_name = 'DTLZ4'
+        self.alpha = 20
+
+    def _evaluate_torch(self, x):
+        xm = x[:, 2:]
+        g = torch.sum(torch.pow(xm - 0.5, 2), dim=1)
+        # alpha = 1
+
+        f1 = torch.cos(x[:, 0] ** self.alpha * np.pi / 2) * torch.cos(x[:, 1] ** self.alpha * np.pi / 2) * (1 + g)
+        f2 = torch.cos(x[:, 0] ** self.alpha * np.pi / 2) * torch.sin(x[:, 1] ** self.alpha * np.pi / 2) * (1 + g)
+        f3 = torch.sin(x[:, 0] ** self.alpha * np.pi / 2) * (1 + g)
+        return torch.stack((f1, f2, f3), dim=1)
+
+    def _evaluate_numpy(self, x):
+        xm = x[:, 2:]
+        g = np.sum(np.power(xm - 0.5, 2), axis=1)
+
+        f1 = np.cos(x[:, 0] ** self.alpha * np.pi / 2) * np.cos(x[:, 1] ** self.alpha * np.pi / 2) * (1 + g)
+        f2 = np.cos(x[:, 0] ** self.alpha * np.pi / 2) * np.sin(x[:, 1] ** self.alpha * np.pi / 2) * (1 + g)
+        f3 = np.sin(x[:, 0] ** self.alpha * np.pi / 2) * (1 + g)
+        return np.stack((f1, f2, f3), axis=1 )
+
+# DTLZ5, DTLZ6.
+# degenerated.
+
+
+# DTLZ7 has disjoint Pareto front.
+
+
+
+if __name__ == '__main__':
+    x = torch.rand(100, 30)
+    problem = DTLZ4()
+
+    y = problem.evaluate(x)
+    print( y )
+
+
+

bike_bench_internal/benchmark_models/libmoon/problem/synthetic/maf.py

@@ -0,0 +1,44 @@
+
+from numpy import array
+import torch
+
+
+
+class MAF1:
+    def __init__(self):
+        '''
+            n_obj can be set as any number. For simlicity, we set it as 3.
+        '''
+        self.n_obj = 3
+        self.n_var = 30
+        self.lb = 0
+        self.ub = 1
+
+    def evaluate(self, x):
+        if type(x) == torch.Tensor:
+
+            g = torch.sum( torch.pow(x[:, 2:] - 0.5, 2), dim=1 )
+
+            f1 = (1 - x[:,0] * x[:,1]) * (1 + g)
+            f2 = (1 - x[:,0] * (1 - x[:,1]) ) * (1 + g)
+            f3 = x[:,0] * (1 + g)
+
+            return torch.stack((f1, f2, f3), dim=1)
+
+        else:
+            assert False
+
+
+
+    def get_pf(self):
+        return array([[0.0, 0.0, 0.0]])
+
+
+
+if __name__ == '__main__':
+    x = torch.rand(100, 30)
+    problem = MAF1()
+
+    y = problem.evaluate(x)
+    print()
+

bike_bench_internal/benchmark_models/libmoon/problem/synthetic/re.py

@@ -0,0 +1,619 @@
+
+import numpy as np
+import torch
+from ..mop import mop
+
+from numpy import array
+
+
+class RE21(mop):
+    def __init__(self, n_var=4, n_obj=2, lbound=np.zeros(4), ubound=np.ones(4)):
+        self.problem_name = 'RE21'
+        self.n_var = n_var
+        self.n_obj = n_obj
+        self.n_cons = 0
+
+        self.n_original_constraints = 0
+
+        self.ideal = array([1237.8414230005742, 0.002761423749158419])
+        # self.nadir = array([2086.36956042, 0.00341421356237])
+        self.nadir = np.array([2886.3695604236013, 0.039999999999998245])
+
+        F = 10.0
+        sigma = 10.0
+        tmp_val = F / sigma
+        self.ubound = np.full(self.n_var, 3 * tmp_val)
+        self.lbound = np.zeros(self.n_var)
+        self.lbound[0] = tmp_val
+        self.lbound[1] = np.sqrt(2.0) * tmp_val
+        self.lbound[2] = np.sqrt(2.0) * tmp_val
+        self.lbound[3] = tmp_val
+
+
+    def _evaluate_numpy(self, x):
+        n_sub = len(x)
+
+        x1 = x[:,0]
+        x2 = x[:,1]
+        x3 = x[:,2]
+        x4 = x[:,3]
+        f = np.zeros((n_sub, self.n_obj) )
+
+        F = 10.0
+        sigma = 10.0
+        E = 2.0 * 1e5
+        L = 200.0
+
+        f[:,0] = L * ((2 * x1) + np.sqrt(2.0) * x2 + np.sqrt(x3) + x4)
+        f[:,1] = ((F * L) / E) * ((2.0 / x1) + (2.0 * np.sqrt(2.0) / x2) - (2.0 * np.sqrt(2.0) / x3) + (2.0 / x4))
+
+        # f_arr = np.stack((f1,f2), axis=1)
+
+        f_arr_norm = (f - self.ideal) / (self.nadir - self.ideal)
+        # f_arr_norm =
+        f_arr_norm[:, 0] = 0.5 * f_arr_norm[:, 0]
+
+        return f_arr_norm
+
+
+    def _evaluate_torch(self, x):
+        x1 = x[:, 0]
+        x2 = x[:, 1]
+        x3 = x[:, 2]
+        x4 = x[:, 3]
+
+        F = 10.0
+        sigma = 10.0
+        E = 2.0 * 1e5
+        L = 200.0
+
+        f1 = L * ( (2 * x1) + np.sqrt(2.0) * x2 + torch.sqrt(x3) + x4 )
+        f2 = ((F * L) / E) * ((2.0 / x1) + (2.0 * np.sqrt(2.0) / x2) - (2.0 * np.sqrt(2.0) / x3) + (2.0 / x4))
+        f_arr = torch.stack((f1, f2), dim=1)
+        f_arr_norm = (f_arr - self.ideal) / (self.nadir - self.ideal)
+        f_arr_norm[:, 0] = 0.5 * f_arr_norm[:, 0]
+        return f_arr_norm
+
+
+class RE22(mop):
+    def __init__(self, n_var=3, n_obj=2, lbound=np.zeros(30),
+                 ubound=np.ones(30)):
+
+
+        self.n_var=n_var
+        self.n_obj=n_obj
+        self.problem_name = 'RE22'
+        self.n_cons = 0
+        self.n_original_constraints = 2
+
+        self.ideal = np.array([5.88, 0.0])
+        self.nadir = np.array([361.262944647, 180.01547])
+
+
+        self.ubound = np.zeros(self.n_var)
+        self.lbound = np.zeros(self.n_var)
+
+        self.lbound[0] = 0.2
+        self.lbound[1] = 0.0
+        self.lbound[2] = 0.0
+        self.ubound[0] = 15
+        self.ubound[1] = 20
+        self.ubound[2] = 40
+
+        self.n_var = n_var
+        self.n_obj = n_obj
+
+        self.feasible_vals = np.array(
+            [0.20, 0.31, 0.40, 0.44, 0.60, 0.62, 0.79, 0.80, 0.88, 0.93, 1.0, 1.20, 1.24, 1.32, 1.40, 1.55, 1.58, 1.60,
+             1.76, 1.80, 1.86, 2.0, 2.17, 2.20, 2.37, 2.40, 2.48, 2.60, 2.64, 2.79, 2.80, 3.0, 3.08, 3, 10, 3.16, 3.41,
+             3.52, 3.60, 3.72, 3.95, 3.96, 4.0, 4.03, 4.20, 4.34, 4.40, 4.65, 4.74, 4.80, 4.84, 5.0, 5.28, 5.40, 5.53,
+             5.72, 6.0, 6.16, 6.32, 6.60, 7.11, 7.20, 7.80, 7.90, 8.0, 8.40, 8.69, 9.0, 9.48, 10.27, 11.0, 11.06, 11.85,
+             12.0, 13.0, 14.0, 15.0])
+
+
+    def _evaluate_numpy(self, x):
+        n_sub = len(x)
+
+        f = np.zeros( (n_sub, self.n_obj) )
+        g = np.zeros( (n_sub, self.n_original_constraints)  )
+        # Reference: getNearestValue_sample2.py (https://gist.github.com/icchi-h/1d0bb1c52ebfdd31f14b3e811328390a)
+        idx_arr = [np.abs(np.asarray(self.feasible_vals) - x0).argmin() for x0 in x[:,0]]
+        x1 = array([self.feasible_vals[idx] for idx in idx_arr])
+        x2 = x[:,1]
+        x3 = x[:,2]
+
+        # First original objective function
+        f[:,0] = (29.4 * x1) + (0.6 * x2 * x3)
+        # Original constraint functions
+        g[:,0] = (x1 * x3) - 7.735 * ((x1 * x1) / x2) - 180.0
+        g[:,1] = 4.0 - (x3 / x2)
+        g = np.where(g < 0, -g, 0)
+        f[:,1] = g[:,0] + g[:,1]
+        f_norm = (f - self.ideal) / (self.nadir - self.ideal)
+        f_norm[:, 0] = 0.5 * f_norm[:, 0]
+
+        return f_norm
+
+    def _evaluate_torch(self, x):
+        pass
+
+
+class RE23(mop):
+    def __init__(self, n_var=4, n_obj=2, lbound=np.zeros(2),
+                 ubound=np.ones(2)):
+        self.problem_name = 'RE23'
+        self.n_obj = n_obj
+        self.n_var = n_var
+        self.n_cons = 0
+        self.n_original_constraints = 3
+
+        self.ideal = array([15.9018007813, 0.0])
+        self.nadir = array([481.608088535, 44.2819047619])
+
+        self.ubound = np.zeros(self.n_var)
+        self.lbound = np.zeros(self.n_var)
+        self.lbound[0] = 1
+        self.lbound[1] = 1
+        self.lbound[2] = 10
+        self.lbound[3] = 10
+        self.ubound[0] = 100
+        self.ubound[1] = 100
+        self.ubound[2] = 200
+        self.ubound[3] = 240
+
+    def _evaluate_numpy(self, x):
+
+        f = np.zeros( (len(x), self.n_obj) )
+        g = np.zeros( (len(x), self.n_original_constraints))
+
+        x1 = 0.0625 * np.round(x[:,0]).astype(np.int32)
+        x2 = 0.0625 * np.round(x[:,1]).astype(np.int32)
+
+        x3 = x[:,2]
+        x4 = x[:,3]
+
+        # First original objective function
+        f[:,0] = (0.6224 * x1 * x3 * x4) + (1.7781 * x2 * x3 * x3) + (3.1661 * x1 * x1 * x4) + (19.84 * x1 * x1 * x3)
+
+        # Original constraint functions
+        g[:,0] = x1 - (0.0193 * x3)
+        g[:,1] = x2 - (0.00954 * x3)
+        g[:,2] = (np.pi * x3 * x3 * x4) + ((4.0 / 3.0) * (np.pi * x3 * x3 * x3)) - 1296000
+        g = np.where(g < 0, -g, 0)
+        f[:,1] = g[:,0] + g[:,1] + g[:,2]
+
+        f_norm = (f - self.ideal) / (self.nadir - self.ideal)
+
+        return f_norm
+
+
+
+
+
+class RE24(mop):
+    def __init__(self, n_var=2, n_obj=2, lbound=np.zeros(2),
+                 ubound=np.ones(2)):
+        super().__init__(n_var=n_var,
+                         n_obj=n_obj,
+                         lbound=lbound,
+                         ubound=ubound, )
+
+        self.problem_name = 'RE24'
+        self.n_obj = 2
+        self.n_var = 2
+
+        self.n_cons = 0
+        self.n_original_constraints = 4
+
+        self.ubound = np.zeros(self.n_var)
+        self.lbound = np.zeros(self.n_var)
+
+        self.lbound[0] = 0.5
+        self.lbound[1] = 0.5
+
+        self.ubound[0] = 4
+        self.ubound[1] = 50
+
+        self.ideal = np.array([60.5, 0.0])
+        self.nadir = np.array([481.608088535, 44.2819047619])
+
+
+
+
+    def _evaluate_numpy(self, x):
+        n_sub = len(x)
+        # f = np.zeros(self.n_objectives)
+        g = np.zeros( (n_sub, self.n_original_constraints) )
+
+        x1 = x[:,0]
+        x2 = x[:,1]
+
+        # First original objective function
+        f1 = x1 + (120 * x2)
+
+        E = 700000
+        sigma_b_max = 700
+        tau_max = 450
+        delta_max = 1.5
+        sigma_k = (E * x1 * x1) / 100
+        sigma_b = 4500 / (x1 * x2)
+        tau = 1800 / x2
+        delta = (56.2 * 10000) / (E * x1 * x2 * x2)
+
+        g[:,0] = 1 - (sigma_b / sigma_b_max)
+        g[:,1] = 1 - (tau / tau_max)
+        g[:,2] = 1 - (delta / delta_max)
+        g[:,3] = 1 - (sigma_b / sigma_k)
+        g = np.where(g < 0, -g, 0)
+        f2 = g[:,0] + g[:,1] + g[:,2] + g[:,3]
+
+        f_arr = np.stack((f1, f2), axis=1)
+        f_norm = (f_arr - self.ideal) / (self.nadir - self.ideal)
+
+        return f_norm
+
+
+    def _evaluate_torch(self, x):
+        pass
+
+
+class RE25(mop):
+    def __init__(self, n_var=3, n_obj=2):
+        self.problem_name = 'RE25'
+        self.n_obj = n_obj
+        self.n_var = n_var
+
+        self.n_cons = 0
+        self.n_original_constraints = 6
+
+        self.ideal = array([0.037591349242869145, 0.0])
+        self.nadir = array([0.40397042546, 2224669.22419])
+
+        self.ubound = np.zeros( self.n_var )
+        self.lbound = np.zeros( self.n_var )
+        self.lbound[0] = 1
+        self.lbound[1] = 0.6
+        self.lbound[2] = 0.09
+        self.ubound[0] = 70
+        self.ubound[1] = 3
+        self.ubound[2] = 0.5
+
+        self.feasible_vals = np.array(
+            [0.009, 0.0095, 0.0104, 0.0118, 0.0128, 0.0132, 0.014, 0.015, 0.0162, 0.0173, 0.018, 0.02, 0.023, 0.025,
+             0.028, 0.032, 0.035, 0.041, 0.047, 0.054, 0.063, 0.072, 0.08, 0.092, 0.105, 0.12, 0.135, 0.148, 0.162,
+             0.177, 0.192, 0.207, 0.225, 0.244, 0.263, 0.283, 0.307, 0.331, 0.362, 0.394, 0.4375, 0.5])
+
+    def _evaluate_numpy(self, x):
+        n_sub = len(x)
+        f = np.zeros( (n_sub, self.n_obj) )
+        g = np.zeros( (n_sub, self.n_original_constraints) )
+        x1 = np.round(x[:,0])
+        x2 = x[:,1]
+
+        # Reference: getNearestValue_sample2.py (https://gist.github.com/icchi-h/1d0bb1c52ebfdd31f14b3e811328390a)
+        idx_array = array([np.abs(np.asarray(self.feasible_vals) - x2).argmin() for x2 in x[:,2]])
+
+        x3 = array( [self.feasible_vals[idx] for idx in idx_array] )
+
+        # first original objective function
+        f[:,0] = (np.pi * np.pi * x2 * x3 * x3 * (x1 + 2)) / 4.0
+
+        # constraint functions
+        Cf = ((4.0 * (x2 / x3) - 1) / (4.0 * (x2 / x3) - 4)) + (0.615 * x3 / x2)
+        Fmax = 1000.0
+        S = 189000.0
+        G = 11.5 * 1e+6
+        K = (G * x3 * x3 * x3 * x3) / (8 * x1 * x2 * x2 * x2)
+        lmax = 14.0
+        lf = (Fmax / K) + 1.05 * (x1 + 2) * x3
+        dmin = 0.2
+        Dmax = 3
+        Fp = 300.0
+        sigmaP = Fp / K
+        sigmaPM = 6
+        sigmaW = 1.25
+
+        g[:,0] = -((8 * Cf * Fmax * x2) / (np.pi * x3 * x3 * x3)) + S
+        g[:,1] = -lf + lmax
+        g[:,2] = -3 + (x2 / x3)
+        g[:,3] = -sigmaP + sigmaPM
+        g[:,4] = -sigmaP - ((Fmax - Fp) / K) - 1.05 * (x1 + 2) * x3 + lf
+        g[:,5] = sigmaW - ((Fmax - Fp) / K)
+
+        g = np.where(g < 0, -g, 0)
+        f[:,1] = g[:,0] + g[:,1] + g[:,2] + g[:,3] + g[:,4] + g[:,5]
+
+        f_norm = (f - self.ideal) / (self.nadir - self.ideal)
+        return f_norm
+
+    def _evaluate_torch(self, x):
+        pass
+
+
+class RE31(mop):
+    def __init__(self, n_obj=3, n_var=3):
+        self.problem_name = 'RE31'
+        self.n_obj = n_obj
+        self.n_var = n_var
+        self.n_cons = 0
+        self.n_original_constraints = 3
+
+        self.ubound = np.zeros(self.n_var)
+        self.lbound = np.zeros(self.n_var)
+        self.lbound[0] = 0.00001
+        self.lbound[1] = 0.00001
+        self.lbound[2] = 1.0
+        self.ubound[0] = 100.0
+        self.ubound[1] = 100.0
+        self.ubound[2] = 3.0
+
+        self.ideal = np.array([5.53731918799e-05, 0.333333333333, 0.0])
+        self.nadir = np.array([500.002668442, 8246211.25124, 19359919.7502])
+
+
+    def _evaluate_numpy(self, x):
+        n_sub = len(x)
+        f = np.zeros( (n_sub, self.n_obj) )
+        g = np.zeros( (n_sub, self.n_original_constraints) )
+
+        x1 = x[:,0]
+        x2 = x[:,1]
+        x3 = x[:,2]
+
+        # First original objective function
+        f[:,0] = x1 * np.sqrt(16.0 + (x3 * x3)) + x2 * np.sqrt(1.0 + x3 * x3)
+        # Second original objective function
+        f[:,1] = (20.0 * np.sqrt(16.0 + (x3 * x3))) / (x1 * x3)
+
+        # Constraint functions
+        g[:,0] = 0.1 - f[:,0]
+        g[:,1] = 100000.0 - f[:,1]
+        g[:,2] = 100000 - ((80.0 * np.sqrt(1.0 + x3 * x3)) / (x3 * x2))
+        g = np.where(g < 0, -g, 0)
+        f[:,2] = g[:,0] + g[:,1] + g[:,2]
+
+        f_norm = (f - self.ideal) / (self.nadir - self.ideal)
+        return f_norm
+
+    def _evaluate_torch(self, x):
+        pass
+
+
+class RE37(mop):
+    def __init__(self, n_obj=3, n_var=4):
+        self.problem_name = 'RE37'
+        self.n_obj = n_obj
+        self.n_var = n_var
+        self.n_cons = 0
+        self.n_original_constraints = 0
+
+        self.lbound = np.full(self.n_var, 0)
+        self.ubound = np.full(self.n_var, 1)
+
+        self.ideal = np.array([0.00889341391106, 0.00488, -0.431499999825])
+        self.nadir = np.array([0.98949120096, 0.956587924661, 0.987530948586])
+
+
+    def _evaluate_numpy(self, x):
+        n_sub = len(x)
+        f = np.zeros( (n_sub, self.n_obj) )
+
+        xAlpha = x[:,0]
+        xHA = x[:,1]
+        xOA = x[:,2]
+        xOPTT = x[:,3]
+
+        # f1 (TF_max)
+        f[:,0] = 0.692 + (0.477 * xAlpha) - (0.687 * xHA) - (0.080 * xOA) - (0.0650 * xOPTT) - (
+                    0.167 * xAlpha * xAlpha) - (0.0129 * xHA * xAlpha) + (0.0796 * xHA * xHA) - (
+                           0.0634 * xOA * xAlpha) - (0.0257 * xOA * xHA) + (0.0877 * xOA * xOA) - (
+                           0.0521 * xOPTT * xAlpha) + (0.00156 * xOPTT * xHA) + (0.00198 * xOPTT * xOA) + (
+                           0.0184 * xOPTT * xOPTT)
+        # f2 (X_cc)
+        f[:,1] = 0.153 - (0.322 * xAlpha) + (0.396 * xHA) + (0.424 * xOA) + (0.0226 * xOPTT) + (
+                    0.175 * xAlpha * xAlpha) + (0.0185 * xHA * xAlpha) - (0.0701 * xHA * xHA) - (
+                           0.251 * xOA * xAlpha) + (0.179 * xOA * xHA) + (0.0150 * xOA * xOA) + (
+                           0.0134 * xOPTT * xAlpha) + (0.0296 * xOPTT * xHA) + (0.0752 * xOPTT * xOA) + (
+                           0.0192 * xOPTT * xOPTT)
+        # f3 (TT_max)
+        f[:,2] = 0.370 - (0.205 * xAlpha) + (0.0307 * xHA) + (0.108 * xOA) + (1.019 * xOPTT) - (
+                    0.135 * xAlpha * xAlpha) + (0.0141 * xHA * xAlpha) + (0.0998 * xHA * xHA) + (
+                           0.208 * xOA * xAlpha) - (0.0301 * xOA * xHA) - (0.226 * xOA * xOA) + (
+                           0.353 * xOPTT * xAlpha) - (0.0497 * xOPTT * xOA) - (0.423 * xOPTT * xOPTT) + (
+                           0.202 * xHA * xAlpha * xAlpha) - (0.281 * xOA * xAlpha * xAlpha) - (
+                           0.342 * xHA * xHA * xAlpha) - (0.245 * xHA * xHA * xOA) + (0.281 * xOA * xOA * xHA) - (
+                           0.184 * xOPTT * xOPTT * xAlpha) - (0.281 * xHA * xAlpha * xOA)
+
+        f_norm = (f - self.ideal) / (self.nadir - self.ideal)
+        return f_norm
+
+
+
+class RE41(mop):
+    def __init__(self, n_obj=4, n_var=7):
+        self.problem_name = 'RE41'
+        self.n_obj = n_obj
+        self.n_var = n_var
+        self.n_cons = 0
+        self.n_original_constraints = 10
+
+        self.lbound = np.zeros(self.n_var)
+        self.ubound = np.zeros(self.n_var)
+        self.lbound[0] = 0.5
+        self.lbound[1] = 0.45
+        self.lbound[2] = 0.5
+        self.lbound[3] = 0.5
+        self.lbound[4] = 0.875
+        self.lbound[5] = 0.4
+        self.lbound[6] = 0.4
+        self.ubound[0] = 1.5
+        self.ubound[1] = 1.35
+        self.ubound[2] = 1.5
+        self.ubound[3] = 1.5
+        self.ubound[4] = 2.625
+        self.ubound[5] = 1.2
+        self.ubound[6] = 1.2
+
+        self.ideal = np.array([15.576004, 3.58525, 10.61064375, 0.0])
+        self.nadir = np.array([39.2905121788, 4.42725, 13.09138125, 9.49401929991])
+
+
+    def _evaluate_numpy(self, x):
+        n_sub = len(x)
+
+        f = np.zeros( (n_sub, self.n_obj) )
+        g = np.zeros( (n_sub, self.n_original_constraints) )
+
+        x1 = x[:,0]
+        x2 = x[:,1]
+        x3 = x[:,2]
+        x4 = x[:,3]
+        x5 = x[:,4]
+        x6 = x[:,5]
+        x7 = x[:,6]
+
+        # First original objective function
+        f[:,0] = 1.98 + 4.9 * x1 + 6.67 * x2 + 6.98 * x3 + 4.01 * x4 + 1.78 * x5 + 0.00001 * x6 + 2.73 * x7
+        # Second original objective function
+        f[:,1] = 4.72 - 0.5 * x4 - 0.19 * x2 * x3
+        # Third original objective function
+        Vmbp = 10.58 - 0.674 * x1 * x2 - 0.67275 * x2
+        Vfd = 16.45 - 0.489 * x3 * x7 - 0.843 * x5 * x6
+        f[:,2] = 0.5 * (Vmbp + Vfd)
+
+        # Constraint functions
+        g[:,0] = 1 - (1.16 - 0.3717 * x2 * x4 - 0.0092928 * x3)
+        g[:,1] = 0.32 - (0.261 - 0.0159 * x1 * x2 - 0.06486 * x1 - 0.019 * x2 * x7 + 0.0144 * x3 * x5 + 0.0154464 * x6)
+        g[:,2] = 0.32 - (
+                    0.214 + 0.00817 * x5 - 0.045195 * x1 - 0.0135168 * x1 + 0.03099 * x2 * x6 - 0.018 * x2 * x7 + 0.007176 * x3 + 0.023232 * x3 - 0.00364 * x5 * x6 - 0.018 * x2 * x2)
+        g[:,3] = 0.32 - (0.74 - 0.61 * x2 - 0.031296 * x3 - 0.031872 * x7 + 0.227 * x2 * x2)
+        g[:,4] = 32 - (28.98 + 3.818 * x3 - 4.2 * x1 * x2 + 1.27296 * x6 - 2.68065 * x7)
+        g[:,5] = 32 - (33.86 + 2.95 * x3 - 5.057 * x1 * x2 - 3.795 * x2 - 3.4431 * x7 + 1.45728)
+        g[:,6] = 32 - (46.36 - 9.9 * x2 - 4.4505 * x1)
+        g[:,7] = 4 - f[:,1]
+        g[:,8] = 9.9 - Vmbp
+        g[:,9] = 15.7 - Vfd
+        g = np.where(g < 0, -g, 0)
+        f[:,3] = g[:,0] + g[:,1] + g[:,2] + g[:,3] + g[:,4] + g[:,5] + g[:,6] + g[:,7] + g[:,8] + g[:,9]
+        f_norm = (f - self.ideal) / (self.nadir - self.ideal)
+        return f_norm
+
+    def _evaluate_torch(self, x):
+        pass
+
+
+class RE42(mop):
+    def __init__(self):
+        self.problem_name = 'RE42'
+
+        self.n_obj = 4
+        self.n_var = 6
+        self.n_cons = 0
+        self.n_original_constraints = 9
+
+        self.lbound = np.zeros(self.n_var )
+        self.ubound = np.zeros(self.n_var )
+        self.lbound[0] = 150.0
+        self.lbound[1] = 20.0
+        self.lbound[2] = 13.0
+        self.lbound[3] = 10.0
+        self.lbound[4] = 14.0
+        self.lbound[5] = 0.63
+        self.ubound[0] = 274.32
+        self.ubound[1] = 32.31
+        self.ubound[2] = 25.0
+        self.ubound[3] = 11.71
+        self.ubound[4] = 18.0
+        self.ubound[5] = 0.75
+
+        self.ideal = np.array([-2756.2590400638524, 3962.557843228888, 1947.880856925791, 0.0])
+        self.nadir = np.array([-1010.5229595219643, 13827.138456300128, 2611.9668107424536, 12.437669929732023 ])
+
+
+
+    def _evaluate_numpy(self, x):
+        n_sub = len(x)
+
+        f = np.zeros( (n_sub, self.n_obj) )
+        # NOT g
+        constraintFuncs = np.zeros( (n_sub, self.n_original_constraints) )
+
+        x_L = x[:,0]
+        x_B = x[:,1]
+        x_D = x[:,2]
+        x_T = x[:,3]
+        x_Vk = x[:,4]
+        x_CB = x[:,5]
+
+        displacement = 1.025 * x_L * x_B * x_T * x_CB
+        V = 0.5144 * x_Vk
+        g = 9.8065
+        Fn = V / np.power(g * x_L, 0.5)
+        a = (4977.06 * x_CB * x_CB) - (8105.61 * x_CB) + 4456.51
+        b = (-10847.2 * x_CB * x_CB) + (12817.0 * x_CB) - 6960.32
+
+        power = (np.power(displacement, 2.0 / 3.0) * np.power(x_Vk, 3.0)) / (a + (b * Fn))
+        outfit_weight = 1.0 * np.power(x_L, 0.8) * np.power(x_B, 0.6) * np.power(x_D, 0.3) * np.power(x_CB, 0.1)
+        steel_weight = 0.034 * np.power(x_L, 1.7) * np.power(x_B, 0.7) * np.power(x_D, 0.4) * np.power(x_CB, 0.5)
+        machinery_weight = 0.17 * np.power(power, 0.9)
+        light_ship_weight = steel_weight + outfit_weight + machinery_weight
+
+        ship_cost = 1.3 * ((2000.0 * np.power(steel_weight, 0.85)) + (3500.0 * outfit_weight) + (
+                    2400.0 * np.power(power, 0.8)))
+        capital_costs = 0.2 * ship_cost
+
+        DWT = displacement - light_ship_weight
+
+        running_costs = 40000.0 * np.power(DWT, 0.3)
+
+        round_trip_miles = 5000.0
+        sea_days = (round_trip_miles / 24.0) * x_Vk
+        handling_rate = 8000.0
+
+        daily_consumption = ((0.19 * power * 24.0) / 1000.0) + 0.2
+        fuel_price = 100.0
+        fuel_cost = 1.05 * daily_consumption * sea_days * fuel_price
+        port_cost = 6.3 * np.power(DWT, 0.8)
+
+        fuel_carried = daily_consumption * (sea_days + 5.0)
+        miscellaneous_DWT = 2.0 * np.power(DWT, 0.5)
+
+        cargo_DWT = DWT - fuel_carried - miscellaneous_DWT
+        port_days = 2.0 * ((cargo_DWT / handling_rate) + 0.5)
+        RTPA = 350.0 / (sea_days + port_days)
+
+        voyage_costs = (fuel_cost + port_cost) * RTPA
+        annual_costs = capital_costs + running_costs + voyage_costs
+        annual_cargo = cargo_DWT * RTPA
+
+        f[:,0] = annual_costs / annual_cargo
+        f[:,1] = light_ship_weight
+        # f_2 is dealt as a minimization problem
+        f[:,2] = -annual_cargo
+
+        # Reformulated objective functions
+        constraintFuncs[:,0] = (x_L / x_B) - 6.0
+        constraintFuncs[:,1] = -(x_L / x_D) + 15.0
+        constraintFuncs[:,2] = -(x_L / x_T) + 19.0
+        constraintFuncs[:,3] = 0.45 * np.power(DWT, 0.31) - x_T
+        constraintFuncs[:,4] = 0.7 * x_D + 0.7 - x_T
+        constraintFuncs[:,5] = 500000.0 - DWT
+        constraintFuncs[:,6] = DWT - 3000.0
+        constraintFuncs[:,7] = 0.32 - Fn
+
+        KB = 0.53 * x_T
+        BMT = ((0.085 * x_CB - 0.002) * x_B * x_B) / (x_T * x_CB)
+        KG = 1.0 + 0.52 * x_D
+        constraintFuncs[:,8] = (KB + BMT - KG) - (0.07 * x_B)
+
+        constraintFuncs = np.where(constraintFuncs < 0, -constraintFuncs, 0)
+        f[:,3] = constraintFuncs[:,0] + constraintFuncs[:,1] + constraintFuncs[:,2] + constraintFuncs[:,3] + constraintFuncs[:,4] + \
+               constraintFuncs[:,5] + constraintFuncs[:,6] + constraintFuncs[:,7] + constraintFuncs[:,8]
+
+        f_norm = (f - self.ideal) / (self.nadir - self.ideal)
+        return f
+
+
+
+
+

bike_bench_internal/benchmark_models/libmoon/problem/synthetic/re_original.py

@@ -0,0 +1,1335 @@
+#!/usr/bin/env python
+"""
+  A real-world multimnist-objective problem suite (the RE benchmark set)
+  Reference:
+  Ryoji Tanabe, Hisao Ishibuchi, "An Easy-to-use Real-world Multi-objective Problem Suite" Applied Soft Computing. 89: 106078 (2020)
+   Copyright (c) 2020 Ryoji Tanabe
+
+  I re-implemented the RE problem set by referring to its C source code (reproblem.c). While variables directly copied from the C source code are written in CamelCase, the other variables are written in snake_case. It is somewhat awkward.
+
+  This program is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  This program is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+"""
+import numpy as np
+
+class RE21():
+    def __init__(self, n_var=4, n_obj=2, lower_bound=np.zeros(30),
+                 upper_bound=np.ones(30)):
+        self.problem_name = 'RE21'
+        self.n_constraints = 0
+        self.n_original_constraints = 0
+
+        F = 10.0
+        sigma = 10.0
+        tmp_val = F / sigma
+
+        self.ubound = np.full(self.n_variables, 3 * tmp_val)
+        self.lbound = np.zeros(self.n_variables)
+        self.lbound[0] = tmp_val
+        self.lbound[1] = np.sqrt(2.0) * tmp_val
+        self.lbound[2] = np.sqrt(2.0) * tmp_val
+        self.lbound[3] = tmp_val
+
+    def evaluate(self, x):
+        f = np.zeros(self.n_objectives)
+        x1 = x[0]
+        x2 = x[1]
+        x3 = x[2]
+        x4 = x[3]
+
+        F = 10.0
+        sigma = 10.0
+        E = 2.0 * 1e5
+        L = 200.0
+
+        f[0] = L * ((2 * x1) + np.sqrt(2.0) * x2 + np.sqrt(x3) + x4)
+        f[1] = ((F * L) / E) * ((2.0 / x1) + (2.0 * np.sqrt(2.0) / x2) - (2.0 * np.sqrt(2.0) / x3) + (2.0 / x4))
+
+        return f
+
+
+class RE22():
+    def __init__(self):
+        self.problem_name = 'RE22'
+        self.n_objectives = 2
+        self.n_variables = 3
+
+        self.n_constraints = 0
+        self.n_original_constraints = 2
+
+        self.ubound = np.zeros(self.n_variables)
+        self.lbound = np.zeros(self.n_variables)
+        self.lbound[0] = 0.2
+        self.lbound[1] = 0.0
+        self.lbound[2] = 0.0
+        self.ubound[0] = 15
+        self.ubound[1] = 20
+        self.ubound[2] = 40
+
+        self.feasible_vals = np.array(
+            [0.20, 0.31, 0.40, 0.44, 0.60, 0.62, 0.79, 0.80, 0.88, 0.93, 1.0, 1.20, 1.24, 1.32, 1.40, 1.55, 1.58, 1.60,
+             1.76, 1.80, 1.86, 2.0, 2.17, 2.20, 2.37, 2.40, 2.48, 2.60, 2.64, 2.79, 2.80, 3.0, 3.08, 3, 10, 3.16, 3.41,
+             3.52, 3.60, 3.72, 3.95, 3.96, 4.0, 4.03, 4.20, 4.34, 4.40, 4.65, 4.74, 4.80, 4.84, 5.0, 5.28, 5.40, 5.53,
+             5.72, 6.0, 6.16, 6.32, 6.60, 7.11, 7.20, 7.80, 7.90, 8.0, 8.40, 8.69, 9.0, 9.48, 10.27, 11.0, 11.06, 11.85,
+             12.0, 13.0, 14.0, 15.0])
+
+    def evaluate(self, x):
+        f = np.zeros(self.n_objectives)
+        g = np.zeros(self.n_original_constraints)
+        # Reference: getNearestValue_sample2.py (https://gist.github.com/icchi-h/1d0bb1c52ebfdd31f14b3e811328390a)
+        idx = np.abs(np.asarray(self.feasible_vals) - x[0]).argmin()
+        x1 = self.feasible_vals[idx]
+        x2 = x[1]
+        x3 = x[2]
+
+        # First original objective function
+        f[0] = (29.4 * x1) + (0.6 * x2 * x3)
+
+        # Original constraint functions
+        g[0] = (x1 * x3) - 7.735 * ((x1 * x1) / x2) - 180.0
+        g[1] = 4.0 - (x3 / x2)
+        g = np.where(g < 0, -g, 0)
+        f[1] = g[0] + g[1]
+
+        return f
+
+
+class RE23():
+    def __init__(self):
+        self.problem_name = 'RE23'
+        self.n_objectives = 2
+        self.n_variables = 4
+        self.n_constraints = 0
+        self.n_original_constraints = 3
+
+        self.ubound = np.zeros(self.n_variables)
+        self.lbound = np.zeros(self.n_variables)
+        self.lbound[0] = 1
+        self.lbound[1] = 1
+        self.lbound[2] = 10
+        self.lbound[3] = 10
+        self.ubound[0] = 100
+        self.ubound[1] = 100
+        self.ubound[2] = 200
+        self.ubound[3] = 240
+
+    def evaluate(self, x):
+        f = np.zeros(self.n_objectives)
+        g = np.zeros(self.n_original_constraints)
+
+        x1 = 0.0625 * int(np.round(x[0]))
+        x2 = 0.0625 * int(np.round(x[1]))
+        x3 = x[2]
+        x4 = x[3]
+
+        # First original objective function
+        f[0] = (0.6224 * x1 * x3 * x4) + (1.7781 * x2 * x3 * x3) + (3.1661 * x1 * x1 * x4) + (19.84 * x1 * x1 * x3)
+
+        # Original constraint functions
+        g[0] = x1 - (0.0193 * x3)
+        g[1] = x2 - (0.00954 * x3)
+        g[2] = (np.pi * x3 * x3 * x4) + ((4.0 / 3.0) * (np.pi * x3 * x3 * x3)) - 1296000
+        g = np.where(g < 0, -g, 0)
+        f[1] = g[0] + g[1] + g[2]
+
+        return f
+
+
+class RE24():
+    def __init__(self):
+        self.problem_name = 'RE24'
+        self.n_objectives = 2
+        self.n_variables = 2
+        self.n_constraints = 0
+        self.n_original_constraints = 4
+
+        self.ubound = np.zeros(self.n_variables)
+        self.lbound = np.zeros(self.n_variables)
+        self.lbound[0] = 0.5
+        self.lbound[1] = 0.5
+        self.ubound[0] = 4
+        self.ubound[1] = 50
+
+    def evaluate(self, x):
+        f = np.zeros(self.n_objectives)
+        g = np.zeros(self.n_original_constraints)
+
+        x1 = x[0]
+        x2 = x[1]
+
+        # First original objective function
+        f[0] = x1 + (120 * x2)
+
+        E = 700000
+        sigma_b_max = 700
+        tau_max = 450
+        delta_max = 1.5
+        sigma_k = (E * x1 * x1) / 100
+        sigma_b = 4500 / (x1 * x2)
+        tau = 1800 / x2
+        delta = (56.2 * 10000) / (E * x1 * x2 * x2)
+
+        g[0] = 1 - (sigma_b / sigma_b_max)
+        g[1] = 1 - (tau / tau_max)
+        g[2] = 1 - (delta / delta_max)
+        g[3] = 1 - (sigma_b / sigma_k)
+        g = np.where(g < 0, -g, 0)
+        f[1] = g[0] + g[1] + g[2] + g[3]
+
+        return f
+
+
+class RE25():
+    def __init__(self):
+        self.problem_name = 'RE25'
+        self.n_objectives = 2
+        self.n_variables = 3
+        self.n_constraints = 0
+        self.n_original_constraints = 6
+
+        self.ubound = np.zeros(self.n_variables)
+        self.lbound = np.zeros(self.n_variables)
+        self.lbound[0] = 1
+        self.lbound[1] = 0.6
+        self.lbound[2] = 0.09
+        self.ubound[0] = 70
+        self.ubound[1] = 3
+        self.ubound[2] = 0.5
+
+        self.feasible_vals = np.array(
+            [0.009, 0.0095, 0.0104, 0.0118, 0.0128, 0.0132, 0.014, 0.015, 0.0162, 0.0173, 0.018, 0.02, 0.023, 0.025,
+             0.028, 0.032, 0.035, 0.041, 0.047, 0.054, 0.063, 0.072, 0.08, 0.092, 0.105, 0.12, 0.135, 0.148, 0.162,
+             0.177, 0.192, 0.207, 0.225, 0.244, 0.263, 0.283, 0.307, 0.331, 0.362, 0.394, 0.4375, 0.5])
+
+    def evaluate(self, x):
+        f = np.zeros(self.n_objectives)
+        g = np.zeros(self.n_original_constraints)
+
+        x1 = np.round(x[0])
+        x2 = x[1]
+        # Reference: getNearestValue_sample2.py (https://gist.github.com/icchi-h/1d0bb1c52ebfdd31f14b3e811328390a)
+        idx = np.abs(np.asarray(self.feasible_vals) - x[2]).argmin()
+        x3 = self.feasible_vals[idx]
+
+        # first original objective function
+        f[0] = (np.pi * np.pi * x2 * x3 * x3 * (x1 + 2)) / 4.0
+
+        # constraint functions
+        Cf = ((4.0 * (x2 / x3) - 1) / (4.0 * (x2 / x3) - 4)) + (0.615 * x3 / x2)
+        Fmax = 1000.0
+        S = 189000.0
+        G = 11.5 * 1e+6
+        K = (G * x3 * x3 * x3 * x3) / (8 * x1 * x2 * x2 * x2)
+        lmax = 14.0
+        lf = (Fmax / K) + 1.05 * (x1 + 2) * x3
+        dmin = 0.2
+        Dmax = 3
+        Fp = 300.0
+        sigmaP = Fp / K
+        sigmaPM = 6
+        sigmaW = 1.25
+
+        g[0] = -((8 * Cf * Fmax * x2) / (np.pi * x3 * x3 * x3)) + S
+        g[1] = -lf + lmax
+        g[2] = -3 + (x2 / x3)
+        g[3] = -sigmaP + sigmaPM
+        g[4] = -sigmaP - ((Fmax - Fp) / K) - 1.05 * (x1 + 2) * x3 + lf
+        g[5] = sigmaW - ((Fmax - Fp) / K)
+
+        g = np.where(g < 0, -g, 0)
+        f[1] = g[0] + g[1] + g[2] + g[3] + g[4] + g[5]
+
+        return f
+
+
+class RE31():
+    def __init__(self):
+        self.problem_name = 'RE31'
+        self.n_objectives = 3
+        self.n_variables = 3
+        self.n_constraints = 0
+        self.n_original_constraints = 3
+
+        self.ubound = np.zeros(self.n_variables)
+        self.lbound = np.zeros(self.n_variables)
+        self.lbound[0] = 0.00001
+        self.lbound[1] = 0.00001
+        self.lbound[2] = 1.0
+        self.ubound[0] = 100.0
+        self.ubound[1] = 100.0
+        self.ubound[2] = 3.0
+
+    def evaluate(self, x):
+        f = np.zeros(self.n_objectives)
+        g = np.zeros(self.n_original_constraints)
+
+        x1 = x[0]
+        x2 = x[1]
+        x3 = x[2]
+
+        # First original objective function
+        f[0] = x1 * np.sqrt(16.0 + (x3 * x3)) + x2 * np.sqrt(1.0 + x3 * x3)
+        # Second original objective function
+        f[1] = (20.0 * np.sqrt(16.0 + (x3 * x3))) / (x1 * x3)
+
+        # Constraint functions
+        g[0] = 0.1 - f[0]
+        g[1] = 100000.0 - f[1]
+        g[2] = 100000 - ((80.0 * np.sqrt(1.0 + x3 * x3)) / (x3 * x2))
+        g = np.where(g < 0, -g, 0)
+        f[2] = g[0] + g[1] + g[2]
+
+        return f
+
+
+class RE32():
+    def __init__(self):
+        self.problem_name = 'RE32'
+        self.n_objectives = 3
+        self.n_variables = 4
+        self.n_constraints = 0
+        self.n_original_constraints = 4
+
+        self.ubound = np.zeros(self.n_variables)
+        self.lbound = np.zeros(self.n_variables)
+        self.lbound[0] = 0.125
+        self.lbound[1] = 0.1
+        self.lbound[2] = 0.1
+        self.lbound[3] = 0.125
+        self.ubound[0] = 5.0
+        self.ubound[1] = 10.0
+        self.ubound[2] = 10.0
+        self.ubound[3] = 5.0
+
+    def evaluate(self, x):
+        f = np.zeros(self.n_objectives)
+        g = np.zeros(self.n_original_constraints)
+
+        x1 = x[0]
+        x2 = x[1]
+        x3 = x[2]
+        x4 = x[3]
+
+        P = 6000
+        L = 14
+        E = 30 * 1e6
+
+        # // deltaMax = 0.25
+        G = 12 * 1e6
+        tauMax = 13600
+        sigmaMax = 30000
+
+        # First original objective function
+        f[0] = (1.10471 * x1 * x1 * x2) + (0.04811 * x3 * x4) * (14.0 + x2)
+        # Second original objective function
+        f[1] = (4 * P * L * L * L) / (E * x4 * x3 * x3 * x3)
+
+        # Constraint functions
+        M = P * (L + (x2 / 2))
+        tmpVar = ((x2 * x2) / 4.0) + np.power((x1 + x3) / 2.0, 2)
+        R = np.sqrt(tmpVar)
+        tmpVar = ((x2 * x2) / 12.0) + np.power((x1 + x3) / 2.0, 2)
+        J = 2 * np.sqrt(2) * x1 * x2 * tmpVar
+
+        tauDashDash = (M * R) / J
+        tauDash = P / (np.sqrt(2) * x1 * x2)
+        tmpVar = tauDash * tauDash + ((2 * tauDash * tauDashDash * x2) / (2 * R)) + (tauDashDash * tauDashDash)
+        tau = np.sqrt(tmpVar)
+        sigma = (6 * P * L) / (x4 * x3 * x3)
+        tmpVar = 4.013 * E * np.sqrt((x3 * x3 * x4 * x4 * x4 * x4 * x4 * x4) / 36.0) / (L * L)
+        tmpVar2 = (x3 / (2 * L)) * np.sqrt(E / (4 * G))
+        PC = tmpVar * (1 - tmpVar2)
+
+        g[0] = tauMax - tau
+        g[1] = sigmaMax - sigma
+        g[2] = x4 - x1
+        g[3] = PC - P
+        g = np.where(g < 0, -g, 0)
+        f[2] = g[0] + g[1] + g[2] + g[3]
+
+        return f
+
+
+class RE33():
+    def __init__(self):
+        self.problem_name = 'RE33'
+        self.n_objectives = 3
+        self.n_variables = 4
+        self.n_constraints = 0
+        self.n_original_constraints = 4
+
+        self.ubound = np.zeros(self.n_variables)
+        self.lbound = np.zeros(self.n_variables)
+        self.lbound[0] = 55
+        self.lbound[1] = 75
+        self.lbound[2] = 1000
+        self.lbound[3] = 11
+        self.ubound[0] = 80
+        self.ubound[1] = 110
+        self.ubound[2] = 3000
+        self.ubound[3] = 20
+
+    def evaluate(self, x):
+        f = np.zeros(self.n_objectives)
+        g = np.zeros(self.n_original_constraints)
+
+        x1 = x[0]
+        x2 = x[1]
+        x3 = x[2]
+        x4 = x[3]
+
+        # First original objective function
+        f[0] = 4.9 * 1e-5 * (x2 * x2 - x1 * x1) * (x4 - 1.0)
+        # Second original objective function
+        f[1] = ((9.82 * 1e6) * (x2 * x2 - x1 * x1)) / (x3 * x4 * (x2 * x2 * x2 - x1 * x1 * x1))
+
+        # Reformulated objective functions
+        g[0] = (x2 - x1) - 20.0
+        g[1] = 0.4 - (x3 / (3.14 * (x2 * x2 - x1 * x1)))
+        g[2] = 1.0 - (2.22 * 1e-3 * x3 * (x2 * x2 * x2 - x1 * x1 * x1)) / np.power((x2 * x2 - x1 * x1), 2)
+        g[3] = (2.66 * 1e-2 * x3 * x4 * (x2 * x2 * x2 - x1 * x1 * x1)) / (x2 * x2 - x1 * x1) - 900.0
+        g = np.where(g < 0, -g, 0)
+        f[2] = g[0] + g[1] + g[2] + g[3]
+
+        return f
+
+
+class RE34():
+    def __init__(self):
+        self.problem_name = 'RE34'
+        self.n_objectives = 3
+        self.n_variables = 5
+        self.n_constraints = 0
+        self.n_original_constraints = 0
+
+        self.lbound = np.full(self.n_variables, 1)
+        self.ubound = np.full(self.n_variables, 3)
+
+    def evaluate(self, x):
+        f = np.zeros(self.n_objectives)
+        g = np.zeros(self.n_original_constraints)
+
+        x1 = x[0]
+        x2 = x[1]
+        x3 = x[2]
+        x4 = x[3]
+        x5 = x[4]
+
+        f[0] = 1640.2823 + (2.3573285 * x1) + (2.3220035 * x2) + (4.5688768 * x3) + (7.7213633 * x4) + (4.4559504 * x5)
+        f[1] = 6.5856 + (1.15 * x1) - (1.0427 * x2) + (0.9738 * x3) + (0.8364 * x4) - (0.3695 * x1 * x4) + (
+                    0.0861 * x1 * x5) + (0.3628 * x2 * x4) - (0.1106 * x1 * x1) - (0.3437 * x3 * x3) + (
+                           0.1764 * x4 * x4)
+        f[2] = -0.0551 + (0.0181 * x1) + (0.1024 * x2) + (0.0421 * x3) - (0.0073 * x1 * x2) + (0.024 * x2 * x3) - (
+                    0.0118 * x2 * x4) - (0.0204 * x3 * x4) - (0.008 * x3 * x5) - (0.0241 * x2 * x2) + (0.0109 * x4 * x4)
+
+        return f
+
+
+class RE35():
+    def __init__(self):
+        self.problem_name = 'RE35'
+        self.n_objectives = 3
+        self.n_variables = 7
+        self.n_constraints = 0
+        self.n_original_constraints = 11
+
+        self.lbound = np.zeros(self.n_variables)
+        self.ubound = np.zeros(self.n_variables)
+        self.lbound[0] = 2.6
+        self.lbound[1] = 0.7
+        self.lbound[2] = 17
+        self.lbound[3] = 7.3
+        self.lbound[4] = 7.3
+        self.lbound[5] = 2.9
+        self.lbound[6] = 5.0
+        self.ubound[0] = 3.6
+        self.ubound[1] = 0.8
+        self.ubound[2] = 28
+        self.ubound[3] = 8.3
+        self.ubound[4] = 8.3
+        self.ubound[5] = 3.9
+        self.ubound[6] = 5.5
+
+    def evaluate(self, x):
+        f = np.zeros(self.n_objectives)
+        g = np.zeros(self.n_original_constraints)
+
+        x1 = x[0]
+        x2 = x[1]
+        x3 = np.round(x[2])
+        x4 = x[3]
+        x5 = x[4]
+        x6 = x[5]
+        x7 = x[6]
+
+        # First original objective function (weight)
+        f[0] = 0.7854 * x1 * (x2 * x2) * (((10.0 * x3 * x3) / 3.0) + (14.933 * x3) - 43.0934) - 1.508 * x1 * (
+                    x6 * x6 + x7 * x7) + 7.477 * (x6 * x6 * x6 + x7 * x7 * x7) + 0.7854 * (x4 * x6 * x6 + x5 * x7 * x7)
+
+        # Second original objective function (stress)
+        tmpVar = np.power((745.0 * x4) / (x2 * x3), 2.0) + 1.69 * 1e7
+        f[1] = np.sqrt(tmpVar) / (0.1 * x6 * x6 * x6)
+
+        # Constraint functions
+        g[0] = -(1.0 / (x1 * x2 * x2 * x3)) + 1.0 / 27.0
+        g[1] = -(1.0 / (x1 * x2 * x2 * x3 * x3)) + 1.0 / 397.5
+        g[2] = -(x4 * x4 * x4) / (x2 * x3 * x6 * x6 * x6 * x6) + 1.0 / 1.93
+        g[3] = -(x5 * x5 * x5) / (x2 * x3 * x7 * x7 * x7 * x7) + 1.0 / 1.93
+        g[4] = -(x2 * x3) + 40.0
+        g[5] = -(x1 / x2) + 12.0
+        g[6] = -5.0 + (x1 / x2)
+        g[7] = -1.9 + x4 - 1.5 * x6
+        g[8] = -1.9 + x5 - 1.1 * x7
+        g[9] = -f[1] + 1300.0
+        tmpVar = np.power((745.0 * x5) / (x2 * x3), 2.0) + 1.575 * 1e8
+        g[10] = -np.sqrt(tmpVar) / (0.1 * x7 * x7 * x7) + 1100.0
+        g = np.where(g < 0, -g, 0)
+        f[2] = g[0] + g[1] + g[2] + g[3] + g[4] + g[5] + g[6] + g[7] + g[8] + g[9] + g[10]
+
+        return f
+
+
+class RE36():
+    def __init__(self):
+        self.problem_name = 'RE36'
+        self.n_objectives = 3
+        self.n_variables = 4
+        self.n_constraints = 0
+        self.n_original_constraints = 1
+
+        self.lbound = np.full(self.n_variables, 12)
+        self.ubound = np.full(self.n_variables, 60)
+
+    def evaluate(self, x):
+        f = np.zeros(self.n_objectives)
+        g = np.zeros(self.n_original_constraints)
+
+        # all the four variables must be inverger values
+        x1 = np.round(x[0])
+        x2 = np.round(x[1])
+        x3 = np.round(x[2])
+        x4 = np.round(x[3])
+
+        # First original objective function
+        f[0] = np.abs(6.931 - ((x3 / x1) * (x4 / x2)))
+        # Second original objective function (the maximum value among the four variables)
+        l = [x1, x2, x3, x4]
+        f[1] = max(l)
+
+        g[0] = 0.5 - (f[0] / 6.931)
+        g = np.where(g < 0, -g, 0)
+        f[2] = g[0]
+
+        return f
+
+
+class RE37():
+    def __init__(self):
+        self.problem_name = 'RE37'
+        self.n_objectives = 3
+        self.n_variables = 4
+        self.n_constraints = 0
+        self.n_original_constraints = 0
+
+        self.lbound = np.full(self.n_variables, 0)
+        self.ubound = np.full(self.n_variables, 1)
+
+    def evaluate(self, x):
+        f = np.zeros(self.n_objectives)
+
+        xAlpha = x[0]
+        xHA = x[1]
+        xOA = x[2]
+        xOPTT = x[3]
+
+        # f1 (TF_max)
+        f[0] = 0.692 + (0.477 * xAlpha) - (0.687 * xHA) - (0.080 * xOA) - (0.0650 * xOPTT) - (
+                    0.167 * xAlpha * xAlpha) - (0.0129 * xHA * xAlpha) + (0.0796 * xHA * xHA) - (
+                           0.0634 * xOA * xAlpha) - (0.0257 * xOA * xHA) + (0.0877 * xOA * xOA) - (
+                           0.0521 * xOPTT * xAlpha) + (0.00156 * xOPTT * xHA) + (0.00198 * xOPTT * xOA) + (
+                           0.0184 * xOPTT * xOPTT)
+        # f2 (X_cc)
+        f[1] = 0.153 - (0.322 * xAlpha) + (0.396 * xHA) + (0.424 * xOA) + (0.0226 * xOPTT) + (
+                    0.175 * xAlpha * xAlpha) + (0.0185 * xHA * xAlpha) - (0.0701 * xHA * xHA) - (
+                           0.251 * xOA * xAlpha) + (0.179 * xOA * xHA) + (0.0150 * xOA * xOA) + (
+                           0.0134 * xOPTT * xAlpha) + (0.0296 * xOPTT * xHA) + (0.0752 * xOPTT * xOA) + (
+                           0.0192 * xOPTT * xOPTT)
+        # f3 (TT_max)
+        f[2] = 0.370 - (0.205 * xAlpha) + (0.0307 * xHA) + (0.108 * xOA) + (1.019 * xOPTT) - (
+                    0.135 * xAlpha * xAlpha) + (0.0141 * xHA * xAlpha) + (0.0998 * xHA * xHA) + (
+                           0.208 * xOA * xAlpha) - (0.0301 * xOA * xHA) - (0.226 * xOA * xOA) + (
+                           0.353 * xOPTT * xAlpha) - (0.0497 * xOPTT * xOA) - (0.423 * xOPTT * xOPTT) + (
+                           0.202 * xHA * xAlpha * xAlpha) - (0.281 * xOA * xAlpha * xAlpha) - (
+                           0.342 * xHA * xHA * xAlpha) - (0.245 * xHA * xHA * xOA) + (0.281 * xOA * xOA * xHA) - (
+                           0.184 * xOPTT * xOPTT * xAlpha) - (0.281 * xHA * xAlpha * xOA)
+
+        return f
+
+
+class RE41():
+    def __init__(self):
+        self.problem_name = 'RE41'
+        self.n_objectives = 4
+        self.n_variables = 7
+        self.n_constraints = 0
+        self.n_original_constraints = 10
+
+        self.lbound = np.zeros(self.n_variables)
+        self.ubound = np.zeros(self.n_variables)
+        self.lbound[0] = 0.5
+        self.lbound[1] = 0.45
+        self.lbound[2] = 0.5
+        self.lbound[3] = 0.5
+        self.lbound[4] = 0.875
+        self.lbound[5] = 0.4
+        self.lbound[6] = 0.4
+        self.ubound[0] = 1.5
+        self.ubound[1] = 1.35
+        self.ubound[2] = 1.5
+        self.ubound[3] = 1.5
+        self.ubound[4] = 2.625
+        self.ubound[5] = 1.2
+        self.ubound[6] = 1.2
+
+    def evaluate(self, x):
+        f = np.zeros(self.n_objectives)
+        g = np.zeros(self.n_original_constraints)
+
+        x1 = x[0]
+        x2 = x[1]
+        x3 = x[2]
+        x4 = x[3]
+        x5 = x[4]
+        x6 = x[5]
+        x7 = x[6]
+
+        # First original objective function
+        f[0] = 1.98 + 4.9 * x1 + 6.67 * x2 + 6.98 * x3 + 4.01 * x4 + 1.78 * x5 + 0.00001 * x6 + 2.73 * x7
+        # Second original objective function
+        f[1] = 4.72 - 0.5 * x4 - 0.19 * x2 * x3
+        # Third original objective function
+        Vmbp = 10.58 - 0.674 * x1 * x2 - 0.67275 * x2
+        Vfd = 16.45 - 0.489 * x3 * x7 - 0.843 * x5 * x6
+        f[2] = 0.5 * (Vmbp + Vfd)
+
+        # Constraint functions
+        g[0] = 1 - (1.16 - 0.3717 * x2 * x4 - 0.0092928 * x3)
+        g[1] = 0.32 - (0.261 - 0.0159 * x1 * x2 - 0.06486 * x1 - 0.019 * x2 * x7 + 0.0144 * x3 * x5 + 0.0154464 * x6)
+        g[2] = 0.32 - (
+                    0.214 + 0.00817 * x5 - 0.045195 * x1 - 0.0135168 * x1 + 0.03099 * x2 * x6 - 0.018 * x2 * x7 + 0.007176 * x3 + 0.023232 * x3 - 0.00364 * x5 * x6 - 0.018 * x2 * x2)
+        g[3] = 0.32 - (0.74 - 0.61 * x2 - 0.031296 * x3 - 0.031872 * x7 + 0.227 * x2 * x2)
+        g[4] = 32 - (28.98 + 3.818 * x3 - 4.2 * x1 * x2 + 1.27296 * x6 - 2.68065 * x7)
+        g[5] = 32 - (33.86 + 2.95 * x3 - 5.057 * x1 * x2 - 3.795 * x2 - 3.4431 * x7 + 1.45728)
+        g[6] = 32 - (46.36 - 9.9 * x2 - 4.4505 * x1)
+        g[7] = 4 - f[1]
+        g[8] = 9.9 - Vmbp
+        g[9] = 15.7 - Vfd
+
+        g = np.where(g < 0, -g, 0)
+        f[3] = g[0] + g[1] + g[2] + g[3] + g[4] + g[5] + g[6] + g[7] + g[8] + g[9]
+
+        return f
+
+
+class RE42():
+    def __init__(self):
+        self.problem_name = 'RE42'
+        self.n_objectives = 4
+        self.n_variables = 6
+        self.n_constraints = 0
+        self.n_original_constraints = 9
+
+        self.lbound = np.zeros(self.n_variables)
+        self.ubound = np.zeros(self.n_variables)
+        self.lbound[0] = 150.0
+        self.lbound[1] = 20.0
+        self.lbound[2] = 13.0
+        self.lbound[3] = 10.0
+        self.lbound[4] = 14.0
+        self.lbound[5] = 0.63
+        self.ubound[0] = 274.32
+        self.ubound[1] = 32.31
+        self.ubound[2] = 25.0
+        self.ubound[3] = 11.71
+        self.ubound[4] = 18.0
+        self.ubound[5] = 0.75
+
+    def evaluate(self, x):
+        f = np.zeros(self.n_objectives)
+        # NOT g
+        constraintFuncs = np.zeros(self.n_original_constraints)
+
+        x_L = x[0]
+        x_B = x[1]
+        x_D = x[2]
+        x_T = x[3]
+        x_Vk = x[4]
+        x_CB = x[5]
+
+        displacement = 1.025 * x_L * x_B * x_T * x_CB
+        V = 0.5144 * x_Vk
+        g = 9.8065
+        Fn = V / np.power(g * x_L, 0.5)
+        a = (4977.06 * x_CB * x_CB) - (8105.61 * x_CB) + 4456.51
+        b = (-10847.2 * x_CB * x_CB) + (12817.0 * x_CB) - 6960.32
+
+        power = (np.power(displacement, 2.0 / 3.0) * np.power(x_Vk, 3.0)) / (a + (b * Fn))
+        outfit_weight = 1.0 * np.power(x_L, 0.8) * np.power(x_B, 0.6) * np.power(x_D, 0.3) * np.power(x_CB, 0.1)
+        steel_weight = 0.034 * np.power(x_L, 1.7) * np.power(x_B, 0.7) * np.power(x_D, 0.4) * np.power(x_CB, 0.5)
+        machinery_weight = 0.17 * np.power(power, 0.9)
+        light_ship_weight = steel_weight + outfit_weight + machinery_weight
+
+        ship_cost = 1.3 * ((2000.0 * np.power(steel_weight, 0.85)) + (3500.0 * outfit_weight) + (
+                    2400.0 * np.power(power, 0.8)))
+        capital_costs = 0.2 * ship_cost
+
+        DWT = displacement - light_ship_weight
+
+        running_costs = 40000.0 * np.power(DWT, 0.3)
+
+        round_trip_miles = 5000.0
+        sea_days = (round_trip_miles / 24.0) * x_Vk
+        handling_rate = 8000.0
+
+        daily_consumption = ((0.19 * power * 24.0) / 1000.0) + 0.2
+        fuel_price = 100.0
+        fuel_cost = 1.05 * daily_consumption * sea_days * fuel_price
+        port_cost = 6.3 * np.power(DWT, 0.8)
+
+        fuel_carried = daily_consumption * (sea_days + 5.0)
+        miscellaneous_DWT = 2.0 * np.power(DWT, 0.5)
+
+        cargo_DWT = DWT - fuel_carried - miscellaneous_DWT
+        port_days = 2.0 * ((cargo_DWT / handling_rate) + 0.5)
+        RTPA = 350.0 / (sea_days + port_days)
+
+        voyage_costs = (fuel_cost + port_cost) * RTPA
+        annual_costs = capital_costs + running_costs + voyage_costs
+        annual_cargo = cargo_DWT * RTPA
+
+        f[0] = annual_costs / annual_cargo
+        f[1] = light_ship_weight
+        # f_2 is dealt as a minimization problem
+        f[2] = -annual_cargo
+
+        # Reformulated objective functions
+        constraintFuncs[0] = (x_L / x_B) - 6.0
+        constraintFuncs[1] = -(x_L / x_D) + 15.0
+        constraintFuncs[2] = -(x_L / x_T) + 19.0
+        constraintFuncs[3] = 0.45 * np.power(DWT, 0.31) - x_T
+        constraintFuncs[4] = 0.7 * x_D + 0.7 - x_T
+        constraintFuncs[5] = 500000.0 - DWT
+        constraintFuncs[6] = DWT - 3000.0
+        constraintFuncs[7] = 0.32 - Fn
+
+        KB = 0.53 * x_T
+        BMT = ((0.085 * x_CB - 0.002) * x_B * x_B) / (x_T * x_CB)
+        KG = 1.0 + 0.52 * x_D
+        constraintFuncs[8] = (KB + BMT - KG) - (0.07 * x_B)
+
+        constraintFuncs = np.where(constraintFuncs < 0, -constraintFuncs, 0)
+        f[3] = constraintFuncs[0] + constraintFuncs[1] + constraintFuncs[2] + constraintFuncs[3] + constraintFuncs[4] + \
+               constraintFuncs[5] + constraintFuncs[6] + constraintFuncs[7] + constraintFuncs[8]
+
+        return f
+
+
+class RE61():
+    def __init__(self):
+        self.problem_name = 'RE61'
+        self.n_objectives = 6
+        self.n_variables = 3
+        self.n_constraints = 0
+        self.n_original_constraints = 7
+
+        self.lbound = np.zeros(self.n_variables)
+        self.ubound = np.zeros(self.n_variables)
+        self.lbound[0] = 0.01
+        self.lbound[1] = 0.01
+        self.lbound[2] = 0.01
+        self.ubound[0] = 0.45
+        self.ubound[1] = 0.10
+        self.ubound[2] = 0.10
+
+    def evaluate(self, x):
+        f = np.zeros(self.n_objectives)
+        g = np.zeros(self.n_original_constraints)
+
+        # First original objective function
+        f[0] = 106780.37 * (x[1] + x[2]) + 61704.67
+        # Second original objective function
+        f[1] = 3000 * x[0]
+        # Third original objective function
+        f[2] = 305700 * 2289 * x[1] / np.power(0.06 * 2289, 0.65)
+        # Fourth original objective function
+        f[3] = 250 * 2289 * np.exp(-39.75 * x[1] + 9.9 * x[2] + 2.74)
+        # Fifth original objective function
+        f[4] = 25 * (1.39 / (x[0] * x[1]) + 4940 * x[2] - 80)
+
+        # Constraint functions
+        g[0] = 1 - (0.00139 / (x[0] * x[1]) + 4.94 * x[2] - 0.08)
+        g[1] = 1 - (0.000306 / (x[0] * x[1]) + 1.082 * x[2] - 0.0986)
+        g[2] = 50000 - (12.307 / (x[0] * x[1]) + 49408.24 * x[2] + 4051.02)
+        g[3] = 16000 - (2.098 / (x[0] * x[1]) + 8046.33 * x[2] - 696.71)
+        g[4] = 10000 - (2.138 / (x[0] * x[1]) + 7883.39 * x[2] - 705.04)
+        g[5] = 2000 - (0.417 * x[0] * x[1] + 1721.26 * x[2] - 136.54)
+        g[6] = 550 - (0.164 / (x[0] * x[1]) + 631.13 * x[2] - 54.48)
+
+        g = np.where(g < 0, -g, 0)
+        f[5] = g[0] + g[1] + g[2] + g[3] + g[4] + g[5] + g[6]
+
+        return f
+
+
+class RE91():
+    def __init__(self):
+        self.problem_name = 'RE91'
+        self.n_objectives = 9
+        self.n_variables = 7
+        self.n_constraints = 0
+        self.n_original_constraints = 0
+
+        self.lbound = np.zeros(self.n_variables)
+        self.ubound = np.zeros(self.n_variables)
+        self.lbound[0] = 0.5
+        self.lbound[1] = 0.45
+        self.lbound[2] = 0.5
+        self.lbound[3] = 0.5
+        self.lbound[4] = 0.875
+        self.lbound[5] = 0.4
+        self.lbound[6] = 0.4
+        self.ubound[0] = 1.5
+        self.ubound[1] = 1.35
+        self.ubound[2] = 1.5
+        self.ubound[3] = 1.5
+        self.ubound[4] = 2.625
+        self.ubound[5] = 1.2
+        self.ubound[6] = 1.2
+
+    def evaluate(self, x):
+        f = np.zeros(self.n_objectives)
+        g = np.zeros(self.n_original_constraints)
+
+        x1 = x[0]
+        x2 = x[1]
+        x3 = x[2]
+        x4 = x[3]
+        x5 = x[4]
+        x6 = x[5]
+        x7 = x[6]
+        # stochastic variables
+        x8 = 0.006 * (np.random.normal(0, 1)) + 0.345
+        x9 = 0.006 * (np.random.normal(0, 1)) + 0.192
+        x10 = 10 * (np.random.normal(0, 1)) + 0.0
+        x11 = 10 * (np.random.normal(0, 1)) + 0.0
+
+        # First function
+        f[0] = 1.98 + 4.9 * x1 + 6.67 * x2 + 6.98 * x3 + 4.01 * x4 + 1.75 * x5 + 0.00001 * x6 + 2.73 * x7
+        # Second function
+        f[1] = max(0.0, (1.16 - 0.3717 * x2 * x4 - 0.00931 * x2 * x10 - 0.484 * x3 * x9 + 0.01343 * x6 * x10) / 1.0)
+        # Third function
+        f[2] = max(0.0, (
+                    0.261 - 0.0159 * x1 * x2 - 0.188 * x1 * x8 - 0.019 * x2 * x7 + 0.0144 * x3 * x5 + 0.87570001 * x5 * x10 + 0.08045 * x6 * x9 + 0.00139 * x8 * x11 + 0.00001575 * x10 * x11) / 0.32)
+        # Fourth function
+        f[3] = max(0.0, (
+                    0.214 + 0.00817 * x5 - 0.131 * x1 * x8 - 0.0704 * x1 * x9 + 0.03099 * x2 * x6 - 0.018 * x2 * x7 + 0.0208 * x3 * x8 + 0.121 * x3 * x9 - 0.00364 * x5 * x6 + 0.0007715 * x5 * x10 - 0.0005354 * x6 * x10 + 0.00121 * x8 * x11 + 0.00184 * x9 * x10 - 0.018 * x2 * x2) / 0.32)
+        # Fifth function
+        f[4] = max(0.0, (
+                    0.74 - 0.61 * x2 - 0.163 * x3 * x8 + 0.001232 * x3 * x10 - 0.166 * x7 * x9 + 0.227 * x2 * x2) / 0.32)
+        # Sixth function
+        tmp = ((
+                           28.98 + 3.818 * x3 - 4.2 * x1 * x2 + 0.0207 * x5 * x10 + 6.63 * x6 * x9 - 7.77 * x7 * x8 + 0.32 * x9 * x10) + (
+                           33.86 + 2.95 * x3 + 0.1792 * x10 - 5.057 * x1 * x2 - 11 * x2 * x8 - 0.0215 * x5 * x10 - 9.98 * x7 * x8 + 22 * x8 * x9) + (
+                           46.36 - 9.9 * x2 - 12.9 * x1 * x8 + 0.1107 * x3 * x10)) / 3
+        f[5] = max(0.0, tmp / 32)
+        # Seventh function
+        f[6] = max(0.0, (
+                    4.72 - 0.5 * x4 - 0.19 * x2 * x3 - 0.0122 * x4 * x10 + 0.009325 * x6 * x10 + 0.000191 * x11 * x11) / 4.0)
+        # EighthEighth function
+        f[7] = max(0.0, (
+                    10.58 - 0.674 * x1 * x2 - 1.95 * x2 * x8 + 0.02054 * x3 * x10 - 0.0198 * x4 * x10 + 0.028 * x6 * x10) / 9.9)
+        # Ninth function
+        f[8] = max(0.0, (
+                    16.45 - 0.489 * x3 * x7 - 0.843 * x5 * x6 + 0.0432 * x9 * x10 - 0.0556 * x9 * x11 - 0.000786 * x11 * x11) / 15.7)
+
+        return f
+
+
+class CRE21():
+    def __init__(self):
+        self.problem_name = 'CRE21'
+        self.n_objectives = 2
+        self.n_variables = 3
+        self.n_constraints = 3
+
+        self.ubound = np.zeros(self.n_variables)
+        self.lbound = np.zeros(self.n_variables)
+        self.lbound[0] = 0.00001
+        self.lbound[1] = 0.00001
+        self.lbound[2] = 1.0
+        self.ubound[0] = 100.0
+        self.ubound[1] = 100.0
+        self.ubound[2] = 3.0
+
+    def evaluate(self, x):
+        f = np.zeros(self.n_objectives)
+        g = np.zeros(self.n_constraints)
+
+        x1 = x[0]
+        x2 = x[1]
+        x3 = x[2]
+
+        # First original objective function
+        f[0] = x1 * np.sqrt(16.0 + (x3 * x3)) + x2 * np.sqrt(1.0 + x3 * x3)
+        # Second original objective function
+        f[1] = (20.0 * np.sqrt(16.0 + (x3 * x3))) / (x1 * x3)
+
+        # Constraint functions
+        g[0] = 0.1 - f[0]
+        g[1] = 100000.0 - f[1]
+        g[2] = 100000 - ((80.0 * np.sqrt(1.0 + x3 * x3)) / (x3 * x2))
+        g = np.where(g < 0, -g, 0)
+
+        return f, g
+
+
+class CRE22():
+    def __init__(self):
+        self.problem_name = 'CRE22'
+        self.n_objectives = 2
+        self.n_variables = 4
+        self.n_constraints = 4
+
+        self.ubound = np.zeros(self.n_variables)
+        self.lbound = np.zeros(self.n_variables)
+        self.lbound[0] = 0.125
+        self.lbound[1] = 0.1
+        self.lbound[2] = 0.1
+        self.lbound[3] = 0.125
+        self.ubound[0] = 5.0
+        self.ubound[1] = 10.0
+        self.ubound[2] = 10.0
+        self.ubound[3] = 5.0
+
+    def evaluate(self, x):
+        f = np.zeros(self.n_objectives)
+        g = np.zeros(self.n_constraints)
+
+        x1 = x[0]
+        x2 = x[1]
+        x3 = x[2]
+        x4 = x[3]
+
+        P = 6000
+        L = 14
+        E = 30 * 1e6
+
+        # // deltaMax = 0.25
+        G = 12 * 1e6
+        tauMax = 13600
+        sigmaMax = 30000
+
+        # First original objective function
+        f[0] = (1.10471 * x1 * x1 * x2) + (0.04811 * x3 * x4) * (14.0 + x2)
+        # Second original objective function
+        f[1] = (4 * P * L * L * L) / (E * x4 * x3 * x3 * x3)
+
+        # Constraint functions
+        M = P * (L + (x2 / 2))
+        tmpVar = ((x2 * x2) / 4.0) + np.power((x1 + x3) / 2.0, 2)
+        R = np.sqrt(tmpVar)
+        tmpVar = ((x2 * x2) / 12.0) + np.power((x1 + x3) / 2.0, 2)
+        J = 2 * np.sqrt(2) * x1 * x2 * tmpVar
+
+        tauDashDash = (M * R) / J
+        tauDash = P / (np.sqrt(2) * x1 * x2)
+        tmpVar = tauDash * tauDash + ((2 * tauDash * tauDashDash * x2) / (2 * R)) + (tauDashDash * tauDashDash)
+        tau = np.sqrt(tmpVar)
+        sigma = (6 * P * L) / (x4 * x3 * x3)
+        tmpVar = 4.013 * E * np.sqrt((x3 * x3 * x4 * x4 * x4 * x4 * x4 * x4) / 36.0) / (L * L)
+        tmpVar2 = (x3 / (2 * L)) * np.sqrt(E / (4 * G))
+        PC = tmpVar * (1 - tmpVar2)
+
+        g[0] = tauMax - tau
+        g[1] = sigmaMax - sigma
+        g[2] = x4 - x1
+        g[3] = PC - P
+        g = np.where(g < 0, -g, 0)
+
+        return f, g
+
+
+class CRE23():
+    def __init__(self):
+        self.problem_name = 'CRE23'
+        self.n_objectives = 2
+        self.n_variables = 4
+        self.n_constraints = 4
+
+        self.ubound = np.zeros(self.n_variables)
+        self.lbound = np.zeros(self.n_variables)
+        self.lbound[0] = 55
+        self.lbound[1] = 75
+        self.lbound[2] = 1000
+        self.lbound[3] = 11
+        self.ubound[0] = 80
+        self.ubound[1] = 110
+        self.ubound[2] = 3000
+        self.ubound[3] = 20
+
+    def evaluate(self, x):
+        f = np.zeros(self.n_objectives)
+        g = np.zeros(self.n_constraints)
+
+        x1 = x[0]
+        x2 = x[1]
+        x3 = x[2]
+        x4 = x[3]
+
+        # First original objective function
+        f[0] = 4.9 * 1e-5 * (x2 * x2 - x1 * x1) * (x4 - 1.0)
+        # Second original objective function
+        f[1] = ((9.82 * 1e6) * (x2 * x2 - x1 * x1)) / (x3 * x4 * (x2 * x2 * x2 - x1 * x1 * x1))
+
+        # Reformulated objective functions
+        g[0] = (x2 - x1) - 20.0
+        g[1] = 0.4 - (x3 / (3.14 * (x2 * x2 - x1 * x1)))
+        g[2] = 1.0 - (2.22 * 1e-3 * x3 * (x2 * x2 * x2 - x1 * x1 * x1)) / np.power((x2 * x2 - x1 * x1), 2)
+        g[3] = (2.66 * 1e-2 * x3 * x4 * (x2 * x2 * x2 - x1 * x1 * x1)) / (x2 * x2 - x1 * x1) - 900.0
+        g = np.where(g < 0, -g, 0)
+
+        return f, g
+
+
+class CRE24():
+    def __init__(self):
+        self.problem_name = 'CRE24'
+        self.n_objectives = 2
+        self.n_variables = 7
+        self.n_constraints = 11
+
+        self.lbound = np.zeros(self.n_variables)
+        self.ubound = np.zeros(self.n_variables)
+
+        self.lbound[0] = 2.6
+        self.lbound[1] = 0.7
+        self.lbound[2] = 17
+        self.lbound[3] = 7.3
+        self.lbound[4] = 7.3
+        self.lbound[5] = 2.9
+        self.lbound[6] = 5.0
+        self.ubound[0] = 3.6
+        self.ubound[1] = 0.8
+        self.ubound[2] = 28
+        self.ubound[3] = 8.3
+        self.ubound[4] = 8.3
+        self.ubound[5] = 3.9
+        self.ubound[6] = 5.5
+
+    def evaluate(self, x):
+        f = np.zeros(self.n_objectives)
+        g = np.zeros(self.n_constraints)
+
+        x1 = x[0]
+        x2 = x[1]
+        x3 = np.round(x[2])
+        x4 = x[3]
+        x5 = x[4]
+        x6 = x[5]
+        x7 = x[6]
+
+        # First original objective function (weight)
+        f[0] = 0.7854 * x1 * (x2 * x2) * (((10.0 * x3 * x3) / 3.0) + (14.933 * x3) - 43.0934) - 1.508 * x1 * (
+                    x6 * x6 + x7 * x7) + 7.477 * (x6 * x6 * x6 + x7 * x7 * x7) + 0.7854 * (x4 * x6 * x6 + x5 * x7 * x7)
+
+        # Second original objective function (stress)
+        tmpVar = np.power((745.0 * x4) / (x2 * x3), 2.0) + 1.69 * 1e7
+        f[1] = np.sqrt(tmpVar) / (0.1 * x6 * x6 * x6)
+
+        # Constraint functions
+        g[0] = -(1.0 / (x1 * x2 * x2 * x3)) + 1.0 / 27.0
+        g[1] = -(1.0 / (x1 * x2 * x2 * x3 * x3)) + 1.0 / 397.5
+        g[2] = -(x4 * x4 * x4) / (x2 * x3 * x6 * x6 * x6 * x6) + 1.0 / 1.93
+        g[3] = -(x5 * x5 * x5) / (x2 * x3 * x7 * x7 * x7 * x7) + 1.0 / 1.93
+        g[4] = -(x2 * x3) + 40.0
+        g[5] = -(x1 / x2) + 12.0
+        g[6] = -5.0 + (x1 / x2)
+        g[7] = -1.9 + x4 - 1.5 * x6
+        g[8] = -1.9 + x5 - 1.1 * x7
+        g[9] = -f[1] + 1300.0
+        tmpVar = np.power((745.0 * x5) / (x2 * x3), 2.0) + 1.575 * 1e8
+        g[10] = -np.sqrt(tmpVar) / (0.1 * x7 * x7 * x7) + 1100.0
+        g = np.where(g < 0, -g, 0)
+
+        return f, g
+
+
+class CRE25():
+    def __init__(self):
+        self.problem_name = 'CRE25'
+        self.n_objectives = 2
+        self.n_variables = 4
+        self.n_constraints = 1
+
+        self.lbound = np.full(self.n_variables, 12)
+        self.ubound = np.full(self.n_variables, 60)
+
+    def evaluate(self, x):
+        f = np.zeros(self.n_objectives)
+        g = np.zeros(self.n_constraints)
+
+        # all the four variables must be inverger values
+        x1 = np.round(x[0])
+        x2 = np.round(x[1])
+        x3 = np.round(x[2])
+        x4 = np.round(x[3])
+
+        # First original objective function
+        f[0] = np.abs(6.931 - ((x3 / x1) * (x4 / x2)))
+        # Second original objective function (the maximum value among the four variables)
+        l = [x1, x2, x3, x4]
+        f[1] = max(l)
+
+        g[0] = 0.5 - (f[0] / 6.931)
+        g = np.where(g < 0, -g, 0)
+
+        return f, g
+
+
+class CRE31():
+    def __init__(self):
+        self.problem_name = 'CRE31'
+        self.n_objectives = 3
+        self.n_variables = 7
+        self.n_constraints = 10
+
+        self.lbound = np.zeros(self.n_variables)
+        self.ubound = np.zeros(self.n_variables)
+        self.lbound[0] = 0.5
+        self.lbound[1] = 0.45
+        self.lbound[2] = 0.5
+        self.lbound[3] = 0.5
+        self.lbound[4] = 0.875
+        self.lbound[5] = 0.4
+        self.lbound[6] = 0.4
+        self.ubound[0] = 1.5
+        self.ubound[1] = 1.35
+        self.ubound[2] = 1.5
+        self.ubound[3] = 1.5
+        self.ubound[4] = 2.625
+        self.ubound[5] = 1.2
+        self.ubound[6] = 1.2
+
+    def evaluate(self, x):
+        f = np.zeros(self.n_objectives)
+        g = np.zeros(self.n_constraints)
+
+        x1 = x[0]
+        x2 = x[1]
+        x3 = x[2]
+        x4 = x[3]
+        x5 = x[4]
+        x6 = x[5]
+        x7 = x[6]
+
+        # First original objective function
+        f[0] = 1.98 + 4.9 * x1 + 6.67 * x2 + 6.98 * x3 + 4.01 * x4 + 1.78 * x5 + 0.00001 * x6 + 2.73 * x7
+        # Second original objective function
+        f[1] = 4.72 - 0.5 * x4 - 0.19 * x2 * x3
+        # Third original objective function
+        Vmbp = 10.58 - 0.674 * x1 * x2 - 0.67275 * x2
+        Vfd = 16.45 - 0.489 * x3 * x7 - 0.843 * x5 * x6
+        f[2] = 0.5 * (Vmbp + Vfd)
+
+        # Constraint functions
+        g[0] = 1 - (1.16 - 0.3717 * x2 * x4 - 0.0092928 * x3)
+        g[1] = 0.32 - (0.261 - 0.0159 * x1 * x2 - 0.06486 * x1 - 0.019 * x2 * x7 + 0.0144 * x3 * x5 + 0.0154464 * x6)
+        g[2] = 0.32 - (
+                    0.214 + 0.00817 * x5 - 0.045195 * x1 - 0.0135168 * x1 + 0.03099 * x2 * x6 - 0.018 * x2 * x7 + 0.007176 * x3 + 0.023232 * x3 - 0.00364 * x5 * x6 - 0.018 * x2 * x2)
+        g[3] = 0.32 - (0.74 - 0.61 * x2 - 0.031296 * x3 - 0.031872 * x7 + 0.227 * x2 * x2)
+        g[4] = 32 - (28.98 + 3.818 * x3 - 4.2 * x1 * x2 + 1.27296 * x6 - 2.68065 * x7)
+        g[5] = 32 - (33.86 + 2.95 * x3 - 5.057 * x1 * x2 - 3.795 * x2 - 3.4431 * x7 + 1.45728)
+        g[6] = 32 - (46.36 - 9.9 * x2 - 4.4505 * x1)
+        g[7] = 4 - f[1]
+        g[8] = 9.9 - Vmbp
+        g[9] = 15.7 - Vfd
+        g = np.where(g < 0, -g, 0)
+
+        return f, g
+
+
+class CRE32():
+    def __init__(self):
+        self.problem_name = 'CRE32'
+        self.n_objectives = 3
+        self.n_variables = 6
+        self.n_constraints = 9
+
+        self.lbound = np.zeros(self.n_variables)
+        self.ubound = np.zeros(self.n_variables)
+        self.lbound[0] = 150.0
+        self.lbound[1] = 20.0
+        self.lbound[2] = 13.0
+        self.lbound[3] = 10.0
+        self.lbound[4] = 14.0
+        self.lbound[5] = 0.63
+        self.ubound[0] = 274.32
+        self.ubound[1] = 32.31
+        self.ubound[2] = 25.0
+        self.ubound[3] = 11.71
+        self.ubound[4] = 18.0
+        self.ubound[5] = 0.75
+
+    def evaluate(self, x):
+        f = np.zeros(self.n_objectives)
+        # NOT g
+        constraintFuncs = np.zeros(self.n_constraints)
+
+        x_L = x[0]
+        x_B = x[1]
+        x_D = x[2]
+        x_T = x[3]
+        x_Vk = x[4]
+        x_CB = x[5]
+
+        displacement = 1.025 * x_L * x_B * x_T * x_CB
+        V = 0.5144 * x_Vk
+        g = 9.8065
+        Fn = V / np.power(g * x_L, 0.5)
+        a = (4977.06 * x_CB * x_CB) - (8105.61 * x_CB) + 4456.51
+        b = (-10847.2 * x_CB * x_CB) + (12817.0 * x_CB) - 6960.32
+
+        power = (np.power(displacement, 2.0 / 3.0) * np.power(x_Vk, 3.0)) / (a + (b * Fn))
+        outfit_weight = 1.0 * np.power(x_L, 0.8) * np.power(x_B, 0.6) * np.power(x_D, 0.3) * np.power(x_CB, 0.1)
+        steel_weight = 0.034 * np.power(x_L, 1.7) * np.power(x_B, 0.7) * np.power(x_D, 0.4) * np.power(x_CB, 0.5)
+        machinery_weight = 0.17 * np.power(power, 0.9)
+        light_ship_weight = steel_weight + outfit_weight + machinery_weight
+
+        ship_cost = 1.3 * ((2000.0 * np.power(steel_weight, 0.85)) + (3500.0 * outfit_weight) + (
+                    2400.0 * np.power(power, 0.8)))
+        capital_costs = 0.2 * ship_cost
+
+        DWT = displacement - light_ship_weight
+
+        running_costs = 40000.0 * np.power(DWT, 0.3)
+
+        round_trip_miles = 5000.0
+        sea_days = (round_trip_miles / 24.0) * x_Vk
+        handling_rate = 8000.0
+
+        daily_consumption = ((0.19 * power * 24.0) / 1000.0) + 0.2
+        fuel_price = 100.0
+        fuel_cost = 1.05 * daily_consumption * sea_days * fuel_price
+        port_cost = 6.3 * np.power(DWT, 0.8)
+
+        fuel_carried = daily_consumption * (sea_days + 5.0)
+        miscellaneous_DWT = 2.0 * np.power(DWT, 0.5)
+
+        cargo_DWT = DWT - fuel_carried - miscellaneous_DWT
+        port_days = 2.0 * ((cargo_DWT / handling_rate) + 0.5)
+        RTPA = 350.0 / (sea_days + port_days)
+
+        voyage_costs = (fuel_cost + port_cost) * RTPA
+        annual_costs = capital_costs + running_costs + voyage_costs
+        annual_cargo = cargo_DWT * RTPA
+
+        f[0] = annual_costs / annual_cargo
+        f[1] = light_ship_weight
+        # f_2 is dealt as a minimization problem
+        f[2] = -annual_cargo
+
+        # Reformulated objective functions
+        constraintFuncs[0] = (x_L / x_B) - 6.0
+        constraintFuncs[1] = -(x_L / x_D) + 15.0
+        constraintFuncs[2] = -(x_L / x_T) + 19.0
+        constraintFuncs[3] = 0.45 * np.power(DWT, 0.31) - x_T
+        constraintFuncs[4] = 0.7 * x_D + 0.7 - x_T
+        constraintFuncs[5] = 500000.0 - DWT
+        constraintFuncs[6] = DWT - 3000.0
+        constraintFuncs[7] = 0.32 - Fn
+
+        KB = 0.53 * x_T
+        BMT = ((0.085 * x_CB - 0.002) * x_B * x_B) / (x_T * x_CB)
+        KG = 1.0 + 0.52 * x_D
+        constraintFuncs[8] = (KB + BMT - KG) - (0.07 * x_B)
+        constraintFuncs = np.where(constraintFuncs < 0, -constraintFuncs, 0)
+
+        return f, constraintFuncs
+
+
+class CRE51():
+    def __init__(self):
+        self.problem_name = 'CRE51'
+        self.n_objectives = 5
+        self.n_variables = 3
+        self.n_constraints = 7
+
+        self.lbound = np.zeros(self.n_variables)
+        self.ubound = np.zeros(self.n_variables)
+        self.lbound[0] = 0.01
+        self.lbound[1] = 0.01
+        self.lbound[2] = 0.01
+        self.ubound[0] = 0.45
+        self.ubound[1] = 0.10
+        self.ubound[2] = 0.10
+
+    def evaluate(self, x):
+        f = np.zeros(self.n_objectives)
+        g = np.zeros(self.n_constraints)
+
+        # First original objective function
+        f[0] = 106780.37 * (x[1] + x[2]) + 61704.67
+        # Second original objective function
+        f[1] = 3000 * x[0]
+        # Third original objective function
+        f[2] = 305700 * 2289 * x[1] / np.power(0.06 * 2289, 0.65)
+        # Fourth original objective function
+        f[3] = 250 * 2289 * np.exp(-39.75 * x[1] + 9.9 * x[2] + 2.74)
+        # Fifth original objective function
+        f[4] = 25 * (1.39 / (x[0] * x[1]) + 4940 * x[2] - 80)
+
+        # Constraint functions
+        g[0] = 1 - (0.00139 / (x[0] * x[1]) + 4.94 * x[2] - 0.08)
+        g[1] = 1 - (0.000306 / (x[0] * x[1]) + 1.082 * x[2] - 0.0986)
+        g[2] = 50000 - (12.307 / (x[0] * x[1]) + 49408.24 * x[2] + 4051.02)
+        g[3] = 16000 - (2.098 / (x[0] * x[1]) + 8046.33 * x[2] - 696.71)
+        g[4] = 10000 - (2.138 / (x[0] * x[1]) + 7883.39 * x[2] - 705.04)
+        g[5] = 2000 - (0.417 * x[0] * x[1] + 1721.26 * x[2] - 136.54)
+        g[6] = 550 - (0.164 / (x[0] * x[1]) + 631.13 * x[2] - 54.48)
+        g = np.where(g < 0, -g, 0)
+
+        return f, g
+
+
+
+
+if __name__ == '__main__':
+    np.random.seed(seed=1)
+    fun = RE21()
+
+    x = fun.lbound + (fun.ubound - fun.lbound) * np.random.rand(fun.n_variables)
+    print("Problem = {}".format(fun.problem_name))
+    print("Number of objectives = {}".format(fun.n_objectives))
+    print("Number of variables = {}".format(fun.n_variables))
+    print("Number of constraints = {}".format(fun.n_constraints))
+    print("Lower bounds = ", fun.lbound)
+    print("Upper bounds = ", fun.ubound)
+    print("x = ", x)
+
+    if 'CRE' in fun.problem_name:
+        f, g = fun.evaluate(x)
+        print("f(x) = {}".format(f))
+        print("g(x) = {}".format(g))
+    else:
+        f = fun.evaluate(x)
+        print("f(x) = {}".format(f))

bike_bench_internal/benchmark_models/libmoon/problem/synthetic/vlmop.py

@@ -0,0 +1,67 @@
+import matplotlib.pyplot as plt
+import torch
+import numpy as np
+
+from ..mop import mop
+
+
+
+class VLMOP1(mop):
+    def __init__( self, n_var=10, n_obj=2, lbound=-np.ones(10), ubound=np.ones(10) ):
+        super().__init__(n_var=n_var,
+                         n_obj=n_obj,
+                         lbound=lbound,
+                         ubound=ubound, )
+        self.problem_name = 'VLMOP1'
+
+    def _evaluate_torch(self, x):
+        f1 = 1 - torch.exp(-1 * torch.sum((x - 1 / np.sqrt(self.n_var))**2, dim=1))
+        f2 = 1 - torch.exp(-1 * torch.sum((x + 1 / np.sqrt(self.n_var))**2, dim=1))
+        return torch.stack((f1, f2), dim=1)
+
+    def _evaluate_numpy(self, x):
+        f1 = 1 - np.exp(-1 * np.sum((x - 1 / np.sqrt(self.n_var)) ** 2, axis=1 ) )
+        f2 = 1 - np.exp(-1 * np.sum((x + 1 / np.sqrt(self.n_var)) ** 2, axis=1 ) )
+        return np.stack((f1, f2), axis=1)
+
+
+    def get_pf(self):
+        x = torch.linspace(-1 / np.sqrt(self.n_var), 1 / np.sqrt(self.n_var), 100)
+        x = torch.tile(x.unsqueeze(1), (1, self.n_var))
+        with torch.no_grad():
+            return self._evaluate_torch(x).numpy()
+
+
+
+
+class VLMOP2(mop):
+    def __init__(self, n_var=10, n_obj=2, lbound=-np.ones(10), ubound=np.ones(10)):
+        super().__init__(n_var=n_var,
+                         n_obj=n_obj,
+                         lbound=lbound,
+                         ubound=ubound, )
+        self.problem_name = 'VLMOP2'
+
+    def _evaluate_torch(self, x):
+        f1 = torch.norm(x - 1 / np.sqrt(self.n_var), dim=1)**2 / 4
+        f2 = torch.norm(x + 1 / np.sqrt(self.n_var), dim=1)**2 / 4
+        return torch.stack((f1, f2), dim=1)
+
+    def _evaluate_numpy(self, x):
+        f1 = np.linalg.norm(x - 1 / np.sqrt(self.n_var), axis=1)**2 / 4
+        f2 = np.linalg.norm(x + 1 / np.sqrt(self.n_var), axis=1)**2 / 4
+        return np.stack((f1, f2), axis=1)
+
+    def get_pf(self):
+        x = torch.linspace(-1 / np.sqrt(self.n_var), 1 / np.sqrt(self.n_var), 100)
+        x = torch.tile(x.unsqueeze(1), (1, self.n_var))
+        with torch.no_grad():
+            return self._evaluate_torch(x).numpy()
+
+
+if __name__ == '__main__':
+    problem = VLMOP2()
+    pf = problem.get_pf()
+    plt.plot(pf[:, 0], pf[:, 1], '-')
+    plt.show()
+    print()

bike_bench_internal/benchmark_models/libmoon/problem/synthetic/wfg.py

@@ -0,0 +1,3 @@
+
+
+

bike_bench_internal/benchmark_models/libmoon/problem/synthetic/zdt.py

@@ -0,0 +1,175 @@
+"""
+ZDT test suite for multi-objective problem
+
+Reference
+----------
+    Zitzler, E., Deb, K., & Thiele, L. (2000). Comparison of multiobjective
+    evolutionary algorithms: Empirical results. Evolutionary computation,
+    8(2), 173-195. DOI: 10.1162/106365600568202
+"""
+import numpy as np
+import torch
+from matplotlib import pyplot as plt
+from ..mop import mop
+
+class ZDT1(mop):
+
+    def __init__(self, n_var=30, n_obj=2, lbound=np.zeros(30),
+                 ubound=np.ones(30)):
+        super().__init__(n_var=n_var,
+                         n_obj=n_obj,
+                         lbound=lbound,
+                         ubound=ubound, )
+        self.problem_name = 'ZDT1'
+
+
+    def _evaluate_torch(self, x: torch.Tensor):
+        f1 = x[:, 0]
+        n = len(x)
+        g = 1 + 9/(n-1) * torch.sum(x[:, 1:], dim=1)
+        h = 1 - torch.sqrt(f1 / g)
+        f2 = g * h
+        return torch.stack((f1, f2), dim=1)
+
+    def _evaluate_numpy(self, x: np.ndarray):
+        n = len(x)
+        f1 = x[:, 0]
+        g = 1 + 9 / (n-1) * np.sum(x[:, 1:], axis=1)
+        f2 = 1 - np.sqrt(f1 / g)
+        return np.stack((f1, f2), axis=1)
+
+    def _get_pf(self, n_points: int = 100):
+        f1 = np.linspace(0, 1, n_points)
+        f2 = 1 - np.sqrt(f1)
+        return np.stack((f1, f2), axis=1)
+
+
+class ZDT2(mop):
+
+    def __init__(self, n_var=30, n_obj=2, lbound=np.zeros(30), ubound=np.ones(30)):
+        super().__init__(n_var=n_var,
+                         n_obj=n_obj,
+                         lbound=lbound,
+                         ubound=ubound)
+        self.problem_name = 'ZDT2'
+
+    def _evaluate_torch(self, x: torch.Tensor):
+        f1 = x[:, 0]
+        g = 1 + 9 * torch.mean(x[:, 1:], dim=1)
+        f2 = g * (1 - (f1 / g) ** 2)
+        return torch.stack((f1, f2), dim=1)
+
+    def _evaluate_numpy(self, x: np.ndarray):
+        f1 = x[:, 0]
+        g = 1 + 9 * np.mean(x[:, 1:], axis=1)
+        f2 = g * (1 - (f1 / g) ** 2)
+        return np.stack((f1, f2), axis=1)
+
+    def _get_pf(self, n_points: int = 100):
+        f1 = np.linspace(0, 1, n_points)
+        f2 = 1 - f1 ** 2
+        return np.stack((f1, f2), axis=1)
+
+
+class ZDT3(mop):
+
+    def __init__(self, n_var=30, n_obj=2, lbound=np.zeros(30), ubound=np.ones(30)):
+        super().__init__(n_var=n_var,
+                         n_obj=n_obj,
+                         lbound=lbound,
+                         ubound=ubound, )
+        self.problem_name = 'ZDT3'
+
+
+    def _evaluate_torch(self, x: torch.Tensor):
+        f1 = x[:, 0]
+        g = 1 + 9 * torch.mean(x[:, 1:], dim=1)
+        f2 = g * (1 - torch.sqrt(f1 / g) - f1 / g * torch.sin(10 * np.pi * f1))
+        return torch.stack((f1, f2), dim=1)
+
+    def _evaluate_numpy(self, x: np.ndarray):
+        f1 = x[:, 0]
+        g = 1 + 9 * np.mean(x[:, 1:], axis=1)
+        f2 = g * (1 - np.sqrt(f1 / g) - f1 / g * np.sin(10 * np.pi * f1))
+        return np.stack((f1, f2), axis=1)
+
+    def _get_pf(self, n_points: int = 100):
+        f1 = np.hstack([np.linspace(0, 0.0830, int(n_points / 5)),
+                        np.linspace(0.1822, 0.2578, int(n_points / 5)),
+                        np.linspace(0.4093, 0.4539, int(n_points / 5)),
+                        np.linspace(0.6183, 0.6525, int(n_points / 5)),
+                        np.linspace(0.8233, 0.8518, n_points - 4 * int(n_points / 5))])
+        f2 = 1 - np.sqrt(f1) - f1 * np.sin(10 * np.pi * f1)
+        return np.stack((f1, f2), axis=1)
+
+
+class ZDT4(mop):
+
+    def __init__(self, n_var=10, n_obj=2, lbound=-5*np.ones(10), ubound=5*np.ones(10)):
+        lbound[0] = 0
+        ubound[0] = 1
+
+        super().__init__(n_var=n_var,
+                         n_obj=n_obj,
+                         lbound=lbound,
+                         ubound=ubound, )
+        self.problem_name = 'ZDT4'
+
+
+    def _evaluate_torch(self, x: torch.Tensor):
+        f1 = x[:, 0]
+        g = 1 + 10 * (self.n_var - 1) + torch.sum(x[:, 1:] ** 2 - 10 * torch.cos(4 * np.pi * x[:, 1:]), dim=1)
+        f2 = g * (1 - torch.sqrt(f1 / g))
+        return torch.stack((f1, f2), dim=1)
+
+    def _evaluate_numpy(self, x: np.ndarray):
+        f1 = x[:, 0]
+        g = 1 + 10 * (self.n_var - 1) + np.sum(x[:, 1:] ** 2 - 10 * np.cos(4 * np.pi * x[:, 1:]), axis=1)
+        f2 = g * (1 - np.sqrt(f1 / g))
+        return np.stack((f1, f2), axis=1)
+
+    def _get_pf(self, n_points: int = 100):
+        f1 = np.linspace(0, 1, n_points)
+        f2 = 1 - np.sqrt(f1)
+        return np.stack((f1, f2), axis=1)
+
+
+class ZDT6(mop):
+
+    def __init__(self, n_var=30, n_obj=2, lbound=np.zeros(30), ubound=np.ones(30) ) -> None:
+        super().__init__(n_var=n_var,
+                         n_obj=n_obj,
+                         lbound=lbound,
+                         ubound=ubound, )
+        self.problem_name = 'ZDT6'
+
+
+    def _evaluate_torch(self, x: torch.Tensor):
+        f1 = 1 - torch.exp(-4 * x[:, 0]) * (torch.sin(6 * np.pi * x[:, 0])) ** 6
+        g = 1 + 9 * (torch.sum(x[:, 1:], dim=1) / (self.n_var - 1)) ** 0.25
+        f2 = g * (1 - (f1 / g) ** 2)
+        return torch.stack((f1, f2), dim=1)
+
+    def _evaluate_numpy(self, x: np.ndarray):
+        f1 = 1 - np.exp(-4 * x[:, 0]) * (np.sin(6 * np.pi * x[:, 0])) ** 6
+        g = 1 + 9 * (np.sum(x[:, 1:], axis=1) / (self.n_var - 1)) ** 0.25
+        f2 = g * (1 - (f1 / g) ** 2)
+        return np.stack((f1, f2), axis=1)
+
+    def _get_pf(self, n_points: int = 100):
+        f1 = np.linspace(0, 1, n_points)
+        f2 = 1 - f1 ** 2
+        return np.stack((f1, f2), axis=1)
+
+
+
+if __name__ == '__main__':
+    problem = ZDT3()
+
+    res = problem.evaluate(torch.rand(10, problem.get_number_variable))
+    pf = problem.get_pf()
+    x = np.random.rand(10, problem.get_number_variable)
+    y = problem.evaluate(x)
+
+    plt.scatter(pf[:, 0], pf[:, 1], c='none', edgecolors='r')
+    plt.show()

bike_bench_internal/benchmark_models/libmoon/solver/__init__.py

@@ -0,0 +1 @@
+# from gradient.epo_solver import EPOSolver

bike_bench_internal/benchmark_models/libmoon/solver/gradient/__init__.py

@@ -0,0 +1,24 @@
+from .base_solver import GradAggSolver
+from .mgda_solver import MGDASolver
+from .epo_solver import EPOSolver
+from .moosvgd import MOOSVGDSolver
+from .gradhv import GradHVSolver
+from .pmtl import PMTLSolver
+
+
+
+from .core_solver import CoreAgg, CoreMGDA, CoreEPO, CoreMOOSVGD, CoreHVGrad
+
+def get_core_solver(args, pref=None):
+    if args.mtd == 'agg':
+        return CoreAgg(pref=pref, agg_mtd=args.agg_mtd)
+    elif args.mtd == 'mgda':
+        return CoreMGDA()
+    elif args.mtd == 'epo':
+        return CoreEPO(pref=pref)
+    elif args.mtd == 'moosvgd':
+        return CoreMOOSVGD(args=args)
+    elif args.mtd == 'hvgrad':
+        return CoreHVGrad(args=args)
+    else:
+        assert False, 'not implemented'

bike_bench_internal/benchmark_models/libmoon/solver/gradient/base_solver.py

@@ -0,0 +1,82 @@
+from numpy import array
+from torch.autograd import Variable
+from torch.optim import SGD
+from torch import Tensor
+from ...util_global.constant import scalar_dict, solution_eps, get_hv_ref_dict
+import torch
+from tqdm import tqdm
+from pymoo.indicators.hv import HV
+
+
+class GradBaseSolver:
+    def __init__(self, step_size, max_iter, tol):
+        self.step_size = step_size
+        self.max_iter = max_iter
+        self.tol = tol
+
+    def solve(self, problem, x, prefs, args):
+        '''
+            :param problem:
+            :param x:
+            :param agg:
+            :return:
+                is a dict with keys: x, y
+        '''
+
+        # The abstract class cannot be implemented directly.
+        raise NotImplementedError
+
+
+
+class GradAggSolver(GradBaseSolver):
+    def __init__(self, step_size, max_iter, tol, device='cpu'):
+        self.device = device
+        super().__init__(step_size, max_iter, tol)
+
+    def solve(self, problem, x, prefs, args, ref_point):
+        x = torch.tensor(x, dtype=torch.float32, requires_grad=True, device=self.device)
+
+        # ref_point = array([2.0, 2.0])
+        # ind = HV(ref_point = get_hv_ref_dict(args.problem_name))
+        # ind = HV(ref_point = array([1.0, 1.0]))
+        
+        # ind = HV(ref_point = ref_point)
+    
+
+
+
+        # hv_arr = []
+        y_arr = []
+
+        if not isinstance(prefs, torch.Tensor):
+            prefs = torch.tensor(prefs, dtype=torch.float32, device=self.device)
+        else:
+            prefs = prefs.to(dtype=torch.float32, device=self.device)
+
+        # prefs = Tensor(prefs)
+        optimizer = SGD([x], lr=self.step_size)
+        agg_func = scalar_dict[args.agg]
+        res = {}
+        for i in tqdm(range(self.max_iter)):
+            y = problem.evaluate(x)
+            
+            # hv_arr.append(ind.do(y.detach().cpu().numpy()))
+
+            agg_val = agg_func(y, prefs)
+            optimizer.zero_grad()
+            torch.sum(agg_val).backward()
+            optimizer.step()
+
+            y_arr.append(y.detach().cpu().numpy())
+
+            if 'lbound' in dir(problem):
+                x.data = torch.clamp(x.data, 
+                     torch.tensor(problem.lbound, device=x.device, dtype=torch.float32) + solution_eps, 
+                     torch.tensor(problem.ubound, device=x.device, dtype=torch.float32) - solution_eps)
+
+
+        res['x'] = x.detach().cpu().numpy()
+        res['y'] = y.detach().cpu().numpy()
+        # res['hv_arr'] = hv_arr
+        res['y_arr'] = y_arr
+        return res

bike_bench_internal/benchmark_models/libmoon/solver/gradient/core_solver.py

@@ -0,0 +1,82 @@
+import numpy as np
+from .mgda_core import solve_mgda
+from .epo_solver import EPO_LP
+import torch
+from .gradhv import HvMaximization
+from ...util_global.constant import get_hv_ref_dict
+
+
+
+class CoreHVGrad:
+    def __init__(self, args):
+        self.args = args
+        self.hv_solver = HvMaximization(args.n_prob, args.n_obj, get_hv_ref_dict(args.problem))
+
+    def get_alpha(self, losses):
+        alpha = self.hv_solver.compute_weights(losses).T
+        return alpha
+
+
+
+class CoreMOOSVGD:
+    def __init__(self):
+        pass
+
+    def get_alpha(self):
+        return 0
+
+
+class CoreMGDA:
+    def __init__(self):
+        pass
+
+    def get_alpha(self, G, losses=None, pref=None):
+        _, alpha = solve_mgda(G, return_coeff=True)
+        return alpha
+
+
+class CoreGrad:
+    def __init__(self):
+        pass
+
+
+
+class CoreEPO(CoreGrad):
+    def __init__(self, pref):
+        self.pref = pref
+        self.epo_lp = EPO_LP(m=len(pref), n=1, r=1/np.array(pref))
+
+
+    def get_alpha(self, G, losses):
+        if type(G) == torch.Tensor:
+            G = G.detach().cpu().numpy().copy()
+        GG = G @ G.T
+
+        alpha = self.epo_lp.get_alpha(losses, G=GG, C=True)
+        return alpha
+
+
+
+
+class CoreAgg(CoreGrad):
+    def __init__(self, pref, agg_mtd='ls'):
+        self.agg_mtd = agg_mtd
+        self.pref = pref
+
+    def get_alpha(self, G, losses):
+        if self.agg_mtd == 'ls':
+            alpha = self.pref
+        elif self.agg_mtd == 'mtche':
+            idx = np.argmax(losses)
+            alpha = np.zeros_like(self.pref )
+            alpha[idx] = 1.0
+        else:
+            assert False
+        return alpha
+
+
+
+if __name__ == '__main__':
+    agg = CoreAgg( pref=np.array([1, 0]) )
+
+

bike_bench_internal/benchmark_models/libmoon/solver/gradient/epo_solver.py

@@ -0,0 +1,198 @@
+import numpy as np
+import cvxpy as cp
+import cvxopt
+
+from .base_solver import GradBaseSolver
+from torch.autograd import Variable
+from tqdm import tqdm
+import torch
+from torch.optim import SGD
+from numpy import array
+from pymoo.indicators.hv import HV
+import warnings
+warnings.filterwarnings("ignore")
+
+from ...util_global.constant import solution_eps
+
+
+
+class EPO_LP(object):
+    def __init__(self, m, n, r, eps=1e-4):
+        cvxopt.glpk.options["msg_lev"] = "GLP_MSG_OFF"
+        self.m = m
+        self.n = n
+        self.r = r
+        self.eps = eps
+        self.last_move = None
+        self.a = cp.Parameter(m)        # Adjustments
+        self.C = cp.Parameter((m, m))   # C: Gradient inner products, G^T G
+        self.Ca = cp.Parameter(m)       # d_bal^TG
+        self.rhs = cp.Parameter(m)      # RHS of constraints for balancing
+
+        self.alpha = cp.Variable(m)     # Variable to optimize
+
+        obj_bal = cp.Maximize(self.alpha @ self.Ca)   # objective for balance
+        constraints_bal = [self.alpha >= 0, cp.sum(self.alpha) == 1,  # Simplex
+                           self.C @ self.alpha >= self.rhs]
+        self.prob_bal = cp.Problem(obj_bal, constraints_bal)  # LP balance
+
+        obj_dom = cp.Maximize(cp.sum(self.alpha @ self.C))  # obj for descent
+        constraints_res = [self.alpha >= 0, cp.sum(self.alpha) == 1,  # Restrict
+                           self.alpha @ self.Ca >= -cp.neg(cp.max(self.Ca)),
+                           self.C @ self.alpha >= 0]
+        constraints_rel = [self.alpha >= 0, cp.sum(self.alpha) == 1,  # Relaxed
+                           self.C @ self.alpha >= 0]
+        self.prob_dom = cp.Problem(obj_dom, constraints_res)  # LP dominance
+        self.prob_rel = cp.Problem(obj_dom, constraints_rel)  # LP dominance
+
+        self.gamma = 0     # Stores the latest Optimum value of the LP problem
+        self.mu_rl = 0     # Stores the latest non-uniformity
+
+
+    def get_alpha(self, l, G, r=None, C=False, relax=False):
+
+        r = self.r if r is None else r
+        assert len(l) == len(G) == len(r) == self.m, "length != m"
+
+        rl, self.mu_rl, self.a.value = adjustments(l, r)
+
+        self.C.value = G if C else G @ G.T
+        self.Ca.value = self.C.value @ self.a.value
+
+
+        if self.mu_rl > self.eps:
+            J = self.Ca.value > 0
+            if len(np.where(J)[0]) > 0:
+                J_star_idx = np.where(rl == np.max(rl))[0]
+                self.rhs.value = self.Ca.value.copy()
+                self.rhs.value[J] = -np.inf     # Not efficient; but works.
+                self.rhs.value[J_star_idx] = 0
+            else:
+                self.rhs.value = np.zeros_like(self.Ca.value)
+            self.gamma = self.prob_bal.solve(solver=cp.GLPK, verbose=False)
+            # self.gamma = self.prob_bal.solve(verbose=False)
+            self.last_move = "bal"
+        else:
+            if relax:
+                self.gamma = self.prob_rel.solve(solver=cp.GLPK, verbose=False)
+            else:
+                self.gamma = self.prob_dom.solve(solver=cp.GLPK, verbose=False)
+            # self.gamma = self.prob_dom.solve(verbose=False)
+            self.last_move = "dom"
+
+        return self.alpha.value
+
+
+def mu(rl, normed=False):
+    if len(np.where(rl < 0)[0]):
+        raise ValueError(f"rl<0 \n rl={rl}")
+        return None
+    m = len(rl)
+    l_hat = rl if normed else rl / rl.sum()
+    eps = np.finfo(rl.dtype).eps
+    l_hat = l_hat[l_hat > eps]
+    return np.sum(l_hat * np.log(l_hat * m))
+
+
+def adjustments(l, r=1):
+    m = len(l)
+    rl = r * l
+    l_hat = rl / rl.sum()
+    mu_rl = mu(l_hat, normed=True)
+
+    eps = 1e-3   # clipping by eps is to avoid log(0), zxy Dec. 5.
+    a = r * ( np.log( np.clip(l_hat * m, eps, np.inf) ) - mu_rl)
+    return rl, mu_rl, a
+
+
+
+def solve_epo(grad_arr, losses, pref, epo_lp):
+
+    '''
+        input: grad_arr: (m,n).
+        losses : (m,).
+        pref: (m,) inv.
+
+        return : gw: (n,).
+    '''
+    if type(pref) == torch.Tensor:
+        pref = pref.numpy()
+
+    pref = np.array(pref)
+    G = grad_arr.detach().clone().numpy()
+
+    if type(losses) == torch.Tensor:
+        losses_np = losses.detach().clone().numpy().squeeze()
+    else:
+        losses_np = losses
+
+    m = G.shape[0]
+    n = G.shape[1]
+    GG = G @ G.T
+
+    # epo_lp = EPO_LP(m=m, n=n, r=np.array(pref))
+
+    alpha = epo_lp.get_alpha(losses_np, G=GG, C=True)
+    if alpha is None:   # A patch for the issue in cvxpy
+        alpha = pref / pref.sum()
+    gw = alpha @ G
+
+    # return torch.Tensor(gw).unsqueeze(0)
+    return torch.Tensor(gw), alpha
+
+
+
+
+
+
+
+class EPOSolver(GradBaseSolver):
+    def __init__(self, step_size, max_iter, tol):
+        super().__init__(step_size, max_iter, tol)
+
+
+    def solve(self, problem, x, prefs, args, ref_point):
+        x = Variable(x, requires_grad=True)
+
+        epo_arr = [  EPO_LP(m=args.n_obj, n=args.n_var, r=np.array( 1/pref )) for pref in prefs ]
+        optimizer = SGD([x], lr=self.step_size)
+
+        ind = HV(ref_point=ref_point)
+        hv_arr = []
+        y_arr = []
+
+
+        for i in tqdm( range(self.max_iter) ):
+
+            # optimizer.zero_grad()
+            y = problem.evaluate(x)
+            y_arr.append(y.detach().numpy() )
+
+            alpha_arr = [0] * args.n_prob
+            for prob_idx in range( args.n_prob ):
+                grad_arr = [0] * args.n_obj
+                for obj_idx in range(args.n_obj):
+                    y[prob_idx][obj_idx].backward(retain_graph=True)
+                    grad_arr[obj_idx] = x.grad[prob_idx].clone()
+                    x.grad.zero_()
+
+                grad_arr = torch.stack(grad_arr)
+                _, alpha = solve_epo(grad_arr, losses=y[prob_idx], pref=prefs[prob_idx], epo_lp=epo_arr[prob_idx])
+                alpha_arr[prob_idx] = alpha
+
+            optimizer.zero_grad()
+            alpha_arr = torch.Tensor( np.array(alpha_arr) )
+            torch.sum(alpha_arr * y).backward()
+            optimizer.step()
+
+            if 'lbound' in dir(problem):
+                x.data = torch.clamp(x.data, torch.Tensor(problem.lbound) + solution_eps, torch.Tensor(problem.ubound)-solution_eps )
+
+
+        res = {}
+        res['x'] = x.detach().numpy()
+        res['y'] = y.detach().numpy()
+        res['hv_arr'] = [0]
+        res['y_arr'] = y_arr
+
+        return res

bike_bench_internal/benchmark_models/libmoon/solver/gradient/functions_evaluation.py

@@ -0,0 +1,94 @@
+"""
+    The function fastNonDominatedSort is based on the sorting algorithm described by
+    Deb, Kalyanmoy, et al.
+    "A fast and elitist multiobjective genetic algorithm: NSGA-II."
+    IEEE transactions on evolutionary computation 6.2 (2002): 182-197.
+"""
+
+
+import numpy as np
+import pdb
+from .functions_hv_python3 import HyperVolume
+
+
+def determine_non_dom_mo_sol(mo_obj_val):
+    # get set of non-dominated solutions, returns indices of non-dominated and booleans of dominated mo_sol
+    n_mo_sol = mo_obj_val.shape[1]
+    domination_rank = fastNonDominatedSort(mo_obj_val)
+    non_dom_indices = np.where(domination_rank == 0)
+    non_dom_indices = non_dom_indices[0] # np.where returns a tuple, so we need to get the array inside the tuple
+    # non_dom_mo_sol = mo_sol[:,non_dom_indices]
+    # non_dom_mo_obj_val = mo_obj_val[:,non_dom_indices]
+    mo_sol_is_non_dominated = np.zeros(n_mo_sol,dtype = bool)
+    mo_sol_is_non_dominated[non_dom_indices] = True
+    mo_sol_is_dominated = np.bitwise_not(mo_sol_is_non_dominated)
+    return(non_dom_indices,mo_sol_is_dominated)
+
+def fastNonDominatedSort(objVal):
+    # Based on Deb et al. (2002) NSGA-II
+    N_OBJECTIVES = objVal.shape[0]
+    N_SOLUTIONS = objVal.shape[1]
+
+    rankIndArray = - 999 * np.ones(N_SOLUTIONS, dtype = int) # -999 indicates unassigned rank
+    solIndices = np.arange(0,N_SOLUTIONS) # array of 0 1 2 ... N_SOLUTIONS
+    ## compute the entire domination matrix
+    # dominationMatrix: (i,j) is True if solution i dominates solution j
+    dominationMatrix = np.zeros((N_SOLUTIONS,N_SOLUTIONS), dtype = bool)
+    for p in solIndices:
+        objValA = objVal[:,p][:,None] # add [:,None] to preserve dimensions
+        # objValArray =  np.delete(objVal, obj = p axis = 1) # dont delete solution p because it messes up indices
+        dominates = checkDomination(objValA,objVal)
+        dominationMatrix[p,:] = dominates
+
+    # count the number of times a solution is dominated
+    dominationCounter = np.sum(dominationMatrix, axis = 0)
+
+    ## find rank 0 solutions to initialize loop
+    isRankZero = (dominationCounter == 0) # column and row binary indices of solutions that are rank 0
+
+    rankZeroRowInd = solIndices[isRankZero]
+    # mark rank 0's solutions by -99 so that they are not considered as members of next rank
+    dominationCounter[rankZeroRowInd] = -99
+    # initialize rank counter at 0
+    rankCounter = 0
+    # assign solutions in rank 0 rankIndArray = 0
+    rankIndArray[isRankZero] = rankCounter
+
+    isInCurRank = isRankZero
+    # while the current rank is not empty
+    while not (np.sum(isInCurRank) == 0):
+        curRankRowInd = solIndices[isInCurRank] # column and row numbers of solutions that are in current rank
+        # for each solution in current rank
+        for p in curRankRowInd:
+            # decrease domination counter of each solution dominated by solution p which is in the current rank
+            dominationCounter[dominationMatrix[p,:]] -= 1 #dominationMatrix[p,:] contains indices of the solutions dominated by p
+        # all solutions that now have dominationCounter == 0, are in the next rank
+        isInNextRank = (dominationCounter == 0)
+        rankIndArray[isInNextRank] = rankCounter + 1
+        # mark next rank's solutions by -99 so that they are not considered as members of future ranks
+        dominationCounter[isInNextRank] = -99
+        # increase front counter
+        rankCounter += 1
+        # check which solutions are in current rank (next rank became current rank)
+        isInCurRank = (rankIndArray == rankCounter)
+        if not np.all(isInNextRank == isInCurRank): # DEBUGGING, if it works fine, replace above assignment
+            pdb.set_trace()
+    return(rankIndArray)
+
+def checkDomination(objValA,objValArray):
+    dominates = ( np.any(objValA < objValArray, axis = 0) & np.all(objValA <= objValArray , axis = 0) )
+    return(dominates)
+
+def compute_hv_in_higher_dimensions(mo_obj_val,ref_point):
+    n_mo_obj = mo_obj_val.shape[0]
+    n_mo_sol = mo_obj_val.shape[1]
+    assert len(ref_point) == n_mo_obj
+    # initialize hv computation instance
+    hv_computation_instance = HyperVolume(tuple(ref_point))
+    # turn numpy array to list of tuples
+    list_of_mo_obj_val = list()
+    for i_mo_sol in range(n_mo_sol):
+        list_of_mo_obj_val.append(tuple(mo_obj_val[:,i_mo_sol]))
+
+    hv = float(hv_computation_instance.compute(list_of_mo_obj_val))
+    return(hv)

bike_bench_internal/benchmark_models/libmoon/solver/gradient/functions_hv_grad_3d.py

@@ -0,0 +1,215 @@
+'''
+The function grad_multi_sweep and its subfunctions in this file are based on the algorithm described in:
+Emmerich, Michael, and André Deutz.
+"Time complexity and zeros of the hypervolume indicator gradient field."
+EVOLVE-a bridge between probability, set oriented numerics,
+and evolutionary computation III. Springer, Heidelberg, 2014. 169-193.
+'''
+import numpy as np
+import copy
+from .functions_hv_python3 import HyperVolume
+from .functions_evaluation import determine_non_dom_mo_sol
+
+def determine_mo_sol_in_exterior(mo_obj_val,ref_point):
+    # select only mo-solutions that are in the exterior
+    ref_point_temp = ref_point[:,None] # add axis so that comparison works
+    exterior_booleans = np.any(mo_obj_val > ref_point_temp, axis = 0)
+    exterior_indices = np.where(exterior_booleans == True)
+    return(exterior_indices,exterior_booleans)
+
+def determine_mo_sol_in_interior(mo_obj_val,ref_point):
+    # select only mo-solutions that are in the interior
+    ref_point_temp = ref_point[:,None] # add axis so that comparison works
+    interior_booleans = np.all(mo_obj_val < ref_point_temp, axis = 0)
+    interior_indices = np.where(interior_booleans == True)
+    return(interior_indices,interior_booleans)
+
+def determine_mo_sol_on_ref_boundary(mo_obj_val,ref_point):
+    # select only mo-solutions that are on the reference boundary
+    ref_point_temp = ref_point[:,None] # add axis so that comparison works
+    boundary_booleans = np.logical_and( np.all(mo_obj_val <= ref_point_temp, axis = 0) , np.any(mo_obj_val == ref_point_temp, axis = 0) )
+    boundary_indices = np.where(boundary_booleans == True)
+    return(boundary_indices,boundary_booleans)
+
+def compute_domination_properties(mo_obj_val):
+    '''
+    compute properties needed for HV gradient computation
+    '''
+    n_mo_sol = mo_obj_val.shape[1]
+    # mo_sol i stricly dominates j (all entries in i  < j ) if strong_domination_matrix[i,j] = True
+    strong_domination_matrix = np.zeros((n_mo_sol,n_mo_sol), dtype = np.bool_)
+    for i in range(0,n_mo_sol):
+        cur_col = mo_obj_val[:,i][:,None]
+        strong_domination_matrix[i,:] = np.all(cur_col < mo_obj_val,axis = 0)
+    # mo_sol i weakly dominates j (all entries in i  <= j and at least one entry i < j ) if weak_domination_matrix[i,j] = True
+    weak_domination_matrix = np.zeros((n_mo_sol,n_mo_sol), dtype = np.bool_)
+    for i in range(0,n_mo_sol):
+        cur_col = mo_obj_val[:,i][:,None]
+        weak_domination_matrix[i,:] = np.logical_and( np.all(cur_col <= mo_obj_val,axis = 0) , np.any(cur_col < mo_obj_val,axis = 0) )
+    # a mo_sol i is strongly dominated if any other solutions j strongly dominates it (any True in the column strong_domination_matrix[:,i] )
+    is_strongly_dominated = np.any(strong_domination_matrix, axis = 0)
+    # a mo_sol i is weakly dominated if any other solutions j weakly dominates it (any True in the column weak_domination_matrix[:,i] )
+    is_weakly_dominated = np.any(weak_domination_matrix, axis = 0)
+    # weakly but not strongly dominated
+    is_weakly_but_not_strongly_dominated = np.logical_and( np.logical_not(is_strongly_dominated) , np.any(weak_domination_matrix, axis = 0) )
+
+    # no other solution weakly dominates it
+    is_not_weakly_dominated = np.logical_not(is_weakly_dominated)
+
+    # mo_sol i shares coordinate with j if at least 1 entry i = j
+    coordinate_sharing_matrix = np.zeros((n_mo_sol,n_mo_sol), dtype = np.bool_)
+    for i in range(0,n_mo_sol):
+        cur_col = mo_obj_val[:,i][:,None]
+        coordinate_sharing_matrix[i,:] = np.any(cur_col == mo_obj_val,axis = 0)
+
+    ## is not weakly dominated but share some coordinate with another not weakly dominated mo_sol
+    # set diagonal entries of coordinate_sharing_matrix to zero to make check work (we do not care of shared coordinates with itself)
+    subset_of_coordinate_sharing_matrix_with_diag_zero = copy.copy(coordinate_sharing_matrix)
+    np.fill_diagonal(subset_of_coordinate_sharing_matrix_with_diag_zero,False) # inplace operation
+    # select subset of columns so that the comparison is for all mo-solutions but only with respect to all non-weakly dominated solutions
+    subset_of_coordinate_sharing_matrix_with_diag_zero = subset_of_coordinate_sharing_matrix_with_diag_zero[:,is_not_weakly_dominated]
+    # check per row if any coordinate is shared with a non-weakly dominated solution
+    has_shared_coordinates = np.any(subset_of_coordinate_sharing_matrix_with_diag_zero,axis = 1)
+    is_not_weakly_dominated_and_shares_coordinates_with_other_not_weakly_dominated = np.logical_and( is_not_weakly_dominated , has_shared_coordinates  )
+
+    ## is not weakly dominated and share no coordinate with another not weakly dominated mo_sol
+    is_not_weakly_dominated_and_shares_no_coordinates_with_other_not_weakly_dominated = np.logical_and( is_not_weakly_dominated , np.logical_not(has_shared_coordinates)  )
+
+
+    return(is_strongly_dominated,is_weakly_dominated,is_weakly_but_not_strongly_dominated,is_not_weakly_dominated,is_not_weakly_dominated_and_shares_coordinates_with_other_not_weakly_dominated,is_not_weakly_dominated_and_shares_no_coordinates_with_other_not_weakly_dominated)
+
+
+def compute_subsets(mo_obj_val,ref_point):
+    is_strongly_dominated,is_weakly_dominated,is_weakly_but_not_strongly_dominated, \
+    is_not_weakly_dominated,is_not_weakly_dominated_and_shares_coordinates_with_other_not_weakly_dominated, \
+    is_not_weakly_dominated_and_shares_no_coordinates_with_other_not_weakly_dominated = compute_domination_properties(mo_obj_val)
+    ## Z: indices of mo-solutions for which all partial derivatives are zero
+    # E: in exterior of reference space
+    E,_ = determine_mo_sol_in_exterior(mo_obj_val,ref_point)
+    # S: that are strictly dominated
+    S = np.where(is_strongly_dominated == True)
+    # Z = E cup S
+    Z = np.union1d(E,S)
+
+    ## U: indices of mo-solutions for which some partial derivatives are undefined
+    # D: that are non-dominated but have duplicate coordinates
+    D = np.where(is_not_weakly_dominated_and_shares_coordinates_with_other_not_weakly_dominated == True)
+    # W: that are weakly dominated, i.e. not strictly dominated but there is a weakly better (Pareto dominating) solution
+    W = np.where(is_weakly_but_not_strongly_dominated == True)
+    # B: that are on the boundary of the reference space
+    B,_ = determine_mo_sol_on_ref_boundary(mo_obj_val,ref_point)
+    # # U = D cup (W \ E) cup (B \ S)
+    # U = np.union1d( np.union1d(D, np.setdiff1d(W,E)) , np.setdiff1d(B,S) )
+    #DEVIATION FROM EMMERICH & DEUTZ PAPER:
+    # U = (D \ E) cup (W \ E) cup (B \ S)
+    U = np.union1d( np.union1d( np.setdiff1d(D,E) , np.setdiff1d(W,E) ) , np.setdiff1d(B,S) )
+
+    ## P: indices of mo-solutions for which all partial derivatives are positive (negative??? which one makes sense in our case)
+    # N: that are not Pareto dominated AND have no duplicate coordinate with any other non-dominated solution
+    N = np.where(is_not_weakly_dominated_and_shares_no_coordinates_with_other_not_weakly_dominated == True)
+    # I: that are in the interior of the reference space
+    I,_ = determine_mo_sol_in_interior(mo_obj_val,ref_point)
+    # P = N intersect I
+    P = np.intersect1d(N,I)
+
+    return(P,U,Z)
+
+def grad_multi_sweep_with_duplicate_handling(mo_obj_val,ref_point):
+    # find unique mo_obj_val (it also sorts columns which is unnecessary but gets fixed in when using mapping_indices)
+    unique_mo_obj_val, mapping_indices = np.unique(mo_obj_val, axis = 1, return_inverse = True)
+    # compute hv_grad for unique mo-solutions
+    unique_hv_grad = grad_multi_sweep(unique_mo_obj_val,ref_point)
+    # assign the same gradients to duplicate mo_obj_val (and undo the unnecessary sorting)
+    hv_grad = unique_hv_grad[:,mapping_indices]
+    return(hv_grad)
+
+def grad_multi_sweep(mo_obj_val,ref_point):
+    '''
+    Based on:
+    Emmerich, Michael, and André Deutz.
+    "Time complexity and zeros of the hypervolume indicator gradient field."
+    EVOLVE-a bridge between probability, set oriented numerics,
+    and evolutionary computation III. Springer, Heidelberg, 2014. 169-193.
+    '''
+    n_obj = mo_obj_val.shape[0]
+    n_mo_sol = mo_obj_val.shape[1]
+    assert n_obj == len(ref_point)
+    hv_grad = np.zeros_like(mo_obj_val)
+
+
+
+    P, U, Z = compute_subsets(mo_obj_val,ref_point)
+    #####
+    if not (len(U) == 0):
+        # raise ValueError("Partial derivatives might be only one-sided in indice" + str(U))
+        print("Partial derivatives might be only one-sided in indices" + str(U))
+        print(mo_obj_val[:,U])
+        hv_grad[:,U] = 0
+
+    for k in range(0,n_obj):
+        temp_ref_point = copy.copy(ref_point)
+        temp_ref_point = np.delete(temp_ref_point,k,axis = 0)
+        temp_ref_point = tuple(temp_ref_point)
+        hv_instance = HyperVolume(temp_ref_point)
+        sorted_P = copy.copy(mo_obj_val[:,P])
+        # descending order sorting
+        sort_order= np.argsort(-sorted_P[k,:])
+        sorted_P = sorted_P[:,sort_order]
+        # remove k-th row
+        sorted_P = np.delete(sorted_P,k,0)
+        # initialize queue by turning array of columns into list of columns
+        Q = sorted_P.T.tolist() # it should be possible to delete this row, Q is overwritten in the next line
+        Q = list()
+        for i in range(sorted_P.shape[1]):
+            Q.append(tuple(sorted_P[:,i]))
+        queue_index = len(P) # this initialization is actually index of last queue entry +1. The +1 is convenient because the while loop will always update the index by -1, so it all matches up in the end
+        T = list()
+        while (len(Q) > 0):
+            # take last element in list
+            q = Q.pop()
+            # compute hypervolume contribution of q when added to T
+
+            if len(T) == 0:
+                T_with_q = list()
+                T_with_q.append(q)
+                hv_contribution = hv_instance.compute(T_with_q)
+            else:
+                T_with_q = copy.copy(T)
+                T_with_q.append(q)
+                hv_contribution = hv_instance.compute(T_with_q) - hv_instance.compute(T)
+            # queue_index is the index of q
+            queue_index = queue_index - 1 # -1 because the counter is always lagging 1 and it was initialized with +1
+            # mo_sol_index is the index of q in mo_obj_val
+            mo_sol_index = P[sort_order[queue_index]]
+            hv_grad[k,mo_sol_index] = hv_contribution
+
+            ## add q to T and remove points dominated by q
+            # initialize T by q in first iteration
+            if len(T) == 0:
+                T.append(q)
+                # T = q
+            else:
+               ## remove columns in T that are dominated by q
+               # loop through T
+                i = 0
+                while (i < len(T)):
+                    # remove entry if dominated by q
+                    if check_weak_domination_in_tuple(q,T[i]):
+                        del T[i]
+                    else:
+                        # if entry not deleted, move to next entry
+                        i = i + 1
+
+                # add q to T
+                T.append(q)
+    return(hv_grad)
+
+def check_weak_domination_in_tuple(tuple_A,tuple_B):
+    assert len(tuple_A) == len(tuple_B)
+    # initialize as True
+    A_weakly_dominates_B = True
+    for i in range(len(tuple_B)):
+        if tuple_A[i] > tuple_B[i]:
+            A_weakly_dominates_B = False
+            break
+    return(A_weakly_dominates_B)

bike_bench_internal/benchmark_models/libmoon/solver/gradient/functions_hv_python3.py

@@ -0,0 +1,274 @@
+"""
+Copyright (C) 2010 Simon Wessing
+TU Dortmund University
+
+This program is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+__author__ = "Simon Wessing"
+"""
+
+
+class HyperVolume:
+    """
+    Hypervolume computation based on variant 3 of the algorithm in the paper:
+    C. M. Fonseca, L. Paquete, and M. Lopez-Ibanez. An improved dimension-sweep
+    algorithm for the hypervolume indicator. In IEEE Congress on Evolutionary
+    Computation, pages 1157-1163, Vancouver, Canada, July 2006.
+    Minimization is implicitly assumed here!
+    """
+
+    def __init__(self, referencePoint):
+        """Constructor."""
+        self.referencePoint = referencePoint
+        self.list = []
+
+    def compute(self, front):
+        """Returns the hypervolume that is dominated by a non-dominated front.
+        Before the HV computation, front and reference point are translated, so
+        that the reference point is [0, ..., 0].
+        """
+
+        def weaklyDominates(point, other):
+            for i in range(len(point)):
+                if point[i] > other[i]:
+                    return False
+            return True
+
+        relevantPoints = []
+        referencePoint = self.referencePoint
+        dimensions = len(referencePoint)
+        for point in front:
+            # only consider points that dominate the reference point
+            if weaklyDominates(point, referencePoint):
+                relevantPoints.append(point)
+        if any(referencePoint):
+            # shift points so that referencePoint == [0, ..., 0]
+            # this way the reference point doesn't have to be explicitly used
+            # in the HV computation
+            for j in range(len(relevantPoints)):
+                relevantPoints[j] = [relevantPoints[j][i] - referencePoint[i] for i in range(dimensions)]
+        self.preProcess(relevantPoints)
+        bounds = [-1.0e308] * dimensions
+        hyperVolume = self.hvRecursive(dimensions - 1, len(relevantPoints), bounds)
+        return hyperVolume
+
+
+    def hvRecursive(self, dimIndex, length, bounds):
+        """Recursive call to hypervolume calculation.
+        In contrast to the paper, the code assumes that the reference point
+        is [0, ..., 0]. This allows the avoidance of a few operations.
+        """
+        hvol = 0.0
+        sentinel = self.list.sentinel
+        if length == 0:
+            return hvol
+        elif dimIndex == 0:
+            # special case: only one dimension
+            # why using hypervolume at all?
+            return -sentinel.next[0].cargo[0]
+        elif dimIndex == 1:
+            # special case: two dimensions, end recursion
+            q = sentinel.next[1]
+            h = q.cargo[0]
+            p = q.next[1]
+            while p is not sentinel:
+                pCargo = p.cargo
+                hvol += h * (q.cargo[1] - pCargo[1])
+                if pCargo[0] < h:
+                    h = pCargo[0]
+                q = p
+                p = q.next[1]
+            hvol += h * q.cargo[1]
+            return hvol
+        else:
+            remove = self.list.remove
+            reinsert = self.list.reinsert
+            hvRecursive = self.hvRecursive
+            p = sentinel
+            q = p.prev[dimIndex]
+            while q.cargo != None:
+                if q.ignore < dimIndex:
+                    q.ignore = 0
+                q = q.prev[dimIndex]
+            q = p.prev[dimIndex]
+            while length > 1 and (q.cargo[dimIndex] > bounds[dimIndex] or q.prev[dimIndex].cargo[dimIndex] >= bounds[dimIndex]):
+                p = q
+                remove(p, dimIndex, bounds)
+                q = p.prev[dimIndex]
+                length -= 1
+            qArea = q.area
+            qCargo = q.cargo
+            qPrevDimIndex = q.prev[dimIndex]
+            if length > 1:
+                hvol = qPrevDimIndex.volume[dimIndex] + qPrevDimIndex.area[dimIndex] * (qCargo[dimIndex] - qPrevDimIndex.cargo[dimIndex])
+            else:
+                qArea[0] = 1
+                qArea[1:dimIndex+1] = [qArea[i] * -qCargo[i] for i in range(dimIndex)]
+            q.volume[dimIndex] = hvol
+            if q.ignore >= dimIndex:
+                qArea[dimIndex] = qPrevDimIndex.area[dimIndex]
+            else:
+                qArea[dimIndex] = hvRecursive(dimIndex - 1, length, bounds)
+                if qArea[dimIndex] <= qPrevDimIndex.area[dimIndex]:
+                    q.ignore = dimIndex
+            while p is not sentinel:
+                pCargoDimIndex = p.cargo[dimIndex]
+                hvol += q.area[dimIndex] * (pCargoDimIndex - q.cargo[dimIndex])
+                bounds[dimIndex] = pCargoDimIndex
+                reinsert(p, dimIndex, bounds)
+                length += 1
+                q = p
+                p = p.next[dimIndex]
+                q.volume[dimIndex] = hvol
+                if q.ignore >= dimIndex:
+                    q.area[dimIndex] = q.prev[dimIndex].area[dimIndex]
+                else:
+                    q.area[dimIndex] = hvRecursive(dimIndex - 1, length, bounds)
+                    if q.area[dimIndex] <= q.prev[dimIndex].area[dimIndex]:
+                        q.ignore = dimIndex
+            hvol -= q.area[dimIndex] * q.cargo[dimIndex]
+            return hvol
+
+
+    def preProcess(self, front):
+        """Sets up the list data structure needed for calculation."""
+        dimensions = len(self.referencePoint)
+        nodeList = MultiList(dimensions)
+        nodes = [MultiList.Node(dimensions, point) for point in front]
+        for i in range(dimensions):
+            # sort by dimension
+            nodes = sorted(nodes, key=lambda node: node.cargo[i])
+            nodeList.extend(nodes, i)
+        self.list = nodeList
+
+
+
+class MultiList:
+    """A special data structure needed by FonsecaHyperVolume.
+
+    It consists of several doubly linked lists that share common nodes. So,
+    every node has multiple predecessors and successors, one in every list.
+    """
+
+    class Node:
+
+        def __init__(self, numberLists, cargo=None):
+            self.cargo = cargo
+            self.next  = [None] * numberLists
+            self.prev = [None] * numberLists
+            self.ignore = 0
+            self.area = [0.0] * numberLists
+            self.volume = [0.0] * numberLists
+
+        def __str__(self):
+            return str(self.cargo)
+
+
+    def __init__(self, numberLists):
+        """Constructor.
+
+        Builds 'numberLists' doubly linked lists.
+        """
+        self.numberLists = numberLists
+        self.sentinel = MultiList.Node(numberLists)
+        self.sentinel.next = [self.sentinel] * numberLists
+        self.sentinel.prev = [self.sentinel] * numberLists
+
+
+    def __str__(self):
+        strings = []
+        for i in range(self.numberLists):
+            currentList = []
+            node = self.sentinel.next[i]
+            while node != self.sentinel:
+                currentList.append(str(node))
+                node = node.next[i]
+            strings.append(str(currentList))
+        stringRepr = ""
+        for string in strings:
+            stringRepr += string + "\n"
+        return stringRepr
+
+
+    def __len__(self):
+        """Returns the number of lists that are included in this MultiList."""
+        return self.numberLists
+
+
+    def getLength(self, i):
+        """Returns the length of the i-th list."""
+        length = 0
+        sentinel = self.sentinel
+        node = sentinel.next[i]
+        while node != sentinel:
+            length += 1
+            node = node.next[i]
+        return length
+
+
+    def append(self, node, index):
+        """Appends a node to the end of the list at the given index."""
+        lastButOne = self.sentinel.prev[index]
+        node.next[index] = self.sentinel
+        node.prev[index] = lastButOne
+        # set the last element as the new one
+        self.sentinel.prev[index] = node
+        lastButOne.next[index] = node
+
+
+    def extend(self, nodes, index):
+        """Extends the list at the given index with the nodes."""
+        sentinel = self.sentinel
+        for node in nodes:
+            lastButOne = sentinel.prev[index]
+            node.next[index] = sentinel
+            node.prev[index] = lastButOne
+            # set the last element as the new one
+            sentinel.prev[index] = node
+            lastButOne.next[index] = node
+
+
+    def remove(self, node, index, bounds):
+        """Removes and returns 'node' from all lists in [0, 'index'[."""
+        for i in range(index):
+            predecessor = node.prev[i]
+            successor = node.next[i]
+            predecessor.next[i] = successor
+            successor.prev[i] = predecessor
+            if bounds[i] > node.cargo[i]:
+                bounds[i] = node.cargo[i]
+        return node
+
+
+    def reinsert(self, node, index, bounds):
+        """
+        Inserts 'node' at the position it had in all lists in [0, 'index'[
+        before it was removed. This method assumes that the next and previous
+        nodes of the node that is reinserted are in the list.
+        """
+        for i in range(index):
+            node.prev[i].next[i] = node
+            node.next[i].prev[i] = node
+            if bounds[i] > node.cargo[i]:
+                bounds[i] = node.cargo[i]
+
+
+
+if __name__ == "__main__":
+
+    # Example:
+    referencePoint = [2, 2, 2]
+    hv = HyperVolume(referencePoint)
+    front = [[1,0,1], [0,1,0]]
+    volume = hv.compute(front)

bike_bench_internal/benchmark_models/libmoon/solver/gradient/gradhv.py

@@ -0,0 +1,122 @@
+"""
+The class HvMaximization is based on the algorithm described by
+Wang, Hao, et al.
+"Hypervolume indicator gradient ascent multimnist-objective optimization."
+International conference on evolutionary multimnist-criterion optimization. Springer, Cham, 2017.
+"""
+
+import numpy as np
+import torch
+
+from .functions_evaluation import fastNonDominatedSort
+from .functions_hv_grad_3d import grad_multi_sweep_with_duplicate_handling
+from .base_solver import GradBaseSolver
+from torch.autograd import Variable
+
+
+from tqdm import tqdm
+from pymoo.indicators.hv import HV
+from ...util_global.constant import solution_eps, get_hv_ref_dict
+
+"""
+The class HvMaximization is based on the algorithm described by
+Wang, Hao, et al.
+"Hypervolume indicator gradient ascent multi-objective optimization."
+International conference on evolutionary multi-criterion optimization. Springer, Cham, 2017.
+"""
+
+import numpy as np
+import torch
+
+from .functions_evaluation import fastNonDominatedSort
+from .functions_hv_grad_3d import grad_multi_sweep_with_duplicate_handling
+
+
+class HvMaximization(object):
+    """
+    Mo optimizer for calculating dynamic weights using higamo style hv maximization
+    based on Hao Wang et al.'s HIGA-MO
+    uses non-dominated sorting to create multiple fronts, and maximize hypervolume of each
+    """
+
+    def __init__(self, n_mo_sol, n_mo_obj, ref_point, obj_space_normalize=True):
+        super(HvMaximization, self).__init__()
+        self.name = 'hv_maximization'
+        self.ref_point = np.array(ref_point)
+        self.n_mo_sol = n_mo_sol
+        self.n_mo_obj = n_mo_obj
+        self.obj_space_normalize = obj_space_normalize
+
+    def compute_weights(self, mo_obj_val):
+        n_mo_obj = self.n_mo_obj
+        n_mo_sol = self.n_mo_sol
+
+        # non-dom sorting to create multiple fronts
+        hv_subfront_indices = fastNonDominatedSort(mo_obj_val)
+        dyn_ref_point = 1.1 * np.max(mo_obj_val, axis=1)
+        for i_obj in range(0, n_mo_obj):
+            dyn_ref_point[i_obj] = np.maximum(self.ref_point[i_obj], dyn_ref_point[i_obj])
+        number_of_fronts = np.max(hv_subfront_indices) + 1  # +1 because of 0 indexing
+
+        obj_space_multifront_hv_gradient = np.zeros((n_mo_obj, n_mo_sol))
+        for i_fronts in range(0, number_of_fronts):
+            # compute HV gradients for current front
+            temp_grad_array = grad_multi_sweep_with_duplicate_handling(mo_obj_val[:, (hv_subfront_indices == i_fronts)],
+                                                                       dyn_ref_point)
+            obj_space_multifront_hv_gradient[:, (hv_subfront_indices == i_fronts)] = temp_grad_array
+
+        # normalize the hv_gradient in obj space (||dHV/dY|| == 1)
+        normalized_obj_space_multifront_hv_gradient = np.zeros((n_mo_obj, n_mo_sol))
+        for i_mo_sol in range(0, n_mo_sol):
+            w = np.sqrt(np.sum(obj_space_multifront_hv_gradient[:, i_mo_sol] ** 2.0))
+            if np.isclose(w, 0):
+                w = 1
+            if self.obj_space_normalize:
+                normalized_obj_space_multifront_hv_gradient[:, i_mo_sol] = obj_space_multifront_hv_gradient[:,
+                                                                           i_mo_sol] / w
+            else:
+                normalized_obj_space_multifront_hv_gradient[:, i_mo_sol] = obj_space_multifront_hv_gradient[:, i_mo_sol]
+
+        dynamic_weights = torch.tensor(normalized_obj_space_multifront_hv_gradient, dtype=torch.float)
+        return (dynamic_weights)
+
+
+
+class GradHVSolver(GradBaseSolver):
+    def __init__(self, step_size, max_iter, tol):
+
+        super().__init__(step_size, max_iter, tol)
+
+    def solve(self, problem, x, prefs, args, ref_point):
+        if args.n_obj != 2:
+            assert False, 'hvgrad only supports 2 obj problem'
+
+        hv_maximizer = HvMaximization(args.n_prob, args.n_obj, ref_point)
+
+
+        x = Variable(x, requires_grad=True)
+        optimizer = torch.optim.SGD([x,], lr=self.step_size)
+        hv_ind = HV(ref_point= ref_point)
+        hv_arr = [0] * self.max_iter
+        y_arr=[]
+        for iter_idx in tqdm(range(self.max_iter)):
+            y = problem.evaluate(x)
+            y_np = y.detach().numpy()
+            y_arr.append(y_np)
+            hv_arr[iter_idx] = hv_ind.do(y_np)
+            weight = hv_maximizer.compute_weights(y_np.T)
+            weight = torch.tensor(weight.T, dtype=torch.float)
+            optimizer.zero_grad()
+            torch.sum(weight*y).backward()
+            optimizer.step()
+            if 'lbound' in dir(problem):
+                x.data = torch.clamp(x.data, torch.Tensor(problem.lbound) + solution_eps, torch.Tensor(problem.ubound)-solution_eps )
+
+
+        res = {}
+        res['x'] = x.detach().numpy()
+        res['y'] = y.detach().numpy()
+        res['hv_arr'] = hv_arr
+        res['y_arr'] = y_arr
+
+        return res

bike_bench_internal/benchmark_models/libmoon/solver/gradient/mgda_core.py

@@ -0,0 +1,72 @@
+import torch
+from numpy import array
+import numpy as np
+from cvxopt import matrix, solvers
+solvers.options['show_progress'] = False
+
+def solve_mgda_analy(grad_1, grad_2, return_coeff = False):
+    '''
+        Noted that, solve_mgda_analy only support 2-objective case.
+        grad_i.shape: (n,).
+        This function support grad_i as both Tensor and numpy.
+    '''
+
+    v1v1 = grad_1 @ grad_1
+    v2v2 = grad_2 @ grad_2
+    v1v2 = grad_1 @ grad_2
+
+    if v1v2 >= v1v1:
+        gamma = 0.999
+    elif v1v2 >= v2v2:
+        gamma = 0.001
+    else:
+        gamma = -1.0 * ((v1v2 - v2v2) / (v1v1 + v2v2 - 2 * v1v2))
+
+    coeff = array([gamma, 1-gamma])
+    gw = coeff[0] * grad_1 + coeff[1] * grad_2
+    if return_coeff:
+        return gw, coeff
+    else:
+        return gw
+
+
+def solve_mgda(G, return_coeff=False):
+    '''
+        input G: (m,n).
+        output gw (n,).
+        comments: This function is used to solve the dual MGDA problem. It can handle m>2.
+    '''
+    if type(G) == torch.Tensor:
+        G = G.detach().cpu().numpy().copy()
+
+
+    m = G.shape[0]
+    if m == 2:
+        return solve_mgda_analy(G[0], G[1], return_coeff=return_coeff)
+    else:
+        Q = G @ G.T
+        Q = matrix(np.float64(Q))
+        p = np.zeros(m)
+        A = np.ones(m)
+
+        A = matrix(A, (1, m))
+        b = matrix(1.0)
+
+        G_cvx = -np.eye(m)
+        h = [0.0] * m
+        h = matrix(h)
+
+        G_cvx = matrix(G_cvx)
+        p = matrix(p)
+        sol = solvers.qp(Q, p, G_cvx, h, A, b)
+
+        res = np.array(sol['x']).squeeze()
+        res = res / sum(res)  # important
+        gw = torch.Tensor( res @ G )
+
+        if return_coeff:
+            return gw, res
+        else:
+            return gw
+
+

bike_bench_internal/benchmark_models/libmoon/solver/gradient/mgda_solver.py

@@ -0,0 +1,71 @@
+import torch.autograd
+
+from .mgda_core import solve_mgda
+
+from .base_solver import GradBaseSolver
+from torch.autograd import Variable
+from torch.optim import SGD
+from tqdm import tqdm
+from torch import Tensor
+import numpy as np
+from ...util_global.constant import solution_eps, get_hv_ref_dict
+from pymoo.indicators.hv import HV
+
+
+'''
+    MGDA solver, published in: 
+    1. Multiple-gradient descent algorithm (MGDA) for multiobjective optimizationAlgorithme de descente à gradients multiples pour lʼoptimisation multiobjectif
+    2. Sener, Ozan, and Vladlen Koltun. "Multi-task learning as multimnist-objective optimization." Advances in neural information processing systems 31 (2018).
+'''
+
+
+class MGDASolver(GradBaseSolver):
+    def __init__(self, step_size, max_iter, tol):
+        super().__init__(step_size, max_iter, tol)
+
+
+    def solve(self, problem, x, prefs, args, ref_point):
+        x = Variable(x, requires_grad=True)
+        optimizer = SGD([x], lr=self.step_size)
+
+        ind = HV(ref_point=ref_point)
+        hv_arr = []
+        y_arr = []
+
+        for i in tqdm(range(self.max_iter)):
+            grad_arr = [0] * args.n_prob
+            y = problem.evaluate(x)
+            y_np = y.detach().numpy()
+            y_arr.append(y_np)
+            hv_arr.append(ind.do(y_np))
+
+            for prob_idx in range( args.n_prob ):
+                grad_arr[prob_idx] = [0] * args.n_obj
+                for obj_idx in range(args.n_obj):
+                    y[prob_idx][obj_idx].backward(retain_graph=True)
+                    grad_arr[prob_idx][obj_idx] = x.grad[prob_idx].clone()
+                    x.grad.zero_()
+                grad_arr[prob_idx] = torch.stack(grad_arr[prob_idx])
+
+            grad_arr = torch.stack(grad_arr)
+            gw_arr = [solve_mgda(G, return_coeff=True) for G in grad_arr]
+            optimizer.zero_grad()
+            weights = Tensor( np.array([gw[1] for gw in gw_arr]) )
+            # weights = Tensor( np.array([1.0, 0.0]) )
+            torch.sum(weights * y).backward()
+            optimizer.step()
+
+            if 'lbound' in dir(problem):
+                x.data = torch.clamp(x.data, torch.Tensor(problem.lbound) + solution_eps, torch.Tensor(problem.ubound) - solution_eps )
+
+        res = {}
+        res['x'] = x.detach().numpy()
+        res['y'] = y.detach().numpy()
+        res['hv_arr'] = hv_arr
+        res['y_arr'] = y_arr
+
+        return res
+
+
+if __name__ == '__main__':
+    print()

bike_bench_internal/benchmark_models/libmoon/solver/gradient/min_norm_solvers_numpy.py

@@ -0,0 +1,386 @@
+# This code is from
+# Multi-Task Learning as Multi-Objective Optimization
+# Ozan Sener, Vladlen Koltun
+# Neural Information Processing Systems (NeurIPS) 2018
+# https://github.com/intel-isl/MultiObjectiveOptimization
+
+import numpy as np
+import torch
+
+class MinNormSolver:
+    MAX_ITER = 250
+    STOP_CRIT = 1e-5
+
+    def _min_norm_element_from2(v1v1, v1v2, v2v2):
+        """
+        Analytical solution for min_{c} |cx_1 + (1-c)x_2|_2^2
+        d is the distance (objective) optimzed
+        v1v1 = <x1,x1>
+        v1v2 = <x1,x2>
+        v2v2 = <x2,x2>
+        """
+        if v1v2 >= v1v1:
+            # Case: Fig 1, third column
+            gamma = 0.999
+            cost = v1v1
+            return gamma, cost
+        if v1v2 >= v2v2:
+            # Case: Fig 1, first column
+            gamma = 0.001
+            cost = v2v2
+            return gamma, cost
+        # Case: Fig 1, second column
+        gamma = -1.0 * ((v1v2 - v2v2) / (v1v1 + v2v2 - 2 * v1v2))
+        cost = v2v2 + gamma * (v1v2 - v2v2)
+        return gamma, cost
+
+    def _min_norm_2d(vecs, dps):
+        """
+        Find the minimum norm solution as combination of two points
+        This is correct only in 2D
+        ie. min_c |\sum c_i x_i|_2^2 st. \sum c_i = 1 , 1 >= c_1 >= 0 for all i, c_i + c_j = 1.0 for some i, j
+        """
+        dmin = 1e8
+        for i in range(len(vecs)):
+            for j in range(i + 1, len(vecs)):
+                if (i, j) not in dps:
+                    dps[(i, j)] = 0.0
+                    for k in range(len(vecs[i])):
+                        dps[(i, j)] += torch.mul(vecs[i][k], vecs[j][k]).sum().data.cpu()
+                    dps[(j, i)] = dps[(i, j)]
+                if (i, i) not in dps:
+                    dps[(i, i)] = 0.0
+                    for k in range(len(vecs[i])):
+                        dps[(i, i)] += torch.mul(vecs[i][k], vecs[i][k]).sum().data.cpu()
+                if (j, j) not in dps:
+                    dps[(j, j)] = 0.0
+                    for k in range(len(vecs[i])):
+                        dps[(j, j)] += torch.mul(vecs[j][k], vecs[j][k]).sum().data.cpu()
+                c, d = MinNormSolver._min_norm_element_from2(dps[(i, i)], dps[(i, j)], dps[(j, j)])
+                if d < dmin:
+                    dmin = d
+                    sol = [(i, j), c, d]
+        return sol, dps
+
+    def _projection2simplex(y):
+        """
+        Given y, it solves argmin_z |y-z|_2 st \sum z = 1 , 1 >= z_i >= 0 for all i
+        """
+        m = len(y)
+        sorted_y = np.flip(np.sort(y), axis=0)
+        tmpsum = 0.0
+        tmax_f = (np.sum(y) - 1.0) / m
+        for i in range(m - 1):
+            tmpsum += sorted_y[i]
+            tmax = (tmpsum - 1) / (i + 1.0)
+            if tmax > sorted_y[i + 1]:
+                tmax_f = tmax
+                break
+        return np.maximum(y - tmax_f, np.zeros(y.shape))
+
+    def _next_point(cur_val, grad, n):
+        proj_grad = grad - (np.sum(grad) / n)
+        tm1 = -1.0 * cur_val[proj_grad < 0] / proj_grad[proj_grad < 0]
+        tm2 = (1.0 - cur_val[proj_grad > 0]) / (proj_grad[proj_grad > 0])
+
+        # tm1 = np.array(tm1)
+        # tm2 = np.array(tm2)
+
+        # skippers = np.sum(tm1 < 1e-7) + np.sum(tm2 < 1e-7)
+        t = 1
+        if len(tm1[tm1 > 1e-7]) > 0:
+            t = np.min(tm1[tm1 > 1e-7])
+        if len(tm2[tm2 > 1e-7]) > 0:
+            t = min(t, np.min(tm2[tm2 > 1e-7]))
+
+        next_point = proj_grad * t + cur_val
+        next_point = MinNormSolver._projection2simplex(next_point)
+        return next_point
+
+    def find_min_norm_element(vecs):
+        """
+            Given a list of vectors (vecs), this method finds the minimum norm element in the convex hull.
+            as min |u|_2 st. u = \sum c_i vecs[i] and \sum c_i = 1.
+            It is quite geometric, and the main idea is the fact that if d_{ij} = min |u|_2 st u = c x_i + (1-c) x_j; the solution lies in (0, d_{i,j}).
+            Hence, we find the best 2-task solution, and then run the projected gradient descent until convergence.
+        """
+        # Solution lying at the combination of two points
+        dps = {}   # What does dps mean?
+
+        init_sol, dps = MinNormSolver._min_norm_2d(vecs, dps)
+
+        n = len(vecs)
+        sol_vec = np.zeros(n)
+        sol_vec[init_sol[0][0]] = init_sol[1]
+        sol_vec[init_sol[0][1]] = 1 - init_sol[1]
+
+        if n < 3:
+            # This is optimal for n=2, so return the solution
+            return sol_vec, init_sol[2]
+
+        iter_count = 0
+
+        grad_mat = np.zeros((n, n))
+        for i in range(n):
+            for j in range(n):
+                grad_mat[i, j] = dps[(i, j)]
+
+        while iter_count < MinNormSolver.MAX_ITER:
+            grad_dir = -1.0 * np.dot(grad_mat, sol_vec)
+            new_point = MinNormSolver._next_point(sol_vec, grad_dir, n)
+            # Re-compute the inner products for line search
+            v1v1 = 0.0
+            v1v2 = 0.0
+            v2v2 = 0.0
+            for i in range(n):
+                for j in range(n):
+                    v1v1 += sol_vec[i] * sol_vec[j] * dps[(i, j)]
+                    v1v2 += sol_vec[i] * new_point[j] * dps[(i, j)]
+                    v2v2 += new_point[i] * new_point[j] * dps[(i, j)]
+            nc, nd = MinNormSolver._min_norm_element_from2(v1v1, v1v2, v2v2)
+            new_sol_vec = nc * sol_vec + (1 - nc) * new_point
+            change = new_sol_vec - sol_vec
+            change = np.array(change)
+
+            if np.sum( np.abs(change) ) < MinNormSolver.STOP_CRIT:
+                return sol_vec, nd
+
+            sol_vec = new_sol_vec
+
+    def find_min_norm_element_FW(vecs):
+        """
+        Given a list of vectors (vecs), this method finds the minimum norm element in the convex hull
+        as min |u|_2 st. u = \sum c_i vecs[i] and \sum c_i = 1.
+        It is quite geometric, and the main idea is the fact that if d_{ij} = min |u|_2 st u = c x_i + (1-c) x_j; the solution lies in (0, d_{i,j})
+        Hence, we find the best 2-task solution, and then run the Frank Wolfe until convergence
+        """
+        # Solution lying at the combination of two points
+        dps = {}
+        init_sol, dps = MinNormSolver._min_norm_2d(vecs, dps)
+
+        n = len(vecs)
+        sol_vec = np.zeros(n)
+        sol_vec[init_sol[0][0]] = init_sol[1]
+        sol_vec[init_sol[0][1]] = 1 - init_sol[1]
+
+        if n < 3:
+            # This is optimal for n=2, so return the solution
+            return sol_vec, init_sol[2]
+
+        iter_count = 0
+
+        grad_mat = np.zeros((n, n))
+        for i in range(n):
+            for j in range(n):
+                grad_mat[i, j] = dps[(i, j)]
+
+        while iter_count < MinNormSolver.MAX_ITER:
+            t_iter = np.argmin(np.dot(grad_mat, sol_vec))
+
+            v1v1 = np.dot(sol_vec, np.dot(grad_mat, sol_vec))
+            v1v2 = np.dot(sol_vec, grad_mat[:, t_iter])
+            v2v2 = grad_mat[t_iter, t_iter]
+
+            nc, nd = MinNormSolver._min_norm_element_from2(v1v1, v1v2, v2v2)
+            new_sol_vec = nc * sol_vec
+            new_sol_vec[t_iter] += 1 - nc
+
+            change = new_sol_vec - sol_vec
+            if np.sum(np.abs(change)) < MinNormSolver.STOP_CRIT:
+                return sol_vec, nd
+            sol_vec = new_sol_vec
+
+
+def gradient_normalizers(grads, losses, normalization_type):
+    gn = {}
+    if normalization_type == 'l2':
+        for t in grads:
+            gn[t] = np.sqrt(np.sum([gr.pow(2).sum().data.cpu() for gr in grads[t]]))
+    elif normalization_type == 'loss':
+        for t in grads:
+            gn[t] = losses[t]
+    elif normalization_type == 'loss+':
+        for t in grads:
+            gn[t] = losses[t] * np.sqrt(np.sum([gr.pow(2).sum().data.cpu() for gr in grads[t]]))
+    elif normalization_type == 'none':
+        for t in grads:
+            gn[t] = 1.0
+    else:
+        print('ERROR: Invalid Normalization Type')
+    return gn
+
+
+class MinNormSolverNumpy:
+    MAX_ITER = 250
+    STOP_CRIT = 1e-6
+
+    def _min_norm_element_from2(v1v1, v1v2, v2v2):
+        """
+        Analytical solution for min_{c} |cx_1 + (1-c)x_2|_2^2
+        d is the distance (objective) optimzed
+        v1v1 = <x1,x1>
+        v1v2 = <x1,x2>
+        v2v2 = <x2,x2>
+        """
+        if v1v2 >= v1v1:
+            # Case: Fig 1, third column
+            gamma = 0.999
+            cost = v1v1
+            return gamma, cost
+        if v1v2 >= v2v2:
+            # Case: Fig 1, first column
+            gamma = 0.001
+            cost = v2v2
+            return gamma, cost
+        # Case: Fig 1, second column
+        gamma = -1.0 * ((v1v2 - v2v2) / (v1v1 + v2v2 - 2 * v1v2))
+        cost = v2v2 + gamma * (v1v2 - v2v2)
+        return gamma, cost
+
+    def _min_norm_2d(vecs, dps):
+        """
+        Find the minimum norm solution as combination of two points
+        This solution is correct if vectors(gradients) lie in 2D
+        ie. min_c |\sum c_i x_i|_2^2 st. \sum c_i = 1 , 1 >= c_1 >= 0 for all i, c_i + c_j = 1.0 for some i, j
+        """
+        dmin = 1e8
+        for i in range(len(vecs)):
+            for j in range(i + 1, len(vecs)):
+                if (i, j) not in dps:
+                    dps[(i, j)] = 0.0
+                    dps[(i, j)] = np.dot(vecs[i], vecs[j])
+                    dps[(j, i)] = dps[(i, j)]
+                if (i, i) not in dps:
+                    dps[(i, i)] = 0.0
+                    dps[(i, i)] = np.dot(vecs[i], vecs[i])
+                if (j, j) not in dps:
+                    dps[(j, j)] = 0.0
+                    dps[(j, j)] = np.dot(vecs[j], vecs[j])
+                c, d = MinNormSolver._min_norm_element_from2(dps[(i, i)], dps[(i, j)], dps[(j, j)])
+                if d < dmin:
+                    dmin = d
+                    sol = [(i, j), c, d]
+        return sol, dps
+
+    def _projection2simplex(y):
+        """
+        Given y, it solves argmin_z |y-z|_2 st \sum z = 1 , 1 >= z_i >= 0 for all i
+        """
+        m = len(y)
+        sorted_y = np.flip(np.sort(y), axis=0)
+        tmpsum = 0.0
+        tmax_f = (np.sum(y) - 1.0) / m
+        for i in range(m - 1):
+            tmpsum += sorted_y[i]
+            tmax = (tmpsum - 1) / (i + 1.0)
+            if tmax > sorted_y[i + 1]:
+                tmax_f = tmax
+                break
+        return np.maximum(y - tmax_f, np.zeros(y.shape))
+
+    def _next_point(cur_val, grad, n):
+        proj_grad = grad - (np.sum(grad) / n)
+        tm1 = -1.0 * cur_val[proj_grad < 0] / proj_grad[proj_grad < 0]
+        tm2 = (1.0 - cur_val[proj_grad > 0]) / (proj_grad[proj_grad > 0])
+
+        skippers = np.sum(tm1 < 1e-7) + np.sum(tm2 < 1e-7)
+        t = 1
+        if len(tm1[tm1 > 1e-7]) > 0:
+            t = np.min(tm1[tm1 > 1e-7])
+        if len(tm2[tm2 > 1e-7]) > 0:
+            t = min(t, np.min(tm2[tm2 > 1e-7]))
+
+        next_point = proj_grad * t + cur_val
+        next_point = MinNormSolver._projection2simplex(next_point)
+        return next_point
+
+    def find_min_norm_element(vecs):
+        """
+        Given a list of vectors (vecs), this method finds the minimum norm element in the convex hull
+        as min |u|_2 st. u = \sum c_i vecs[i] and \sum c_i = 1.
+        It is quite geometric, and the main idea is the fact that if d_{ij} = min |u|_2 st u = c x_i + (1-c) x_j; the solution lies in (0, d_{i,j})
+        Hence, we find the best 2-task solution, and then run the projected gradient descent until convergence
+        """
+        # Solution lying at the combination of two points
+        dps = {}
+        init_sol, dps = MinNormSolver._min_norm_2d(vecs, dps)
+
+        n = len(vecs)
+        sol_vec = np.zeros(n)
+        sol_vec[init_sol[0][0]] = init_sol[1]
+        sol_vec[init_sol[0][1]] = 1 - init_sol[1]
+
+        if n < 3:
+            # This is optimal for n=2, so return the solution
+            return sol_vec, init_sol[2]
+
+        iter_count = 0
+
+        grad_mat = np.zeros((n, n))
+        for i in range(n):
+            for j in range(n):
+                grad_mat[i, j] = dps[(i, j)]
+
+        while iter_count < MinNormSolver.MAX_ITER:
+            grad_dir = -1.0 * np.dot(grad_mat, sol_vec)
+            new_point = MinNormSolver._next_point(sol_vec, grad_dir, n)
+            # Re-compute the inner products for line search
+            v1v1 = 0.0
+            v1v2 = 0.0
+            v2v2 = 0.0
+            for i in range(n):
+                for j in range(n):
+                    v1v1 += sol_vec[i] * sol_vec[j] * dps[(i, j)]
+                    v1v2 += sol_vec[i] * new_point[j] * dps[(i, j)]
+                    v2v2 += new_point[i] * new_point[j] * dps[(i, j)]
+            nc, nd = MinNormSolver._min_norm_element_from2(v1v1, v1v2, v2v2)
+            new_sol_vec = nc * sol_vec + (1 - nc) * new_point
+            change = new_sol_vec - sol_vec
+            if np.sum(np.abs(change)) < MinNormSolver.STOP_CRIT:
+                return sol_vec, nd
+            sol_vec = new_sol_vec
+        return sol_vec, nd
+
+    def find_min_norm_element_FW(vecs):
+        """
+        Given a list of vectors (vecs), this method finds the minimum norm element in the convex hull
+        as min |u|_2 st. u = \sum c_i vecs[i] and \sum c_i = 1.
+        It is quite geometric, and the main idea is the fact that if d_{ij} = min |u|_2 st u = c x_i + (1-c) x_j; the solution lies in (0, d_{i,j})
+        Hence, we find the best 2-task solution, and then run the Frank Wolfe until convergence
+        """
+        # Solution lying at the combination of two points
+        dps = {}
+        init_sol, dps = MinNormSolver._min_norm_2d(vecs, dps)
+
+        n = len(vecs)
+        sol_vec = np.zeros(n)
+        sol_vec[init_sol[0][0]] = init_sol[1]
+        sol_vec[init_sol[0][1]] = 1 - init_sol[1]
+
+        if n < 3:
+            # This is optimal for n=2, so return the solution
+            return sol_vec, init_sol[2]
+
+        iter_count = 0
+
+        grad_mat = np.zeros((n, n))
+        for i in range(n):
+            for j in range(n):
+                grad_mat[i, j] = dps[(i, j)]
+
+        while iter_count < MinNormSolver.MAX_ITER:
+            t_iter = np.argmin(np.dot(grad_mat, sol_vec))
+
+            v1v1 = np.dot(sol_vec, np.dot(grad_mat, sol_vec))
+            v1v2 = np.dot(sol_vec, grad_mat[:, t_iter])
+            v2v2 = grad_mat[t_iter, t_iter]
+
+            nc, nd = MinNormSolver._min_norm_element_from2(v1v1, v1v2, v2v2)
+            new_sol_vec = nc * sol_vec
+            new_sol_vec[t_iter] += 1 - nc
+
+            change = new_sol_vec - sol_vec
+            if np.sum(np.abs(change)) < MinNormSolver.STOP_CRIT:
+                return sol_vec, nd
+            sol_vec = new_sol_vec
+        return sol_vec, nd

bike_bench_internal/benchmark_models/libmoon/solver/gradient/moosvgd.py

@@ -0,0 +1,101 @@
+from .mgda_core import solve_mgda
+from .base_solver import GradBaseSolver
+
+import torch
+import math
+from torch.autograd import Variable
+from torch.optim import SGD
+import sys
+from ...util_global.constant import solution_eps
+from tqdm import tqdm
+
+def kernel_functional_rbf(losses):
+    '''
+        input losses: (n_prob, n_obj)
+        output kernel_matrix: (n_prob, n_prob)
+        comments: This function is used to compute the kernel matrix for SVGD.
+    '''
+    n = losses.shape[0]
+    # losses shape : (10,) * (3,)
+    pairwise_distance = torch.norm(losses[:, None] - losses, dim=2).pow(2)
+    h = median(pairwise_distance) / math.log(n)
+    A = 5e-6  # Noted, this bandwith parameter is important.
+    kernel_matrix = torch.exp(-pairwise_distance / A*h)  # 5e-6 for zdt1,2,3, zxy, Dec 5, 2023
+    return kernel_matrix
+
+def median(tensor):
+    """
+        torch.median() acts differently from np.median(). We want to simulate numpy implementation.
+    """
+    tensor = tensor.detach().flatten()
+    tensor_max = tensor.max()[None]
+    return (torch.cat((tensor, tensor_max)).median() + tensor.median()) / 2.
+
+def get_svgd_gradient(G, inputs, losses):
+    '''
+    :param G.shape: (n_prob, n_obj, n_var)
+    :param inputs.shape: (n_prob, n_var)
+    :param losses.shape: (n_prob, n_obj)
+    :return:
+    '''
+    n_prob = inputs.size(0)
+    # G shape (n_prob, n_obj, n_var)
+    g_w = [0] * n_prob
+
+    for idx in range(n_prob):
+        g_w[idx] = torch.Tensor( solve_mgda(G[idx], return_coeff=False) )
+
+    g_w = torch.stack(g_w)  # (n_prob, n_var)
+    # See https://github.com/activatedgeek/svgd/issues/1#issuecomment-649235844 for why there is a factor -0.5
+
+    kernel = kernel_functional_rbf(losses)
+    kernel_grad = -0.5 * torch.autograd.grad(kernel.sum(), inputs, allow_unused=True)[0]   # (n_prob, n_var)
+    gradient = (kernel.mm(g_w) - kernel_grad) / n_prob
+
+    return gradient
+
+
+
+class MOOSVGDSolver(GradBaseSolver):
+    def __init__(self, step_size, max_iter, tol):
+        super().__init__(step_size, max_iter, tol)
+
+    def solve(self, problem, x, prefs, args, ref_point):
+        x = Variable(x, requires_grad=True)
+        optimizer = SGD([x], lr=self.step_size)
+        for i in tqdm(range(self.max_iter)):
+            y = problem.evaluate(x)
+            grad_arr = [0] * args.n_prob
+            for prob_idx in range(args.n_prob):
+                grad_arr[prob_idx] = [0] * args.n_obj
+                for obj_idx in range(args.n_obj):
+                    y[prob_idx][obj_idx].backward(retain_graph=True)
+                    grad_arr[prob_idx][obj_idx] = x.grad[prob_idx].clone()
+                    x.grad.zero_()
+                grad_arr[prob_idx] = torch.stack(grad_arr[prob_idx])
+
+            grad_arr = torch.stack(grad_arr).detach()
+            gw = get_svgd_gradient(grad_arr, x, y)
+            optimizer.zero_grad()
+            x.grad = gw
+            optimizer.step()
+
+            if 'lbound' in dir(problem):
+                x.data = torch.clamp(x.data, torch.Tensor(problem.lbound) + solution_eps, torch.Tensor(problem.ubound) - solution_eps)
+
+        res={}
+        res['x'] = x.detach().numpy()
+        res['y'] = y.detach().numpy()
+        res['hv_arr'] = [0]
+        return res
+
+
+
+if __name__ == '__main__':
+    losses = torch.rand(10, 3)
+    kernel = kernel_functional_rbf(losses)
+
+
+
+
+

bike_bench_internal/benchmark_models/libmoon/solver/gradient/pmgda.py

@@ -0,0 +1,14 @@
+from .base_solver import GradBaseSolver
+
+
+
+class PMGDASolver(GradBaseSolver):
+    def __init__(self, step_size, max_iter, tol):
+        print('pmgda solver')
+        print()
+
+        super().__init__(step_size, max_iter, tol)
+
+    def solve(self, problem, x, prefs, args):
+        print()
+

bike_bench_internal/benchmark_models/libmoon/solver/gradient/pmtl.py

@@ -0,0 +1,147 @@
+import torch
+
+from .base_solver import GradBaseSolver
+from matplotlib import pyplot as plt
+from .min_norm_solvers_numpy import MinNormSolver
+import numpy as np
+
+
+
+from torch.autograd import Variable
+from tqdm import tqdm
+from ...util_global.constant import solution_eps
+
+from .mgda_core import solve_mgda
+
+
+
+
+def get_d_moomtl(grads):
+    """
+        calculate the gradient direction for MOO-MTL
+    """
+    nobj, dim = grads.shape
+    sol, nd = MinNormSolver.find_min_norm_element(grads)
+    return sol
+
+
+def get_d_paretomtl(grads, value, weights, i):
+    # calculate the gradient direction for Pareto MTL
+    nobj, dim = grads.shape
+
+    # check active constraints
+    normalized_current_weight = weights[i] / np.linalg.norm(weights[i])
+    normalized_rest_weights = np.delete(weights, (i), axis=0) / np.linalg.norm(np.delete(weights, (i), axis=0), axis=1,
+                                                                     keepdims=True)
+    # shape: (9, 2)
+    w = normalized_rest_weights - normalized_current_weight
+    # solve QP
+    gx = np.dot(w, value / np.linalg.norm(value))
+    idx = gx > 0
+    vec = np.concatenate((grads, np.dot(w[idx], grads)), axis=0)
+    # use MinNormSolver to solve QP
+    # vec.shape:
+
+    # sol, nd = MinNormSolver.find_min_norm_element( vec )
+    _, sol = solve_mgda( torch.Tensor(vec), return_coeff=True)
+
+    # reformulate ParetoMTL as linear scalarization method, return the weights
+    weight0 = sol[0] + np.sum(np.array([sol[j] * w[idx][j - 2, 0] for j in np.arange(2, 2 + np.sum(idx))]))
+    weight1 = sol[1] + np.sum(np.array([sol[j] * w[idx][j - 2, 1] for j in np.arange(2, 2 + np.sum(idx))]))
+    weight = np.stack([weight0, weight1])
+
+    return weight
+
+
+
+
+def get_d_paretomtl_init(grads, value, weights, i):
+    # calculate the gradient direction for Pareto MTL initialization
+    nobj, dim = grads.shape
+
+    # check active constraints
+    normalized_current_weight = weights[i] / np.linalg.norm(weights[i])
+    normalized_rest_weights = np.delete(weights, (i), axis=0) / np.linalg.norm(np.delete(weights, (i), axis=0), axis=1,
+                                                                               keepdims=True)
+    w = normalized_rest_weights - normalized_current_weight
+    gx = np.dot(w, value / np.linalg.norm(value))
+    idx = gx > 0
+
+
+    if np.sum(idx) <= 0:
+        return np.zeros(nobj)
+    if np.sum(idx) == 1:
+        sol = np.ones(1)
+    else:
+        vecs = np.dot(w[idx], grads)
+        _, sol = solve_mgda( torch.Tensor(vecs), return_coeff=True)
+        # print()
+
+    # calculate the weights
+    weight0 = np.sum(np.array([sol[j] * w[idx][j, 0] for j in np.arange(0, np.sum(idx))]))
+    weight1 = np.sum(np.array([sol[j] * w[idx][j, 1] for j in np.arange(0, np.sum(idx))]))
+    weight = np.stack([weight0, weight1])
+
+    return weight
+
+
+def circle_points(r, n):
+    # generate evenly distributed preference vector
+    circles = []
+    for r, n in zip(r, n):
+        t = np.linspace(0, 0.5 * np.pi, n)
+        x = r * np.cos(t)
+        y = r * np.sin(t)
+        circles.append(np.c_[x, y])
+    return circles
+
+
+
+
+class PMTLSolver(GradBaseSolver):
+    def __init__(self, step_size, max_iter, tol):
+        super().__init__(step_size, max_iter, tol)
+
+    def solve(self, problem, x, prefs, args, ref_point):
+        if args.n_obj != 2:
+            assert False, 'hvgrad only supports 2 obj problem'
+
+        x = Variable(x, requires_grad=True)
+        warmup_iter = self.max_iter // 5
+        optimizer = torch.optim.SGD([x], lr=self.step_size)
+
+        y_arr = []
+        for iter_idx in tqdm( range(self.max_iter) ):
+            y = problem.evaluate(x)
+            y_np = y.detach().numpy()
+            y_arr.append(y_np)
+
+            grad_arr = [0] * args.n_prob
+            for prob_idx in range(args.n_prob):
+                grad_arr[prob_idx] = [0] * args.n_obj
+                for obj_idx in range(args.n_obj):
+                    y[prob_idx][obj_idx].backward(retain_graph=True)
+                    grad_arr[prob_idx][obj_idx] = x.grad[prob_idx].clone()
+                    x.grad.zero_()
+                grad_arr[prob_idx] = torch.stack(grad_arr[prob_idx])
+
+            grad_arr = torch.stack(grad_arr)
+            grad_arr_np = grad_arr.detach().numpy()
+            if iter_idx < warmup_iter:
+                weights = [ get_d_paretomtl_init(grad_arr_np[i], y_np[i], prefs, i) for i in range(args.n_prob) ]
+            else:
+                weights = [ get_d_paretomtl(grad_arr_np[i], y_np[i], prefs, i) for i in range(args.n_prob) ]
+
+            optimizer.zero_grad()
+            torch.sum(torch.tensor(weights) * y).backward()
+            optimizer.step()
+
+            if 'lbound' in dir(problem):
+                x.data = torch.clamp(x.data, torch.Tensor(problem.lbound) + solution_eps, torch.Tensor(problem.ubound) - solution_eps)
+
+        res={}
+        res['x'] = x.detach().numpy()
+        res['y'] = y.detach().numpy()
+        res['hv_arr'] = [0]
+        res['y_arr'] = y_arr
+        return res

bike_bench_internal/benchmark_models/libmoon/solver/gradient/run/run_grad.py

@@ -0,0 +1,99 @@
+import numpy as np
+import os
+import sys
+print(f"vscode current run path is {os.getcwd()}")
+os.chdir(sys.path[0])
+print(f"set  py path as current path ")
+print(f"vscode current run path is {os.getcwd()}")
+
+from libmoon.solver.gradient import MGDASolver, GradAggSolver, EPOSolver, MOOSVGDSolver, GradHVSolver, PMTLSolver
+from libmoon.util_global.weight_factor.funs import uniform_pref
+from libmoon.util_global.constant import problem_dict
+from libmoon.visulization.view_res import vis_res, vedio_res
+import argparse
+import torch
+from matplotlib import pyplot as plt
+import pickle
+import time
+
+
+
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser( description= 'example script' )
+    parser.add_argument( '--n-partition', type=int, default=10 )
+    parser.add_argument( '--agg', type=str, default='tche')  # If solve is agg, then choose a specific agg method.
+    parser.add_argument('--solver', type=str, default='hvgrad')
+    # ['agg', 'epo', 'moosvgd', 'hvgrad', 'pmtl', 'mgda']
+    parser.add_argument( '--problem-name', type=str, default='VLMOP2')
+    parser.add_argument('--iter', type=int, default=2000)
+    parser.add_argument('--step-size', type=float, default=0.1)
+    parser.add_argument('--tol', type=float, default=1e-6)
+    parser.add_argument('--plt-pref-flag', type=str, default='N')
+
+    args = parser.parse_args()
+    problem = problem_dict[args.problem_name]
+    args.n_obj, args.n_var = problem.n_obj, problem.n_var
+    root_name = os.path.dirname(os.path.dirname(__file__))
+
+
+    if args.solver == 'mgda':
+        solver = MGDASolver(args.step_size, args.iter, args.tol)
+    elif args.solver == 'agg':
+        solver = GradAggSolver(args.step_size, args.iter, args.tol)
+    elif args.solver == 'epo':
+        solver = EPOSolver(args.step_size, args.iter, args.tol)
+    elif args.solver == 'moosvgd':
+        solver = MOOSVGDSolver(args.step_size, args.iter, args.tol)
+    elif args.solver == 'hvgrad':
+        solver = GradHVSolver(args.step_size, args.iter, args.tol)
+    elif args.solver == 'pmtl':
+        solver = PMTLSolver(args.step_size, args.iter, args.tol)
+    elif args.solver=='pmgda':
+        assert False, 'will be implemented soon'
+    else:
+        raise Exception('solver not supported')
+
+    if args.solver == 'agg':
+        args.folder_name = os.path.join(root_name, 'output', args.problem_name, '{}_{}'.format(args.solver, args.agg))
+    else:
+        args.folder_name = os.path.join(root_name, 'output', args.problem_name, args.solver)
+
+    os.makedirs(args.folder_name, exist_ok=True)
+    prefs = uniform_pref( args.n_partition, problem.n_obj, clip_eps=1e-2)
+
+    args.n_prob = len(prefs)
+    if 'lbound' in dir(problem):
+        if args.problem_name == 'VLMOP1':
+            x0 = torch.rand(args.n_prob, problem.n_var) * 2 / np.sqrt(problem.n_var) - 1 / np.sqrt(problem.n_var)
+        else:
+            x0 = torch.rand(args.n_prob, problem.n_var)
+    else:
+        x0 = torch.rand( args.n_prob, problem.n_var ) * 20 - 10
+
+
+    ts = time.time()
+    res = solver.solve( problem, x=x0, prefs=prefs, args=args)
+
+    elapsed = time.time() - ts
+    res['elapsed'] = elapsed
+
+
+    use_fig=False
+    if use_fig:
+        vis_res(res, problem, prefs, args)
+        fig_name = os.path.join(args.folder_name, 'res.svg')
+        plt.savefig(fig_name)
+        print('Save fig to %s' % fig_name)
+        plt.show()
+
+
+    use_vedio=True
+    if use_vedio:
+        vedio_res(res, problem, prefs, args)
+
+
+    pickle_name = os.path.join(args.folder_name, 'res.pkl')
+    with open(pickle_name, 'wb') as f:
+        pickle.dump(res, f)
+
+    print('Save pickle to %s' % pickle_name)