This page describes a strategy for calculating pi one digit at a time using a streaming algorithm that does not require floating point arithmetic. Although there are more exotic strategies for calculating pi, this one is surprisingly simple to implement, and it streams the digits of pi infinitely. Before we dive in, give it a try in your browser!
It has long been known that pi can be represented as the product of an infinite series. The Wallis Product published in 1656 represents Pi as:
Which can be expanded to:
And alternatively expressed as:
The paper Spigot Algorithm for the Digits of Pi (1995) by
Stanley Rabinowitz and Stan Wagon describes how this infinite product series can be calculated in base 10 to stream the output digit by digit. The algorithm is is described in the paper on page 5. All mathematical operations are achieved using integer arithmetic, completely avoiding issues associated with floating point error accumulation. It is worth noting that the integers are large, so overflows with typical integer data types will occur quickly. It is necessary to use arbitrary length integer data types (BigInt) or to store large integers in an array pre-sized according to the number of digits to be calculated (a little over 3 times the number of digits) which is O(N) space complexity. Because the algorithm requires iterating the array fully for every digit to be calculated, this strategy runs in O(N2) time.
Appendix 2 has a useful bit of code attributed to “Macalester student Simeon Simeonov” implementing the algorithm in Pascal. It notes that this code makes use of the fact that the queue of predigits always has a pile of 9s to the right of its leftmost member, and so only this leftmost predigit and the number of 9s need be remembered.
Program PiSpigot;
const n = 1000;
len = 10 * n div 3;
var i, j, k, q, x, nines, predigit: integer;
a: array[1..len] of longint;
begin
for j := 1 to len do a[j] := 2; {Start with 2s}
nines := 0;
predigit := 0; {First predigit is a 0}
for j := 1 to n do begin
q := 0;
for i := len downto 1do {Work backwards}
begin
x := 10 * a[i] + q * i;
a[i] := x mod (2 * i - 1);
q := x div (2 * i - 1);
end;
a[1] := q mod 10;
q := q div 10;
if q = 9 then
nines := nines + 1elseif q = 10 then
begin
write(predigit + 1);
for k := 1 to nines do write(0); {zeros}
predigit := 0;
nines := 0;
end
else begin
write(predigit);
predigit := q;
if nines <> 0 then
begin
for k := 1 to nines do write(9);
nines := 0;
end;
end;
end;
writeln(predigit);
end.
int maxDigits = 1000;
int len = (10 * maxDigits) / 3;
int[] a = newint[len + 1];
int nines = 0;
int predigit = 0;
for (int j = 1; j <= len; j++)
a[j] = 2;
for (int j = 1; j <= maxDigits; j++)
{
int q = 0;
for (int i = len; i >= 1; i--)
{
int x = 10 * a[i] + q * i;
a[i] = x % (2 * i - 1);
q = x / (2 * i - 1);
}
a[1] = q % 10;
q /= 10;
if (q == 9)
{
nines++;
}
elseif (q == 10)
{
Console.Write(predigit + 1);
for (int k = 1; k <= nines; k++)
Console.Write(0);
predigit = 0;
nines = 0;
}
else {
if (j > 1)
Console.Write(predigit);
if (j == 2)
Console.Write(".");
predigit = q;
if (nines != 0)
{
for (int k = 1; k <= nines; k++)
Console.Write(9);
nines = 0;
}
}
}
It generates the first thousand digits of pi in about 60 milliseconds. The first 10,000 took 3.3 sec, and additional performance benchmarks can be found at the bottom of this page. The output is accurate according to some online tables I found of pi to a large number of decimal places.
Performance of the JavaScript implementation is surprisingly good. You can generate virtually infinite digits of pi using your desktop or mobile device. You never know when you’re going to need a few thousand digits of pi, and it is always good to be prepared.
This a bak of the envelope test just using my moderate Desktop PC and the C# implementation above, but performance is pretty good for several thousand digits and balloons after that as expected for an algorithm with O(N2) time complexity.
I spent the last few years experimenting with designs for internet-enabled temperature sensors so I could monitor the relationship between my indoor, outdoor, garage, and attic temperatures over long time spans. This page documents my favorite implementations using NodeMCU, a USB powered and programmable Arduino-compatible IoT development board with an ESP8266 WiFi module that can be purchased on Amazon in small quantities for about $5 each. The design presented on this page uses a DS18B20 1-Wire thermometer packaged in a weatherproof enclosure sold in small quantities for about $2 each.
Connecting the 1-Wire sensor only requires adding an external resistor to pull the data pin high. I found it convenient to solder the resistor and sensor wires directly to pins of the NodeMCU development board. It is worth noting that the DS18B20 is available in several different form factors, but the communication protocol and wiring requirements are the same.
The NodeMCU can be programmed via its micro USB connection using the Arduino IDE. Full source code for this project is available on GitHub, but the following sketch is a trimmed-down version which is likely sufficient for others to replicate this project.
#include<ESP8266WiFi.h>#include<ESP8266HTTPClient.h>#include<OneWire.h>#include<DallasTemperature.h>OneWire oneWire(D2); // Data wire connected to D2
DallasTemperature sensors(&oneWire);
voidsetup() {
sensors.begin();
WiFi.begin("WIFI_SSID", "WIFI_PASSWORD"); // Customize for your network
while (WiFi.status() != WL_CONNECTED) {delay(1000);}
WiFi.setAutoReconnect(true);
WiFi.persistent(true);
}
voidloop() {
// Read temperature
sensors.requestTemperatures();
float temp = sensors.getTempFByIndex(0);
// Prepare the endpoint URL
String resource = String("https://swharden.com/weather/v1/write/");
// Create JSON body for the POST request
String data = String("{")
+ String("\"key\": \"MY_SECRET_API_KEY\",") // Customize for security
+ String("\"sensor\": \"3\",") // Sensor ID for multi-sensor environments
+ String("\"temperature\": \"") + String(temp) + String("\",")
+ String("}");
// Create a WiFi client that does not use SSL
WiFiClientSecure client;
client.setInsecure();
client.connect("swharden.com", 443);
// Make the POST request including our custom JSON body
HTTPClient http;
http.begin(client, resource);
http.addHeader("Content-Type", "application/json");
http.POST(data);
http.end();
// Wait one minute before retrying
int seconds =60;
while (seconds--) {delay(1000);}
}
Users will want to customize their web server architecture, but I achieved a live temperature dashboard using a PHP endpoint that interprets the body of the POST request and logs the time/sensor/temperature as a single line appended to a CSV file. The interactive graphs are achieved using Dygraph, a JavaScript library ideal for interactively displaying time-series data on webpages.
I am fortunate to already have easy access to an indoor/outdoor pass-through created to run HF and VHF coax cabling from indoor radios to outdoor antennas. This made it easy to house the weather-sensitive NodeMCU module inside while snaking the weather-insensitive temperature sensor outside. This is on the North side of my house too, so it’s always in shadow and it should be a pretty good location for direct air temperature measurement.
The sensor in my garage is perhaps the cleanest build. I kept the USB cable short and only had the tip of the temperature sensor stick out the top. The NodeMCU has blue LEDs on it already, so I drilled a hole in the enclosure (sealed off with clear plastic glued to the underside of the hole) allowing the blinking light to be observed from outside. Short flashes indicate it is waiting until the next measurement, and it turns solid while it’s measuring temperature and submitting the HTTP request.
Originally I created a “weatherproof” enclosure to house the NodeMCU microcontroller and temperature sensor, powered using low voltage delivered through landscaping wire typically used for lard lighting. This design used a BMP280 sensor (described below) to measure temperature and barometric pressure. I put “weatherproof” in quotes because after almost two years it failed. The experience was an interesting (and humbling) venture into designing IoT gear intended to stand up against the elements. For what it’s worth, it did make it through two Florida hurricanes.
Previous WiFi temperature sensor designs used the BMP280 I2C temperature and pressure sensor which I had on hand from previous precision pressure measurement project, but the software required to read temperature and pressure was surprisingly complex due to the calibration procedure described in the datasheet. In contrast, the design using the SD18B20 digital temperature sensor was extremely easy to implement. There are Arduino sketches on this project’s GitHub repository demonstrating how to use both sensors with NodeMCU.
I recently had the need to convert a segmented line (polyline) with arbitrary node points to one with evenly spaced points. After exploring a few solutions I settled on the following two strategies because of their simplicity and performance. Code on this page avoids trigonometric functions and square root operations to enhance performance.
Simple linear interpolation may be employed to place points at specific distances along the length of the line. This method minimizes computations but it may place points closer together in 2D space than the specified distance, a phenomenon which is exacerbated when the original polyline contains very sharp angles.
// Start with a segmented line with arbitrary pointsList<Vector2> polyline1 = [new(1, 1), new(3, 8), new(5, 2), new(9, 9)];
// Generate a new segmented line with evenly spaced pointsdouble spacing = 1;
List<Vector2> polyline2 = ResampleAlongLine(polyline1, spacing);
List<Vector2> ResampleAlongLine(List<Vector2> polyline1, float spacing)
{
List<Vector2> polyline2 = [];
// start at the same first point as the original line polyline2.Add(new(polyline1[0].X, polyline1[0].Y));
float distanceToNextPoint2 = spacing;
for (int i = 1; i < polyline1.Count; i++)
{
Vector2 p0 = polyline1[i - 1];
Vector2 p1 = polyline1[i];
float dx = p0.X - p1.X;
float dy = p0.Y - p1.Y;
float distanceToNextPoint1 = (float)Math.Sqrt(dx * dx + dy * dy);
while (distanceToNextPoint2 <= distanceToNextPoint1)
{
float t = distanceToNextPoint2 / distanceToNextPoint1;
float x = p0.X + t * (p1.X - p0.X);
float y = p0.Y + t * (p1.Y - p0.Y);
polyline2.Add(new Vector2(x, y));
// consider the point we just placed p0 and re-evaluate p0 = polyline2.Last();
distanceToNextPoint1 -= distanceToNextPoint2;
distanceToNextPoint2 = spacing;
}
// move backwards as needed before advancing to the next given pair distanceToNextPoint2 -= distanceToNextPoint1;
}
return polyline2;
}
This implementation walks along interpolated line segments in short steps and places new points as the target distance from the previous point measured in 2D space has been met. This strategy can be significantly slower than the linear interpolation only method described above, but this implementation ensures that sequential points are always the target distance apart in 2D space (equal Euclidean spacing on a Cartesian coordinate system) and may be preferred for some applications. The accuracy can be increased by reducing the epsilon at the expense of increased computing time. Performance may be improved by starting with a large epsilon and decreasing it dynamically as the calculated distance between the test point and previous point approaches the target separation.
// Start with a segmented line with arbitrary pointsList<Vector2> polyline1 = [new(1, 1), new(3, 8), new(5, 2), new(9, 9)];
// Generate a new segmented line with evenly spaced pointsdouble spacing = 1;
List<Vector2> polyline2 = ResampleIn2D(polyline1, spacing);
List<Vector2> ResampleIn2D(List<Vector2> polyline, float spacing, float epsilon = 0.0001f)
{
List<Vector2> polyline2 = [];
// start at the same first point as the original line polyline2.Add(new(polyline[0].X, polyline[0].Y));
// compare this because calculating square roots is very costlyfloat spacingSquared = spacing * spacing;
// consider each line segment between pairs of points in the given polylinefor (int i = 1; i < polyline.Count; i++)
{
Vector2 p1 = polyline[i - 1];
Vector2 p2 = polyline[i];
// use epsilon to determine how many points along the line to tryfloat sdx = p1.X - p2.X;
float sdy = p1.X - p2.X;
float segLength = (float)Math.Sqrt(sdx * sdx + sdy * sdy);
int stepCount = (int)(segLength / epsilon);
// check the distance at each step along the segmentfor (int j = 0; j <= stepCount; j++)
{
// decide the test point using linear interpolationfloat fraction = (float)j / stepCount;
float lx = p1.X + (p2.X - p1.X) * fraction;
float ly = p1.Y + (p2.Y - p1.Y) * fraction;
Vector2 pTest = new(lx, ly);
// add the test point if it's the right distance from the last onefloat dx = polyline2.Last().X - pTest.X;
float dy = polyline2.Last().Y - pTest.Y;
float pTestDistanceSquared = (dx * dx + dy * dy);
if (pTestDistanceSquared >= spacingSquared)
{
polyline2.Add(pTest);
}
}
}
return polyline2;
}
Below is the single file .NET 9 Console Application I used to generate the figures on this page.ScottPlot.net was used to generate the images. Note that Vector2 has a Lerp() method for linear interpolation and Distance() and DistanceSquared() methods for measuring distances between points, but I chose to implement these functions discretely to make the code on this page easier to translate to other languages.
View the full source code
usingSystem.Numerics;
List<Vector2> polyline1 = [new(1, 1), new(3, 8), new(5, 2), new(9, 9)];
for (int i = 1; i <= 10; i++)
{
List<Vector2> polyline2 = ResampleAlongLine(polyline1, i);
Plotting.Plot(polyline1, polyline2, i, "line", "Even Spacing Along Polyline", ScottPlot.Colors.Black.WithAlpha(.5));
List<Vector2> polyline3 = ResampleIn2D(polyline1, i);
Plotting.Plot(polyline1, polyline3, i, "2d", "Even Spacing in 2D Space", ScottPlot.Colors.Magenta.WithAlpha(.5));
ShowDistances(i, polyline3);
}
staticvoid ShowDistances(float expected, List<Vector2> polyline)
{
var distances = Enumerable.Range(0, polyline.Count - 1)
.Select(i => Vector2.Distance(polyline[i], polyline[i + 1]));
Console.WriteLine($"[{expected}] " + string.Join(", ", distances));
}
static List<Vector2> ResampleAlongLine(List<Vector2> polyline1, float spacing)
{
if (polyline1.Count < 2)
thrownew InvalidDataException("need at least 2 points");
if (spacing <= 0)
thrownew InvalidDataException("spacing must be positive");
List<Vector2> polyline2 = []; // will be filled with evenly-spaced points polyline2.Add(new(polyline1[0].X, polyline1[0].Y)); // start at the same first pointfloat distanceToNextPoint2 = spacing;
for (int i = 1; i < polyline1.Count; i++)
{
Vector2 p0 = polyline1[i - 1];
Vector2 p1 = polyline1[i];
float dx = p0.X - p1.X;
float dy = p0.Y - p1.Y;
float distanceToNextPoint1 = (float)Math.Sqrt(dx * dx + dy * dy);
while (distanceToNextPoint2 <= distanceToNextPoint1)
{
float t = distanceToNextPoint2 / distanceToNextPoint1;
float x = p0.X + t * (p1.X - p0.X);
float y = p0.Y + t * (p1.Y - p0.Y);
polyline2.Add(new Vector2(x, y));
// consider the point we just placed p0 and re-evaluate p0 = polyline2.Last();
distanceToNextPoint1 -= distanceToNextPoint2;
distanceToNextPoint2 = spacing;
}
// move backwards as needed before advancing to the next given pair distanceToNextPoint2 -= distanceToNextPoint1;
}
return polyline2;
}
static List<Vector2> ResampleIn2D(List<Vector2> polyline, float spacing, float epsilon = 0.0001f)
{
List<Vector2> polyline2 = [];
// start at the same first point as the original line polyline2.Add(new(polyline[0].X, polyline[0].Y));
// compare this because calculating square roots is very costlyfloat spacingSquared = spacing * spacing;
// consider each pair of points in the given polylinefor (int i = 1; i < polyline.Count; i++)
{
Vector2 p1 = polyline[i - 1];
Vector2 p2 = polyline[i];
// use epsilon to determine how many points along the line to tryint stepCount = (int)(Vector2.Distance(p1, p2) / epsilon);
for (int j = 0; j <= stepCount; j++)
{
// decide the test point using linear interpolationfloat fraction = (float)j / stepCount;
float lx = p1.X + (p2.X - p1.X) * fraction;
float ly = p1.Y + (p2.Y - p1.Y) * fraction;
Vector2 pTest = new(lx, ly);
// add the test point if it's the right distance from the last onefloat dx = polyline2.Last().X - pTest.X;
float dy = polyline2.Last().Y - pTest.Y;
float pTestDistanceSquared = (dx * dx + dy * dy);
if (pTestDistanceSquared >= spacingSquared)
{
polyline2.Add(pTest);
}
}
}
return polyline2;
}
publicstaticclassExtensions{
publicstatic ScottPlot.Coordinates[] ToCoordinates(this List<Vector2> points)
{
return [.. points.Select(x => new ScottPlot.Coordinates(x.X, x.Y))];
}
}
staticclassPlotting{
publicstaticvoid Plot(List<Vector2> points1, List<Vector2> points2, int spacing, string name, string legend, ScottPlot.Color color)
{
ScottPlot.Plot plot = new();
var sp1 = plot.Add.Scatter(points1.ToCoordinates());
sp1.Color = ScottPlot.Colors.C1;
sp1.LineWidth = 2;
sp1.MarkerSize = 7;
sp1.LegendText = "Original Polyline";
var markers = plot.Add.Markers(points2.ToCoordinates());
markers.Color = color;
markers.MarkerShape = ScottPlot.MarkerShape.OpenCircle;
markers.MarkerSize = 13;
markers.MarkerLineWidth = 2;
markers.LegendText = legend;
plot.Axes.SetLimits(0, 10, 0, 10);
plot.SavePng($"spacing-{name}-{spacing}.png", 700, 400);
}
}
Spline Interpolation with ScottPlot - Demonstrates additional types of interpolation: Bezier, Catmull-Rom, Chaikin, Cubic, etc. The project is open source under a MIT license.
The simple moving average (SMA) for an array of length N using a window of size W can be efficiently calculated using a running sum maintained in memory.
As the window moves along the data, new values are shifted into it by adding their value to the running sum, and old values are shifted out by subtracting their value from the running sum.
This strategy reduces the number of operations from approximately O(N*W) to O(N) which offers considerable performance enhancement for applications requiring large smoothing windows.
double[] SmoothForward(double[] data, int windowSize)
{
double[] smooth = newdouble[data.Length];
double runningSum = 0;
int pointsInSum = 0;
for (int i = 0; i < smooth.Length; i++)
{
runningSum += data[i];
if (pointsInSum < windowSize)
{
pointsInSum++;
smooth[i] += runningSum / pointsInSum;
continue;
}
runningSum -= data[i - windowSize];
smooth[i] += runningSum / windowSize;
}
return smooth;
}
Applying two successive moving window averages in opposite directions is a quick way to achieve bidirectional smoothing and prevents the “trailing” effect caused by a lagging running average which can be seen in the example above.
Bidirectional smoothing doubles the number of operations from O(N) to O(N*2) but this strategy used here is still far more performant than O(N*W) realized by convolution with a triangular kernel which is expected to produce a similar result.
double[] SmoothBidirectional(double[] data, int windowSize)
{
double[] smooth = newdouble[data.Length];
// smooth from left to rightdouble runningSum = 0;
int pointsInSum = 0;
for (int i = 0; i < smooth.Length; i++)
{
runningSum += data[i];
if (pointsInSum < windowSize)
{
pointsInSum++;
smooth[i] += runningSum / pointsInSum;
continue;
}
runningSum -= data[i - windowSize];
smooth[i] += runningSum / windowSize;
}
// smooth from right to left runningSum = 0;
pointsInSum = 0;
for (int i = smooth.Length - 1; i >= 0; i--)
{
runningSum += data[i];
if (pointsInSum < windowSize)
{
pointsInSum++;
smooth[i] += runningSum / pointsInSum;
continue;
}
runningSum -= data[i + windowSize];
smooth[i] += runningSum / windowSize;
}
// average the two directionsfor (int i = 0; i < smooth.Length; i++)
smooth[i] /= 2;
return smooth;
}
Successively smoothing data using repeatedly halved window sizes is a strategy I have recently found to be extremely useful for aggressively smoothing large amounts of signal data over long timescales. This HBSMA strategy is especially good at smoothing discontinuous data with “jumps” in the signal. With a performance approximating O(N*log(W)*2) it offers vastly superior performance over convolution with a Gaussian kernel, while producing roughly similar output.
I think of this strategy has having three primary functions that each act according to the size of the window:
Large windows: smooth jumps in discontinuous data
Small windows: smooth high frequency noise
Intermediate windows: smooth transitions between jumps
A looping strategy is demonstrated here, but a recursive algorithm could also be realized to achieve the same effect.
Let’s compare results of the halving bidirectional simple moving average (HBSMA) with a traditional convolution strategy using a Gaussian kernel.
In this strategy our rectangular window is replaced with a bell-shaped window (a kernel) which is used to produce a weighted moving average for each point in the original dataset.
A Gaussian window provides bidirectional smoothing,
but since it never reaches zero on the edges I prefer using
a Hanning window for applications like this.
The HBSMA strategy produces a result strikingly similar with the Gaussian convolution, but in less than one tenth of the time.
Performance benchmarking on my system revealed a 10x improvement in speed for the data above,
but demonstrated over 130x performance improvement for larger signals (10,000 points) smoothed with a bigger window (100 points).
It is worth noting that moving window averaging strategies must account for the situation at the ends of the original data where
there are not enough data points to fill the window. In these cases a “padding” strategy may be used to fill the window to ensure
the same number of data points are averaged to create each value on the smoothed curve.
The examples in this article do not use fancy padding (their behavior is best described by the “None” option below),
but it is worth being aware of common padding strategies:
None - Calculated the average from fewer values at the edge of the data
Constant - Pad with a repeated user-defined value
Mirror - Pad with data values mirrored at the ege of the data
Nearest - Pad with repeated copies of the value at the edge of the data
Wrap - Wrap data values around so one end is averaged with the opposite end
The Halving Bidirectional Simple Moving Average (HBSMA) strategy enables data smoothing with results similar to Gaussian convolution
but in a small fraction of the time. The performance of HBSMA is approximately O(N*log(M)*2) whereas the convolution strategy is approximately
O(N*M), and the larger the smoothing window size is the more gains may be had by this data smoothing strategy.
This article describes how to fit arbitrary curves with any number of parameters to X/Y data points using C#. The strategy described here uses particle swarm optimization to iteratively work toward an optimal solution. Unlike other curve fitting strategies, the first derivative of the error function does not need to be calculated mathematically, making it particularly useful for problems where the solution landscape is noisy, nonlinear, or completely unknown. Users only need to provide a function and the limits for each parameter, and the swarm will converge on the optimal solution.
SwarmFit is a free and open-source NuGet package that makes it easy to fit arbitrary curves to data. I created SwarmFit after researching the topic of particle swarm optimization in preparation for this article. Although this page details how particle swarm optimization works under the hood, SwarmFit be used without any knowledge of the topic.
Add the SwarmFit NuGet package: dotnet add SwarmFit
Add the following to Program.cs then dotnet run
// data points to fitdouble[] xs = [1, 2, 3, 4, 5];
double[] ys = [304, 229, 174, 134, 111];
// a custom function to fit to the data using any number of parametersstaticdouble MyFitFunc(double x, double[] parameters)
{
double a = parameters[0];
double b = parameters[1];
double c = parameters[2];
return a + b * Math.Exp(x * c);
}
// the minimum and maximum value for each parameterdouble[] minVars = [-100, -5000, -10];
double[] maxVars = [100, 5000, 10];
// perform the fit with general purpose settingsdouble[] solution = QuickFit.Solve(xs, ys, MyFitFunc, minVars, maxVars);
// display the solutiondouble a = solution[0];
double b = solution[1];
double c = solution[2];
Console.WriteLine($"Y = {a} + {b} * e^(x * {c})");
Advanced users may wish to use F12 to peek inside the QuickFit.Solve() method to see how it creates a SwarmFitter, initializes it, and runs the solver. The SwarmFitter has many properties which may be adjusted to change the behavior of the fitter for enhanced logging, improved performance, or potentially better fits. See the SwarmFit GitHub repository for additional information.
Particle swarm optimization is a type of swarm intelligence algorithm that represents candidate solutions as particles in a swarm that all work together to minimize the error function. In this strategy, a set of parameter values is represented by a particle, and a small swarm of particles iteratively works toward a solution by each of their parameters to minimize the total difference between the data points and the curve generated by those parameters.
Particles move toward the optimal solution by combining their own values and recent experience with information gathered from neighboring particles. The combination of self-guided and socially guided behavior allows the swarm to explore a broad range of solutions while progressively converging onto the best one. Randomness is added to the system by randomly killing particles and spawning new ones, increasing the chances that particles will find the global minimum and not get trapped in a local one.
Here’s a summary of how particle swarm optimization may be implemented in C#. Full source code may be reviewed from the SwarmFit GitHub repository and additional information is available in the resources section at the bottom of the page.
Data points: Start with a collection of X/Y points represented by xs[] and ys[] (two arrays of equal length).
Curve formula: The user must create a function which returns Y given an X position and a collection of parameters. In this example we’ll assume parameters are in an array named vars[]. In C# a static function can be represented by Func<> making it easy to pass into a fitter class as a parameter.
Error function: The fitter must be able to evaluate how well a set of vars[] performs, so predict the Y for every X in the dataset. This can be achieved by returning the sum of the differences between each X’s actual Y value and the predicted Y value obtained from the formula given that X and the current parameters.
Parameter limits: The minimum and maximum possible limits for each parameter must be prepared by the user. These limits are used when initializing the swarm (placing particles randomly across the whole field) and also when adding randomness to the system (randomly displacing particles by a distance relative to the field size).
Initialize the swarm: Create a population of particles randomly placed on the field and given a random movement direction.
Evaluate fitness: For each particle, calculate how well the curve defined by its parameters fits the data.
Track best positions: Each particle keeps track of its own best known position (the one with the smallest error), and the optimizer keeps track of the best known position across all particles.
Update movement vector: Adjust the direction and speed of each particle based on its own best solution and the best solution found by the swarm.
Move particles: Move each particle forward according to its new direction and speed. This new set of particles represents new and potentially improved solutions.
Add randomness: Randomly kill particles and add new ones. This step is optional, but I find it particularly interesting. It is recommended to place new particles a random distance from the population mean which is proportional to the variance of the population.
Iterate: Repeat the process of adjusting movement, updating positions, and evaluating fitness over several iterations, allowing the swarm to progressively improve its collective solution.
Return the best solution: Once a certain number of iterations is reached, the best set of curve parameters is returned as the final solution.
This article describes how to perform a simple exponential fit using pure C#. Exponential curves which pass through 0 can be described using two variables, so data can be translated into linear space, fitted using a linear least squares strategy, then translated back into exponential space.
// These are the data points we wish to fitdouble[] xs = [1, 2, 3, 4, 5, 6, 7];
double[] ys = [258, 183, 127, 89, 65, 48, 35];
// Log-transform values, then perform a linear fitdouble[] logYs = ys.Select(x => Math.Log(x)).ToArray();
(double slope, double intercept) = LeastSquaresFit(xs, logYs);
// Convert back to exponential form and display the resultdouble a = Math.Exp(intercept);
Console.WriteLine($"y = {a}*e^({slope}*x)");
The code above is ideal for calculating two terms to generate a fitted curve according to the formula:
Y = A * eB * x
…but what about fitting to a curve with third term describing a Y offset?
Y = C + A * eB * x
…or a forth term describing an X offset too?
Y = C + A * eB * x + D
The linear least squares fitting strategy described above is not sufficient for fitting data to these complex curves. The linear least squares strategy can be applied in higher dimensional space using gradient descent to seek optimized curves for more advanced equations. This strategy is summarized by The Math Coffeeshop in a fantastic YouTube video about using Linear Least Squares to Solve Nonlinear Problems
. However, formulas representing the partial derivative of each unknown in the error function must be calculated which may be difficult. An alternative strategy is to use Particle swarm optimization to iteratively work toward an optimal solution for a fitted curve even if the derivative of the error function cannot be calculated. This strategy is described in my follow-up article, Curve Fitting in C# using Particle Swarm Optimization
Adding a permissive license to your open source projects makes it easy for others to benefit from reusing your code. I typically recommend the MIT license for most projects because it’s simple, well understood, and allows the code to be used for essentially any purpose. Using code from websites and GitHub repositories without a license is difficult, because as I understand it the default terms are “all rights reserved” by the copyright holder unless explicitly stated otherwise. Since not having a license puts restrictions on how others may use your code, I created the GitHub License Checker tool to help identify GitHub repositories with missing license files.
It’s a Vanilla JavaScript app that works by hitting the paginated “repos” GitHub API for a given username. Because a page has up to 100 repositories and most users are unlikely to have more than a few hundred repositories, the API rate limiting is rarely engaged so a GitHub API key is not required.
The GitHub repositories API returns a lot of useful information. Because it’s so easy to use client-side without an API key, I’m curious what other useful tools could be created using it. Perhaps ScottPlot could be used to generate some interesting charts. Ideas that jump out are plotting every repo’s star count vs size, or date vs number of commits. I’ll leave that for another day, but C# developers should note that the Octokit NuGet Package is an official .NET package for interacting with the GitHub API.
