module adjcvmod
	use basicmodule
	implicit none
	
	integer, parameter	:: nax=3
	real				:: cvcoefs(nax)
	real				:: linexpar=1,aNh
	real				:: level

	contains
	
	function stat_adj(x,stats) result(val)
		real, intent(in)	:: x(nax), stats(:)
		real				:: val
		val=stats(1)*exp(x(1)+x(2)*stats(2)+x(3)*stats(2)**2)
	end function
	
	function aNquants(x) result(val)
		real	:: x(nax), val(size(allstats,3))
		integer	:: i,j,l
		real	:: c
		val=0
!$omp parallel do private(c,l)
		do j=1,size(allstats,3)
			c=0
			do l=1,size(allstats,2)
				c=c+gausscdf((1-stat_adj(x,allstats(:,l,j)))/aNh)
			enddo
			val(j)=c/size(allstats,2)
		enddo
	end function

	function anSize(x) result(val)
		real	:: x(nax), val(size(allstats,3))
		real	:: stats(size(allstats,2)),quants(size(allstats,3))
		integer	:: i,j,l
		real	:: c
		c=0
!$omp parallel do private(stats,l)
		do j=1,size(allstats,3)
			do l=1,size(allstats,2)
				stats(l)=stat_adj(x,allstats(:,l,j))
			enddo
			quants(j)=quantile_s(stats,1-level)
		enddo


		c=maxval(quants)		
		x(1)=x(1)-log(c)
!$omp parallel do private(c,l)
		do j=1,size(allstats,3)
			c=0
			do l=1,size(allstats,2)
				if(stat_adj(x,allstats(:,l,j))>1) c=c+1
			enddo
			val(j)=c/size(allstats,2)
		enddo
	end function
	
	elemental function linex(x) result(val)
		real, intent(in)	:: x
		real				:: val
		val=exp(x)-x-1
	end function
	
	function adobj(x) result(val)
		real	:: x(nax), val
		val=sum(linex(-linexpar*(invlogit(aNquants(x))-invlogit(1-level))))
		val=val+.001*sum(x(2:)**2)
	end function
					
	subroutine findajdcv
		use uminf_int
		use U4INF_int
		external	::  ex_adjNyb
		real	:: x(nax),x0(nax),f,rparam(7),bestf,bestx(nax),qs(2)
		integer	:: iparam(7),j
		real, allocatable	:: dum(:,:)
		bestf=1E100
		qs=.3*quantile_v(allstats(1,:,size(allstats,3)/2),[.93,.97])
		aNh=(qs(2)-qs(1))
		print *,"bandwidth for quantile estimation",aNh

		call u4inf(iparam,rparam)
		iparam(3:5)=1000
		call erset(0,0,0)
		x0=0
		linexpar=.3
		call uminf(ex_adjNyb,x,xguess=x0,iparam=iparam,rparam=rparam)
		linexpar=3
		call u4inf(iparam,rparam)
		iparam(3:5)=1000
		x0=x
		call uminf(ex_adjNyb,x,xguess=x0,iparam=iparam,rparam=rparam)
		linexpar=12
		call u4inf(iparam,rparam)
		iparam(3:5)=1000
		x0=x
		call uminf(ex_adjNyb,x,xguess=x0,iparam=iparam,rparam=rparam)
		call mdisp(aNsize(x))
		call ex_adjNyb(nax,x,f)
		call mdisp([f,x])
		cvcoefs=x
		call savecvcoefs
	end subroutine
	
	subroutine savecvcoefs
		call savecsv(trim(dir)//"cvx_"//getfilesuffix()//".csv",reshape(cvcoefs,[1,nax]))
	end subroutine
	
	subroutine loadcvcoefs
		real, allocatable :: mdata(:,:)
		mdata=loadcsv(trim(dir)//"cvx_"//getfilesuffix()//".csv")
		cvcoefs=[mdata]
	end subroutine
end module
	
subroutine ex_adjNyb(n0,x,f)	
	use adjcvmod
	implicit none
	integer	:: n0
	real	:: x(n0),f
	
	f=adobj(x)
	if((f.ne.f) .or. f>1E50) f=1E100
end subroutine	
	
	