Skip to content

Commit

Permalink
Merge pull request #743 from galacticusorg/featConstraintTools
Browse files Browse the repository at this point in the history 
Add new feature that are useful for MCMC (and similar) constraints work
  • Loading branch information
@abensonca
abensonca authored Nov 27, 2024
2 parents c5185af + 2b4513c commit ddb74ce
Show file tree
Hide file tree
Showing 4 changed files with 320 additions and 22 deletions.
2 changes: 1 addition & 1 deletion parameters/tutorials/mcmcBase.xml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -139,7 +139,7 @@
<redshiftInterval value="1" />
<computeScatter value="false" />
<systematicErrorPolynomialCoefficient value="0.0 0.0"/>
<likelihoodBin value="11" />
<likelihoodBins value="9" />
</outputAnalysis>

</parameters>
109 changes: 109 additions & 0 deletions source/merger_trees.filter.base_node.F90
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,109 @@
!! Copyright 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018,
!! 2019, 2020, 2021, 2022, 2023, 2024
!! 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 filter which passes if the base node in the tree passes the given galactic filter.
!!}

use :: Galactic_Filters, only : galacticFilterClass

!![
<mergerTreeFilter name="mergerTreeFilterBaseNode">
<description>A merger tree filter which passes if the base node in the tree passes the given galactic filter.</description>
</mergerTreeFilter>
!!]
type, extends(mergerTreeFilterClass) :: mergerTreeFilterBaseNode
!!{
A merger tree filter class which passes if the base node in the tree passes the given galactic filter.
!!}
private
class (galacticFilterClass), pointer :: galacticFilter_ => null()
contains
final :: baseNodeDestructor
procedure :: passes => baseNodePasses
end type mergerTreeFilterBaseNode

interface mergerTreeFilterBaseNode
!!{
Constructors for the ``baseNode'' merger tree filter class.
!!}
module procedure baseNodeConstructorParameters
module procedure baseNodeConstructorInternal
end interface mergerTreeFilterBaseNode

contains

function baseNodeConstructorParameters(parameters) result(self)
!!{
Constructor for the ``baseNode'' merger tree filter class which takes a parameter set as input.
!!}
use :: Input_Parameters, only : inputParameters
implicit none
type (mergerTreeFilterBaseNode) :: self
type (inputParameters ), intent(inout) :: parameters
class(galacticFilterClass ), pointer :: galacticFilter_
!![
<objectBuilder class="galacticFilter" name="galacticFilter_" source="parameters"/>
!!]
self=mergerTreeFilterBaseNode(galacticFilter_)
!![
<inputParametersValidate source="parameters"/>
<objectDestructor name="galacticFilter_"/>
!!]
return
end function baseNodeConstructorParameters

function baseNodeConstructorInternal(galacticFilter_) result(self)
!!{
Internal constructor for the ``baseNode'' merger tree filter class.
!!}
implicit none
type (mergerTreeFilterBaseNode) :: self
class(galacticFilterClass ), intent(in ), target :: galacticFilter_
!![
<constructorAssign variables="*galacticFilter_"/>
!!]

return
end function baseNodeConstructorInternal

subroutine baseNodeDestructor(self)
!!{
Destructor for the {\normalfont \ttfamily baseNode} merger tree filter class.
!!}
implicit none
type(mergerTreeFilterBaseNode), intent(inout) :: self

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

logical function baseNodePasses(self,tree) result(passes)
!!{
Implement a merger tree filter which passes if the base node in the tree passes the given merger tree filter.
!!}
implicit none
class(mergerTreeFilterBaseNode), intent(inout) :: self
type (mergerTree ), intent(in ) :: tree

passes=self%galacticFilter_%passes(tree%nodeBase)
return
end function baseNodePasses
170 changes: 170 additions & 0 deletions source/posterior_sampling.state.initialize.Gaussian_sphere.F90
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,170 @@
!! Copyright 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018,
!! 2019, 2020, 2021, 2022, 2023, 2024
!! 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/>.

!!{
Implementation of a posterior sampling state initializor class which initializes states using random draws from a Gaussian
sphere centered on the mode of the prior distribution.
!!}

!![
<posteriorSampleStateInitialize name="posteriorSampleStateInitializeGaussianSphere">
<description>
A posterior sampling state initialization class which initializes states using random draws from a Gaussian sphere centered on
the mode of the prior distribution.
</description>
</posteriorSampleStateInitialize>
!!]
type, extends(posteriorSampleStateInitializeClass) :: posteriorSampleStateInitializeGaussianSphere
!!{
Implementation of a posterior sampling state initialization class which initializes states using random draws from a Gaussian
sphere centered on the mode of the prior distribution.
!!}
private
class (randomNumberGeneratorClass), pointer :: randomNumberGenerator_ => null()
double precision :: radiusSphere
logical :: radiusIsRelative
contains
final :: gaussianSphereDestructor
procedure :: initialize => gaussianSphereInitialize
end type posteriorSampleStateInitializeGaussianSphere

interface posteriorSampleStateInitializeGaussianSphere
!!{
Constructors for the {\normalfont \ttfamily gaussianSphere} posterior sampling state initialization class.
!!}
module procedure gaussianSphereConstructorParameters
module procedure gaussianSphereConstructorInternal
end interface posteriorSampleStateInitializeGaussianSphere

contains

function gaussianSphereConstructorParameters(parameters) result(self)
!!{
Constructor for the {\normalfont \ttfamily gaussianSphere} posterior sampling state initialization class.
!!}
use :: Input_Parameters, only : inputParameters
implicit none
type (posteriorSampleStateInitializeGaussianSphere) :: self
type (inputParameters ), intent(inout) :: parameters
class (randomNumberGeneratorClass ), pointer :: randomNumberGenerator_
double precision :: radiusSphere
logical :: radiusIsRelative

!![
<inputParameter>
<name>radiusSphere</name>
<description>The radius of the Gaussian sphere.</description>
<source>parameters</source>
</inputParameter>
<inputParameter>
<name>radiusIsRelative</name>
<description>If true, the radius of the sphere is assumed to be relative to the extent of the prior, otherwise it is assumed to be an absolute radius.</description>
<source>parameters</source>
</inputParameter>
<objectBuilder class="randomNumberGenerator" name="randomNumberGenerator_" source="parameters"/>
!!]
self=posteriorSampleStateInitializeGaussianSphere(radiusSphere,radiusIsRelative,randomNumberGenerator_)
!![
<inputParametersValidate source="parameters"/>
<objectDestructor name="randomNumberGenerator_"/>
!!]
return
end function gaussianSphereConstructorParameters

function gaussianSphereConstructorInternal(radiusSphere,radiusIsRelative,randomNumberGenerator_) result(self)
!!{
Internal constructor for the {\normalfont \ttfamily gaussianSphere} posterior sampling state initialization class.
!!}
implicit none
type (posteriorSampleStateInitializeGaussianSphere) :: self
double precision , intent(in ) :: radiusSphere
logical , intent(in ) :: radiusIsRelative
class (randomNumberGeneratorClass ), intent(in ), target :: randomNumberGenerator_
!![
<constructorAssign variables="radiusSphere, radiusIsRelative, *randomNumberGenerator_"/>
!!]

return
end function gaussianSphereConstructorInternal

subroutine gaussianSphereDestructor(self)
!!{
Destructor for the {\normalfont \ttfamily gaussianSphere} posterior sampling state initialization class.
!!}
implicit none
type(posteriorSampleStateInitializeGaussianSphere), intent(inout) :: self

!![
<objectDestructor name="self%randomNumberGenerator_"/>
!!]
return
end subroutine gaussianSphereDestructor

subroutine gaussianSphereInitialize(self,simulationState,modelParameters_,modelLikelihood,timeEvaluatePrevious,logLikelihood,logPosterior)
!!{
Initialize simulation state by drawing at random from the parameter priors.
!!}
use :: Models_Likelihoods_Constants, only : logImpossible
implicit none
class (posteriorSampleStateInitializeGaussianSphere), intent(inout) :: self
class (posteriorSampleStateClass ), intent(inout) :: simulationState
class (posteriorSampleLikelihoodClass ), intent(inout) :: modelLikelihood
type (modelParameterList ), intent(inout), dimension(: ) :: modelParameters_
double precision , intent( out) :: timeEvaluatePrevious, logLikelihood , &
& logPosterior
double precision , dimension(size(modelParameters_)) :: state
double precision :: distributionMinimum , distributionMaximum, &
& distributionMedian , radius
integer :: j
logical :: first
!$GLC attributes unused :: self, modelLikelihood

! No knowledge of evaluation time.
timeEvaluatePrevious=-1.0d0
! We have no information about the likelihood of this state.
logLikelihood=logImpossible
logPosterior =logImpossible
! Initialize chain to some state vector.
state=0.0d0
do j=1,simulationState%dimension()
! Get the median, minimum, and maximum of the prior.
distributionMinimum=modelParameters_(j)%modelParameter_%map(modelParameters_(j)%modelParameter_%priorMinimum( ))
distributionMaximum=modelParameters_(j)%modelParameter_%map(modelParameters_(j)%modelParameter_%priorMaximum( ))
distributionMedian =modelParameters_(j)%modelParameter_%map(modelParameters_(j)%modelParameter_%priorInvert (0.5d0))
! Determine the radius of the Gaussian.
radius=self%radiusSphere
if (self%radiusIsRelative) radius=radius*(distributionMaximum-distributionMinimum)
! Sample values until a value within the allowed prior range is found.
first=.true.
do while ( &
& first &
& .or. &
& state(j) < distributionMinimum &
& .or. &
& state(j) > distributionMaximum &
& )
first =.false.
state(j)=+self%randomNumberGenerator_%standardNormalSample() &
& * radius &
& + distributionMedian
end do
end do
call simulationState%update(state,.false.,.false.)
return
end subroutine gaussianSphereInitialize
Loading

0 comments on commit ddb74ce

Please sign in to comment.