Merge pull request #10 from Libensemble/add/summer21_sims

Add matching sim_fs to summer21 gens section

README.md

@@ -72,7 +72,7 @@ and has [API documentation available online](https://libensemble.readthedocs.io/
 6. #### Distributed Consensus-based Optimization Methods
    *Generator and Allocator Functions*
 
-   Four generator functions, a matching allocator function, a consensus-subroutines module, and tests for each. 
+   Four generator functions, a matching allocator function, a consensus-subroutines module, a handful of test simulation functions, and tests for each. 
    All created by Caleb Ju at ANL as Given's associate Summer 2021. See the docstrings in each of the modules for more information.
 
    Last included in libEnsemble v0.9.3. See [here](https://github.com/Libensemble/libensemble/tree/v0.9.3/libensemble/gen_funcs).

consensus/sims/alt_rosenbrock.py

@@ -0,0 +1,51 @@
+import numpy as np
+
+const = 1500
+
+
+def EvaluateFunction(x, component=np.nan):
+    """
+    Evaluates the chained Rosenbrock function
+    """
+
+    assert not np.isnan(component), "Must give a component"
+
+    i = component
+    x1 = x[i]
+    x2 = x[i + 1]
+    f = 100 * (x1**2 - x2) ** 2 + (x1 - 1) ** 2
+
+    return f
+
+
+def EvaluateJacobian(x, component=np.nan):
+    """
+    Evaluates the chained Rosenbrock Jacobian
+    """
+
+    df = np.zeros(len(x), dtype=float)
+
+    assert not np.isnan(component), "Must give a component"
+
+    i = component
+    x1 = x[i]
+    x2 = x[i + 1]
+
+    df[i] = 400 * x1 * (x1**2 - x2) + 2 * (x1 - 1)
+    df[i + 1] = -200 * (x1**2 - x2)
+
+    return 1.0 / const * df
+
+
+def alt_rosenbrock_eval(H, persis_info, sim_specs, _):
+    batch = len(H["x"])
+    H_o = np.zeros(batch, dtype=sim_specs["out"])
+
+    for i, x in enumerate(H["x"]):
+        obj_component = H["obj_component"][i]  # which f_i
+        if H[i]["get_grad"]:
+            H_o["gradf_i"][i] = EvaluateJacobian(x, obj_component)
+        else:
+            H_o["f_i"][i] = EvaluateFunction(x, obj_component)
+
+    return H_o, persis_info

consensus/sims/geomedian.py

@@ -0,0 +1,48 @@
+__all__ = ["geomedian_eval"]
+import numpy as np
+import numpy.linalg as la
+
+
+def EvaluateFunction(x, component, B):
+    """
+    Evaluates the sum of squares-variant chwirut function
+    """
+    m = B.shape[0]
+    assert B.shape[1] == len(x)
+    assert 0 <= component <= m - 1
+    i = component
+    b_i = B[i]
+
+    f_i = 1.0 / m * la.norm(x - b_i)
+    return f_i
+
+
+def EvaluateJacobian(x, component, B):
+    """
+    Evaluates the sum of squares-variant chwirut Jacobian
+    """
+    m = B.shape[0]
+    assert B.shape[1] == len(x)
+    assert 0 <= component <= m - 1
+    i = component
+    b_i = B[i]
+
+    df_i = 1.0 / m * (x - b_i) / la.norm(x - b_i)
+    return df_i
+
+
+def geomedian_eval(H, persis_info, sim_specs, _):
+    B = persis_info["params"]["B"]
+
+    num_xs = len(H["x"])  # b==1 always?
+    H_o = np.zeros(num_xs, dtype=sim_specs["out"])
+
+    for k, x in enumerate(H["x"]):
+        i = H[k]["obj_component"]  # f_i
+
+        if H[k]["get_grad"]:
+            H_o["gradf_i"][k] = EvaluateJacobian(x, i, B)
+        else:
+            H_o["f_i"][k] = EvaluateFunction(x, i, B)
+
+    return H_o, persis_info

consensus/sims/linear_regression.py

@@ -0,0 +1,70 @@
+import numpy as np
+
+
+def EvaluateFunction(theta, component, X, y, c, reg):
+    """
+    Evaluates linear regression with l2 regularization
+    """
+    i = component
+    m = len(y)
+
+    y_i = y[i]
+    X_i = X[:, i]
+    XT_theta = np.dot(X_i, theta)
+
+    f_i = (1 / m) * (np.dot(y_i, y_i) - 2 * y_i * XT_theta + XT_theta**2)
+
+    assert reg == "l2", "Only l2 regularization allowed"
+    # if reg is None:
+    #     reg_val = 0
+    # elif reg == 'l1':
+    #     reg_val = (c/m) * np.sum(np.abs(theta))
+    reg_val = (c / m) * np.dot(theta, theta)
+
+    return f_i + reg_val
+
+
+def EvaluateJacobian(theta, component, X, y, c, reg):
+    """
+    Evaluates linear regression with l2 regularization
+    """
+
+    i = component
+    m = len(y)
+
+    y_i = y[i]
+    X_i = X[:, i]
+
+    df_i = (2 / m) * (-y_i + np.dot(X_i, theta)) * X_i
+
+    assert reg == "l2", "Only l2 regularization allowed"
+
+    # if reg is None:
+    #     reg_val = 0
+    # elif reg == 'l1':
+    #     reg_val = (c/m) * np.sign(theta)
+    reg_val = (2 * c / m) * theta
+
+    return df_i + reg_val
+
+
+def linear_regression_eval(H, persis_info, sim_specs, _):
+    X = persis_info["params"]["X"]
+    y = persis_info["params"]["y"]
+    c = persis_info["params"]["c"]
+    reg = persis_info["params"].get("reg", None)
+
+    assert (reg is None) or (reg == "l1") or (reg == "l2"), f"Incompatible regularization {reg}"
+
+    batch = len(H["x"])
+    H_o = np.zeros(batch, dtype=sim_specs["out"])
+
+    for i, x in enumerate(H["x"]):
+        obj_component = H["obj_component"][i]  # which f_i
+
+        if H[i]["get_grad"]:
+            H_o["gradf_i"][i] = EvaluateJacobian(x, obj_component, X, y, c, reg)
+        else:
+            H_o["f_i"][i] = EvaluateFunction(x, obj_component, X, y, c, reg)
+
+    return H_o, persis_info

consensus/sims/logistic_regression.py

@@ -0,0 +1,72 @@
+import numpy as np
+
+
+def EvaluateFunction(theta, component, X, y, c, reg):
+    """
+    Evaluates linear regression with l2 regularization
+    """
+    i = component
+    m = len(y)
+
+    y_i = y[i]
+    X_i = X[:, i]
+
+    base = np.exp(-y_i * np.dot(X_i, theta))
+
+    f_i = (1 / m) * np.log(1 + base)
+
+    assert reg == "l2", "Only l2 regularization allowed"
+    # if reg is None:
+    #     reg_val = 0
+    # elif reg == 'l1':
+    #     reg_val = (c/m) * np.sum(np.abs(theta))
+    reg_val = (c / m) * np.dot(theta, theta)
+
+    return f_i + reg_val
+
+
+def EvaluateJacobian(theta, component, X, y, c, reg):
+    """
+    Evaluates linear regression with l2 regularization
+    """
+
+    i = component
+    m = len(y)
+
+    y_i = y[i]
+    X_i = X[:, i]
+
+    base = np.exp(-y_i * np.dot(X_i, theta))
+
+    df_i = (1 / m) * (-y_i * base) / (1 + base) * X_i
+
+    assert reg == "l2", "Only l2 regularization allowed"
+    # if reg is None:
+    #     reg_val = 0
+    # elif reg == 'l1':
+    #     reg_val = (c/m) * np.sign(theta)
+    reg_val = (2 * c / m) * theta
+
+    return df_i + reg_val
+
+
+def logistic_regression_eval(H, persis_info, sim_specs, _):
+    X = persis_info["params"]["X"]
+    y = persis_info["params"]["y"]
+    c = persis_info["params"]["c"]
+    reg = persis_info["params"].get("reg", None)
+
+    assert (reg is None) or (reg == "l1") or (reg == "l2"), f"Incompatible regularization {reg}"
+
+    batch = len(H["x"])
+    H_o = np.zeros(batch, dtype=sim_specs["out"])
+
+    for i, x in enumerate(H["x"]):
+        obj_component = H["obj_component"][i]  # which f_i
+
+        if H[i]["get_grad"]:
+            H_o["gradf_i"][i] = EvaluateJacobian(x, obj_component, X, y, c, reg)
+        else:
+            H_o["f_i"][i] = EvaluateFunction(x, obj_component, X, y, c, reg)
+
+    return H_o, persis_info

consensus/sims/nesterov_quadratic.py

@@ -0,0 +1,65 @@
+import numpy as np
+
+
+def EvaluateFunction(x, component):
+    """
+    Evaluates the chained Rosenbrock function
+    """
+    n = len(x)
+
+    i = component
+    assert 0 <= i <= n
+
+    if i == 0:
+        x_1 = x[i]
+        f_i = 0.5 * x_1**2 - x_1
+    elif i == n:
+        x_n = x[-1]
+        f_i = 0.5 * x_n**2
+    else:
+        x_1 = x[i - 1]
+        x_2 = x[i]
+        f_i = 0.5 * (x_2 - x_1) ** 2
+
+    return f_i
+
+
+def EvaluateJacobian(x, component=np.nan):
+    """
+    Evaluates the chained Rosenbrock Jacobian
+    """
+    n = len(x)
+    df = np.zeros(n, dtype=float)
+
+    i = component
+    assert 0 <= i <= n
+
+    if i == 0:
+        x_1 = x[i]
+        df[0] = x_1 - 1
+    elif i == n:
+        x_n = x[-1]
+        df[-1] = x_n
+    else:
+        x_1 = x[i - 1]
+        x_2 = x[i]
+
+        df[i - 1] = x_1 - x_2
+        df[i] = x_2 - x_1
+
+    return df
+
+
+def nesterov_quadratic_eval(H, persis_info, sim_specs, _):
+    batch = len(H["x"])
+    H_o = np.zeros(batch, dtype=sim_specs["out"])
+
+    for i, x in enumerate(H["x"]):
+        obj_component = H["obj_component"][i]  # which f_i
+
+        if H[i]["get_grad"]:
+            H_o["gradf_i"][i] = EvaluateJacobian(x, obj_component)
+        else:
+            H_o["f_i"][i] = EvaluateFunction(x, obj_component)
+
+    return H_o, persis_info

consensus/sims/svm.py

@@ -0,0 +1,77 @@
+import numpy as np
+
+
+def EvaluateFunction(theta, component, X, b, c, reg):
+    """
+    Evaluates svm with l1 regularization
+    """
+    i = component
+    m = len(b)
+
+    b_i = b[i]
+    X_i = X[:, i]
+
+    f_i = max(0, 1 - b_i * np.dot(X_i, theta))
+
+    assert reg == "l1", "Only l1 regularization allowed"
+
+    # if reg is None:
+    #     reg_val = 0
+    # elif reg == 'l1':
+    reg_val = (c / m) * np.sum(np.abs(theta))
+    # else:
+    #     reg_val = (c/m) * np.dot(theta, theta)
+
+    return f_i + reg_val
+
+
+def EvaluateJacobian(theta, component, X, b, c, reg):
+    """
+    Evaluates svm with l1 regularization
+    """
+
+    i = component
+    m = len(b)
+
+    b_i = b[i]
+    X_i = X[:, i]
+
+    score = b_i * np.dot(X_i, theta)
+
+    if score >= 1:
+        df_i = np.zeros(len(theta))
+    else:
+        df_i = -b_i * X_i
+
+    assert reg == "l1", "Only l1 regularization allowed"
+
+    # if reg is None:
+    #     reg_val = 0
+    # elif reg == 'l1':
+    reg_val = (c / m) * np.sign(theta)
+    # else:
+    #     reg_val = (2*c/m) * theta
+
+    return df_i + reg_val
+
+
+def svm_eval(H, persis_info, sim_specs, _):
+    X = persis_info["params"]["X"]
+    b = persis_info["params"]["b"]
+    c = persis_info["params"]["c"]
+    reg = persis_info["params"].get("reg", None)
+
+    assert (reg is None) or (reg == "l1") or (reg == "l2"), f"Incompatible regularization {reg}"
+
+    batch = len(H["x"])
+    H_o = np.zeros(batch, dtype=sim_specs["out"])
+
+    for i, x in enumerate(H["x"]):
+        obj_component = H["obj_component"][i]  # which f_i
+
+        if H[i]["get_grad"]:
+            H_o["gradf_i"][i] = EvaluateJacobian(x, obj_component, X, b, c, reg)
+        else:
+            H_o["f_i"][i] = EvaluateFunction(x, obj_component, X, b, c, reg)
+
+    return H_o, persis_info