Skip to content

Instantly share code, notes, and snippets.

@amitjamadagni

amitjamadagni/pricing_benchmark.py

Created September 3, 2016 17:58
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save amitjamadagni/c253be215746fcf5c919f27b923fef87 to your computer and use it in GitHub Desktop.
Save amitjamadagni/c253be215746fcf5c919f27b923fef87 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
pricing_benchmark.py
import numpy as np
import random as rand
from numpy.linalg import inv
from numpy.linalg import norm
from scipy.optimize import minimize
from scipy.optimize import linprog
import warnings
import time
import sys
from memory_profiler import profile
sys.setrecursionlimit(100000)
# The initial inputs should be the
# 1. products higher bound on the prices - array - p_h
# 2. products lower bound on the prices - array - p_l
# 3. L_1 function which is either (t*log(t))^(1/2), c*t^(2/3)
# 4. Choose `n` linearly independent vectors each vector lying between p_l < p < p_h
class Pricing:
"""
Pricing is the starting object for all the computations. The inputs for creating
an Pricing object is the following :
1. n_time : The number of time steps, given `n` time steps, we compute the optimal price for `n+1`.
- Type : array
2. k_products : The number of products.
- Type : array
3. demand_mat : At each time step, each product has a particular demand, with time steps as columns
and products as rows, construct the demand matrix.
- Type : matrix
4. p_h : For each product, there is a higher limit on the prices for each product.
- Type : array
5. p_l : For each product, there is a lower limit on the prices for each product.
- Type : array
"""
def __init__(self, n_time, k_products, demand_mat, p_h, p_l):
self._ntime = n_time
self._kproducts = k_products
self._demand_mat = demand_mat
self._ph = p_h
self._pl = p_l
# generating linearly independent price vectors
def independent_price_vectors(self):
"""
At each instant for the `n` time steps, given there are `k` products. For each time instant, we randomly
generate a price, say `p` which lies between the highest price and the lowest price, for every of the `k` products.
Therefore, at every time instant we generate a column price vector. Stacking column vectors together we end up
with a `k x n` matrix. Each vector is an independent linear vector randomly generated.
"""
mat = np.ones(shape=(self._kproducts, self._ntime))
for product in range(1, self._kproducts):
for time in range(self._ntime):
mat[product][time] = rand.uniform(self._pl[product], self._ph[product])
return mat
# Helper functions for normal distribution, refer notes for the equation.
def P_helper_func_normal(self, price_vectors_mat, k):
"""
Helper function to compute the outer product of the above generated column stack.
Also price vector column used for computation of beta_matrix as mentioned in the notes.
"""
coeff_mat = np.zeros(shape=price_vectors_mat.shape)
const_col = np.zeros(shape=(1, self._ntime))
for ncol in range(self._ntime):
col = price_vectors_mat[:,ncol]
coeff_mat = coeff_mat + np.outer(col, col)
const_col = const_col + col*(self._demand_mat[k, ncol])
return coeff_mat, const_col
# Helper function for poisson distribution, refer notes for the equation.
def P_helper_func_poisson(self, price_vectors_mat, k):
"""
Helper function to compute the outer product of the above generated column stack.
Also price vector column used for computation of beta_matrix as mentioned in the notes.
"""
coeff_mat = np.zeros(shape=price_vectors_mat.shape)
const_col = np.zeros(shape=(1, self._ntime))
for ncol in range(self._ntime):
col = price_vectors_mat[:,ncol]
coeff_mat = coeff_mat + np.outer(col, col)
const_col = const_col + col*(self._demand_mat[k, ncol]-1)
return coeff_mat, const_col
# given a product, computing the beta column vector using the above
# generated linearly independent price vectors for normal distribution
def beta_k_normal(self, price_vectors_mat, k):
"""
Generating beta column got by solving A*X = B equation, for how we arrived at this
refer the notes.
"""
P_help = self.P_helper_func_normal(price_vectors_mat, k)
P_t_outersum = P_help[0]
P_t_sum = P_help[1]
beta = np.linalg.solve(np.array(P_t_outersum), np.reshape(P_t_sum, (self._ntime, 1)))
return beta
# given a product, computing the beta column vector using the above
# generated linearly independent price vectors for poisson distribution
def beta_k_poisson(self, price_vectors_mat, k):
"""
Generating beta column got by solving A*X = B equation, for how we arrived at this
refer the notes.
"""
P_help = self.P_helper_func_poisson(price_vectors_mat, k)
P_t_outersum = P_help[0]
P_t_sum = P_help[1]
beta = np.linalg.solve(np.array(P_t_outersum), np.reshape(P_t_sum, (self._ntime, 1)))
return beta
def beta_mat_normal(self, price_vectors_mat):
"""
Stacking above columns to generate the beta matrix.
"""
x = []
for product in range(self._kproducts):
x.append(self.beta_k_normal(price_vectors_mat, product))
return np.column_stack(x)
def beta_mat_poisson(self, price_vectors_mat):
"""
Stacking above columns to generate the beta matrix.
"""
x = []
for product in range(self._kproducts):
x.append(self.beta_k_poisson(price_vectors_mat, product))
return np.column_stack(x)
###################################################################################################
# The functions defined from here are specific to Normal Distribution. For any other distribution
# the following functions are to be redefined :
#
# 1. Optimization objective function - For normal distribution it is `optimal_function`
# defined below
# 2. Optimization function determinant - For normal distribution it is `optimal_func_der`
# defined below
# 3. Optimization function hessian mat - For normal distribution it is `optimal_func_hess`
# defined below
#
# Reference on how determinant and hessian matrix can be developed from optimization objective
# function :
#
# 1. Documentation from SciPy -
# http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html#scipy.optimize.minimize
#
# 2. Development tutorial usage in SciPy -
# http://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html
###################################################################################################
def optimal_function_normal(p, beta_m):
row_sum = 0
for k in range(1, len(p)):
row_sum = row_sum - p[k]*np.dot(p, beta_m[:,k])
return row_sum
def optimal_func_der_normal(p, beta_m):
der = np.zeros_like(p)
for k in range(len(p)):
s = 0
for i in range(1, len(p)):
s = s - p[i]*(beta_m[k, i] + beta_m[i, k])
der[k] = s
# der[k] = p[k]*(beta_m[:,k].sum() + beta_m[k,:].sum())
# der[k] = -np.dot(p, beta_m[:,k]) - p[k]*beta_m[k, k]
return der
def optimal_func_hess_normal(p, beta_m):
H = np.zeros(shape=beta_m.shape)
for i in range(1, beta_m.shape[0]):
for j in range(1, beta_m.shape[1]):
H[i, j] = -beta_m[i, j] - beta_m[j, i]
# dia = [2 for i in range(beta_m.shape[0])]
# np.fill_diagonal(H, dia)
return H
# def optimal_func_hess_normal(p, beta_m):
# H = np.ones(shape=beta_m.shape)
# for i in range(beta_m.shape[0]):
# for j in range(beta_m.shape[1]):
# H[i, j] = -beta_m[i, j] - beta_m[j, i]
# dia = [beta_m[i,i] for i in range(beta_m.shape[0])]
# np.fill_diagonal(H, dia)
# return H
###################################################################################################
# The functions defined from here are specific to Poisson Distribution. For any other distribution
# the following functions are to be redefined :
#
# 1. Optimization objective function - For normal distribution it is `optimal_function`
# defined below
# 2. Optimization function determinant - For normal distribution it is `optimal_func_der`
# defined below
# 3. Optimization function hessian mat - For normal distribution it is `optimal_func_hess`
# defined below
#
# Reference on how determinant and hessian matrix can be developed from optimization objective
# function :
#
# 1. Documentation from SciPy -
# http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html#scipy.optimize.minimize
#
# 2. Development tutorial usage in SciPy -
# http://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html
###################################################################################################
def optimal_function_poisson(p, beta_m):
row_sum = 0
for k in range(1, len(p)):
row_sum = row_sum - p[k]*np.exp(np.dot(p, beta_m[:,k]))
return row_sum
def optimal_func_der_poisson(p, beta_m):
der = np.zeros_like(p)
for k in range(len(p)):
prod = 0
for m in range(1, beta_m.shape[0]):
exp_sum = 0
for i in range(1, beta_m.shape[0]):
exp_sum = exp_sum + p[i]*beta_m[m, i]
prod = prod - p[m]*beta_m[m, k]*np.exp(exp_sum)
exp_sum_2 = 0
for j in range(1, beta_m.shape[0]):
exp_sum_2 = exp_sum_2 + p[j]*beta_m[k, j]
der[k] = prod - np.exp(exp_sum_2)
return der
def optimal_func_hess_poisson(p, beta_m):
H = np.zeros(shape=beta_m.shape)
for k in range(1, beta_m.shape[0]):
for m in range(1, beta_m.shape[0]):
prod = 0
for i in range(1, beta_m.shape[0]):
exp_sum = 0
for j in range(1, beta_m.shape[0]):
exp_sum = exp_sum + p[i]*beta_m[i, j]
prod = prod - p[i]*beta_m[i, k]*beta_m[i, m]*np.exp(exp_sum)
exp_sum_2 = 0
exp_sum_3 = 0
for i1 in range(1, beta_m.shape[0]):
exp_sum_2 = exp_sum_2 + p[i1]*beta_m[k, i1]
exp_sum_2 = beta_m[k, m]*np.exp(exp_sum_2)
exp_sum_3 = beta_m[m, k]*np.exp(exp_sum_2)
H[k, m] = prod - exp_sum_2 - exp_sum_3
return H
####################################################################################################
# Given the above functions are defined for a particular distribution, the above functions are to be
# passed to the `optimal_solution` as parameters to obtain the optimal solution.
####################################################################################################
def optimal_solution_normal(pricing_object, opt_func, opt_func_der, opt_func_hess, initial_guess, beta_mat, bound_tuples):
res = minimize(opt_func, initial_guess, args=(beta_mat), method='Newton-CG', jac=opt_func_der, hess=opt_func_hess, bounds=bound_tuples)
for bound_index in range(1, len(bound_tuples)):
if bound_tuples[bound_index][0] > res.x[bound_index] or res.x[bound_index] > bound_tuples[bound_index][1] or np.isnan(res.x[bound_index]) or res.x[0] < 0.0 :
pv = pricing_object.independent_price_vectors()
beta_mat = pricing_object.beta_mat_normal(pv)
return optimal_solution_normal(pricing_object, opt_func, opt_func_der, opt_func_hess, initial_guess, beta_mat, bound_tuples)
else:
return res.x
def optimal_solution_poisson(pricing_object, opt_func, opt_func_der, opt_func_hess, initial_guess, beta_mat, bound_tuples):
res = minimize(opt_func, initial_guess, args=(beta_mat), method='Newton-CG', jac=opt_func_der, hess=opt_func_hess, bounds=bound_tuples)
print res.x
for bound_index in range(1, len(bound_tuples)):
if bound_tuples[bound_index][0] > res.x[bound_index] or res.x[bound_index] > bound_tuples[bound_index][1] or np.isnan(res.x[bound_index]): # or res.x[0] > 2.0 or res.x[0] < 0.9:
pv = pricing_object.independent_price_vectors()
beta_mat = pricing_object.beta_mat_poisson(pv)
return optimal_solution_poisson(pricing_object, opt_func, opt_func_der, opt_func_hess, initial_guess, beta_mat, bound_tuples)
else:
return res.x
#####################################################################################################
# Depending on different cases referring to 4.10, 4.11 the below functions are to be used which are
# general in nature and can be extended to any distribution. The bounds are passed in manually
# depending on function L(t).
# 1. optimal_sol_const_1 aims at providing a solution to 4.10
# 2. optimal_sol_const_2 aims at providing a solution to 4.11
# TODO : Detect the bounds automatically and automate the cases to directly obtain the optimal price.
# As of now the check on the cases are manual.
#####################################################################################################
def optimal_sol_const_1_normal(pricing_object, opt_func, opt_func_der, opt_func_hess, initial_guess, beta_mat, bound_tuples, price_vectors_mat, bound):
opt_sol = optimal_solution_normal(pricing_object, opt_func, opt_func_der, opt_func_hess, initial_guess, beta_mat, bound_tuples)
P_t_mat = pricing_object.P_helper_func(price_vectors_mat, 1)[0]
if np.trace(inv(P_t_mat + np.dot(np.transpose(opt_sol), opt_sol)))**-1 < bound:
return optimal_sol_const_1(pricing_object, opt_func, opt_func_der, opt_func_hess, initial_guess, beta_mat, bound_tuples, price_vectors_mat, bound)
else:
return opt_sol
def optimal_sol_const_1_poisson(pricing_object, opt_func, opt_func_der, opt_func_hess, initial_guess, beta_mat, bound_tuples, price_vectors_mat, bound):
opt_sol = optimal_solution_poisson(pricing_object, opt_func, opt_func_der, opt_func_hess, initial_guess, beta_mat, bound_tuples)
P_t_mat = pricing_object.P_helper_func(price_vectors_mat, 1)[0]
if np.trace(inv(P_t_mat + np.dot(np.transpose(opt_sol), opt_sol)))**-1 < bound:
return optimal_sol_const_1(pricing_object, opt_func, opt_func_der, opt_func_hess, initial_guess, beta_mat, bound_tuples, price_vectors_mat, bound)
else:
return opt_sol
def optimal_sol_const_2_normal(pricing_object, opt_func, opt_func_der, opt_func_hess, initial_guess, beta_mat, bound_tuples, price_vectors_mat, bound):
opt_sol = optimal_solution_normal(pricing_object, opt_func, opt_func_der, opt_func_hess, initial_guess, beta_mat, bound_tuples)
P_t_mat = pricing_object.P_helper_func(price_vectors_mat, 1)[0]
if norm(np.dot(inv(P_t_mat), np.transpose(opt_sol)))**2/(1 + np.dot(np.dot(opt_sol, inv(P_t_mat)), np.transpose(opt_sol))) < bound:
return optimal_sol_const_2(pricing_object, opt_func, opt_func_der, opt_func_hess, initial_guess, beta_mat, bound_tuples, price_vectors_mat, bound)
else:
return opt_sol
def optimal_sol_const_2_poisson(pricing_object, opt_func, opt_func_der, opt_func_hess, initial_guess, beta_mat, bound_tuples, price_vectors_mat, bound):
opt_sol = optimal_solution_poisson(pricing_object, opt_func, opt_func_der, opt_func_hess, initial_guess, beta_mat, bound_tuples)
P_t_mat = pricing_object.P_helper_func(price_vectors_mat, 1)[0]
if norm(np.dot(inv(P_t_mat), np.transpose(opt_sol)))**2/(1 + np.dot(np.dot(opt_sol, inv(P_t_mat)), np.transpose(opt_sol))) < bound:
return optimal_sol_const_2(pricing_object, opt_func, opt_func_der, opt_func_hess, initial_guess, beta_mat, bound_tuples, price_vectors_mat, bound)
else:
return opt_sol
########################################################################################################
# CUI for various inputs
# number of products, time steps
# example of pricing object
# x = Pricing(4, 4, np.matrix([[11., 22., 31., 41.], [33., 20., 23., 43.], [23., 18., 28., 45.], [28., 16., 26., 47.]]), [1., 50., 50., 50.], [1., 10., 10., 10.])
# np.matrix([[11., 22., 31., 41., 22.], [33., 20., 23., 43., 23.], [23., 18., 28., 45., 28.], [28., 16., 26., 47., 29.], [22., 19., 30., 48., 30.]])
# [1., 50., 50., 50., 50.]
# [1., 10., 10., 10., 10.]
# initial_guess = [10., 23., 39., 24., 28.]
@profile
def func():
n_products = input("Enter the number of products (that is add +1 to the actual number) : ")
time_steps = input("Enter the number of time steps (should be equal to the above number of products) : ")
demand_mat = input("Enter the demand matrix with time instance as columns and products as rows : ")
max_price = input("Enter the maximum price for each product as an array : ")
min_price = input("Enter the minimum price for each product as an array : ")
print "Creating the pricing object "
pricing_obj = Pricing(time_steps, n_products, demand_mat, max_price, min_price)
print "Generating the independent vectors once for all the cases to follow"
pv = pricing_obj.independent_price_vectors()
# generating the respective beta matrix
distribution = raw_input("Enter the distribution, the supported distributions are :\n1. Normal (N) \n2. Poisson (P)\n ")
print "Generating the beta matrix for the above independent vectors using the above distribution :",distribution
if distribution == "N":
print "Choosen distribution is Normal"
beta_matrix = pricing_obj.beta_mat_normal(pv)
elif distribution == "P":
print "Choosen distribution is Poisson"
beta_matrix = pricing_obj.beta_mat_poisson(pv)
else:
print "Error : choose N for Normal, P for Poisson"
bound_tuple = tuple((pricing_obj._pl[i], pricing_obj._ph[i]) for i in range(len(pricing_obj._ph)))
initial_guess = input("Enter an initial guess for the optimal solution generation : ")
# generating different optimal solutions depending on various conditions
# TODO : this is to be automated to generate the final solution, as the bounds are unknown this
# that is a partial blocker or may be a full blocker.
# y = optimal_solution(x, optimal_function, optimal_func_der, optimal_func_hess, [23, 15.5, 36.4, 22.3], beta_matrix, bound_tuple)
print "Generating the optimal solution to be compared with the bound in order to predict the pricing at future times ..."
warnings.filterwarnings("ignore")
if distribution == "N":
optimal_set = []
t_end = time.time() + 60 * 1
while time.time() < t_end:
y = optimal_solution_normal(pricing_obj, optimal_function_normal, optimal_func_der_normal, optimal_func_hess_normal, initial_guess, beta_matrix, bound_tuple)
optimal_set.append(y)
least = np.argsort([opt[0] for opt in optimal_set])
opt = []
for ind in least:
opt.append(optimal_set[ind])
for i in range(len(opt)):
if np.array_equal(opt[i], np.array(initial_guess)):
continue
else:
print " optimal solution for normal distribution: ", opt[i]
break
else:
print " optimal solution for normal distribution: ", opt[0]
elif distribution == "P":
y = optimal_solution_poisson(pricing_obj, optimal_function_poisson, optimal_func_der_poisson, optimal_func_hess_poisson, initial_guess, beta_matrix, bound_tuple)
print " optimal solution for poisson distribution: ", y
else:
print "Error : choose N for Normal, P for Poisson"
# opt_sol = optimal_sol_const_2(x, optimal_function, optimal_func_der, optimal_func_hess, [23, 15.5, 36.4, 22.3], beta_matrix, bound_tuple, pv, 228.6)
# print opt_sol
if __name__ == '__main__':
func()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment