[oorb] r276 committed - Updates on MCMC ranging and random-walk ranging
2 views
Skip to first unread message
oo...@googlecode.com
Mar 11, 2015, 5:37:25 AM3/11/15
to oo...@googlegroups.com
Revision: 276
Author: jxvir...@gmail.com
Date: Wed Mar 11 09:37:07 2015 UTC
Log: Updates on MCMC ranging and random-walk ranging
https://code.google.com/p/oorb/source/detail?r=276
Modified:
/trunk/classes/StochasticOrbit_class.f90
=======================================
--- /trunk/classes/StochasticOrbit_class.f90 Tue Mar 10 07:57:31 2015 UTC
+++ /trunk/classes/StochasticOrbit_class.f90 Wed Mar 11 09:37:07 2015 UTC
@@ -4006,7 +4006,7 @@
ASSOCIATED(this%orb_arr_cmp)) THEN
indx_ml = MAXLOC(this%pdf_arr_cmp,dim=1)
getNominalOrbit = copy(this%orb_arr_cmp(indx_ml))
- print*,indx_ml
+ print*,indx_ml, this%pdf_arr_cmp(indx_ml)
ELSE
getNominalOrbit = copy(this%orb_ml_cmp)
END IF
@@ -5986,7 +5986,10 @@
+
SUBROUTINE autoMCMCRanging3(this)
+ !! Calls MCMCRanging (earlier: MCMCRanging3). Stable version.
+ !! NOTE: Combine with autoMCMCRanging (take parts from here e.g.
stdevs).
IMPLICIT NONE
TYPE(StochasticOrbit), INTENT(inout) :: this
@@ -6013,8 +6016,8 @@
storb%sor_ntrial_prm, "TRIALS..."
WRITE(stdout,"(A)") "***************"
END IF
- CALL MCMCRanging3(storb)
-! CALL MCMCRanging4(storb)
+ CALL MCMCRanging(storb)
+! CALL MCMCRanging3(storb)
IF (error) THEN
CALL errorMessage("StochasticOrbit / autoMCMCRanging3", &
"TRACE BACK (10)", 1)
@@ -6056,8 +6059,8 @@
storb%sor_norb_prm, "SAMPLE ORBITS..."
WRITE(stdout,"(A)") "***************"
END IF
- CALL MCMCRanging3(storb)
-! CALL MCMCRanging4(storb)
+ CALL MCMCRanging(storb)
+! CALL MCMCRanging3(storb)
IF (error) THEN
CALL errorMessage("StochasticOrbit / autoMCMCRanging3", &
"TRACE BACK (30)", 1)
@@ -6111,8 +6114,8 @@
this%sor_rho_prm(1,1:2) = sor_rho_init(1:2)
this%sor_rho_prm(2,1:2) = sor_rho_init(3:4)
- CALL MCMCRanging3(this)
- !CALL MCMCRanging4(storb)
+ CALL MCMCRanging(this)
+ !CALL MCMCRanging3(storb)
IF (error) THEN
CALL errorMessage("StochasticOrbit / autoMCMCRanging3", &
"TRACE BACK (50)", 1)
@@ -6163,10 +6166,11 @@
- SUBROUTINE autoMCMCRanging4(this)
- ! Automated version for either mcmc or random-walk
+ SUBROUTINE autoMCMCRanging(this)
+ ! Automated version for mcmc. Stable version (pre-2015 version, see
autoMCMCRanging2)
+ ! Calls MCMCRanging (earlier: MCMCRanging3).
! Using for initialization either
- ! - chains with uniform proposal
+ ! - no a priori orbit information: Ranging (i.e. uniform proposal
! - orbital distribution from input file
IMPLICIT NONE
@@ -6174,8 +6178,7 @@
TYPE(StochasticOrbit) :: storb
REAL(bp), DIMENSION(4) :: sor_rho_init
- ! Only initialize if not already done (in oorb.f90 with input file)
- IF (ALL(this%sor_iterate_bounds_prm)) THEN
+ IF
(.NOT.ASSOCIATED(this%orb_arr_cmp) .OR. .NOT.ASSOCIATED(this%pdf_arr_cmp))
THEN
this%sor_niter_cmp = 0
! First iteration with Ranging and uniform sampling
@@ -6191,7 +6194,7 @@
END IF
CALL statisticalRanging(storb)
IF (error .OR. storb%sor_norb_cmp <= 1) THEN
- CALL errorMessage("StochasticOrbit / autoMCMCRanging4", &
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging", &
"First iteration failed.", 1)
RETURN
END IF
@@ -6212,60 +6215,300 @@
END IF
CALL statisticalRanging(storb)
IF (error .OR. storb%sor_norb_cmp <= 1) THEN
- CALL errorMessage("StochasticOrbit / autoMCMCRanging4", &
- "First iteration failed.", 1)
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging", &
+ "Second iteration failed.", 1)
RETURN
END IF
this%sor_niter_cmp = 2
CALL updateRanging(storb, automatic=.TRUE.)
- sor_rho_init=getRangeBounds(storb)
- this%sor_rho_prm(1,1:2) = sor_rho_init(1:2)
- this%sor_rho_prm(2,1:2) = sor_rho_init(3:4)
- ! Use ml orbit as the first trial in the MCMC chain
- this%orb_ml_prm=copy(storb%orb_ml_cmp)
- call nullify(storb)
+ ELSE
+ storb = copy(this)
+ storb%sor_norb_prm = 50
+ storb%sor_ntrial_prm = this%sor_ntrial_sw_prm
+ ! Compute range distribution from orbit distribution
+ CALL constrainRangeDistributions(storb, storb%obss)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging", &
+ "TRACE BACK (15)", 1)
+ STOP
+ END IF
+ CALL setRangeBounds(storb)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging", &
+ "TRACE BACK (20)", 1)
+ STOP
+ END IF
+ storb%orb_ml_prm=getNominalOrbit(storb, ml_orbit=.true.)
+! storb%orb_ml_prm=copy(storb%orb_ml_cmp)
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(2X,A)") "***************"
+ WRITE(stdout,"(2X,A)") "FIRST ITERATION"
+ WRITE(stdout,"(2X,A)") "***************"
+ END IF
+ CALL MCMCRanging(storb)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging", &
+ "TRACE BACK (10)", 1)
+ STOP
+ END IF
+ this%sor_niter_cmp = 1
+
+ ! Is this needed, since burn in in MCMC discards orbits based on
their rms-values?
+ !CALL constrainRangeDistributions(storb, storb%obss)
+ !IF (error) THEN
+ ! CALL errorMessage("StochasticOrbit / autoMCMCRanging", &
+ ! "TRACE BACK (15)", 1)
+ ! STOP
+ !END IF
+ CALL setRangeBounds(storb)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging", &
+ "TRACE BACK (20)", 1)
+ STOP
+ END IF
+ storb%orb_ml_prm=copy(storb%orb_ml_cmp)
+ !storb%orb_ml_prm=getNominalOrbit(storb, ml_orbit=.true.)
+
+ !Second iterarion
+ storb%sor_norb_prm = this%sor_norb_sw_prm
+ storb%sor_ntrial_prm = this%sor_ntrial_sw_prm
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(1X)")
+ WRITE(stdout,"(2X,A)") "*****************"
+ WRITE(stdout,"(2X,A)") "SECOND ITERATION "
+ WRITE(stdout,"(2X,A)") "*****************"
+ END IF
+ CALL MCMCRanging(storb)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging", &
+ "TRACE BACK (10)", 1)
+ STOP
+ END IF
+ this%sor_niter_cmp = 2
+ CALL setRangeBounds(storb)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging", &
+ "TRACE BACK (20)", 1)
+ STOP
+ END IF
+
ENDIF
+ sor_rho_init=getRangeBounds(storb)
+ this%sor_rho_prm(1,1:2) = sor_rho_init(1:2)
+ this%sor_rho_prm(2,1:2) = sor_rho_init(3:4)
+ print*,sor_rho_init
+ ! Use ml orbit as the first trial in the MCMC chain
+ this%orb_ml_prm=copy(storb%orb_ml_cmp)
+ call nullify(storb)
! Start actual Markov chain
IF (info_verb >= 2) THEN
WRITE(stdout,"(2X,A)") "*****************"
- WRITE(stdout,"(A,1X,I0,1X,A)") "FINAL SAMPLING FOR",
this%sor_norb_prm, "SAMPLE ORBITS..."
+ WRITE(stdout,"(2X,A,1X,I0,1X,A)") "FINAL SAMPLING FOR",
this%sor_norb_prm, "SAMPLE ORBITS..."
WRITE(stdout,"(2X,A)") "*****************"
END IF
- CALL MCMCRanging3(this)
-
+ CALL MCMCRanging(this)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / autoMCMCRanging4", &
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging", &
"TRACE BACK (50)", 1)
STOP
END IF
+ this%sor_niter_cmp = 3
+ storb=copy(this)
+ CALL setRangeBounds(storb)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging", &
+ "TRACE BACK (20)", 1)
+ STOP
+ END IF
+ call nullify(storb)
- END SUBROUTINE autoMCMCRanging4
+ END SUBROUTINE autoMCMCRanging
+ SUBROUTINE autoRandomWalkRanging(this)
+ ! Automated version for random walk.
+ ! Calls randomWalkRanging.
+ ! Using for initialization either
+ ! - no a priori orbit information: Ranging (i.e. uniform proposal
+ ! - orbital distribution from input file
+ IMPLICIT NONE
+ TYPE(StochasticOrbit), INTENT(inout) :: this
+ TYPE(StochasticOrbit) :: storb
+ REAL(bp), DIMENSION(4) :: sor_rho_init
+
+ IF
(.NOT.ASSOCIATED(this%orb_arr_cmp) .OR. .NOT.ASSOCIATED(this%pdf_arr_cmp))
THEN
+
+ this%sor_niter_cmp = 0
+ ! First iteration with Ranging and uniform sampling
+ storb = copy(this)
+ storb%sor_norb_prm = 20
+ storb%sor_ntrial_prm = this%sor_ntrial_sw_prm
+ storb%dchi2_rejection_prm = .FALSE.
+
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(2X,A)") "***************"
+ WRITE(stdout,"(2X,A)") "FIRST ITERATION"
+ WRITE(stdout,"(2X,A)") "***************"
+ END IF
+ CALL statisticalRanging(storb)
+ IF (error .OR. storb%sor_norb_cmp <= 1) THEN
+ CALL errorMessage("StochasticOrbit / autoRandomWalkRanging", &
+ "First iteration failed.", 1)
+ RETURN
+ END IF
+ this%sor_niter_cmp = 1
+ CALL updateRanging(storb, automatic=.TRUE.)
+
+ !CALL NULLIFY(storb)
+
+ ! Second iteration
+ !storb = copy(this)
+ storb%sor_norb_prm = this%sor_norb_sw_prm
+ storb%sor_ntrial_prm = this%sor_ntrial_sw_prm
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(1X)")
+ WRITE(stdout,"(2X,A)") "*****************"
+ WRITE(stdout,"(2X,A)") "SECOND ITERATION "
+ WRITE(stdout,"(2X,A)") "*****************"
+ END IF
+ CALL statisticalRanging(storb)
+ IF (error .OR. storb%sor_norb_cmp <= 1) THEN
+ CALL errorMessage("StochasticOrbit / autoRandomWalkRanging", &
+ "Second iteration failed.", 1)
+ RETURN
+ END IF
+ this%sor_niter_cmp = 2
+ CALL updateRanging(storb, automatic=.TRUE.)
+
+ ELSE
+ storb = copy(this)
+ storb%sor_norb_prm = 50
+ storb%sor_ntrial_prm = this%sor_ntrial_sw_prm
+ ! Compute range distribution from orbit distribution
+ CALL constrainRangeDistributions(storb, storb%obss)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoRandomWalkRanging", &
+ "TRACE BACK (15)", 1)
+ STOP
+ END IF
+ CALL setRangeBounds(storb)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoRandomWalkRanging", &
+ "TRACE BACK (20)", 1)
+ STOP
+ END IF
+ storb%orb_ml_prm=getNominalOrbit(storb, ml_orbit=.true.)
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(2X,A)") "***************"
+ WRITE(stdout,"(2X,A)") "FIRST ITERATION"
+ WRITE(stdout,"(2X,A)") "***************"
+ END IF
+ CALL randomWalkRanging(storb)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoRandomWalkRanging", &
+ "TRACE BACK (10)", 1)
+ STOP
+ END IF
+ this%sor_niter_cmp = 1
+
+ ! Is this needed, since burn in in MCMC discards orbits based on
their rms-values?
+ !CALL constrainRangeDistributions(storb, storb%obss)
+ !IF (error) THEN
+ ! CALL errorMessage("StochasticOrbit / autoRandomWalkRanging", &
+ ! "TRACE BACK (15)", 1)
+ ! STOP
+ !END IF
+ CALL setRangeBounds(storb)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoRandomWalkRanging", &
+ "TRACE BACK (20)", 1)
+ STOP
+ END IF
+ if (storb%chi2_min_prm > storb%chi2_min_cmp) storb%chi2_min_prm =
storb%chi2_min_cmp
+ storb%orb_ml_prm=copy(storb%orb_ml_cmp)
+
+ !Second iterarion
+ storb%sor_norb_prm = this%sor_norb_sw_prm
+ storb%sor_ntrial_prm = this%sor_ntrial_sw_prm
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(1X)")
+ WRITE(stdout,"(2X,A)") "*****************"
+ WRITE(stdout,"(2X,A)") "SECOND ITERATION "
+ WRITE(stdout,"(2X,A)") "*****************"
+ END IF
+ CALL randomWalkRanging(storb)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoRandomWalkRanging", &
+ "TRACE BACK (10)", 1)
+ STOP
+ END IF
+ this%sor_niter_cmp = 2
+ CALL setRangeBounds(storb)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoRandomWalkRanging", &
+ "TRACE BACK (20)", 1)
+ STOP
+ END IF
+
+ ENDIF
+ sor_rho_init=getRangeBounds(storb)
+ this%sor_rho_prm(1,1:2) = sor_rho_init(1:2)
+ this%sor_rho_prm(2,1:2) = sor_rho_init(3:4)
+ ! Use ml orbit as the first trial in the MCMC chain
+ this%orb_ml_prm=copy(storb%orb_ml_cmp)
+ print*, 'Chi2min prm, cmp', storb%chi2_min_prm, storb%chi2_min_cmp
+ this%chi2_min_prm = storb%chi2_min_cmp
+ call nullify(storb)
+
+ ! Start actual Markov chain
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(2X,A)") "*****************"
+ WRITE(stdout,"(2X,A,1X,I0,1X,A)") "FINAL SAMPLING FOR",
this%sor_norb_prm, "SAMPLE ORBITS..."
+ WRITE(stdout,"(2X,A)") "*****************"
+ END IF
+ CALL randomWalkRanging(this)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoRandomWalkRanging", &
+ "TRACE BACK (50)", 1)
+ STOP
+ END IF
+ this%sor_niter_cmp = 3
+ print*, 'Chi2min prm, cmp', this%chi2_min_prm, this%chi2_min_cmp
+ storb=copy(this)
+ CALL setRangeBounds(storb)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoRandomWalkRanging", &
+ "TRACE BACK (20)", 1)
+ STOP
+ END IF
+ call nullify(storb)
+
+ END SUBROUTINE autoRandomWalkRanging
+
+
+
!! *Description:*
!!
- !! Markov Chain Monte Carlo ranging method.
- !! - Burn in phase: discard accepted orbits until (R.A., Dec.)
residual rms are reasonable.
- !! - Random-walk option (ran_walk) NOTE: parameter should be made
accessible from conf!
+ !! Random-walk ranging method.
+ !! - Burn in phase: discard accepted orbits until first chi2 in the
phase-space are
+ !! defined by dchi2.
+ !! - Random-walk option (ran_walk) can be turned off, which reverts
the routine to
+ !! a version of MCMC ranging
!!
!! Returns error.
!!
- SUBROUTINE MCMCRanging4(this)
+ SUBROUTINE randomWalkRanging(this)
IMPLICIT NONE
TYPE(StochasticOrbit), INTENT(inout) :: this
TYPE(Orbit) :: orb
- TYPE(SphericalCoordinates), DIMENSION(:), POINTER :: obs_scoords,
comp_scoords !observed
- !and
- !computed
- !sky
- !positions
+ TYPE(SphericalCoordinates), DIMENSION(:), POINTER :: obs_scoords,
comp_scoords
TYPE(SphericalCoordinates) :: obs_scoord1, obs_scoord2, obs_help
TYPE(CartesianCoordinates), DIMENSION(:), POINTER :: obsy_ccoords
TYPE(CartesianCoordinates) :: obs_ccoord_helio1, &
@@ -6294,7 +6537,7 @@
LOGICAL, DIMENSION(:,:), POINTER :: obs_masks
LOGICAL :: accepted, first = .TRUE., burn_in = .TRUE., frst_found
= .FALSE., ran_walk = .false., burnin_case=.false.
- ! SWITCH TO RANDOM WALK RANGING (-> move later to conf)
+ ! SWITCH TO RANDOM WALK RANGING (remnant: no need now when separated
from mcmc)
ran_walk = .true.
IF (info_verb >= 2) THEN
@@ -6304,7 +6547,7 @@
IF (.NOT. this%is_initialized_prm) THEN
error = .TRUE.
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"Object has not yet been initialized.", 1)
RETURN
END IF
@@ -6314,7 +6557,7 @@
DEALLOCATE(this%orb_arr_cmp, stat=err)
IF (err /= 0) THEN
error = .TRUE.
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"Could not deallocate memory (5).", 1)
RETURN
END IF
@@ -6322,7 +6565,7 @@
ALLOCATE(this%orb_arr_cmp(this%sor_norb_prm), stat=err)
IF (err /= 0) THEN
error = .TRUE.
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"Could not allocate memory (5).", 1)
RETURN
END IF
@@ -6340,7 +6583,7 @@
stat=err)
IF (err /= 0) THEN
error = .TRUE.
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"Could not allocate memory (10).", 1)
RETURN
END IF
@@ -6350,7 +6593,7 @@
DEALLOCATE(this%pdf_arr_cmp, stat=err)
IF (err /= 0) THEN
error = .TRUE.
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"Could not deallocate memory (10).", 1)
RETURN
END IF
@@ -6358,7 +6601,7 @@
ALLOCATE(this%pdf_arr_cmp(this%sor_norb_prm), stat=err)
IF (err /= 0) THEN
error = .TRUE.
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"Could not allocate memory (15).", 1)
RETURN
END IF
@@ -6368,7 +6611,7 @@
! Get observed sky positions (original observational data)
obs_scoords => getObservationSCoords(this%obss)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"TRACE BACK (5)", 1)
RETURN
END IF
@@ -6377,7 +6620,7 @@
CALL rotateToEquatorial(obs_scoords(i))
obs_coords(i,:) = getCoordinates(obs_scoords(i))
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"TRACE BACK (10)", 1)
CALL NULLIFY(obs_scoords(i))
RETURN
@@ -6388,7 +6631,7 @@
obs_masks => getObservationMasks(this%obss)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"TRACE BACK (15)", 1)
RETURN
END IF
@@ -6401,7 +6644,7 @@
information_matrix_obs => getBlockDiagInformationMatrix(this%obss)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"TRACE BACK (20)", 1)
DEALLOCATE(obsy_ccoords, stat=err)
DEALLOCATE(residuals, stat=err)
@@ -6413,7 +6656,7 @@
obsy_ccoords => getObservatoryCCoords(this%obss)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"TRACE BACK (25)", 1)
RETURN
END IF
@@ -6423,7 +6666,7 @@
obs_pair = RESHAPE(this%sor_pair_arr_prm, (/ 2 /))
ELSE
error = .TRUE.
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"Use of random observation pairs not accepted.", 1)
RETURN
END IF
@@ -6431,7 +6674,7 @@
obs_scoord1 = copy(obs_scoords(obs_pair(1)))
obs_scoord2 = copy(obs_scoords(obs_pair(2)))
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"TRACE BACK (30)", 1)
RETURN
END IF
@@ -6458,19 +6701,24 @@
state_(1:3) = getPosition(obs_scoord1)
state_(4:6) = getPosition(obs_scoord2)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"TRACE BACK (50)", 1)
DEALLOCATE(obs_masks, stat=err)
RETURN
END IF
! Center the first trial to the ml orbit if available
- !if (exist(this%orb_ml_cmp)) then
- !CALL getEphemerides(this%orb_ml_cmp, obsy_ccoords, comp_scoords)
- !state_(1)=comp_scoords(obs_pair(1),1)
- !state_(4)=comp_scoords(obs_pair(2),1)
- ! else
- state_(1) = (this%sor_rho_prm(1,2) + this%sor_rho_prm(1,1))/2.0_bp
- state_(4) = state_(1) + (this%sor_rho_prm(2,2) +
this%sor_rho_prm(2,1))/2.0_bp
+ IF (exist(this%orb_ml_prm)) then
+ CALL getEphemerides(this%orb_ml_prm, obsy_ccoords, comp_scoords)
+ state_(1)=getDistance(comp_scoords(obs_pair(1)))
+ state_(4)=getDistance(comp_scoords(obs_pair(2)))
+ print*,state_(1), state_(4)
+ else
+ state_(1) = (this%sor_rho_prm(1,2) + this%sor_rho_prm(1,1))/2.0_bp
+ state_(4) = state_(1) + (this%sor_rho_prm(2,2) +
this%sor_rho_prm(2,1))/2.0_bp
+ endif
+ print*,"...centered around ", state_(1), state_(4)
+ !state_(1) = (this%sor_rho_prm(1,2) + this%sor_rho_prm(1,1))/2.0_bp
+ !state_(4) = state_(1) + (this%sor_rho_prm(2,2) +
this%sor_rho_prm(2,1))/2.0_bp
! set counter some parameters
this%sor_norb_cmp = 0
@@ -6523,7 +6771,7 @@
CALL NEW(obs_scoord1, state(1), state(2), state(3),
getTime(obsy_ccoords(obs_pair(1))))
CALL NEW(obs_scoord2, state(4), state(5), state(6),
getTime(obsy_ccoords(obs_pair(2))))
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"TRACE BACK (65)", 1)
RETURN
END IF
@@ -6534,7 +6782,7 @@
CALL NEW(obs_ccoord_topo1, obs_scoord1)
CALL NEW(obs_ccoord_topo2, obs_scoord2)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"TRACE BACK (70)", 1)
RETURN
END IF
@@ -6545,7 +6793,7 @@
CALL rotateToEcliptic(obsy_ccoords(obs_pair(1)))
CALL rotateToEcliptic(obsy_ccoords(obs_pair(2)))
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"TRACE BACK (75)", 1)
RETURN
END IF
@@ -6556,7 +6804,7 @@
obs_ccoord_helio1 = copy(obsy_ccoords(obs_pair(1)) +
obs_ccoord_topo1)
obs_ccoord_helio2 = copy(obsy_ccoords(obs_pair(2)) +
obs_ccoord_topo2)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"TRACE BACK (80)", 1)
RETURN
END IF
@@ -6564,7 +6812,7 @@
CALL estimateLightTime(obs_ccoord_helio1, state(1))! changed from 3
and 6
CALL estimateLightTime(obs_ccoord_helio2, state(1)+state(4))
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"TRACE BACK (85)", 1)
RETURN
END IF
@@ -6580,7 +6828,7 @@
integrator=this%integrator_prm, &
integration_step=this%integration_step_prm)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"TRACE BACK (86)", 1)
error = .FALSE.
CYCLE
@@ -6612,7 +6860,7 @@
END IF
END IF
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"TRACE BACK (xx)", 1)
error = .FALSE.
CYCLE
@@ -6687,7 +6935,7 @@
IF (.NOT. equal(this%t_inv_prm,getTime(orb))) THEN
CALL propagate(orb, this%t_inv_prm)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"TRACE BACK (90)", 1)
error = .FALSE.
CYCLE
@@ -6706,7 +6954,7 @@
! get computed sky positions - spherical coordinates
CALL getEphemerides(orb, obsy_ccoords, comp_scoords,
partials_arr=partials_arr)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"TRACE BACK (45)", 1)
DEALLOCATE(comp_scoords, stat=err)
DEALLOCATE(partials_arr, stat=err)
@@ -6727,7 +6975,7 @@
CALL rotateToEquatorial(comp_scoords(i))
comp_coords(i,:) = getCoordinates(comp_scoords(i))
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"TRACE BACK (55)", 1)
DEALLOCATE(obsy_ccoords, stat=err)
DEALLOCATE(obs_masks, stat=err)
@@ -6759,7 +7007,7 @@
chi2 = chi_square(residuals, information_matrix_obs, errstr=errstr)
IF (LEN_TRIM(errstr) /= 0) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"TRACE BACK (60)", 1)
DEALLOCATE(obs_masks, stat=err)
RETURN
@@ -6777,7 +7025,7 @@
DEALLOCATE(partials_arr, stat=err)
jac_sph_inv = ABS(determinant(jacobian_matrix, errstr))
IF (LEN_TRIM(errstr) /= 0) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"Unsuccessful computation of determinant of orbital
element " // &
"jacobian matrix: " // TRIM(errstr), 1)
CALL matrix_print(jacobian_matrix, stderr, errstr)
@@ -6919,7 +7167,7 @@
! From ranging:
! "- generated distance to object at second epoch (NB!
Relative
! to the distance at the first epoch.)"
- this%sor_rho_arr_cmp(this%sor_norb_cmp,2) = state(1) -
state(4)
+ this%sor_rho_arr_cmp(this%sor_norb_cmp,2) = state(4) ! -
state(4)
! ELSE IF (count_dchi2 > this%count_dchi2_prm) then
! burn_in = .FALSE.
END IF
@@ -6971,6 +7219,7 @@
CALL propagate(this%orb_arr_cmp, this%t_inv_prm)
i = MAXLOC(this%pdf_arr_cmp,dim=1)
this%orb_ml_cmp = copy(this%orb_arr_cmp(i))
+ this%chi2_min_cmp = chi2min
IF (info_verb >= 2) THEN
WRITE(stdout,"(3(A,1X))") "Element types for covariance matrix and
ML orbit:", &
this%cov_type_prm, getElementType(this%orb_ml_cmp)
@@ -6982,7 +7231,7 @@
END IF
ELSE
error = .TRUE.
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"No sample orbits found!", 1)
! RETURN
END IF
@@ -6996,12 +7245,12 @@
information_matrix_obs, stat=err)
IF (err /= 0) THEN
error = .TRUE.
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"Could not deallocate memory (20).", 1)
RETURN
END IF
- END SUBROUTINE MCMCRanging4
+ END SUBROUTINE randomWalkRanging
@@ -7009,12 +7258,13 @@
!! *Description:*
!!
- !! Markov Chain Monte Carlo ranging method.
+ !! Markov Chain Monte Carlo ranging method.
+ !! - Currently stable version! (Pre-2015 version, see MCMCRanging2)
!! - Burn in phase: discard accepted orbits until (R.A., Dec.)
residual rms are reasonable.
!!
!! Returns error.
!!
- SUBROUTINE MCMCRanging3(this)
+ SUBROUTINE MCMCRanging(this)
IMPLICIT NONE
TYPE(StochasticOrbit), INTENT(inout) :: this
@@ -7057,7 +7307,7 @@
IF (.NOT. this%is_initialized_prm) THEN
error = .TRUE.
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"Object has not yet been initialized.", 1)
RETURN
END IF
@@ -7067,7 +7317,7 @@
DEALLOCATE(this%orb_arr_cmp, stat=err)
IF (err /= 0) THEN
error = .TRUE.
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"Could not deallocate memory (5).", 1)
RETURN
END IF
@@ -7075,7 +7325,7 @@
ALLOCATE(this%orb_arr_cmp(this%sor_norb_prm), stat=err)
IF (err /= 0) THEN
error = .TRUE.
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"Could not allocate memory (5).", 1)
RETURN
END IF
@@ -7093,7 +7343,7 @@
stat=err)
IF (err /= 0) THEN
error = .TRUE.
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"Could not allocate memory (10).", 1)
RETURN
END IF
@@ -7103,7 +7353,7 @@
DEALLOCATE(this%pdf_arr_cmp, stat=err)
IF (err /= 0) THEN
error = .TRUE.
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"Could not deallocate memory (10).", 1)
RETURN
END IF
@@ -7111,7 +7361,7 @@
ALLOCATE(this%pdf_arr_cmp(this%sor_norb_prm), stat=err)
IF (err /= 0) THEN
error = .TRUE.
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"Could not allocate memory (15).", 1)
RETURN
END IF
@@ -7121,7 +7371,7 @@
! Get observed sky positions (original observational data)
obs_scoords => getObservationSCoords(this%obss)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"TRACE BACK (5)", 1)
RETURN
END IF
@@ -7130,7 +7380,7 @@
CALL rotateToEquatorial(obs_scoords(i))
obs_coords(i,:) = getCoordinates(obs_scoords(i))
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"TRACE BACK (10)", 1)
CALL NULLIFY(obs_scoords(i))
RETURN
@@ -7141,7 +7391,7 @@
obs_masks => getObservationMasks(this%obss)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"TRACE BACK (15)", 1)
RETURN
END IF
@@ -7154,7 +7404,7 @@
information_matrix_obs => getBlockDiagInformationMatrix(this%obss)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"TRACE BACK (20)", 1)
DEALLOCATE(obsy_ccoords, stat=err)
DEALLOCATE(residuals, stat=err)
@@ -7166,7 +7416,7 @@
obsy_ccoords => getObservatoryCCoords(this%obss)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"TRACE BACK (25)", 1)
RETURN
END IF
@@ -7176,7 +7426,7 @@
obs_pair = RESHAPE(this%sor_pair_arr_prm, (/ 2 /))
ELSE
error = .TRUE.
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"Use of random observation pairs not accepted.", 1)
RETURN
END IF
@@ -7184,7 +7434,7 @@
obs_scoord1 = copy(obs_scoords(obs_pair(1)))
obs_scoord2 = copy(obs_scoords(obs_pair(2)))
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"TRACE BACK (30)", 1)
RETURN
END IF
@@ -7211,12 +7461,12 @@
state_(1:3) = getPosition(obs_scoord1)
state_(4:6) = getPosition(obs_scoord2)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"TRACE BACK (50)", 1)
DEALLOCATE(obs_masks, stat=err)
RETURN
END IF
- ! Center the first trial to the ml orbit if available -> this could be
taken outside mcmcranging?
+ ! Center the first trial to the ml orbit if available
IF (exist(this%orb_ml_prm)) then
CALL getEphemerides(this%orb_ml_prm, obsy_ccoords, comp_scoords)
state_(1)=getDistance(comp_scoords(obs_pair(1)))
@@ -7279,7 +7529,7 @@
CALL NEW(obs_scoord1, state(1), state(2), state(3),
getTime(obsy_ccoords(obs_pair(1))))
CALL NEW(obs_scoord2, state(4), state(5), state(6),
getTime(obsy_ccoords(obs_pair(2))))
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"TRACE BACK (65)", 1)
RETURN
END IF
@@ -7290,7 +7540,7 @@
CALL NEW(obs_ccoord_topo1, obs_scoord1)
CALL NEW(obs_ccoord_topo2, obs_scoord2)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"TRACE BACK (70)", 1)
RETURN
END IF
@@ -7301,7 +7551,7 @@
CALL rotateToEcliptic(obsy_ccoords(obs_pair(1)))
CALL rotateToEcliptic(obsy_ccoords(obs_pair(2)))
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"TRACE BACK (75)", 1)
RETURN
END IF
@@ -7312,7 +7562,7 @@
obs_ccoord_helio1 = copy(obsy_ccoords(obs_pair(1)) +
obs_ccoord_topo1)
obs_ccoord_helio2 = copy(obsy_ccoords(obs_pair(2)) +
obs_ccoord_topo2)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"TRACE BACK (80)", 1)
RETURN
END IF
@@ -7320,7 +7570,7 @@
CALL estimateLightTime(obs_ccoord_helio1, state(1))! changed from 3
and 6
CALL estimateLightTime(obs_ccoord_helio2, state(1)+state(4))
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"TRACE BACK (85)", 1)
RETURN
END IF
@@ -7336,7 +7586,7 @@
integrator=this%integrator_prm, &
integration_step=this%integration_step_prm)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"TRACE BACK (86)", 1)
error = .FALSE.
CYCLE
@@ -7437,7 +7687,7 @@
IF (.NOT. equal(this%t_inv_prm,getTime(orb))) THEN
CALL propagate(orb, this%t_inv_prm)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"TRACE BACK (90)", 1)
error = .FALSE.
CYCLE
@@ -7456,7 +7706,7 @@
! get computed sky positions - spherical coordinates
CALL getEphemerides(orb, obsy_ccoords, comp_scoords,
partials_arr=partials_arr)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"TRACE BACK (45)", 1)
DEALLOCATE(comp_scoords, stat=err)
DEALLOCATE(partials_arr, stat=err)
@@ -7477,7 +7727,7 @@
CALL rotateToEquatorial(comp_scoords(i))
comp_coords(i,:) = getCoordinates(comp_scoords(i))
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"TRACE BACK (55)", 1)
DEALLOCATE(obsy_ccoords, stat=err)
DEALLOCATE(obs_masks, stat=err)
@@ -7509,7 +7759,7 @@
chi2 = chi_square(residuals, information_matrix_obs, errstr=errstr)
IF (LEN_TRIM(errstr) /= 0) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"TRACE BACK (60)", 1)
DEALLOCATE(obs_masks, stat=err)
RETURN
@@ -7527,7 +7777,7 @@
DEALLOCATE(partials_arr, stat=err)
jac_sph_inv = ABS(determinant(jacobian_matrix, errstr))
IF (LEN_TRIM(errstr) /= 0) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"Unsuccessful computation of determinant of orbital
element " // &
"jacobian matrix: " // TRIM(errstr), 1)
CALL matrix_print(jacobian_matrix, stderr, errstr)
@@ -7646,7 +7896,7 @@
END IF
- IF (this%sor_norb_cmp > 1) THEN
+ IF (this%sor_norb_cmp >= 1) THEN
this%repetition_arr_cmp(this%sor_norb_cmp) =
this%repetition_arr_cmp(this%sor_norb_cmp) + 1
!print*,this%sor_norb_cmp,
this%repetition_arr_cmp(this%sor_norb_cmp)
END IF
@@ -7692,7 +7942,7 @@
END IF
ELSE
error = .TRUE.
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"No sample orbits found!", 1)
! RETURN
END IF
@@ -7706,17 +7956,19 @@
information_matrix_obs, stat=err)
IF (err /= 0) THEN
error = .TRUE.
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"Could not deallocate memory (20).", 1)
RETURN
END IF
- END SUBROUTINE MCMCRanging3
+ END SUBROUTINE MCMCRanging
+
+
***The diff for this file has been truncated for email.***
Author: jxvir...@gmail.com
Date: Wed Mar 11 09:37:07 2015 UTC
Log: Updates on MCMC ranging and random-walk ranging
https://code.google.com/p/oorb/source/detail?r=276
Modified:
/trunk/classes/StochasticOrbit_class.f90
=======================================
--- /trunk/classes/StochasticOrbit_class.f90 Tue Mar 10 07:57:31 2015 UTC
+++ /trunk/classes/StochasticOrbit_class.f90 Wed Mar 11 09:37:07 2015 UTC
@@ -4006,7 +4006,7 @@
ASSOCIATED(this%orb_arr_cmp)) THEN
indx_ml = MAXLOC(this%pdf_arr_cmp,dim=1)
getNominalOrbit = copy(this%orb_arr_cmp(indx_ml))
- print*,indx_ml
+ print*,indx_ml, this%pdf_arr_cmp(indx_ml)
ELSE
getNominalOrbit = copy(this%orb_ml_cmp)
END IF
@@ -5986,7 +5986,10 @@
+
SUBROUTINE autoMCMCRanging3(this)
+ !! Calls MCMCRanging (earlier: MCMCRanging3). Stable version.
+ !! NOTE: Combine with autoMCMCRanging (take parts from here e.g.
stdevs).
IMPLICIT NONE
TYPE(StochasticOrbit), INTENT(inout) :: this
@@ -6013,8 +6016,8 @@
storb%sor_ntrial_prm, "TRIALS..."
WRITE(stdout,"(A)") "***************"
END IF
- CALL MCMCRanging3(storb)
-! CALL MCMCRanging4(storb)
+ CALL MCMCRanging(storb)
+! CALL MCMCRanging3(storb)
IF (error) THEN
CALL errorMessage("StochasticOrbit / autoMCMCRanging3", &
"TRACE BACK (10)", 1)
@@ -6056,8 +6059,8 @@
storb%sor_norb_prm, "SAMPLE ORBITS..."
WRITE(stdout,"(A)") "***************"
END IF
- CALL MCMCRanging3(storb)
-! CALL MCMCRanging4(storb)
+ CALL MCMCRanging(storb)
+! CALL MCMCRanging3(storb)
IF (error) THEN
CALL errorMessage("StochasticOrbit / autoMCMCRanging3", &
"TRACE BACK (30)", 1)
@@ -6111,8 +6114,8 @@
this%sor_rho_prm(1,1:2) = sor_rho_init(1:2)
this%sor_rho_prm(2,1:2) = sor_rho_init(3:4)
- CALL MCMCRanging3(this)
- !CALL MCMCRanging4(storb)
+ CALL MCMCRanging(this)
+ !CALL MCMCRanging3(storb)
IF (error) THEN
CALL errorMessage("StochasticOrbit / autoMCMCRanging3", &
"TRACE BACK (50)", 1)
@@ -6163,10 +6166,11 @@
- SUBROUTINE autoMCMCRanging4(this)
- ! Automated version for either mcmc or random-walk
+ SUBROUTINE autoMCMCRanging(this)
+ ! Automated version for mcmc. Stable version (pre-2015 version, see
autoMCMCRanging2)
+ ! Calls MCMCRanging (earlier: MCMCRanging3).
! Using for initialization either
- ! - chains with uniform proposal
+ ! - no a priori orbit information: Ranging (i.e. uniform proposal
! - orbital distribution from input file
IMPLICIT NONE
@@ -6174,8 +6178,7 @@
TYPE(StochasticOrbit) :: storb
REAL(bp), DIMENSION(4) :: sor_rho_init
- ! Only initialize if not already done (in oorb.f90 with input file)
- IF (ALL(this%sor_iterate_bounds_prm)) THEN
+ IF
(.NOT.ASSOCIATED(this%orb_arr_cmp) .OR. .NOT.ASSOCIATED(this%pdf_arr_cmp))
THEN
this%sor_niter_cmp = 0
! First iteration with Ranging and uniform sampling
@@ -6191,7 +6194,7 @@
END IF
CALL statisticalRanging(storb)
IF (error .OR. storb%sor_norb_cmp <= 1) THEN
- CALL errorMessage("StochasticOrbit / autoMCMCRanging4", &
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging", &
"First iteration failed.", 1)
RETURN
END IF
@@ -6212,60 +6215,300 @@
END IF
CALL statisticalRanging(storb)
IF (error .OR. storb%sor_norb_cmp <= 1) THEN
- CALL errorMessage("StochasticOrbit / autoMCMCRanging4", &
- "First iteration failed.", 1)
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging", &
+ "Second iteration failed.", 1)
RETURN
END IF
this%sor_niter_cmp = 2
CALL updateRanging(storb, automatic=.TRUE.)
- sor_rho_init=getRangeBounds(storb)
- this%sor_rho_prm(1,1:2) = sor_rho_init(1:2)
- this%sor_rho_prm(2,1:2) = sor_rho_init(3:4)
- ! Use ml orbit as the first trial in the MCMC chain
- this%orb_ml_prm=copy(storb%orb_ml_cmp)
- call nullify(storb)
+ ELSE
+ storb = copy(this)
+ storb%sor_norb_prm = 50
+ storb%sor_ntrial_prm = this%sor_ntrial_sw_prm
+ ! Compute range distribution from orbit distribution
+ CALL constrainRangeDistributions(storb, storb%obss)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging", &
+ "TRACE BACK (15)", 1)
+ STOP
+ END IF
+ CALL setRangeBounds(storb)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging", &
+ "TRACE BACK (20)", 1)
+ STOP
+ END IF
+ storb%orb_ml_prm=getNominalOrbit(storb, ml_orbit=.true.)
+! storb%orb_ml_prm=copy(storb%orb_ml_cmp)
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(2X,A)") "***************"
+ WRITE(stdout,"(2X,A)") "FIRST ITERATION"
+ WRITE(stdout,"(2X,A)") "***************"
+ END IF
+ CALL MCMCRanging(storb)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging", &
+ "TRACE BACK (10)", 1)
+ STOP
+ END IF
+ this%sor_niter_cmp = 1
+
+ ! Is this needed, since burn in in MCMC discards orbits based on
their rms-values?
+ !CALL constrainRangeDistributions(storb, storb%obss)
+ !IF (error) THEN
+ ! CALL errorMessage("StochasticOrbit / autoMCMCRanging", &
+ ! "TRACE BACK (15)", 1)
+ ! STOP
+ !END IF
+ CALL setRangeBounds(storb)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging", &
+ "TRACE BACK (20)", 1)
+ STOP
+ END IF
+ storb%orb_ml_prm=copy(storb%orb_ml_cmp)
+ !storb%orb_ml_prm=getNominalOrbit(storb, ml_orbit=.true.)
+
+ !Second iterarion
+ storb%sor_norb_prm = this%sor_norb_sw_prm
+ storb%sor_ntrial_prm = this%sor_ntrial_sw_prm
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(1X)")
+ WRITE(stdout,"(2X,A)") "*****************"
+ WRITE(stdout,"(2X,A)") "SECOND ITERATION "
+ WRITE(stdout,"(2X,A)") "*****************"
+ END IF
+ CALL MCMCRanging(storb)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging", &
+ "TRACE BACK (10)", 1)
+ STOP
+ END IF
+ this%sor_niter_cmp = 2
+ CALL setRangeBounds(storb)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging", &
+ "TRACE BACK (20)", 1)
+ STOP
+ END IF
+
ENDIF
+ sor_rho_init=getRangeBounds(storb)
+ this%sor_rho_prm(1,1:2) = sor_rho_init(1:2)
+ this%sor_rho_prm(2,1:2) = sor_rho_init(3:4)
+ print*,sor_rho_init
+ ! Use ml orbit as the first trial in the MCMC chain
+ this%orb_ml_prm=copy(storb%orb_ml_cmp)
+ call nullify(storb)
! Start actual Markov chain
IF (info_verb >= 2) THEN
WRITE(stdout,"(2X,A)") "*****************"
- WRITE(stdout,"(A,1X,I0,1X,A)") "FINAL SAMPLING FOR",
this%sor_norb_prm, "SAMPLE ORBITS..."
+ WRITE(stdout,"(2X,A,1X,I0,1X,A)") "FINAL SAMPLING FOR",
this%sor_norb_prm, "SAMPLE ORBITS..."
WRITE(stdout,"(2X,A)") "*****************"
END IF
- CALL MCMCRanging3(this)
-
+ CALL MCMCRanging(this)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / autoMCMCRanging4", &
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging", &
"TRACE BACK (50)", 1)
STOP
END IF
+ this%sor_niter_cmp = 3
+ storb=copy(this)
+ CALL setRangeBounds(storb)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging", &
+ "TRACE BACK (20)", 1)
+ STOP
+ END IF
+ call nullify(storb)
- END SUBROUTINE autoMCMCRanging4
+ END SUBROUTINE autoMCMCRanging
+ SUBROUTINE autoRandomWalkRanging(this)
+ ! Automated version for random walk.
+ ! Calls randomWalkRanging.
+ ! Using for initialization either
+ ! - no a priori orbit information: Ranging (i.e. uniform proposal
+ ! - orbital distribution from input file
+ IMPLICIT NONE
+ TYPE(StochasticOrbit), INTENT(inout) :: this
+ TYPE(StochasticOrbit) :: storb
+ REAL(bp), DIMENSION(4) :: sor_rho_init
+
+ IF
(.NOT.ASSOCIATED(this%orb_arr_cmp) .OR. .NOT.ASSOCIATED(this%pdf_arr_cmp))
THEN
+
+ this%sor_niter_cmp = 0
+ ! First iteration with Ranging and uniform sampling
+ storb = copy(this)
+ storb%sor_norb_prm = 20
+ storb%sor_ntrial_prm = this%sor_ntrial_sw_prm
+ storb%dchi2_rejection_prm = .FALSE.
+
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(2X,A)") "***************"
+ WRITE(stdout,"(2X,A)") "FIRST ITERATION"
+ WRITE(stdout,"(2X,A)") "***************"
+ END IF
+ CALL statisticalRanging(storb)
+ IF (error .OR. storb%sor_norb_cmp <= 1) THEN
+ CALL errorMessage("StochasticOrbit / autoRandomWalkRanging", &
+ "First iteration failed.", 1)
+ RETURN
+ END IF
+ this%sor_niter_cmp = 1
+ CALL updateRanging(storb, automatic=.TRUE.)
+
+ !CALL NULLIFY(storb)
+
+ ! Second iteration
+ !storb = copy(this)
+ storb%sor_norb_prm = this%sor_norb_sw_prm
+ storb%sor_ntrial_prm = this%sor_ntrial_sw_prm
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(1X)")
+ WRITE(stdout,"(2X,A)") "*****************"
+ WRITE(stdout,"(2X,A)") "SECOND ITERATION "
+ WRITE(stdout,"(2X,A)") "*****************"
+ END IF
+ CALL statisticalRanging(storb)
+ IF (error .OR. storb%sor_norb_cmp <= 1) THEN
+ CALL errorMessage("StochasticOrbit / autoRandomWalkRanging", &
+ "Second iteration failed.", 1)
+ RETURN
+ END IF
+ this%sor_niter_cmp = 2
+ CALL updateRanging(storb, automatic=.TRUE.)
+
+ ELSE
+ storb = copy(this)
+ storb%sor_norb_prm = 50
+ storb%sor_ntrial_prm = this%sor_ntrial_sw_prm
+ ! Compute range distribution from orbit distribution
+ CALL constrainRangeDistributions(storb, storb%obss)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoRandomWalkRanging", &
+ "TRACE BACK (15)", 1)
+ STOP
+ END IF
+ CALL setRangeBounds(storb)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoRandomWalkRanging", &
+ "TRACE BACK (20)", 1)
+ STOP
+ END IF
+ storb%orb_ml_prm=getNominalOrbit(storb, ml_orbit=.true.)
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(2X,A)") "***************"
+ WRITE(stdout,"(2X,A)") "FIRST ITERATION"
+ WRITE(stdout,"(2X,A)") "***************"
+ END IF
+ CALL randomWalkRanging(storb)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoRandomWalkRanging", &
+ "TRACE BACK (10)", 1)
+ STOP
+ END IF
+ this%sor_niter_cmp = 1
+
+ ! Is this needed, since burn in in MCMC discards orbits based on
their rms-values?
+ !CALL constrainRangeDistributions(storb, storb%obss)
+ !IF (error) THEN
+ ! CALL errorMessage("StochasticOrbit / autoRandomWalkRanging", &
+ ! "TRACE BACK (15)", 1)
+ ! STOP
+ !END IF
+ CALL setRangeBounds(storb)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoRandomWalkRanging", &
+ "TRACE BACK (20)", 1)
+ STOP
+ END IF
+ if (storb%chi2_min_prm > storb%chi2_min_cmp) storb%chi2_min_prm =
storb%chi2_min_cmp
+ storb%orb_ml_prm=copy(storb%orb_ml_cmp)
+
+ !Second iterarion
+ storb%sor_norb_prm = this%sor_norb_sw_prm
+ storb%sor_ntrial_prm = this%sor_ntrial_sw_prm
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(1X)")
+ WRITE(stdout,"(2X,A)") "*****************"
+ WRITE(stdout,"(2X,A)") "SECOND ITERATION "
+ WRITE(stdout,"(2X,A)") "*****************"
+ END IF
+ CALL randomWalkRanging(storb)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoRandomWalkRanging", &
+ "TRACE BACK (10)", 1)
+ STOP
+ END IF
+ this%sor_niter_cmp = 2
+ CALL setRangeBounds(storb)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoRandomWalkRanging", &
+ "TRACE BACK (20)", 1)
+ STOP
+ END IF
+
+ ENDIF
+ sor_rho_init=getRangeBounds(storb)
+ this%sor_rho_prm(1,1:2) = sor_rho_init(1:2)
+ this%sor_rho_prm(2,1:2) = sor_rho_init(3:4)
+ ! Use ml orbit as the first trial in the MCMC chain
+ this%orb_ml_prm=copy(storb%orb_ml_cmp)
+ print*, 'Chi2min prm, cmp', storb%chi2_min_prm, storb%chi2_min_cmp
+ this%chi2_min_prm = storb%chi2_min_cmp
+ call nullify(storb)
+
+ ! Start actual Markov chain
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(2X,A)") "*****************"
+ WRITE(stdout,"(2X,A,1X,I0,1X,A)") "FINAL SAMPLING FOR",
this%sor_norb_prm, "SAMPLE ORBITS..."
+ WRITE(stdout,"(2X,A)") "*****************"
+ END IF
+ CALL randomWalkRanging(this)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoRandomWalkRanging", &
+ "TRACE BACK (50)", 1)
+ STOP
+ END IF
+ this%sor_niter_cmp = 3
+ print*, 'Chi2min prm, cmp', this%chi2_min_prm, this%chi2_min_cmp
+ storb=copy(this)
+ CALL setRangeBounds(storb)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoRandomWalkRanging", &
+ "TRACE BACK (20)", 1)
+ STOP
+ END IF
+ call nullify(storb)
+
+ END SUBROUTINE autoRandomWalkRanging
+
+
+
!! *Description:*
!!
- !! Markov Chain Monte Carlo ranging method.
- !! - Burn in phase: discard accepted orbits until (R.A., Dec.)
residual rms are reasonable.
- !! - Random-walk option (ran_walk) NOTE: parameter should be made
accessible from conf!
+ !! Random-walk ranging method.
+ !! - Burn in phase: discard accepted orbits until first chi2 in the
phase-space are
+ !! defined by dchi2.
+ !! - Random-walk option (ran_walk) can be turned off, which reverts
the routine to
+ !! a version of MCMC ranging
!!
!! Returns error.
!!
- SUBROUTINE MCMCRanging4(this)
+ SUBROUTINE randomWalkRanging(this)
IMPLICIT NONE
TYPE(StochasticOrbit), INTENT(inout) :: this
TYPE(Orbit) :: orb
- TYPE(SphericalCoordinates), DIMENSION(:), POINTER :: obs_scoords,
comp_scoords !observed
- !and
- !computed
- !sky
- !positions
+ TYPE(SphericalCoordinates), DIMENSION(:), POINTER :: obs_scoords,
comp_scoords
TYPE(SphericalCoordinates) :: obs_scoord1, obs_scoord2, obs_help
TYPE(CartesianCoordinates), DIMENSION(:), POINTER :: obsy_ccoords
TYPE(CartesianCoordinates) :: obs_ccoord_helio1, &
@@ -6294,7 +6537,7 @@
LOGICAL, DIMENSION(:,:), POINTER :: obs_masks
LOGICAL :: accepted, first = .TRUE., burn_in = .TRUE., frst_found
= .FALSE., ran_walk = .false., burnin_case=.false.
- ! SWITCH TO RANDOM WALK RANGING (-> move later to conf)
+ ! SWITCH TO RANDOM WALK RANGING (remnant: no need now when separated
from mcmc)
ran_walk = .true.
IF (info_verb >= 2) THEN
@@ -6304,7 +6547,7 @@
IF (.NOT. this%is_initialized_prm) THEN
error = .TRUE.
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"Object has not yet been initialized.", 1)
RETURN
END IF
@@ -6314,7 +6557,7 @@
DEALLOCATE(this%orb_arr_cmp, stat=err)
IF (err /= 0) THEN
error = .TRUE.
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"Could not deallocate memory (5).", 1)
RETURN
END IF
@@ -6322,7 +6565,7 @@
ALLOCATE(this%orb_arr_cmp(this%sor_norb_prm), stat=err)
IF (err /= 0) THEN
error = .TRUE.
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"Could not allocate memory (5).", 1)
RETURN
END IF
@@ -6340,7 +6583,7 @@
stat=err)
IF (err /= 0) THEN
error = .TRUE.
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"Could not allocate memory (10).", 1)
RETURN
END IF
@@ -6350,7 +6593,7 @@
DEALLOCATE(this%pdf_arr_cmp, stat=err)
IF (err /= 0) THEN
error = .TRUE.
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"Could not deallocate memory (10).", 1)
RETURN
END IF
@@ -6358,7 +6601,7 @@
ALLOCATE(this%pdf_arr_cmp(this%sor_norb_prm), stat=err)
IF (err /= 0) THEN
error = .TRUE.
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"Could not allocate memory (15).", 1)
RETURN
END IF
@@ -6368,7 +6611,7 @@
! Get observed sky positions (original observational data)
obs_scoords => getObservationSCoords(this%obss)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"TRACE BACK (5)", 1)
RETURN
END IF
@@ -6377,7 +6620,7 @@
CALL rotateToEquatorial(obs_scoords(i))
obs_coords(i,:) = getCoordinates(obs_scoords(i))
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"TRACE BACK (10)", 1)
CALL NULLIFY(obs_scoords(i))
RETURN
@@ -6388,7 +6631,7 @@
obs_masks => getObservationMasks(this%obss)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"TRACE BACK (15)", 1)
RETURN
END IF
@@ -6401,7 +6644,7 @@
information_matrix_obs => getBlockDiagInformationMatrix(this%obss)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"TRACE BACK (20)", 1)
DEALLOCATE(obsy_ccoords, stat=err)
DEALLOCATE(residuals, stat=err)
@@ -6413,7 +6656,7 @@
obsy_ccoords => getObservatoryCCoords(this%obss)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"TRACE BACK (25)", 1)
RETURN
END IF
@@ -6423,7 +6666,7 @@
obs_pair = RESHAPE(this%sor_pair_arr_prm, (/ 2 /))
ELSE
error = .TRUE.
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"Use of random observation pairs not accepted.", 1)
RETURN
END IF
@@ -6431,7 +6674,7 @@
obs_scoord1 = copy(obs_scoords(obs_pair(1)))
obs_scoord2 = copy(obs_scoords(obs_pair(2)))
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"TRACE BACK (30)", 1)
RETURN
END IF
@@ -6458,19 +6701,24 @@
state_(1:3) = getPosition(obs_scoord1)
state_(4:6) = getPosition(obs_scoord2)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"TRACE BACK (50)", 1)
DEALLOCATE(obs_masks, stat=err)
RETURN
END IF
! Center the first trial to the ml orbit if available
- !if (exist(this%orb_ml_cmp)) then
- !CALL getEphemerides(this%orb_ml_cmp, obsy_ccoords, comp_scoords)
- !state_(1)=comp_scoords(obs_pair(1),1)
- !state_(4)=comp_scoords(obs_pair(2),1)
- ! else
- state_(1) = (this%sor_rho_prm(1,2) + this%sor_rho_prm(1,1))/2.0_bp
- state_(4) = state_(1) + (this%sor_rho_prm(2,2) +
this%sor_rho_prm(2,1))/2.0_bp
+ IF (exist(this%orb_ml_prm)) then
+ CALL getEphemerides(this%orb_ml_prm, obsy_ccoords, comp_scoords)
+ state_(1)=getDistance(comp_scoords(obs_pair(1)))
+ state_(4)=getDistance(comp_scoords(obs_pair(2)))
+ print*,state_(1), state_(4)
+ else
+ state_(1) = (this%sor_rho_prm(1,2) + this%sor_rho_prm(1,1))/2.0_bp
+ state_(4) = state_(1) + (this%sor_rho_prm(2,2) +
this%sor_rho_prm(2,1))/2.0_bp
+ endif
+ print*,"...centered around ", state_(1), state_(4)
+ !state_(1) = (this%sor_rho_prm(1,2) + this%sor_rho_prm(1,1))/2.0_bp
+ !state_(4) = state_(1) + (this%sor_rho_prm(2,2) +
this%sor_rho_prm(2,1))/2.0_bp
! set counter some parameters
this%sor_norb_cmp = 0
@@ -6523,7 +6771,7 @@
CALL NEW(obs_scoord1, state(1), state(2), state(3),
getTime(obsy_ccoords(obs_pair(1))))
CALL NEW(obs_scoord2, state(4), state(5), state(6),
getTime(obsy_ccoords(obs_pair(2))))
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"TRACE BACK (65)", 1)
RETURN
END IF
@@ -6534,7 +6782,7 @@
CALL NEW(obs_ccoord_topo1, obs_scoord1)
CALL NEW(obs_ccoord_topo2, obs_scoord2)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"TRACE BACK (70)", 1)
RETURN
END IF
@@ -6545,7 +6793,7 @@
CALL rotateToEcliptic(obsy_ccoords(obs_pair(1)))
CALL rotateToEcliptic(obsy_ccoords(obs_pair(2)))
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"TRACE BACK (75)", 1)
RETURN
END IF
@@ -6556,7 +6804,7 @@
obs_ccoord_helio1 = copy(obsy_ccoords(obs_pair(1)) +
obs_ccoord_topo1)
obs_ccoord_helio2 = copy(obsy_ccoords(obs_pair(2)) +
obs_ccoord_topo2)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"TRACE BACK (80)", 1)
RETURN
END IF
@@ -6564,7 +6812,7 @@
CALL estimateLightTime(obs_ccoord_helio1, state(1))! changed from 3
and 6
CALL estimateLightTime(obs_ccoord_helio2, state(1)+state(4))
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"TRACE BACK (85)", 1)
RETURN
END IF
@@ -6580,7 +6828,7 @@
integrator=this%integrator_prm, &
integration_step=this%integration_step_prm)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"TRACE BACK (86)", 1)
error = .FALSE.
CYCLE
@@ -6612,7 +6860,7 @@
END IF
END IF
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"TRACE BACK (xx)", 1)
error = .FALSE.
CYCLE
@@ -6687,7 +6935,7 @@
IF (.NOT. equal(this%t_inv_prm,getTime(orb))) THEN
CALL propagate(orb, this%t_inv_prm)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"TRACE BACK (90)", 1)
error = .FALSE.
CYCLE
@@ -6706,7 +6954,7 @@
! get computed sky positions - spherical coordinates
CALL getEphemerides(orb, obsy_ccoords, comp_scoords,
partials_arr=partials_arr)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"TRACE BACK (45)", 1)
DEALLOCATE(comp_scoords, stat=err)
DEALLOCATE(partials_arr, stat=err)
@@ -6727,7 +6975,7 @@
CALL rotateToEquatorial(comp_scoords(i))
comp_coords(i,:) = getCoordinates(comp_scoords(i))
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"TRACE BACK (55)", 1)
DEALLOCATE(obsy_ccoords, stat=err)
DEALLOCATE(obs_masks, stat=err)
@@ -6759,7 +7007,7 @@
chi2 = chi_square(residuals, information_matrix_obs, errstr=errstr)
IF (LEN_TRIM(errstr) /= 0) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"TRACE BACK (60)", 1)
DEALLOCATE(obs_masks, stat=err)
RETURN
@@ -6777,7 +7025,7 @@
DEALLOCATE(partials_arr, stat=err)
jac_sph_inv = ABS(determinant(jacobian_matrix, errstr))
IF (LEN_TRIM(errstr) /= 0) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"Unsuccessful computation of determinant of orbital
element " // &
"jacobian matrix: " // TRIM(errstr), 1)
CALL matrix_print(jacobian_matrix, stderr, errstr)
@@ -6919,7 +7167,7 @@
! From ranging:
! "- generated distance to object at second epoch (NB!
Relative
! to the distance at the first epoch.)"
- this%sor_rho_arr_cmp(this%sor_norb_cmp,2) = state(1) -
state(4)
+ this%sor_rho_arr_cmp(this%sor_norb_cmp,2) = state(4) ! -
state(4)
! ELSE IF (count_dchi2 > this%count_dchi2_prm) then
! burn_in = .FALSE.
END IF
@@ -6971,6 +7219,7 @@
CALL propagate(this%orb_arr_cmp, this%t_inv_prm)
i = MAXLOC(this%pdf_arr_cmp,dim=1)
this%orb_ml_cmp = copy(this%orb_arr_cmp(i))
+ this%chi2_min_cmp = chi2min
IF (info_verb >= 2) THEN
WRITE(stdout,"(3(A,1X))") "Element types for covariance matrix and
ML orbit:", &
this%cov_type_prm, getElementType(this%orb_ml_cmp)
@@ -6982,7 +7231,7 @@
END IF
ELSE
error = .TRUE.
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"No sample orbits found!", 1)
! RETURN
END IF
@@ -6996,12 +7245,12 @@
information_matrix_obs, stat=err)
IF (err /= 0) THEN
error = .TRUE.
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / randomWalkRanging", &
"Could not deallocate memory (20).", 1)
RETURN
END IF
- END SUBROUTINE MCMCRanging4
+ END SUBROUTINE randomWalkRanging
@@ -7009,12 +7258,13 @@
!! *Description:*
!!
- !! Markov Chain Monte Carlo ranging method.
+ !! Markov Chain Monte Carlo ranging method.
+ !! - Currently stable version! (Pre-2015 version, see MCMCRanging2)
!! - Burn in phase: discard accepted orbits until (R.A., Dec.)
residual rms are reasonable.
!!
!! Returns error.
!!
- SUBROUTINE MCMCRanging3(this)
+ SUBROUTINE MCMCRanging(this)
IMPLICIT NONE
TYPE(StochasticOrbit), INTENT(inout) :: this
@@ -7057,7 +7307,7 @@
IF (.NOT. this%is_initialized_prm) THEN
error = .TRUE.
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"Object has not yet been initialized.", 1)
RETURN
END IF
@@ -7067,7 +7317,7 @@
DEALLOCATE(this%orb_arr_cmp, stat=err)
IF (err /= 0) THEN
error = .TRUE.
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"Could not deallocate memory (5).", 1)
RETURN
END IF
@@ -7075,7 +7325,7 @@
ALLOCATE(this%orb_arr_cmp(this%sor_norb_prm), stat=err)
IF (err /= 0) THEN
error = .TRUE.
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"Could not allocate memory (5).", 1)
RETURN
END IF
@@ -7093,7 +7343,7 @@
stat=err)
IF (err /= 0) THEN
error = .TRUE.
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"Could not allocate memory (10).", 1)
RETURN
END IF
@@ -7103,7 +7353,7 @@
DEALLOCATE(this%pdf_arr_cmp, stat=err)
IF (err /= 0) THEN
error = .TRUE.
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"Could not deallocate memory (10).", 1)
RETURN
END IF
@@ -7111,7 +7361,7 @@
ALLOCATE(this%pdf_arr_cmp(this%sor_norb_prm), stat=err)
IF (err /= 0) THEN
error = .TRUE.
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"Could not allocate memory (15).", 1)
RETURN
END IF
@@ -7121,7 +7371,7 @@
! Get observed sky positions (original observational data)
obs_scoords => getObservationSCoords(this%obss)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"TRACE BACK (5)", 1)
RETURN
END IF
@@ -7130,7 +7380,7 @@
CALL rotateToEquatorial(obs_scoords(i))
obs_coords(i,:) = getCoordinates(obs_scoords(i))
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"TRACE BACK (10)", 1)
CALL NULLIFY(obs_scoords(i))
RETURN
@@ -7141,7 +7391,7 @@
obs_masks => getObservationMasks(this%obss)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"TRACE BACK (15)", 1)
RETURN
END IF
@@ -7154,7 +7404,7 @@
information_matrix_obs => getBlockDiagInformationMatrix(this%obss)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"TRACE BACK (20)", 1)
DEALLOCATE(obsy_ccoords, stat=err)
DEALLOCATE(residuals, stat=err)
@@ -7166,7 +7416,7 @@
obsy_ccoords => getObservatoryCCoords(this%obss)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"TRACE BACK (25)", 1)
RETURN
END IF
@@ -7176,7 +7426,7 @@
obs_pair = RESHAPE(this%sor_pair_arr_prm, (/ 2 /))
ELSE
error = .TRUE.
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"Use of random observation pairs not accepted.", 1)
RETURN
END IF
@@ -7184,7 +7434,7 @@
obs_scoord1 = copy(obs_scoords(obs_pair(1)))
obs_scoord2 = copy(obs_scoords(obs_pair(2)))
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"TRACE BACK (30)", 1)
RETURN
END IF
@@ -7211,12 +7461,12 @@
state_(1:3) = getPosition(obs_scoord1)
state_(4:6) = getPosition(obs_scoord2)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"TRACE BACK (50)", 1)
DEALLOCATE(obs_masks, stat=err)
RETURN
END IF
- ! Center the first trial to the ml orbit if available -> this could be
taken outside mcmcranging?
+ ! Center the first trial to the ml orbit if available
IF (exist(this%orb_ml_prm)) then
CALL getEphemerides(this%orb_ml_prm, obsy_ccoords, comp_scoords)
state_(1)=getDistance(comp_scoords(obs_pair(1)))
@@ -7279,7 +7529,7 @@
CALL NEW(obs_scoord1, state(1), state(2), state(3),
getTime(obsy_ccoords(obs_pair(1))))
CALL NEW(obs_scoord2, state(4), state(5), state(6),
getTime(obsy_ccoords(obs_pair(2))))
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"TRACE BACK (65)", 1)
RETURN
END IF
@@ -7290,7 +7540,7 @@
CALL NEW(obs_ccoord_topo1, obs_scoord1)
CALL NEW(obs_ccoord_topo2, obs_scoord2)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"TRACE BACK (70)", 1)
RETURN
END IF
@@ -7301,7 +7551,7 @@
CALL rotateToEcliptic(obsy_ccoords(obs_pair(1)))
CALL rotateToEcliptic(obsy_ccoords(obs_pair(2)))
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"TRACE BACK (75)", 1)
RETURN
END IF
@@ -7312,7 +7562,7 @@
obs_ccoord_helio1 = copy(obsy_ccoords(obs_pair(1)) +
obs_ccoord_topo1)
obs_ccoord_helio2 = copy(obsy_ccoords(obs_pair(2)) +
obs_ccoord_topo2)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"TRACE BACK (80)", 1)
RETURN
END IF
@@ -7320,7 +7570,7 @@
CALL estimateLightTime(obs_ccoord_helio1, state(1))! changed from 3
and 6
CALL estimateLightTime(obs_ccoord_helio2, state(1)+state(4))
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"TRACE BACK (85)", 1)
RETURN
END IF
@@ -7336,7 +7586,7 @@
integrator=this%integrator_prm, &
integration_step=this%integration_step_prm)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"TRACE BACK (86)", 1)
error = .FALSE.
CYCLE
@@ -7437,7 +7687,7 @@
IF (.NOT. equal(this%t_inv_prm,getTime(orb))) THEN
CALL propagate(orb, this%t_inv_prm)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"TRACE BACK (90)", 1)
error = .FALSE.
CYCLE
@@ -7456,7 +7706,7 @@
! get computed sky positions - spherical coordinates
CALL getEphemerides(orb, obsy_ccoords, comp_scoords,
partials_arr=partials_arr)
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"TRACE BACK (45)", 1)
DEALLOCATE(comp_scoords, stat=err)
DEALLOCATE(partials_arr, stat=err)
@@ -7477,7 +7727,7 @@
CALL rotateToEquatorial(comp_scoords(i))
comp_coords(i,:) = getCoordinates(comp_scoords(i))
IF (error) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"TRACE BACK (55)", 1)
DEALLOCATE(obsy_ccoords, stat=err)
DEALLOCATE(obs_masks, stat=err)
@@ -7509,7 +7759,7 @@
chi2 = chi_square(residuals, information_matrix_obs, errstr=errstr)
IF (LEN_TRIM(errstr) /= 0) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"TRACE BACK (60)", 1)
DEALLOCATE(obs_masks, stat=err)
RETURN
@@ -7527,7 +7777,7 @@
DEALLOCATE(partials_arr, stat=err)
jac_sph_inv = ABS(determinant(jacobian_matrix, errstr))
IF (LEN_TRIM(errstr) /= 0) THEN
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"Unsuccessful computation of determinant of orbital
element " // &
"jacobian matrix: " // TRIM(errstr), 1)
CALL matrix_print(jacobian_matrix, stderr, errstr)
@@ -7646,7 +7896,7 @@
END IF
- IF (this%sor_norb_cmp > 1) THEN
+ IF (this%sor_norb_cmp >= 1) THEN
this%repetition_arr_cmp(this%sor_norb_cmp) =
this%repetition_arr_cmp(this%sor_norb_cmp) + 1
!print*,this%sor_norb_cmp,
this%repetition_arr_cmp(this%sor_norb_cmp)
END IF
@@ -7692,7 +7942,7 @@
END IF
ELSE
error = .TRUE.
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"No sample orbits found!", 1)
! RETURN
END IF
@@ -7706,17 +7956,19 @@
information_matrix_obs, stat=err)
IF (err /= 0) THEN
error = .TRUE.
- CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ CALL errorMessage("StochasticOrbit / MCMCRanging", &
"Could not deallocate memory (20).", 1)
RETURN
END IF
- END SUBROUTINE MCMCRanging3
+ END SUBROUTINE MCMCRanging
+
+
***The diff for this file has been truncated for email.***
Reply all
Reply to author
Forward
0 new messages