pricing_new_parallel.py

import numpy as np
import random as rand
import multiprocessing
import logging
import warnings
import os
import time
import sys
from numpy.linalg import inv
from numpy.linalg import norm
from scipy.optimize import minimize
from scipy.optimize import linprog
from multiprocessing import Queue, Process, cpu_count
sys.setrecursionlimit(100000000)
# 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

    def P_helper_func(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

    # given a product, computing the beta column vector using the above
    # generated linearly independent price vectors
    def beta_k(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(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(self, price_vectors_mat):
        """
        Stacking above columns to generate the beta matrix.
        """
        x = []
        for product in range(self._kproducts):
            x.append(self.beta_k(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(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(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(beta_m.shape[0]):
    #     for j in range(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 -1*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(pricing_object, opt_func, opt_func_der, opt_func_hess, initial_guess, beta_mat, bound_tuples, queue):
    t_end = time.time() + 60 * 1
    while time.time() < t_end:
        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] > 1.75 or res.x[0] < 0.9:
                pv = pricing_object.independent_price_vectors()
                beta_mat = pricing_object.beta_mat(pv)
                return optimal_solution(pricing_object, opt_func, opt_func_der, opt_func_hess, initial_guess, beta_mat, bound_tuples, queue)
        else:
            # print res.x
            # return res.x
            queue.put(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(pricing_object, opt_func, opt_func_der, opt_func_hess, initial_guess, beta_mat, bound_tuples, price_vectors_mat, bound):
    opt_sol = optimal_solution(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(pricing_object, opt_func, opt_func_der, opt_func_hess, initial_guess, beta_mat, bound_tuples, price_vectors_mat, bound):
    opt_sol = optimal_solution(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.5, 22.5, 31.2, 41.2], [33.1, 45.4, 23.1, 34.3], [15.2, 56.2, 35.4, 41.4], [26.2, 41.4 , 36.4, 41.8]]), [1., 50., 50., 50.],4)
# initial_guess = [23, 15.5, 36.4, 22.3]
########################################################################################################

# 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)

# Reference : http://stackoverflow.com/questions/28530705/how-to-use-pythons-multiprocessing-module?noredirect=1&lq=1

if __name__ == '__main__':
    print "number of cores : ", cpu_count()
    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()

    print "Generating the beta matrix for the above independent vectors"
    beta_matrix = pricing_obj.beta_mat(pv)

    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 : ")
    cpu_counts = input("Enter the number of processes : ")
    queue = Queue()
    processes = []
    # multiprocessing.log_to_stderr(logging.DEBUG) # uncomment this for multi processing logs.
    for i in range(cpu_counts):
        warnings.filterwarnings("ignore")
        processes.append(Process(target=optimal_solution, args=(pricing_obj, optimal_function_normal, optimal_func_der_normal, optimal_func_hess_normal, initial_guess, beta_matrix, bound_tuple, queue,)))
        processes[-1].start()
        results = []
    for i in range(cpu_counts):
        results.append(queue.get())
    for process in processes:
        process.join()
    print results
    print np.mean(results, axis=0)
    # final_result = []
    # for result in results:
    #     if 0.95 < result[0] < 1.05 :
    #         final_result.append(result)
    # print final_result
    # print np.mean(final_result, axis=0)


 #np.matrix([[1.52, 5.62, 3.54, 4.14, 1.23 , 1.46 , 1.78 ], [2.62, 4.14 , 2.66, 5.18, 1.58 , 1.86 , 2.27] , [ 2.34 , 3.42 , 3.44 , 5.68 , 1.96 , 2.89 , 1.64], [2.62, 2.85 , 3.86, 5.06 , 2.86 , 4.25 , 2.45] , [2.45 , 3.56 , 4.32 , 4.24 , 3.67 , 3.87, 2.56] , [2.76 , 3.78 ,4.56, 4.05 , 4.06 , 3.56 , 3.64], [ 2.89 , 4.34 , 4.23 ,3.86 , 4.46, 3.45 ,3.89] ])

#Enter the maximum price for each product as an array : [1., 300., 300., 300. ,300. ,300., 300.]
#Enter the minimum price for each product as an array : [1., 100., 100., 100., 150. ,150. ,130. ]
#Creating the pricing object
#Generating the independent vectors once for all the cases to follow
#Generating the beta matrix for the above independent vectors
#Enter an initial guess for the optimal solution generation : [130., 230.4, 280.6, 280.9 , 270.6 , 290.3 , 240.1]
# [1., 230.4, 280.6, 280.9 , 270.6 , 290.3 , 240.1]
# what to do when nan is being returned
#in case the solution is not changing, what is to be done
#emergency situations, what to do with the price
# to calculate the initial prices , we need to use a particular day , and put in the demands as per the time stamps , with the initial guess as the
# price at which the restaurant is selling and then keep updating on the prices based on the output pricing
# [ 24.3 , 110.3  , 254.8 , 164.4 ,  169.6 , 204.1 , 257.9]