convex_comb_agg_log -- add piecewise relation with a logarithmic number of binary variables
    using the convex combination formulation -- non-disaggregated.
    Parameters:
        - model: a model where to include the piecewise linear relation
        - a[k]: x-coordinate of the k-th point in the piecewise linear relation
        - b[k]: y-coordinate of the k-th point in the piecewise linear relation
    Returns the model with the piecewise linear relation on added variables X, Y, and z.

def convex_comb_agg_log(model,a,b):
    """convex_comb_agg_log -- add piecewise relation with a logarithmic number of binary variables
    using the convex combination formulation -- non-disaggregated.
    Parameters:
        - model: a model where to include the piecewise linear relation
        - a[k]: x-coordinate of the k-th point in the piecewise linear relation
        - b[k]: y-coordinate of the k-th point in the piecewise linear relation
    Returns the model with the piecewise linear relation on added variables X, Y, and z.
    """
    K = len(a)-1
    G = int(math.ceil((math.log(K)/math.log(2))))     # number of required bits
    w,g = {},{}
    for k in range(K+1):
        w[k] = model.addVar(lb=0, ub=1, vtype="C")
    for j in range(G):
        g[j] = model.addVar(vtype="B")
    X = model.addVar(lb=a[0], ub=a[K], vtype="C")
    Y = model.addVar(lb=-model.infinity(), vtype="C")

    model.addCons(X == quicksum(a[k]*w[k]  for k in range(K+1)))
    model.addCons(Y == quicksum(b[k]*w[k]  for k in range(K+1)))
    model.addCons(quicksum(w[k] for k in range(K+1)) == 1)

    # binary variables setup
    for j in range(G):
        zeros,ones = [0],[]
        # print(j,"\tinit zeros:",zeros,"ones:",ones
        for k in range(1,K+1):
            # print(j,k,"\t>zeros:",zeros,"ones:",ones
            if (1 & gray(k)>>j) == 1 and (1 & gray(k-1)>>j) == 1:
                ones.append(k)
            if (1 & gray(k)>>j) == 0 and (1 & gray(k-1)>>j) == 0:
                zeros.append(k)
            # print(j,k,"\tzeros>:",zeros,"ones:",ones

        # print(j,"\tzeros:",zeros,"ones:",ones
        model.addCons(quicksum(w[k] for k in ones) <= g[j])
        model.addCons(quicksum(w[k] for k in zeros) <= 1-g[j])

    return X,Y,w