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

1 files changed, 114 insertions, 24 deletions

diff --git a/src/grouped_polynomial.c b/src/grouped_polynomial.c
index 0d4d75d..9e245a8 100644
--- a/src/grouped_polynomial.c
+++ b/src/grouped_polynomial.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.
@@ -205,9 +205,9 @@ int group_polynomial(Polynomial polynomial, Grouped_Polynomial* grouped_polynomi
 
     if(index==-2){
       fprintf(stderr,"error: monomial (");
-      for(j=0;j<polynomial.monomials[i].length;j++){
-	fprintf(stderr,"%d", polynomial.monomials[i].values[j]);
-	if(j<polynomial.monomials[i].length-1){
+      for(j=0;j<remainder.monomials[i].length;j++){
+	fprintf(stderr,"%d", remainder.monomials[i].values[j]);
+	if(j<remainder.monomials[i].length-1){
 	  fprintf(stderr,",");
 	}
       }
@@ -436,7 +436,7 @@ int simplify_grouped_polynomial(Grouped_Polynomial* polynomial){
 }
 
 
-// derive a flow equation with respect to an unknown variable
+// differentiate a flow equation with respect to an unknown variable
 // equivalent to DB.dl where dl are symbols for the derivatives of the indices in the flow equation with respect to the unknown variable
 // indices specifies the list of indices that depend on the variable
 int flow_equation_derivx(Grouped_Polynomial flow_equation, Int_Array indices, Grouped_Polynomial* dflow){
@@ -487,7 +487,7 @@ int flow_equation_derivx(Grouped_Polynomial flow_equation, Int_Array indices, Gr
 
 
 /*
-// derive a flow equation with respect to an index
+// differentiate a flow equation with respect to an index
 int flow_equation_deriv(Grouped_Polynomial flow_equation, int index, Grouped_Polynomial* output){
   int i,k;
   // temp list of indices
@@ -603,7 +603,7 @@ int grouped_polynomial_print(Grouped_Polynomial grouped_polynomial, char lhs_pre
     init_Char_Array(&buffer, STR_SIZE);
     coefficient_sprint(grouped_polynomial.coefs[i],&buffer,9,rhs_pre);
     if(buffer.length>0){
-      printf("%s",buffer.str);
+      printf("%s",char_array_to_str_noinit(&buffer));
     }
     free_Char_Array(buffer);
 
@@ -732,29 +732,33 @@ int char_array_to_Grouped_Polynomial(Char_Array str, Grouped_Polynomial* output)
 }
 
 
-// eValuate an equation on a vector
-int evaleq(RCC* rccs, Grouped_Polynomial poly){
+// evaluate an equation on a vector
+int evaleq(RCC out, RCC in, Grouped_Polynomial poly){
   int i;
-  long double* res=calloc((*rccs).length,sizeof(long double));
+  long double* res=calloc(out.length,sizeof(long double));
 
-  if((*rccs).length!=poly.length){
-    fprintf(stderr, "error: trying to evaluate an flow equation with %d components on an rcc with %d\n",poly.length,(*rccs).length);
+  if(in.length!=poly.length){
+    fprintf(stderr, "error: trying to evaluate a flow equation with %d components on an rcc with %d\n",poly.length,in.length);
+    exit(-1);
+  }
+  if(out.length!=poly.length){
+    fprintf(stderr, "error: trying to write the output of a flow equation with %d components on an rcc with %d\n",poly.length,out.length);
     exit(-1);
   }
 
-  // initialize vectors to 0
-  for(i=0;i<(*rccs).length;i++){
+  // initialize vectors to 0 in an auxiliary vector (to allow for out=in without interference)
+  for(i=0;i<in.length;i++){
     res[i]=0.;
   }
 
   // for each equation
   for(i=0;i<poly.length;i++){
-    evalcoef(*rccs, poly.coefs[i], res+i);
+    evalcoef(in, poly.coefs[i], res+i);
   }
 
   // copy res to rccs
-  for(i=0;i<(*rccs).length;i++){
-    (*rccs).values[i]=res[i];
+  for(i=0;i<out.length;i++){
+    out.values[i]=res[i];
   }
 
   // free memory
@@ -763,25 +767,29 @@ int evaleq(RCC* rccs, Grouped_Polynomial poly){
 
 }
 // evaluate an equation on a vector (using mpfr floats)
-int evaleq_mpfr(RCC_mpfr* rccs, Grouped_Polynomial poly){
+int evaleq_mpfr(RCC_mpfr out, RCC_mpfr in, Grouped_Polynomial poly){
   int i;
   mpfr_t* res;
 
-  if((*rccs).length!=poly.length){
-    fprintf(stderr, "error: trying to evaluate an flow equation with %d components on an rcc with %d\n",poly.length,(*rccs).length);
+  if(in.length!=poly.length){
+    fprintf(stderr, "error: trying to evaluate a flow equation with %d components on an rcc with %d\n",poly.length,in.length);
+    exit(-1);
+  }
+  if(out.length!=poly.length){
+    fprintf(stderr, "error: trying to write the output of a flow equation with %d components on an rcc with %d\n",poly.length,out.length);
     exit(-1);
   }
 
-  res=calloc((*rccs).length,sizeof(mpfr_t));
+  res=calloc(out.length,sizeof(mpfr_t));
 
   // for each equation
   for(i=0;i<poly.length;i++){
-    evalcoef_mpfr(*rccs, poly.coefs[i], res[i]);
+    evalcoef_mpfr(in, poly.coefs[i], res[i]);
   }
 
   // copy res to rccs
-  for(i=0;i<(*rccs).length;i++){
-    mpfr_set((*rccs).values[i], res[i], MPFR_RNDN);
+  for(i=0;i<out.length;i++){
+    mpfr_set(out.values[i], res[i], MPFR_RNDN);
     mpfr_clear(res[i]);
   }
 
@@ -791,3 +799,85 @@ int evaleq_mpfr(RCC_mpfr* rccs, Grouped_Polynomial poly){
 }
 
 
+// compose two flow equations (replace the rcc's of flow1 by the right hand side of flow2)
+int compose_flow_equations(Grouped_Polynomial flow1, Grouped_Polynomial flow2, Grouped_Polynomial* out){
+  if(flow1.length!=flow2.length){
+    fprintf(stderr, "error: trying to compose two flow equations of different size\n");
+    exit(-1);
+  }
+  int i,j,k;
+  Coefficient constant;
+
+  // init
+  init_Grouped_Polynomial(out, flow1.length);
+  (*out).length=flow1.length;
+
+  // init constant (so we can tell when the constant was not found)
+  constant.length=0;
+
+  // loop over rcc's
+  for(i=0;i<flow1.length;i++){
+    // set indices
+    (*out).indices[i]=flow1.indices[i];
+
+    // passthrough constant terms
+    if((*out).indices[i]<0){
+      int index=intlist_find_err(flow2.indices,flow2.length,(*out).indices[i]);
+      coefficient_cpy(flow2.coefs[index], (*out).coefs+i);
+      constant=flow2.coefs[index];
+      continue;
+    }
+
+    // init
+    init_Coefficient((*out).coefs+i, COEF_SIZE);
+
+    // loop over terms
+    for(j=0;j<flow1.coefs[i].length;j++){
+      Coefficient tmp_coef;
+
+      // init
+      init_Coefficient(&tmp_coef, COEF_SIZE);
+      // init factor
+      Int_Array tmp_factor;
+      init_Int_Array(&tmp_factor, MONOMIAL_SIZE);
+      // init denom
+      coef_denom denom;
+      // index should be that appearing in flow2
+      if(flow2.coefs[i].length<1){
+	fprintf(stderr,"error: composing two flow equations: the %d-th term in the flow equation is empty\n",flow1.indices[i]);
+	exit(-1);
+      }
+      denom.index=flow2.coefs[i].denoms[0].index;
+      denom.power=0;
+      // init num
+      Number tmp_num;
+      number_cpy(flow1.coefs[i].nums[j], &tmp_num);
+
+      // init coefficient with numerical prefactor
+      coefficient_append_noinit(tmp_factor, tmp_num, denom, &tmp_coef);
+
+      // loop over factors
+      for(k=0;k<flow1.coefs[i].factors[j].length;k++){
+	// multiply factors together
+	coefficient_prod_chain(flow2.coefs[intlist_find_err(flow2.indices,flow2.length,flow1.coefs[i].factors[j].values[k])], &tmp_coef);
+      }
+
+      // add to out
+      coefficient_concat_noinit(tmp_coef, (*out).coefs+i);
+    }
+
+  }
+
+  // simplify fractions
+  if(constant.length!=0){
+    for(i=0;i<(*out).length;i++){
+      if((*out).indices[i]>=0){
+	// reduce them to a common denominator (not much is gained from trying to simplify them)
+	coefficient_common_denominator(constant, (*out).coefs+i);
+	//coefficient_simplify_rational(constant, (*out).coefs+i);
+      }
+    }
+  }
+
+  return(0);
+}