#!/usr/bin/env luajit

--[[

   Copyright (c) 2016 Idiap Research Institute, http://www.idiap.ch/
   Written by Francois Fleuret <francois.fleuret@idiap.ch>

   This file is free software: you can redistribute it and/or modify
   it under the terms of the GNU General Public License version 3 as
   published by the Free Software Foundation.

   It 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 this file.  If not, see <http://www.gnu.org/licenses/>.

]]--

require 'torch'
require 'nn'

require 'dagnn'

function checkGrad(model, criterion, input, target)
   local params, gradParams = model:getParameters()

   local epsilon = 1e-5

   local output = model:forward(input)
   local loss = criterion:forward(output, target)
   local gradOutput = criterion:backward(output, target)
   gradParams:zero()
   model:backward(input, gradOutput)
   local analyticalGradParam = gradParams:clone()

   for i = 1, params:size(1) do
      local x = params[i]

      params[i] = x - epsilon
      local output0 = model:forward(input)
      local loss0 = criterion:forward(output0, target)

      params[i] = x + epsilon
      local output1 = model:forward(input)
      local loss1 = criterion:forward(output1, target)

      params[i] = x

      local ana = analyticalGradParam[i]
      local num = (loss1 - loss0) / (2 * epsilon)
      local err

      if num == ana then
         err = 0
      else
         err = torch.abs(num - ana) / torch.abs(num)
      end

      print(
         'CHECK '
            .. err
            .. ' checkGrad ' .. i
            .. ' analytical ' .. ana
            .. ' numerical ' .. num
      )
   end

end

function printTensorTable(t)
   if torch.type(t) == 'table' then
      for i, t in pairs(t) do
         print('-- ELEMENT [' .. i .. '] --')
         printTensorTable(t)
      end
   else
      print(tostring(t))
   end
end

-- torch.setnumthreads(params.nbThreads)
torch.setdefaulttensortype('torch.DoubleTensor')
torch.manualSeed(2)

--                     +--> c ----> e --+
--                    /            /     \
--                   /            /       \
--  input --> a --> b ---> d ----+         g --> output
--                          \             /
--                           \           /
--                            +--> f ---+

a = nn.Linear(10, 10)
b = nn.ReLU()
c = nn.Linear(10, 3)
d = nn.Linear(10, 3)
e = nn.CMulTable()
f = nn.Linear(3, 3)
g = nn.CAddTable()

model = nn.DAG()

model:addEdge(a, b)
model:addEdge(b, c)
model:addEdge(b, d)
model:addEdge(c, e)
model:addEdge(d, e)
model:addEdge(d, f)
model:addEdge(e, g)
model:addEdge(f, g)

model:setInput(a)
model:setOutput(g)

input = torch.Tensor(3, 10):uniform()

print('******************************************************************')
print('** updateOutput **************************************************')
print('******************************************************************')

output = model:updateOutput(input):clone()

printTensorTable(output)

print('******************************************************************')
print('** updateGradInput ***********************************************')
print('******************************************************************')

gradInput = model:updateGradInput(input, output)

printTensorTable(gradInput)

print('******************************************************************')
print('** checkGrad *****************************************************')
print('******************************************************************')

output:uniform()

checkGrad(model, nn.MSECriterion(), input, output)