//----------------------------------------------------------------------------
// Copyright (C) 2013-2015 Fabrice HARROUET (ENIB)
//
// Permission to use, copy, modify, distribute and sell this software
// and its documentation for any purpose is hereby granted without fee,
// provided that the above copyright notice appear in all copies and
// that both that copyright notice and this permission notice appear
// in supporting documentation.
// The author makes no representations about the suitability of this
// software for any purpose.
// It is provided "as is" without express or implied warranty.
//----------------------------------------------------------------------------

// nb: DTW applies to time series
//     using it on these random data is a total nonsense
// #define USE_DTW

using System;
using System.Diagnostics;

using ParallelRunner;

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

public class Work
{
public
Work(int inputWidth,
     int tableLength,
     int threadCount)
  {
  _input=new float[inputWidth];
  _table=new float[tableLength][];
  _threadCount=threadCount;
  _rand=new Random();
  for(int j=0;j<_table.Length;++j)
    {
    _table[j]=new float[_input.Length];
    for(int i=0;i<_input.Length;++i)
      {
      _table[j][i]=(float)_rand.NextDouble();
      }
    }
  _threadEuclidian=new Result[threadCount];
  _euclidian=new Result();
  _threadCosine=new Result[threadCount];
  _cosine=new Result();
  _threadTanimoto=new Result[threadCount];
  _tanimoto=new Result();
  _threadSorensen=new Result[threadCount];
  _sorensen=new Result();
  _threadJaccard=new Result[threadCount];
  _jaccard=new Result();
#if USE_DTW
  _threadDtw=new Result[threadCount];
  _dtw=new Result();
#endif
  for(int i=0;i<_threadCount;++i)
    {
    _threadEuclidian[i]=new Result();
    _threadCosine[i]=new Result();
    _threadTanimoto[i]=new Result();
    _threadSorensen[i]=new Result();
    _threadJaccard[i]=new Result();
#if USE_DTW
    _threadDtw[i]=new Result();
#endif
    }
  }

public class Result
  {
  public const int SIZE=8;
  public int[] line;
  public float[] similarity;
  public float min;
  public float max;
  public
  Result()
    {
    line=new int[SIZE];
    similarity=new float[SIZE];
    reset();
    }
  public
  void
  reset()
    {
    for(int i=0;i<SIZE;++i)
      {
      line[i]=-1;
      similarity[i]=float.MinValue;
      }
    min=float.MaxValue;
    max=float.MinValue;
    }
  public
  void
  insert(int line,
         float similarity)
    {
    for(int i=0;i<SIZE;++i)
      {
      if(similarity>this.similarity[i])
        {
        insertAt(i,line,similarity);
        break;
        }
      }
    }
  public
  void
  insertAt(int index,
           int line,
           float similarity)
    {
    for(int i=SIZE-1;i>index;--i)
      {
      this.line[i]=this.line[i-1];
      this.similarity[i]=this.similarity[i-1];
      }
    this.line[index]=line;
    this.similarity[index]=similarity;
    min=(float)Math.Min(min,similarity);
    max=(float)Math.Max(max,similarity);
    }
  public
  void
  merge(Result other)
    {
    int firstPos=0;
    for(int i=0;i<SIZE;++i)
      {
      int pos;
      for(pos=firstPos;pos<SIZE;++pos)
        {
        if(other.similarity[i]>similarity[pos])
          {
          insertAt(pos,other.line[i],other.similarity[i]);
          firstPos=pos+1;
          break;
          }
        }
      if(pos==SIZE) break; // not inserted so next ones will not
      }
    min=(float)Math.Min(min,other.min);
    max=(float)Math.Max(max,other.max);
    }
  public override
  string
  ToString()
    {
    string s=String.Format("[{0:F3}%;{1:F3}%]",100.0f*min,100.0f*max);
    for(int i=0;i<SIZE;++i)
      {
      s+=String.Format("({0}:{1:F3}%)",line[i],100.0f*similarity[i]);
      }
    return s;
    }
  }

public
float // iteration per seconds
run(float duration,
    bool spinning)
  {
  int iterationCount=0;
  float elapsed=0.0f;
  // ensure cleanup is done through Dispose() 
  using(Runner runner=Runner.create(_threadCount,spinning))
    {
    Stopwatch sw=Stopwatch.StartNew();
    while(elapsed<duration)
      {
      for(int i=0;i<_input.Length;++i)
        {
        _input[i]=(float)_rand.NextDouble();
        }
      runner.run(_similarityTask);
      _euclidian.reset();
      _cosine.reset();
      _tanimoto.reset();
      _sorensen.reset();
      _jaccard.reset();
#if USE_DTW
      _dtw.reset();
#endif
      for(int i=0;i<_threadCount;++i)
        {
        _euclidian.merge(_threadEuclidian[i]);
        _cosine.merge(_threadCosine[i]);
        _tanimoto.merge(_threadTanimoto[i]);
        _sorensen.merge(_threadSorensen[i]);
        _jaccard.merge(_threadJaccard[i]);
#if USE_DTW
        _dtw.merge(_threadDtw[i]);
#endif
        }
      ++iterationCount;
      elapsed=(float)sw.Elapsed.TotalSeconds;
      }
    }
  // display last results to check correctness
  Console.WriteLine("euclidian: {0}",_euclidian);
  Console.WriteLine("   cosine: {0}",_cosine);
  Console.WriteLine(" tanimoto: {0}",_tanimoto);
  Console.WriteLine(" sorensen: {0}",_sorensen);
  Console.WriteLine("  jaccard: {0}",_jaccard);
#if USE_DTW
  Console.WriteLine("      dtw: {0}",_dtw);
#endif
  return iterationCount/elapsed;
  }

protected
void
_similarityTask(int threadIndex)
  {
  Result euclidian=_threadEuclidian[threadIndex];
  euclidian.reset();
  Result cosine=_threadCosine[threadIndex];
  cosine.reset();
  Result tanimoto=_threadTanimoto[threadIndex];
  tanimoto.reset();
  Result sorensen=_threadSorensen[threadIndex];
  sorensen.reset();
  Result jaccard=_threadJaccard[threadIndex];
  jaccard.reset();
#if USE_DTW
  Result dtw=_threadDtw[threadIndex];
  dtw.reset();
#endif
  int beginIndex=_table.Length*threadIndex/_threadCount;
  int endIndex=Math.Min(_table.Length*(threadIndex+1)/_threadCount,
                        _table.Length);
  for(int j=beginIndex;j<endIndex;++j)
    {
    float[] line=_table[j];
    float sumDelta2=0.0f;
    float sumAB=0.0f;
    float sumA2=0.0f;
    float sumB2=0.0f;
    float sumMin=0.0f;
    float sumMax=0.0f;
    for(int i=0;i<_input.Length;++i)
      {
      float a=line[i];
      float b=_input[i];
      float delta=a-b;
      sumDelta2+=delta*delta;
      sumAB+=a*b;
      sumA2+=a*a;
      sumB2+=b*b;
      sumMin+=Math.Min(a,b);
      sumMax+=Math.Max(a,b);
      }
    // max distance is sqrt(inputLength)
    euclidian.insert(j,1.0f-(float)Math.Sqrt(sumDelta2/_input.Length));
    // ignore negative cosine
    cosine.insert(j,(float)Math.Max(0.0,sumAB/Math.Sqrt(sumA2*sumB2)));
    tanimoto.insert(j,sumAB/(sumA2+sumB2-sumAB));
    sorensen.insert(j,2.0f*sumAB/(sumA2+sumB2));
    jaccard.insert(j,sumMin/sumMax);
#if USE_DTW
    // max DTW similarity is inputLength
    dtw.insert(j,_computeDtw(line)/_input.Length);
#endif
    }
  }

#if USE_DTW
protected
float
_computeDtw(float[] line)
  {
  // no need to allocate a huge array, only the previous row is needed while
  // computing a row (this prevents from unnecessary cache-trashing)
  float[,] array=new float[2,_input.Length];
  int previousIndex=0,currentIndex=1;
  array[currentIndex,0]=0.0f;
  for(int j=1;j<_input.Length;++j)
    {
    array[currentIndex,j]=float.MaxValue;
    }
  float previousResult=0.0f; // any value suitable here, will be overwritten
  for(int i=0;i<line.Length;++i)
    {
    previousIndex=currentIndex;
    currentIndex=1-currentIndex;
    previousResult=float.MaxValue;
    float refValue=line[i];
    for(int j=0;j<_input.Length;++j)
      {
      float nextResult=Math.Abs(refValue-_input[j])+
                       Math.Min(previousResult,array[previousIndex,j]);
      // precompute min here in order to save one reread in next row
      array[currentIndex,j]=Math.Min(previousResult,nextResult);
      previousResult=nextResult;
      }
    }
  return previousResult;
  }
protected
float
_computeDtwOriginal(float[] line)
  {
  float[,] array=new float[1+line.Length,1+_input.Length];
  array[0,0]=0.0f;
  for(int j=0;j<_input.Length;++j)
    {
    array[0,j+1]=float.MaxValue;
    }
  for(int i=0;i<line.Length;++i)
    {
    float value=line[i];
    array[i+1,0]=float.MaxValue;
    for(int j=0;j<_input.Length;++j)
      {
      array[i+1,j+1]=Math.Abs(value-_input[j])+
                     Math.Min(array[i+1,j],
                              Math.Min(array[i,j],
                                       array[i,j+1]));
      }
    }
  return array[line.Length,_input.Length];
  }
#endif

protected float[] _input;
protected float[][] _table;
protected int _threadCount;
protected Random _rand;
protected Result[] _threadEuclidian;
protected Result _euclidian;
protected Result[] _threadCosine;
protected Result _cosine;
protected Result[] _threadTanimoto;
protected Result _tanimoto;
protected Result[] _threadSorensen;
protected Result _sorensen;
protected Result[] _threadJaccard;
protected Result _jaccard;
#if USE_DTW
protected Result[] _threadDtw;
protected Result _dtw;
#endif
}

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

public class Test
{
static
void
taskA(int iterationIndex,
      int threadIndex)
  {
  Console.WriteLine("task A: iteration {0}, thread {1}",
                    iterationIndex,threadIndex);
  }

static
void
taskB(int iterationIndex,
      int threadIndex)
  {
  Console.WriteLine("task B: iteration {0}, thread {1}",
                    iterationIndex,threadIndex);
  }

public static
void
Main(string[] args)
  {
  int threadCount=Environment.ProcessorCount;
  bool simple=false;
  bool spinning=false;
  for(int i=0;i<args.Length;++i)
    {
    int tmp;
    if(args[i]=="simple")
      {
      simple=true;
      }
    else if(args[i]=="spin")
      {
      spinning=true;
      }
    else if(int.TryParse(args[i],out tmp))
      {
      threadCount=tmp;
      }
    }
  if(simple)
    {
    Console.WriteLine("simpleTest: threadCount={0}, spinning={1}",
                      threadCount,spinning);
    // ensure cleanup is done through Dispose() 
    using(Runner runner=Runner.create(threadCount,spinning))
      {
      for(int i=0;i<2;++i)
        {
        runner.run(delegate(int t) { taskA(i,t); });
        runner.run(delegate(int t) { taskB(i,t); });
        }
      }
    }
  else
    {
    Console.WriteLine("workTest: threadCount={0}, spinning={1}",
                      threadCount,spinning);
    int inputWidth=200;
    int tableLength=1000;
    Work w=new Work(inputWidth,tableLength,threadCount);
    float perf=w.run(2.0f,spinning);
    Console.WriteLine("{0} iteration/second",
                      perf);
    }
  }
}

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