#!/usr/bin/env ruby
require 'pp'

$:<<".."

require 'tools.rb'
require 'matrix_math.rb'


module PCA
  class Input
    def initialize(i=nil)
      @vecs=[]
      if i
        if i.is_a?(Array)
          i.each{|k|self<<k}
        end
      end
    end
    def <<(v)
      assert {v.is_a?(Array)}
      if @vecs.length==0
        @vecs<<v
      else
        assert {@vecs[0].length == v.length }
        @vecs<<v
      end
    end
    
    def trainCount
      @vecs.length
    end
    
    def featureCount
      @vecs[0].length
    end
    
    def dmatrix
      Linalg::DMatrix[*@vecs]
    end
  end
  
  class Base
    attr_reader :eigenvalues
    def initialize(input)
      assert { input.is_a?(Input) }
      @input=input
      
      computePCA
      @input=nil
    end
    
    def dimension
      @base.hsize
    end
    
    def reduce(size)
      assert{size<=@eigenvalues.length}
      return if size==@eigenvalues.length
      @eigenvalues=@eigenvalues[0...size]
      
      @base=@base.cut(size,@base.vsize)
      @baseInv=@base.pseudo_inverse
    end
    
    def compress(v)
      v=Linalg::DMatrix[v]
      v2=v-@middle
      #pp "---",v2,@baseInv,"---"
      @baseInv*v2.transpose
      
      #pp v2
    end
    
    def decompress(v)
      #pp "---",v,@base,"---"
      (@base*v).transpose+@middle
    end
    
    def measureErrors(input)
      errs=input.map{|i|measureSingleError(i)}
      
      mean=errs.inject(0){|old,new|
       # pp "#{old} #{new}"
        old+new
      }/input.length
      [errs.max,errs.min,mean]
    end
    
    def measureSingleError(input)
      inter=compress(input)
      output=decompress(inter)
        pp inter
      input=Linalg::DMatrix[input]
      diff=input-output
      (diff*diff.transpose)[0,0]
    end

    
    private
    def computePCA
      #return doOrthoPCA if @input.trainCount<@input.FeatureCount
      
      dm=@input.dmatrix
      @middle=dm.middle_row
      puts "COV..."
      cov=dm.covariance
      puts "SVD..."
      u, s, vt = cov.singular_value_decomposition
      @eigenvalues=(0...s.vsize).map{|x|s[x,x]}
      @base=u
      
      #sort
      puts "INV"
      @baseInv=@base.pseudo_inverse
      puts "READY"
    end
    def sort
      pp self
      a=@base.to_a
      x=
      (0...a.length).to_a.map{|i|
        [a[i],@eigenvalues[i]]
      }
      y=x.sort{|a,b| a[1]<=>b[1]}.reverse
      
      z=y.map{|a|a[0]}
      #pp x,y,z
      
      @base=Linalg::DMatrix[*@base]
      @eigenvalues.sort!.reverse!
      pp self
      exit
    end
  end
  
  def doOrthoPCA

  trainCount=@input.trainCount
  featureCount=@input.featureCount
  
  #middle=@input.middle
  d=@input.dmatrix
  @middle=d.middle_row
  dt=d.transpose
  
  t=d*dt
  
  t=t*(1.0/trainCount)
  
  eigenValues=eigenVectors=nil
  t.eigensystem {|_eigenVectors,_eigenValues,img|
    eigenValues=_eigenValues
    eigenVectors=_eigenVectors
  }
  
  nvr=dT*eigenVectors

  #Tools::Math::Matrix nvr(trainCount,featureCount);

  #nvr=dT*vr;


#  PCAOutput o;

  
  
#  for(size_t y=0;y<ev.size();y++)
#    {
#      PCAVector v;
#      double scale;

#      if(ev[y]>0)
#  scale=sqrt(ev[y]*trainCount);
#      else
#  scale=0;

#      if(scale<=0.0f)
#  scale=0;
#      else
#  scale=1.0/scale;


#      for(size_t x=0;x<nvr.nrows();x++)
#  v.push_back(nvr(x,y)*scale);
#      o.vectors.push_back(v);
#      o.values.push_back(ev[y]>0?ev[y]:0);
#    }
  end
  
end

if false
i=PCA::Input.new
i<<[1,1,0]
i<<[0,1,1]
i<<[0,1,2]
i<<[9,1,1]
i<<[0,1,0]

base=PCA::Base.new(i)

#pp base

c=base.compress([5,1,0])
#pp c
r=base.decompress(c)
pp r

    
end