#!/usr/bin/env python2.3
################################################################################
#
#       This file is part of the GQL (Graphical Query Language) Toolkit
#
#       file:   Plotting
#       author: Ivan Costa (filho@molgen.mpg.de)
#
#       Copyright (C) 2003-2009 lexander Schliep and Ivan Costa
#
#       Contact: filho@molgen.mpg.de
#
#       Information: http://ghmm.org/gql
#
#   GQL is free software; you can redistribute it and/or modify
#   it under the terms of the GNU General Public License as published by
#   the Free Software Foundation; either version 2 of the License, or
#   (at your option) any later version.
#
#   GQL is distributed in the hope that it will be useful,
#   but WITHOUT ANY WARRANTY; without even the implied warranty of
#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#   GNU General Public License for more details.
#
#   You should have received a copy of the GNU General Public License
#   along with GQL; if not, write to the Free Software
#   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
#
#
#
#       This file is version $Revision: 2125 $
#                       from $Date: 2009-10-22 17:32:35 -0300 (Thu, 22 Oct 2009) $
#             last change by $Author: filho $.
#
################################################################################


import GQL
import pylab
import numpy
import scipy

################################################################################
#      Auxiliary functions for creating time series plots
################################################################################

def plotTimeSeries(profileSet,clusters):
  colors =  ['b','g','r','c','m','y','k','#808080']
  series = []

  #excluding end symbols
  for s in profileSet.profile:
      s = numpy.array(s)
      series.append(s[range(len(s)-1)].tolist())
  series = numpy.array(series)
  clusters = numpy.array(clusters)
  noMin = min(clusters)
  if noMin != 0:
      clusters=clusters-noMin
  noClusters = max(clusters)+1

  # defining ranges
  maxV = numpy.max(series)
  minV = numpy.min(series)
  
  auxseries = []
  fig = pylab.figure()
  means = []
  stds = [] 
  for i in range(noClusters):
      auxseries = numpy.transpose(series[clusters==i])      
      fig.add_subplot(331 + i)
      pylab.plot(range(len(auxseries)),auxseries,color=colors[i])
      pylab.title('Cluster '+str(noClusters))
      
      pylab.axis([0,len(auxseries)-1,minV,maxV])
      means.append(scipy.mean(series[clusters==i],axis=0))
      stds.append(scipy.std(series[clusters==i],axis=0))
  fig.add_subplot(331 + noClusters)
  print means
  pylab.plot(range(len(means[0])),
                       numpy.transpose(means))
  pylab.title('Mean of Clusters') 
  
  #pylab.show()


if __name__ == '__main__':

   import sys
   import GQLCluster
   
   fileName = sys.argv[1]

   profileSet = GQL.ProfileSet()
   profileSet.ReadDataFromCaged(fileName,end=0)
 
   plotTimeSeries(profileSet,profileSet.seq_classes)
   pylab.savefig(fileName[:-4]+'.pdf')
