/*
 * Title: evaluate.c
 * Author: Jed Marti KI7NNP
 * Description: Evaluate an inductor group using a result
 *   structure. We compute the RMS error by percent of the measured
 *   inductance to avoid over committment to bigger coils. Apologies
 *   for all the macros used here. These work with the equations
 *   defined in equations.h
 * Revision History: (Created Tue Jan 24 13:35:21 2023)
 *   Sat Jan 28 15:16:00 2023 - JBM quintic for diameter.
 *   Sun Jan 29 17:03:54 2023 - JBM New evaluation strategy with flags.
 *   Tue Feb  6 17:22:00 2024 - Fix RMSERROR calculation.
 */


#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>

#include "equations.h"
#include "indequ.h"


// The range table.
extern range R[];

// The symbol table.
extern symtab S[];

// Swapper. We reorder how we fix the number of tests so that the tail gets done as much as
// the start.
uint8_t swapper = 0;


// Doubtless (hopefully) gets open coded here.
double sqrd(double x) { return x * x; }
double cubed(double x) { return x * x * x; }
double quartd(double x) { return x * x * x * x; }
double quintd(double x) { return x * x * x * x * x; }



double evalInductorAbs(anInductor *l, result *r)
/*     We have a measured inductor l. Use its dimensions and data with the equation r to predict
 *   inductance. Then return the squared percentage error. For all the measured inductors, we
 *   will have the sum  |comp - L| / L.
 * Parameters:
 *   l: a measured inductor data.
 *   r: The coefficients of our polynomial.
 * External References: none
 * Returned values: Inductance in Henries or NORESULT (a super large value) indicating a value
 *   of 0.0 or division by 0.0.
 */
{
  double den;
  double num;

  // Compute denominator in anticipation of it being 0. Equations are defined in equations.h
  den = DEN;
  if (den == 0.0 || den == NAN) return NORESULT;
  num = NUM;
  if (num == 0.0 || num == NAN) return NORESULT;  
  return num / den;
}



double evalInductor(anInductor *l, result *r)
/*     Evaluate the current equation against that measured and return the percentage error.
 * Parameters:
 *   l: The measured inductance value for this coil.
 *   r: The parameter values for the equation.
 * External References:
 *   NORESULT: A constant double that's used to indicate equation failure.
 * Returned value: percentage error, greater than or equal to 0.0.
 * Errors detected: none
 */
{
  double v;
  v = evalInductorAbs(l, r);
  if (v == NORESULT) return NORESULT;
  #ifdef RMSERROR
    return sqrd((v - l->L) / l->L);
  #else
    return fabs(v - l->L) / l->L;
  #endif
}



void showResult(FILE *fout, result *r)
/*     Show the current values in a form that can be used in a program.
 * Parameters:
 *   fout: File to write to.
 *   r: The result to show.
 * External references: none
 * Returned values: none
 * Errors detected: none
 */
{
  fprintf(fout, "// Goodness = %g\n", r->goodness);
#define X(nme) fprintf(fout, "%5s = %20.14f\n", #nme, r->nme);
  VARS
#undef X
}



void evalInductors(anInductor *l, result *r)
/*     For a list of inductors and polynomial coefficients compute the "goodness" of the equation
 *   against this inductor. The values can range from 0.0 - perfect to pretty bad.
 * Parameters:
 *   l: A list of measured inductors.
 *   r: The polynomial coefficients.
 * External References: none
 * Returned values: none - modified the goodness value in r.
 * Errors detected: Only counts values that don't have 0 denominators. If no results returned
 *   then NORESULT is value.
 */
{
  double err = 0.0;
  double v;
  int32_t cnt = 0;

  // Now sum up the errors in all inductors.
  for (; l != NULL; l = l->next)
  {
    v = evalInductor(l, r);
    if (v != NORESULT) { err += v; ++cnt; }
    // printf("v=%g\n", v);
  }

  // If no results, then return a large bad value otherwise the mean of the errors.
  if (cnt == 0) r->goodness = NORESULT;
  #ifdef RMSERROR
    else r->goodness = sqrt(err / (double) cnt);
  #else
    else r->goodness = err / (double) cnt;
  #endif
}



static int64_t lowp(int32_t *lcnt,  double *lp, double r, int32_t i, int64_t cnt)
/*     Set the loop values for testing parameter minors.
 * Parameters:
 *   lcnt: Address of where to store count.
 *   lp: Address of where to store low value.
 *   r: Base value.
 *   i: index of variable in range table.
 *   cnt: Total count done so far.
 * External References:
 *   R[]: The range table.
 *   S[]: Symbol table. Looking for how many values to test to keep loop limits small.
 * Returned values: Values stored include the count of minor values and the lowest value. The
 *   increment comes from the range table.
 * Errors detected: none
 */
{
  int32_t v2;
  double  lw;
  double  hg;
  v2 = R[i].mcount / 2;
  lw = fmax(R[i].low, r - v2*R[i].minor);
  hg = fmin(R[i].high, r + v2*R[i].minor);
  *lcnt = (int32_t) ((hg - lw) / R[i].minor);
  *lp = lw;
  if (cnt * *lcnt >= S[s_minormax].value.iv || *lcnt == 0)
  {
    *lcnt = 1;
    *lp = r;
    return cnt;
  }
  return cnt * *lcnt;
}

// Earlier version only modifies those parameters flagged. Probably not a good idea even though
// it seems to work.
//  if ((flag & (1 << s_##p)) && cnt < S[s_minormax].value.iv)        \
//    cnt = lowp(&ilim##p, &lowParam->p, r->p, s_##p, cnt); \
//else { lowParam->p = r->p;  ilim##p = 1; }

// Macro for setting the limits and low value for a parameter. Also decides if there are too many
// to investigate and sets the limits to leave at just the major value.
#define setlim(p)                                \
if (cnt < S[s_minormax].value.iv) cnt = lowp(&ilim##p, &lowParam->p, r->p, s_##p, cnt); \
    else { lowParam->p = r->p;  ilim##p = 1; }

result *runMajor(anInductor *l, result *r)
/*     We have a major result. Scan all the parameters and their minor ranges collecting the
 *   best result. is on. Keep track of only the best solution for all the inductors. The library
 *   has stored r, not the new best result as a repeat of r would only give the same value.
 * Parameters:
 *   l: The list of inductors to test against.
 *   r: The set of major parameters to test.
 * External References: The macro above for setting limits.
 */
{
  result *bestResult;
  result *tempResult;
  result *lowParam;
  uint8_t anyMinor = 0;
  uint32_t flag;

  // The equation values computed whilst searching the limits.
#define X(nme) double nme;
  VARS
#undef X

  // The index for each parameter minor values.
#define X(nme) int32_t i##nme;
VARS
#undef X
    
  // Limits.
#define X(nme) int32_t ilim##nme;
    VARS
#undef X

  // Total count to do.
  int64_t cnt = 1;

  bestResult = newResult();
  bestResult->goodness = NORESULT;
  
  tempResult = newResult();
  lowParam = newResult();
  flag = r->flags;

  // Do the top values first.
  if (swapper & 0x01)
  {
#define X(nme) setlim(nme);
  VARS
  }

  // Otherwise do in reverse order.
  else
  {
    VARSI
    #undef X
  }
  ++swapper;
    

  cnt = 0;
#define X(nme) \
  for (i##nme = 0; i##nme < ilim##nme; ++i##nme)	\
  { \
    nme = lowParam->nme + ((double) i##nme) * R[s_##nme].minor; \
    tempResult->nme = nme;
VARS
#undef X
    // Now get the goodness into tempResult.
    evalInductors(l, tempResult);
    ++cnt;

    // If we have a good value then test.
    if (tempResult->goodness != NORESULT)
    {
      // If the best result is good and ours is good, keep the best.
      if (bestResult->goodness != NORESULT)
      {
        if (bestResult->goodness > tempResult->goodness) copyResult(bestResult, tempResult);
      }

      // bestResult is bad, use the new one which is good.
      else copyResult(bestResult, tempResult);
    }
#define X(nme) }
    VARS
#undef X
  // free the ones for the next round.
  freeResult(tempResult);
  freeResult(lowParam);
  //  printf("cnt = %ld\n", cnt);
  //  printf("*************************%g*************************\n", bestResult->goodness);
  return bestResult;
}



// EOF
