LUT.cc

Go to the documentation of this file.

// ----------------------------------------------------------------------------
// openCARP is an open cardiac electrophysiology simulator.
//
// Copyright (C) 2020 openCARP project
//
// This program is licensed under the openCARP Academic Public License (APL)
// v1.0: You can use and redistribute it and/or modify it in non-commercial
// academic environments under the terms of APL as published by the openCARP
// project v1.0, or (at your option) any later version. Commercial use requires
// a commercial license (info@opencarp.org).
//
// This program is distributed without any warranty; see the openCARP APL for
// more details.
//
// You should have received a copy of the openCARP APL along with this program
// and can find it online: http://www.opencarp.org/license
// ----------------------------------------------------------------------------

#include "LUT.h"
#include "basics.h"
#include "MULTI_ION_IF.h"

#include <string.h>
#include <assert.h>
#include <cmath>

namespace limpet {

using ::opencarp::dupstr;
using ::opencarp::f_open;
using ::opencarp::FILE_SPEC;
using ::opencarp::is_big_endian;
using ::opencarp::log_msg;

/*  create and initialize LUT
 *
 * \param p pointer to lookup table structure
 * \param rows  number of lookup variables
 * \param mn  lower boundary of LUT
 * \param mx  upper boundary of LUT
 * \param res resolution of LUT
 * \param name  table name
 *
 */
void LUT_alloc(LUT *plut,int cols,float mn,float mx,float res,const char *name, Target target)
{
  plut->name = dupstr(name);

  // table bounds
  plut->cols   = cols;
  plut->mn     = mn;
  plut->mx     = mx;
  plut->res    = res;
  plut->step   = 1/res;

  // in terms of indices
  plut->mn_ind = (int)(mn * plut->step);
  plut->mx_ind = (int)(mx * plut->step);
  plut->rows   = plut->mx_ind - plut->mn_ind + 1;

  // build a table and offset it
  plut->tab  = build_matrix_ns<LUT_data_t>( plut->rows, plut->cols, sizeof(LUT_data_t), target);
  plut->tab -= plut->mn_ind;
}


/*  dump LUT to file
 *
 * \param p pointer to lookup table structure
 * \param fname filename to store LUT dump
 *
 * \return 0 if successull, -1 otherwise
 */
int LUT_dump( LUT *plut, const char *fname )
{
  FILE_SPEC fdump = f_open(fname, "wb");

  if (!fdump)
    return -1;

  // write endianness and table dimensions first
  int endian_flg = is_big_endian();
  fwrite(&endian_flg, sizeof(int), 1, fdump->fd);
  fwrite(&plut->rows, sizeof(plut->rows), 1, fdump->fd);
  fwrite(&plut->cols, sizeof(plut->cols), 1, fdump->fd);

  for (int i=plut->mn_ind; i<=plut->mx_ind; i++)
    for (int j=0; j<plut->cols; j++)
      fwrite(&plut->tab[i][j], sizeof(LUT_data_t), 1, fdump->fd);

  f_close(fdump);

  return 0;
}


void destroy_lut( LUT *plut, Target target )
{
  if ( plut != NULL && plut->tab != NULL )  {
    if (plut->tab != NULL) {
      plut->tab += plut->mn_ind;
      deallocate_on_target<LUT_data_t>(target, plut->tab[0]);
      deallocate_on_target(target, plut->tab);
    }
    free(plut->name);
  }
  if( plut )
    memset( plut, 0, sizeof(LUT) );
}




#define MAX_LUT_NUMINF 10       

int check_LUT( LUT* lut )
{
  int retval = 0;
  int numinf = 0;

  for ( int i=lut->mn_ind; i<=lut->mx_ind; i++ )
    for ( int j=0; j<lut->cols; j++ )
      if ( ! std::isfinite(lut->tab[i][j])) {
        log_msg(0, 2, 0, "LUT WARNING: %s=%g produces %g in entry number %d!\n",lut->name,i*lut->res,lut->tab[i][j],j);
        retval = 1;
    if( numinf++ > MAX_LUT_NUMINF ) {
          log_msg(0, 3, 0, "suppressing further errors!!!\n" );
      return retval;
    }
      }

  return retval;
}

void IIF_warn(const int wv, const char error[]) {
  // log_msg(_nc_logf, 2, 0, "danger [tm = %f]: %s : local node = %d, global node = %d\n",
  //         current_global_time(), error, wv, current_global_node(wv) ); //FIXME
}

void LUT_problem( LUT *lt, double val, int wv, const char *tabname )
{
  // struct exception_type except =  { .exit=0 };
  char error[5000]; //FIXME, replace me with a C++ string.

  if ( std::isfinite(val) ) {
#ifdef CHATTY_LUT
    snprintf(error, sizeof error, "bounds exceeded for %s-table = %g (limits: %g - %g)",
            tabname, val, lt->mn, lt->mx);
    IIF_warn(wv, error);
#endif
    return;
  } else if (std::isinf(val)) {
    snprintf(error, sizeof error, "inf passed to LUT_index() for %s-table", tabname);
    IIF_warn(wv, error);
  } else if ( std::isnan(val) ) {
    snprintf(error, sizeof error, "NaN passed to LUT_index() for %s-table", tabname);
    IIF_warn(wv, error);
  }
}

#ifndef IMP_FAST
int LUT_index( LUT *tab, GlobalData_t val, int locind )
{
  int indx =  (int)(tab->step*val);

  if (indx < tab->mn_ind) {
    indx = tab->mn_ind;
  } else if (indx > tab->mx_ind) {
    indx = tab->mx_ind;
  } else {
    return indx;
  }

  LUT_problem(tab, val, locind, tab->name );
  return indx;
}
#endif

int LUT_out_of_bounds(LUT *tab, GlobalData_t val) {
  return (val < tab->mn || val > tab->mx) ;
}

inline LUT_data_t LUT_interp( LUT *t, int i, int j, GlobalData_t x )
{
  return (1.-x)*t->tab[i][j] + x*t->tab[i+1][j];
}


inline LUT_data_t LUT_derror( LUT *t, int idx, GlobalData_t x )
{
  return  (x-idx*t->res)/t->res;
}

LUT_data_t LUT_interpRow(LUT *const tab,GlobalData_t val,int i,LUT_data_t* row)
{
  int    idx  = LUT_index(tab, val, i);
  GlobalData_t derr = LUT_derror(tab, idx, val);
  if (LUT_out_of_bounds(tab, val)) {
    for (int j=0;j<tab->cols;j++)
      row[j] = tab->tab[idx][j];
  } else {
    for (int j=0;j<tab->cols;j++)
      row[j] = LUT_interp(tab, idx, j, derr );
  }
  return derr;
}


LUT_data_t*
LUT_row( LUT *lut, GlobalData_t val, int locind )
{
  return lut->tab[LUT_index( lut, val, locind )];
}

#ifdef __cplusplus
extern "C"
{
#endif  // ifdef __cplusplus


void LUT_interpRow_n_elements(char *table, char *val_ptr, int offset,
                              int distance, int index, int n, char* row_ptr,
                              int lut_numelements) {
  if (n <= 0)
    return;
  LUT *tab = (LUT *)table;
  GlobalData_t *val = (GlobalData_t *)(val_ptr + offset);
  LUT_data_t* row = (LUT_data_t*)row_ptr;
  if (distance == 1) {
    for (unsigned i = 0; i < n; ++i)
      LUT_interpRow(tab, *(val + i), i + index*n,
                    row + (i * (lut_numelements)));
  } else {
    for (unsigned i = 0; i < n; ++i)
      LUT_interpRow(tab, *(val + i*distance + index*n), i + index*n,
                    row + (i * (lut_numelements)));
  }
}

#ifdef HAS_MLIR_CPU_MODEL
void compute_LUT_interpRow_mlir_8xf64(double , double, double, double, int, double, double*, double *, bool, int, double, char*);
void compute_LUT_interpRow_mlir_4xf64(double , double, double, double, int, double, double*, double *, bool, int, double, char*);
void compute_LUT_interpRow_mlir_2xf64(double , double, double, double, int, double, double*, double *, bool, int, double, char*);

void LUT_interpRow_mlir(char *table, int i, char* row_ptr,
                        int lut_num_elements, int vector_size)
{
  LUT *const tab = (LUT *const)table;
    LUT_data_t* row = (LUT_data_t*)row_ptr;
    
    if (lut_num_elements > 0)
    {
    #ifndef IMP_FAST
      if (vector_size == 8)
      {
        compute_LUT_interpRow_mlir_8xf64((double) tab->step, (double) tab->mn, (double) tab->mx, (double) tab->res, (int) tab->cols, (double) tab->mn_ind, (double*)tab->tab[tab->mn_ind], (double *) row, true, i, (double) tab->mx_ind, table);
      }
      else if (vector_size == 4)
      {
        compute_LUT_interpRow_mlir_4xf64((double) tab->step, (double) tab->mn, (double) tab->mx, (double) tab->res, (int) tab->cols, (double) tab->mn_ind, (double*)tab->tab[tab->mn_ind], (double *) row, true, i, (double) tab->mx_ind, table);
      }
      else if (vector_size == 2)
      {
        compute_LUT_interpRow_mlir_2xf64((double) tab->step, (double) tab->mn, (double) tab->mx, (double) tab->res, (int) tab->cols, (double) tab->mn_ind, (double*)tab->tab[tab->mn_ind], (double *) row, true, i, (double) tab->mx_ind, table);
      }
       
    #else
      if (vector_size == 8)
      {
        compute_LUT_interpRow_mlir_8xf64((double) tab->step, (double) tab->mn, (double) tab->mx, (double) tab->res, (int) tab->cols, (double) tab->mn_ind, (double*)tab->tab[tab->mn_ind], (double *) row, false, i, (double) tab->mx_ind, table);
      }
      else if (vector_size == 4)
      {
        compute_LUT_interpRow_mlir_4xf64((double) tab->step, (double) tab->mn, (double) tab->mx, (double) tab->res, (int) tab->cols, (double) tab->mn_ind, (double*)tab->tab[tab->mn_ind], (double *) row, false, i, (double) tab->mx_ind, table);
      }
      else if (vector_size == 2)
      {
        compute_LUT_interpRow_mlir_2xf64((double) tab->step, (double) tab->mn, (double) tab->mx, (double) tab->res, (int) tab->cols, (double) tab->mn_ind, (double*)tab->tab[tab->mn_ind], (double *) row, false, i, (double) tab->mx_ind, table);
      }
    #endif      
    }
    
}


void LUT_problem_mlir( char *table, GlobalData_t val, int locind)
{
  LUT *const tab = (LUT *const)table;
  LUT_problem(tab, val, locind, tab->name );
}
#endif

#ifdef __cplusplus
}
#endif  // ifdef __cplusplus

}  // namespace limpet