1 module matplotd.hist;
2 
3 import matplotd.base;
4 
5 struct HistogramConfig
6 {
7     /// title of the plot
8     string title;
9 
10     /// number of bins
11     int numBins = 100;
12 
13     /// step size of the points in the plot
14     double stepSize = 0.005;
15 
16     /// whether the histogram should be plotted with cumulative probabilities
17     bool cumulative = false;
18 
19     /// histogram type to plot ('bar', 'step', 'stepfilled')
20     string histType = "bar";
21 
22     /// color of the histogram lines
23     string color = "blue";
24 
25     /// whether the reference pdf should be plotted
26     bool plotReference = true;
27 }
28 
29 /**
30 Plot distribution histogram
31 
32 Params:
33     pdf = probability density function
34     values = sampled values from a distribution
35     fileName = path where the file should be saved
36 */
37 
38 void histogram(S, Pdf)(Pdf pdf, S[] values, string fileName, HistogramConfig config = HistogramConfig())
39 in
40 {
41     with(config)
42     assert(histType == "bar" || histType == "step" || histType == "stepfilled",
43                 "Invalid histogram type");
44 }
45 body
46 {
47     import pyd.embedded : InterpContext;
48     import pyd.extra : d_to_python_numpy_ndarray;
49     import std.algorithm.searching : maxPos, minPos;
50     import std.algorithm.iteration : sum;
51     import std.array : array;
52     import std.math : isFinite;
53     import std.range : front, iota;
54 
55     import std.stdio;
56     writeln("starting python hist plot");
57 
58     auto pythonContext = new InterpContext();
59     pythonContext.num_bins = config.numBins;
60     pythonContext.fileName = fileName;
61     pythonContext.title = config.title;
62     pythonContext.histType = config.histType;
63     pythonContext.py_stmts(`
64         import matplotlib.pyplot as plt
65         import numpy as np
66     `);
67 
68     // double is needed for NumPy
69     double[] npValues;
70     // apply filtering due to weird output errors
71     foreach (v; values)
72     {
73         if (v.isFinite)
74             npValues ~= v;
75     }
76     pythonContext.sample = npValues.d_to_python_numpy_ndarray;
77     pythonContext.cumulative = config.cumulative;
78     pythonContext.nMax = 0;
79     pythonContext.color = config.color;
80     pythonContext.py_stmts(`
81         n, bins, patches = plt.hist(sample, num_bins, normed=1,
82             cumulative=cumulative, histtype=histType, color=color)
83         nMax = np.max(n)
84     `);
85 
86     static if (!is(Pdf == typeof(null)))
87     {
88         if (config.plotReference)
89         {
90 
91             // plot actual density function
92             double[] xs = iota(values.minPos.front, values.maxPos.front, config.stepSize).array;
93             double[] ys = new double[xs.length];
94             foreach (i, x; xs)
95                 ys[i] = pdf(x);
96 
97             if (config.cumulative)
98             {
99                 // normalize
100                 auto total = ys.sum();
101                 foreach (ref y; ys)
102                     y /= total;
103                 foreach (i, ref y; ys[1..$])
104                     y += ys[i];
105             }
106             else
107             {
108                 // we try to scale the pdf according to highest bar of the histogram
109                 // this is not a 100% perfect solution
110                 auto nMax = pythonContext.nMax.to_d!double;
111                 auto factor = nMax / ys.maxPos.front;
112                 foreach (ref y; ys)
113                     y *= factor;
114             }
115 
116             pythonContext.xs = xs.d_to_python_numpy_ndarray;
117             pythonContext.ys = ys.d_to_python_numpy_ndarray;
118             pythonContext.py_stmts(`plt.plot(xs, ys, color='black')`);
119         }
120     }
121 
122     // save file
123     pythonContext.py_stmts(`
124         plt.title(title)
125         plt.savefig(fileName, bbox_inches='tight')
126         plt.close()
127     `);
128 }
129 
130 /// ditto
131 void histogram(S)(S[] values, string fileName, HistogramConfig config = HistogramConfig())
132 {
133     histogram(null, values, fileName, config);
134 }