module taugrid
	use basicmodule
	use compute
	implicit none
	
	real		:: lllb,msxw(3)
	
	contains
	
	subroutine loadtaudata
		real, allocatable :: mdata(:,:)
		mdata=loadcsv("c:/Dropbox/TVExtremes/data/damages/Disasters_All_Variables.csv",.true.)
		ys(1:6,:)=transpose(mdata(:,3:8))
		tau=1
		yks=mdata(:,2)
		ys(k,:)=tau
	end subroutine
	
	subroutine loadtaudatanorm
		real, allocatable :: mdata(:,:)
		mdata=loadcsv("c:/Dropbox/TVExtremes/data/damages/Disasters_Normalized_All_Variables.csv",.true.)
		ys(1:6,:)=transpose(mdata(:,4:9))
		tau=mdata(:,2)
		yks=mdata(:,3)
		ys(k,:)=tau
	end subroutine

	
	subroutine settaugrid()
		integer, parameter	:: ngrid=50
		real	:: msx(3),stats(nsim),cv0,cv,msxbounds(2,3),u(3)
		real	:: ll0,lltst
		real, allocatable	:: ys0(:,:)
		integer, allocatable	:: yks0(:)
		integer	:: l,i,j
		integer	:: tstind0

		ys0=ys; yks0=yks
		if(all(tau==1)) then
			msx=maxll([0.45,0.88,0.84],.false.)
		else
			msx=maxll([0.33,0.32,0.99],.false.)
		endif
		call mdisp(msx)
		ll0=getconstll(msx)
		do l=1,nsim
			call draw_Ys(msx,l)
			stats(l)=getconstll(maxll(msx,.false.))-getconstll(msx)
		enddo
		cv0=quantile_s(stats,0.995)
		print *,"initial cv=",cv0
		lllb=ll0-cv0
		ys=ys0; yks=yks0
		do i=1,2
			do j=1,3
				msxw=0
				msxw(j)=merge(1,-1,i==1)
				msxbounds(i,j)=exmsx(msx)
			enddo
		enddo
		print *,"bounds on msx based on intial cv"
		call mdisp(msxbounds)
		if(allocated(taumsxgrid)) deallocate(taumsxgrid) 
		allocate(taumsxgrid(3,ngrid))
		j=1
		do 
			call rnun(u)
			msx=msxbounds(1,:)+u*(msxbounds(2,:)-msxbounds(1,:))
			lltst=ll0-getconstll(msx)
			if(lltst>1.2*cv0) cycle
			do l=1,nsim
				call draw_Ys(msx,l)
				stats(l)=getconstll(maxll(msx,.false.))-getconstll(msx)
			enddo
			cv=quantile_s(stats,0.99)
			ys=ys0; yks=yks0
			if(lltst>cv) cycle
			taumsxgrid(:,j)=msx
			j=j+1
			if(j>ngrid) exit
		enddo
		call mdisp(taumsxgrid)
		call savemat(trim(dir)//"taumsxgrid.mat",taumsxgrid)
		call printtime
	end subroutine
	
	subroutine tauobj(msx,iact,f,ierr)
		real	:: msx(3),f
		integer	:: iact
		logical	:: ierr
		ierr=.false.
		if(iact==0) then
			f=sum(msxw*msx)
			return
		endif
		f=getconstll(msx)-lllb
	end subroutine		

	subroutine dtauobj(msx,iact,g)
		real	:: msx(3),g(3)
		real	:: xsi
		integer	:: iact
		if(iact==0) then
			g=msxw
			return
		endif
		g=getdconstll(msx)
	end subroutine		

	function exmsx(msx0) result(val)
		use nnlpg_int
		real		:: msx0(3),val,msx(3)
		real		:: lbounds(3),ubounds(3)
		integer		:: m
		external	:: extauobj, exdtauobj
		
		lbounds=[-1E3,.001,-0.99]
		ubounds=[1E3,1E3,5.0]
		m=1
		call erset(0,0,0)
		call nnlpg(extauobj,exdtauobj,m,m,0,lbounds,ubounds,msx,xguess=msx0,iprint=0)
		val=sum(abs(msxw)*msx)
	end function
	

end module


	subroutine extauobj(msx,iact,f,ierr)
		use taugrid
		real	:: msx(3),f
		integer	:: iact
		logical	:: ierr	
		call tauobj(msx,iact,f,ierr)
	end subroutine
	
	subroutine exdtauobj(msx,iact,g)
		use taugrid
		real	:: msx(3),g(3)
		integer	:: iact
		call dtauobj(msx,iact,g)
	end subroutine
