wdk.h

#ifndef WDK_H
#define WDK_H


#include <stdexcept> 
#include <iostream>
#include <sstream>
#include <vector>
#include <list>
#include <map>
#include <algorithm>
#include <numeric>
#include <cassert>
#include <cmath>
#include <cctype>

using namespace std;


typedef float RealType;
typedef long int IndexConcreteType;
typedef int ValueConcreteType;


enum KernelType {DOT_PRODUCT_TYPE=0,HISTOGRAM_INTERSECTION_TYPE};
class UserKernelParametersClass
{
 public:
  UserKernelParametersClass():
    mKernelType(HISTOGRAM_INTERSECTION_TYPE),
    mExpGamma(0),
    mPolyDegree(0),
    mPolyOffset(1),
    mPolyCoefficient(1),
    mTotalNormalization(true),
    mProductNormalization(true),
    mAttributeNormalization(true)
    {}

  void Input(const char* aString) throw(exception)
    {
      vector<string> options;
      stringstream ss;
      ss<<(string)aString;
      while (!ss.eof())
        {
          string data;
          ss>>data;
          if (data!="") options.push_back(data);
        }
      for (vector<string>::iterator it=options.begin();it!=options.end();++it)
        {
          if ((*it)=="-Kd") mPolyDegree=atof((*(++it)).c_str());
          else if ((*it)=="-Kr") mPolyOffset=atof((*(++it)).c_str());
          else if ((*it)=="-Ks") mPolyCoefficient=atof((*(++it)).c_str());
          else if ((*it)=="-Kg") mExpGamma=atof((*(++it)).c_str());
          else if ((*it)=="-KDot") mKernelType=DOT_PRODUCT_TYPE;
          else if ((*it)=="-KHistInt") mKernelType=HISTOGRAM_INTERSECTION_TYPE;
          else if ((*it)=="-KNoTotNorm") mTotalNormalization=false;
          else if ((*it)=="-KNoProdNorm") mProductNormalization=false;
          else if ((*it)=="-KNoAttrNorm") mAttributeNormalization=false;
          else throw(invalid_argument("Argument:"+(*it)+" is not valid"));
        }
      if(mExpGamma*mPolyDegree)//Mutual exclusion between polinomial and exponential kernel
        throw(invalid_argument("Error: exponential and polinomial cannot be set active at the same time"));
      ComputeInversePolinomialNormalizationCoefficient();
    }
  float ComputeInversePolinomialNormalizationCoefficient()
    {
      mInversePolinomialNormalizationCoefficient=1/pow(mPolyCoefficient+mPolyOffset,mPolyDegree);
    }
  KernelType mKernelType; 
  float mExpGamma; 
  float mPolyDegree; 
  float mPolyOffset; 
  float mPolyCoefficient; 
  bool mTotalNormalization; 
  bool mProductNormalization; 
  bool mAttributeNormalization; 
  float mInversePolinomialNormalizationCoefficient; 
};

static UserKernelParametersClass user_kernel_param; 
template <class OutType, class InType>
OutType stream_cast(const InType & t)
{
 stringstream ss;
 ss << t; // first insert value to stream
 OutType result; // value will be converted to OutType
 ss >> result; // write value to result
 return result;
}

pair<string,string> Split(const string& aInput,char sep=':') throw(exception)
{
  string::size_type pos=aInput.find(sep);
  if (pos==string::npos) throw domain_error("Non tokenizable element:"+aInput);
  return make_pair(aInput.substr(0,pos),aInput.substr(pos+1,aInput.size()-(pos+1)));
}

bool IsNumber(const string& aInput)
{
  for (unsigned i=0;i<aInput.size();++i)
    if (!isdigit(aInput[i])) return false;
  return true;
}

template<typename IndexType, typename ValueType>
  struct LessPair{
    inline bool operator()(const pair<IndexType,ValueType>& A,const pair<IndexType,ValueType>& B)
    {
      return A.first<B.first;
    }
  };

template<typename IndexType, typename ValueType>
  struct EqualPair{
    EqualPair(const IndexType& aItem){mItem=aItem;}
    inline bool operator()(const pair<IndexType,ValueType>& A)
    {
      return A.first==mItem;
    }
    IndexType mItem;
  };

//-----------------------------------------------------------------------------------------------

//Note:friend declaration is a usage of function not a proper declaration
//We first declare operator<<, then qualify it as friend, then define it and define its specialization
template<typename IndexType, typename ValueType> class SparseVectorClass;
template<typename IndexType, typename ValueType> ostream& operator<< (ostream& ss, const SparseVectorClass<IndexType,ValueType>& data);


template<typename IndexType, typename ValueType>
class SparseVectorClass
{
  friend ostream& operator<< <IndexType,ValueType>(ostream& ss, const SparseVectorClass<IndexType,ValueType>& data);
 public:
  SparseVectorClass(){}

  void Insert(const IndexType& aIndex, const ValueType& aValue)
    {
      mData.push_back(make_pair(aIndex,aValue));
    }

  void Insert(const pair<IndexType,ValueType>& aItem)
    {
      mData.push_back(aItem);
    }

  int Size()const{return mData.size();}
 public:
  list<pair<IndexType,ValueType> > mData; 
};

//Note: generally we want to characterize the semantic of the sparse vector by prepending a description before its data, i.e. we do not want to output the index but only the value
template<typename IndexType, typename ValueType>
ostream& operator<< (ostream& out, const SparseVectorClass<IndexType,ValueType>& data)
{
  typename list<pair<IndexType,ValueType> >::const_iterator it=data.mData.begin();
  for (;it!=data.mData.end();++it)
    out<<" "<<it->second;
  return out;
}

//Note: in the special case of an <int:int> sparse vector the output includes both the index and the value
template<>
ostream& operator<< <IndexConcreteType,ValueConcreteType>(ostream& out, const SparseVectorClass<IndexConcreteType,ValueConcreteType>& data)
    {
      list<pair<IndexConcreteType,ValueConcreteType> >::const_iterator it=data.mData.begin();
      for (;it!=data.mData.end();++it)
        out<<" "<<it->first<<":"<<it->second;
      return out;
    }


//-----------------------------------------------------------------------------------------------

template<typename IndexType, typename ValueType>
class DotSparseVectorClass:public SparseVectorClass<IndexType,ValueType>
{
 public:
  DotSparseVectorClass(){}

  void Sort()
    {
      this->mData.sort(LessPair<IndexType,ValueType>());
    }

  RealType operator*(const DotSparseVectorClass<IndexType,ValueType>& aInstance)const
    {
      RealType result=0;
      typename list<pair<IndexType,ValueType> >::const_iterator it,jt;
      it=this->mData.begin();
      jt=aInstance.mData.begin();
      while(it!=this->mData.end() && jt!=aInstance.mData.end())
        {
          if (it->first==jt->first)
            {
              result+=it->second * jt->second;
              ++it;
              ++jt;
            }
          else if (it->first<jt->first) ++it;
          else //(it->first==jt->first)
            ++jt;
        }
      return result;
    }
};


//-----------------------------------------------------------------------------------------------

template<typename IndexType, typename ValueType>
class HistogramIntersectionOptimizedSparseVectorClass:public SparseVectorClass<IndexType,ValueType>
{
  friend ostream& operator<<(ostream& out, const HistogramIntersectionOptimizedSparseVectorClass& data)
    {
      out<<" dim:"<<data.mData.size();
      for (unsigned i=0;i<data.mData.size();++i)
        {
          out<<" bin:"<<data.mData[i].first;
          out<<" dim:"<<data.mData[i].second.size();
          for (unsigned j=0;j<data.mData[i].second.size();++j)
            out<<" "<<data.mData[i].second[j];
        }
      return out;
    }
  
 public:
  void Insert(const IndexType& aIndex, const ValueType& aValue)
    {
      //find index
      EqualPair<IndexType,vector<ValueType> > equal(aIndex);
      typename vector<pair<IndexType,vector<ValueType> > >::iterator it=find_if(mData.begin(),mData.end(),equal);
      //if it exists
      if (it!=mData.end())
        {
          //add value to corresponding vector
          it->second.push_back(aValue);
        }
      //else add new vector containing value
      else
        {
          vector<ValueType> tmp_vec;
          tmp_vec.push_back(aValue);
          mData.push_back(make_pair(aIndex,tmp_vec));
        }
    }
  void Insert(const pair<IndexType,ValueType>& aItem)
    {
      Insert(aItem.first,aItem.second);
    }
  void Sort()
    {
      sort(mData.begin(),mData.end(),LessPair<IndexType,vector<ValueType> >());
      for (typename vector<pair<IndexType,vector<ValueType> > >::iterator it=mData.begin();it!=mData.end();++it)
        {
          sort(it->second.begin(),it->second.end());
        }
    }
  RealType operator*(const HistogramIntersectionOptimizedSparseVectorClass<IndexType,ValueType>& aInstance)const
    {
      typename vector<pair<IndexType,vector<ValueType> > >::const_iterator it=mData.begin();
      typename vector<pair<IndexType,vector<ValueType> > >::const_iterator jt=aInstance.mData.begin();

      ValueType result=0;
      while(it!=mData.end() && jt!=aInstance.mData.end())
        {
          if (it->first==jt->first)
            {
              result+=OptimizedHistogramIntersection(it->second,jt->second);
              it++;
              jt++;
            }
          else if (it->first < jt->first) it++;
          else //it->first > jt->first
            jt++;
        }
      return result;
    }
  ValueType OptimizedHistogramIntersectionBaseVersion(const vector<ValueType>& aVecA, const vector<ValueType>& aVecB)const
    {
      ValueType result=0;
      for (unsigned i=0;i<aVecA.size();++i)
        {
          unsigned j;
          for (j=0;j<aVecB.size() && aVecB[j]<=aVecA[i];++j)
            result+=aVecB[j];
          result+=aVecA[i]*(aVecB.size()-j);
        }
      return result;
    }
  ValueType OptimizedHistogramIntersection(const vector<ValueType>& aVecA, const vector<ValueType>& aVecB)const
    {
      ValueType result=0;
      ValueType cumulative=0;
      unsigned i=0,j=0;
      while (i<aVecA.size() && j<aVecB.size())
        {
          if (aVecA[i]<aVecB[j])
            {
              result+=aVecA[i]*(aVecB.size()-j)+cumulative;
              i++;
            }
          else if (aVecA[i]>=aVecB[j])
            {
              cumulative+=aVecB[j];
              j++;
            }
        }
      //if there are still elements in aVecA (i.e. i<aVewcA.size()) but not in aVecB then it means that they are all greater and therefore:
      result+=(aVecA.size()-i)*cumulative;
      return result;
    }
 public:
  vector<pair<IndexType,vector<ValueType> > > mData;
};

//-----------------------------------------------------------------------------------------------

class AttributeClass
{
  friend ostream& operator<<(ostream& out, const AttributeClass& data)
    {
      out<<"attribute:"<<data.mAttributeType<<" ";
      out<<"dim:"<<data.mValue.Size();
      out<<data.mValue;
      out<<data.mOptimizedValue;
      return out;
    }
  friend istream& operator>>(istream& in, AttributeClass& data) throw(exception)
    {
      string token;
      in>>token;pair<string,string> key_value_token=Split(token);
      if(key_value_token.first!="attribute") throw domain_error("Unexpected element:"+key_value_token.first);
        
      data.mAttributeType=stream_cast<IndexConcreteType>(key_value_token.second);

      in>>token;key_value_token=Split(token);
      if (key_value_token.first!="dim") throw domain_error("Unexpected element:"+key_value_token.first);
      if (!IsNumber(key_value_token.second)) throw domain_error("Expected a number:"+key_value_token.second);
      int dim=stream_cast<int>(key_value_token.second);

      for (int i=0;i<dim;i++)
        {
          IndexConcreteType index;
          ValueConcreteType value;
          in>>token;pair<string,string> key_value_token=Split(token);
          if (!IsNumber(key_value_token.first)) throw domain_error("Expected a number:"+key_value_token.first);
          if (!IsNumber(key_value_token.second)) throw domain_error("Expected a number:"+key_value_token.second);
          index=stream_cast<IndexConcreteType>(key_value_token.first);
          value=stream_cast<ValueConcreteType>(key_value_token.second);
          data.mValue.Insert(index,value);
          data.mOptimizedValue.Insert(index,value);
        }
      return in;
    }
 public:
  AttributeClass():
    mAttributeType(0),
    mNorm(0),
    mInverseSquaredNorm(0){}
  void Merge(AttributeClass& aAttribute) throw(exception)
    {
      switch(user_kernel_param.mKernelType)
        {
        case DOT_PRODUCT_TYPE:
          MergeDotProduct(aAttribute);
          break;
        case HISTOGRAM_INTERSECTION_TYPE:
          MergeHistogramIntersection(aAttribute);
          break;
        default: throw(invalid_argument("ERROR:Unknown kernel"));
        }
    }
  RealType operator*(const AttributeClass& aInstance)const throw(exception)
    {
      if (user_kernel_param.mAttributeNormalization)
        {
          //Note: lazy normalization: this is needed since the
          //optimization procedure alters the attribute data types, so
          //that no normalization can be computed just after parsing
          //the data in input
          if (mNorm==0) ComputeNormalizationFactor();
          if (aInstance.mNorm==0) aInstance.ComputeNormalizationFactor();       
          switch(user_kernel_param.mKernelType)
            {
            case DOT_PRODUCT_TYPE:
              return (mValue*aInstance.mValue)*mInverseSquaredNorm*aInstance.mInverseSquaredNorm;
            case HISTOGRAM_INTERSECTION_TYPE:
              return (mOptimizedValue*aInstance.mOptimizedValue)*mInverseSquaredNorm*aInstance.mInverseSquaredNorm;
            default: throw(invalid_argument("ERROR:Unknown kernel"));
            }
        }
      else
        {
          switch(user_kernel_param.mKernelType)
            {
            case DOT_PRODUCT_TYPE:
              return (mValue*aInstance.mValue);
            case HISTOGRAM_INTERSECTION_TYPE:
              return (mOptimizedValue*aInstance.mOptimizedValue);
            default: throw(invalid_argument("ERROR:Unknown kernel"));
            }
        }
    }
 protected:
  void MergeDotProduct(AttributeClass& aAttribute)
    {
      //invoke add on matching histogram objects
      list<pair<IndexConcreteType,ValueConcreteType> >::iterator it,jt;
      it=mValue.mData.begin();
      jt=aAttribute.mValue.mData.begin();

      while(it!=mValue.mData.end() && jt!=aAttribute.mValue.mData.end())
        {
          if (it->first==jt->first)
            {
              it->second+=jt->second;
              ++it;
              ++jt;
            }
          else if (it->first < jt->first) ++it;
          //case when there is an element in aAttribute that is not present in current object
          else 
            {
              mValue.Insert(*jt);
              ++jt;
            }
        }
      //extremal cases: copy remaining aAttribute.mValue elements
      while (jt!=aAttribute.mValue.mData.end())
        {
          mValue.Insert(*jt);
          ++jt;
        }
      //sort for later use
      mValue.Sort();
    }
  void MergeHistogramIntersection(AttributeClass& aAttribute)
    {
      for (list<pair<IndexConcreteType,ValueConcreteType> >::iterator it=aAttribute.mValue.mData.begin();it!=aAttribute.mValue.mData.end();++it)
        mOptimizedValue.Insert(it->first,it->second);
     
      //sort for later use
      mOptimizedValue.Sort();
    }
  void ComputeNormalizationFactor()const throw(exception)
    {
      switch(user_kernel_param.mKernelType)
        {
        case DOT_PRODUCT_TYPE:
          mNorm=(mValue*mValue);
          mInverseSquaredNorm=1/sqrt(mNorm);
          break;
        case HISTOGRAM_INTERSECTION_TYPE:
          mNorm=(mOptimizedValue*mOptimizedValue);
          mInverseSquaredNorm=1/sqrt(mNorm);
          break;
        default:
          throw(invalid_argument("ERROR:Unknown kernel"));
        }
    }

 public:
  IndexConcreteType mAttributeType;
  mutable RealType mNorm;
  mutable RealType mInverseSquaredNorm;
  DotSparseVectorClass<IndexConcreteType,ValueConcreteType> mValue; 
  HistogramIntersectionOptimizedSparseVectorClass<IndexConcreteType,ValueConcreteType> mOptimizedValue;
};

//-----------------------------------------------------------------------------------------------
class PartClass
{
 public:
  friend ostream& operator<<(ostream& out, const PartClass& data)
    {
      out<<"part:"<<data.mPartType<<" ";
      out<<"dim:"<<data.mAttr.Size();
      out<<data.mAttr;
      return out;
    }
  friend istream& operator>>(istream& in, PartClass& data) throw(exception)
    {
      string token;
      in>>token;pair<string,string> key_value_token=Split(token);
      if (key_value_token.first!="part") throw domain_error("Unexpected element:"+key_value_token.first);
      data.mPartType=stream_cast<IndexConcreteType>(key_value_token.second);
      in>>token;key_value_token=Split(token);
      if (key_value_token.first!="dim") throw domain_error("Unexpected element:"+key_value_token.first);
      if (!IsNumber(key_value_token.second)) throw domain_error("Expected a number:"+key_value_token.second);
      int dim=stream_cast<int>(key_value_token.second);
      for (int i=0;i<dim;i++)
        {
          AttributeClass attribute;
          in>>attribute;
          data.mAttr.Insert(attribute.mAttributeType,attribute);
        }
      data.mAttr.Sort();
      return in;
    }
 public:
  PartClass(){}
  void Merge(PartClass& aPart)
    {
      //invoke add on matching attribute objects
      list<pair<IndexConcreteType,AttributeClass> >::iterator it,jt;
      it=mAttr.mData.begin();
      jt=aPart.mAttr.mData.begin();

      while(it!=mAttr.mData.end() && jt!=aPart.mAttr.mData.end())
        {
          if (it->first==jt->first)
            {
              it->second.Merge(jt->second);
              ++it;
              ++jt;
            }
          else if (it->first < jt->first) ++it;
          //case when there is an element in aPart that is not present in current object
          else 
            {
              mAttr.Insert(*jt);
              ++jt;
            }
        }
      //extremal cases: copy remaining aPart.mAttr elements
      while (jt!=aPart.mAttr.mData.end())
        {
          mAttr.Insert(*jt);
          ++jt;
        }
      //sort for later use
      mAttr.Sort();
    }
  RealType operator*(const PartClass& aInstance)const
    {
      return mAttr*aInstance.mAttr;
    }
 public:
  IndexConcreteType mPartType;
  DotSparseVectorClass<IndexConcreteType,AttributeClass> mAttr; 
};

//-----------------------------------------------------------------------------------------------

class WDKDataClass
{
 public: 
  friend ostream& operator<<(ostream& out, const WDKDataClass& data)
    {
      out<<"dim:"<<data.mPart.Size();
      out<<data.mPart;
      return out;
    }
  friend istream& operator>>(istream& in, WDKDataClass& data) throw(exception)
    {
      string token;
      int dim;
      in>>token;pair<string,string> key_value_token=Split(token);
      if (key_value_token.first!="dim") throw domain_error("Unexpected element:"+key_value_token.first);
      if (!IsNumber(key_value_token.second)) throw domain_error("Expected a number:"+key_value_token.second);
      dim=stream_cast<int>(key_value_token.second);
      for (int i=0;i<dim;i++)
        {
          PartClass part;
          in>>part;
          data.mPart.Insert(part.mPartType,part);
        }
      data.Optimize();
      return in;
    }
 public: 
  WDKDataClass():mNorm(0),mInverseSquaredNorm(0){}
  WDKDataClass(const char* aS):mNorm(1),mInverseSquaredNorm(1)
    {
      mSerialized=(string)(aS);
      stringstream ss(aS);
      ss>>(*this);
      ComputeNormalizationFactor();
    }
  string Serialize()
    {
      return mSerialized;
    }
  RealType operator*(const WDKDataClass& aInstance)const
    {
      if (mNorm==0 || aInstance.mNorm==0) return 0;//Case of at least one null element
      if (user_kernel_param.mProductNormalization)
        return Product(aInstance)*mInverseSquaredNorm*aInstance.mInverseSquaredNorm;
      else
        return Product(aInstance);
    }
 protected:
  void Optimize()
    {
      //sort by part type (selector)
      mPart.Sort();
      //optimize attributes: 
      //in dot product case compress attributes 
      //in histogram intersection kernel compute attribute bin vectors 
      for (list<pair<IndexConcreteType,PartClass> >::iterator it=mPart.mData.begin();it!=mPart.mData.end();)
        {
          list<pair<IndexConcreteType,PartClass> >::iterator jt=it;
          jt++;
          if (jt!=mPart.mData.end())
            {
              //if part types (selectors) are equal then optimize histograms
              if (it->first==jt->first)
                {
                  it->second.Merge(jt->second);
                  mPart.mData.erase(jt);
                }
              //else advance index
              else ++it;
            }
          else ++it;//we are one element before the end of data structure => just advance to next element and terminate
        }
    }
  void ComputeNormalizationFactor()
    {
      mNorm=Product(*this);
      mInverseSquaredNorm=1/sqrt(mNorm);
    }
 protected:
  RealType Product(const WDKDataClass& aInstance)const
    {
      return mPart*aInstance.mPart;
    }
 public:
  RealType mNorm;
  RealType mInverseSquaredNorm;
 protected:
  DotSparseVectorClass<IndexConcreteType,PartClass> mPart; 
  string mSerialized;
};

class UserKernelClass
{
 public:
  RealType K(const WDKDataClass& aX,const WDKDataClass& aZ, const UserKernelParametersClass& aParam)
    {
      try
        {
          if (aParam.mTotalNormalization)
            {
              if (aParam.mPolyDegree!=0) 
                {
                  if (aParam.mProductNormalization) return k(aX,aZ,aParam)*aParam.mInversePolinomialNormalizationCoefficient;
                  else return k(aX,aZ,aParam)/sqrt(k(aX,aX,aParam)*k(aZ,aZ,aParam));
                }
              else //if (aParam.mExpGamma) or linear case 
                return k(aX,aZ,aParam); //Note:normalization is implicit in exponential case
            }
          else return k(aX,aZ,aParam);
        }
      catch(exception e)
        {
          cerr<<e.what()<<endl;
          exit(1);
        }
    }
  RealType k(const WDKDataClass& aX,const WDKDataClass& aZ, const UserKernelParametersClass& aParam)
    {
      if (aParam.mPolyDegree!=0)
        {
          return pow(aParam.mPolyCoefficient*(aX*aZ)+aParam.mPolyOffset,aParam.mPolyDegree);
        }
      else if (aParam.mExpGamma)
        {
          if (aParam.mProductNormalization) //if data product is normalized then squared norm of data is 1
            return exp(-aParam.mExpGamma*(2 -2*(aX*aZ)));
          else
            return exp(-aParam.mExpGamma*(aX.mNorm -2*(aX*aZ) + aZ.mNorm));
        }
      else return (aX*aZ);
    }
};

#endif

---