Least-squares homework

homework/least-squares/Makefile

+CFLAGS = -Wall -O3 -std=gnu11 $(shell gsl-config --cflags)
+LDLIBS = -lm $(shell gsl-config --libs)
+
+LS_fit.png: data.txt fit.txt Makefile
+	echo '\
+	set terminal png;\
+	set output "$@";\
+	set key top right;\
+	set tics out;\
+	set xlabel "Time (days)";\
+	set ylabel "ln(activity)";\
+	set title "Radioactive decay";\
+	plot\
+	"data.txt" using 1:2:3 with errorbars title "Data points",\
+	"fit.txt" using 1:2 with line title "Least squares fit",\
+	"fit.txt" using 1:3 with line title "Least squares fit lower estimate",\
+	"fit.txt" using 1:4 with line title "Least squares fit high estimate",\
+	' | gnuplot
+
+data.txt: main
+	./$<
+
+fit.txt: main
+	./$<
+
+main: main.c functions.c
+
+.PHONEY: clean
+
+clean:
+	$(RM) main fit.txt data.txt LS_fit.png values.txt
\ No newline at end of file

homework/least-squares/functions.c

+#include "functions.h"
+
+void print_vector(char s[], gsl_vector* v) {
+	printf("%s\n",s);
+	for(int i=0;i<v->size;i++) {
+		printf("%15g",gsl_vector_get(v,i));
+		printf("\n");
+	}
+}
+
+void print_matrix(FILE* stream,gsl_matrix* M) {
+    for(int i = 0; i < M->size1; i++) {
+        for(int j = 0; j < M->size2; j++) {
+            fprintf(stream,"%15g\t", gsl_matrix_get(M, i, j));
+        }
+        fprintf(stream,"\n");
+    }
+}
+
+void GS_decomp(gsl_matrix* A, gsl_matrix* R) {
+    int n = A -> size1;
+    int m = A -> size2;
+    for(int i=0;i<m;i++) {
+        gsl_vector_view ai = gsl_matrix_column(A,i);
+        double Rii = gsl_blas_dnrm2(&(ai.vector));
+        gsl_matrix_set(R,i,i,Rii);
+        gsl_vector_scale(&(ai.vector),1/Rii);
+        for(int j=i+1;j<m;j++) {
+            gsl_vector_view aj = gsl_matrix_column(A,j);
+            double Rij;
+            gsl_blas_ddot(&(ai.vector),&(aj.vector),&Rij);
+            gsl_matrix_set(R,i,j,Rij);
+            gsl_vector* qiRij = gsl_vector_alloc(n);
+            gsl_vector_memcpy(qiRij,&(ai.vector));
+            gsl_vector_scale(qiRij,Rij);
+            gsl_vector_sub(&(aj.vector),qiRij);
+            gsl_vector_free(qiRij);
+        }
+    }
+}
+
+void back_sub(gsl_matrix* R, gsl_vector* x, gsl_vector* b) {
+    int m = x->size;
+    for(int i=m-1;i>=0;i--) {
+        double xi = gsl_vector_get(b,i);
+        for(int j=i+1;j<m;j++) {
+            xi += -gsl_matrix_get(R,i,j)*gsl_vector_get(x,j);
+        }
+        xi /= gsl_matrix_get(R,i,i);
+        gsl_vector_set(x,i,xi);
+    }
+}
+
+void GS_solve(gsl_matrix* Q, gsl_matrix* R, gsl_vector* b, gsl_vector* x) {
+    assert(Q->size1==Q->size2);
+    gsl_blas_dgemv(CblasTrans,1,Q,b,0,x);
+    back_sub(R,x,x);
+}
+
+void GS_inverse(gsl_matrix* Q, gsl_matrix* R, gsl_matrix* B) {
+    for(int j=0;j<Q->size2;j++) {
+        gsl_vector* ej = gsl_vector_alloc(Q->size1);
+        gsl_vector_set_basis(ej,j);
+        gsl_vector_view Bj = gsl_matrix_column(B,j);
+        GS_solve(Q,R,ej,&(Bj.vector));
+        gsl_vector_free(ej);
+    }
+}
+
+void LS_fit(data data, double (*f)(int,double), gsl_vector* c, gsl_vector* dc, gsl_matrix* cova) {
+    int nr_data = (data.xs)->size;
+    int nr_params = c->size;
+    gsl_matrix* A = gsl_matrix_alloc(nr_data,nr_params);
+    gsl_matrix* AtA = gsl_matrix_alloc(nr_params,nr_params);
+    gsl_matrix* R = gsl_matrix_alloc(nr_params,nr_params);
+    gsl_vector* b = gsl_vector_alloc(nr_data);
+    for(int i=0; i<nr_data; i++) {
+        for(int k=0; k<nr_params; k++) {
+            double Aik = f(k,gsl_vector_get(data.xs,i))/gsl_vector_get(data.dy,i);
+            gsl_matrix_set(A,i,k,Aik);
+        }
+        double bi = gsl_vector_get(data.ys,i)/gsl_vector_get(data.dy,i);
+        gsl_vector_set(b,i,bi);
+    }
+    gsl_blas_dgemm(CblasTrans,CblasNoTrans,1,A,A,0,AtA);
+
+    GS_decomp(A,R);
+    gsl_blas_dgemv(CblasTrans,1,A,b,0,c);
+    back_sub(R,c,c);
+    GS_decomp(AtA,R);
+    GS_inverse(AtA,R,cova);
+    for(int i=0;i<nr_params;i++) {gsl_vector_set(dc,i,sqrt(gsl_matrix_get(cova,i,i)));}
+
+    gsl_matrix_free(A);
+    gsl_matrix_free(AtA);
+    gsl_matrix_free(R);
+    gsl_vector_free(b);
+}
\ No newline at end of file

homework/least-squares/functions.h

+#ifndef FUNCTIONS_H
+#define FUNCTIONS_H
+
+#include<math.h>
+#include<time.h>
+#include<stdio.h>
+#include<assert.h>
+#include<gsl/gsl_vector.h>
+#include<gsl/gsl_matrix.h>
+#include<gsl/gsl_blas.h>
+#include<gsl/gsl_linalg.h>
+#define RND (double)rand()/10000
+
+typedef struct {gsl_vector* xs; gsl_vector* ys; gsl_vector* dy;} data;
+
+void print_vector(char s[], gsl_vector* v);
+void print_matrix(FILE* stream, gsl_matrix* M);
+void GS_decomp(gsl_matrix* A, gsl_matrix* R);
+void back_sub(gsl_matrix* R, gsl_vector* x, gsl_vector* b);
+void GS_solve(gsl_matrix* Q, gsl_matrix* R, gsl_vector* b, gsl_vector* x);
+void GS_inverse(gsl_matrix* Q, gsl_matrix* R, gsl_matrix* B);
+void LS_fit(data data, double (*f)(int,double), gsl_vector* c, gsl_vector* dc, gsl_matrix* cova);
+
+#endif
\ No newline at end of file

homework/least-squares/main.c

+#include "functions.h"
+
+double funcs(int i, double x) {
+    switch (i) {
+    case 0: return 1; break;
+    case 1: return x; break;
+    default: return NAN; break;
+    }
+}
+
+int main() {
+    gsl_vector* x = gsl_vector_alloc(9);
+    gsl_vector* y = gsl_vector_alloc(9);
+    gsl_vector* dy = gsl_vector_alloc(9);
+    gsl_vector* c = gsl_vector_alloc(2);
+    gsl_vector* dc = gsl_vector_alloc(2);
+    gsl_matrix* cova = gsl_matrix_alloc(2,2);
+    gsl_vector_set(x,0,1);gsl_vector_set(y,0,log(117.0));gsl_vector_set(dy,0,1.0/20.0);
+    gsl_vector_set(x,1,2);gsl_vector_set(y,1,log(100.0));gsl_vector_set(dy,1,1.0/20.0);
+    gsl_vector_set(x,2,3);gsl_vector_set(y,2,log(88.0));gsl_vector_set(dy,2,1.0/20.0);
+    gsl_vector_set(x,3,4);gsl_vector_set(y,3,log(72.0));gsl_vector_set(dy,3,1.0/20.0);
+    gsl_vector_set(x,4,6);gsl_vector_set(y,4,log(53.0));gsl_vector_set(dy,4,1.0/20.0);
+    gsl_vector_set(x,5,9);gsl_vector_set(y,5,log(29.5));gsl_vector_set(dy,5,1.0/20.0);
+    gsl_vector_set(x,6,10);gsl_vector_set(y,6,log(25.2));gsl_vector_set(dy,6,1.0/20.0);
+    gsl_vector_set(x,7,13);gsl_vector_set(y,7,log(15.2));gsl_vector_set(dy,7,1.0/20.0);
+    gsl_vector_set(x,8,15);gsl_vector_set(y,8,log(11.1));gsl_vector_set(dy,8,1.0/20.0);
+    data data = {.xs=x,.ys=y,.dy=dy};
+
+    LS_fit(data,&funcs,c,dc,cova);
+
+    FILE* data_points_stream = fopen("data.txt","w");
+    for(int i=0;i<9;i++) {fprintf(data_points_stream,"%g %g %g\n",gsl_vector_get(x,i),gsl_vector_get(y,i),gsl_vector_get(dy,i));}
+    fclose(data_points_stream);
+
+    FILE* fit_stream = fopen("fit.txt","w");
+    double fit_0 = gsl_vector_get(c,0);
+    double fit_1 = gsl_vector_get(c,1);
+    for(double k=gsl_vector_get(x,0);k<gsl_vector_get(x,8);k+=1.0/36) {
+        double fit = fit_0+fit_1*k;
+        double fit_minus = (fit_0-gsl_vector_get(dc,0))+(fit_1-gsl_vector_get(dc,0))*k;
+        double fit_plus = (fit_0+gsl_vector_get(dc,0))+(fit_1+gsl_vector_get(dc,0))*k;
+        fprintf(fit_stream,"%g %g %g %g\n",k,fit,fit_minus,fit_plus);
+    }
+    fclose(fit_stream);
+
+    FILE* values_stream = fopen("values.txt","w");
+    double half = -log(2)/gsl_vector_get(c,1);
+    fprintf(values_stream,"The half life of ThX from a least squares fit: %g days\nCompared to the modern value: 3.6319 days\n",half);
+    fprintf(values_stream,"Covariance matrix of fitting parameters:\n");
+    print_matrix(values_stream,cova);
+    fprintf(values_stream,"Meaning that the uncertainty in the lambda parameter is: %g\n",gsl_vector_get(dc,1));
+    double half_min = -log(2)/(gsl_vector_get(c,1)-gsl_vector_get(dc,1));
+    double half_max = -log(2)/(gsl_vector_get(c,1)+gsl_vector_get(dc,1));
+    fprintf(values_stream,"Thus the half life is between %g and %g days\n",half_min,half_max);
+    fclose(values_stream);
+
+
+    gsl_vector_free(x);
+    gsl_vector_free(y);
+    gsl_vector_free(dy);
+    gsl_vector_free(c);
+    gsl_vector_free(dc);
+    gsl_matrix_free(cova);
+return 0;
+}
\ No newline at end of file

homework/least-squares/values.txt

+The half life of ThX from a least squares fit: 4.04452 days
+Compared to the modern value: 3.6319 days
+Covariance matrix of fitting parameters:
+    0.000890278	      -8.75e-05	
+      -8.75e-05	       1.25e-05	
+Meaning that the uncertainty in the lambda parameter is: 0.00353553
+Thus the half life is between 3.96277 and 4.12971 days