.gitignore

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

Makefile

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

README.md

@@ -2,13 +2,5 @@
 
 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.
-
-`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.
+- `linreg.py` fit a polynomial to the data.

derive.py

@@ -22,7 +22,6 @@ 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 
@@ -47,7 +46,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:w:")
+    opt,files = getopt.getopt(sys.argv[1:], "PKx:y:s:lLvd:p:hf:N:F:S:W:")
 except:
     usage()
     raise
@@ -64,7 +63,6 @@ fmt="%.10g"
 N=2
 fit_range=5.
 minstep=0.
-xstep=None
 
 i_base = 1
 
@@ -98,9 +96,6 @@ for o,v in opt:
         if o=="-W":
             minstep = float(eval(v))
 
-        if o=="-w":
-            xstep = float(eval(v))
-
         if o=="-l":
             logscale = True
         
@@ -180,16 +175,14 @@ fit_range *= sigma
 sigma = 1/(math.sqrt(2)*sigma)
 
 p = numpy.arange(2*P+1).reshape((-1,1,1))
-Pr = numpy.arange(P+1)
-Ai = Pr+Pr.reshape((-1,1))
+Ai = numpy.arange(P+1)
+Ai = Ai+Ai.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:
@@ -226,11 +219,7 @@ 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)
@@ -245,22 +234,6 @@ 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,52 +109,7 @@ static v_index_t max_size = 0x1000000;
 static bool binary_out = false;
 static bool xyz_out = false;
 
-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 const char **files = NULL;
 
 static void parse_size(const char *arg, struct axis_para *a)
 {
@@ -282,7 +237,12 @@ 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 = optarg ? parse_column(optarg) : 1; break;
+		case 'K': column0 = 0;
+			if (optarg)
+				column0 = parse_column(optarg);
+			else
+				column0 = 1;
+			break;
 		case 'w': weights = parse_column(optarg); break;
 		case 'B': binary_out = true; break;
 		case '3': xyz_out = true; break;
@@ -304,12 +264,38 @@ 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);