minimization homework

homework/minimization/Makefile

+CFLAGS = -Wall -O3 -std=gnu11 $(shell gsl-config --cflags)
+LDLIBS = -lm $(shell gsl-config --libs)
+
+default: higgs.png simplex.png #animation_files
+
+higgs.png: BW_fit.txt Makefile
+	echo '\
+	set terminal png;\
+	set output "$@";\
+	set key top right;\
+	set tics out;\
+	set grid;\
+	set xlabel "Energy [GeV]";\
+	set ylabel "Sigma";\
+	set title "Higgs production signal";\
+	plot\
+	"BW_fit.txt" using 1:2 with line title "Fit",\
+	"CERN.txt" using 1:2:3 with yerrorbar title "CERN data",\
+	' | gnuplot
+
+simplex.png: path.txt contour.txt Makefile
+	echo '\
+	set terminal png;\
+	set output "$@";\
+	set key top right;\
+	set tics out;\
+	set grid;\
+	set xlabel "x";\
+	set ylabel "y";\
+	set title "Simplex steps taken while minimizing Rosenbrock";\
+	plot\
+	"path.txt" using 1:2 with line lw 3 title "Simplex",\
+	' | gnuplot
+
+NUMBERS = 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33\
+34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55
+animation_files: path.txt contour.txt folder Makefile
+	$(foreach var,$(NUMBERS), \
+		echo '\
+		set terminal png;\
+		set output "./animation/animation_file$(var).png";\
+		set key top right;\
+		set tics out;\
+		set grid;\
+		set xrange [-1.5:1.5];\
+		set yrange [-3:3];\
+		set xlabel "x";\
+		set ylabel "y";\
+		set title "Simplex steps taken while minimizing Rosenbrock";\
+		plot\
+		"path.txt" every :::$(var)::$(var) using 1:2 with line lw 3 title "Simplex",\
+		' | gnuplot ;)
+
+BW_fit.txt: main
+	./$< >out.txt
+
+path.txt: main
+	./$< >out.txt
+
+contour.txt: main
+	./$< >out.txt
+
+main: main.c functions.c
+
+folder:
+	mkdir animation
+
+.PHONEY: clean
+
+clean:
+	$(RM) main *.o *.txt *.png
+	$(RM) -r animation
\ No newline at end of file

homework/minimization/functions.c

+#include "functions.h"
+static double energy[30] = {101,103,105,107,109,111,113,115,117,119,121,123,125,127,129,\
+                            131,133,135,137,139,141,143,145,147,149,151,153,155,157,159};
+static double sigma[30] = {-0.25,-0.30,-0.15,-1.71,-0.81,0.65,-0.91,0.91,0.96,-2.52,-1.01,2.01,4.83,4.58,1.26,\
+                           1.01,-1.26,0.45,0.15,-0.91,-0.81,-1.41,1.36,0.50,-0.45,1.61,-2.21,-1.86,1.76,-0.50};
+static double dsigma[30] = {2.0,2.0,1.9,1.9,1.9,1.9,1.9,1.9,1.6,1.6,1.6,1.6,1.6,1.6,1.3,\
+                            1.3,1.3,1.3,1.3,1.3,1.1,1.1,1.1,1.1,1.1,1.1,1.1,0.9,0.9,0.9};
+
+double RB(gsl_vector* x) {
+    double x0 = gsl_vector_get(x,0);
+    double x1 = gsl_vector_get(x,1);
+    return pow((1-x0),2)+100*pow((x1-pow(x0,2)),2);
+}
+
+double HB(gsl_vector* x) {
+    double x0 = gsl_vector_get(x,0);
+    double x1 = gsl_vector_get(x,1);
+    return pow(pow(x0,2)+x1-11,2)+pow(x0+pow(x1,2)-7,2);
+}
+
+double F(double e, double m, double G, double A) {
+    return A/(pow(e-m,2)+pow(G,2)/4);
+}
+
+double BWd(gsl_vector* x) {
+    double m = gsl_vector_get(x,0);
+    double G = gsl_vector_get(x,1);
+    double A = gsl_vector_get(x,2);
+    double result = 0;
+    for(int i=0;i<30;i++) {
+        result += pow(F(energy[i],m,G,A)-sigma[i],2)/pow(dsigma[i],2);
+    }
+    return result;
+}
+
+void data(int i, double* e, double* sig, double* err) {
+    *e = energy[i]; *sig = sigma[i]; *err = dsigma[i];
+}
+
+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(gsl_matrix* M) {
+    for(int i = 0; i < M->size1; i++) {
+        for(int j = 0; j < M->size2; j++) {
+            printf("%10g\t", gsl_matrix_get(M, i, j));
+        }
+        printf("\n");
+    }
+}
+
+int qnewton(double f(gsl_vector* x), gsl_vector* x, double eps) {
+    int dim = x->size;
+    gsl_matrix* B = gsl_matrix_alloc(dim,dim);
+    gsl_vector* nx = gsl_vector_alloc(dim);
+    gsl_vector* nxps = gsl_vector_alloc(dim);
+    gsl_vector* dx = gsl_vector_alloc(dim);
+    gsl_matrix_set_identity(B);
+    double new_gradient;
+    int reset = 0, steps = 0;
+    do{
+        if(reset>0) {gsl_matrix_set_identity(B); reset=0;}
+        double fx = f(x);
+        //Calculate numerical derivatives at x
+        for(int i=0;i<dim;i++) {
+            double* xi = gsl_vector_ptr(x,i);
+            *xi+=DX; double dfx = f(x); *xi-=DX;
+            gsl_vector_set(nx,i,(dfx-fx)/DX);
+        }
+        //Calculate optimal delta x
+        gsl_blas_dgemv(CblasNoTrans,-1,B,nx,0,dx);
+        gsl_blas_daxpy(1,dx,x);
+        double l=1;
+        double a_cond;
+        gsl_blas_ddot(dx,nx,&a_cond); a_cond *= 1e-4; a_cond += fx;
+        while(f(x)>a_cond && l>DX) {
+            l /= 2;
+            gsl_vector_scale(dx,1.0/2);
+            gsl_blas_ddot(dx,nx,&a_cond); a_cond *= 1e-4; a_cond += fx;
+            gsl_blas_daxpy(-1,dx,x);
+        }
+        if(l<DX) {reset++;}
+        //Update Hessian matrix
+        for(int i=0;i<dim;i++) {
+            double* xi = gsl_vector_ptr(x,i);
+            double ndfx = f(x); *xi+=DX;
+            double dfx = f(x); *xi-=DX;
+            gsl_vector_set(nxps,i,(dfx-ndfx)/DX);
+        }
+        new_gradient = gsl_blas_dnrm2(nxps);
+        gsl_blas_daxpy(-1,nx,nxps); gsl_vector_memcpy(nx,dx);
+        gsl_blas_dgemv(CblasNoTrans,-1,B,nxps,1,nx);
+        gsl_blas_ddot(dx,nxps,&a_cond); //a_cond is used because it isnt currently in use
+        if(a_cond<1e-6) {a_cond=0; reset++;} //only update if dot product is reasonable
+        gsl_blas_dger(1.0/a_cond,nx,dx,B);
+        fx = f(x); steps++;
+        if(steps>=1e6) {printf("Maximum number of steps (1e6) reached.\n"); break;}
+    } while(new_gradient>eps);
+    gsl_matrix_free(B);
+    gsl_vector_free(nx);
+    gsl_vector_free(nxps);
+    gsl_vector_free(dx);
+return steps;
+}
+
+int simplex(double f(gsl_vector* x), gsl_matrix* S, gsl_vector* result, double eps) {
+    int dim = S->size1, verts = dim+1, steps = 0;
+    long unsigned int min_index, max_index;
+    double min_val, max_val, size = 0;
+    gsl_vector* vert_vals = gsl_vector_alloc(verts);
+    gsl_vector* pce = gsl_vector_alloc(dim);
+    gsl_vector* change = gsl_vector_alloc(dim);
+    //FILE* path_stream = fopen("path.txt","w");
+    do {
+        /*Printing current vertices in a form that Gnuplot can plot
+        for(int j=0;j<verts;j++) {
+            for(int i=0;i<dim;i++) {
+                fprintf(path_stream,"%g ",gsl_matrix_get(S,i,j));
+            } fprintf(path_stream,"\n");
+        } 
+        for(int i=0;i<dim;i++) {fprintf(path_stream,"%g ",gsl_matrix_get(S,i,0));} fprintf(path_stream,"\n\n");
+        */
+        steps++;
+        gsl_vector_set_zero(pce);
+        for(int j=0;j<verts;j++) {
+            gsl_vector_view temp = gsl_matrix_column(S,j);
+            gsl_vector* vertj = &(temp.vector);
+            gsl_blas_daxpy(1,vertj,pce); //Sums all vertices! Still needs to subtract phi and /n to get final pce
+            gsl_vector_set(vert_vals,j,f(vertj)); //Function value at j'th vertex
+        }
+        //printf("S before step:\n"); //For debugging
+        //print_matrix(S); //For debugging
+        gsl_vector_minmax_index(vert_vals,&min_index,&max_index); 
+        min_val = gsl_vector_get(vert_vals,min_index); max_val = gsl_vector_get(vert_vals,max_index);
+        gsl_vector_view plo_view = gsl_matrix_column(S,min_index); gsl_vector* plo = &(plo_view.vector);
+        gsl_vector_view phi_view = gsl_matrix_column(S,max_index); gsl_vector* phi = &(phi_view.vector);
+        gsl_blas_daxpy(-1,phi,pce); gsl_vector_scale(pce,1.0/dim); //We now got plo, phi and pce
+
+        gsl_vector_memcpy(change,pce); gsl_blas_daxpy(-1,phi,change); gsl_blas_daxpy(2,change,phi); //reflection
+        double reflec_val = f(phi);
+        if(reflec_val<min_val) {
+            gsl_blas_daxpy(1,change,phi); //Expansion
+            if(f(phi)>reflec_val) {gsl_blas_daxpy(-1,change,phi);} //Accept reflection
+        }
+        else {
+            if(reflec_val>max_val) {
+                gsl_blas_daxpy(-1.0/2,change,phi); //Contraction
+                if(f(phi)>max_val) {
+                    gsl_blas_daxpy(-3.0/2,change,phi); //Reset phi
+                    for(int i = 0;i < verts; i++) { //Reduction
+                        if(i!=min_index) {
+                            gsl_vector_view Si_view = gsl_matrix_column(S,i); gsl_vector* Si = &(Si_view.vector);
+                            gsl_vector_memcpy(change,plo); gsl_blas_daxpy(-1,Si,change); //Change is vector Si->plo
+                            gsl_blas_daxpy(1.0/2,change,Si);
+                        }
+                    }
+                }
+            }
+        }
+        //printf("S after step:\n"); //For debugging
+        //print_matrix(S); //For debugging
+        size=0; //Reset size
+        for(int i=0;i<verts;i++) {
+            for(int j=i+1;j<verts;j++) {
+                gsl_vector_view Si_view = gsl_matrix_column(S,i); gsl_vector* Si = &(Si_view.vector);
+                gsl_vector_view Sj_view = gsl_matrix_column(S,j); gsl_vector* Sj = &(Sj_view.vector);
+                gsl_vector* norm_calc = gsl_vector_alloc(dim); gsl_vector_memcpy(norm_calc,Si);
+                gsl_blas_daxpy(-1,Sj,norm_calc); double sizei = gsl_blas_dnrm2(norm_calc); //Calculate norm of Si->Sj
+                if(sizei>size) {size=sizei;} //Update size, so that size equals largest distance between 2 vertices
+                gsl_vector_free(norm_calc);
+            }
+        }
+        //printf("Size: %g\n",size); //For debugging
+    } while(size>eps && steps<1e6);
+    //fclose(path_stream);
+    gsl_vector_view min_view = gsl_matrix_column(S,min_index); gsl_vector* min_vec = &(min_view.vector);
+    gsl_vector_memcpy(result,min_vec); //Puts final min vector in result
+    gsl_vector_free(vert_vals);
+    gsl_vector_free(pce);
+    gsl_vector_free(change);
+return steps;
+}
\ No newline at end of file

homework/minimization/functions.h

+#ifndef FUNCTIONS_H
+#define FUNCTIONS_H
+
+#include<math.h>
+#include<time.h>
+#include<stdio.h>
+#include<float.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
+#define DX sqrt(DBL_EPSILON)
+
+double RB(gsl_vector* x);
+double HB(gsl_vector* x);
+double F(double e, double m, double G, double A);
+double BWd(gsl_vector* x);
+void data(int i, double* e, double* sig, double* err);
+void print_vector(char s[], gsl_vector* v);
+void print_matrix(gsl_matrix* M);
+int qnewton(double f(gsl_vector* x), gsl_vector* x, double eps);
+int simplex(double f(gsl_vector* x), gsl_matrix* S, gsl_vector* result, double eps);
+
+#endif
\ No newline at end of file

homework/minimization/main.c

+#include "functions.h"
+
+int main() {
+
+    gsl_vector* x = gsl_vector_alloc(2);
+    gsl_vector_set(x,0,15);
+    gsl_vector_set(x,1,18);
+    int steps;
+
+    print_vector("Initial x for minimization of Rosenbrock:",x);
+    steps = qnewton(RB,x,1e-4);
+    printf("Minimization took %i steps. ",steps);
+    print_vector("X after minimization:",x); printf("\n\n");
+
+    gsl_vector_set(x,0,0);
+    gsl_vector_set(x,1,0);
+    print_vector("Initial x for minimization of Himmelblau:",x);
+    steps = qnewton(HB,x,1e-6);
+    printf("Minimization took %i steps. ",steps);
+    print_vector("X after minimization:",x); printf("\n\n");
+
+    FILE* plot_stream = fopen("BW_fit.txt","w");
+    gsl_vector* mGA = gsl_vector_alloc(3);
+    gsl_vector_set(mGA,0,110);
+    gsl_vector_set(mGA,1,10);
+    gsl_vector_set(mGA,2,10);
+    print_vector("Initial m G and A for minimization of deviation function:",mGA);
+    steps = qnewton(BWd,mGA,1e-3);
+    printf("Minimization took %i steps. ",steps);
+    print_vector("m G and A after minimization:",mGA); printf("\n\n");
+    for(double e=101;e<=159;e+=1.0/32) {
+        fprintf(plot_stream,"%g %g\n",e,F(e,gsl_vector_get(mGA,0),gsl_vector_get(mGA,1),gsl_vector_get(mGA,2)));
+    }
+    gsl_vector_free(x);
+    gsl_vector_free(mGA);
+    fclose(plot_stream);
+
+    FILE* data_stream = fopen("CERN.txt","w");
+    double e, sig, err;
+    for(int i=0;i<30;i++) {
+        data(i,&e,&sig,&err);
+        fprintf(data_stream,"%g %g %g\n",e,sig,err);
+    }
+    fclose(data_stream);
+
+    gsl_matrix* S = gsl_matrix_alloc(2,3);
+    gsl_vector* result = gsl_vector_alloc(2);
+    gsl_matrix_set(S,0,0,0); gsl_matrix_set(S,0,1,-3.0/2); gsl_matrix_set(S,0,2,3.0/2);
+    gsl_matrix_set(S,1,0,3); gsl_matrix_set(S,1,1,-1.0/2); gsl_matrix_set(S,1,2,-1.0/2);
+    printf("Initial simplex for downhill finding of minimum in Rosenbrock:\n"); print_matrix(S);
+    steps = simplex(RB,S,result,1e-6);
+    printf("Minimization took %i steps. ",steps);
+    print_vector("Found minimum:",result); printf("\n\n");
+    gsl_matrix_free(S);
+    gsl_vector_free(result);
+    
+    S = gsl_matrix_alloc(3,4);
+    result = gsl_vector_alloc(3);
+    gsl_matrix_set(S,0,0,0); gsl_matrix_set(S,0,1,200); gsl_matrix_set(S,0,2,0); gsl_matrix_set(S,0,3,200);
+    gsl_matrix_set(S,1,0,0); gsl_matrix_set(S,1,1,0); gsl_matrix_set(S,1,2,200); gsl_matrix_set(S,1,3,200);
+    gsl_matrix_set(S,2,0,0); gsl_matrix_set(S,2,1,0); gsl_matrix_set(S,2,2,0); gsl_matrix_set(S,2,3,200);
+    printf("Initial simplex for downhill finding of minimum in deviation function:\n"); print_matrix(S);
+    steps = simplex(BWd,S,result,1e-4);
+    printf("Minimization took %i steps. ",steps);
+    print_vector("Found minimum:",result); printf("\n\n");
+    gsl_matrix_free(S);
+    gsl_vector_free(result);
+    
+return 0;
+}
\ No newline at end of file

homework/minimization/out.txt

+Initial x for minimization of Rosenbrock:
+             15
+             18
+Minimization took 129 steps. X after minimization:
+       0.999995
+       0.999991
+
+
+Initial x for minimization of Himmelblau:
+              0
+              0
+Minimization took 26 steps. X after minimization:
+              3
+              2
+
+
+Initial m G and A for minimization of deviation function:
+            110
+             10
+             10
+Minimization took 888 steps. m G and A after minimization:
+        125.973
+        2.06333
+        9.76338
+
+
+Initial simplex for downhill finding of minimum in Rosenbrock:
+         0	      -1.5	       1.5	
+         3	      -0.5	      -0.5	
+Minimization took 100 steps. Found minimum:
+              1
+              1
+Ive made a gif of this process :)
+
+
+Initial simplex for downhill finding of minimum in deviation function:
+         0	       200	         0	       200	
+         0	         0	       200	       200	
+         0	         0	         0	       200	
+Minimization took 169 steps. Found minimum:
+        125.973
+        2.05282
+        9.71236
+
+