text:
/****************************************************************************/
/* Copyright 2005-2006, Francis Russell                                     */
/*                                                                          */
/* Licensed under the Apache License, Version 2.0 (the License);            */
/* you may not use this file except in compliance with the License.         */
/* You may obtain a copy of the License at                                  */
/*                                                                          */
/*     http://www.apache.org/licenses/LICENSE-2.0                           */
/*                                                                          */
/* Unless required by applicable law or agreed to in writing, software      */
/* distributed under the License is distributed on an AS IS BASIS,          */
/* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. */
/* See the License for the specific language governing permissions and      */
/* limitations under the License.                                           */
/*                                                                          */
/****************************************************************************/

#ifndef DESOLA_TG_OBJECTS_HPP
#define DESOLA_TG_OBJECTS_HPP

#include <typeinfo>
#include <cassert>
#include <cstddef>
#include <cstring>
#include <string>
#include <map>
#include <utility>
#include <algorithm>
#include <functional>
#include <TaskGraph>
#include <boost/functional/hash.hpp>
#include <boost/function.hpp>
#include "Desola_tg_fwd.hpp"
#include "TaskGraphWrappers.hpp"

namespace desola
{

namespace detail
{

template<typename T_element>
class TGScalar
{
public:
  virtual const TGScalarExpr<T_element> getExpression() const = 0;
  virtual void setExpression(const TGScalarExpr<T_element>& e) = 0;
  virtual void addExpression(const TGScalarExpr<T_element>& e) = 0;
  virtual InternalScalar<T_element>* createInternalRep() const = 0;
  virtual bool isParameter() const = 0;
  virtual void addParameterMappings(InternalScalar<T_element>& internal, ParameterHolder& params) const = 0;
  virtual std::size_t hashValue() const = 0;
  virtual bool matches(const TGScalar& s) const = 0;
  virtual void createTaskGraphVariable() = 0;
  virtual ~TGScalar() {}
};

template<typename T_element>
class TGVector
{
public:
  virtual const TGScalarExpr<T_element> getExpression(const tg::TaskExpression& row) const = 0;
  virtual void setExpression(const tg::TaskExpression& row, const TGScalarExpr<T_element>& e) = 0;
  virtual void addExpression(const tg::TaskExpression& row, const TGScalarExpr<T_element>& e) = 0;
  virtual std::size_t getRows() const = 0;
  virtual InternalVector<T_element>* createInternalRep() const = 0;
  virtual bool isParameter() const = 0;
  virtual void addParameterMappings(InternalVector<T_element>& internal, ParameterHolder& params) const = 0;
  virtual std::size_t hashValue() const = 0;
  virtual bool matches(const TGVector& v) const = 0;
  virtual void createTaskGraphVariable() = 0;
  virtual ~TGVector() {}
};

template<typename T_element>
class TGMatrix
{
public:
  typedef boost::function<void (NameGenerator& generator, const tg::TaskExpression& row, const tg::TaskExpression& col, const TGScalarExpr<T_element>& value)> MatrixIterationCallback;

  virtual const TGScalarExpr<T_element> getExpression(NameGenerator& generator, const tg::TaskExpression& row, const tg::TaskExpression& col) const = 0;
  virtual void setExpression(NameGenerator& generator, const tg::TaskExpression& row, const tg::TaskExpression& col, const TGScalarExpr<T_element>& e) = 0;
  virtual void addExpression(NameGenerator& generator, const tg::TaskExpression& row, const tg::TaskExpression& col, const TGScalarExpr<T_element>& e) = 0;
  virtual std::size_t getRows() const = 0;
  virtual std::size_t getCols() const = 0;
  virtual void iterateDense(NameGenerator& generator, MatrixIterationCallback& callback) const = 0;
  virtual void iterateSparse(NameGenerator& generator, MatrixIterationCallback& callback) const = 0;
  virtual InternalMatrix<T_element>* createInternalRep() const = 0;
  virtual bool isParameter() const = 0;
  virtual void addParameterMappings(InternalMatrix<T_element>& internal, ParameterHolder& params) const = 0;
  virtual std::size_t hashValue() const = 0;
  virtual bool matches(const TGMatrix& m) const = 0;
  virtual void createTaskGraphVariable() = 0;
  virtual ~TGMatrix() {}
};

template<typename T_element>
class TGScalarExpr
{
  const tg::TaskExpression expr;

public:
  TGScalarExpr() : expr(T_element())
  {
  }

  TGScalarExpr(const tg::TaskExpression& e) : expr(e)
  {
  }

  inline const tg::TaskExpression getExpression() const
  {
    return expr;
  }

  inline const TGScalarExpr mul(const TGScalarExpr& right) const
  {
    return TGScalarExpr(expr*right.expr);
  }

  inline const TGScalarExpr div(const TGScalarExpr& right) const
  {
    return TGScalarExpr(expr/right.expr);
  }

  inline const TGScalarExpr add(const TGScalarExpr& right) const
  {
    return TGScalarExpr(expr+right.expr);
  }

  inline const TGScalarExpr sub(const TGScalarExpr& right) const
  {
    return TGScalarExpr(expr-right.expr);
  }

  inline const TGScalarExpr sqrt() const
  {
    return TGScalarExpr(tSqrt(expr));
  }

  inline const TGScalarExpr abs() const
  {
    return TGScalarExpr(tAbs(expr));
  }

  inline const TGScalarExpr negate() const
  {
    return TGScalarExpr(-expr);
  }
};

template<typename T_element>
class TGConventionalScalar : public TGScalar<T_element>
{
  friend class Mapper;

private:
  const bool parameter;
  const std::string name;
  TaskScalarVariableWrapper<T_element> value;

  static inline const std::string getPrefix()
  {
    return std::string("convScalar");
  }

public:
  TGConventionalScalar(NameGenerator& generator, const ConventionalScalar<T_element>& internal, const bool hasData) :
    parameter(true), name(generator.getName(getPrefix())), value(true, name)
  {
  }

  TGConventionalScalar(const bool param, NameGenerator& generator, const ExprNode<scalar, T_element>& s) : parameter(param),
    name(generator.getName(getPrefix())), value(param, name)
  {
  }

  const TGScalarExpr<T_element> getExpression() const
  {
    return TGScalarExpr<T_element>(*value);
  }

  void setExpression(const TGScalarExpr<T_element>& e)
  {
    *value = e.getExpression();
  }

  void addExpression(const TGScalarExpr<T_element>& e)
  {
    *value += e.getExpression();
  }

  virtual InternalScalar<T_element>* createInternalRep() const
  {
    assert(parameter);
    return new ConventionalScalar<T_element>();
  }

  virtual bool isParameter() const
  {
    return parameter;
  }

  virtual std::size_t hashValue() const
  {
    const char* nodeTypeString = typeid(*this).name();
    std::size_t seed = boost::hash_range(nodeTypeString, nodeTypeString+strlen(nodeTypeString));
    boost::hash_combine(seed, parameter);
    boost::hash_combine(seed, name);
    return seed;
  }

  virtual bool matches(const TGScalar<T_element>& s) const
  {
    if (typeid(s) == typeid(*this))
    {
      const TGConventionalScalar& right = static_cast<const TGConventionalScalar&>(s);
      return parameter == right.parameter &&
	     name == right.name;
    }
    else
    {
      return false;
    }
  }

  virtual void createTaskGraphVariable()
  {
    value.instantiate();
  }

  virtual void addParameterMappings(InternalScalar<T_element>& internal, ParameterHolder& params) const
  {
    assert(parameter);
    Mapper mapper(*this, params);
    internal.accept(mapper);
  }

  class Mapper : public InternalScalarVisitor<T_element>
  {
  private:
    const TGConventionalScalar<T_element>& internal;
    ParameterHolder& parameterHolder;

  public:
    Mapper(const TGConventionalScalar<T_element>& i, ParameterHolder& p) : internal(i), parameterHolder(p)
    {
    }

    virtual void visit(ConventionalScalar<T_element>& s)
    {
      assert(s.isAllocated());
      parameterHolder.addParameter(internal.name, s.getValue());
    }
  };
};

template<typename T_element>
class TGConventionalVector : public TGVector<T_element>
{
private:
  const bool parameter;
  const std::string name;
  const std::size_t rows;
  TaskArrayWrapper<T_element, 1> value;

  static inline const std::string getPrefix()
  {
    return std::string("convVector");
  }

public:
  TGConventionalVector(NameGenerator& generator, const ConventionalVector<T_element>& internal, const bool hasData) :
    parameter(true), name(generator.getName(getPrefix())), rows(internal.getRowCount()), value(true, name, rows)
  {
  }

  TGConventionalVector(const bool param, NameGenerator& generator, const ExprNode<vector, T_element>& v) :
    parameter(param), name(generator.getName(getPrefix())), rows(v.getRowCount()), value(param, name, rows)
  {
  }

  const TGScalarExpr<T_element> getExpression(const tg::TaskExpression& row) const
  {
    return TGScalarExpr<T_element>((*value)[row]);
  }

  void setExpression(const tg::TaskExpression& row, const TGScalarExpr<T_element>& e)
  {
    (*value)[row] = e.getExpression();
  }

  void addExpression(const tg::TaskExpression& row, const TGScalarExpr<T_element>& e)
  {
    (*value)[row] += e.getExpression();
  }

  virtual std::size_t getRows() const
  {
    return rows;
  }

  virtual InternalVector<T_element>* createInternalRep() const
  {
    assert(parameter);
    return new ConventionalVector<T_element>(rows);
  }

  virtual bool isParameter() const
  {
    return parameter;
  }

  virtual void addParameterMappings(InternalVector<T_element>& internal, ParameterHolder& params) const
  {
    assert(parameter);
    Mapper mapper(*this, params);
    internal.accept(mapper);
  }

  virtual std::size_t hashValue() const
  {
    const char* nodeTypeString = typeid(*this).name();
    std::size_t seed = boost::hash_range(nodeTypeString, nodeTypeString+strlen(nodeTypeString));
    boost::hash_combine(seed, parameter);
    boost::hash_combine(seed, rows);
    boost::hash_combine(seed, name);
    return seed;
  }

  virtual bool matches(const TGVector<T_element>& v) const
  {
    if (typeid(v) == typeid(*this))
    {
      const TGConventionalVector& right = static_cast<const TGConventionalVector&>(v);
      return parameter == right.parameter &&
             name == right.name &&
	     rows == right.rows;
    }
    else
    {
      return false;
    }
  }

  virtual void createTaskGraphVariable()
  {
    value.instantiate();
  }

  class Mapper : public InternalVectorVisitor<T_element>
  {
  private:
    const TGConventionalVector<T_element>& internal;
    ParameterHolder& parameterHolder;

  public:
    Mapper(const TGConventionalVector<T_element>& i, ParameterHolder& p) : internal(i), parameterHolder(p)
    {
    }

    virtual void visit(ConventionalVector<T_element>& v)
    {
      assert(v.isAllocated());
      parameterHolder.addParameter(internal.name, v.getValue());
    }
  };

};

template<typename T_element>
class TGConventionalMatrix : public TGMatrix<T_element>
{
private:
  const bool parameter;
  const std::string name;
  const std::size_t rows;
  const std::size_t cols;
  TaskArrayWrapper<T_element, 2> value;

  static inline const std::string getPrefix()
  {
    return std::string("convMatrix");
  }

public:
  typedef typename TGMatrix<T_element>::MatrixIterationCallback MatrixIterationCallback;

  TGConventionalMatrix(NameGenerator& generator, ConventionalMatrix<T_element>& internal, const bool hasData) :
    parameter(true), name(generator.getName(getPrefix())),  rows(internal.getRowCount()), cols(internal.getColCount()), value(true, name, rows, cols)
  {
  }

  TGConventionalMatrix(const bool param, NameGenerator& generator, const ExprNode<matrix, T_element>& m) :
    parameter(param), name(generator.getName(getPrefix())), rows(m.getRowCount()), cols(m.getColCount()), value(param, name, rows, cols)
  {
  }

  const TGScalarExpr<T_element> getExpression(NameGenerator& generator, const tg::TaskExpression& row, const tg::TaskExpression& col) const
  {
    return TGScalarExpr<T_element>((*value)[row][col]);
  }

  void setExpression(NameGenerator& generator, const tg::TaskExpression& row, const tg::TaskExpression& col, const TGScalarExpr<T_element>& e)
  {
    (*value)[row][col] = e.getExpression();
  }

  void addExpression(NameGenerator& generator, const tg::TaskExpression& row, const tg::TaskExpression& col, const TGScalarExpr<T_element>& e)
  {
    (*value)[row][col] += e.getExpression();
  }

  virtual std::size_t getRows() const
  {
    return rows;
  }

  virtual std::size_t getCols() const
  {
    return cols;
  }

  virtual void iterateDense(NameGenerator& generator, MatrixIterationCallback& callback) const
  {
    using namespace tg;

    tVarNamed(unsigned, i, generator.getName("ConvMatrix_row").c_str());
    tVarNamed(unsigned, j, generator.getName("ConvMatrix_col").c_str());

    tFor(i, 0u, rows-1)
    {
      tFor(j, 0u, cols-1)
      {
        callback(generator, i, j, getExpression(generator, i, j));
      }
    }
  }

  virtual void iterateSparse(NameGenerator& generator, MatrixIterationCallback& callback) const
  {
    iterateDense(generator, callback);
  }

  virtual InternalMatrix<T_element>* createInternalRep() const
  {
    assert(parameter);
    return new ConventionalMatrix<T_element>(rows, cols);
  }

  virtual bool isParameter() const
  {
    return parameter;
  }

  virtual void addParameterMappings(InternalMatrix<T_element>& internal, ParameterHolder& params) const
  {
    assert(parameter);
    Mapper mapper(*this, params);
    internal.accept(mapper);
  }

  virtual std::size_t hashValue() const
  {
    const char* nodeTypeString = typeid(*this).name();
    std::size_t seed = boost::hash_range(nodeTypeString, nodeTypeString+strlen(nodeTypeString));
    boost::hash_combine(seed, parameter);
    boost::hash_combine(seed, rows);
    boost::hash_combine(seed, cols);
    boost::hash_combine(seed, name);
    return seed;
  }

  virtual bool matches(const TGMatrix<T_element>& m) const
  {
    if (typeid(m) == typeid(*this))
    {
      const TGConventionalMatrix& right = static_cast<const TGConventionalMatrix&>(m);
      return parameter == right.parameter &&
             name == right.name &&
             rows == right.rows &&
	     cols == right.cols;
    }
    else
    {
      return false;
    }
  }

  virtual void createTaskGraphVariable()
  {
    value.instantiate();
  }

  class Mapper : public InternalMatrixVisitor<T_element>
  {
  private:
    const TGConventionalMatrix<T_element>& internal;
    ParameterHolder& parameterHolder;

  public:
    Mapper(const TGConventionalMatrix<T_element>& i, ParameterHolder& p) : internal(i), parameterHolder(p)
    {
    }

    virtual void visit(ConventionalMatrix<T_element>& m)
    {
      assert(m.isAllocated());
      parameterHolder.addParameter(internal.name, m.getValue());
    }

    virtual void visit(CRSMatrix<T_element>& m)
    {
      // FIXME: Handle this error better
      assert(0 && "Attempted to map parameters for CRSMatrix instead of a ConventionalMatrix");
    }
  };
};

template<typename T_element>
class TGCRSMatrix : public TGMatrix<T_element>
{
public:
  typedef typename TGMatrix<T_element>::MatrixIterationCallback MatrixIterationCallback;

private:
  static inline std::string getColIndPrefix()
  {
    return std::string("colInd");
  }

  static inline std::string getRowPtrPrefix()
  {
    return std::string("rowPtr");
  }

  static inline std::string getValPrefix()
  {
    return std::string("val");
  }

  const bool parameter;
  const std::string col_ind_name;
  const std::string row_ptr_name;
  const std::string val_name;
  const std::size_t nnz;
  const std::size_t rows;
  const std::size_t cols;
  TaskArrayWrapper<int, 1> col_ind;
  TaskArrayWrapper<int, 1> row_ptr;
  TaskArrayWrapper<T_element, 1> val;
  const CRSMatrix<T_element>* possibleData;

  void specialisedIterateSparse(NameGenerator& generator, MatrixIterationCallback& callback) const
  {
    using namespace tg;

    assert(possibleData != NULL);
    RowLengthStatistics stats(*possibleData);

    std::vector< std::pair<std::size_t, std::size_t> > freqToRowLengths;

    BOOST_FOREACH(const RowLengthStatistics::value_type& rowFreq, std::make_pair(stats.begin(), stats.end()))
    {
      freqToRowLengths.push_back(std::make_pair(rowFreq.second, rowFreq.first));
    }

    std::sort(freqToRowLengths.begin(), freqToRowLengths.end(), std::greater< std::pair<std::size_t, std::size_t> >());

    tVarNamed(unsigned, row, generator.getName("row").c_str());
    tVarNamed(unsigned, rowLength, generator.getName("rowLength").c_str());
    tVarNamed(unsigned, valOffset, generator.getName("valOffset").c_str());
    tVarNamed(unsigned, valPtr, generator.getName("valPtr").c_str());
    tVarNamed(unsigned, valPtrStart, generator.getName("valPtrStart").c_str());
    tVarNamed(unsigned, valPtrEnd, generator.getName("valPtrEnd").c_str());

    tFor(row, 0u, rows-1)
    {
      valPtrStart = (*row_ptr)[row];
      valPtrEnd = (*row_ptr)[row+1];
      rowLength = valPtrEnd - valPtrStart;

      const double fractionToSpecialise = 0.95;
      std::size_t rowsSpecialised = 0;

      for(std::size_t index = 0; index < freqToRowLengths.size(); ++index)
      {
        rowsSpecialised += freqToRowLengths[index].first;
        const std::size_t constRowLength = freqToRowLengths[index].second;
        tIf(rowLength == constRowLength)
        {
          tFor(valOffset, 0, constRowLength - 1)
          {
            callback(generator, row, (*col_ind)[valPtrStart + valOffset], TGScalarExpr<T_element>((*val)[valPtrStart + valOffset]));
          }
          tContinue;
        }

        if (static_cast<double>(rowsSpecialised)/stats.numRows() >= fractionToSpecialise)
          break;
      }

      tFor(valPtr, valPtrStart, valPtrEnd - 1)
      {
        callback(generator, row, (*col_ind)[valPtr], TGScalarExpr<T_element>((*val)[valPtr]));
      }
    }
  }

public:
  TGCRSMatrix(NameGenerator& generator, CRSMatrix<T_element>& internal, const bool hasData) :
    parameter(true), col_ind_name(generator.getName(getColIndPrefix())), row_ptr_name(generator.getName(getRowPtrPrefix())),
    val_name(generator.getName(getValPrefix())), nnz(internal.nnz()), rows(internal.getRowCount()),
    cols(internal.getColCount()), col_ind(true, col_ind_name, nnz), row_ptr(true, row_ptr_name, internal.row_ptr_size()),
    val(true, val_name, nnz), possibleData(hasData ? &internal : NULL)
  {
  }

  // TaskGraph code shouldn't generate CRS intermediates or returns so comment out this constructor for now
  //TGCRSMatrix(bool param, NameGenerator& generator, const ExprNode<matrix, T_element>& m) :  parameter(param), col_ind_name(generator.getName(getColIndPrefix())),
  //                                                                        row_ptr_name(generator.getName(getRowPtrPrefix())), val_name(generator.getName(getValPrefix())),
  //                                                                        nnz(0), rows(m.getRowCount()), cols(m.getColCount()), col_ind(param, col_ind_name, 0),
  //                                                                        row_ptr(param, row_ptr_name, 0), val(param, val_name, 0)
  //{
  //  //FIXME: Handle this error better
  //  assert(0 && "CRS sparse matrices should not be generated as intermediates or return values by TaskGraph code");
  //}

  const TGScalarExpr<T_element> getExpression(NameGenerator& generator, const tg::TaskExpression& row, const tg::TaskExpression& col) const
  {
    using namespace tg;

    tVarNamed(int, valPtr, generator.getName("valPtr").c_str());
    tVarTemplateTypeNamed(T_element, result, generator.getName("elementVal").c_str());

    result = T_element();

    tFor(valPtr, (*row_ptr)[row], (*row_ptr)[row+1]-1)
    {
      tIf((*col_ind)[valPtr] == col)
        result = (*val)[valPtr];
    }

    return TGScalarExpr<T_element>(result);
  }

  void setExpression(NameGenerator& generator, const tg::TaskExpression& row, const tg::TaskExpression& col, const TGScalarExpr<T_element>& e)
  {
    //FIXME: Handle this problem better
    assert(0 && "Editing a CRS Matrix inside TaskGraph code is not implemented.");
  }

  void addExpression(NameGenerator& generator, const tg::TaskExpression& row, const tg::TaskExpression& col, const TGScalarExpr<T_element>& e)
  {

    //FIXME: Handle this problem better
    assert(0 && "Editing a CRS Matrix inside TaskGraph code is not implemented.");
  }

  virtual std::size_t getRows() const
  {
    return rows;
  }

  virtual std::size_t getCols() const
  {
    return cols;
  }

  virtual InternalMatrix<T_element>* createInternalRep() const
  {
    assert(parameter);
    return new CRSMatrix<T_element>(rows, cols);
  }

  virtual bool isParameter() const
  {
    return parameter;
  }

  virtual void iterateDense(NameGenerator& generator, MatrixIterationCallback& callback) const
  {
    //FIXME: Implement me
    assert(0 && "Not yet implemented dense iteration over CRS matrix.");
  }

  virtual void iterateSparse(NameGenerator& generator, MatrixIterationCallback& callback) const
  {
    using namespace tg;
    const ConfigurationManager& configurationManager(ConfigurationManager::getConfigurationManager());

    if (configurationManager.singleForLoopSparseIterationEnabled())
    {
      tVarNamed(unsigned, valPtr, generator.getName("valPtr").c_str());
      tVarNamed(unsigned, currentRow, generator.getName("currentRow").c_str());

      currentRow = 0u;

      tFor(valPtr, 0u, nnz-1)
      {
        tWhile((*row_ptr)[currentRow+1] == valPtr)
          ++currentRow;

        callback(generator, currentRow, (*col_ind)[valPtr], TGScalarExpr<T_element>((*val)[valPtr]));
      }
    }
    else if (possibleData != NULL && configurationManager.sparseSpecialisationEnabled())
    {
      specialisedIterateSparse(generator, callback);
    }
    else
    {
      tVarNamed(unsigned, row, generator.getName("row").c_str());
      tVarNamed(unsigned, valPtr, generator.getName("valPtr").c_str());
      tVarNamed(unsigned, valPtrStart, generator.getName("valPtrStart").c_str());
      tVarNamed(unsigned, valPtrEnd, generator.getName("valPtrEnd").c_str());

      tFor(row, 0u, rows-1)
      {
        valPtrStart = (*row_ptr)[row];
        valPtrEnd = (*row_ptr)[row+1];

        tFor(valPtr, valPtrStart, valPtrEnd - 1)
        {
          callback(generator, row, (*col_ind)[valPtr], TGScalarExpr<T_element>((*val)[valPtr]));
        }
      }
    }
  }

  virtual void addParameterMappings(InternalMatrix<T_element>& internal, ParameterHolder& params) const
  {
    assert(parameter);
    Mapper mapper(*this, params);
    internal.accept(mapper);
  }

  virtual std::size_t hashValue() const
  {
    const char* nodeTypeString = typeid(*this).name();
    std::size_t seed = boost::hash_range(nodeTypeString, nodeTypeString+strlen(nodeTypeString));
    boost::hash_combine(seed, parameter);
    boost::hash_combine(seed, rows);
    boost::hash_combine(seed, cols);
    boost::hash_combine(seed, col_ind_name);
    boost::hash_combine(seed, row_ptr_name);
    boost::hash_combine(seed, val_name);
    return seed;
  }

  virtual bool matches(const TGMatrix<T_element>& m) const
  {
    if (typeid(m) == typeid(*this))
    {
      const TGCRSMatrix& right = static_cast<const TGCRSMatrix&>(m);
      return parameter == right.parameter &&
             col_ind_name == right.col_ind_name &&
             row_ptr_name == right.row_ptr_name &&
             val_name == right.val_name &&
             rows == right.rows &&
	     cols == right.cols;
    }
    else
    {
      return false;
    }
  }

  virtual void createTaskGraphVariable()
  {
    col_ind.instantiate();
    row_ptr.instantiate();
    val.instantiate();
  }

  class Mapper : public InternalMatrixVisitor<T_element>
  {
  private:
    const TGCRSMatrix<T_element>& internal;
    ParameterHolder& parameterHolder;

  public:
    Mapper(const TGCRSMatrix<T_element>& i, ParameterHolder& p) : internal(i), parameterHolder(p)
    {
    }

    virtual void visit(CRSMatrix<T_element>& m)
    {
      assert(m.isAllocated());
      parameterHolder.addParameter(internal.col_ind_name, m.get_col_ind());
      parameterHolder.addParameter(internal.row_ptr_name, m.get_row_ptr());
      parameterHolder.addParameter(internal.val_name, m.get_val());
    }

    virtual void visit(ConventionalMatrix<T_element>& m)
    {
      // FIXME: Handle this error better
      assert(0 && "Attempted to map parameters for ConventionalMatrix instead of a CRSMatrix");
    }
  };
};

template<typename T_element>
class TGMatrixGen : public InternalMatrixVisitor<T_element>
{
private:
  NameGenerator& generator;
  TGMatrix<T_element>* result;
  const bool hasData;

public:
  TGMatrixGen(NameGenerator& g, const bool _hasData) : generator(g), hasData(_hasData)
  {
  }

  void visit(ConventionalMatrix<T_element>& m)
  {
    result = new TGConventionalMatrix<T_element>(generator, m, hasData);
  }

  void visit(CRSMatrix<T_element>& m)
  {
    result = new TGCRSMatrix<T_element>(generator, m, hasData);
  }

  TGMatrix<T_element>* getResult() const
  {
    assert(result != NULL);
    return result;
  }
};

template<typename T_element>
class TGVectorGen : public InternalVectorVisitor<T_element>
{
private:
  NameGenerator& generator;
  TGVector<T_element>* result;
  const bool hasData;

public:
  TGVectorGen(NameGenerator& g, const bool _hasData) : generator(g), hasData(_hasData)
  {
  }

  void visit(ConventionalVector<T_element>& m)
  {
    result = new TGConventionalVector<T_element>(generator, m, hasData);
  }

  TGVector<T_element>* getResult() const
  {
    assert(result != NULL);
    return result;
  }
};

template<typename T_element>
class TGScalarGen : public InternalScalarVisitor<T_element>
{
private:
  NameGenerator& generator;
  TGScalar<T_element>* result;
  const bool hasData;

public:
  TGScalarGen(NameGenerator& g, const bool _hasData) : generator(g), hasData(_hasData)
  {
  }

  void visit(ConventionalScalar<T_element>& m)
  {
    result = new TGConventionalScalar<T_element>(generator, m, hasData);
  }

  TGScalar<T_element>* getResult() const
  {
    assert(result != NULL);
    return result;
  }
};

}

}
#endif

max_stars_repo_name: FrancisRussell/desola

max_issues_repo_name: FrancisRussell/desola

max_forks_repo_name: FrancisRussell/desola


text:
/*
Metaheuristic Optimization Using Population-Based Simulated Annealing.

Copyright (c) 2021 Gabriele Gilardi

    func            Function to minimize
    LB              Lower boundaries of the search space
    UB              Upper boundaries of the search space
    nPop            Number of agents (population)
    epochs          Number of iterations
    nMove           Number of neighbours of a state evaluated at each epoch
    T0              Initial temperature
    alphaT          Temperature reduction rate
    sigma0          Initial standard deviation used to search the neighboroud
                    of a state (given as a fraction of the search space)
    alphaS          Standard deviation reduction rate
    prob            Probability the dimension of a state is changed
    IntVar          List of indexes specifying which variable should be treated
                    as integer
    normalize       Specifies if the search space should be normalized
    args            Tuple containing any parameter that needs to be passed to
                    the function

    Dimensions:
    (nVar)          LB, UB, LB_orig, UB_orig, sigma, best_pos
    (nPop, nVar)    agent_pos, neigh_pos, agent_pos_orig, neigh_pos_orig
    (nPop)          agent_cost, neigh_cost
    (epochs)        F
    (0-nVar)        IntVar
*/

#include <random>
#include <Eigen/Dense>
#include <iostream>

using namespace std;
using namespace Eigen;

/* Structure used to pass the parameters */
struct Parameters {
    int nPop;
    int epochs;
    int nMove;
    double T0;
    double alphaT;
    double sigma0;
    double alphaS;
    double prob;
    bool normalize;
    ArrayXi IntVar;
    ArrayXXd args;
    int seed;
};

/* Structure used to return the results */
struct Results {
    double best_cost;
    ArrayXXd best_pos;
    ArrayXd F;
    double T;
    ArrayXXd sigma;
};

/* Minimize a function using simulated annealing */
Results sa(ArrayXd (*func)(ArrayXXd, ArrayXXd), ArrayXXd LB, ArrayXXd UB,
           Parameters p)
{
    int nVar, nIntVar, idx;
    double T, prob_swap, delta, best_cost;
    Results res;

    /*Eigen declarations*/
    Index r_min, c_min;
    ArrayXd agent_cost, neigh_cost, F;
    ArrayXXd LBe, UBe, LBe_orig, UBe_orig, sigma, agent_pos, neigh_pos,
             agent_pos_orig, neigh_pos_orig, best_pos, rn, flips;

    /* Random generator and probability distributions */
    mt19937 generator(p.seed);
    uniform_real_distribution<double> unif(0.0, 1.0);
    normal_distribution<double> norm(0.0, 1.0);

    nVar = LB.size();
    nIntVar = p.IntVar.size();

    /* Create boundaries for all agents*/
    LBe = LB.replicate(p.nPop, 1);
    UBe = UB.replicate(p.nPop, 1);

    /* Normalize search space */
    if (p.normalize) {
        LBe_orig = LBe;
        UBe_orig = UBe;
        LBe.setZero(p.nPop, nVar);
        UBe.setOnes(p.nPop, nVar);
    }

    T = p.T0;                           // Temperature
    sigma = p.sigma0 * (UBe - LBe);     // Standard deviation of each dimension

    /* Initial position of each agent */
    rn.setZero(p.nPop, nVar);
    for (int i=0; i<p.nPop; i++) {
        for (int j=0; j<nVar; j++) {
            rn(i, j) = unif(generator);
        }
    }
    agent_pos = LBe + rn * (UBe - LBe);

    /* Correct for any integer variable */
    for (int j=0; j<nIntVar; j++) {
        idx = p.IntVar(j);
        agent_pos.col(idx) = round(agent_pos.col(idx));
    }

    /* Initial cost of each agent */
    if (p.normalize) {
        agent_pos_orig = LBe_orig + agent_pos * (UBe_orig - LBe_orig);
        agent_cost = func(agent_pos_orig, p.args);
    }
    else {
        agent_cost = func(agent_pos, p.args);
    }

    /* Initial (overall) best position/cost */
    best_cost = agent_cost.minCoeff(&r_min, &c_min);
    best_pos = agent_pos.row(r_min);

    /* Main loop (T = const) */
    neigh_pos.setZero(p.nPop, nVar);
    F.setZero(p.epochs);
    for (int epoch=0; epoch<p.epochs; epoch++) {

        /* Sub-loop (search the neighboroud of a state) */
        for (int move=0; move<p.nMove; move++) {

            /* Randomly decide in which dimension to search */
            for (int i=0; i<p.nPop; i++) {
                for (int j=0; j<nVar; j++) {
                    rn(i, j) = unif(generator);
                }
            }
            flips = (rn <= p.prob).cast<double>();

            /* Create each agent's neighbours */
            for (int i=0; i<p.nPop; i++) {
                for (int j=0; j<nVar; j++) {
                    rn(i, j) = norm(generator);
                }
            }
            neigh_pos = agent_pos + flips * rn * sigma;

            /* Correct for any integer variable */
            for (int j=0; j<nIntVar; j++) {
                idx = p.IntVar(j);
                neigh_pos.col(idx) = round(neigh_pos.col(idx));
            }

            /* Impose position boundaries */
            neigh_pos = neigh_pos.max(LBe);
            neigh_pos = neigh_pos.min(UBe);

            for (int j=0; j<nIntVar; j++) {
                idx = p.IntVar(j);
                neigh_pos.col(idx) = neigh_pos.col(idx).max(ceil(LBe.col(idx)));
                neigh_pos.col(idx) = neigh_pos.col(idx).min(floor(UBe.col(idx)));
            }

            /* Evaluate the cost of each agent's neighbour */
            if (p.normalize) {
                neigh_pos_orig = LBe_orig + neigh_pos * (UBe_orig - LBe_orig);
                neigh_cost = func(neigh_pos_orig, p.args);
            }
            else {
                neigh_cost = func(neigh_pos, p.args);
            }

            /* Decide if each agent will change its state */
            for (int i=0; i<p.nPop; i++) {

                /* Swap states if the neighbour state is better ... */
                if (neigh_cost(i) <= agent_cost(i)) {
                    agent_cost(i) = neigh_cost(i);
                    agent_pos.row(i) = neigh_pos.row(i);
                }

                /* ... or decide probabilistically */
                else {

                    /* Acceptance probability */
                    delta = (neigh_cost(i) - agent_cost(i)) / agent_cost(i);
                    prob_swap = exp(-delta / T);

                    /* Randomly swap states */
                    if (unif(generator) <= prob_swap) {
                        agent_cost(i) = neigh_cost(i);
                        agent_pos.row(i) = neigh_pos.row(i);
                    }
                }

                /* Update the (overall) best position/cost */
                best_cost = agent_cost.minCoeff(&r_min, &c_min);
                best_pos = agent_pos.row(r_min);
            }
        }

        /* Save the best cost for this epoch */
        F(epoch) = best_cost;

        /* Cooling scheduling */
        T = p.alphaT * T;

        /* Random neighboroud search schedule */
        sigma = p.alphaS * sigma;
    }

    /* De-normalize */
    if (p.normalize) {
        best_pos = LB + best_pos * (UB - LB);
        sigma = sigma * (UBe_orig - LBe_orig);
    }

    /* Copy solution */
    res.best_cost = best_cost;
    res.T = T;
    res.best_pos = best_pos;
    res.sigma = sigma.row(0);
    res.F = F;

    return res;
}

ext: cpp

max_stars_repo_name: gabrielegilardi/SimulatedAnnealing

max_issues_repo_name: gabrielegilardi/SimulatedAnnealing

max_forks_repo_name: gabrielegilardi/SimulatedAnnealing


text:
#include "vtkboneFaimVersion5InputReader.h"
#include "vtkboneFiniteElementModel.h"
#include "vtkboneMaterialTable.h"
#include "vtkboneLinearIsotropicMaterial.h"
#include "vtkboneLinearOrthotropicMaterial.h"
#include "vtkboneConstraint.h"
#include "vtkboneConstraintCollection.h"
#include "vtkboneSolverParameters.h"
#include "vtkbone_version.h"
#include "vtkObjectFactory.h"
#include "vtkIntArray.h"
#include "vtkFloatArray.h"
#include "vtkDoubleArray.h"
#include "vtkCharArray.h"
#include "vtkCellArray.h"
#include "vtkPointData.h"
#include "vtkCellData.h"
#include "vtkStreamingDemandDrivenPipeline.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkInformationIntegerKey.h"
#include "vtkInformationDoubleKey.h"
#include "vtkInformationStringVectorKey.h"
#include <boost/format.hpp>
#include <boost/tokenizer.hpp>
#include <boost/lexical_cast.hpp>
#include <vector>
#include <string>
#include <sstream>
#include <ios>
#include <cassert>

using boost::format;
using boost::lexical_cast;
using std::string;
using std::ios_base;

#ifdef read
#warning That is weird. Why is read defined?
#undef read
#endif

vtkStandardNewMacro(vtkboneFaimVersion5InputReader);

class FAIMInputReaderHelper
{
public:
  FAIMInputReaderHelper(const char* FileName, istream& in);
  int IsIgnorable(const string& l);
  // Returns the next line, that is non-ignorable
  string GetNextLine();
  string& GetLine();
  void SplitLine(const string& l, std::vector<string>& tokens);
  void GetLineAndParse(int N, std::vector<string>& tokens);
  string GetNextToken();
  void GetNTokens(int N, std::vector<string>& tokens);
  template <typename T> void ReadArray(vtkIdType N, T* data);

  const char* FileName;
  istream& in;
  long line_count;
  std::vector<string> UnusedTokens;
  string line;
  string next_line;
};

//-----------------------------------------------------------------------
FAIMInputReaderHelper::FAIMInputReaderHelper(
  const char* FileName_,
  istream& in_)
:
  FileName (FileName_),
  in (in_),
  line_count (-1)
{}

//-----------------------------------------------------------------------
int FAIMInputReaderHelper::IsIgnorable(const string& l)
{
  if (l.empty())
    {
    return 1;
    }
  if ((l[0] == '\n') || (l[0] == '#'))
    {
    return 1;
    }
  return 0;
}

//-----------------------------------------------------------------------
string FAIMInputReaderHelper::GetNextLine()
{
  string l;
  while (1)
    {
    getline(in, l);
    this->line_count += 1;
    if (in.eof())
      {
      return string();  // Return blank string
      }
    if (in.fail())
      {
      throw ios_base::failure((format("IO Error reading %s at line %d")
                               % this->FileName % this->line_count).str());
      }
    if (!this->IsIgnorable(l))
      {
      break;
      }
    }
  return l;
}

//-----------------------------------------------------------------------
string& FAIMInputReaderHelper::GetLine()
{
  if (this->UnusedTokens.size() > 0)
    {
    throw ios_base::failure((format("Left over values at end of line %d")
                              % this->line_count).str());
    }
  if (this->line_count == -1)
    {
    this->next_line = this->GetNextLine();
    }
  this->line = this->next_line;
  if (this->line.empty())
    {
    throw ios_base::failure((format("Unexpected end of file at line %d")
                            % this->line_count).str());
    }
  this->next_line = this->GetNextLine();
  return this->line;
}

//-----------------------------------------------------------------------
void FAIMInputReaderHelper::SplitLine(
  const string& l,
  std::vector<string>& tokens)
{
  typedef boost::tokenizer<boost::char_separator<char> > tokenizer;
  boost::char_separator<char> separators(" \t,");

  tokens.clear();
  tokenizer tok(l, separators);
  for(tokenizer::iterator beg=tok.begin(); beg!=tok.end();++beg)
    {
    tokens.push_back(*beg);
    }
}

//-----------------------------------------------------------------------
void FAIMInputReaderHelper::GetLineAndParse(
  int N,
  std::vector<string>& tokens)
{
  this->GetLine();
  this->SplitLine(this->line, tokens);
  if (N != -1 && tokens.size() != N)
    {
    throw ios_base::failure((format("Unexpected number of tokens at line %d")
                             % this->line_count).str());
    }
}

//-----------------------------------------------------------------------
string FAIMInputReaderHelper::GetNextToken()
{
  if (this->UnusedTokens.size() == 0)
    {
    this->GetLine();
    this->SplitLine(this->line, this->UnusedTokens);
    if (this->UnusedTokens.size() == 0)
      {
      throw ios_base::failure((format("Insufficient number of values at line %d")
                              % this->line_count).str());
      }
    }
  string token = this->UnusedTokens[0];
  this->UnusedTokens.erase(this->UnusedTokens.begin());
  return token;
}

//-----------------------------------------------------------------------
void FAIMInputReaderHelper::GetNTokens(int N, std::vector<string>& tokens)
{
  tokens.clear();
  while (tokens.size() != N)
    {
    tokens.push_back(this->GetNextToken());
    }
}

//-----------------------------------------------------------------------
template <typename T>
void FAIMInputReaderHelper::ReadArray(vtkIdType N, T* data)
{
  for (vtkIdType i=0; i<N; i++)
    {
    data[i] = lexical_cast<T>(this->GetNextToken());
    }
}

//-----------------------------------------------------------------------
vtkboneFaimVersion5InputReader::vtkboneFaimVersion5InputReader()
{
  this->FileName = NULL;
  this->SetNumberOfInputPorts(0);
}

//-----------------------------------------------------------------------
vtkboneFaimVersion5InputReader::~vtkboneFaimVersion5InputReader()
{
  this->SetFileName(0);
}

//----------------------------------------------------------------------------
void vtkboneFaimVersion5InputReader::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os,indent);
  os << indent << "File Name: "
     << (this->FileName ? this->FileName : "(none)") << "\n";
}

//----------------------------------------------------------------------------
int vtkboneFaimVersion5InputReader::RequestData(
                               vtkInformation *,
                               vtkInformationVector **,
                               vtkInformationVector *outputVector)
{

  vtkInformation *outInfo = outputVector->GetInformationObject(0);

  vtkboneFiniteElementModel *output = vtkboneFiniteElementModel::SafeDownCast(
    outInfo->Get(vtkDataObject::DATA_OBJECT()));
  if (!output)
  {
    vtkErrorMacro("No output object.");
    return 0;
  }

  // Make sure we have a file to read.
  if (!this->FileName)
  {
    vtkErrorMacro("FileName not set.");
    return 0;
  }

  // Open the input file.
  std::ifstream fin (this->FileName);
  if (!fin)
    {
    vtkErrorMacro("Error opening file " << this->FileName);
    return 0;
    }
  vtkDebugMacro(<<"\n Reading " << this->FileName << " ...");

  std::ostringstream history;
  history << "Model read from Faim version 5 input file \"" << this->FileName << "\" using vtkbone version " << VTKBONE_VERSION;
  output->AppendHistory(history.str().c_str());

  // Variables that we are going to read.
  vtkIdType numberOfElements;
  vtkIdType numberOfNodes;
  double spacing[3];
  int numberOfMaterials;

  vtkInformation* info = output->GetInformation();

  FAIMInputReaderHelper reader(this->FileName, fin);

  try
    {

    // --------------------------------------------------------------
    // Read Model Parameters

    std::vector<string> tokens;
    reader.GetLineAndParse(8, tokens);
    numberOfElements = lexical_cast<vtkIdType>(tokens[1]);
    numberOfNodes = lexical_cast<vtkIdType>(tokens[2]);
    int nip = lexical_cast<int>(tokens[3]);
    int nodof = lexical_cast<int>(tokens[4]);
    if (nodof != 3)
      {
      vtkErrorMacro("Unsupported number of degrees of freedom in " << this->FileName);
      return 0;
      }
    int nod = lexical_cast<int>(tokens[5]);
    if (nod != 8)
      {
      vtkErrorMacro("Unsupported number of nodes per cell in " << this->FileName);
      return 0;
      }
    int nst = lexical_cast<int>(tokens[6]);
    int ndim = lexical_cast<int>(tokens[7]);
    if (ndim != 3)
      {
      vtkErrorMacro("Unsupported number of dimensions in " << this->FileName);
      return 0;
      }

    reader.GetLineAndParse(3, tokens);
    for (int i=0; i<3; i++)
      {
      spacing[i] = lexical_cast<double>(tokens[i]);
      }

    reader.GetLineAndParse(2, tokens);
    int nprops = lexical_cast<int>(tokens[0]);
    numberOfMaterials = lexical_cast<int>(tokens[1]);
    if (nprops != 9)
      {
      vtkErrorMacro(<<"Only material property type 9 (orthotropic) currently handled.");
      return 0;
      }

    // --------------------------------------------------------------
    // Read Material Table

    vtkSmartPointer<vtkboneMaterialTable> materialTable =
                               vtkSmartPointer<vtkboneMaterialTable>::New();
    for (int i=0; i<numberOfMaterials; i++)
      {
      reader.GetLineAndParse(10, tokens);
      int materialIndex = lexical_cast<int>(tokens[0]);
      if (materialIndex < 1)
        {
        vtkErrorMacro(<<"Invalid index in material table line " << reader.line_count);
        return 0;
        }
      if (tokens[2] == tokens[1] && tokens[3] == tokens[1] &&
          tokens[6] == tokens[4] && tokens[5] == tokens[4])
        {
        vtkSmartPointer<vtkboneLinearIsotropicMaterial> material = vtkSmartPointer<vtkboneLinearIsotropicMaterial>::New();
        try
          {
          material->SetYoungsModulus(lexical_cast<double>(tokens[1]));
          material->SetPoissonsRatio(lexical_cast<double>(tokens[4]));
          }
        catch(boost::bad_lexical_cast &)
          {
          vtkErrorMacro(<<"Unable to parse material properties at line " << reader.line_count);
          return 0;
          }
        materialTable->AddMaterial(materialIndex, material);
        }
      else
        {
        vtkSmartPointer<vtkboneLinearOrthotropicMaterial> material = vtkSmartPointer<vtkboneLinearOrthotropicMaterial>::New();
        try
          {
          material->SetYoungsModulusX(lexical_cast<double>(tokens[1]));
          material->SetYoungsModulusY(lexical_cast<double>(tokens[2]));
          material->SetYoungsModulusZ(lexical_cast<double>(tokens[3]));
          material->SetPoissonsRatioXY(lexical_cast<double>(tokens[4]));
          material->SetPoissonsRatioYZ(lexical_cast<double>(tokens[5]));
          material->SetPoissonsRatioZX(lexical_cast<double>(tokens[6]));
          material->SetShearModulusXY(lexical_cast<double>(tokens[7]));
          material->SetShearModulusYZ(lexical_cast<double>(tokens[8]));
          material->SetShearModulusZX(lexical_cast<double>(tokens[9]));
          }
        catch(boost::bad_lexical_cast &)
          {
          vtkErrorMacro(<<"Unable to parse material properties at line " << reader.line_count);
          return 0;
          }
        materialTable->AddMaterial(materialIndex, material);
        }
      }
    output->SetMaterialTable(materialTable);

    // --------------------------------------------------------------
    // Read Solver Parameters

    reader.GetLineAndParse(6, tokens);
    double convergenceTolerance = lexical_cast<double>(tokens[0]);
    vtkboneSolverParameters::CONVERGENCE_TOLERANCE()->Set (info, convergenceTolerance);
    int iterationLimit = lexical_cast<int>(tokens[1]);
    vtkboneSolverParameters::MAXIMUM_ITERATIONS()->Set (info, iterationLimit);

    // --------------------------------------------------------------
    // Read Nodes (Points)

    vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
    points->SetDataTypeToFloat();
    points->SetNumberOfPoints(numberOfNodes);
    reader.ReadArray(numberOfNodes*3,
                     reinterpret_cast<float*>(points->GetVoidPointer(0)));
    output->SetPoints(points);

    // --------------------------------------------------------------
    // Read Elements (Cells)

    vtkSmartPointer<vtkCellArray> cells = vtkSmartPointer<vtkCellArray>::New();
    // The following allocation is exact
    const int pointsPerCell = 8;
    cells->Allocate(cells->EstimateSize(numberOfElements, pointsPerCell));
    {  // scope
    vtkSmartPointer<vtkIdList> pointIds = vtkSmartPointer<vtkIdList>::New();
    pointIds->SetNumberOfIds(pointsPerCell);
    for (vtkIdType newCellId=0; newCellId < numberOfElements; newCellId++)
      {
      reader.GetLineAndParse(pointsPerCell, tokens);
      // Convert from FAIM topology to VTK_VOXEL topolgy.
      // Note: Also convert from 1-indexed to 0-indexed.
      pointIds->SetId(0, lexical_cast<vtkIdType>(tokens[0])-1);
      pointIds->SetId(1, lexical_cast<vtkIdType>(tokens[3])-1);
      pointIds->SetId(2, lexical_cast<vtkIdType>(tokens[4])-1);
      pointIds->SetId(3, lexical_cast<vtkIdType>(tokens[7])-1);
      pointIds->SetId(4, lexical_cast<vtkIdType>(tokens[1])-1);
      pointIds->SetId(5, lexical_cast<vtkIdType>(tokens[2])-1);
      pointIds->SetId(6, lexical_cast<vtkIdType>(tokens[5])-1);
      pointIds->SetId(7, lexical_cast<vtkIdType>(tokens[6])-1);
      cells->InsertNextCell(pointIds);
      }
    }
    output->SetCells(VTK_VOXEL, cells);

    // --------------------------------------------------------------
    // Read Material Indices (Cell Scalars)

    vtkSmartPointer<vtkIntArray> scalars = vtkSmartPointer<vtkIntArray>::New();
    scalars->SetNumberOfTuples(numberOfElements);
    reader.ReadArray(numberOfElements,
                     reinterpret_cast<int*>(scalars->GetVoidPointer(0)));
    output->GetCellData()->SetScalars(scalars);

    // --------------------------------------------------------------
    // Read Fixed Nodes Constraint Set

    reader.GetLineAndParse(1, tokens);
    vtkIdType numberOfFixedFreedoms = lexical_cast<vtkIdType>(tokens[0]);
    if (numberOfFixedFreedoms > 0)
      {
      vtkSmartPointer<vtkIdTypeArray> ids = vtkSmartPointer<vtkIdTypeArray>::New();
      ids->Allocate(numberOfFixedFreedoms);
      vtkSmartPointer<vtkCharArray> senses = vtkSmartPointer<vtkCharArray>::New();
      senses->Allocate(numberOfFixedFreedoms);
      vtkSmartPointer<vtkDoubleArray> values = vtkSmartPointer<vtkDoubleArray>::New();
      values->Allocate(numberOfFixedFreedoms);
      for (vtkIdType i=0; i<numberOfFixedFreedoms; i++)
        {
        reader.GetLineAndParse(4, tokens);
        vtkIdType nodeId = lexical_cast<vtkIdType>(tokens[0]) - 1;
        if (lexical_cast<int>(tokens[1]) == 0)
          {
          ids->InsertNextValue(nodeId);
          senses->InsertNextValue(0);
          values->InsertNextValue(0.0);
          }
        if (lexical_cast<int>(tokens[2]) == 0)
          {
          ids->InsertNextValue(nodeId);
          senses->InsertNextValue(1);
          values->InsertNextValue(0.0);
          }
        if (lexical_cast<int>(tokens[3]) == 0)
          {
          ids->InsertNextValue(nodeId);
          senses->InsertNextValue(2);
          values->InsertNextValue(0.0);
          }
        }
      if (ids->GetNumberOfTuples() != senses->GetNumberOfTuples() ||
          ids->GetNumberOfTuples() != values->GetNumberOfTuples())
        {
        vtkErrorMacro(<<"Internal Error");
        return 0;
        }
      output->ApplyBoundaryCondition(ids, senses, values, "FixedConstraints");
      }

    // --------------------------------------------------------------
    // Read Force Constraints (on nodes)

    reader.GetLineAndParse(1, tokens);
    vtkIdType numberOfForceConstraints = lexical_cast<vtkIdType>(tokens[0]);
    if (numberOfForceConstraints > 0)
      {
      vtkSmartPointer<vtkIdTypeArray> ids = vtkSmartPointer<vtkIdTypeArray>::New();
      ids->Allocate(numberOfForceConstraints);
      vtkSmartPointer<vtkCharArray> senses = vtkSmartPointer<vtkCharArray>::New();
      senses->SetName("SENSE");
      senses->Allocate(numberOfForceConstraints);
      vtkSmartPointer<vtkFloatArray> values = vtkSmartPointer<vtkFloatArray>::New();
      values->SetName("VALUE");
      values->Allocate(numberOfForceConstraints);
      for (vtkIdType i=0; i<numberOfForceConstraints; i++)
        {
        reader.GetLineAndParse(3, tokens);
        ids->InsertNextValue(lexical_cast<vtkIdType>(tokens[0]) - 1);
        senses->InsertNextValue(lexical_cast<vtkIdType>(tokens[1]) - 1);
        values->InsertNextValue(lexical_cast<float>(tokens[2]));
        }
      vtkSmartPointer<vtkboneConstraint> constraint = vtkSmartPointer<vtkboneConstraint>::New();
      constraint->SetName("ForceConstraints");
      constraint->SetIndices(ids);
      constraint->SetConstraintType(vtkboneConstraint::FORCE);
      constraint->SetConstraintAppliedTo(vtkboneConstraint::NODES);
      constraint->GetAttributes()->AddArray(senses);
      constraint->GetAttributes()->AddArray(values);
      output->GetConstraints()->AddItem(constraint);
      }

    // --------------------------------------------------------------
    // Read Displacement Constraints

    reader.GetLineAndParse(1, tokens);
    vtkIdType numberOfDisplacementConstraints = lexical_cast<vtkIdType>(tokens[0]);
    if (numberOfDisplacementConstraints > 0)
      {
      vtkSmartPointer<vtkIdTypeArray> ids = vtkSmartPointer<vtkIdTypeArray>::New();
      ids->Allocate(numberOfDisplacementConstraints);
      vtkSmartPointer<vtkCharArray> senses = vtkSmartPointer<vtkCharArray>::New();
      senses->SetName("SENSE");
      senses->Allocate(numberOfDisplacementConstraints);
      vtkSmartPointer<vtkFloatArray> values = vtkSmartPointer<vtkFloatArray>::New();
      values->SetName("VALUE");
      values->Allocate(numberOfDisplacementConstraints);
      for (vtkIdType i=0; i<numberOfDisplacementConstraints; i++)
        {
        reader.GetLineAndParse(3, tokens);
        ids->InsertNextValue(lexical_cast<vtkIdType>(tokens[0]) - 1);
        senses->InsertNextValue(lexical_cast<vtkIdType>(tokens[1]) - 1);
        values->InsertNextValue(lexical_cast<float>(tokens[2]));
        }
      vtkSmartPointer<vtkboneConstraint> constraint = vtkSmartPointer<vtkboneConstraint>::New();
      constraint->SetName("DisplacementConstraints");
      constraint->SetIndices(ids);
      constraint->SetConstraintType(vtkboneConstraint::DISPLACEMENT);
      constraint->SetConstraintAppliedTo(vtkboneConstraint::NODES);
      constraint->GetAttributes()->AddArray(senses);
      constraint->GetAttributes()->AddArray(values);
      output->GetConstraints()->AddItem(constraint);
      }

    // --------------------------------------------------------------
    // Read Node Sets

    reader.GetLineAndParse(1, tokens);
    vtkIdType numberOfNodeSets = lexical_cast<vtkIdType>(tokens[0]);
    for (int nodeSetId=0; nodeSetId<numberOfNodeSets; nodeSetId++)
      {
      reader.GetLineAndParse(1, tokens);
      vtkIdType numberOfIds = lexical_cast<vtkIdType>(tokens[0]);
      string selectionName = (format("NodeSet%d") % (nodeSetId+1)).str();
      vtkSmartPointer<vtkIdTypeArray> ids = vtkSmartPointer<vtkIdTypeArray>::New();
      ids->SetName (selectionName.c_str());
      ids->SetNumberOfTuples(numberOfIds);
      reader.ReadArray(numberOfIds,
                       reinterpret_cast<vtkIdType*>(ids->GetVoidPointer(0)));
      // Convert to 0-indexed from 1-indexed.
      for (vtkIdType i=0; i<numberOfIds; i++)
        {
        ids->SetValue(i,ids->GetValue(i)-1);
        }
      output->AddNodeSet(ids);
      vtkboneSolverParameters::POST_PROCESSING_NODE_SETS()->Append (
                                                info, selectionName.c_str());
      }

    // --------------------------------------------------------------
    // Read Element Sets

    reader.GetLineAndParse(1, tokens);
    vtkIdType numberOfElementSets = lexical_cast<vtkIdType>(tokens[0]);
    for (int elementSetId=0; elementSetId<numberOfElementSets; elementSetId++)
      {
      reader.GetLineAndParse(1, tokens);
      vtkIdType numberOfIds = lexical_cast<vtkIdType>(tokens[0]);
      string selectionName = (format("ElementSet%d") % (elementSetId+1)).str();
      vtkSmartPointer<vtkIdTypeArray> ids = vtkSmartPointer<vtkIdTypeArray>::New();
      ids->SetName (selectionName.c_str());
      ids->SetNumberOfTuples(numberOfIds);
      reader.ReadArray(numberOfIds,
                       reinterpret_cast<vtkIdType*>(ids->GetVoidPointer(0)));
      // Convert to 0-indexed from 1-indexed.
      for (vtkIdType i=0; i<numberOfIds; i++)
        {
        ids->SetValue(i,ids->GetValue(i)-1);
        }
      output->AddElementSet(ids);
      vtkboneSolverParameters::POST_PROCESSING_ELEMENT_SETS()->Append (
                                                info, selectionName.c_str());
      }

    }
  catch (boost::bad_lexical_cast&)
    {
    vtkErrorMacro(<<"Error reading file " << this->FileName
                  << ". Parse error line " << reader.line_count);
    return 0;
    }
  catch (std::exception& e)
    {
    vtkErrorMacro(<<"Error reading file " << this->FileName << ". " << e.what());
    return 0;
    }

  return 1;
}

ext: cxx


text:
/**
 ** Copyright (c) 2014 Illumina, Inc.
 **
 ** This file is part of Illumina's Enhanced Artificial Genome Engine (EAGLE),
 ** covered by the "BSD 2-Clause License" (see accompanying LICENSE file)
 **/

#include <iostream>
#include <string>
#include <vector>
#include <boost/assign.hpp>

using namespace std;
using boost::assign::list_of;

#include "Helpers.hh"

#include "RegistryName.hh"
#include "testBcl.hh"

using eagle::io::BclTile;

CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( TestBcl, registryName("Bcl"));

void TestBcl::setUp()
{
}

void TestBcl::tearDown()
{
}

void TestBcl::testBclTile()
{
    // Create a tile
    const unsigned int expectedReadCount = 5;
    const unsigned int clusterLength     = 101;
    const std::string filenameTemplate      = "";
    const std::string statsFilenameTemplate = "";
    const std::string filterFilename        = "";
    const std::string posFilename           = "";
    const std::string controlFilename       = "";

    BclTile tile( expectedReadCount, clusterLength, filenameTemplate, statsFilenameTemplate, filterFilename, posFilename, controlFilename, false );

    // Test that we can fill it
    const char bufCluster[102] = "12345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901";
    CPPUNIT_ASSERT_NO_THROW( tile.addClusterToRandomLocation( bufCluster ) );
    CPPUNIT_ASSERT_NO_THROW( tile.addClusterToRandomLocation( bufCluster ) );
    CPPUNIT_ASSERT_NO_THROW( tile.addClusterToRandomLocation( bufCluster ) );
    CPPUNIT_ASSERT_NO_THROW( tile.addClusterToRandomLocation( bufCluster ) );
    CPPUNIT_ASSERT_NO_THROW( tile.addClusterToRandomLocation( bufCluster ) ); // expectedReadCount == 5 => tile is now full

    // Test that we cannot overfill it
    CPPUNIT_ASSERT_THROW( tile.addClusterToRandomLocation( bufCluster ), eagle::common::OutOfLimitsException ); // too many clusters added to tile
}

max_stars_repo_name: sequencing/EAGLE

max_issues_repo_name: sequencing/EAGLE

max_forks_repo_name: sequencing/EAGLE


text:
/*
 * Copyright 2018,
 * Mirabel Joseph
 *
 * CNRS/AIST
 *
 */

#ifndef SOT_CORE_VARIADICOP_HH
#define SOT_CORE_VARIADICOP_HH

/* --------------------------------------------------------------------- */
/* --- INCLUDE --------------------------------------------------------- */
/* --------------------------------------------------------------------- */

/* Matrix */
#include <dynamic-graph/linear-algebra.h>

/* SOT */
#include <dynamic-graph/all-signals.h>
#include <dynamic-graph/entity.h>
#include <sot/core/flags.hh>
#include <sot/core/matrix-geometry.hh>
#include <sot/core/pool.hh>

/* STD */
#include <string>

#include <boost/function.hpp>

namespace dynamicgraph {
namespace sot {

/* --------------------------------------------------------------------- */
/* --- CLASS ----------------------------------------------------------- */
/* --------------------------------------------------------------------- */

template <typename Tin, typename Tout, typename Time>
class VariadicAbstract : public Entity {
public: /* --- CONSTRUCTION --- */
  static std::string getTypeInName(void);
  static std::string getTypeOutName(void);

  VariadicAbstract(const std::string &name, const std::string &className)
      : Entity(name), SOUT(className + "(" + name + ")::output(" +
                           getTypeOutName() + ")::sout"),
        baseSigname(className + "(" + name + ")::input(" + getTypeInName() +
                    ")::") {
    signalRegistration(SOUT);
  }

  virtual ~VariadicAbstract(void) {
    for (std::size_t i = 0; i < signalsIN.size(); ++i) {
      _removeSignal(i);
    }
  };

public: /* --- SIGNAL --- */
  typedef SignalPtr<Tin, int> signal_t;
  SignalTimeDependent<Tout, int> SOUT;

  std::size_t addSignal() {
    std::ostringstream oss;
    oss << "sin" << signalsIN.size();
    return addSignal(oss.str());
  }

  std::size_t addSignal(const std::string &name) {
    signal_t *sig = new signal_t(NULL, baseSigname + name);
    try {
      _declareSignal(sig);
      signalsIN.push_back(sig);
      // names.push_back (name);
      return signalsIN.size() - 1;
    } catch (const ExceptionAbstract &) {
      delete sig;
      throw;
    }
  }

  void removeSignal() {
    assert(signalsIN.size() > 0);
    _removeSignal(signalsIN.size() - 1);
    // names.pop_back();
    signalsIN.pop_back();
  }

  void setSignalNumber(const int &n) {
    assert(n >= 0);
    const std::size_t oldSize = signalsIN.size();
    for (std::size_t i = n; i < oldSize; ++i)
      _removeSignal(i);
    signalsIN.resize(n, NULL);
    // names.resize(n);

    for (std::size_t i = oldSize; i < (std::size_t)n; ++i) {
      assert(signalsIN[i] == NULL);
      std::ostringstream oss;
      oss << baseSigname << "sin" << i;
      // names[i] = oss.str();
      // signal_t* s = new signal_t (NULL,names[i]);
      signal_t *s = new signal_t(NULL, oss.str());
      signalsIN[i] = s;
      _declareSignal(s);
    }
    updateSignalNumber(n);
  }

  int getSignalNumber() const { return (int)signalsIN.size(); }

  signal_t *getSignalIn(int i) {
    if (i < 0 || i >= (int)signalsIN.size())
      throw std::out_of_range("Wrong signal index");
    return signalsIN[i];
  }

protected:
  std::vector<signal_t *> signalsIN;
  // Use signal->shortName instead
  // std::vector< std::string > names;

  virtual void updateSignalNumber(int n) { (void)n; };

private:
  void _removeSignal(const std::size_t i) {
    // signalDeregistration(names[i]);
    signalDeregistration(signalsIN[i]->shortName());
    SOUT.removeDependency(*signalsIN[i]);
    delete signalsIN[i];
  }
  void _declareSignal(signal_t *s) {
    signalRegistration(*s);
    SOUT.addDependency(*s);
  }
  const std::string baseSigname;
};

template <typename Operator>
class VariadicOp : public VariadicAbstract<typename Operator::Tin,
                                           typename Operator::Tout, int> {
  typedef typename Operator::Tin Tin;
  typedef typename Operator::Tout Tout;
  typedef VariadicOp<Operator> Self;

public: /* --- CONSTRUCTION --- */
  Operator op;

  typedef VariadicAbstract<Tin, Tout, int> Base;

  // static std::string getTypeInName ( void ) { return Operator::nameTypeIn ();
  // } static std::string getTypeOutName( void ) { return
  // Operator::nameTypeOut(); }
  static const std::string CLASS_NAME;
  virtual const std::string &getClassName() const { return CLASS_NAME; }
  std::string getDocString() const { return op.getDocString(); }

  VariadicOp(const std::string &name) : Base(name, CLASS_NAME) {
    this->SOUT.setFunction(boost::bind(&Self::computeOperation, this, _1, _2));
    op.initialize(this, this->commandMap);
  }

  virtual ~VariadicOp(void){};

protected:
  Tout &computeOperation(Tout &res, int time) {
    std::vector<const Tin *> in(this->signalsIN.size());
    for (std::size_t i = 0; i < this->signalsIN.size(); ++i) {
      const Tin &x = this->signalsIN[i]->access(time);
      in[i] = &x;
    }
    op(in, res);
    return res;
  }

  inline void updateSignalNumber(int n) { op.updateSignalNumber(n); }
};
} // namespace sot
} // namespace dynamicgraph

#endif // #ifndef SOT_CORE_VARIADICOP_HH

ext: hh

max_stars_repo_name: machines-in-motion/sot-core

max_issues_repo_name: machines-in-motion/sot-core

max_forks_repo_name: machines-in-motion/sot-core


text:
#include "businessobjectlist.h"
#include "bitcoinunits.h"
#include "ui_businessobjectlist.h"
#include "secdialog.h"
#include "ui_secdialog.h"
#include "walletmodel.h"
#include <QPainter>
#include <QTableWidget>
#include <QGridLayout>
#include <QUrl>
#include <univalue.h>
#include <boost/algorithm/string/case_conv.hpp> // for to_lower()

QString ToQstring(std::string s);
std::string FromQStringW(QString qs);
std::string RoundToString(double d, int place);
std::string GetBusinessObjectList(std::string sType, std::string sFields);
UniValue GetBusinessObjectByFieldValue(std::string sType, std::string sFieldName, std::string sSearchValue);
UniValue GetBusinessObject(std::string sType, std::string sPrimaryKey, std::string& sError);
double cdbl(std::string s, int place);
std::string ObjectType = "";
extern int GetUrlColumn(std::string sTarget);

QStringList BusinessObjectList::GetHeaders()
{
	QStringList pHeaders;
	UniValue aBO = GetBusinessObjectByFieldValue("object", "object_name", ObjectType);
	std::string sFields = "id," + aBO["fields"].getValStr();
	std::vector<std::string> vFields = Split(sFields.c_str(),",");
	for (int i = 0; i < (int)vFields.size(); i++)
	{
		std::string sFieldName = vFields[i];
		pHeaders << ToQstring(sFieldName);
	}
	return pHeaders;
}

BusinessObjectList::BusinessObjectList(const PlatformStyle *platformStyle, QWidget *parent) : ui(new Ui::BusinessObjectList)
{
    ui->setupUi(this);
}

BusinessObjectList::~BusinessObjectList()
{
    delete ui;
}

void BusinessObjectList::setModel(WalletModel *model)
{
    this->model = model;
}

void BusinessObjectList::UpdateObject(std::string objType)
{
	ObjectType = objType;

	LogPrintf(" objType changed to ObjType %s ",ObjectType.c_str());

	UniValue aBO = GetBusinessObjectByFieldValue("object", "object_name", ObjectType);
	std::string sFields = aBO["fields"].getValStr();
    QString pString = ToQstring(GetBusinessObjectList(ObjectType, sFields));

    QStringList pHeaders = GetHeaders();
    this->createUI(pHeaders, pString);
}

void BusinessObjectList::createUI(const QStringList &headers, const QString &pStr)
{
    ui->tableWidget->setShowGrid(true);
    ui->tableWidget->setSelectionMode(QAbstractItemView::SingleSelection);
    ui->tableWidget->setSelectionBehavior(QAbstractItemView::SelectRows);
    ui->tableWidget->horizontalHeader()->setStretchLastSection(true);

    QVector<QVector<QString> > pMatrix;
	if (pStr == "") return;

    pMatrix = SplitData(pStr);
	int iAmountCol = GetUrlColumn("Amount");
	int iFooterRow = (iAmountCol > -1) ? 1 : 0;
    int rows = pMatrix.size();
    ui->tableWidget->setRowCount(rows + iFooterRow);
    int cols = pMatrix[0].size() - 1;
    ui->tableWidget->setColumnCount(cols);
    ui->tableWidget->setHorizontalHeaderLabels(headers);
    QString s;
	double dGrandTotal = 0;
    for(int i=0; i < rows; i++)
	{
        for(int j=0; j < cols; j++)
		{
            ui->tableWidget->setItem(i,j, new QTableWidgetItem(pMatrix[i][j]));
		}
		if (iAmountCol > -1)
		{
			dGrandTotal += cdbl(FromQStringW(pMatrix[i][iAmountCol]), 2);
		}
	}
	if (iFooterRow > 0)
	{
		ui->tableWidget->setItem(rows, 0, new QTableWidgetItem("Grand Total:"));
		ui->tableWidget->setItem(rows, iAmountCol, new QTableWidgetItem(ToQstring(RoundToString(dGrandTotal, 2))));
	}

    ui->tableWidget->setEditTriggers(QAbstractItemView::NoEditTriggers);
    ui->tableWidget->resizeRowsToContents();
    ui->tableWidget->resizeColumnsToContents();
    ui->tableWidget->setContextMenuPolicy(Qt::CustomContextMenu);

	// Column widths should be set
	for (int j=0; j < cols; j++)
	{
		ui->tableWidget->setColumnWidth(j, 128);
	}

    // Connect SLOT to context menu
    connect(ui->tableWidget, SIGNAL(customContextMenuRequested(QPoint)), this, SLOT(slotCustomMenuRequested(QPoint)));
}

void BusinessObjectList::slotCustomMenuRequested(QPoint pos)
{
    /* Create an object context menu */
    QMenu * menu = new QMenu(this);
    //  Create, Connect and Set the actions to the menu
    menu->addAction(tr("Navigate To"), this, SLOT(slotNavigateTo()));
	menu->addAction(tr("List"), this, SLOT(slotList()));
	menu->addAction(tr("View"), this, SLOT(slotView()));

    menu->popup(ui->tableWidget->viewport()->mapToGlobal(pos));
}

std::string GetHtmlForm(std::string PK)
{
	UniValue aBO = GetBusinessObjectByFieldValue("object", "object_name", ObjectType);
	std::string sFields = aBO["fields"].getValStr();
	std::vector<std::string> vFields = Split(sFields.c_str(), ",");
	std::vector<std::string> vPK = Split(PK.c_str(), "-");
	std::string sError = "";
	if (vPK.size() > 2)
	{
		aBO = GetBusinessObject(vPK[0], vPK[1], sError);
		std::string sHTML = "
";
		for (int i = 0; i < (int)vFields.size(); i++)
		{
			std::string sFieldName = vFields[i];
			std::string sFV = aBO[sFieldName].getValStr();
			std::string sRow = "
" + sFieldName + ":	" + sFV + "
";
			sHTML += sRow;
		}
		sHTML += "
";
		return sHTML;
	}
	return "Business object not found";
}

int GetUrlColumn(std::string sTarget)
{
	boost::to_upper(sTarget);
	UniValue aBO = GetBusinessObjectByFieldValue("object", "object_name", ObjectType);
	std::string sFields = "id," + aBO["fields"].getValStr();
	std::vector<std::string> vFields = Split(sFields.c_str(), ",");
	for (int i = 0; i < (int)vFields.size(); i++)
	{
		std::string sFieldName = vFields[i];
		boost::to_upper(sFieldName);
		if (sFieldName == sTarget) return i;
	}
	return -1;
}

void BusinessObjectList::slotList()
{
    int row = ui->tableWidget->selectionModel()->currentIndex().row();
    if(row >= 0)
    {
		// Navigate to the List of the Object
        std::string sID = FromQStringW(ui->tableWidget->item(row, 0)->text()); // PK-2PK-IPFS Hash of business object
		int iCol = GetUrlColumn("object_name");
		if (iCol > -1)
		{
			std::string sTarget = FromQStringW(ui->tableWidget->item(row, iCol)->text());
			UpdateObject(sTarget);
		}
    }
}

void BusinessObjectList::slotView()
{
    int row = ui->tableWidget->selectionModel()->currentIndex().row();
    if(row >= 0)
    {
        QMessageBox msgBox;
        std::string id = FromQStringW(ui->tableWidget->item(row, 0)->text()); // PK-2PK-IPFS Hash of business object
		msgBox.setWindowTitle(ToQstring(id));
	    msgBox.setText(ToQstring(GetHtmlForm(id)));
        msgBox.setStandardButtons(QMessageBox::Ok);
        msgBox.setDefaultButton(QMessageBox::Ok);
		msgBox.exec();
    }
}

void BusinessObjectList::slotNavigateTo()
{
    int row = ui->tableWidget->selectionModel()->currentIndex().row();
    if(row >= 0)
    {
		// Open the URL
        QMessageBox msgBox;
		int iURLCol = GetUrlColumn("URL");
		if (iURLCol > -1)
		{
			QString Url = ui->tableWidget->item(row, iURLCol)->text();
			QUrl pUrl(Url);
			QDesktopServices::openUrl(pUrl);
		}
    }
}

QVector<QVector<QString> > BusinessObjectList::SplitData(const QString &pStr)
{
	QStringList proposals = pStr.split(QRegExp("<object>"),QString::SkipEmptyParts);
    int nProposals = proposals.size();
    QVector<QVector<QString> > proposalMatrix;
    for (int i=0; i < nProposals; i++)
    {
        proposalMatrix.append(QVector<QString>());
        QStringList proposalDetail = proposals[i].split(QRegExp("<col>"));
        int detailSize = proposalDetail.size();
        for (int j = 0; j < detailSize; j++)
		{
			QString sData = proposalDetail[j];
			/*  Reserved for BitcoinUnits
				sData = BitcoinUnits::format(2, cdbl(FromQStringW(sData), 2) * 100, false, BitcoinUnits::separatorAlways);
			*/
			proposalMatrix[i].append(sData);
		}
    }
	return proposalMatrix;
}