import numpy
import copy
from model import *
from triangle import *
from setmask import *
from parameterize import *
from setflowequation import *
from independent import *
from dependent import *
from SetIceShelfBC import *
from EnumDefinitions import *
from solve import *
from MatlabFuncs import *

#This test runs test3005 with autodiff on, and checks that 
#the value of the scalar forward difference match a step-wise differential

#First configure
md=triangle(model(),'../Exp/Square.exp',50000.)
md=setmask(md,'all','')
md=parameterize(md,'../Par/SquareShelfConstrained.py')
md=setflowequation(md,'macayeal','all')
md.cluster=generic('name',oshostname(),'np',3)
md.prognostic.requested_outputs=IceVolumeEnum()
md.verbose=verbose('autodiff',True)

#setup autodiff parameters
index=1 #this is the scalar component we are checking against
md.autodiff.independents=[\
	independent('name','Thickness','type','vertex','nods',md.mesh.numberofvertices,'fos_forward_index',index)
	]

md.autodiff.dependents=[\
	dependent('name','IceVolume','type','scalar')\
	]
md.autodiff.driver='fos_forward'

#parameters for the step-wise derivative
delta=0.001
h1=md.geometry.thickness[index]
h0=h1*(1.-delta)
h2=h1*(1.+delta)
deltaH=(h2-h0)

#save model
md2=copy.deepcopy(md)

#evaluate derivative by forward and backward stepping 
#forward
md=copy.deepcopy(md2)
md.autodiff.isautodiff=False
md.geometry.thickness[index]=h0
md.geometry.bed=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness
md.geometry.surface=md.geometry.bed+md.geometry.thickness
md=SetIceShelfBC(md)

md=solve(md,PrognosticSolutionEnum())
V0=md.results.PrognosticSolution.IceVolume

#backward
md=copy.deepcopy(md2)
md.autodiff.isautodiff=False
md.geometry.thickness[index]=h2
md.geometry.bed=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness
md.geometry.surface=md.geometry.bed+md.geometry.thickness
md=SetIceShelfBC(md)

md=solve(md,PrognosticSolutionEnum())
V2=md.results.PrognosticSolution.IceVolume

#compute resulting derivative
dVdh_an=(V2-V0)/deltaH

#evaluate derivative using ADOLC 
md=md2
md.autodiff.isautodiff=True
md.geometry.thickness[index]=h1
md.geometry.bed=-md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness
md.geometry.surface=md.geometry.bed+md.geometry.thickness
md=SetIceShelfBC(md)

md=solve(md,PrognosticSolutionEnum())
#retrieve directly
dVdh_ad=md.results.PrognosticSolution.AutodiffJacobian

print "dV/dh: analytical:  #16.16g\n       using adolc:  #16.16g\n" % (dVdh_an,dVdh_ad)

#Fields and tolerances to track changes
field_names     =['dV/dh-dV/dh0']
field_tolerances=[1e-13]
field_values=[dVdh_ad-dVdh_an]
