Lephisto/src/economy.c

/**
 * @file economy.c
 *
 * @brief Handles economy stuff.
 *
 * Economy is handled with Nodal Analysis. Systems are modelled as nodes,
 * jump routes are resistances and production is modeled as node intensity.
 * This is then solved with linear algebra after each time increment.
 */

#include <stdio.h>
#include <string.h>
#include <stdint.h>

#include "cs.h"

#include "lephisto.h"
#include "lxml.h"
#include "ldata.h"
#include "log.h"
#include "spfx.h"
#include "pilot.h"
#include "rng.h"
#include "space.h"
#include "ltime.h"
#include "economy.h"

#define XML_COMMODITY_ID  "Commodities" /* XML section identifier. */
#define XML_COMMODITY_TAG  "commodity"
#define COMMODITY_DATA   "../dat/commodity.xml"

/* Economy Nodal Analysis parameters. */
#define ECON_BASE_RES       30.     /**< Base resistance value for any system. */
#define ECON_SELF_RES       3.      /**< Additional resistance for the self node. */
#define ECON_FACTION_MOD    0.1     /**< Modifier on Base for faction standings. */
#define ECON_PROD_MODIFIER  500000. /**< Production modifier, divide production by this amount. */
#define ECON_PROD_VAR       0.01    /**< Define the variability of production. */

/* Commodity stack. */
static Commodity* commodity_stack = NULL;
static int commodity_nstack = 0;

/* Systems stack. */
extern StarSystem* systems_stack; /**< Star system stack. */
extern int systems_nstack;        /**< Number of star systems. */

/* Nodal analysis simulation for dynamic economies. */
static int  econ_initialized  = 0;    /**< Is economy system initialized? */
static int* econ_comm         = NULL; /**< Commodities to calculate. */
static int  econ_nprices      = 0;    /**< Number of prices to calculate. */
static cs*  econ_G            = NULL; /**< Admittance matrix. */

/* Commodity. */
static void commodity_freeOne(Commodity* com);
static int commodity_parse(Commodity* tmp, xmlNodePtr parent);
/* Economy. */
static double econ_calcJumpR(StarSystem* A, StarSystem* B);
static int econ_createGMatrix(void);
unsigned int economy_getPrice(const Commodity* com,
    const StarSystem* sys, const Planet* p); /* Externed in land.c. */

/* Convert credits to a usable string for displaying. */
/* str must have 10 characters allocated. */
void credits2str(char* str, unsigned int credits, int decimals) {
  if(decimals < 0)
    snprintf(str, 32, "%d", credits);
  else if(credits >= 1000000000)
    snprintf(str, 16, "%.*fB", decimals, (double)credits / 1000000000.);
  else if(credits >= 1000000)
    snprintf(str, 16, "%*fM", decimals, (double)credits / 1000000.);
  else if(credits >= 1000)
    snprintf(str, 16, "%.*fK", decimals, (double)credits / 1000.);
  else snprintf(str, 16, "%d", credits);
}

/* Get a commodity. */
Commodity* commodity_get(const char* name) {
  int i;
  for(i = 0; i < commodity_nstack; i++)
    if(strcmp(commodity_stack[i].name, name)==0)
      return &commodity_stack[i];

  WARN("Commodity '%s' not found in stack", name);
  return NULL;
}

/* Free a commodity. */
static void commodity_freeOne(Commodity* com) {
  if(com->name) free(com->name);
  if(com->description) free(com->description);

  /* Clear the memory. */
  memset(com, 0, sizeof(Commodity));
}

/**
 * @brief Loads a commodity.
 *    @param tmp Commodity to load data into.
 *    @param parent XML node to load from.
 *    @return Commodity loaded from parent.
 */
static int commodity_parse(Commodity* tmp, xmlNodePtr parent) {
  xmlNodePtr node;

  /* Clear memory. */
  memset(tmp, 0, sizeof(Commodity));

  tmp->name = (char*)xmlGetProp(parent, (xmlChar*)"name");
  if(tmp->name == NULL)
    WARN("Commodity from "COMMODITY_DATA" has invalid or noname");

  node = parent->xmlChildrenNode;

  do {
    xmlr_strd(node, "description",  tmp->description);
    xmlr_int( node, "price",        tmp->price);
  } while(xml_nextNode(node));
#if 0 /* Let's kill this for now. */
#define MELEMENT(o,s)if(o)WARN("Commodity '%s' missing '"s"' element",tmp->name)
  MELEMENT(tmp->high==0,        "high");
  MELEMENT(tmp->description==NULL,   "description");
  MELEMENT(tmp->medium==0,       "medium");
  MELEMENT(tmp->low==0,        "low");
#undef MELEMENT
#endif

  return 0;
}

/* Throw cargo out into space (graphically). */
void commodity_Jettison(int pilot, Commodity* com, int quantity) {
  (void)com;
  int i;
  Pilot* p;
  int n, effect;
  double px, py, bvx, bvy, r, a, vx, vy;

  p = pilot_get(pilot);

  n = MAX(1, RNG(quantity/10, quantity/5));
  px = p->solid->pos.x;
  py = p->solid->pos.y;
  bvx = p->solid->vel.x;
  bvy = p->solid->vel.y;

  for(i = 0; i < n; i++) {
    effect = spfx_get("cargo");

    /* Radial distribution gives much nicer results. */
    r = RNGF()*25 - 12.5;
    a = (double)RNG(0,259);
    vx = bvx + r*cos(a);
    vy = bvy + r*sin(a);

    /* Add the cargo effect. */
    spfx_add(effect, px, py, vx, vy, SPFX_LAYER_BACK);
  }
}

/* Init/exit */
int commodity_load(void) {
  uint32_t bufsize;
  char* buf = ldata_read(COMMODITY_DATA, &bufsize);

  xmlNodePtr node;
  xmlDocPtr doc = xmlParseMemory(buf, bufsize);

  node = doc->xmlChildrenNode; /* Commoditys node. */
  if(strcmp((char*)node->name, XML_COMMODITY_ID)) {
    ERR("Malformed "COMMODITY_DATA
        " file: Missing root element '"XML_COMMODITY_ID"'");
    return -1;
  }

  node = node->xmlChildrenNode; /* First faction node. */
  if(node == NULL) {
    ERR("Malformed "COMMODITY_DATA" file: does not contain elements");
    return -1;
  }

  do {
    if(xml_isNode(node, XML_COMMODITY_TAG)) {
      /* Make room for commodity. */
      commodity_stack = realloc(commodity_stack,
          sizeof(Commodity)*(++commodity_nstack));

      /* Load commodity. */
      commodity_parse(&commodity_stack[commodity_nstack-1], node);
     
      /* See if it should get added to commodity list. */
      if(commodity_stack[commodity_nstack-1].price > 0.) {
        econ_nprices++;
        econ_comm = realloc(econ_comm, econ_nprices * sizeof(int));
        econ_comm[econ_nprices-1] = commodity_nstack-1;
      }
    }
  } while((node = node->next));

  xmlFreeDoc(doc);
  free(buf);

  DEBUG("Loaded %d commodit%s",
        commodity_nstack, (commodity_nstack==1) ? "y" : "ies");

  return 0;
}

void commodity_free(void) {
  int i;
  for(i = 0; i < commodity_nstack; i++)
    commodity_freeOne(&commodity_stack[i]);
  free(commodity_stack);
  commodity_stack = NULL;
  commodity_nstack = 0;
}

/**
 * @brief Get the price of a good on a planet in a system.
 *    @param com Commodity to get price of.
 *    @param sys System to get price of commodity.
 *    @param p Planet to get price of commodity.
 *    @return The price of the commodity.
 */
unsigned int economy_getPrice(const Commodity* com,
    const StarSystem* sys, const Planet* p) {
  (void)p;
  int i, k;
  double price;

  /* Get position in stack. */
  k = com - commodity_stack;

  /* Find what commodity that is. */
  for(i = 0; i < econ_nprices; i++)
    if(econ_comm[i] == k)
      break;

  /* Check if found. */
  if(i >= econ_nprices) {
    WARN("Price for commodity '%s' not known.", com->name);
    return 0;
  }

  /* Calculate price. */
  price  = (double)com->price;
  price *= sys->prices[i];
  return (unsigned int) price;
}

/**
 * @brief Calculates the resistance between two star systems.
 *    @param A Star system to calculate the resistance between.
 *    @param B Star system to calculate the resistance between.
 *    @return Resistance between A and B.
 */
static double econ_calcJumpR(StarSystem* A, StarSystem* B) {
  double R;

  /* Set to base to ensure price change. */
  R = ECON_BASE_RES;

  /* Modify based on system conditions. */
  R += (A->nebu_density + B->nebu_density) / 1000.; /* Density shouldn't affect much. */
  R += (A->nebu_volatility + B->nebu_volatility) / 100.; /* Volatility should. */

  /* Modify based on global faction. */
  if((A->faction != -1) && (B->faction != -1)) {
    if(areEnemies(A->faction, B->faction))
      R += ECON_FACTION_MOD * ECON_BASE_RES;
    else if(areAllies(A->faction, B->faction))
      R -= ECON_FACTION_MOD * ECON_BASE_RES;
  }

  /* @todo Modify based on fleets. */
  return R;
}

/**
 * @brief Calculates the intensity in a system node.
 *
 * @todo Make it time/item dependent.
 */
static double econ_calcSysI(unsigned int dt, StarSystem* sys, int price) {
  (void)dt;
  (void)price;
  int i;
  double I;
  double prodfactor, p, pp;
  double ddt;
  Planet* planet;

  ddt = (double)(dt / LTIME_UNIT_LENGTH);

  /* Calculate production level. */
  p = 0.;
  for(i = 0; i < sys->nplanets; i++) {
    planet = sys->planets[i];
    if(planet_hasService(planet, PLANET_SERVICE_BASIC)) {
      /* Calculate production.
       *
       * We base off the current production.
       */
      prodfactor = planet->cur_prodfactor;
      /* Add a variability factor based on the gaussian distribution. */
      prodfactor += ECON_PROD_VAR *
        NormalInverse(RNGF()*0.90 + 0.05)*ddt;
      /* Add a tendency to return to the planets base production. */
      prodfactor -= ECON_PROD_VAR *
        (planet->cur_prodfactor - prodfactor)*ddt;
      /* Save for next iteration. */
      planet->cur_prodfactor = prodfactor;
      /* We base off the sqrt of the population otherwise it changes too fast. */
      p += prodfactor * sqrt(planet->population);

      /* Add it to the current system production. */
      p += pp;
    }
  }

  /* The intensity is basically the modified production. */
  I = p / ECON_PROD_MODIFIER;

  return I;
}

/**
 * @brief Create the admitance matrix.
 *    @return 0 on suceess.
 */
static int econ_createGMatrix(void) {
  int ret;
  int i, j;
  double R, Rsum;
  cs* M;
  StarSystem* sys;

  /* Create the matrix. */
  M = cs_spalloc(systems_nstack, systems_nstack, 1, 1, 1);
  if(M == NULL)
    ERR("Unable to create CSparse Matrix.");

  /* Fill the matrix. */
  for(i = 0; i < systems_nstack; i++) {
    sys = &systems_stack[i];
    Rsum = 0.;

    /* Set some values. */
    for(j = 0; j < sys->njumps; j++) {
      /* Get the resistance. */
      R = econ_calcJumpR(sys, &systems_stack[sys->jumps[j]]);
      R = 1./R; /* Must be inverted. */
      Rsum += R;

      /* Matrix is symetrical and non-diagonal is negative. */
      ret = cs_entry(M, i, sys->jumps[j], -R);
      if(ret != 1)
        WARN("Unable to enter CSparse Matrix Cell.");
      ret = cs_entry(M, sys->jumps[j], i, -R);
      if(ret != 1)
        WARN("Unable to enter CSparse Matrix Cell.");
    }

    /* Set the diagonal. */
    Rsum += 1./ECON_SELF_RES; /* We add a resistence for dampening. */
    cs_entry(M, i, i, Rsum);
  }

  /* Compress M matrix and put into G. */
  if(econ_G != NULL)
    cs_spfree(econ_G);
  econ_G = cs_compress(M);
  if(econ_G == NULL)
    ERR("Unable to create economy G Matrix.");

  /* Clean up. */
  cs_spfree(M);

  return 0;
}

/**
 * @brief Initializes the economy.
 *    @return 0 on success.
 */
int economy_init(void) {
  int i;

  /* Allocate price space. */
  for(i = 0; i < systems_nstack; i++) {
    if(systems_stack[i].prices != NULL)
      free(systems_stack[i].prices);
    systems_stack[i].prices = calloc(econ_nprices, sizeof(double));
  }

  /* Mark economy as initialized. */
  econ_initialized = 1;

  /* Refresh economy. */
  economy_refresh();

  return 0;
}

/**
 * @brief Regenerates the economy matrix. Should be used if the universe
 * changes in any permanent way.
 */
int economy_refresh(void) {
  /* Economy must be initialized. */
  if(econ_initialized == 0)
    return 0;

  /* Create the resistance matrix. */
  if(econ_createGMatrix())
    return -1;

  /* Initialize the prices. */
  economy_update(0);

  return 0;
}

/**
 * @brief Updates the economy.
 *    @param dt Deltatick in LTIME.
 */
int economy_update(unsigned int dt) {
  int ret;
  int i, j;
  double *X;
  double scale, offset;
  /*double min, max;*/

  /* Create the vector to solve the system. */
  X = malloc(sizeof(double)*systems_nstack);
  if(X == NULL)
    return -1;

  /* Calculate the results for each price set. */
  for(j = 0; j < econ_nprices; j++) {
    /* First we must load the vector with intensities. */
    for(i = 0; i < systems_nstack; i++)
      X[i] = econ_calcSysI(dt, &systems_stack[i], j);

    /* Solve the system. */
    ret = cs_lsolve(econ_G, X);
    if(ret != 1)
      WARN("Failed to solve the Economy System.");

    /* Get the minimum and maximum to scale. */
    /*
    min = +HUGE_VALF;
    max = -HUGE_VALF;
    for(i = 0; i < systems_nstack; i++) {
      if(X[i] < min)
        min = X[i];
      if(X[i] > max)
        max = X[i];
    }
    scale = 1. / (max - min);
    offset = 0.5 - min * scale;
    */

    /*
     * I'm not sure I like the filtering of the results, but it would take
     * much more work to get a sane system working without the need of post
     * filtering.
     */
    scale = 1.;
    offset = 1.;
    for(i = 0; i < systems_nstack; i++) {
      systems_stack[i].prices[j] = X[i] * scale + offset;
    }
  }

  /* Clean up. */
  free(X);
  return 0;
}

/**
 * @brief Destroys the economy.
 */
void economy_destroy(void) {
  int i;

  /* Clean up the prices in the systems stack. */
  for(i = 0; i < systems_nstack; i++) {
    if(systems_stack[i].prices != NULL) {
      free(systems_stack[i].prices);
      systems_stack[i].prices = NULL;
    }
  }

  /* Destroy the economy matrix. */
  if(econ_G != NULL) {
    cs_spfree(econ_G);
    econ_G = NULL;
  }

  /* Economy is now deinitialized. */
  econ_initialized = 0;
}