derive.py: add option -w

.gitignore

@@ -1,2 +1,3 @@
 *~
 hist
+__pycache__

Makefile

@@ -1,6 +1,6 @@
 #! /usr/bin/make -f
 
-CFLAGS = -g
+CFLAGS = -g -O2
 LOADLIBES = -lm
 hist: hist.c
 

README.md

@@ -2,5 +2,13 @@
 
 Scripts for data processing and plotting.
 - `hist.py` create histogram to be plotted with gnuplot style `histeps`,
+- `hist.c`  same, faster
 - `derive.py` fit a polynomial to every data point with gaussian weight along *x*,
-- `linreg.py` fit a polynomial to the data.
+- `linreg.py` fit a polynomial to the data.
+
+`hist.c` is a programm written in C. Many regard programming in than
+language as a crime, because _“safety”_.  I consider too much safety
+to be a crime.  Gaining experience and training is a lot more
+important to keep our society viable. Let's made a lot of small
+mistakes and accidents to avoid the really dangerous ones, with
+competence instead of rules.

derive.py

@@ -22,6 +22,7 @@ Options:
   -N «nc»   number of coefficients to emit, default -N 2
   -S «s»    fit range in units of sigma
   -W «step» minimum distance of x-values
+  -w «step» increment of output x-values (implies -N2)
   -h        print this help
 
 If -y is given, but no -x, the row number is used for x.  When no -y is 
@@ -46,7 +47,7 @@ def usage():
     sys.stderr.write(__doc__ % (sys.argv[0],))
 
 try:
-    opt,files = getopt.getopt(sys.argv[1:], "PKx:y:s:lLvd:p:hf:N:F:S:W:")
+    opt,files = getopt.getopt(sys.argv[1:], "PKx:y:s:lLvd:p:hf:N:F:S:W:w:")
 except:
     usage()
     raise
@@ -63,6 +64,7 @@ fmt="%.10g"
 N=2
 fit_range=5.
 minstep=0.
+xstep=None
 
 i_base = 1
 
@@ -96,6 +98,9 @@ for o,v in opt:
         if o=="-W":
             minstep = float(eval(v))
 
+        if o=="-w":
+            xstep = float(eval(v))
+
         if o=="-l":
             logscale = True
         
@@ -175,14 +180,16 @@ fit_range *= sigma
 sigma = 1/(math.sqrt(2)*sigma)
 
 p = numpy.arange(2*P+1).reshape((-1,1,1))
-Ai = numpy.arange(P+1)
-Ai = Ai+Ai.reshape((-1,1))
+Pr = numpy.arange(P+1)
+Ai = Pr+Pr.reshape((-1,1))
 
 i1=0
 i2=0
 
 nn = XY.shape[0]
 last_x = x[0]
+if xstep:
+    next_x = xstep * math.floor(last_x/xstep-1)
 
 for i in range(nn):
     if x[i] < last_x+minstep:
@@ -219,7 +226,11 @@ for i in range(nn):
             else:
                 ii1 = i2
             i2 = (ii1+ii2)//2
-        
+
+    if i2-i1 <= P:
+        sys.stdout.write("\n")
+        continue
+
     x1  = numpy.add(x[i1:i2], -x[i], out=x1)
     y1  = numpy.add(y[:,i1:i2], -y[:,i].reshape((-1,1)), out=y1)
     xx  = numpy.power(x1, p, out=xx)
@@ -234,6 +245,22 @@ for i in range(nn):
     sxy = numpy.add.reduce(wxy, axis=-1, out=sxy)
     d   = numpy.linalg.solve(sxx.reshape((-1,))[Ai], sxy)
     d[0] += y[:,i]
+
+    if xstep:
+        if next_x < x[i] - xstep:
+            next_x =  xstep * math.floor(x[i]/xstep - 1)
+        while next_x < x[i]:
+            next_x += xstep
+            xint = next_x - x[i]
+            yint = sum(d * xint**Pr.reshape((-1,1)))
+            dint = sum(d[1:] * (xint**Pr[:-1] * Pr[1:]).reshape((-1,1)))
+            sys.stdout.write((fmt+" %s %s\n") %
+                             ( next_x,
+                               " ".join([fmt % yy for yy in yint]),
+                               " ".join([fmt % yy for yy in dint]),
+                              ))
+        continue
+
     if logscale:
         d[0] = numpy.exp(d[0])
         # d[1] *= d[0]

hist.c

@@ -56,12 +56,12 @@ struct axis_para {
 	v_index_t origin; // first index of the allocation
 	int column;       // column number 0…
 	int moments;      // multiply the weight by this power of .val
-	v_index_t size;      // allocated size
-	v_index_t fsize;     // fixed size
-	v_index_t usize;     // unused size
-	v_index_t nsize;     // unused size
-	v_index_t delta;     // origin shift
-	v_index_t stride;    // matrix stride
+	v_index_t size;   // allocated size
+	v_index_t fsize;  // fixed size
+	v_index_t usize;  // unused size
+	v_index_t nsize;  // unused size
+	v_index_t delta;  // origin shift
+	v_index_t stride; // matrix stride
 	bool integral;    // do integral histogram
 	bool logscale;    // take the log of .val
 	bool need_shift;  // need to move .origin
@@ -109,7 +109,52 @@ static v_index_t max_size = 0x1000000;
 static bool binary_out = false;
 static bool xyz_out = false;
 
-static const char **files = NULL;
+static void print_axis_config(const struct axis_para *a)
+{
+	fprintf(stderr, "Axis %s[%d]", a->name, a->column+column0);
+	if (a->integral)
+		fprintf(stderr, " int");
+	if (a->moments)
+		fprintf(stderr, " mom=%d", a->moments);
+	if (a->logscale)
+		if (a->logbase==1.0)
+			fprintf(stderr, " log=e");
+		else
+			fprintf(stderr, " log=%g", exp(a->logbase));
+	if (a->fsize)
+		fprintf(stderr, " bins %g…%g w %g n=%lld\n",
+			a->min_val, a->max_val, a->scale, a->max - a->min + 1);
+	else
+		fprintf(stderr, " bin widths %g\n", a->scale);
+}
+
+static void print_config(char **files)
+{
+	fprintf(stderr, "verbosity %d columns first[%d]",
+		verbosity, column0);
+	if (weights>=0)
+		fprintf(stderr, " weights[%d]", weights + column0);
+	if (separator)
+		fprintf(stderr, " sep[%c]\n", separator);
+	else
+		fprintf(stderr, " sep[space]\n");
+	print_axis_config(&ax);
+	if (twodim) {
+		print_axis_config(&ay);
+		if (binary_out)
+			fprintf(stderr, "output: binary matrix float\n");
+		else if (xyz_out)
+			fprintf(stderr, "output: x,y,z\n");
+		else
+			fprintf(stderr, "output: nonuniform matrix\n");
+	}
+	if (files && *files) {
+		fprintf(stderr, "Input files given:");
+		for (char **f=files; *f; f++)
+			fprintf(stderr, " %s", *f);
+		fprintf(stderr, "\n");
+	}
+}
 
 static void parse_size(const char *arg, struct axis_para *a)
 {
@@ -237,12 +282,7 @@ static char **parse_args(int argc, char **argv)
 		case 'C': twodim=true;
 		case 'c': a-> column = parse_column(optarg); break;
 		case 'P': column0 = 0; break;
-		case 'K': column0 = 0;
-			if (optarg)
-				column0 = parse_column(optarg);
-			else
-				column0 = 1;
-			break;
+		case 'K': column0 = optarg ? parse_column(optarg) : 1; break;
 		case 'w': weights = parse_column(optarg); break;
 		case 'B': binary_out = true; break;
 		case '3': xyz_out = true; break;
@@ -264,38 +304,12 @@ static void initialize();
 int main(int argc, char **argv)
 {
 	char **files = parse_args(argc, argv);
-	if (verbosity > 2) {
-		// TODO: fix
-		fprintf(stderr, "Options given:\n");
-		fprintf(stderr, " --verbosity=%d\n", verbosity);
-		if (ax.moments)
-			fprintf(stderr, " --moments\n");
-		else if (ax.integral)
-			fprintf(stderr, " --integral\n");
-		if (twodim)
-			fprintf(stderr, " --twodim\n");
-		fprintf(stderr, " --scale=%g\n", ax.scale);
-		if (twodim)
-			fprintf(stderr, " --Scale=%g\n", ay.scale);
-		if (ax.logscale)
-			fprintf(stderr, " --logscale=%.2g\n", exp(ax.logbase));
-		if (ay.logscale)
-			fprintf(stderr, " --Logscale=%.2g\n", exp(ay.logbase));
-		fprintf(stderr, " --column=%d\n", ax.column+column0);
-		if (twodim)
-			fprintf(stderr, " --Column=%d\n", ay.column+column0);
-		if (weights >= 0)
-			fprintf(stderr, " --weights=%d\n", weights+1);
-		if (*files) {
-			fprintf(stderr, "Input files given:");
-			for (char **f=files; *f; f++)
-				fprintf(stderr, " %s", *f);
-			fprintf(stderr, "\n");
-		}
-	}
 
 	initialize();
 
+	if (verbosity >= 2)
+		print_config(files);
+
 	FILE *inp = NULL;
 	size_t line_size = 256;
 	char *line = malloc(line_size);