MODULE minimize
        USE nrtype
        USE interfaces, ONLY: mrqmin_TGM,mrqmin_ABD
        USE error
        USE share, ONLY: convABD,convTGM
        CONTAINS
!####################################
        SUBROUTINE minimTGM_Q(w_sp,f_sp_norm,sig_noise,TGM,TGM_mask&
        &,covarTGM,alpha_covarTGM,chisq)
        IMPLICIT NONE
        REAL(DP), DIMENSION(:), INTENT(IN) :: w_sp,f_sp_norm,sig_noise
        REAL(DP), DIMENSION(:), INTENT(INOUT) :: TGM
        REAL(DP), DIMENSION(SIZE(TGM)) :: TGM_old
        REAL(DP), DIMENSION(:,:), INTENT(INOUT) :: covarTGM,alpha_covarTGM
        REAL(DP), INTENT(OUT) :: chisq
        REAL(DP) :: alamda,tol,conv,chisq_old
        LOGICAL(LGT), DIMENSION(:), INTENT(IN) :: TGM_mask

        INTEGER(I1B) :: i
        alamda=-1.
        conv=1.
        tol=0.01_dp
        chisq=1e9_dp
        chisq_old=1e9_dp


        DO i=1,101
        TGM_old=TGM
!        write(*,*) 'start mrqmin_TGM'
        CALL mrqmin_TGM(w_sp,f_sp_norm,sig_noise,TGM,TGM_mask&
        &,covarTGM,alpha_covarTGM,chisq,func_TGM_Q,alamda)
        !x=w_sp;y=f_sp_norm;sig=sig_noise;a=TGM;maska=TGM_mask;covar=covar;alpha=alpha_covar;
        !chisq=chisq;alamda;

!        write(*,*) TGM
!        write(*,*) 'chi2', chisq_old,chisq,conv,alamda,TGM
        IF(chisq<chisq_old) THEN
!          IF(chisq<10e-30) chisq=10e-30
          conv=(chisq_old-chisq)/(chisq+chisq_old)
          chisq_old=chisq
          IF(conv<=tol) EXIT
        END IF
        IF(ALL(TGM_old==TGM)) THEN
!         write(*,*) 'TGM minim makes no progress'
!         convTGM=.FALSE.
         EXIT
        END IF
        IF(i>100) THEN
!         CALL warning_msg('reached max number of iterations in mrqmin_TGM')
         write(*,*) 'reached max number of iterations in mrqmin_TGM_Q'
         convTGM=.FALSE.
         EXIT
        END IF
        END DO
        alamda=0
        CALL mrqmin_TGM(w_sp,f_sp_norm,sig_noise,TGM,TGM_mask&
        &,covarTGM,alpha_covarTGM,chisq,func_TGM_Q,alamda)
        END SUBROUTINE minimTGM_Q
!###################################

        SUBROUTINE minimTGM(w_sp,f_sp_norm,sig_noise,TGM,TGM_mask&
        &,covarTGM,alpha_covarTGM,chisq)
        IMPLICIT NONE
        REAL(DP), DIMENSION(:), INTENT(IN) :: w_sp,f_sp_norm,sig_noise
        REAL(DP), DIMENSION(:), INTENT(INOUT) :: TGM
        REAL(DP), DIMENSION(SIZE(TGM)) :: TGM_old
        REAL(DP), DIMENSION(:,:), INTENT(INOUT) :: covarTGM,alpha_covarTGM
        REAL(DP), INTENT(OUT) :: chisq
        REAL(DP) :: alamda,tol,conv,chisq_old
        LOGICAL(LGT), DIMENSION(:), INTENT(IN) :: TGM_mask

        INTEGER(I1B) :: i
        alamda=-1.
        conv=1.
        tol=0.01_dp
        chisq=1e9_dp
        chisq_old=1e9_dp


        DO i=1,101
        TGM_old=TGM
!        write(*,*) 'start mrqmin_TGM'
        CALL mrqmin_TGM(w_sp,f_sp_norm,sig_noise,TGM,TGM_mask&
        &,covarTGM,alpha_covarTGM,chisq,func_TGM,alamda)
        !x=w_sp;y=f_sp_norm;sig=sig_noise;a=TGM;maska=TGM_mask;covar=covar;alpha=alpha_covar;
        !chisq=chisq;alamda;

!        write(*,*) TGM
!        write(*,*) 'chi2', chisq_old,chisq,conv,alamda,TGM
        IF(chisq<chisq_old) THEN
!          IF(chisq<10e-30) chisq=10e-30
          conv=(chisq_old-chisq)/(chisq+chisq_old)
          chisq_old=chisq
          IF(conv<=tol) EXIT
        END IF
        IF(ALL(TGM_old==TGM)) THEN
!         write(*,*) 'TGM minim makes no progress'
!         convTGM=.FALSE.
         EXIT
        END IF
        IF(i>100) THEN
!         CALL warning_msg('reached max number of iterations in mrqmin_TGM')
         write(*,*) 'reached max number of iterations in mrqmin_TGM'
         convTGM=.FALSE.
         EXIT
        END IF
        END DO
        alamda=0
        CALL mrqmin_TGM(w_sp,f_sp_norm,sig_noise,TGM,TGM_mask&
        &,covarTGM,alpha_covarTGM,chisq,func_TGM,alamda)
        END SUBROUTINE minimTGM
!###############        
        SUBROUTINE minimABD(w_sp,f_sp_norm,sig_noise,ABD,ABD_mask&
        &,covarABD,alpha_covarABD,chisq)
        IMPLICIT NONE
        REAL(DP), DIMENSION(:), INTENT(IN) :: w_sp,f_sp_norm,sig_noise
        REAL(DP), DIMENSION(:), INTENT(INOUT) :: ABD
        REAL(DP), DIMENSION(SIZE(ABD)) :: ABD_old
        REAL(DP), DIMENSION(:,:), INTENT(INOUT) :: covarABD,alpha_covarABD
        REAL(DP), INTENT(OUT) :: chisq
        REAL(DP) :: alamda,tol,conv,chisq_old
        LOGICAL(LGT), DIMENSION(:), INTENT(IN) :: ABD_mask
        INTEGER(I1B) :: i
        alamda=-1.
        tol=0.1_dp
        chisq=1e9_dp
        chisq_old=1e9_dp

        DO i=1,101
!        write(*,*) i,'minimize.f95 1', ABD
        CALL mrqmin_ABD(w_sp,f_sp_norm,sig_noise,ABD,ABD_mask&
        &,covarABD,alpha_covarABD,chisq,func_ABD,alamda)
!        write(*,*) 'chi2 ABD', chisq_old,chisq,conv,alamda
!        write(*,*) 'minimize.f95 2', chisq
        IF(chisq<chisq_old) THEN
!          IF(chisq<10e-30) chisq=10e-30
          conv=(chisq_old-chisq)/(chisq+chisq_old)
          chisq_old=chisq
          IF(conv<=tol) EXIT
        END IF
        IF(ALL(ABD_old==ABD)) THEN
!         write(*,*) 'ABD minim makes no progress',chisq
         convABD=.FALSE.
         EXIT
        END IF
        IF(i>100) THEN
!         CALL warning_msg('reached max number of iterations in mrqmin_ABD')
         write(*,*) 'reached max number of iterations in mrqmin_ABD'
         convABD=.FALSE.
         EXIT
        END IF
        END DO
        alamda=0.
        CALL mrqmin_ABD(w_sp,f_sp_norm,sig_noise,ABD,ABD_mask&
        &,covarABD,alpha_covarABD,chisq,func_ABD,alamda)
        END SUBROUTINE minimABD
!#########################
SUBROUTINE func_TGM_Q(x,a,yfit,dyda)
        USE make_model, ONLY: make_model_TGM_quick
        USE share, ONLY: TGM_mask
        IMPLICIT NONE
        REAL(DP), DIMENSION(:), INTENT(IN) :: x,a
        REAL(DP), DIMENSION(SIZE(a)) :: a_new
        REAL(DP), DIMENSION(:), INTENT(OUT) :: yfit
        REAL(DP), DIMENSION(SIZE(x)) :: yfit_new
        REAL(DP), DIMENSION(:,:), INTENT(OUT) :: dyda
        REAL(DP),DIMENSION(SIZE(a)) :: incr
        REAL(DP), PARAMETER :: eps=5e-3_dp
        INTEGER(I1B) :: j

        dyda=0.
        incr=ABS(a)*eps
        WHERE(incr<=0.) incr=eps        
        CALL make_model_TGM_quick(yfit,a)

        a_new=a
        DO j=1,INT(SIZE(a,1),I1B)
         IF(TGM_mask(j)) THEN
           a_new(j)=a(j)+incr(j)
!        write(*,*) j,a(j),a_new(j),incr(j)
           CALL make_model_TGM_quick(yfit_new,a_new)
           dyda(:,j)=(yfit_new-yfit)/incr(j)
        !sometimes dyda(j,j)=0 because the derivatives are zero when one
        !parameter reaches the grid border. This makes the matrix singular. 
        !to avoid this, give a fake value and continue.
!        if(ALL(dyda(:,j)==0)) then
!        write(*,*) j,'TGM zero'
!        dyda(j,j)=1.
!        end if
         a_new(j)=a(j)
         END IF
        END DO

END SUBROUTINE func_TGM_Q
!#########################
SUBROUTINE func_TGM(x,a,yfit,dyda)
        USE make_model, ONLY: make_model_TGM
        USE share, ONLY: TGM_mask
        IMPLICIT NONE
        REAL(DP), DIMENSION(:), INTENT(IN) :: x,a
        REAL(DP), DIMENSION(SIZE(a)) :: a_new
        REAL(DP), DIMENSION(:), INTENT(OUT) :: yfit
        REAL(DP), DIMENSION(SIZE(x)) :: yfit_new
        REAL(DP), DIMENSION(:,:), INTENT(OUT) :: dyda
        REAL(DP),DIMENSION(SIZE(a)) :: incr
        REAL(DP), PARAMETER :: eps=5e-3_dp
        INTEGER(I1B) :: j

        dyda=0.
        incr=ABS(a)*eps
        WHERE(incr<=0.) incr=eps        
        CALL make_model_TGM(yfit,a)

        a_new=a
        DO j=1,INT(SIZE(a,1),I1B)
         IF(TGM_mask(j)) THEN
           a_new(j)=a(j)+incr(j)
!        write(*,*) j,a(j),a_new(j),incr(j)
           CALL make_model_TGM(yfit_new,a_new)
           dyda(:,j)=(yfit_new-yfit)/incr(j)
        !sometimes dyda(j,j)=0 because the derivatives are zero when one
        !parameter reaches the grid border. This makes the matrix singular. 
        !to avoid this, give a fake value and continue.
!        if(ALL(dyda(:,j)==0)) then
!        write(*,*) j,'TGM zero'
!        dyda(j,j)=1.
!        end if
         a_new(j)=a(j)
         END IF
        END DO

END SUBROUTINE func_TGM

SUBROUTINE func_ABD(x,a,yfit,dyda)
        USE make_model, ONLY: make_model_ABD
        USE utils, ONLY: ABD_mask
        IMPLICIT NONE
        REAL(DP), DIMENSION(:), INTENT(IN) :: x,a
        REAL(DP), DIMENSION(SIZE(a)) :: a_new
        REAL(DP), DIMENSION(:), INTENT(OUT) :: yfit
        REAL(DP), DIMENSION(SIZE(x)) :: yfit_new
        REAL(DP), DIMENSION(:,:), INTENT(OUT) :: dyda
        REAL(DP),DIMENSION(SIZE(a)) :: incr
        REAL(DP), PARAMETER :: eps=5e-3_dp
        INTEGER(I1B) :: j

        incr=ABS(a)*eps
        WHERE(incr<=0.) incr=eps

        CALL make_model_ABD(yfit,a)

        a_new=a
        
        DO j=1,INT(SIZE(a,1),I1B)
         IF(ABD_mask(j)) THEN
         a_new(j)=a(j)+incr(j)
!        write(*,*) j,a(j),a_new(j)
           CALL make_model_ABD(yfit_new,a_new)
           dyda(:,j)=(yfit_new-yfit)/incr(j)
        !sometimes dyda(j,j)=0 because the derivatives are zero when one
        !parameter reaches the grid border. This makes the matrix singular. 
        !to avoid this, give a fake value and continue.
!        if(ALL(dyda(:,j)==0)) then
!        write(*,*) j,a(j),a_new(j),'ABD zero'
!         dyda(j,j)=1
!        end if

         a_new(j)=a(j)
        END IF
        END DO

END SUBROUTINE func_ABD

END MODULE minimize