From 167980ea437881ec56186332370afcc169f2e4dd Mon Sep 17 00:00:00 2001
From: Ian Jauslin <ian@jauslin.org>
Date: Tue, 14 Jun 2022 09:26:07 +0200
Subject: Update to v1.5.

  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.

  * Multiple bug fixes
---
 src/mean.c | 141 +++++++++++++++++++++++++++++++++++++++++++++++--------------
 1 file changed, 109 insertions(+), 32 deletions(-)

(limited to 'src/mean.c')

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
-- 
cgit v1.2.3-70-g09d2