Ian Jauslin
summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorIan Jauslin <ian@jauslin.org>2022-06-14 09:26:07 +0200
committerIan Jauslin <ian@jauslin.org>2022-06-14 09:46:36 +0200
commit3f0510629e422e979b57d3f93791937912a4183a (patch)
treebf2589b2689044261b0cd4d9e6b3082194fdd9e9 /src/mean.c
parent469bdc80712dbf9c12562059dc4594620b59a076 (diff)
Update to v1.5.HEADmaster
The update to version 1.5 is rather substantial, and introduces some minor backward-incompatibilities: * The header "#!symbols" has been replaced by "#!virtual_fields" * Multiplying polynomials using the '*' symbol is no longer supported (or, rather, the symbolic capabilities of meankondo were enhanced, and the syntax has been changed). * 'meantools exp' has been removed (its functionality is now handled by other means) * 'meantoolds derive' has been replaced by 'meantools differentiate' * The symbolic capabilities were enhanced: polynomials can now be multiplied, added, exponentiated, and their logarithms can be taken directly in the configuration file. * The flow equation can now be processed after being computed using the various "#!postprocess_*" entries. * Deprecated kondo_preprocess. * Compute the mean using an LU decomposition if possible. * More detailed checks for syntax errors in configuration file. * Check that different '#!group' entries are indeed uncorrelated. * New flags in meankondo: '-p' and '-A'. * New tool: meantools expand. * Improve conversion to LaTeX using meantools-convert * Assign terms randomly to different threads. * Created vim files to implement syntax highlighting for configuration files. * Multiple bug fixes
Diffstat (limited to 'src/mean.c')
-rw-r--r--src/mean.c141
1 files changed, 109 insertions, 32 deletions
diff --git a/src/mean.c b/src/mean.c
index 39ece56..0c5707e 100644
--- a/src/mean.c
+++ b/src/mean.c
@@ -1,5 +1,5 @@
/*
-Copyright 2015 Ian Jauslin
+Copyright 2015-2022 Ian Jauslin
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
@@ -30,6 +30,7 @@ As of version 1.0, the mean of a monomial is computed directly
#include "array.h"
#include "fields.h"
#include "number.h"
+#include "determinant.h"
// mean of a monomial
int mean(Int_Array monomial, Polynomial* out, Fields_Table fields, Polynomial_Matrix propagator){
@@ -61,10 +62,79 @@ int mean(Int_Array monomial, Polynomial* out, Fields_Table fields, Polynomial_Ma
// compute the mean of a monomial of internal fields (with split + and -)
int mean_internal(Int_Array internal_plus, Int_Array internal_minus, Polynomial* out, Polynomial_Matrix propagator, Fields_Table fields){
+ int ret;
+ Number num;
+
if(internal_plus.length!=internal_minus.length){
fprintf(stderr,"error: monomial contains unmatched +/- fields\n");
exit(-1);
}
+
+ ret=mean_determinant(internal_plus, internal_minus, &num, propagator, fields);
+ // cannot compute the mean as a determinant, use permutations
+ // can be because some fields are not Fermions
+ // can be because the propagator has non-numeric values (inverting polynomials is not implemented, and would be required for the computation of the determinant)
+ if(ret==-1){
+ mean_permutations(internal_plus, internal_minus, out, propagator, fields);
+ }
+ else{
+ polynomial_multiply_scalar(*out, num);
+ free_Number(num);
+ }
+
+ return(0);
+}
+
+// compute the mean of a monomial by computing a determinant
+// can only be used if all of the propagators are numbers
+int mean_determinant(Int_Array internal_plus, Int_Array internal_minus, Number* out, Polynomial_Matrix propagator, Fields_Table fields){
+ Number_Matrix M;
+ int n=internal_minus.length;
+ int i,j;
+ int a,b;
+ int sign;
+
+ init_Number_Matrix(&M,n);
+
+ // extra sign: the monomial is sorted in such a way that minus fields are on the left of plus fields, but the determinant formula requires the fields to be alternated +-
+ if((n+1)/2%2==1){
+ sign=-1;
+ }
+ else{
+ sign=1;
+ }
+
+ // construct matrix
+ for(i=0;i<n;i++){
+ a=intlist_find_err(propagator.indices, propagator.length, internal_plus.values[i]);
+ for(j=0;j<n;j++){
+ b=intlist_find_err(propagator.indices, propagator.length, -internal_minus.values[j]);
+ // ignore 0
+ if(propagator.matrix[a][b].length!=0){
+ // check whether the fields are Fermions, and whether the entry is a number
+ if(is_fermion(internal_plus.values[i], fields)==0 || is_fermion(internal_minus.values[j], fields)==0 || polynomial_is_number(propagator.matrix[a][b])==0){
+ free_Number_Matrix(M);
+ return(-1);
+ }
+
+ number_add_chain(propagator.matrix[a][b].nums[0], M.matrix[i]+j);
+ }
+ }
+ }
+
+ // compute determinant
+ determinant_inplace(M, out);
+
+ number_Qprod_chain(quot(sign,1), out);
+
+ free_Number_Matrix(M);
+
+ return(0);
+}
+
+
+// compute the mean of a monomial by summing over permutations
+int mean_permutations(Int_Array internal_plus, Int_Array internal_minus, Polynomial* out, Polynomial_Matrix propagator, Fields_Table fields){
int n=internal_minus.length;
// pairing as an array of positions
int* pairing=calloc(n,sizeof(int));
@@ -118,7 +188,7 @@ int mean_internal(Int_Array internal_plus, Int_Array internal_minus, Polynomial*
pairing_sign=permutation_signature(pairing,n);
}
- // only simplify in mean_symbols
+ // only simplify in mean_virtual_fields
polynomial_prod_chain_nosimplify(num_summed,out,fields);
free_Polynomial(num_summed);
free(pairing);
@@ -346,10 +416,10 @@ int get_internals(Int_Array monomial, Int_Array* internal_plus, Int_Array* inter
}
-// compute the mean of a monomial containing symbolic expressions
+// compute the mean of a monomial containing virtual fields
// keep track of which means were already computed
-int mean_symbols(Int_Array monomial, Polynomial* output, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, Identities* computed){
- Int_Array symbol_list;
+int mean_virtual_fields(Int_Array monomial, Polynomial* output, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, Identities* computed){
+ Int_Array virtual_field_list;
int i;
int power;
int* current_term;
@@ -375,43 +445,43 @@ int mean_symbols(Int_Array monomial, Polynomial* output, Fields_Table fields, Po
}
}
- init_Int_Array(&symbol_list, monomial.length);
+ init_Int_Array(&virtual_field_list, monomial.length);
init_Int_Array(&base_monomial, monomial.length);
- // generate symbols list
+ // generate virtual_fields list
for(i=0;i<monomial.length;i++){
- if(field_type(monomial.values[i], fields)==FIELD_SYMBOL){
- int_array_append(intlist_find_err(fields.symbols.indices, fields.symbols.length, monomial.values[i]), &symbol_list);
+ if(field_type(monomial.values[i], fields)==FIELD_VIRTUAL){
+ int_array_append(intlist_find_err(fields.virtual_fields.indices, fields.virtual_fields.length, monomial.values[i]), &virtual_field_list);
}
else{
int_array_append(monomial.values[i], &base_monomial);
}
}
- power=symbol_list.length;
+ power=virtual_field_list.length;
// trivial case
if(power==0){
mean(monomial, &mean_num, fields, propagator);
polynomial_concat_noinit(mean_num, output);
- free_Int_Array(symbol_list);
+ free_Int_Array(virtual_field_list);
free_Int_Array(base_monomial);
return(0);
}
else{
// initialize current term to a position that has no repetitions
current_term=calloc(power,sizeof(int));
- exists_next=init_prod(current_term, symbol_list, fields, power, base_monomial)+1;
+ exists_next=init_prod(current_term, virtual_field_list, fields, power, base_monomial)+1;
}
- // loop over terms; the loop stops when all the pointers are at the end of the first symbol
+ // loop over terms; the loop stops when all the pointers are at the end of the first virtual field
while(exists_next==1){
// construct monomial
int_array_cpy(base_monomial, &tmp_monomial);
tmp_num=number_one();
for(i=0;i<power;i++){
- int_array_concat(fields.symbols.expr[symbol_list.values[i]].monomials[current_term[i]], &tmp_monomial);
- number_prod_chain(fields.symbols.expr[symbol_list.values[i]].nums[current_term[i]], &tmp_num);
+ int_array_concat(fields.virtual_fields.expr[virtual_field_list.values[i]].monomials[current_term[i]], &tmp_monomial);
+ number_prod_chain(fields.virtual_fields.expr[virtual_field_list.values[i]].nums[current_term[i]], &tmp_num);
}
// check whether the monomial vanishes
if(check_monomial_match(tmp_monomial, fields)==1){
@@ -432,7 +502,7 @@ int mean_symbols(Int_Array monomial, Polynomial* output, Fields_Table fields, Po
free_Int_Array(tmp_monomial);
// next term
- exists_next=next_prod(current_term, symbol_list, fields, power, base_monomial)+1;
+ exists_next=next_prod(current_term, virtual_field_list, fields, power, base_monomial)+1;
// simplfiy every 25 steps (improves both memory usage and performance)
@@ -451,13 +521,13 @@ int mean_symbols(Int_Array monomial, Polynomial* output, Fields_Table fields, Po
// free memory
free(current_term);
- free_Int_Array(symbol_list);
+ free_Int_Array(virtual_field_list);
free_Int_Array(base_monomial);
return(0);
}
// first term in product with no repetitions
-int init_prod(int* current_term, Int_Array symbol_list, Fields_Table fields, int power, Int_Array base_monomial){
+int init_prod(int* current_term, Int_Array virtual_field_list, Fields_Table fields, int power, Int_Array base_monomial){
// index we want to increment
int move=0;
// tmp monomial
@@ -474,7 +544,7 @@ int init_prod(int* current_term, Int_Array symbol_list, Fields_Table fields, int
// loop until move is out of range
while(move>=0 && move<power){
// move
- current_term[move]=next_term_norepeat(current_term[move], fields.symbols.expr[symbol_list.values[move]], &monomial, fields);
+ current_term[move]=next_term_norepeat(current_term[move], fields.virtual_fields.expr[virtual_field_list.values[move]], &monomial, fields);
// if the next term does not exist, then move previous index
if(current_term[move]==-1){
move--;
@@ -495,7 +565,7 @@ int init_prod(int* current_term, Int_Array symbol_list, Fields_Table fields, int
}
// next term in product with no repetitions
-int next_prod(int* current_term, Int_Array symbol_list, Fields_Table fields, int power, Int_Array base_monomial){
+int next_prod(int* current_term, Int_Array virtual_field_list, Fields_Table fields, int power, Int_Array base_monomial){
// index we want to increment
int move=power-1;
// tmp monomial
@@ -506,13 +576,13 @@ int next_prod(int* current_term, Int_Array symbol_list, Fields_Table fields, int
init_Int_Array(&monomial, base_monomial.length+5*power);
int_array_cpy_noinit(base_monomial, &monomial);
for(i=0;i<=move;i++){
- int_array_concat(fields.symbols.expr[symbol_list.values[i]].monomials[current_term[i]],&monomial);
+ int_array_concat(fields.virtual_fields.expr[virtual_field_list.values[i]].monomials[current_term[i]],&monomial);
}
// loop until move is out of range
while(move>=0 && move<power){
// move
- current_term[move]=next_term_norepeat(current_term[move], fields.symbols.expr[symbol_list.values[move]], &monomial, fields);
+ current_term[move]=next_term_norepeat(current_term[move], fields.virtual_fields.expr[virtual_field_list.values[move]], &monomial, fields);
// if the next term does not exist, then move previous index
if(current_term[move]==-1){
move--;
@@ -623,13 +693,13 @@ int mean_groups(Int_Array monomial, Polynomial* output, Fields_Table fields, Pol
int group=-2;
int next_group=-2;
Polynomial tmp_poly;
- int sign;
+ int sign=1;
init_Polynomial(output, MONOMIAL_SIZE);
- // check whether there are symbols
- // IMPORTANT: the symbols must be at the end of the monomial
- if(monomial.length==0 || field_type(monomial.values[monomial.length-1], fields)!=FIELD_SYMBOL){
+ // check whether there are virtual fields
+ // IMPORTANT: the virtual fields must be at the end of the monomial
+ if(monomial.length==0 || field_type(monomial.values[monomial.length-1], fields)!=FIELD_VIRTUAL){
// mean
mean(monomial, &num_mean, fields, propagator);
// add to output
@@ -650,7 +720,7 @@ int mean_groups(Int_Array monomial, Polynomial* output, Fields_Table fields, Pol
}
// if group changes, take mean
if((i>0 && next_group!=group) || i==monomial.length){
- mean_symbols(tmp_monomial, &tmp_poly, fields, propagator, groups, computed);
+ mean_virtual_fields(tmp_monomial, &tmp_poly, fields, propagator, groups, computed);
// if zero
if(polynomial_is_zero(tmp_poly)==1){
// set output to 0
@@ -702,10 +772,11 @@ struct mean_args{
Fields_Table fields;
Polynomial_Matrix propagator;
Groups groups;
+ int print_progress;
};
// multithreaded
-int polynomial_mean_multithread(Polynomial* polynomial, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, int threads){
+int polynomial_mean_multithread(Polynomial* polynomial, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, int threads, int print_progress){
int i;
Polynomial thread_polys[threads];
pthread_t thread_ids[threads];
@@ -720,11 +791,15 @@ int polynomial_mean_multithread(Polynomial* polynomial, Fields_Table fields, Pol
args[i].fields=fields;
args[i].propagator=propagator;
args[i].groups=groups;
+ args[i].print_progress=print_progress;
}
// split polynomial
+ // randomly choose the thread
+ // see random number generator
+ srand(time(NULL));
for(i=0;i<len;i++){
- polynomial_append((*polynomial).monomials[i], (*polynomial).factors[i], (*polynomial).nums[i], thread_polys+(i % threads));
+ polynomial_append((*polynomial).monomials[i], (*polynomial).factors[i], (*polynomial).nums[i], thread_polys+(rand() % threads));
}
// start threads
@@ -750,12 +825,12 @@ int polynomial_mean_multithread(Polynomial* polynomial, Fields_Table fields, Pol
// mean for one of the threads
void* polynomial_mean_thread(void* mean_args){
struct mean_args *args=mean_args;
- polynomial_mean((*args).polynomial,(*args).fields,(*args).propagator,(*args).groups);
+ polynomial_mean((*args).polynomial,(*args).fields,(*args).propagator,(*args).groups, (*args).print_progress);
return(NULL);
}
// single threaded version
-int polynomial_mean(Polynomial* polynomial, Fields_Table fields, Polynomial_Matrix propagator, Groups groups){
+int polynomial_mean(Polynomial* polynomial, Fields_Table fields, Polynomial_Matrix propagator, Groups groups, int print_progress){
int i,j;
Polynomial output;
Polynomial tmp_poly;
@@ -769,7 +844,9 @@ int polynomial_mean(Polynomial* polynomial, Fields_Table fields, Polynomial_Matr
// mean of each monomial
for(i=0;i<(*polynomial).length;i++){
- fprintf(stderr,"computing %d of %d means\n",i,(*polynomial).length-1);
+ if(print_progress==1){
+ fprintf(stderr,"computing %d of %d means\n",i,(*polynomial).length-1);
+ }
mean_groups((*polynomial).monomials[i], &tmp_poly, fields, propagator, groups, &computed);
// write factors