module cymod			! routines for Campbell Yogo application
	use compute
	implicit none
	
	contains
	
	function getparafromx(x) result(para)
		real	:: x(:)
		type(tpara)	:: para
		para%hAR(1:p)=hlb+x(1:p)*hub
		para%hMA(1:p-1)=hlb+x(p+1:2*p-1)*hub
	end function
	
	subroutine maxLL(LL)		! maximize (profiled) likelihood and store maximum
		use bconf_int
		real	:: LL
		external		:: evalll_ext
		integer, parameter	:: nruns=100
		real		:: x(2*p-1),lb(2*p-1), ub(2*p-1)
		type(tpara)	:: para
		real		:: rparam(7)
		integer	:: iparam(7)
		integer	:: i
		real		:: bestx(2*p-1,nruns), cll, bestll(nruns),om2s(nruns)
		real	,save	:: outtab(2,5)=-1

		call du4inf(iparam,rparam)
		iparam(3)=200

		lb=1E-5
		ub=1		
		do i=1,nruns
			call rnun(x)
			call erset(0,0,0)
			call bconf(evalll_ext,0,lb,ub,bestx(:,i),xguess=x)
			call evalll_ext(2*p-1,bestx(:,i),bestll(i))
			bestll(i)=-bestll(i)
			om2s(i)=om2
		enddo
		call mdisp(bestll)
		i=maxloc(bestll,dim=1)
		
		print *,"maximum likelihood for p=",p
		print *,bestll(i)
		para=getparafromx(bestx(:,i))
		print *,"h"
		call mdisp([para%har(1:p),para%hma(1:p-1)])
		print *,"om2"
		print *,om2s(i)
		print *,"implied c and g values (real part in top row, imaginary part in bottom)"
		call mdisp(real([getcfromh(para%har(1:p)),getcfromh(para%hma(1:p-1))]).cvr.aimag([getcfromh(para%har(1:p)),getcfromh(para%hma(1:p-1))]))
		LL=bestll(i)
	end subroutine
	
	subroutine genpoints(LL)		! generate random points within 2 log points of maximum likelihood
		real		:: LL
		integer, parameter	:: npoints=1000
		real		:: x(2*p-1),paralist(0:2*p,npoints)
		real		:: cll
		type(tpara)	:: para
		integer	:: i,j
		logical	:: exitflag
		
		i=0
		exitflag=.false.
		do j=1,10*npoints
			if(exitflag) cycle
			call rnun(x)
			para=getparafromx(x)
			cll=getKalman_ll(para)
			cll=cll-LL
			if(abs(cll)<2) then
				i=i+1
				i=min(i,npoints)
				if(i==npoints) exitflag=.true.
				paralist(:,i)=[cll,om2,para%hAR(1:p),para%hMA(1:p-1)]
				if(i==1) then
					print *,"c and g for first column"
					call mdisp(real([getcfromh(para%har(1:p)),getcfromh(para%hma(1:p-1))]).cvr.aimag([getcfromh(para%har(1:p)),getcfromh(para%hma(1:p-1))]))
				endif		
			endif
		enddo
		print *,i
		call savemat(trim(dir)//"CYpoints.txt",paralist)
		call storeML(paralist,"paralist")
		print *,"done"
	end subroutine


end module
	
subroutine evalll_ext(n,x,val)
	use CYmod
	implicit none
	integer	:: n
	real		:: x(2*p-1),val
	val=-getKalman_ll(getparafromx(x))
	if(.not. isfinite(val)) val=1E100
end subroutine	