Skip to content

Add overflow protection for MaxEnt kernel functions #632

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
fassaad merged 3 commits into from
Jun 3, 2026
Merged

Add overflow protection for MaxEnt kernel functions #632

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
0859709
Add overflow protection tests for MaxEnt kernel functions and update …
johanneshofmann87 Apr 12, 2026
f3dde0c
Fix check subroutine: use ieee_is_finite to detect NaN and Inf before…
Copilot Apr 12, 2026
3b69003
Merge branch 'master' into 631-nans-for-large-betaomega-in-maxent
johanneshofmann87 Jun 3, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
91 changes: 81 additions & 10 deletions Analysis/maxent_wrapper_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,13 @@ Real (Kind=Kind(0.d0)) function XKER_ph(tau,om, beta)

pi = 3.1415927

XKER_ph = (exp(-tau*om) + exp(-( beta - tau )*om ) )/( pi*(1.d0 + exp( - beta * om ) ) )
if ( beta*om > 600 ) then
XKER_ph = (exp(-tau*om) + exp(-( beta - tau )*om )) / pi
elseif ( beta*om < -600 ) then
XKER_ph = (exp((beta-tau)*om) + exp(tau*om)) / pi
else
XKER_ph = (exp(-tau*om) + exp(-( beta - tau )*om ) )/( pi*(1.d0 + exp( - beta * om ) ) )
endif
Comment thread
johanneshofmann87 marked this conversation as resolved.

end function XKER_ph

Expand All @@ -54,7 +60,13 @@ Real (Kind=Kind(0.d0)) function XKER_ph_c(tau,om, beta)

pi = 3.1415927

XKER_ph_c = (exp(-tau*om) + exp(-( beta - tau )*om ) )/( pi*(1.d0 + exp( - beta * om ) ) )
if ( beta*om > 600 ) then
XKER_ph_c = (exp(-tau*om) + exp(-( beta - tau )*om )) / pi
elseif ( beta*om < -600 ) then
XKER_ph_c = (exp((beta-tau)*om) + exp(tau*om)) / pi
else
XKER_ph_c = (exp(-tau*om) + exp(-( beta - tau )*om ) )/( pi*(1.d0 + exp( - beta * om ) ) )
endif

end function XKER_ph_c

Expand All @@ -65,7 +77,13 @@ Real (Kind=Kind(0.d0)) function XKER_pp(tau,om, beta)

pi = 3.1415927

XKER_pp = exp(-tau*om) / ( pi*(1.d0 + exp( - beta * om ) ) )
if ( beta*om > 600 ) then
XKER_pp = exp(-tau*om) / pi
elseif ( beta*om < -600 ) then
XKER_pp = exp((beta-tau)*om) / pi
else
XKER_pp = exp(-tau*om) / ( pi*(1.d0 + exp( - beta * om ) ) )
endif

end function XKER_pp

Expand All @@ -76,7 +94,13 @@ Real (Kind=Kind(0.d0)) function XKER_p(tau,om, beta)

pi = 3.1415927

XKER_p = exp(-tau*om) / ( pi*(1.d0 + exp( - beta * om ) ) )
if ( beta*om > 600 ) then
XKER_p = exp(-tau*om) / pi
elseif ( beta*om < -600 ) then
XKER_p = exp( (beta-tau)*om) / pi
else
XKER_p = exp(-tau*om) / ( pi*(1.d0 + exp( - beta * om ) ) )
endif

end function XKER_p

Expand All @@ -102,7 +126,11 @@ Real (Kind=Kind(0.d0)) function F_QFI_ph(om, beta)
real (Kind=Kind(0.d0)) :: om, beta
real (Kind=Kind(0.d0)) :: pi
pi = 3.1415927
F_QFI_ph = (4.d0/pi) * ( (exp(beta*om) - 1.d0)/( exp(beta*om) + 1.d0 ) )**2
if ( abs(beta*om) > 600 ) then
F_QFI_ph = 4.d0/pi
else
F_QFI_ph = (4.d0/pi) * ( (exp(beta*om) - 1.d0)/( exp(beta*om) + 1.d0 ) )**2
endif

end function F_QFI_ph

Expand All @@ -112,7 +140,11 @@ Real (Kind=Kind(0.d0)) function F_QFI_ph_c(om, beta)
real (Kind=Kind(0.d0)) :: om, beta
real (Kind=Kind(0.d0)) :: pi
pi = 3.1415927
F_QFI_ph_c = (4.d0/pi) * ( (exp(beta*om) - 1.d0)/( exp(beta*om) + 1.d0 ) )**2
if ( abs(beta*om) > 600 ) then
F_QFI_ph_c = 4.d0/pi
else
F_QFI_ph_c = (4.d0/pi) * ( (exp(beta*om) - 1.d0)/( exp(beta*om) + 1.d0 ) )**2
endif

end function F_QFI_ph_c

Expand All @@ -121,7 +153,11 @@ Real (Kind=Kind(0.d0)) function F_QFI_pp(om, beta)
real (Kind=Kind(0.d0)) :: om, beta
real (Kind=Kind(0.d0)) :: pi
pi = 3.1415927
F_QFI_pp = (4.d0/pi) * ( (exp(beta*om) - 1.d0)/( exp(beta*om) + 1.d0 ) )**2
if ( abs(beta*om) > 600 ) then
F_QFI_pp = 4.d0/pi
else
F_QFI_pp = (4.d0/pi) * ( (exp(beta*om) - 1.d0)/( exp(beta*om) + 1.d0 ) )**2
endif

end function F_QFI_pp

Expand Down Expand Up @@ -262,7 +298,13 @@ Real (Kind=Kind(0.d0)) function Back_trans_ph(Aom, om, beta )
Implicit None
real (Kind=Kind(0.d0)), intent(in) :: Aom, beta, om

Back_trans_ph = Aom/(1.d0 + exp(-beta*om) )
if ( beta*om > 600 ) then
Back_trans_ph = Aom
elseif ( beta*om < -600 ) then
Back_trans_ph = Aom * exp(beta*om)
else
Back_trans_ph = Aom/(1.d0 + exp(-beta*om) )
endif
! This gives S(q,om) = chi(q,om)/(1 - e^(-beta om))

end function BACK_TRANS_PH
Expand All @@ -276,6 +318,10 @@ Real (Kind=Kind(0.d0)) function Back_trans_ph_c(Aom, om, beta)
Zero = 1.D-8
if ( abs(om) < zero ) then
Back_trans_ph_c = beta * Aom/2.d0
elseif ( beta*om > 600 ) then
Back_trans_ph_c = Aom / om
elseif ( beta*om < -600 ) then
Back_trans_ph_c = -Aom / om
else
Back_trans_ph_c = Aom * (1.d0 - exp(-beta*om) ) / (om *( 1.d0 + exp(-beta*om) ) )
! This gives sigma'(q,om) = A_c(q,om)*(1 - e^(-beta om))/(1 + e^(-beta om))/om
Expand All @@ -296,6 +342,10 @@ Real (Kind=Kind(0.d0)) function Back_trans_pp(Aom, om, beta)
Zero = 1.D-8
if ( abs(om) < Zero ) then
Back_trans_pp = beta * Aom/2.d0
elseif ( beta*om > 600 ) then
Back_trans_pp = Aom / om
elseif ( beta*om < -600 ) then
Back_trans_pp = -Aom / om
else
Back_trans_pp = Aom * (1.d0 - exp(-beta*om) ) / (om *( 1.d0 + exp(-beta*om) ) )
endif
Expand All @@ -313,7 +363,14 @@ Real (Kind=Kind(0.d0)) function XKER_p_ph(tau,om, beta)
real (Kind=Kind(0.d0)) :: tau, om, pi, beta

pi = 3.1415927
XKER_p_ph = (exp(-tau*om) + exp(-(beta-tau)*om)) / (pi*(1.d0 + exp( -beta * om )) )

if ( beta*om > 600 ) then
XKER_p_ph = (exp(-tau*om) + exp(-(beta-tau)*om)) / pi
elseif ( beta*om < -600 ) then
XKER_p_ph = (exp((beta-tau)*om) + exp(tau*om)) / pi
else
XKER_p_ph = (exp(-tau*om) + exp(-(beta-tau)*om)) / (pi*(1.d0 + exp( -beta * om )) )
endif

end function XKER_p_ph

Expand Down Expand Up @@ -442,10 +499,16 @@ Subroutine Set_default(Default,beta,Channel, OM_st, Om_en, xmom1,Default_model_e
case("PH")
If (.not. Default_model_exists ) Default = 1.d0/(Om_en - Om_st) ! Flat default
!Compute sum rule for A(om)
! Note: PH channel uses symmetric spectrum with Om_st = 0 (enforced in Max_SAC.F90),
! so om > 0 always holds here and the beta*om < -600 branch is never reached.
X = 0.d0
Do nw = 1, Ndis
Om = Om_st + dble(nw)*dom
Default(nw) = (1.d0 + exp(-beta*om)) * Default(nw)
if ( beta*om > 600 ) then
Default(nw) = Default(nw)
else
Default(nw) = (1.d0 + exp(-beta*om)) * Default(nw)
endif
X = X + Default(nw)
enddo
X = X*dom
Expand All @@ -467,6 +530,10 @@ Subroutine Set_default(Default,beta,Channel, OM_st, Om_en, xmom1,Default_model_e
Om = Om_st + dble(nw)*dom
if ( abs(om) < zero ) then
Default(nw) = Default(nw)*2.d0/ beta
elseif ( beta*om > 600 ) then
Default(nw) = Default(nw) * om
elseif ( beta*om < -600 ) then
Default(nw) = Default(nw) * (-om)
else
Default(nw) = Default(nw) * (om *( 1.d0 + exp(-beta*om) ) )/ (1.d0 - exp(-beta*om) )
endif
Expand All @@ -492,6 +559,10 @@ Subroutine Set_default(Default,beta,Channel, OM_st, Om_en, xmom1,Default_model_e
Om = Om_st + dble(nw)*dom
if ( abs(om) < zero ) then
Default(nw) = Default(nw)*2.d0/ beta
elseif ( beta*om > 600 ) then
Default(nw) = Default(nw) * om
elseif ( beta*om < -600 ) then
Default(nw) = Default(nw) * (-om)
else
Default(nw) = Default(nw) * (om *( 1.d0 + exp(-beta*om) ) )/ (1.d0 - exp(-beta*om) )
endif
Expand Down
Loading
Loading