Compare commits

..

2 Commits

2 changed files with 64 additions and 22 deletions

View File

@@ -9,6 +9,7 @@ import logging
import argparse
from pathlib import PurePath##https://docs.python.org/3/library/pathlib.html#module-pathlib
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from scipy.stats import norm
import pandas as pd
@@ -36,16 +37,24 @@ parser_sim.add_argument('mean', type=float,
parser_sim.add_argument('stddev', type=float,
help="sample standard deviation")
## General Arguments
parser.add_argument('minvalue', type=float,
parser.add_argument('--lowerbound', type=float,
help='Tolerance low limit')
parser.add_argument('maxvalue', type=float,
parser.add_argument('--upperbound', 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('--legend', action="store_true",
help='Put legend on the PDF graph')
parser.add_argument('--graphinfo', action="store_true",
help='Put information on the PDF graph')
parser.add_argument('--xlabel',
help='X-axis label, if needed')
parser.add_argument('--outfile',
help="output graph to PDF file")
parser.add_argument('--fontsize', default=14, type=int,
help="Adjust font size")
args = parser.parse_args()
@@ -72,10 +81,7 @@ fh.setFormatter(spamformatter)
logger.addHandler(ch)
logger.addHandler(fh)
logger.info("Creating infocalc log file %s", logpath)
lowerbound = args.minvalue
upperbound = args.maxvalue
logger.debug("Creating infocalc log file %s", logpath)
# seed values for variable scoping
mean = 0
@@ -96,14 +102,26 @@ elif args.mode == "SIM":
stddev = args.stddev
samplesize = args.samplesize
# time to deal with the bounds
# 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)
prob = 0
if args.upperbound and args.lowerbound:
prob = norm.cdf(args.upperbound, mean, stddev) - norm.cdf(args.lowerbound, mean, stddev)
elif args.upperbound:
prob = norm.cdf(args.upperbound, mean, stddev)
elif args.lowerbound:
prob = 1 - norm.cdf(args.lowerbound, mean, stddev)
else:
prob = 1# no bounds set!
##TODO!!!!
#print("probability: %f", prob)
info = -np.emath.log2(prob)
#print("information content: %f bits", info)
## set default fontsize
matplotlib.rcParams['font.size']=args.fontsize
## place text on plot: https://matplotlib.org/3.3.4/gallery/recipes/placing_text_boxes.html
fig, ax = plt.subplots()
if args.graphinfo:#put info on corner of graph
textstr = '\n'.join((
r'$n=%d$' % (samplesize),
@@ -115,26 +133,49 @@ if args.graphinfo:#put info on corner of graph
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,
ax.text(0.05, 0.95, textstr, transform=ax.transAxes, fontsize=args.fontsize,
verticalalignment='top', bbox=props)
xgraphlimits = {"min": mean-3*stddev, "max": mean+3*stddev}
if args.lowerbound and xgraphlimits["min"] > args.lowerbound:
xgraphlimits["min"] = args.lowerbound
if args.upperbound and xgraphlimits["max"] < args.upperbound:
xgraphlimits["max"] = args.upperbound
x = np.linspace(mean-3*stddev, mean+3*stddev, 500)
x = np.linspace(xgraphlimits["min"], xgraphlimits["max"], 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")
if args.lowerbound:
plt.axvline(args.lowerbound, color="red")
if args.upperbound:
plt.axvline(args.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)
# Filter for which region to fill
coloredregion = x#default fill all
if args.lowerbound and args.upperbound:
coloredregion = (x >= args.lowerbound) & ( x <= args.upperbound )
elif args.upperbound:
coloredregion = x <= args.upperbound
elif args.lowerbound:
coloredregion = x >= args.lowerbound
plt.fill_between(x, 0, y, where=coloredregion, color="grey", alpha=0.5, label="Design range",)
if args.xlabel:
plt.xlabel(args.xlabel)
plt.ylabel('Probability density')
if args.legend:
plt.legend()
#plt.grid(True)
top = plt.ylim()[1]
if args.outfile:
logger.info(f"Graph output to {args.outfile}")
plt.savefig(args.outfile,bbox_inches='tight')
else:
plt.show()

View File

@@ -1,4 +1,5 @@
#!/bin/bash
# Get infocalc.py from https://gitea.cs.ru.is/AxiomaticDesign/adcalc
echo "Loading data from file"
./infocalc.py DATA testdata.csv data1 0.9 1.1 --graphinfo
echo "Creating simulated curve from parameters"