#!/usr/bin/env python
## Axiomatic Design Information Calculator (and plotter)
## Author: Joseph Timothy foley <foley AT RU.IS>
## Start Date:  2026-02-27
## Input:  data in csv file
## Output:  information calculation and PDF for report/presentation
import os
import logging
import argparse
from pathlib import PurePath##https://docs.python.org/3/library/pathlib.html#module-pathlib
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
import pandas as pd

#Main program loop
print("""Axiomatic Design Information Calculator by Joseph. T. Foley<foley AT ru DOT is>
From https://gitea.cs.ru.is/AxiomaticDesign/adcalc/""")
parser = argparse.ArgumentParser(
    description="Axiomatic Design Information Calculator.")
parser.add_argument('csvfile',
                    help="CSV file with data and headers")
parser.add_argument('column',
                    help='Which column header to take data from')
parser.add_argument('minvalue', type=float,
                    help='Tolerance low limit')
parser.add_argument('maxvalue', type=float,
                    help='Tolerance high limit')
parser.add_argument('--normalizey', action="store_true",
                    help='Set y-axis to normalized probability density')
parser.add_argument('--log', default="INFO",
                    help='Console log level:  Number or DEBUG, INFO, WARNING, ERROR')
parser.add_argument('--graphinfo',
                    help='Put information on the PDF graph')
args = parser.parse_args()

## Set up logging
numeric_level = getattr(logging, args.log.upper(), None)
if not isinstance(numeric_level, int):
    raise ValueError(f'Invalid log level: {args.log}')
#print(f"Log level:  {numeric_level}")
logger = logging.getLogger("app")
logger.setLevel(numeric_level)
# log everything to file
logpath = os.path.splitext(args.csvfile)[0]+".log"
fh = logging.FileHandler(logpath)
fh.setLevel(logging.DEBUG)
# log to console
ch = logging.StreamHandler()
ch.setLevel(numeric_level)
# create formatter and add to handlers
consoleformatter = logging.Formatter('%(message)s')
ch.setFormatter(consoleformatter)
spamformatter = logging.Formatter('%(asctime)s %(name)s[%(levelname)s] %(message)s')
fh.setFormatter(spamformatter)
# add the handlers to logger
logger.addHandler(ch)
logger.addHandler(fh)

logger.info("Creating infocalc log file %s", logpath)

# filename pre-processing for output
inpath = PurePath(args.csvfile)
print(f"Input: {inpath}")

# grab the data and process
data = np.array(pd.read_csv(inpath)[args.column])
lowerbound = args.minvalue
upperbound = args.maxvalue
logger.debug(f"data:{data}, lower:{lowerbound}, upper:{upperbound}")

mean = data.mean()
stddev = data.std(ddof=1)
# Delta Degrees of Freedom: ddof=0 for population, ddof=1 for sample std dev
prob = norm.cdf(upperbound, mean, stddev) - norm.cdf(lowerbound, mean, stddev)
#print("probability: %f", prob)
info = -np.emath.log2(prob)
#print("information content: %f bits", info)
## place text on plot:  https://matplotlib.org/3.3.4/gallery/recipes/placing_text_boxes.html
fig, ax = plt.subplots()
textstr = '\n'.join((
    r'$n=%d$' % (len(data)),
    r'$\mu=%.2f$' % (mean, ),
    r'$\sigma=%.2f$' % (stddev, ),
    r'$P=%.2f$' % (prob, ),
    r'$I=%.2f$ bits' % (info, )))
# these are matplotlib.patch.Patch properties
props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)

# place a text box in upper left in axes coords
ax.text(0.05, 0.95, textstr, transform=ax.transAxes, fontsize=14,
        verticalalignment='top', bbox=props)

x = np.linspace(mean-3*stddev, mean+3*stddev, 500)
y = norm.pdf(x, loc=mean, scale=stddev)
if args.normalizey:
    y = y * stddev#rescale back to unity area
plt.axvline(x=mean, color="green", linestyle="dashed", label="mean")
plt.axvline(lowerbound, color="red")
plt.axvline(upperbound, color="red")

plt.plot(x, y, 'b-', label='Normal distribution')
#yt = scipy.stats.t.pdf(x, len(data)-1, mean, stddev)
#plt.plot(x, yt, 'g-', label='T Distribution')
coloredregion = (x >= lowerbound) & ( x <= upperbound ) #select x values
plt.fill_between(x, 0, y, where=coloredregion, color="grey", alpha=0.5, label="Design range")
plt.xlabel('X')
plt.ylabel('Probability density')
plt.legend()
plt.grid(True)

top = plt.ylim()[1]
plt.show()