Skip to content

optimizer_utils — Foundation Utilities

Shared utilities used by all optimizers in QuantumExpectedValueFunctionProject/. Provides VariableRegister (qubit-encoded integer with operator circuits) and PowerSystem_1Bus (problem specification container).


VariableRegister

QuantumExpectedValueFunctionProject.optimizer_utils.VariableRegister

VariableRegister(max_value, encoding='binary')

class VariableRegister: Use this object to store an integer variable in a register of qubits, and to perform interactions between registers of qubits. We use a register of qubits 's', which represent binary variables 'y' to encode the integer variable 'v'. Depending on the encoding, this integer variable 'v' is created by a linear combination of binary variables 'sum b*y', and 'b' depends on the type of encoding.

Source code in QuantumExpectedValueFunctionProject/optimizer_utils.py
def __init__(self, max_value, encoding='binary'):
    try:
        assert(encoding in ['binary', 'unary'])
    except:
        print("ERROR: Encoding type must be 'binary' or 'unary', got '{}'".format(encoding))
    self.encoding = encoding
    self.max_value = max_value
    if self.max_value > 0:
        self.width = int(np.log2(self.max_value)) + 1 if self.encoding == 'binary' else self.max_value
    elif self.max_value == 0:
        self.width = 0
    else:
        print("ERROR: max_value must be non-negative")

    if 2**self.width - 1 > self.max_value and self.encoding == 'binary':
        print("WARNING: Possible variable range ({}) is larger than variable maximum value ({})".format(2**self.width-1, self.max_value))

Functions

__unaryLessThanOperator

__unaryLessThanOperator(other, amplitude)

__unaryLessThanOperator If this > other, return e^i(this - other). Otherwise, Identity

Source code in QuantumExpectedValueFunctionProject/optimizer_utils.py
def __unaryLessThanOperator(self, other, amplitude):
    ''' __unaryLessThanOperator
        If this > other, return e^i(this - other). Otherwise, Identity
    '''
    # TODO: implement with a real circuit
    qc = QuantumCircuit(self.width + other.width)
    length = self.width+other.width
    dim = 2**(length)
    op = np.zeros(shape=(dim,dim), dtype=np.complex128)
    for i in range(dim):
        bstr = ''.join(list(int_to_bstr_N(i, length)))
        my_bstr = bstr[:self.width]
        other_bstr = bstr[self.width:]
        my_val = self.getValue(my_bstr)
        other_val = self.getValue(other_bstr)
        if my_val > other_val:
            op[i,i] = np.exp(1j*amplitude*(my_val-other_val)) 
        else:
            op[i,i] = 1.
    qc.append(Operator(op), range(length))
    return qc

__binaryLessThanOperator

__binaryLessThanOperator(other, amplitude)

__binaryLessThanOperator If this > other, return e^i(this - other). Otherwise, Identity

Source code in QuantumExpectedValueFunctionProject/optimizer_utils.py
def __binaryLessThanOperator(self, other, amplitude):
    ''' __binaryLessThanOperator
        If this > other, return e^i(this - other). Otherwise, Identity
    '''
    # TODO: implement with a real circuit
    qc = QuantumCircuit(self.width + other.width)
    length = self.width+other.width
    dim = 2**(length)
    op = np.zeros(shape=(dim,dim), dtype=np.complex128)
    for i in range(dim):
        bstr = ''.join(list(int_to_bstr_N(i, length)))
        my_bstr = bstr[:self.width]
        other_bstr = bstr[self.width:]
        my_val = self.getValue(my_bstr[::-1])
        other_val = other.getValue(other_bstr[::-1])
        if my_val > other_val:
            op[i,i] = np.exp(1j*amplitude*(my_val-other_val))
        else:
            op[i,i] = 1.
    qc.append(Operator(op), range(length))
    return qc

__binaryLessThanValue

__binaryLessThanValue(value, amplitude)

__binaryLessThanValue

Source code in QuantumExpectedValueFunctionProject/optimizer_utils.py
def __binaryLessThanValue(self, value, amplitude): 
    ''' __binaryLessThanValue
    '''
    # TODO: implement with a real circuit
    qc = QuantumCircuit(self.width)
    dim = 2**self.width
    op = np.zeros(shape=(dim,dim), dtype=np.complex128)
    for i in range(dim):
        bstr = ''.join(list(int_to_bstr_N(i, self.width)))
        my_val = self.getValue(bstr[::-1])
        if my_val > value:
            op[i,i] = np.exp(1j*amplitude*(my_val-value))
        else:
            op[i,i] = 1.
    qc.append(Operator(op), range(self.width))
    return qc

__binaryNumberOperator

__binaryNumberOperator(amplitude)

__binaryNumberOperator Return amplitude * v in Little Endian

Source code in QuantumExpectedValueFunctionProject/optimizer_utils.py
def __binaryNumberOperator(self, amplitude):
    ''' __binaryNumberOperator
        Return amplitude * v in Little Endian
    '''
    qc = QuantumCircuit(self.width)
    for i in range(self.width):
        qc.p(amplitude * 2**i, i)
    return qc

__unaryNumberOperator

__unaryNumberOperator(amplitude)

__unaryNumberOperator Return amplitude * v encoded as a unary operator

Source code in QuantumExpectedValueFunctionProject/optimizer_utils.py
def __unaryNumberOperator(self, amplitude):
    ''' __unaryNumberOperator
        Return amplitude * v encoded as a unary operator
    '''
    qc = QuantumCircuit(self.width)
    for i in range(self.width):
        qc.p(amplitude, i)
    return qc

__unarySwapOperator

__unarySwapOperator(other, amplitude)

__unarySwapOperator Return integer encoding as a unary operator For now, the swap interactions are all-to-all, single Trotter step

Source code in QuantumExpectedValueFunctionProject/optimizer_utils.py
def __unarySwapOperator(self, other, amplitude):
    ''' __unarySwapOperator
        Return integer encoding as a unary operator
        For now, the swap interactions are all-to-all, single Trotter step
    ''' 
    qc = QuantumCircuit(self.width+other.width, name='swap({})'.format(np.round(amplitude)))
    for i in range(self.width):
        for j in range(self.width, self.width+other.width):
            qc.append(SwapGate().power(amplitude), [i,j])
    return qc

__binarySwapOperator

__binarySwapOperator(other, amplitude)

__binarySwapOperator Return integer swap operation for the binary encoding

Source code in QuantumExpectedValueFunctionProject/optimizer_utils.py
def __binarySwapOperator(self, other, amplitude):
    ''' __binarySwapOperator
        Return integer swap operation for the binary encoding
    '''
    ''' kind of works
    length = self.width+other.width
    qc = QuantumCircuit(length)
    # TODO: find a real circuit
    dim = 2**length
    ham = np.zeros(shape=(dim,dim), dtype=np.complex128)
    for i in range(dim):
        # decode first number
        bstr_i = int_to_bstr_N(i, length)[::-1]
        my_int_i = self.getValue(bstr_i[:self.width])
        other_int_i = other.getValue(bstr_i[self.width:])
        for j in range(dim):
            # decode second number
            bstr_j = int_to_bstr_N(j, length)[::-1]
            my_int_j = self.getValue(bstr_j[:self.width])
            other_int_j = other.getValue(bstr_j[self.width:])
            if my_int_i + other_int_i == my_int_j + other_int_j and i!=j:
                ham[i,j] = np.abs(my_int_i - my_int_j)*self.max_value#/self.max_value #1.
    #ham[0,0] = 1.
    #ham[dim-1,dim-1] = 1.
    #ham /= sp.linalg.norm(ham)
    op = sp.linalg.expm(-1j *np.pi/2* amplitude * ham)
    qc.append(Operator(op), list(range(length)))
    '''
    length = self.width+other.width
    ancilla = length
    qc = QuantumCircuit(length)
    # TODO: account for different register width, remove ancilla
    assert(self.width == other.width)
    for i in range(self.width):
        j = i+self.width
        #qc.append(SwapGate().power(amplitude), [i,j]) # swap bits of equal value
        # promote two bits of equal value - manual matrix, bit order (i+1),i,j
        op = [[1,0,0,0,0,0,0,0],#111
              [0,1,0,0,0,0,0,0],#110
              [0,0,1,0,0,0,0,0],#101
              [0,0,0,np.cos(1/2*amplitude*np.pi/2),1j*np.sin(1/2*amplitude*np.pi/2),0,0,0],#100 - 011
              [0,0,0,1j*np.sin(1/2*amplitude*np.pi/2),np.cos(1/2*amplitude*np.pi/2),0,0,0],#011 - 100
              [0,0,0,0,0,1,0,0],#010
              [0,0,0,0,0,0,1,0],#001
              [0,0,0,0,0,0,0,1]]#000
        if i < self.width-1:
            qc.append(Operator(op), [i,j,i+1])
            qc.append(Operator(op), [i,j,j+1])
        # swap bits of equal significance
        qc.append(SwapGate().power(amplitude), [i,j]) # swap bits of equal value
    return qc

__unaryProductOperator

__unaryProductOperator(other, amplitude)

__unaryProductOperator

Source code in QuantumExpectedValueFunctionProject/optimizer_utils.py
def __unaryProductOperator(self, other, amplitude):
    ''' __unaryProductOperator
    '''
    qc = QuantumCircuit(self.width + other.width)
    for j in range(self.width):
        for k in range(self.width, self.width + other.width):
            qc.cp(amplitude, j, k)
    return qc

__binaryProductOperator

__binaryProductOperator(other, amplitude)

__binaryProductOperator

Source code in QuantumExpectedValueFunctionProject/optimizer_utils.py
def __binaryProductOperator(self, other, amplitude):
    ''' __binaryProductOperator
    '''
    qc = QuantumCircuit(self.width + other.width)
    for j in range(self.width):
        for k in range(other.width):
            qc.cp(2**j * 2**k * amplitude, j, k+self.width)
    return qc

__binarySquaredOperator

__binarySquaredOperator(amplitude)

__binarySquaredOperator

Source code in QuantumExpectedValueFunctionProject/optimizer_utils.py
def __binarySquaredOperator(self, amplitude):
    ''' __binarySquaredOperator
    '''
    qc = QuantumCircuit(self.width)
    for j in range(self.width):
        for k in range(j+1,self.width):
            qc.cp(2*2**j * 2**k * amplitude, j, k)
    return qc

__unarySquaredOperator

__unarySquaredOperator(amplitude)

__unarySquaredOperator

Source code in QuantumExpectedValueFunctionProject/optimizer_utils.py
def __unarySquaredOperator(self, amplitude):
    ''' __unarySquaredOperator
    '''
    qc = QuantumCircuit(self.width)
    for j in range(self.width):
        for k in range(j+1,self.width):
            qc.cp(2*amplitude, j, k)
    return qc

swapOperator

swapOperator(other, amplitude)

swapOperator Return a circuit which does an amplitude * "integer swap" between two registers.

Source code in QuantumExpectedValueFunctionProject/optimizer_utils.py
def swapOperator(self, other, amplitude):
    ''' swapOperator
        Return a circuit which does an amplitude * "integer swap" between two registers.
    '''
    if self.encoding == 'binary':
        if other.encoding != 'binary':
            print("ERROR: I have encoding '{}' and do not support SWAP with a register of encoding '{}'".format(self.encoding, other.encoding))
            exit(1)
        return self.__binarySwapOperator(other, amplitude)
    elif self.encoding == 'unary':
        if other.encoding != 'unary':
            print("ERROR: I have encoding '{}' and do not support SWAP with a register of encoding '{}'".format(self.encoding, other.encoding))
            exit(1)
        return self.__unarySwapOperator(other, amplitude)
    else:
        print("ERROR: Unimplemented")
        exit(1)

numberOperator

numberOperator(amplitude)

numberOperator Return a circuit which does amplitude*v using the corresponding encoding.

Source code in QuantumExpectedValueFunctionProject/optimizer_utils.py
def numberOperator(self, amplitude):
    ''' numberOperator
        Return a circuit which does amplitude*v using the corresponding encoding.
    '''
    if self.encoding == 'binary':
        return self.__binaryNumberOperator(amplitude)
    elif self.encoding == 'unary':
        return self.__unaryNumberOperator(amplitude)
    else:
        print("ERROR: Unimplemented")
        exit(1)

productOperator

productOperator(other, amplitude)

productOperator Return a circuit which, given two registers, computes amplitudethisother

Source code in QuantumExpectedValueFunctionProject/optimizer_utils.py
def productOperator(self, other, amplitude):
    ''' productOperator
        Return a circuit which, given two registers, computes amplitude*this*other
    '''
    if self.encoding == 'binary':
        if other.encoding != 'binary':
            print("ERROR: I have encoding '{}' and do not support N_i*N_j with a register of encoding '{}'".format(self.encoding, other.encoding))
            exit(1)
        return self.__binaryProductOperator(other, amplitude)
    elif self.encoding == 'unary':
        if other.encoding != 'unary':
            print("ERROR: I have encoding '{}' and do not support N_i*N_j with a register of encoding '{}'".format(self.encoding, other.encoding))
            exit(1)
        return self.__unaryProductOperator(other, amplitude)
    else:
        print("ERROR: Unimplemented")
        exit(1)

squaredOperator

squaredOperator(amplitude)

squaredOperator Return a circuit which does amplitudethisthis

Source code in QuantumExpectedValueFunctionProject/optimizer_utils.py
def squaredOperator(self, amplitude):
    ''' squaredOperator
        Return a circuit which does amplitude*this*this
    '''
    if self.encoding == 'binary':
        return self.__binarySquaredOperator(amplitude)
    elif self.encoding == 'unary':
        return self.__unarySquaredOperator(amplitude)
    else:
        print("ERROR: Unimplemented")
        exit(1)

lessThanOperator

lessThanOperator(other, amplitude)

lessThanOperator Return a circuit which does phase amplitude*(this-other) if this > other, otherwise 1

Source code in QuantumExpectedValueFunctionProject/optimizer_utils.py
def lessThanOperator(self, other, amplitude):
    ''' lessThanOperator
        Return a circuit which does phase amplitude*(this-other) if this > other, otherwise 1
    '''
    if self.encoding == 'binary':
        if other.encoding != 'binary':
            print("ERROR: I have encoding '{}' and do not support LTPhase with a register of encoding '{}'".format(self.encoding, other.encoding))
            exit(1)
        return self.__binaryLessThanOperator(other, amplitude)
    elif self.encoding == 'unary':
        if other.encoding != 'unary':
            print("ERROR: I have encoding '{}' and do not support LTPhase with a register of encoding '{}'".format(self.encoding, other.encoding))
            exit(1)
        return self.__unaryLessThanOperator(other, amplitude)
    else:
        print("ERROR: Unimplemented")
        exit(1)

oracleLessThanOperator

oracleLessThanOperator(other, amplitude)

Return a circuit which rotates an ancilla by (this-other) if this > other

Source code in QuantumExpectedValueFunctionProject/optimizer_utils.py
def oracleLessThanOperator(self, other, amplitude):
    ''' 
        Return a circuit which rotates an ancilla by (this-other) if this > other
    '''
    if self.encoding == 'binary':
        print("ERROR: oracleLessThanOperator not implemented for encoding 'binary'")
        exit(1)
    elif self.encoding == 'unary':
        if other.encoding != 'unary':
            print("ERROR: I have encoding {} but other has encoding {}".format('unary', 'binary'))
        return self.__unaryOracleLessThanOperator(other, amplitude)
    else:
        print("ERROR: Unimplemented")
        exit(1)

lessThanValue

lessThanValue(value, amplitude)

lessThanValue Return a circuit which does phase amplitude*(this-value) if this > value, otherwise Identity

Source code in QuantumExpectedValueFunctionProject/optimizer_utils.py
def lessThanValue(self, value, amplitude):
    ''' lessThanValue
        Return a circuit which does phase amplitude*(this-value) if this > value, otherwise Identity
    '''
    if self.encoding == 'binary':
        return self.__binaryLessThanValue(value, amplitude)
    elif self.encoding == 'unary':
        print("ERROR, unimplemented lessThanValue for unary encoding")
        # TODO do we need to implement this?
    else:
        print("ERROR: Unimplemented")
        exit(1)

getValue

getValue(bstr)

Given a register's bitstring, return the integer it encodes

Source code in QuantumExpectedValueFunctionProject/optimizer_utils.py
def getValue(self, bstr):
    ''' Given a register's bitstring, return the integer it encodes
    '''
    assert(len(bstr) == self.width)
    if self.encoding == 'binary':
        v = 0
        revbstr = bstr[::-1] # reverse the bitstring to do Little Endian encoding
        for i in range(self.width):
            v += int(revbstr[i])*2**i
        return v
    elif self.encoding == 'unary':
        return sum([int(b) for b in bstr])
    else:
        print("ERROR: Unimplemented")
        exit(1)

setValue

setValue(integer)

Given an integer, return a circuit which creates that integer

Source code in QuantumExpectedValueFunctionProject/optimizer_utils.py
def setValue(self, integer):
    ''' Given an integer, return a circuit which creates that integer
    '''
    assert(integer <= self.max_value)
    qc = QuantumCircuit(self.width)
    if self.encoding == 'unary':
        for i in range(integer):
            qc.x(i)
    elif self.encoding == 'binary':
        bstr = int_to_bstr_N(integer, self.width)
        for i,val in enumerate(bstr):
            if val == '1':
                qc.x(i)
    return qc    

PowerSystem_1Bus

QuantumExpectedValueFunctionProject.optimizer_utils.PowerSystem_1Bus

PowerSystem_1Bus(gas_costs, wind_costs, decision_levels, undersatisfied_cost, demand, pdf, normalization=None)

All-to-all lossless single-bus power system specification.

Stores all parameters of the two-stage stochastic unit commitment problem and provides utilities for normalizing costs, evaluating the classical solution, and post-processing quantum measurement results.

The single-bus approximation ignores transmission losses, so every generator contributes directly to meeting total demand \(d\).

Parameters:

Name Type Description Default
gas_costs

List of per-unit costs for gas generators.

required
wind_costs

List of per-unit operational costs for wind turbines.

required
decision_levels

Number of discrete levels per variable (integer range 0 to decision_levels-1).

required
undersatisfied_cost

Penalty cost per unit of unmet demand (recourse cost \(c_r\)).

required
demand

Total demand \(d\) that must be met.

required
pdf

Dict mapping scenario tuples \((\xi_0, \xi_1, \ldots)\) to probabilities.

required
normalization

Optional (max_cost, max_amplitude) pair for scaling.

None
Note

Discretization constraint: cost * (decision_levels - 1) == max_cost.

Initialize the power system with costs, discretization, and PDF.

Source code in QuantumExpectedValueFunctionProject/optimizer_utils.py
def __init__(self, gas_costs, wind_costs, decision_levels, undersatisfied_cost, demand, pdf, normalization=None):
    '''Initialize the power system with costs, discretization, and PDF.'''
    assert(np.isclose(sum(pdf.values()),1.))

    # number of each generator type and cost
    self.num_gas_generators = len(gas_costs)
    self.num_wind_turbines = len(wind_costs)
    self.gas_costs = gas_costs
    self.wind_costs = wind_costs
    self.undersatisfied_cost = undersatisfied_cost

    self.normalization = normalization
    if normalization is not None: # sometimes, we want to normalize the operator amplitudes; this parameter is something like (max cost, max op amplitude)
        obj = normalization[0]
        su2 = normalization[1]
        for i in range(len(self.gas_costs)):
            self.gas_costs[i] *= su2/obj
        for i in range(len(self.wind_costs)):
            self.wind_costs[i] *= su2/obj
        self.undersatisfied_cost *= su2/obj

    # discretization
    # TODO: make this variable specific
    self.decision_levels = decision_levels
    # demand
    self.demand = demand

    # probability distriution function, or scenarios
    self.pdf = pdf
    self.scenarios = list(self.pdf.keys()) # index is scenario 'number'

    # decompose pdf into samples
    '''
    self.sample_list = []
    self.sample_hist = {}
    denominators = [Fraction(Decimal(str(self.pdf[scene]))).denominator for scene in self.scenarios] # iterate over 'scenarios' to preserve order of vars
    # get least common multiple of the Pr's denominators
    lcm = 1 
    for denom in denominators:
        lcm = lcm*denom//gcd(lcm,denom)
    # convert Pr to frequency of occurency in SAA
    for scene, prob in self.pdf.items():
        reps = prob * lcm
        self.sample_hist[scene] = int(reps)
        for _ in range(int(np.around(reps))):
            self.sample_list.append(scene)
    '''
    #print(self.sample_hist)
    # expanded costs
    self.variable_costs = copy.deepcopy(self.gas_costs)
    for scenario, amp in self.pdf.items():
        [self.variable_costs.append(amp*w) for w in self.wind_costs]
        self.variable_costs.append(amp*self.undersatisfied_cost)

Attributes

scenarios instance-attribute

scenarios = list(keys())

self.sample_list = [] self.sample_hist = {} denominators = [Fraction(Decimal(str(self.pdf[scene]))).denominator for scene in self.scenarios] # iterate over 'scenarios' to preserve order of vars

get least common multiple of the Pr's denominators

lcm = 1 for denom in denominators: lcm = lcm*denom//gcd(lcm,denom)

convert Pr to frequency of occurency in SAA

for scene, prob in self.pdf.items(): reps = prob * lcm self.sample_hist[scene] = int(reps) for _ in range(int(np.around(reps))): self.sample_list.append(scene)

Functions

getFirstStageCosts

getFirstStageCosts(gas_decisions)

Given a list of decisions over the gas variables, return the cost and wind decision variables for each

Source code in QuantumExpectedValueFunctionProject/optimizer_utils.py
def getFirstStageCosts(self, gas_decisions):
    ''' Given a list of decisions over the gas variables, return the cost and wind decision variables for each
    '''
    decision_to_cost = {}
    decision_to_decision = {}
    for gas_dec in gas_decisions:
        result = self.cobylaSolve(gas_values=gas_dec)
        if result.fun is not None:
            decision_to_cost[gas_dec] = result.fun
        decision_to_decision[gas_dec] = result.x
    return decision_to_cost,decision_to_decision

cobylaSolve

cobylaSolve(discrete=True, gas_values=None)

Use cobyla to solve the optimization problem. Either solve the entire two-stage problem via direct expansion or, given a gas decision, solve only the second stage

Source code in QuantumExpectedValueFunctionProject/optimizer_utils.py
def cobylaSolve(self, discrete=True, gas_values=None):
    ''' Use cobyla to solve the optimization problem. 
        Either solve the entire two-stage problem via direct expansion or, 
        given a gas decision, solve _only_ the second stage 
    '''
    num_vars = self.num_gas_generators + len(self.scenarios)*(self.num_wind_turbines + 1)
    gas_upper_bounds = [self.decision_levels - 1] * self.num_gas_generators

    # Get the second stage upper bounds
    second_stage_upper_bounds = []
    for scenario_id, scenario in enumerate(self.scenarios):
        for wind in scenario:
            second_stage_upper_bounds.append(wind)
        second_stage_upper_bounds.append(self.demand)
    if gas_values is None:
        bounds = Bounds([0.]*num_vars, gas_upper_bounds + second_stage_upper_bounds)
    else:
        assert(len(gas_values) == self.num_gas_generators)
        gas_values = list(gas_values)
        other_zeros = [0.]*(num_vars - self.num_gas_generators)
        bounds = Bounds(gas_values + other_zeros, gas_values + second_stage_upper_bounds)

    # Set the linear constraint
    A = []
    for scenario_id, scenario in enumerate(self.scenarios):
        constraint_row = [1]*self.num_gas_generators + [0]*scenario_id*(self.num_wind_turbines + 1)
        for wind in scenario:
            constraint_row.append(1)
        constraint_row.append(1)
        constraint_row += [0]*(num_vars - len(constraint_row))
        A.append(constraint_row)
    lin_constraint = LinearConstraint(A, self.demand, self.demand)

    # QUICK AND DIRTY force non-anticipativity among wind
    #B = []
    #for w in range(self.num_wind_turbines):
    #    for scenario_i, scenario in enumerate(self.scenarios):
    #        for scenario_j, scenario in enumerate(self.scenarios):
    #            if scenario_j <= scenario_i:
    #                continue
    #            constraint_row = [0]*self.num_gas_generators + [0]*scenario_i*(self.num_wind_turbines + 1)
    #            # si
    #            constraint_row += [0]*w + [1] + [0]*(self.num_wind_turbines-w)
    #            constraint_row += [0]*(scenario_j-scenario_i-1)*(self.num_wind_turbines + 1)
    #            # sj
    #            constraint_row += [0]*w + [-1] + [0]*(self.num_wind_turbines-w)
    #            constraint_row += [0]*(num_vars - len(constraint_row))
    #            #print(num_vars, len(constraint_row), constraint_row)
    #            B.append(constraint_row)
    #nonant_constraint = LinearConstraint(B, 0, 0)


    # Optimize
    # If we use discrete variables 
    if discrete:
        #res = milp(self.variable_costs, constraints=[lin_constraint, nonant_constraint], bounds=bounds, integrality=1)
        res = milp(self.variable_costs, constraints=[lin_constraint,], bounds=bounds, integrality=1)
    else:
        # TODO: initial condition 
        x0 = [0] * num_vars
        x0[0] = self.demand
        res = minimize(self._price, x0, constraints=[lin_constraint], bounds=bounds)

    return res

plotMeasurementsVExpectedCost

plotMeasurementsVExpectedCost(measured_decisions, pdf=False)

plotMeasurementsVExpectedCost Given a dictionary of gas generator measured results, plot the measurement count vs. expected cost of that decision It is the quantum optimizer's job to create a histogram of gas decisions.

Source code in QuantumExpectedValueFunctionProject/optimizer_utils.py
def plotMeasurementsVExpectedCost(self, measured_decisions, pdf=False):
    ''' plotMeasurementsVExpectedCost 
        Given a dictionary of gas generator measured results, plot the measurement count vs. expected cost of that decision
        It is the quantum optimizer's job to create a histogram of gas decisions.
    '''
    #plt.cla()
    first_index = 2 if not pdf else 3
    ex_costs,_ = self.getFirstStageCosts(measured_decisions.keys())

    #plt.figure()
    plt.figure().set_figheight(10)

    plt.subplot(first_index, 1, 1)
    measured_decisions = {k:measured_decisions[k] for k in sorted(measured_decisions.keys())}
    plt.bar([str(a) for a in measured_decisions.keys()], measured_decisions.values())
    plt.ylabel(r"Pr$(x)$")
    plt.xlabel(r"$x$")

    plt.subplot(first_index, 1, 2)
    ex_costs = {k:ex_costs[k] for k in sorted(ex_costs.keys())}
    if self.normalization is not None:
        #plt.bar([str(a) for a in val_of_decisions.keys()], [system.normalization[0]/system.normalization[1]*v for v in val_of_decisions.values()])
        #plt.bar([str(a) for a in ex_costs.keys()], [v for v in ex_costs.values()])
        plt.bar([str(a) for a in ex_costs.keys()], [self.unNormalize(v) for v in ex_costs.values()])
    else:
        plt.bar([str(a) for a in ex_costs.keys()], ex_costs.values())
    plt.ylabel(r'$\sum_g c_gx_g + \mathbb{E}_\xi [L(x,\xi)]$')
    plt.xlabel(r"$x$")

    if pdf:
        plt.subplot(3, 1, 3)
        pdf = {k:self.pdf[k] for k in sorted(self.pdf.keys())}
        plt.bar([str(a) for a in pdf.keys()], pdf.values())
        plt.ylabel(r"Pr$(\xi)$")
        plt.xlabel(r"$\xi$")
    plt.tight_layout()
    return plt