Skip to content

Commit

Permalink
Merge pull request #487 from galacticusorg/brownianBridgeFix
Browse files Browse the repository at this point in the history
Fix mathematical implementation of Brownian bridge solutions
  • Loading branch information
abensonca authored Oct 3, 2023
2 parents 03a6865 + c6f4d6d commit 82f6406
Show file tree
Hide file tree
Showing 11 changed files with 911 additions and 125 deletions.
11 changes: 11 additions & 0 deletions doc/Galacticus.bib
Original file line number Diff line number Diff line change
Expand Up @@ -4441,6 +4441,17 @@ @article{kummer_effective_2018
pages = {388--399}
}

@misc{kiwiakos_answer_2014,
title = {Answer to "{Brownian} {Bridge}'s first passage time distribution"},
url = {https://quant.stackexchange.com/a/14219},
urldate = {2023-09-28},
journal = {Quantitative Finance Stack Exchange},
author = {Kiwiakos},
month = jul,
year = {2014},
file = {Snapshot:/home/abensonca/.mozilla/firefox/f54gqgdx.default/zotero/storage/RQ6Q725B/14219.html:text/html},
}

@article{lacey_merger_1993,
title = {Merger rates in hierarchical models of galaxy formation},
volume = {262},
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -71,10 +71,10 @@
type (rootFinder ) :: finder
type (integrator ) :: integrator_ , integratorSubresolution_
! Parent halo shared variables.
double precision :: parentDTimeDDeltaCritical , parentDelta , &
& parentHaloMass , parentSigma , &
& parentSigmaSquared , parentTime , &
& probabilityMinimumMass , probabilitySeek , &
double precision :: parentDTimeDDeltaCritical , parentDelta , &
& parentHaloMass , parentSigma , &
& parentSigmaSquared , parentTime , &
& probabilityMinimumMass , probabilitySeek , &
& normalization

! Record of mass resolution.
Expand All @@ -92,7 +92,7 @@
! Record of issued warnings.
logical :: subresolutionFractionIntegrandFailureWarned
! Option controlling whether only lower-half of the distribution function should be used.
logical :: distributionFunctionLowerHalfOnly
logical :: distributionFunctionLowerHalfOnly , distributionFunctionNormalize
! Minimum mass to which subresolution fractions will be integrated.
double precision :: massMinimum
! Current epoch.
Expand Down Expand Up @@ -150,7 +150,8 @@ function generalizedPressSchechterConstructorParameters(parameters) result(self)
class (excursionSetFirstCrossingClass ), pointer :: excursionSetFirstCrossing_
class (mergerTreeBranchingProbabilityModifierClass ), pointer :: mergerTreeBranchingProbabilityModifier_
double precision :: deltaStepMaximum , massMinimum
logical :: smoothAccretion , distributionFunctionLowerHalfOnly
logical :: smoothAccretion , distributionFunctionLowerHalfOnly, &
& distributionFunctionNormalize

!![
<inputParameter>
Expand All @@ -177,13 +178,24 @@ function generalizedPressSchechterConstructorParameters(parameters) result(self)
<description>If true, only the lower half ($M &lt; M_0/2$) of the branching rate distribution function is used, as per the algorithm of \cite{cole_hierarchical_2000}.</description>
<source>parameters</source>
</inputParameter>
<inputParameter>
<name>distributionFunctionNormalize</name>
<defaultValue>.true.</defaultValue>
<description>
If using the full range ($M &lt; M_0$) of the branching rate distribution, if this parameter is {\normalfont \ttfamily
true} then divide the branching rate by 2. This is appropriate if two progenitors are to be sampled (i.e. a binary
split). If the branching rate applies to only a single branch is it more appropriate to set this parameter to be
{\normalfont \ttfamily true} in which case this normalization by a factor 2 is \emph{not} applied.
</description>
<source>parameters</source>
</inputParameter>
<objectBuilder class="criticalOverdensity" name="criticalOverdensity_" source="parameters"/>
<objectBuilder class="cosmologicalMassVariance" name="cosmologicalMassVariance_" source="parameters"/>
<objectBuilder class="cosmologyFunctions" name="cosmologyFunctions_" source="parameters"/>
<objectBuilder class="excursionSetFirstCrossing" name="excursionSetFirstCrossing_" source="parameters"/>
<objectBuilder class="mergerTreeBranchingProbabilityModifier" name="mergerTreeBranchingProbabilityModifier_" source="parameters"/>
!!]
self=mergerTreeBranchingProbabilityGnrlzdPrssSchchtr(deltaStepMaximum,massMinimum,smoothAccretion,distributionFunctionLowerHalfOnly,cosmologyFunctions_,criticalOverdensity_,cosmologicalMassVariance_,excursionSetFirstCrossing_,mergerTreeBranchingProbabilityModifier_)
self=mergerTreeBranchingProbabilityGnrlzdPrssSchchtr(deltaStepMaximum,massMinimum,smoothAccretion,distributionFunctionLowerHalfOnly,distributionFunctionNormalize,cosmologyFunctions_,criticalOverdensity_,cosmologicalMassVariance_,excursionSetFirstCrossing_,mergerTreeBranchingProbabilityModifier_)
!![
<inputParametersValidate source="parameters"/>
<objectDestructor name="criticalOverdensity_" />
Expand All @@ -195,7 +207,7 @@ <description>If true, only the lower half ($M &lt; M_0/2$) of the branching rate
return
end function generalizedPressSchechterConstructorParameters

function generalizedPressSchechterConstructorInternal(deltaStepMaximum,massMinimum,smoothAccretion,distributionFunctionLowerHalfOnly,cosmologyFunctions_,criticalOverdensity_,cosmologicalMassVariance_,excursionSetFirstCrossing_,mergerTreeBranchingProbabilityModifier_) result(self)
function generalizedPressSchechterConstructorInternal(deltaStepMaximum,massMinimum,smoothAccretion,distributionFunctionLowerHalfOnly,distributionFunctionNormalize,cosmologyFunctions_,criticalOverdensity_,cosmologicalMassVariance_,excursionSetFirstCrossing_,mergerTreeBranchingProbabilityModifier_) result(self)
!!{
Internal constructor for the \cite{cole_hierarchical_2000} merger tree building class.
!!}
Expand All @@ -208,10 +220,11 @@ function generalizedPressSchechterConstructorInternal(deltaStepMaximum,massMinim
class (excursionSetFirstCrossingClass ), intent(in ), target :: excursionSetFirstCrossing_
class (mergerTreeBranchingProbabilityModifierClass ), intent(in ), target :: mergerTreeBranchingProbabilityModifier_
double precision , intent(in ) :: deltaStepMaximum , massMinimum
logical , intent(in ) :: smoothAccretion , distributionFunctionLowerHalfOnly
logical , intent(in ) :: smoothAccretion , distributionFunctionLowerHalfOnly , &
& distributionFunctionNormalize
double precision , parameter :: toleranceAbsolute =0.0d+0, toleranceRelative =1.0d-9
!![
<constructorAssign variables="deltaStepMaximum, massMinimum, smoothAccretion, distributionFunctionLowerHalfOnly, *criticalOverdensity_, *cosmologicalMassVariance_, *cosmologyFunctions_, *excursionSetFirstCrossing_, *mergerTreeBranchingProbabilityModifier_"/>
<constructorAssign variables="deltaStepMaximum, massMinimum, smoothAccretion, distributionFunctionLowerHalfOnly, distributionFunctionNormalize, *criticalOverdensity_, *cosmologicalMassVariance_, *cosmologyFunctions_, *excursionSetFirstCrossing_, *mergerTreeBranchingProbabilityModifier_"/>
!!]

self%excursionSetsTested =.false.
Expand Down Expand Up @@ -307,12 +320,16 @@ double precision function generalizedPressSchechterMassBranch(self,haloMass,delt
call self%computeCommonFactors(node,haloMass,deltaCritical,time)
! Determine the upper mass limit to use.
if (self%distributionFunctionLowerHalfOnly) then
massUpper =+0.5d0*haloMass
self%normalization=+1.0d0
massUpper =+0.5d0*haloMass
self%normalization =+1.0d0
else
massUpper =+ haloMass
self%normalization=+0.5d0
end if
massUpper =+ haloMass
if (self%distributionFunctionNormalize) then
self%normalization=+0.5d0
else
self%normalization=+1.0d0
end if
end if
! Check that the root is bracketed.
rootFunctionLower=generalizedPressSchechterMassBranchRoot(massResolution)
rootFunctionUpper=generalizedPressSchechterMassBranchRoot(massUpper )
Expand Down Expand Up @@ -451,11 +468,15 @@ double precision function generalizedPressSchechterProbability(self,haloMass,del
call self%computeCommonFactors(node,haloMass,deltaCritical,time)
massMinimum=massResolution
if (self%distributionFunctionLowerHalfOnly) then
massMaximum =+0.5d0*self%parentHaloMass
normalization=+1.0d0
massMaximum =+0.5d0*self%parentHaloMass
normalization =+1.0d0
else
massMaximum =+ self%parentHaloMass
normalization=+0.5d0
massMaximum =+ self%parentHaloMass
if (self%distributionFunctionNormalize) then
normalization=+0.5d0
else
normalization=+1.0d0
end if
end if
self_ => self
generalizedPressSchechterProbability = +normalization &
Expand Down
2 changes: 1 addition & 1 deletion source/merger_trees.build.controller.constrained.F90
Original file line number Diff line number Diff line change
Expand Up @@ -110,7 +110,7 @@ function constrainedConstructorParameters(parameters) result(self)
<name>label</name>
<source>parameters</source>
<defaultValue>var_str(' ')</defaultValue>
<description>A label to apply to constrained node.</description>
<description>A label to apply to the constrained node.</description>
</inputParameter>
!!]
if (label == '') label=' '
Expand Down
157 changes: 157 additions & 0 deletions source/merger_trees.build.controller.filtered.F90
Original file line number Diff line number Diff line change
@@ -0,0 +1,157 @@
!! Copyright 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018,
!! 2019, 2020, 2021, 2022, 2023
!! Andrew Benson <[email protected]>
!!
!! This file is part of Galacticus.
!!
!! Galacticus is free software: you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! Galacticus is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with Galacticus. If not, see <http://www.gnu.org/licenses/>.

!!{
Contains a module which implements a merger tree build controller class which follows branches only if they pass a filter.
!!}

use :: Galactic_Filters, only : galacticFilterClass

!![
<mergerTreeBuildController name="mergerTreeBuildControllerFiltered">
<description>A merger tree build controller class which builds filtered trees.</description>
</mergerTreeBuildController>
!!]
type, extends(mergerTreeBuildControllerClass) :: mergerTreeBuildControllerFiltered
!!{
A merger tree build controller class which builds filtered trees. Specifically, only nodes which pass a provided filtered are
followed.
!!}
private
class(mergerTreeBranchingProbabilityClass), pointer :: mergerTreeBranchingProbability_ => null()
class(galacticFilterClass ), pointer :: galacticFilter_ => null()
contains
final :: filteredDestructor
procedure :: control => filteredControl
procedure :: branchingProbabilityObject => filteredBranchingProbabilityObject
procedure :: nodesInserted => filteredNodesInserted
end type mergerTreeBuildControllerFiltered

interface mergerTreeBuildControllerFiltered
!!{
Constructors for the ``filtered'' merger tree build controller class.
!!}
module procedure filteredConstructorParameters
module procedure filteredConstructorInternal
end interface mergerTreeBuildControllerFiltered

contains

function filteredConstructorParameters(parameters) result(self)
!!{
Constructor for the ``filtered'' merger tree build controller class which takes a parameter set as input.
!!}
use :: Input_Parameters, only : inputParameter, inputParameters
implicit none
type (mergerTreeBuildControllerFiltered ) :: self
type (inputParameters ), intent(inout) :: parameters
class(mergerTreeBranchingProbabilityClass), pointer :: mergerTreeBranchingProbability_
class(galacticFilterClass ), pointer :: galacticFilter_

!![
<objectBuilder class="mergerTreeBranchingProbability" name="mergerTreeBranchingProbability_" source="parameters"/>
<objectBuilder class="galacticFilter" name="galacticFilter_" source="parameters"/>
!!]
self=mergerTreeBuildControllerFiltered(mergerTreeBranchingProbability_,galacticFilter_)
!![
<inputParametersValidate source="parameters"/>
<objectDestructor name="mergerTreeBranchingProbability_"/>
!!]
return
end function filteredConstructorParameters

function filteredConstructorInternal(mergerTreeBranchingProbability_,galacticFilter_) result(self)
!!{
Internal constructor for the ``filtered'' merger tree build controller class .
!!}
implicit none
type (mergerTreeBuildControllerFiltered ) :: self
class(mergerTreeBranchingProbabilityClass), intent(in ), target :: mergerTreeBranchingProbability_
class(galacticFilterClass ), intent(in ), target :: galacticFilter_
!![
<constructorAssign variables="*mergerTreeBranchingProbability_, *galacticFilter_"/>
!!]

return
end function filteredConstructorInternal

subroutine filteredDestructor(self)
!!{
Destructor for the {\normalfont \ttfamily filtered} merger tree build controller class.
!!}
implicit none
type(mergerTreeBuildControllerFiltered), intent(inout) :: self

!![
<objectDestructor name="self%mergerTreeBranchingProbability_"/>
<objectDestructor name="self%galacticFilter_" />
!!]
return
end subroutine filteredDestructor

logical function filteredControl(self,node,treeWalker_)
!!{
Skip side branches of a tree under construction.
!!}
implicit none
class(mergerTreeBuildControllerFiltered), intent(inout) :: self
type (treeNode ), intent(inout), pointer :: node
class(mergerTreeWalkerClass ), intent(inout), optional :: treeWalker_
!$GLC attributes unused :: self

filteredControl=.true.
! Move to the next node in the tree while such exists, and the current node does not pass the filter.
do while (filteredControl.and..not.self%galacticFilter_%passes(node))
if (present(treeWalker_)) then
filteredControl=treeWalker_%next(node)
else
filteredControl=.false.
end if
end do
return
end function filteredControl

function filteredBranchingProbabilityObject(self,node) result(mergerTreeBranchingProbability_)
!!{
Return a pointer the the merger tree branching probability object to use.
!!}
implicit none
class(mergerTreeBranchingProbabilityClass), pointer :: mergerTreeBranchingProbability_
class(mergerTreeBuildControllerFiltered ), intent(inout) :: self
type (treeNode ), intent(inout) :: node
!$GLC attributes unused :: node

mergerTreeBranchingProbability_ => self%mergerTreeBranchingProbability_
return
end function filteredBranchingProbabilityObject

subroutine filteredNodesInserted(self,nodeCurrent,nodeProgenitor1,nodeProgenitor2,didBranch)
!!{
Act on the insertion of nodes into the merger tree.
!!}
implicit none
class (mergerTreeBuildControllerFiltered), intent(inout) :: self
type (treeNode ), intent(inout) :: nodeCurrent , nodeProgenitor1
type (treeNode ), intent(inout), optional :: nodeProgenitor2
logical , intent(in ), optional :: didBranch
!$GLC attributes unused :: self, nodeCurrent, nodeProgenitor1, nodeProgenitor2, didBranch

! Nothing to do.
return
end subroutine filteredNodesInserted
Loading

0 comments on commit 82f6406

Please sign in to comment.