#####################################################
# Conjoint Analysis                                 #
# Chapter 12 - Chestnut Ridge Smart Watch Example   #
#####################################################

# This script reproduces the full Chestnut Ridge smart watch conjoint analysis example in Python:
#   Step 1. Conjoint study design (full and fractional factorial designs)
#   Step 2. Estimate the part-worths for each respondent
#   Step 3. Part-worths, attribute importance, willingness-to-pay, and market simulations
# Run it from top to bottom. Each step prints its result or shows a plot.
# The "# %%" lines split the script into cells you can run one at a time
# in VS Code, Spyder, or PyCharm.

# %%
## Install Packages (if needed - only installs packages that are missing)
## Each entry is {name used with import: name used with pip install}
import sys, subprocess, importlib.util
packages = {"pandas": "pandas", "numpy": "numpy", "statsmodels": "statsmodels", "matplotlib": "matplotlib"}
new_packages = [pip_name for import_name, pip_name in packages.items() if importlib.util.find_spec(import_name) is None]
if new_packages: subprocess.check_call([sys.executable, "-m", "pip", "install", *new_packages])

## Load Packages and Set Seed
import itertools
from tkinter import Tk, filedialog

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.api as sm

rng = np.random.default_rng(1)

## Show all columns when printing tables
pd.set_option("display.max_columns", None)
pd.set_option("display.width", 200)


## Function to choose a file with a dialog box
## (the Python version of file.choose() in R - no need to set the working directory)
def file_choose(title="Choose a file"):
    root = Tk()
    root.withdraw()                     # hide the empty Tk window
    root.attributes("-topmost", True)   # bring the dialog to the front
    path = filedialog.askopenfilename(title=title, filetypes=[("CSV files", "*.csv"), ("All files", "*.*")])
    root.destroy()
    if not path:
        raise SystemExit("No file was chosen.")
    return path


# %%
#####################################################
# Step 1: Conjoint Study Design Setup               #
#####################################################

## Set up attributes and levels as a dictionary
## The first level of each attribute is the baseline level
attrib_level = {"brand": ["CR", "Apple", "Samsung", "FitBit"],
                "ship": ["$0", "$10", "$20"],
                "restock": ["0%", "5%", "10%", "15%"],
                "retdays": ["7 days", "14 days", "21 days"],
                "price": ["$150", "$200", "$250", "$300"]}


## Function to set the attributes as categories with the levels in the order above
def set_levels(df):
    df = df.copy()
    for a, levels in attrib_level.items():
        df[a] = pd.Categorical(df[a], categories=levels)
    return df


## Function to create dummy variables (one per non-baseline level) plus an intercept
def dummy_code(df):
    X = pd.get_dummies(set_levels(df), drop_first=True, dtype=float)
    return sm.add_constant(X)


## Create the full factorial design (576 product bundles)
experiment = pd.DataFrame(list(itertools.product(*attrib_level.values())),
                          columns=list(attrib_level.keys()))
print(len(experiment))

## Create the fractional factorial design (30 product profiles)
## A D-optimal design picks the 30 profiles from the full factorial that give the most
## precise regression estimates. Start from 30 random profiles and keep swapping in
## profiles from the full factorial as long as the swap improves the design.
## Note: This will NOT produce the same 30 profiles as the R version. R's caFactorialDesign()
## and this Python code use different random number generators and different search steps,
## so they start from different random profiles and end at a different design.
## The correlations in the Python design are therefore also different from R's
## (e.g., brand and shipping correlate at about 0.13 here vs. 0.01 in R).
## Setting the same seed (1) in both languages does not make them match.
X_cand = dummy_code(experiment).to_numpy()


def d_criterion(rows):
    sign, logdet = np.linalg.slogdet(X_cand[rows].T @ X_cand[rows])
    return logdet if sign > 0 else -np.inf


cards = 30
rows = rng.choice(len(X_cand), size=cards, replace=False)
current = d_criterion(rows)
improved = True
while improved:
    improved = False
    for i in range(cards):
        best, best_j = current, None
        for j in range(len(X_cand)):
            if j in rows:
                continue
            trial = rows.copy()
            trial[i] = j
            value = d_criterion(trial)
            if value > best + 1e-9:
                best, best_j = value, j
        if best_j is not None:
            rows[i] = best_j
            current = best
            improved = True
design = experiment.iloc[np.sort(rows)]
print(design)


## Function to check for correlation in a fractional factorial design
## Each attribute is coded 1, 2, 3, ... by level
def design_cor(df):
    return set_levels(df).apply(lambda x: x.cat.codes + 1).corr().round(3)


print(design_cor(design))

## Run the conjoint analysis study

## Note: Python and R do not produce the same fractional factorial design.
## The survey in this example was run using the design created in R,
## so read in that design for the rest of the analysis.
design = pd.read_csv(file_choose("Choose conjoint_profiles.csv"))  ## Choose the file named conjoint_profiles.csv
print(design_cor(design))


# %%
#####################################################
# Step 2: Estimate the Part-Worths                  #
#####################################################

## Read in the survey preference results
pref = pd.read_csv(file_choose("Choose conjoint_preferences.csv"))  ## Choose the file named conjoint_preferences.csv

## Set up the product profiles as dummy variables
## Base Case: Brand CR, Shipping $0, Restock 0%, Retdays 7 days, Price $150
X = dummy_code(design)
print(X)

## Estimate the part-worths for each respondent using OLS regression
part_worths = []
for r in pref.columns:
    fit = sm.OLS(pref[r], X).fit()
    temp = {"intercept": fit.params["const"]}
    for a, levels in attrib_level.items():
        for level in levels:
            ## Part-worths are 0 for each baseline level
            temp[level] = 0.0 if level == levels[0] else fit.params[f"{a}_{level}"]
    part_worths.append(temp)
part_worths = pd.DataFrame(part_worths, index=pref.columns).round(3)
print(part_worths)


# %%
#####################################################
# Step 3: Part-Worth Analysis, Willingness-to-Pay,  #
#         and Market Simulations                    #
#####################################################

## Function to create a horizontal bar chart with labels
## Bars are shown in the order they are listed
def plot_bars(values, title, labels=None):
    if labels is None:
        labels = values.round(3).astype(str)
    fig, ax = plt.subplots(figsize=(7, 5))
    bars = ax.barh(values.index, values, color="steelblue")
    ax.bar_label(bars, labels=list(labels), padding=3)
    ax.axvline(0, color="black", linewidth=1)
    ax.invert_yaxis()
    ax.margins(x=0.2)
    ax.set_title(title, loc="left")
    plt.tight_layout()
    plt.show()


## Average Attribute Part-Worths
## Exclude the baseline levels since their part-worths are always 0
baseline = [levels[0] for levels in attrib_level.values()]
avg_pw = part_worths.mean().drop(baseline)
print(avg_pw.round(3))
plot_bars(avg_pw, "Average Attribute Part-Worths")

## Average Attribute Importance
## Importance is the range (max - min) of the part-worths within each attribute
importance = pd.DataFrame({a: part_worths[levels].max(axis=1) - part_worths[levels].min(axis=1)
                           for a, levels in attrib_level.items()})
importance.columns = ["Brand", "Shipping", "Restock", "Return Days", "Price"]
avg_importance = importance.mean()
print(avg_importance.round(3))
plot_bars(avg_importance, "Average Attribute Importance")

## Percentage Average Attribute Importance
pct_importance = avg_importance / avg_importance.sum()
pct_importance.index = pct_importance.index + " %"
print((pct_importance * 100).round(2))
plot_bars(pct_importance, "Percentage Average Attribute Importance",
          labels=[f"{v:.2%}" for v in pct_importance])

## Average Willingness to Pay for a Feature
## The change in price from $150 to $300 ($150) tells us how many dollars one util is worth
avg_all = part_worths.mean()
dollars_per_util = 150 / (avg_all["$150"] - avg_all["$300"])
features = ["Apple", "Samsung", "FitBit", "$10", "$20", "5%", "10%", "15%", "14 days", "21 days"]
wtp = avg_all[features] * dollars_per_util
wtp.index = ["WTP - Brand Apple", "WTP - Brand Samsung", "WTP - Brand FitBit",
             "WTP - Shipping $10", "WTP - Shipping $20",
             "WTP - Restock 5%", "WTP - Restock 10%", "WTP - Restock 15%",
             "WTP - Return Days 14", "WTP - Return Days 21"]
print(wtp.round(2))
plot_bars(wtp, "Average Willingness to Pay for a Feature",
          labels=[f"-${abs(v):.2f}" if v < 0 else f"${v:.2f}" for v in wtp])


## Function to calculate each respondent's utility for a set of products
## Each row of products is one product profile
def calc_util(products):
    util = pd.DataFrame(index=part_worths.index)
    for _, p in products.iterrows():
        u = part_worths["intercept"].copy()
        for a in attrib_level:
            u = u + part_worths[p[a]]
        util[p["product"]] = u
    return util


## Function to estimate market shares using three choice rules
def calc_shares(util):
    ## Maximum utility (first choice) rule: each respondent chooses the product with the highest utility
    first_choice = util.eq(util.max(axis=1), axis=0).astype(float)
    ## Share of preference rule: choice probability is proportional to utility
    ## Negative utilities are set to 0 so that no product receives a negative probability
    util_pos = util.clip(lower=0)
    share_pref = util_pos.div(util_pos.sum(axis=1), axis=0)
    ## Logit share of preference rule: choice probability is proportional to exp(utility)
    logit_share = np.exp(util).div(np.exp(util).sum(axis=1), axis=0)
    return pd.DataFrame({"First Choice": first_choice.mean(),
                         "Share of Preference": share_pref.mean(),
                         "Logit Share": logit_share.mean()})


## Function to plot market shares by choice rule
def plot_shares(shares, title):
    fig, axes = plt.subplots(3, 1, figsize=(7, 7), sharex=True)
    for ax, rule in zip(axes, shares.columns):
        bars = ax.barh(shares.index, shares[rule], color="steelblue")
        ax.bar_label(bars, labels=[f"{v:.2%}" for v in shares[rule]], padding=3)
        ax.invert_yaxis()
        ax.set_title(rule, fontsize=10)
    axes[-1].set_xlim(0, shares.to_numpy().max() * 1.2)
    axes[-1].xaxis.set_major_formatter(lambda x, pos: f"{x:.0%}")
    fig.suptitle(title, x=0.02, ha="left")
    plt.tight_layout()
    plt.show()


## Market Share - Current Products
current = pd.DataFrame({"product": ["Apple", "Samsung", "FitBit"],
                        "brand": ["Apple", "Samsung", "FitBit"],
                        "ship": ["$0", "$20", "$10"],
                        "restock": ["15%", "0%", "10%"],
                        "retdays": ["7 days", "14 days", "14 days"],
                        "price": ["$200", "$300", "$250"]})
util_current = calc_util(current)
print(util_current.round(3))
shares_current = calc_shares(util_current)
print(shares_current.round(4))
plot_shares(shares_current, "Market Share - Current Products")

## Market Share - Proposed Products (add the CR smart watch)
cr = pd.DataFrame({"product": ["CR"], "brand": ["CR"], "ship": ["$10"], "restock": ["10%"],
                   "retdays": ["14 days"], "price": ["$250"]})
proposed = pd.concat([cr, current], ignore_index=True)
util_proposed = calc_util(proposed)
print(util_proposed.round(3))
shares_proposed = calc_shares(util_proposed)
print(shares_proposed.round(4))
plot_shares(shares_proposed, "Market Share - Proposed Products")