Skip to content

Commit 45420e4

Browse files
devriesdevries
committed
Initial conversion from CVS to Mercurial
0 parents  commit 45420e4

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

51 files changed

+4816
-0
lines changed

‎AUTHORS

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1 @@
1+
Christopher De Vries ([email protected])

‎ChangeLog

Whitespace-only changes.

‎Makefile.ORIG

Lines changed: 100 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,100 @@
1+
CC = gcc
2+
LIBS = -lgsl -lgslcblas -lm
3+
4+
OBJS = grow_array.o devo2.o read_data.o
5+
6+
CFLAGS = -g -O -Wall
7+
8+
DESTDIR = /home/cdevries/bin
9+
10+
all: twolayer5_devo twolayer5_simplex twolayer5_hybrid hill5_devo hill5_simplex hill5_hybrid twolayer6_devo twolayer6_simplex twolayer6_hybrid hill7_devo hill7_simplex hill7_hybrid hill6core_devo hill6core_simplex hill6core_hybrid hill6_devo hill6_simplex hill6_hybrid
11+
cp noisescript.py $(DESTDIR)
12+
13+
test: ga_test devo2_test read_data_test
14+
15+
ga_test: ga_test.o $(OBJS)
16+
$(CC) $(CFLAGS) ga_test.o $(OBJS) $(LIBS) -o $@
17+
18+
devo2_test: devo2_test.o $(OBJS)
19+
$(CC) $(CFLAGS) devo2_test.o $(OBJS) $(LIBS) -o $@
20+
21+
read_data_test: read_data_test.o $(OBJS)
22+
$(CC) $(CFLAGS) read_data_test.o $(OBJS) $(LIBS) -o $@
23+
24+
twolayer5_devo: twolayer5_devo.o twolayer5.o $(OBJS)
25+
$(CC) $(CFLAGS) twolayer5_devo.o twolayer5.o $(OBJS) $(LIBS) -o $@
26+
mv $@ $(DESTDIR)
27+
28+
twolayer5_simplex: twolayer5_simplex.o twolayer5.o $(OBJS)
29+
$(CC) $(CFLAGS) twolayer5_simplex.o twolayer5.o $(OBJS) $(LIBS) -o $@
30+
mv $@ $(DESTDIR)
31+
32+
twolayer5_hybrid: twolayer5_hybrid.o twolayer5.o $(OBJS)
33+
$(CC) $(CFLAGS) twolayer5_hybrid.o twolayer5.o $(OBJS) $(LIBS) -o $@
34+
mv $@ $(DESTDIR)
35+
36+
hill5_devo: hill5_devo.o hill5.o $(OBJS)
37+
$(CC) $(CFLAGS) hill5_devo.o hill5.o $(OBJS) $(LIBS) -o $@
38+
mv $@ $(DESTDIR)
39+
40+
hill5_simplex: hill5_simplex.o hill5.o $(OBJS)
41+
$(CC) $(CFLAGS) hill5_simplex.o hill5.o $(OBJS) $(LIBS) -o $@
42+
mv $@ $(DESTDIR)
43+
44+
hill5_hybrid: hill5_hybrid.o hill5.o $(OBJS)
45+
$(CC) $(CFLAGS) hill5_hybrid.o hill5.o $(OBJS) $(LIBS) -o $@
46+
mv $@ $(DESTDIR)
47+
48+
twolayer6_devo: twolayer6_devo.o twolayer6.o $(OBJS)
49+
$(CC) $(CFLAGS) twolayer6_devo.o twolayer6.o $(OBJS) $(LIBS) -o $@
50+
mv $@ $(DESTDIR)
51+
52+
twolayer6_simplex: twolayer6_simplex.o twolayer6.o $(OBJS)
53+
$(CC) $(CFLAGS) twolayer6_simplex.o twolayer6.o $(OBJS) $(LIBS) -o $@
54+
mv $@ $(DESTDIR)
55+
56+
twolayer6_hybrid: twolayer6_hybrid.o twolayer6.o $(OBJS)
57+
$(CC) $(CFLAGS) twolayer6_hybrid.o twolayer6.o $(OBJS) $(LIBS) -o $@
58+
mv $@ $(DESTDIR)
59+
60+
hill7_devo: hill7_devo.o hill7.o $(OBJS)
61+
$(CC) $(CFLAGS) hill7_devo.o hill7.o $(OBJS) $(LIBS) -o $@
62+
mv $@ $(DESTDIR)
63+
64+
hill7_simplex: hill7_simplex.o hill7.o $(OBJS)
65+
$(CC) $(CFLAGS) hill7_simplex.o hill7.o $(OBJS) $(LIBS) -o $@
66+
mv $@ $(DESTDIR)
67+
68+
hill7_hybrid: hill7_hybrid.o hill7.o $(OBJS)
69+
$(CC) $(CFLAGS) hill7_hybrid.o hill7.o $(OBJS) $(LIBS) -o $@
70+
mv $@ $(DESTDIR)
71+
72+
hill6core_devo: hill6core_devo.o hill6core.o $(OBJS)
73+
$(CC) $(CFLAGS) hill6core_devo.o hill6core.o $(OBJS) $(LIBS) -o $@
74+
mv $@ $(DESTDIR)
75+
76+
hill6core_simplex: hill6core_simplex.o hill6core.o $(OBJS)
77+
$(CC) $(CFLAGS) hill6core_simplex.o hill6core.o $(OBJS) $(LIBS) -o $@
78+
mv $@ $(DESTDIR)
79+
80+
hill6core_hybrid: hill6core_hybrid.o hill6core.o $(OBJS)
81+
$(CC) $(CFLAGS) hill6core_hybrid.o hill6core.o $(OBJS) $(LIBS) -o $@
82+
mv $@ $(DESTDIR)
83+
84+
hill6_devo: hill6_devo.o hill6.o $(OBJS)
85+
$(CC) $(CFLAGS) hill6_devo.o hill6.o $(OBJS) $(LIBS) -o $@
86+
mv $@ $(DESTDIR)
87+
88+
hill6_simplex: hill6_simplex.o hill6.o $(OBJS)
89+
$(CC) $(CFLAGS) hill6_simplex.o hill6.o $(OBJS) $(LIBS) -o $@
90+
mv $@ $(DESTDIR)
91+
92+
hill6_hybrid: hill6_hybrid.o hill6.o $(OBJS)
93+
$(CC) $(CFLAGS) hill6_hybrid.o hill6.o $(OBJS) $(LIBS) -o $@
94+
mv $@ $(DESTDIR)
95+
96+
.c.o:
97+
$(CC) $(CFLAGS) -c $*.c
98+
99+
clean:
100+
rm -f *.o ga_test devo2_test read_data_test

‎Makefile.am

Lines changed: 23 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,23 @@
1+
bin_PROGRAMS = twolayer5_devo twolayer5_simplex twolayer5_hybrid hill5_devo hill5_simplex hill5_hybrid twolayer6_devo twolayer6_simplex twolayer6_hybrid hill6_devo hill6_simplex hill6_hybrid hill6core_devo hill6core_simplex hill6core_hybrid hill7_devo hill7_simplex hill7_hybrid hill5map hill5_plot
2+
twolayer5_devo_SOURCES = twolayer5_devo.c twolayer5.c twolayer5.h devo2.h devo2.c grow_array.h grow_array.c read_data.c
3+
twolayer5_simplex_SOURCES = twolayer5_simplex.c twolayer5.c twolayer5.h grow_array.h grow_array.c read_data.c
4+
twolayer5_hybrid_SOURCES = twolayer5_hybrid.c twolayer5.c twolayer5.h devo2.h devo2.c grow_array.h grow_array.c read_data.c
5+
hill5_devo_SOURCES = hill5_devo.c hill5.c hill5.h devo2.h devo2.c grow_array.h grow_array.c read_data.c
6+
hill5_simplex_SOURCES = hill5_simplex.c hill5.c hill5.h grow_array.h grow_array.c read_data.c
7+
hill5_hybrid_SOURCES = hill5_hybrid.c hill5.c devo2.h devo2.c grow_array.h grow_array.c read_data.c
8+
twolayer6_devo_SOURCES = twolayer6_devo.c twolayer6.c twolayer6.h devo2.h devo2.c grow_array.h grow_array.c read_data.c
9+
twolayer6_simplex_SOURCES = twolayer6_simplex.c twolayer6.c twolayer6.h grow_array.h grow_array.c read_data.c
10+
twolayer6_hybrid_SOURCES = twolayer6_hybrid.c twolayer6.c twolayer6.h devo2.h devo2.c grow_array.h grow_array.c read_data.c
11+
hill6_devo_SOURCES = hill6_devo.c hill6.c hill6.h devo2.h devo2.c grow_array.h grow_array.c read_data.c
12+
hill6_simplex_SOURCES = hill6_simplex.c hill6.c hill6.h grow_array.h grow_array.c read_data.c
13+
hill6_hybrid_SOURCES = hill6_hybrid.c hill6.c hill6.h devo2.h devo2.c grow_array.h grow_array.c read_data.c
14+
hill6core_devo_SOURCES = hill6core_devo.c hill6core.c hill6core.h devo2.h devo2.c grow_array.h grow_array.c read_data.c
15+
hill6core_simplex_SOURCES = hill6core_simplex.c hill6core.c hill6core.h grow_array.h grow_array.c read_data.c
16+
hill6core_hybrid_SOURCES = hill6core_hybrid.c hill6core.c hill6core.h devo2.h devo2.c grow_array.h grow_array.c read_data.c
17+
hill7_devo_SOURCES = hill7_devo.c hill7.c hill7.h devo2.h devo2.c grow_array.h grow_array.c read_data.c
18+
hill7_simplex_SOURCES = hill7_simplex.c hill7.c hill7.h grow_array.h grow_array.c read_data.c
19+
hill7_hybrid_SOURCES = hill7_hybrid.c hill7.c hill7.h devo2.h devo2.c grow_array.h grow_array.c read_data.c
20+
hill5map_SOURCES = hill5map.c hill5.c hill5.h devo2.h devo2.c grow_array.h grow_array.c read_data.c
21+
hill5_plot_SOURCES = hill5_plot.c hill5.c hill5.h grow_array.h grow_array.c read_data.c
22+
bin_SCRIPTS = noisescript.py noisescript2.py
23+
EXTRA_DIST = $(bin_SCRIPTS)

‎NEWS

Whitespace-only changes.

‎README

Lines changed: 53 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,53 @@
1+
Thank you for downloading the analytic infall models. These models are
2+
discussed in detail in De Vries & Myers (2005). More documentation and updates
3+
are available from http://cfa-www.harvard.edu/~cdevries/analytic_infall.html
4+
. Although this paper includes all the analytic models discussed in the paper,
5+
I would restrict myself to the HILL5 and TWOLAYER6 models, which were found to
6+
be the best models.
7+
8+
In order to fit a spectrum you need to place it in a two column text format
9+
(the first column should be velocity in kilmeters per second, and the second
10+
should be brightness temperature in Kelvins). Lines that begin with a hash mark
11+
'#' will be ignored. Use the hybrid version of the model you want to run, this
12+
will allow for some slow differential evolution minimization which will
13+
(hopefully) find the global minimum well followed by Nelder-Mead simplex
14+
minimization which will quickly find the bottom of that well. The hybrid
15+
programs of interest are "hill5_hybrid" and "twolayer6_hybrid". The arguments
16+
to both programs are:
17+
18+
1. input filename --- the file with your spectrum
19+
2. frequency --- to perform the J(T) conversion. A frequency of 0 will work in
20+
units of J(T) instead of T.
21+
3. vmin --- minimum velocity of the line profile.
22+
4. vmax --- maximum velocity of the line profile. Be sure to get the entire
23+
line full width as well as a little baseline. This is the region
24+
over which chi-squared is calculated, so artifacts of the fit can
25+
sometimes appear outside this region.
26+
5. population in generation --- The number of solutions to calculate each
27+
generation of the differential evolution. I
28+
usually use about 200-300 here.
29+
6. generations per check --- The number of generations to run before checking
30+
for convergence in the differential evolution. I
31+
usually pick 300 here.
32+
7. checks to convergence --- The number of checks to make before deciding the
33+
differential evolution algorithm has converged. I
34+
usually pick 3 to 5 here.
35+
8. output file --- Where to place the fit and parameters.
36+
37+
The output file will have a header line the following (a HILL 5 example is
38+
shown):
39+
40+
# Tau: 4.46394
41+
# Vlsr: 0.000296732
42+
# Vin: 0.0993918
43+
# sigma: 0.0937655
44+
# Tpeak: 10.5371
45+
# Attained Chisq: 0.539567
46+
47+
giving the parameters of the fit and Chi Squared (not a reduced
48+
chi-squared). Followed by the fit in a three column data file. The first two
49+
columns are the two columns of the input file, while the third column is the
50+
fit of the analytic model. Always check that the fit is good.
51+
52+
Please email me at [email protected] if you have any questions or
53+
comments.

‎configure.in

Lines changed: 25 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,25 @@
1+
dnl Process this file with autoconf to produce a configure script.
2+
AC_INIT(devo2.c)
3+
AM_INIT_AUTOMAKE(analytic_infall,1.1)
4+
5+
dnl Checks for programs.
6+
AC_PROG_CC
7+
8+
dnl Checks for libraries.
9+
AC_CHECK_LIB(m,main)
10+
AC_CHECK_LIB(cfitsio,main)
11+
AC_CHECK_LIB(gslcblas,main)
12+
AC_CHECK_LIB(gsl,main,[AC_DEFINE([HAVE_LIBGSL]) LIBS="-lgsl $LIBS"],AC_ERROR([Gnu Scientific Library not found]))
13+
14+
AM_PATH_GSL(1.3, ,AC_ERROR([Version 1.3 or greater of Gnu Scientific Library Required]))
15+
16+
dnl Checks for header files.
17+
AC_HEADER_STDC
18+
19+
dnl Checks for typedefs, structures, and compiler characteristics.
20+
AC_C_CONST
21+
AC_TYPE_SIZE_T
22+
23+
dnl Checks for library functions.
24+
AC_PROG_INSTALL
25+
AC_OUTPUT(Makefile)

‎devo2.c

Lines changed: 124 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,124 @@
1+
#include "devo2.h"
2+
#include <stdlib.h>
3+
#include <stdio.h>
4+
#include <math.h>
5+
#include <float.h>
6+
7+
void devo2_init(devo2_struct *st, int num_param, double *min_param, double *max_param, int num_vecs, double mutation_prob, double diff_factor, double (*evaluate)(double *param_vector)) {
8+
int i, j;
9+
10+
/* Check num_vecs Argument */
11+
if(num_vecs<4) {
12+
fprintf(stderr, "Too few parameter vectors to evolve.\n");
13+
exit(1);
14+
}
15+
16+
/* Assign some basic values */
17+
st->num_param=num_param;
18+
st->num_vecs=num_vecs;
19+
st->mutation_prob = mutation_prob;
20+
st->diff_factor = diff_factor;
21+
st->gen_count=1;
22+
st->minimization_function = evaluate;
23+
st->best_vector=NULL;
24+
st->best_score=DBL_MAX;
25+
26+
/* Allocate memory and assign pointers */
27+
st->param_matrix1=(double*)malloc(num_param*num_vecs*sizeof(double));
28+
st->param_matrix2=(double*)malloc(num_param*num_vecs*sizeof(double));
29+
st->scores = (double*)malloc(num_vecs*sizeof(double));
30+
st->trial = (double*)malloc(num_param*sizeof(double));
31+
st->prev_gen=st->param_matrix1;
32+
st->next_gen=st->param_matrix2;
33+
34+
/* Check to see if memory allocation succeeded */
35+
if(st->param_matrix1==NULL || st->param_matrix2==NULL || st->scores==NULL || st->trial==NULL) {
36+
fprintf(stderr, "ERROR in devo2_init: Unable to allocate sufficient memory\n");
37+
exit(1);
38+
}
39+
40+
/* Initialize Random Number Generator */
41+
gsl_rng_env_setup();
42+
st->rng = gsl_rng_alloc(gsl_rng_default);
43+
44+
/* Initialize parent generation */
45+
for(i=0;i<num_vecs;i++) {
46+
for(j=0;j<num_param;j++) {
47+
st->prev_gen[i*num_param+j]=gsl_rng_uniform(st->rng)*(max_param[j]-min_param[j])+min_param[j];
48+
}
49+
st->scores[i]=(*st->minimization_function)(&(st->prev_gen[i*num_param]));
50+
if(st->scores[i]<st->best_score) {
51+
st->best_score = st->scores[i];
52+
st->best_vector = &(st->prev_gen[i*num_param]);
53+
}
54+
}
55+
56+
return;
57+
}
58+
59+
void devo2_step(devo2_struct *st) {
60+
/* Begin Differential Evolution */
61+
int i, j, k;
62+
int a, b, c;
63+
double score;
64+
double *temp;
65+
66+
for(i=0;i<st->num_vecs;i++) {
67+
/* Select other parent vectors for mutation */
68+
do a=(int)gsl_rng_uniform_int(st->rng,(long int)st->num_vecs); while(a==i);
69+
do b=(int)gsl_rng_uniform_int(st->rng,(long int)st->num_vecs); while(b==i || b==a);
70+
do c=(int)gsl_rng_uniform_int(st->rng,(long int)st->num_vecs); while(c==i || c==a || c==b);
71+
72+
/* Randomly pick parameter to adjust */
73+
j=(int)gsl_rng_uniform_int(st->rng,(long int)st->num_param);
74+
75+
/* Create a Child Parameter Vector */
76+
for(k=1;k<=st->num_param;k++) {
77+
if(gsl_rng_uniform(st->rng) < st->mutation_prob || k==st->num_param) {
78+
/* Mutate parameter */
79+
st->trial[j]=st->prev_gen[c*st->num_param+j]+st->diff_factor*(st->prev_gen[a*st->num_param+j]-st->prev_gen[b*st->num_param+j]);
80+
}
81+
else {
82+
/* Copy parameter */
83+
st->trial[j]=st->prev_gen[i*st->num_param+j];
84+
}
85+
j=(j+1)%st->num_param; /* Started on random param, must use mod */
86+
}
87+
88+
/* Evaluate Child, and retain either child or parent or exit */
89+
score=(*st->minimization_function)(st->trial);
90+
91+
/* Evaluate whether to pass parent or child to next generation */
92+
if(score<=st->scores[i]) {
93+
for(j=0;j<st->num_param;j++) {
94+
st->next_gen[i*st->num_param+j]=st->trial[j];
95+
}
96+
st->scores[i]=score;
97+
/* Check if it is best score */
98+
if(score<st->best_score) {
99+
st->best_score = score;
100+
st->best_vector = &(st->next_gen[i*st->num_param]);
101+
}
102+
}
103+
else {
104+
for(j=0;j<st->num_param;j++) {
105+
st->next_gen[i*st->num_param+j]=st->prev_gen[i*st->num_param+j];
106+
}
107+
}
108+
}
109+
110+
/* Swap Generations */
111+
temp=st->next_gen;
112+
st->next_gen=st->prev_gen;
113+
st->prev_gen=temp;
114+
st->gen_count++;
115+
return;
116+
}
117+
118+
void devo2_free(devo2_struct *st) {
119+
free(st->param_matrix1);
120+
free(st->param_matrix2);
121+
free(st->scores);
122+
free(st->trial);
123+
gsl_rng_free(st->rng);
124+
}

0 commit comments

Comments
 (0)