Mixed scheme implementation

[Giannis1993]

This is to discuss how to implement/assemble mixed schemes. I am a fan of dictionaries, so here is how I implemented it. I remember however that you guys had some other thoughts with a list/tuple of subproblems.

---

In the main function I would write (if n=4):

def test_mixed_square(n):

    # Instantiate problem
    prob = Square(n)

    # Define variable and response sets of a mixed approximation scheme as dictionaries
    var_set = {0: np.array([0]),
               1: np.array([1]),
               2: np.arange(2, prob.n)}
    resp_set = {0: np.array([0]),
                1: np.array([1])}

    # Instantiate subproblem objects for a mixed approximation scheme
    subprob_map =  {(0, 0): Subproblem(intervening=MMA(prob.xmin[var_set[0]], prob.xmax[var_set[0]]),
                                       approximation=Taylor1(),
                                       ml=MoveLimitIntervening(xmin=prob.xmin[var_set[0]],
                                                               xmax=prob.xmax[var_set[0]])),
                    (0, 1): Subproblem(intervening=MMA(prob.xmin[var_set[1]], prob.xmax[var_set[1]]),
                                       approximation=Taylor1(),
                                       ml=MoveLimitIntervening(xmin=prob.xmin[var_set[1]],
                                                               xmax=prob.xmax[var_set[1]])),
                    (0, 2): Subproblem(intervening=MMA(prob.xmin[var_set[2]], prob.xmax[var_set[2]]),
                                       approximation=Taylor1(),
                                       ml=MoveLimitIntervening(xmin=prob.xmin[var_set[2]],
                                                               xmax=prob.xmax[var_set[2]])),
                    (1, 0): Subproblem(intervening=Linear(),
                                       approximation=Taylor1(),
                                       ml=MoveLimitIntervening(xmin=prob.xmin[var_set[0]],
                                                               xmax=prob.xmax[var_set[0]])),
                    (1, 1): Subproblem(intervening=Linear(),
                                       approximation=Taylor1(),
                                       ml=MoveLimitIntervening(xmin=prob.xmin[var_set[1]],
                                                               xmax=prob.xmax[var_set[1]])),
                    (1, 2): Subproblem(intervening=Linear(),
                                       approximation=Taylor1(),
                                       ml=MoveLimitIntervening(xmin=prob.xmin[var_set[2]],
                                                               xmax=prob.xmax[var_set[2]]))}

    # Instantiate a mixed scheme
    subprob = Mixed(subprob_map, var_set, resp_set)

    # Instantiate solver
    solver = SvanbergIP(prob.n, prob.m)

    # Initialize iteration counter and design
    itte = 0
    x_k = prob.x0.copy()

    # Optimization loop
    while not converged:

        # Evaluate responses and sensitivities at current point, i.e. g(X^(k)), dg(X^(k))
        f = prob.g(x_k)
        df = prob.dg(x_k)

        # Print current iteration and x_k
        print('iter: {:^4d}  |  x: {:<20s}  |  obj: {:^9.3f}  |  constr: {:^6.3f}'.format(
            itte, np.array2string(x_k[0:2]), f[0], f[1]))

        # Build approximate sub-problem at X^(k)
        subprob.build(x_k, f, df)

        # Call solver (x_k, g and dg are within approx instance)
        x, y, z, lam, xsi, eta, mu, zet, s = solver.subsolv(subprob)
        x_k = x.copy()

        itte += 1

---

Whereas the Mixed class would be as:

from sao.problems.subproblem import Subproblem
import numpy as np

## Mixed scheme assembly class
class Mixed(Subproblem):

    def __init__(self, subprob_map, var_set, resp_set):                 # no need to run the parent constructor
        self.num_of_var_sets = len(var_set.keys())                      # number of variable sets
        self.num_of_resp_sets = len(resp_set.keys())                    # number of response sets

        # Calculate n, m
        self.n = 0
        self.m = -1
        for l in range(0, self.num_of_var_sets):
            self.n += len(var_set[l])
        for p in range(0, self.num_of_resp_sets):
            self.m += len(resp_set[p])

        # Initialize arrays
        self.alpha = np.zeros(self.n, dtype=float)
        self.beta = np.ones(self.n, dtype=float)
        self.f = np.zeros(self.m + 1, dtype=float)
        self.df = np.zeros((self.m + 1, self.n), dtype=float)
        self.ddf = np.zeros((self.m + 1, self.n), dtype=float)

        # Store dictionaries
        self.subprob_map = subprob_map                                  # dictionary of subproblems
        self.var_set = var_set                                          # dictionary of variable sets
        self.resp_set = resp_set                                        # dictionary of response sets

    def build(self, x, f, df, ddf=None):
        self.n, self.m = len(x), len(f) - 1             # might be unnecessary
        self.f, self.df, self.ddf = f, df, ddf

       # Build constituent subproblems
        for p in range(0, self.num_of_resp_sets):
            for l in range(0, self.num_of_var_sets):
                if ddf is not None:
                    self.subprob_map[p, l].build(x[self.var_set[l]],
                                                 f[self.resp_set[p]],
                                                 df[np.ix_(self.resp_set[p], self.var_set[l])],
                                                 ddf[np.ix_(self.resp_set[p], self.var_set[l])])
                else:
                    self.subprob_map[p, l].build(x[self.var_set[l]],
                                                 f[self.resp_set[p]],
                                                 df[np.ix_(self.resp_set[p], self.var_set[l])])

        # Get mixed subproblem bounds
        self.get_bounds()

    # This is to find the most conservative bounds for each variable
    def get_bounds(self):
        self.alpha[:] = -1.
        self.beta[:] = 1.
        for p in range(0, self.num_of_resp_sets):
            for l in range(0, self.num_of_var_sets):
                self.alpha[self.var_set[l]] = np.maximum(self.alpha[self.var_set[l]],
                                                         self.subprob_map[p, l].alpha)
                self.beta[self.var_set[l]] = np.minimum(self.beta[self.var_set[l]],
                                                        self.subprob_map[p, l].beta)

    # Assembly of response vector for a mixed scheme
    def g(self, x):
        g_value = - (self.num_of_var_sets - 1) * self.f      # cuz each time a var_set is added, so is the const. term
        for p in range(0, self.num_of_resp_sets):
            for l in range(0, self.num_of_var_sets):
                g_value[self.resp_set[p]] += self.subprob_map[p, l].approx.g(
                    self.subprob_map[p, l].inter.y(x[self.var_set[l]]).T)
        return g_value

    # Assembly of sensitivity array for a mixed scheme
    def dg(self, x):
        dg_value = np.zeros_like(self.df)
        for p in range(0, self.num_of_resp_sets):
            for l in range(0, self.num_of_var_sets):
                dg_value[[self.resp_set[p]], [self.var_set[l]]] = self.subprob_map[p, l].approx.dg(
                    self.subprob_map[p, l].inter.y(x[self.var_set[l]]).T,
                    self.subprob_map[p, l].inter.dydx(x[self.var_set[l]]))
        return dg_value

    # Assembly of 2nd-order sensitivity array for a mixed scheme
    def ddg(self, x):
        ddg_value = np.zeros_like(self.df)
        for p in range(0, self.num_of_resp_sets):
            for l in range(0, self.num_of_var_sets):
                ddg_value[[self.resp_set[p]], [self.var_set[l]]] = self.subprob_map[p, l].approx.ddg(
                    self.subprob_map[p, l].inter.y(x[self.var_set[l]]).T,
                    self.subprob_map[p, l].inter.dydx(x[self.var_set[l]]),
                    self.subprob_map[p, l].inter.ddyddx(x[self.var_set[l]]))
        return ddg_value

[MaxvdKolk]

I have not tried to run the code, so just some comments based on reading these functions:

• It would be great to have a default option available, or some other way to specify the application of subproblems to a bunch of combinations directly. For instance, in the example the combinations (1,0), (1,1), (1,2) are all Linear & Taylor1. Maybe we can wrap this subproblem map in a (very small) class where we can pass a such a map and a default instance. Inside that class, we can then call self.map.get(key, self.default) as fallback to the default.

• I am wondering if there is a way to avoid the duplication of code for accessing the individual variables, notice how get_bounds, g, dg, ddg all perform the same nested for loop across the combinations of variable and response sets. There are of course many ways to resolve that, we could store all indices on class initialisation as the outer product of var_set.keys() and resp_set.keys() Note this does not necessarily have to be a prebuild matrix, that can also be done by some iterator, you could think of

  def problems(self):
    for p in range(...):
        for l in range(...):
            yield (p, l)
  
  # inside functions
  for (p, l) in self.problems:
    # apply operation to each problem 

• Do we want to have var_set inside, or outside? What I mean is prob.xmin[var_set[1]] versus problem_set[index].xmin? Maybe if we flip the current structure inside out, we have less nesting of the indexing for the problems?

I would have to take a bit more in depth look though, to really suggest other implementations. Maybe it's good to discuss the current structure and possible alternatives this Friday? I doubt I can take more detailed look before, but perhaps in the weekend I can

[aatmdelissen]

My suggestion for the initialisation is as follows

# Default initialize all to MMA, Taylor1
mix = Mixed(var_sets, resp_sets, default_intervening=MMA, default_approximation=Taylor1, default_movelimit=...)

# Specify a few sets
mix.set_group(0, 1, intervening=ConLin, approximation=Taylor2, ...)
mix.set_group(1, range(5), approximation=Taylor2)
mix.set_group(1, 2, intervening=Linear)

In this way a default uniform mixed scheme is initialized, and by using the set_group method, specific groups are altered. Intervening variables, approximation, and move limit should be able to set independently, e.g.:

# One variable set, three response sets
mix = Mixed(var_sets, resp_sets, default_intervening=MMA, default_approximation=Taylor1, default_movelimit=...)
mix.set_group(0, [0, 1], intervening=ConLin)
mix.set_group(0, 1, approximation=Taylor2)

Which results in (0,0 - ConLin, Taylor1) and (0,1 - ConLin, Taylor2) and (0,2 - MMA, Taylor1).

[Giannis1993]

Let's continue this discussion on #52