Compare commits

...

2 Commits

2 changed files with 64 additions and 22 deletions

View File

@@ -9,6 +9,7 @@ import logging
import argparse import argparse
from pathlib import PurePath##https://docs.python.org/3/library/pathlib.html#module-pathlib from pathlib import PurePath##https://docs.python.org/3/library/pathlib.html#module-pathlib
import numpy as np import numpy as np
import matplotlib
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from scipy.stats import norm from scipy.stats import norm
import pandas as pd import pandas as pd
@@ -36,16 +37,24 @@ parser_sim.add_argument('mean', type=float,
parser_sim.add_argument('stddev', type=float, parser_sim.add_argument('stddev', type=float,
help="sample standard deviation") help="sample standard deviation")
## General Arguments ## General Arguments
parser.add_argument('minvalue', type=float, parser.add_argument('--lowerbound', type=float,
help='Tolerance low limit') help='Tolerance low limit')
parser.add_argument('maxvalue', type=float, parser.add_argument('--upperbound', type=float,
help='Tolerance high limit') help='Tolerance high limit')
parser.add_argument('--normalizey', action="store_true", parser.add_argument('--normalizey', action="store_true",
help='Set y-axis to normalized probability density') help='Set y-axis to normalized probability density')
parser.add_argument('--log', default="INFO", parser.add_argument('--log', default="INFO",
help='Console log level: Number or DEBUG, INFO, WARNING, ERROR') 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", parser.add_argument('--graphinfo', action="store_true",
help='Put information on the PDF graph') 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() args = parser.parse_args()
@@ -72,10 +81,7 @@ fh.setFormatter(spamformatter)
logger.addHandler(ch) logger.addHandler(ch)
logger.addHandler(fh) logger.addHandler(fh)
logger.info("Creating infocalc log file %s", logpath) logger.debug("Creating infocalc log file %s", logpath)
lowerbound = args.minvalue
upperbound = args.maxvalue
# seed values for variable scoping # seed values for variable scoping
mean = 0 mean = 0
@@ -96,14 +102,26 @@ elif args.mode == "SIM":
stddev = args.stddev stddev = args.stddev
samplesize = args.samplesize samplesize = args.samplesize
# time to deal with the bounds
# Delta Degrees of Freedom: ddof=0 for population, ddof=1 for sample std dev # 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) #print("probability: %f", prob)
info = -np.emath.log2(prob) info = -np.emath.log2(prob)
#print("information content: %f bits", info) #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 ## place text on plot: https://matplotlib.org/3.3.4/gallery/recipes/placing_text_boxes.html
fig, ax = plt.subplots() fig, ax = plt.subplots()
if args.graphinfo:#put info on corner of graph if args.graphinfo:#put info on corner of graph
textstr = '\n'.join(( textstr = '\n'.join((
r'$n=%d$' % (samplesize), 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) props = dict(boxstyle='round', facecolor='wheat', alpha=0.5)
# place a text box in upper left in axes coords # 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) 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) y = norm.pdf(x, loc=mean, scale=stddev)
if args.normalizey: if args.normalizey:
y = y * stddev#rescale back to unity area y = y * stddev#rescale back to unity area
plt.axvline(x=mean, color="green", linestyle="dashed", label="mean") plt.axvline(x=mean, color="green", linestyle="dashed", label="mean")
plt.axvline(lowerbound, color="red") if args.lowerbound:
plt.axvline(upperbound, color="red") plt.axvline(args.lowerbound, color="red")
if args.upperbound:
plt.axvline(args.upperbound, color="red")
plt.plot(x, y, 'b-', label='Normal distribution') plt.plot(x, y, 'b-', label='Normal distribution')
#yt = scipy.stats.t.pdf(x, len(data)-1, mean, stddev) #yt = scipy.stats.t.pdf(x, len(data)-1, mean, stddev)
#plt.plot(x, yt, 'g-', label='T Distribution') #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] 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() plt.show()

View File

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