diff --git a/src/diffevolution.c b/src/diffevolution.c index 81c7b80..b9420aa 100644 --- a/src/diffevolution.c +++ b/src/diffevolution.c @@ -10,7 +10,8 @@ #include #include #include -#include "omp.h" +#include + #include "eem.h" #include "kappa.h" #include "neemp.h" @@ -38,55 +39,58 @@ void run_diff_evolution(struct subset * const ss) { printf("DE Generating population of size %d\n", s.population_size); /* Set bounds for each parameters in kappa_data */ - float *bounds = (float*) malloc((ts.atom_types_count*2+1)*2*sizeof(float)); - //compute bounds, 0 means set them to fixed numbers taken from Tomas's full scan, 1 means try to find them with broad search + float *bounds = (float *) malloc((ts.atom_types_count * 2 + 1) * 2 * sizeof(float)); + /* Compute bounds, 0 means set them to fixed numbers taken from Tomas's full scan, 1 means try to find them with broad search */ compute_parameters_bounds(bounds, 0); - //generate population + /* Generate population */ fill_ss(ss, s.population_size); generate_random_population(ss, bounds, s.population_size); de_ss = ss; - //evaluate the fitness function for all points + /* Evaluate the fitness function for all points */ if (s.verbosity >= VERBOSE_KAPPA) printf("DE Calculating charges and evaluating fitness function for whole population\n"); int i = 0; -#pragma omp parallel for num_threads(s.de_threads) default(shared) private(i) + #pragma omp parallel for num_threads(s.de_threads) default(shared) private(i) for (i = 0; i < ss->kappa_data_count; i++) { calculate_charges(ss, &ss->data[i]); calculate_statistics(ss, &ss->data[i]); } - //minimize part of population + /* Minimize part of population */ int minimized_initial = 0; - int* good_indices = (int*) calloc(s.population_size, sizeof(int)); + int *good_indices = (int *) calloc(s.population_size, sizeof(int)); for (i = 0; i < s.population_size; i++) good_indices[i] = -1; + if (s.polish > 2) { if (s.verbosity >= VERBOSE_KAPPA) printf("DE minimizing part of population\n"); minimized_initial = minimize_part_of_population(ss, good_indices); } - //if we minimized zero data or polish <= 2, use all kappa_data instead of minimized + + /* If we minimized zero data or polish <= 2, use all kappa_data instead of minimized */ if (minimized_initial == 0) { for (i = 0; i < s.population_size; i++) good_indices[i] = i; + minimized_initial = s.population_size; } - //find the best kappa_data + /* find the best kappa_data */ set_the_best(ss); /* Run the optimization for max iterations */ - //TODO include iters_max for DE in limits or use one already there (that is used for discard) - //TODO can we tell we have converged? if yes, include to while condition + /* TODO include iters_max for DE in limits or use one already there (that is used for discard) */ + /* TODO can we tell we have converged? if yes, include to while condition */ int iter = 0; - struct kappa_data *trial = (struct kappa_data*)malloc(sizeof(struct kappa_data)); - struct kappa_data *so_far_best = (struct kappa_data*)malloc(sizeof(struct kappa_data)); + struct kappa_data *trial = (struct kappa_data *) malloc(sizeof(struct kappa_data)); + struct kappa_data *so_far_best = (struct kappa_data *) malloc(sizeof(struct kappa_data)); kd_init(trial); trial->parent_subset = ss; kd_init(so_far_best); - //copy ss->best into so_far_best + /* copy ss->best into so_far_best */ kd_copy_parameters(ss->best, so_far_best); calculate_charges(ss, so_far_best); calculate_statistics(ss, so_far_best); @@ -96,17 +100,17 @@ void run_diff_evolution(struct subset * const ss) { calculate_statistics(ss, so_far_best); kd_print_results(so_far_best); float mutation_constant = s.mutation_constant; - int iters_with_evolution=0; - int condition=1; + int iters_with_evolution = 0; + int condition = 1; int minimized = 0; -#pragma omp parallel num_threads(s.de_threads) default(shared) + #pragma omp parallel num_threads(s.de_threads) default(shared) { while (condition) { -#pragma omp master + #pragma omp master condition = (iter < s.limit_de_iters); { -#pragma omp master + #pragma omp master { iter++; if (s.verbosity >= VERBOSE_KAPPA) { @@ -124,17 +128,20 @@ void run_diff_evolution(struct subset * const ss) { if (s.dither) mutation_constant = get_random_float(0.5, 1); + /* Recombine parts of best, a and b to obtain new trial structure */ -#pragma omp critical + #pragma omp critical evolve_kappa(trial, so_far_best, a, b, bounds, mutation_constant, s.recombination_constant); iters_with_evolution++; + /* Evaluate the new trial structure */ calculate_charges(ss, trial); calculate_statistics(ss, trial); if (s.verbosity >= VERBOSE_KAPPA) printf("Trial stats %f %f %d\n", trial->full_stats.R_w, trial->full_stats.R2, is_quite_good(trial)); + /* If the new structure is better than what we have before, reassign */ -#pragma omp critical + #pragma omp critical { if (compare_and_set(trial, so_far_best)) { calculate_charges(ss, so_far_best); @@ -146,25 +153,30 @@ void run_diff_evolution(struct subset * const ss) { } } } + /* All other threads do this */ if (s.polish > 1 && is_quite_good(trial) && (s.de_threads == 0 || omp_get_thread_num() != 0)) { minimized++; if (s.verbosity >= VERBOSE_KAPPA) printf("\nDE min thread %d\n", omp_get_thread_num()); + /* Copy trial into private structure */ - struct kappa_data *min_trial = (struct kappa_data*)malloc(sizeof(struct kappa_data)); + struct kappa_data *min_trial = (struct kappa_data *) malloc(sizeof(struct kappa_data)); kd_init(min_trial); min_trial->parent_subset = ss; -#pragma omp critical + + #pragma omp critical kd_copy_parameters(trial, min_trial); + /* Run local minimization */ minimize_locally(min_trial, 500); calculate_charges(de_ss, min_trial); calculate_statistics(de_ss, min_trial); int cond = 0; + /* If better, swap for so_far_best */ -#pragma omp critical + #pragma omp critical cond = (kd_sort_by_is_better(min_trial, trial) && compare_and_set(min_trial, so_far_best)); if (cond) { calculate_charges(ss, so_far_best); @@ -207,26 +219,26 @@ void generate_random_population(struct subset* ss, float *bounds, int size) { assert(bounds != NULL); /* Get random numbers by Latin Hypercube Sampling */ - int dimensions_count = ts.atom_types_count*2+1; + int dimensions_count = ts.atom_types_count * 2 + 1; int points_count = size; int seed = rand(); double* random_lhs = latin_random_new(dimensions_count, points_count, &seed); /* Redistribute random_lhs[dim_num, point_num] to ss->data */ - //every row contains values of one dimension for all population members - //start with alpha and beta for all atom types + /* Every row contains values of one dimension for all population members */ + /* Start with alpha and beta for all atom types */ for (int i = 0; i < points_count; i++) { for (int j = 0; j < ts.atom_types_count; j++) { - ss->data[i].parameters_alpha[j] = interpolate_to_different_bounds((float)random_lhs[j*points_count + i], bounds[2+j*4], bounds[2+j*4+1]); - ss->data[i].parameters_beta[j] = interpolate_to_different_bounds((float)random_lhs[(j+1)*points_count + i], bounds[2+j*4+2], bounds[2+j*4+3]); + ss->data[i].parameters_alpha[j] = interpolate_to_different_bounds((float)random_lhs[j * points_count + i], bounds[2 + j * 4], bounds[2 + j * 4 + 1]); + ss->data[i].parameters_beta[j] = interpolate_to_different_bounds((float)random_lhs[(j + 1) * points_count + i], bounds[2 + j * 4 + 2], bounds[2 + j * 4 + 3]); } } - //the last row of random_lhs is kappa for all population members + /* The last row of random_lhs is kappa for all population members */ for (int i = 0; i < points_count; i++) { if (s.fixed_kappa < 0) - ss->data[i].kappa = interpolate_to_different_bounds((float)random_lhs[(dimensions_count-2)*points_count + i], bounds[0], bounds[1]); + ss->data[i].kappa = interpolate_to_different_bounds((float)random_lhs[(dimensions_count - 2) * points_count + i], bounds[0], bounds[1]); else ss->data[i].kappa = s.fixed_kappa; } @@ -235,23 +247,24 @@ void generate_random_population(struct subset* ss, float *bounds, int size) { } /* Run local minimization on part of population */ -int minimize_part_of_population(struct subset* ss, int* good_indices) { +int minimize_part_of_population(struct subset * ss, int *good_indices) { assert(ss != NULL); assert(good_indices != NULL); int quite_good = 0; int i = 0; - //we minimize all with R2>0.2 && R>0 -#pragma omp parallel for num_threads(s.de_threads) shared(ss, quite_good, good_indices) private(i) + + /* We minimize all with R2>0.2 && R>0 */ + #pragma omp parallel for num_threads(s.de_threads) shared(ss, quite_good, good_indices) private(i) for (i = 0; i < ss->kappa_data_count; i++) { if (ss->data[i].full_stats.R2 > 0.2 && ss->data[i].full_stats.R > 0) { -#pragma omp critical + #pragma omp critical { good_indices[quite_good] = i; quite_good++; } - struct kappa_data* m = (struct kappa_data*) malloc (sizeof(struct kappa_data)); + struct kappa_data* m = (struct kappa_data *) malloc (sizeof(struct kappa_data)); kd_init(m); m->parent_subset = ss; kd_copy_parameters(&ss->data[i], m); @@ -269,7 +282,7 @@ int minimize_part_of_population(struct subset* ss, int* good_indices) { } /* Evolve kappa_data, i.e. create a new trial structure */ -int evolve_kappa(struct kappa_data* trial, struct kappa_data* x, struct kappa_data* a, struct kappa_data *b, float *bounds, float mutation_constant, float recombination_constant) { +int evolve_kappa(struct kappa_data *trial, struct kappa_data *x, struct kappa_data *a, struct kappa_data *b, float *bounds, float mutation_constant, float recombination_constant) { assert(trial != NULL); assert(x != NULL); @@ -280,36 +293,36 @@ int evolve_kappa(struct kappa_data* trial, struct kappa_data* x, struct kappa_da int changed = 0; kd_copy_parameters(x, trial); - //evolve alpha parameters + /* Evolve alpha parameters */ for (int i = 0; i < ts.atom_types_count; i++) - // if random number is higher than the recombination constant, we will combine i-th atom type - if (get_random_float(0,1) < recombination_constant) { + /* If random number is higher than the recombination constant, we will combine i-th atom type */ + if (get_random_float(0,1) < recombination_constant) { changed++; - trial->parameters_alpha[i] += mutation_constant*(a->parameters_alpha[i] - b->parameters_alpha[i]); - // check bounds, if the evolved parameters are out of bounds, discard changes - if (bounds[2 + i*4] > trial->parameters_alpha[i] || bounds[2 + i*4 + 1] < trial->parameters_alpha[i]) { + trial->parameters_alpha[i] += mutation_constant * (a->parameters_alpha[i] - b->parameters_alpha[i]); + /* Check bounds, if the evolved parameters are out of bounds, discard changes */ + if (bounds[2 + i * 4] > trial->parameters_alpha[i] || bounds[2 + i * 4 + 1] < trial->parameters_alpha[i]) { trial->parameters_alpha[i] = x->parameters_alpha[i]; changed--; } } - //evolve beta parameters + /* Evolve beta parameters */ for (int i = 0; i < ts.atom_types_count; i++) - if (get_random_float(0,1) < recombination_constant) { + if (get_random_float(0,1) < recombination_constant) { changed++; trial->parameters_beta[i] += mutation_constant*(a->parameters_beta[i] - b->parameters_beta[i]); - if (bounds[2 + i*4 + 2] > trial->parameters_beta[i] || bounds[2 + i*4 + 3] < trial->parameters_beta[i]) { + if (bounds[2 + i * 4 + 2] > trial->parameters_beta[i] || bounds[2 + i * 4 + 3] < trial->parameters_beta[i]) { trial->parameters_beta[i] = x->parameters_beta[i]; changed--; } } - //always change kappa - if (s.fixed_kappa < 0) { - trial->kappa += mutation_constant*(a->kappa - b->kappa)/(bounds[1]-bounds[0]); + /* Always change kappa */ + if (s.fixed_kappa < 0) { + trial->kappa += mutation_constant * (a->kappa - b->kappa) / (bounds[1] - bounds[0]); changed++; if (bounds[0] > trial->kappa || bounds[1] < trial->kappa) { - trial->kappa -= mutation_constant*(a->kappa - b->kappa)/(bounds[1]-bounds[0]); + trial->kappa -= mutation_constant * (a->kappa - b->kappa) / (bounds[1] - bounds[0]); } } else @@ -336,29 +349,32 @@ void minimize_locally(struct kappa_data* t, int max_calls) { assert(t != NULL); - int n = 2*ts.atom_types_count + 1; //number of variables - int npt = 2*n + 1; //number of interpolation conditions - double* x = (double*) malloc(n*sizeof(double)); + int n = 2 * ts.atom_types_count + 1; /* Number of variables */ + int npt = 2 * n + 1; /* Number of interpolation conditions */ + double* x = (double *) malloc(n * sizeof(double)); kappa_data_to_double_array(t, x); + double rhobeg = 0.2; double rhoend = 0.0001; int iprint = 0; int maxfun = max_calls; - double* w = (double*) malloc(((npt+13)*(npt+n) + 3*n*(n+3)/2)*sizeof(double)); - //call fortran code NEWUOA for local minimization + double *w = (double *) malloc(((npt + 13) * (npt + n) + 3 * n * (n + 3) / 2) * sizeof(double)); + + /* Call fortran code NEWUOA for local minimization */ newuoa_(&n, &npt, x, &rhobeg, &rhoend, &iprint, &maxfun, w); double_array_to_kappa_data(x, t); + free(w); free(x); } /* Used by NEWUOA algorithm. Evaluates the vector in the local minimization: converts it to kappa_data, computes charges, computes statistics and return the fitness score that should be minimized */ -extern void calfun_(int n, double*x, double* f) { +extern void calfun_(int n, double *x, double *f) { assert(x != NULL); assert(f != NULL); - struct kappa_data* t = (struct kappa_data*) malloc (sizeof(struct kappa_data)); + struct kappa_data *t = (struct kappa_data *) malloc (sizeof(struct kappa_data)); kd_init(t); double_array_to_kappa_data(x, t); calculate_charges(de_ss, t); @@ -369,7 +385,7 @@ extern void calfun_(int n, double*x, double* f) { case SORT_R2: case SORT_RW: case SORT_SPEARMAN: - *f = 1 - (double)(result) + n - n; //+n-n just to get rid of compilaiton warning + *f = 1 - (double)(result) + n - n; /* +n-n just to get rid of compilaiton warning */ break; default: *f = (double) (result) + n -n; @@ -379,33 +395,33 @@ extern void calfun_(int n, double*x, double* f) { } /* Convert kappa_data into an array of doubles, used in local minimization */ -void kappa_data_to_double_array(struct kappa_data* t, double* x) { +void kappa_data_to_double_array(struct kappa_data *t, double *x) { assert(t != NULL); assert(x != NULL); x[0] = t->kappa; for (int i = 0; i < ts.atom_types_count; i++) { - x[i*2 + 1] = t->parameters_alpha[i]; - x[i*2 + 2] = t->parameters_beta[i]; + x[i * 2 + 1] = t->parameters_alpha[i]; + x[i * 2 + 2] = t->parameters_beta[i]; } } /* Convert array of doubles into kappa_data, used in local minimization */ -void double_array_to_kappa_data(double* x, struct kappa_data* t) { +void double_array_to_kappa_data(double *x, struct kappa_data *t) { assert(x != NULL); assert(t != NULL); t->kappa = (float)x[0]; for (int i = 0; i < ts.atom_types_count; i++) { - t->parameters_alpha[i] = (float)x[i*2 + 1]; - t->parameters_beta[i] = (float)x[i*2 + 2]; + t->parameters_alpha[i] = (float) x[i * 2 + 1]; + t->parameters_beta[i] = (float) x[i * 2 + 2]; } } /* Returns true if R2 is above 0.6, used in decision whether to minimize kappa_data */ -int is_quite_good(struct kappa_data* t) { +int is_quite_good(const struct kappa_data * const t) { assert(t != NULL); @@ -419,32 +435,28 @@ void compute_parameters_bounds(float* bounds, int by_atom_type) { assert(bounds != NULL); - //returns bounds[kappa_low, kappa_high, alpha_1_low, alpha_1_high, beta_1_low, beta_1_high, alpha_2_low, ...] + /* Returns bounds[kappa_low, kappa_high, alpha_1_low, alpha_1_high, beta_1_low, beta_1_high, alpha_2_low, ...] */ float toH = 0.036749309; - bounds[0] = 0.0005; //kappa_low - bounds[1] = 3.5; //kappa_high + bounds[0] = 0.0005; /* kappa_low */ + bounds[1] = 3.5; /* kappa_high */ for (int j = 0; j < ts.atom_types_count; j++) { - //set bounds to values corresponding with their affinity and ionenergies + /* Set bounds to values corresponding with their affinity and ionenergies */ if (by_atom_type == 1) { - //set bounds for particular atom type + /* Set bounds for particular atom type */ int atom_number = ts.atom_types[j].Z; - bounds[2 + j*4] = (((float)(ionenergies[atom_number]) + (float)(affinities[atom_number])*5)/2)*toH; //alpha_low - bounds[2 + j*4 + 1] = (float) (bounds[2 + j*4] + 0.1); //alpha_high - bounds[2 + j*4] -= 0.1; - //bounds[2 + j*4] *= toEV; - //bounds[2 + j*4 + 1] *= toEV; - bounds[2 + j*4 + 2] = ((float)(ionenergies[atom_number]) - (float)(affinities[atom_number])*5)/2*toH; //beta_low - bounds[2 + j*4 + 3] = (float) (bounds[2 + j*4 + 2] + 0.1); //beta_high - bounds[2 + j*4 + 2] -= 0.1; - //bounds[2 + j*4 + 2] *= toEV; - //bounds[2 + j*4 + 3] *= toEV; + bounds[2 + j * 4] = (((float)(ionenergies[atom_number]) + (float)(affinities[atom_number]) * 5) / 2) * toH; /* alpha_low */ + bounds[2 + j * 4 + 1] = (float) (bounds[2 + j * 4] + 0.1); /* alpha_high */ + bounds[2 + j * 4] -= 0.1; + bounds[2 + j * 4 + 2] = ((float)(ionenergies[atom_number]) - (float)(affinities[atom_number]) * 5) / 2* toH; /* beta_low */ + bounds[2 + j * 4 + 3] = (float) (bounds[2 + j * 4 + 2] + 0.1); /* beta_high */ + bounds[2 + j * 4 + 2] -= 0.1; } - //set bounds to constant values taken from previous experience + /* Set bounds to constant values taken from previous experience */ else if (by_atom_type == 0) { - bounds[2 + j*4] = 1.8; - bounds[2 + j*4 + 1] = 3.2; - bounds[2 + j*4 + 2] = 0; - bounds[2 + j*4 + 3] = 1.0; + bounds[2 + j * 4] = 1.8; + bounds[2 + j * 4 + 1] = 3.2; + bounds[2 + j * 4 + 2] = 0; + bounds[2 + j * 4 + 3] = 1.0; } } if (s.verbosity >= VERBOSE_KAPPA) { @@ -452,7 +464,8 @@ void compute_parameters_bounds(float* bounds, int by_atom_type) { for (int i = 0; i < ts.atom_types_count; i++) { char buff[10]; at_format_text(&ts.atom_types[i], buff); - printf("%s %5.3f - %5.3f, %5.3f - %5.3f\n", buff, bounds[2+i*4], bounds[2+i*4+1], bounds[2+i*4+2], bounds[2+i*4+3]); + printf("%s %5.3f - %5.3f, %5.3f - %5.3f\n", buff,\ + bounds[2 + i * 4], bounds[2 + i * 4 + 1], bounds[2 + i * 4 + 2], bounds[2 + i * 4 + 3]); } } } @@ -460,15 +473,15 @@ void compute_parameters_bounds(float* bounds, int by_atom_type) { /* Get random float within the bounds */ float get_random_float(float low, float high) { - //better random number generator would be nice, but this is sufficient - float n = low + (float)(rand())/((float)(RAND_MAX/(high-low))); + /* Better random number generator would be nice, but this is sufficient */ + float n = low + (float) (rand()) / ((float) (RAND_MAX / (high - low))); return n; } /* Interpolate the number from [0,1] to [low, high] */ float interpolate_to_different_bounds(float x, float low, float high) { - return low + x*(high-low); + return low + x * (high - low); } /* Compute the sum of the vector */ @@ -476,6 +489,6 @@ int sum(int* vector, int size) { int sum = 0; for (int i = 0; i < size; i++) - sum+=vector[i]; + sum += vector[i]; return sum; } diff --git a/src/diffevolution.h b/src/diffevolution.h index 971ca90..2bbf2d7 100644 --- a/src/diffevolution.h +++ b/src/diffevolution.h @@ -25,7 +25,7 @@ extern void newuoa_(int* n, int* npt, double* x, double* rhobeg, double* rhoend, extern void calfun_(int n, double*x, double* f); void kappa_data_to_double_array(struct kappa_data* trial, double* x); void double_array_to_kappa_data(double* x, struct kappa_data* trial); -int is_quite_good(struct kappa_data* t); +int is_quite_good(const struct kappa_data * const t); struct subset * de_ss; #endif /* __DIFFEVOLUTION_H__ */ diff --git a/src/guidedmin.c b/src/guidedmin.c index e5f6284..1d4dcaa 100644 --- a/src/guidedmin.c +++ b/src/guidedmin.c @@ -38,41 +38,45 @@ void run_guided_min(struct subset * const ss) { printf("GM Generating %d vectors\n", s.gm_size); /* Set bounds for each parameters in kappa_data */ - float *bounds = (float*) malloc((ts.atom_types_count*2+1)*2*sizeof(float)); - //compute bounds, 0 means set them to fixed numbers taken from Tomas's full scan, 1 means try to find them with broad search + float *bounds = (float *) malloc((ts.atom_types_count * 2 + 1) * 2 * sizeof (float)); + + /* Compute bounds, 0 means set them to fixed numbers taken from Tomas's full scan, 1 means try to find them with broad search */ compute_parameters_bounds(bounds, 0); - //generate population + + /* Generate population */ fill_ss(ss, s.gm_size); generate_random_population(ss, bounds, s.gm_size); - de_ss = ss; - //evaluate the fitness function for all points + de_ss = ss; + + /* Evaluate the fitness function for all points */ if (s.verbosity >= VERBOSE_KAPPA) printf("GM Calculating charges and evaluating fitness function for whole set\n"); + int i = 0; -#pragma omp parallel for num_threads(s.gm_threads) default(shared) private(i) + #pragma omp parallel for num_threads(s.gm_threads) default(shared) private(i) for (i = 0; i < ss->kappa_data_count; i++) { calculate_charges(ss, &ss->data[i]); calculate_statistics(ss, &ss->data[i]); } - //minimize part of population that has R>0.3 + /* Minimize part of population that has R > 0.3 */ int minimized_initial = 0; if (s.verbosity >= VERBOSE_KAPPA) printf("GM minimizing population\n"); minimized_initial = minimize_part_of_gm_set(ss, s.gm_iterations_beg); - //if we minimized zero data, inform user and suggest larger set + /* If we minimized zero data, inform user and suggest larger set */ if (minimized_initial == 0) { EXIT_ERROR(RUN_ERROR, "No vector in set was worth minimizing. Please choose larger set than %d (option --gm-size)\n.", s.gm_size); } - //find the best kappa_data + /* Find the best kappa_data */ set_the_best(ss); - struct kappa_data *so_far_best = (struct kappa_data*)malloc(sizeof(struct kappa_data)); + struct kappa_data *so_far_best = (struct kappa_data *) malloc(sizeof(struct kappa_data)); kd_init(so_far_best); - //copy ss->best into so_far_best + /* Copy ss->best into so_far_best */ kd_copy_parameters(ss->best, so_far_best); calculate_charges(ss, so_far_best); calculate_statistics(ss, so_far_best); @@ -96,15 +100,15 @@ int minimize_part_of_gm_set(struct subset* ss, int min_iterations) { int quite_good = 0; int i = 0; - //we minimize all with R2>0.2 && R>0 -#pragma omp parallel for num_threads(s.gm_threads) shared(ss, quite_good) private(i) + /* We minimize all with R2 > 0.2 && R > 0 */ + #pragma omp parallel for num_threads(s.gm_threads) shared(ss, quite_good) private(i) for (i = 0; i < ss->kappa_data_count; i++) { if (ss->data[i].full_stats.R2 > 0.2 && ss->data[i].full_stats.R > 0) { -#pragma omp critical + #pragma omp critical { quite_good++; } - struct kappa_data* m = (struct kappa_data*) malloc (sizeof(struct kappa_data)); + struct kappa_data* m = (struct kappa_data *) malloc (sizeof(struct kappa_data)); kd_init(m); m->parent_subset = ss; kd_copy_parameters(&ss->data[i], m); diff --git a/src/guidedmin.h b/src/guidedmin.h index 211859a..b456a2b 100644 --- a/src/guidedmin.h +++ b/src/guidedmin.h @@ -13,17 +13,7 @@ #include "diffevolution.h" void run_guided_min(struct subset * const ss); -//void generate_random_population(struct subset* ss, float *bounds, int size); -int minimize_part_of_gm_set(struct subset* ss, int min_iterations); -//void compute_parameters_bounds(float* bounds, int by_atom_type); -//float get_random_float(float low, float high); -//float interpolate_to_different_bounds(float x, float low, float high); -//int sum(int* vector, int size); -//void minimize_locally(struct kappa_data* trial, int max_calls); -//extern void newuoa_(int* n, int* npt, double* x, double* rhobeg, double* rhoend, int* iprint, int* maxfun, double* w); -//extern void calfun_(int n, double*x, double* f); -//void kappa_data_to_double_array(struct kappa_data* trial, double* x); -//void double_array_to_kappa_data(double* x, struct kappa_data* trial); +int minimize_part_of_gm_set(struct subset *ss, int min_iterations); -struct subset * de_ss; +struct subset *de_ss; #endif /* __GUIDEDMIN_H__ */ diff --git a/src/kappa.c b/src/kappa.c index 5b97108..92c8c3d 100644 --- a/src/kappa.c +++ b/src/kappa.c @@ -227,11 +227,11 @@ void find_the_best_parameters_for_subset(struct subset * const ss) { } if (s.params_method == PARAMS_DE) { - //runs a differential evolution algorithm to find the best parameters, ss->best is set after the call + /* Runs a differential evolution algorithm to find the best parameters, ss->best is set after the call */ run_diff_evolution(ss); } if (s.params_method == PARAMS_GM) { - //runs a guided minimization algorithm, ss->best is set after the call + /* Runs a guided minimization algorithm, ss->best is set after the call */ run_guided_min(ss); } @@ -245,7 +245,7 @@ void find_the_best_parameters_for_subset(struct subset * const ss) { ss->best = &ss->data[ss->kappa_data_count - 1]; } else if (s.params_method == PARAMS_DE || s.params_method == PARAMS_GM) { - //well, nothing, the best structure has been already set + /* well, nothing, the best structure has been already set */ } } diff --git a/src/settings.c b/src/settings.c index 6a84231..c41b3af 100644 --- a/src/settings.c +++ b/src/settings.c @@ -98,7 +98,7 @@ void s_init(void) { s.de_threads = 1; s.limit_de_iters = NO_LIMIT_ITERS; s.limit_de_time = NO_LIMIT_TIME; - s.polish = -1; //0 off, 1 only result, 2 result + during evolve, 3 result, evolve and some structures in initial population + s.polish = -1; /* 0 off, 1 only result, 2 result + during evolve, 3 result, evolve and some structures in initial population */ s.gm_size = 100; s.gm_iterations_beg = 1000; s.gm_iterations_end = 2000; @@ -367,7 +367,7 @@ void parse_options(int argc, char **argv) { case 173: s.rw = (float) atof(optarg); break; - //DE settings + /* DE settings */ case 180: s.population_size = atoi(optarg); break; @@ -411,7 +411,7 @@ void parse_options(int argc, char **argv) { case 190: s.polish = atoi(optarg); break; - //GM settings + /* GM settings */ case 191: s.gm_size = atoi(optarg); break; @@ -454,7 +454,7 @@ void check_settings(void) { if(s.mode == MODE_PARAMS) { if(s.chg_file[0] == '\0') EXIT_ERROR(ARG_ERROR, "%s", "No .chg file provided.\n"); - //if user did not specify the optimization method for parameters calculation, set linear regression + /* If user did not specify the optimization method for parameters calculation, set linear regression */ if (s.params_method == PARAMS_NOT_SET) s.params_method = PARAMS_LR_FULL; @@ -478,12 +478,13 @@ void check_settings(void) { } } - if (s.params_method == PARAMS_DE) { //all settings are optional, so check for mistakes and set defaults + if (s.params_method == PARAMS_DE) { + /* All settings are optional, so check for mistakes and set defaults */ if (s.population_size == 0) - s.population_size = 500; //1.2*(ts.atom_types_count*2+1); + s.population_size = 500; /* 1.2 * (ts.atom_types_count * 2 + 1); */ if (s.limit_de_iters == NO_LIMIT_ITERS && s.limit_de_time == NO_LIMIT_TIME) s.limit_de_iters = 1000; - if (s.mutation_constant < 0) // if not set + if (s.mutation_constant < 0) /* If not set */ s.mutation_constant = 0.75; if (s.recombination_constant < 0) s.recombination_constant = 0.7; @@ -492,7 +493,8 @@ void check_settings(void) { } - if (s.params_method == PARAMS_GM) { //all settings are optional, so check for mistakes + if (s.params_method == PARAMS_GM) { + /* All settings are optional, so check for mistakes */ if (s.gm_size < 1) EXIT_ERROR(ARG_ERROR, "%s", "Size of GM set has to be positive.\n"); if (s.gm_iterations_beg < 1 || s.gm_iterations_end < 1) @@ -512,7 +514,8 @@ void check_settings(void) { if(s.limit_time != NO_LIMIT_TIME && s.limit_time > 36000 * 1000) EXIT_ERROR(ARG_ERROR, "%s", "Maximum time should not be higher than 1000 hours.\n"); - //TODO verify with Tomas if this is the intended behavior + + /* TODO verify with Tomas if this is the intended behavior */ if(s.params_method == PARAMS_LR_FULL_BRENT /*!s.full_scan_only*/ && (s.sort_by != SORT_R && s.sort_by != SORT_R2 && s.sort_by != SORT_SPEARMAN && s.sort_by != SORT_RW)) EXIT_ERROR(ARG_ERROR, "%s", "Full scan must be used for sort-by other than R, R2 or Spearman.\n"); @@ -728,7 +731,6 @@ void print_settings(void) { printf("\t - iterations for set at the beginning %d\n", s.gm_iterations_beg); printf("\t - iterations for the result at the end %d\n", s.gm_iterations_end); printf("\t - threads used for minimization %d\n", s.gm_threads); - } } diff --git a/src/settings.h b/src/settings.h index 2930d99..bb735da 100644 --- a/src/settings.h +++ b/src/settings.h @@ -27,7 +27,7 @@ enum params_calc_method { PARAMS_LR_FULL, PARAMS_LR_FULL_BRENT, PARAMS_DE, - PARAMS_GM, //guided minimization + PARAMS_GM, PARAMS_NOT_SET }; @@ -80,29 +80,29 @@ struct settings { enum atom_type_customization at_customization; enum discarding_mode discard; - //settings regarding PARAMS_LR_FULL* parameters' calculation method + /* Settings regarding PARAMS_LR_FULL* parameters' calculation method */ float kappa_max; float kappa_set; float full_scan_precision; - //settings regarding PARAMS_DE optimization method + /* Settings regarding PARAMS_DE optimization method */ int population_size; float mutation_constant; - int dither; //set mutation constant to random value from [0.5, 1] each iteration + int dither; /* Set mutation constant to random value from [0.5, 1] each iteration */ float recombination_constant; - float fixed_kappa; //fix kappa to given value - int de_threads; //number of threads used to paralellize DE + float fixed_kappa; /* Fix kappa to given value */ + int de_threads; /* Number of threads used to paralellize DE */ int limit_de_iters; time_t limit_de_time; - int polish; //use NEWUOA minimization to polish trial or results + int polish; /* Use NEWUOA minimization to polish trial or results */ - //settings regarding PARAMS_GM optimization method + /* Settings regarding PARAMS_GM optimization method */ int gm_size; int gm_iterations_beg; int gm_iterations_end; int gm_threads; - //other settings + /* Other settings */ int random_seed; enum verbosity_levels verbosity; diff --git a/src/statistics.c b/src/statistics.c index 497dcee..decd87c 100644 --- a/src/statistics.c +++ b/src/statistics.c @@ -92,8 +92,8 @@ static void set_total_R_w(struct kappa_data * const kd) { weighted_corr_sum += weight * kd->per_at_stats[i].R; } - //consider total R - weighted_corr_sum += 3 * (pow(s.rw, kd->full_stats.R2)) * kd->full_stats.R2 - kd->full_stats.RMSD/3; + /* consider total R */ + weighted_corr_sum += 3 * (pow(s.rw, kd->full_stats.R2)) * kd->full_stats.R2 - kd->full_stats.RMSD / 3; /* Normalize the results */ kd->full_stats.R_w = (float) (weighted_corr_sum / (ts.atom_types_count * s.rw + 3 * s.rw));