[oorb] r267 committed - Improvements to MCMCRanging. Updates to stdout from Ranging and Constr...
2 views
Skip to first unread message
oo...@googlecode.com
Feb 25, 2015, 3:00:44 AM2/25/15
to oo...@googlegroups.com
Revision: 267
Author: jxvir...@gmail.com
Date: Wed Feb 25 08:00:17 2015 UTC
Log: Improvements to MCMCRanging. Updates to stdout from Ranging and
ConstrainRangeDistributions.
https://code.google.com/p/oorb/source/detail?r=267
Modified:
/trunk/classes/StochasticOrbit_class.f90
=======================================
--- /trunk/classes/StochasticOrbit_class.f90 Mon Sep 15 18:11:32 2014 UTC
+++ /trunk/classes/StochasticOrbit_class.f90 Wed Feb 25 08:00:17 2015 UTC
@@ -168,6 +168,8 @@
INTEGER :: sor_niter_cmp = -1
INTEGER :: sor_niter_prm = -1
INTEGER :: sor_rho_histo_cmp = 1
+ ! Define number of burn-in orbits required in MCMC ranging
+ INTEGER :: sor_burnin_prm = 1
LOGICAL, DIMENSION(4) :: sor_iterate_bounds_prm
= .TRUE.
LOGICAL :: sor_random_obs_prm
= .FALSE.
! Allow generation of rho from Gaussian p.d.f.
@@ -768,6 +770,7 @@
this%sor_niter_cmp = -1
this%sor_niter_prm = -1
this%sor_rho_histo_cmp = 1
+ this%sor_burnin_prm = 1
this%sor_random_obs_prm = .FALSE.
this%sor_gaussian_pdf_prm = .FALSE.
this%dchi2_rejection_prm = .TRUE.
@@ -1059,6 +1062,7 @@
copy_SO%sor_niter_cmp = this%sor_niter_cmp
copy_SO%sor_niter_prm = this%sor_niter_prm
copy_SO%sor_rho_histo_cmp = this%sor_rho_histo_cmp
+ copy_SO%sor_burnin_prm = this%sor_burnin_prm
copy_SO%sor_random_obs_prm = this%sor_random_obs_prm
copy_SO%sor_gaussian_pdf_prm = this%sor_gaussian_pdf_prm
copy_SO%dchi2_rejection_prm = this%dchi2_rejection_prm
@@ -5960,8 +5964,176 @@
END SUBROUTINE includeObservations
+ SUBROUTINE autoMCMCRanging3(this)
+ IMPLICIT NONE
+ TYPE(StochasticOrbit), INTENT(inout) :: this
+ TYPE(StochasticOrbit) :: storb
+ REAL(bp), DIMENSION(6) :: stdev, stdev_
+ REAL(bp), DIMENSION(4) :: sor_rho_init
+ REAL(bp) :: dchi2, ddchi2
+ INTEGER :: i, stdev_flag, niter_, burn_in_
+
+ ! Not used (Move this to option file?)
+ !this%sor_burnin_prm = 10
+ ! First iteration
+ storb = copy(this)
+ storb%sor_norb_prm = 20
+ storb%sor_ntrial_prm = this%sor_ntrial_sw_prm
+ !burn_in_=this%sor_burnin_prm
+ !storb%sor_burnin_prm = 1
+
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(A)") "***************"
+ WRITE(stdout,"(A,2(1X,I0,1X,A))") "FIRST SAMPLING FOR", &
+ storb%sor_norb_prm, "SAMPLE ORBITS WITH MAX", &
+ storb%sor_ntrial_prm, "TRIALS..."
+ WRITE(stdout,"(A)") "***************"
+ END IF
+ CALL MCMCRanging3(storb)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging3", &
+ "TRACE BACK (10)", 1)
+ STOP
+ END IF
+ CALL constrainRangeDistributions(storb, storb%obss)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging3", &
+ "TRACE BACK (15)", 1)
+ STOP
+ END IF
+ CALL setRangeBounds(storb)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging3", &
+ "TRACE BACK (20)", 1)
+ STOP
+ END IF
+ sor_rho_init(1:4) = getRangeBounds(storb)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging3", &
+ "TRACE BACK (25)", 1)
+ STOP
+ END IF
+ CALL NULLIFY(storb)
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(A)") "Updated range bounds:"
+ WRITE(stdout,*) sor_rho_init
+ END IF
+
+ ! Second iteration
+ storb = copy(this)
+ storb%sor_norb_prm = this%sor_norb_sw_prm
+ storb%sor_ntrial_prm = this%sor_ntrial_sw_prm
+ storb%sor_rho_prm(1,1:2) = sor_rho_init(1:2)
+ storb%sor_rho_prm(2,1:2) = sor_rho_init(3:4)
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(A)") "***************"
+ WRITE(stdout,"(A,1X,I0,1X,A)") "SECOND SAMPLING FOR", &
+ storb%sor_norb_prm, "SAMPLE ORBITS..."
+ WRITE(stdout,"(A)") "***************"
+ END IF
+ CALL MCMCRanging3(storb)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging3", &
+ "TRACE BACK (30)", 1)
+ STOP
+ END IF
+ CALL constrainRangeDistributions(storb, storb%obss)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging3", &
+ "TRACE BACK (35)", 1)
+ STOP
+ END IF
+ CALL setRangeBounds(storb)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging3", &
+ "TRACE BACK (40)", 1)
+ STOP
+ END IF
+ sor_rho_init(1:4) = getRangeBounds(storb)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging3", &
+ "TRACE BACK (45)", 1)
+ STOP
+ END IF
+
+ ! Compute moments
+ stdev_=getStandardDeviations(storb,storb%element_type_prm)
+ CALL NULLIFY(storb)
+
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(A)") "Updated range bounds:"
+ WRITE(stdout,*) sor_rho_init
+ END IF
+
+ ! Final iterations
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(A)") "***************"
+ WRITE(stdout,"(A,1X,I0,1X,A)") "FINAL SAMPLING FOR",
this%sor_norb_prm, "SAMPLE ORBITS..."
+ WRITE(stdout,"(A)") "***************"
+ END IF
+
+ this%sor_niter_cmp = 2
+ niter_ = this%sor_niter_cmp
+ stdev_flag = 1
+ DO WHILE (stdev_flag > 0) !.OR. (this%dchi2_rejection_prm .AND.
ddchi2 > 2.0_bp))
+ !print*,this%sor_niter_cmp,this%sor_niter_prm + niter_
+ IF (this%sor_niter_cmp >= (this%sor_niter_prm + niter_)) THEN
+ EXIT
+ END IF
+ this%sor_niter_cmp = this%sor_niter_cmp + 1
+ !print*, this%sor_niter_cmp
+
+ 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)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging3", &
+ "TRACE BACK (50)", 1)
+ STOP
+ END IF
+
+ ! Check moments
+ stdev=getStandardDeviations(this,this%element_type_prm)
+ stdev_flag = 0
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(A)") "CHECKING orbital element standard
deviations"
+ WRITE(stdout,"(A)") " from previous and current round"
+ END IF
+ DO i=1,6
+ IF (ABS(stdev(i) - stdev_(i)) > 0.1*stdev_(i))
stdev_flag=stdev_flag + 1
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(2(F10.6,1X),I0)") stdev_(i), stdev(i),
stdev_flag
+ END IF
+ END DO
+ stdev_=stdev
+
+ dchi2 = MAXVAL(this%rchi2_arr_cmp) - MINVAL(this%rchi2_arr_cmp)
+ ddchi2 = ABS(dchi2 - this%dchi2_prm)
+ print*,'ddchi2',ddchi2
+
+ storb=copy(this)
+ CALL constrainRangeDistributions(storb, storb%obss)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging3", &
+ "TRACE BACK (55)", 1)
+ STOP
+ END IF
+ CALL setRangeBounds(storb)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging3", &
+ "TRACE BACK (60)", 1)
+ STOP
+ END IF
+ sor_rho_init(1:4) = getRangeBounds(storb)
+ CALL NULLIFY(storb)
+
+ END DO
+
+ END SUBROUTINE autoMCMCRanging3
+
+
SUBROUTINE autoMCMCRanging2(this)
@@ -6097,7 +6269,682 @@
+ !! *Description:*
+ !!
+ !! Markov Chain Monte Carlo ranging method.
+ !! - Burn in phase: discard accepted orbits until (R.A., Dec.)
residual rms are reasonable.
+ !!
+ !! Returns error.
+ !!
+ SUBROUTINE MCMCRanging3(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) :: obs_scoord1, obs_scoord2, obs_help
+ TYPE(CartesianCoordinates), DIMENSION(:), POINTER :: obsy_ccoords
+ TYPE(CartesianCoordinates) :: obs_ccoord_helio1, &
+ obs_ccoord_helio2, obs_ccoord_topo1, obs_ccoord_topo2
+ TYPE(Time) :: tt
+ REAL(bp), DIMENSION(:,:,:), POINTER :: information_matrix_obs, &
+ partials_arr
+ REAL(bp), DIMENSION(:,:), ALLOCATABLE :: residuals, &
+ comp_coords, &
+ obs_coords, &
+ jacobians
+ REAL(bp), DIMENSION(:), ALLOCATABLE :: cosdec0
+ REAL(bp), DIMENSION(6,6):: information_matrix_elem, &
+ jacobian_matrix
+ REAL(bp), DIMENSION(6) :: elements, rans, state, state_, &
+ proposal_density
+ REAL(bp), DIMENSION(3) :: pos, rms
+ REAL(bp) :: chi2, ran, t1, t2, obs_coord, jac_sph_inv, &
+ jac_sph_inv_, chi2_, chi2min, avalue, a, q
+ INTEGER, DIMENSION(:,:), POINTER :: obs_pair_arr
+ INTEGER, DIMENSION(6) :: n0, n0_
+ INTEGER, DIMENSION(2) :: obs_pair
+ INTEGER :: ndof, i, itrial, err, nobs, j, count_rms, norb_acc
+ INTEGER, parameter :: rms_prm = 10 ! Make user-defined parameter?
this%rms_prm
+ LOGICAL, DIMENSION(:,:), POINTER :: obs_masks
+ LOGICAL :: accepted, first = .TRUE., burn_in = .TRUE., rms_found
= .FALSE.
+
+ IF (info_verb >= 2) THEN
+ WRITE(*,*)"Starting MCMC ranging"
+ END IF
+
+ IF (.NOT. this%is_initialized_prm) THEN
+ error = .TRUE.
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "Object has not yet been initialized.", 1)
+ RETURN
+ END IF
+
+ ! Allocate memory for solution storing - orbits
+ IF (ASSOCIATED(this%orb_arr_cmp)) THEN
+ DEALLOCATE(this%orb_arr_cmp, stat=err)
+ IF (err /= 0) THEN
+ error = .TRUE.
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "Could not deallocate memory (5).", 1)
+ RETURN
+ END IF
+ END IF
+ ALLOCATE(this%orb_arr_cmp(this%sor_norb_prm), stat=err)
+ IF (err /= 0) THEN
+ error = .TRUE.
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "Could not allocate memory (5).", 1)
+ RETURN
+ END IF
+ nobs = getNrOfObservations(this%obss)
+ ALLOCATE(this%res_arr_cmp(this%sor_norb_prm,nobs,6), &
+ this%pdf_arr_cmp(this%sor_norb_prm), &
+ residuals(nobs,6), obs_coords(nobs,6), &
+ comp_coords(nobs, 6), cosdec0(nobs), &
+ this%reg_apr_arr_cmp(this%sor_norb_prm), &
+ this%jac_arr_cmp(this%sor_norb_prm,3), &
+ this%rchi2_arr_cmp(this%sor_norb_prm), &
+ this%rms_arr_cmp(this%sor_norb_prm,6), &
+ this%repetition_arr_cmp(this%sor_norb_prm), &
+ stat=err)
+ IF (err /= 0) THEN
+ error = .TRUE.
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "Could not allocate memory (10).", 1)
+ RETURN
+ END IF
+ this%repetition_arr_cmp = 0
+ ! Allocate memory for solution storing - pdf
+ IF (ASSOCIATED(this%pdf_arr_cmp)) THEN
+ DEALLOCATE(this%pdf_arr_cmp, stat=err)
+ IF (err /= 0) THEN
+ error = .TRUE.
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "Could not deallocate memory (10).", 1)
+ RETURN
+ END IF
+ END IF
+ ALLOCATE(this%pdf_arr_cmp(this%sor_norb_prm), stat=err)
+ IF (err /= 0) THEN
+ error = .TRUE.
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "Could not allocate memory (15).", 1)
+ RETURN
+ END IF
+
+ ! Get observed sky positions (original observational data)
+ obs_scoords => getObservationSCoords(this%obss)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "TRACE BACK (5)", 1)
+ RETURN
+ END IF
+ obs_coords = 0.0_bp
+ DO i=1,nobs
+ CALL rotateToEquatorial(obs_scoords(i))
+ obs_coords(i,:) = getCoordinates(obs_scoords(i))
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "TRACE BACK (10)", 1)
+ CALL NULLIFY(obs_scoords(i))
+ RETURN
+ END IF
+ ! CALL NULLIFY(obs_scoords(i))
+ cosdec0(i) = COS(obs_coords(i,3))
+ END DO
+
+ obs_masks => getObservationMasks(this%obss)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "TRACE BACK (15)", 1)
+ RETURN
+ END IF
+ ! Observation number counter (observation mask must be up-to-date!),
+ ! construct cosine array:
+ DO i=1,6
+ n0(i) = COUNT(obs_masks(:,i))
+ END DO
+ ndof = COUNT(obs_masks) - 6
+
+ information_matrix_obs => getBlockDiagInformationMatrix(this%obss)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "TRACE BACK (20)", 1)
+ DEALLOCATE(obsy_ccoords, stat=err)
+ DEALLOCATE(residuals, stat=err)
+ DEALLOCATE(information_matrix_obs, stat=err)
+ DEALLOCATE(obs_coords, stat=err)
+ DEALLOCATE(comp_coords, stat=err)
+ RETURN
+ END IF
+
+ obsy_ccoords => getObservatoryCCoords(this%obss)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "TRACE BACK (25)", 1)
+ RETURN
+ END IF
+
+ ! Use predefined observation pair if not using random selection
+ IF (.NOT. this%sor_random_obs_prm) THEN
+ obs_pair = RESHAPE(this%sor_pair_arr_prm, (/ 2 /))
+ ELSE
+ error = .TRUE.
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "Use of random observation pairs not accepted.", 1)
+ RETURN
+ END IF
+ ! Copy observation pair
+ obs_scoord1 = copy(obs_scoords(obs_pair(1)))
+ obs_scoord2 = copy(obs_scoords(obs_pair(2)))
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "TRACE BACK (30)", 1)
+ RETURN
+ END IF
+
+ ! proposal density; assuming 3sigma rho bounds
+ proposal_density(1) = (this%sor_rho_prm(1,2) -
this%sor_rho_prm(1,1))/6.0_bp
+ proposal_density(2) = 1.0_bp /
SQRT(information_matrix_obs(obs_pair(1),2,2))
+ proposal_density(3) = 1.0_bp /
SQRT(information_matrix_obs(obs_pair(1),3,3))
+ proposal_density(4) = (this%sor_rho_prm(2,2) -
this%sor_rho_prm(2,1))/6.0_bp
+ proposal_density(5) = 1.0_bp /
SQRT(information_matrix_obs(obs_pair(2),2,2))
+ proposal_density(6) = 1.0_bp /
SQRT(information_matrix_obs(obs_pair(2),3,3))
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(A)") "Proposal density for the first observation
(y_l [au], RA [as], Dec [as]): "
+ WRITE(stdout,"(F15.7)") proposal_density(1)
+ WRITE(stdout,"(F15.7)") proposal_density(2)/rad_asec
+ WRITE(stdout,"(F15.7)") proposal_density(3)/rad_asec
+ WRITE(stdout,"(A)") "Proposal density for the second observation
(y_r [au], RA [as], Dec [as]): "
+ WRITE(stdout,"(F15.7)") proposal_density(4)
+ WRITE(stdout,"(F15.7)") proposal_density(5)/rad_asec
+ WRITE(stdout,"(F15.7)") proposal_density(6)/rad_asec
+ END IF
+
+ ! rho1, RA1, Dec1, drho12, RA2, Dec2 of the proposed orbit
+ state_(1:3) = getPosition(obs_scoord1)
+ state_(4:6) = getPosition(obs_scoord2)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "TRACE BACK (50)", 1)
+ DEALLOCATE(obs_masks, stat=err)
+ RETURN
+ END IF
+ 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
+ this%sor_ntrial_cmp = -1 ! first loop is setting the parameters
+ chi2min = HUGE(chi2min)
+ first = .TRUE.
+ itrial = 0
+
+ norb_acc=0
+ count_rms=0
+ burn_in = .true.
+ DO WHILE (this%sor_norb_cmp < this%sor_norb_prm .AND. &
+ this%sor_ntrial_cmp < this%sor_ntrial_prm)
+
+ ! Get next state(topocentric distances and sky positions for
+ ! the two observations) around the previous computed state
+ IF (itrial == 0 .AND. info_verb >= 3) THEN
+ WRITE(stdout,"(A)") "Generating new state..."
+ END IF
+
+ this%sor_ntrial_cmp = this%sor_ntrial_cmp + 1
+ itrial = itrial + 1
+ IF (info_verb >= 4) THEN
+ WRITE(stdout,"(A,I0)") "Trial state #", itrial
+ END IF
+
+ IF (first) THEN
+ state = state_
+ ELSE
+ CALL randomGaussian(rans)
+ state(1) = state_(1) + rans(1) * proposal_density(1) + rans(4) *
proposal_density(4)
+ state(2:3) = state_(2:3) + rans(2:3) * proposal_density(2:3)
+ state(4) = state_(4) + rans(1) * proposal_density(1) - rans(4) *
proposal_density(4)
+ state(5:6) = state_(5:6) + rans(5:6) * proposal_density(5:6)
+ END IF
+
+ ! Check if the distances are negative
+ IF (state(1) <= 0.0_bp .OR. state(4) <= 0.0_bp) THEN
+ IF (info_verb >= 5) THEN
+ WRITE(stdout,"(2X,A,F10.7,A)") &
+ "Failed (One or both topocentric distances smaller than
the Earth radius.)"
+ END IF
+ CYCLE
+ END IF
+
+ ! Create new spherical coordinates with generated states in
equatorial frame
+ CALL NULLIFY(obs_scoord1)
+ CALL NULLIFY(obs_scoord2)
+ 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", &
+ "TRACE BACK (65)", 1)
+ RETURN
+ END IF
+
+ ! Create new topocentric cartesian coordinates of the observations
in equatorial frame
+ CALL NULLIFY(obs_ccoord_topo1)
+ CALL NULLIFY(obs_ccoord_topo2)
+ CALL NEW(obs_ccoord_topo1, obs_scoord1)
+ CALL NEW(obs_ccoord_topo2, obs_scoord2)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "TRACE BACK (70)", 1)
+ RETURN
+ END IF
+
+ ! Rotate to ecliptic
+ CALL rotateToEcliptic(obs_ccoord_topo1)
+ CALL rotateToEcliptic(obs_ccoord_topo2)
+ CALL rotateToEcliptic(obsy_ccoords(obs_pair(1)))
+ CALL rotateToEcliptic(obsy_ccoords(obs_pair(2)))
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "TRACE BACK (75)", 1)
+ RETURN
+ END IF
+
+ ! Cartesian heliocentric coordinates
+ CALL NULLIFY(obs_ccoord_helio1)
+ CALL NULLIFY(obs_ccoord_helio2)
+ 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", &
+ "TRACE BACK (80)", 1)
+ RETURN
+ END IF
+
+ 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", &
+ "TRACE BACK (85)", 1)
+ RETURN
+ END IF
+
+ ! Find orbit candidate at the epoch of the first observation by
+ ! using the chosen method to solve the 2-point boundary value
+ ! problem:
+ CALL NULLIFY(orb)
+ CALL NEW(orb, obs_ccoord_helio1, obs_ccoord_helio2, &
+ this%sor_2point_method_prm, this%apriori_a_max_prm, &
+ center=this%center_prm, &
+ perturbers=this%perturbers_prm, &
+ integrator=this%integrator_prm, &
+ integration_step=this%integration_step_prm)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "TRACE BACK (86)", 1)
+ error = .FALSE.
+ CYCLE
+ END IF
+
+ CALL setParameters(orb, &
+ dyn_model=this%dyn_model_prm, &
+ perturbers=this%perturbers_prm, &
+ integrator=this%integrator_prm, &
+ integration_step=this%integration_step_prm)
+
+ ! checking if the epochs are the same
+ IF (.NOT. equal(this%t_inv_prm,getTime(orb))) THEN
+ CALL propagate(orb, this%t_inv_prm)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "TRACE BACK (90)", 1)
+ error = .FALSE.
+ CYCLE
+ END IF
+ END IF
+
+ IF (this%informative_apriori_prm .AND. .NOT. first) THEN
+ ! Semimajor axis:
+ a = -1.0_bp
+ IF (this%apriori_a_min_prm >= 0.0_bp .OR. &
+ this%apriori_a_max_prm >= 0.0_bp) THEN
+ a = getSemimajorAxis(orb)
+ IF (this%apriori_a_min_prm >= 0.0_bp .AND. &
+ a < this%apriori_a_min_prm) THEN
+ ! Semimajor axis too small
+ IF (info_verb >= 4) THEN
+ WRITE(stdout,"(2X,A,F13.7,A)") &
+ "Failed (semimajor axis too small: ", a, " au)"
+ END IF
+ CYCLE
+ END IF
+ IF (this%apriori_a_max_prm >= 0.0_bp .AND. &
+ a > this%apriori_a_max_prm) THEN
+ ! Semimajor axis too large
+ IF (info_verb >= 4) THEN
+ WRITE(stdout,"(2X,A,F10.7,A)") &
+ "Failed (semimajor axis too large: ", a, " au)"
+ END IF
+ CYCLE
+ END IF
+ END IF
+ ! Periapsis distance:
+ IF (this%apriori_periapsis_min_prm >= 0.0_bp .OR. &
+ this%apriori_periapsis_max_prm >= 0.0_bp) THEN
+ IF (a >= 0.0_bp) THEN
+ CALL getPeriapsisDistance(orb, q, a)
+ ELSE
+ CALL getPeriapsisDistance(orb, q)
+ END IF
+ IF (error) THEN
+ error = .FALSE.
+ CYCLE
+ END IF
+ ! Periapsis distance too small:
+ IF (this%apriori_periapsis_min_prm >= 0.0_bp .AND. &
+ q < this%apriori_periapsis_min_prm) THEN
+ IF (info_verb >= 4) THEN
+ WRITE(stdout,"(2X,A,F13.7,A)") &
+ "Failed (periapsis distance too small: ", a, " au)"
+ END IF
+ CYCLE
+ END IF
+ ! Periapsis distance too large:
+ IF (this%apriori_periapsis_max_prm >= 0.0_bp .AND. &
+ q > this%apriori_periapsis_max_prm) THEN
+ IF (info_verb >= 4) THEN
+ WRITE(stdout,"(2X,A,F10.7,A)") &
+ "Failed (periapsis distance too large: ", a, " au)"
+ END IF
+ CYCLE
+ END IF
+ END IF
+ ! Apoapsis distance:
+ IF (this%apriori_apoapsis_min_prm >= 0.0_bp .OR. &
+ this%apriori_apoapsis_max_prm >= 0.0_bp) THEN
+ IF (a >= 0.0_bp) THEN
+ CALL getApoapsisDistance(orb, Q, a)
+ ELSE
+ CALL getApoapsisDistance(orb, Q)
+ END IF
+ ! Apoapsis distance too small:
+ IF (this%apriori_apoapsis_min_prm >= 0.0_bp .AND. &
+ Q < this%apriori_apoapsis_min_prm) THEN
+ IF (info_verb >= 4) THEN
+ WRITE(stdout,"(2X,A,F13.7,A)") &
+ "Failed (apoapsis distance too small: ", a, " au)"
+ END IF
+ CYCLE
+ END IF
+ ! Apoapsis distance too large:
+ IF (this%apriori_apoapsis_max_prm >= 0.0_bp .AND. &
+ Q > this%apriori_apoapsis_max_prm) THEN
+ IF (info_verb >= 4) THEN
+ WRITE(stdout,"(2X,A,F10.7,A)") &
+ "Failed (apoapsis distance too large: ", a, " au)"
+ END IF
+ CYCLE
+ END IF
+ END IF
+ END IF
+
+ IF (info_verb >= 3) THEN
+ WRITE(stdout,"(A,1X,I0,1X,A)") "New state generated using",
itrial, "trials."
+ END IF
+ itrial = 0
+
+ !!
+ !! COMPUTE CHI2 AND PDF FOR PROPOSED ORBIT
+ !!
+
+ ! get computed sky positions - spherical coordinates
+ CALL getEphemerides(orb, obsy_ccoords, comp_scoords,
partials_arr=partials_arr)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "TRACE BACK (45)", 1)
+ DEALLOCATE(comp_scoords, stat=err)
+ DEALLOCATE(partials_arr, stat=err)
+ !DEALLOCATE(obsy_ccoords, stat=err)
+ !DEALLOCATE(obs_masks, stat=err)
+ !DEALLOCATE(residuals, stat=err)
+ !DEALLOCATE(information_matrix_obs, stat=err)
+ !DEALLOCATE(obs_coords, stat=err)
+ !DEALLOCATE(comp_coords, stat=err)
+ !RETURN
+ error = .FALSE.
+ CYCLE
+ END IF
+
+ ! Residuals
+ comp_coords = 0.0_bp
+ DO i=1,nobs
+ CALL rotateToEquatorial(comp_scoords(i))
+ comp_coords(i,:) = getCoordinates(comp_scoords(i))
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "TRACE BACK (55)", 1)
+ DEALLOCATE(obsy_ccoords, stat=err)
+ DEALLOCATE(obs_masks, stat=err)
+ DEALLOCATE(residuals, stat=err)
+ DEALLOCATE(information_matrix_obs, stat=err)
+ DEALLOCATE(obs_coords, stat=err)
+ DEALLOCATE(comp_coords, stat=err)
+ RETURN
+ END IF
+ CALL NULLIFY(comp_scoords(i))
+ END DO
+ DEALLOCATE(comp_scoords)
+ residuals(1:nobs,1:6) = obs_coords(1:nobs,1:6) - &
+ comp_coords(1:nobs,1:6)
+ residuals(1:nobs,2) = residuals(1:nobs,2) * &
+ COS(obs_coords(1:nobs,3))
+ DO i=1,nobs
+ IF (ABS(residuals(i,2)) > pi) THEN
+ obs_coord = obs_coords(i,2)
+ IF (obs_coord < comp_coords(i,2)) THEN
+ obs_coord = obs_coord + two_pi
+ ELSE
+ comp_coords(i,2) = comp_coords(i,2) + two_pi
+ END IF
+ residuals(i,2) = (obs_coord - &
+ comp_coords(i,2)) * COS(obs_coords(i,3))
+ END IF
+ END DO
+
+ chi2 = chi_square(residuals, information_matrix_obs, errstr=errstr)
+ IF (LEN_TRIM(errstr) /= 0) THEN
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "TRACE BACK (60)", 1)
+ DEALLOCATE(obs_masks, stat=err)
+ RETURN
+ END IF
+
+ ! Jacobians
+ ! Determinant of Jacobian between topocentric
+ ! coordinates (inverse problem coordinates)
+ ! and orbital parameters required for output
+ ! ("Topocentric Wrt Cartesian/Keplerian"):
+ jacobian_matrix(1:3,:) = partials_arr(1:3,:,obs_pair(1)) / &
+ cosdec0(obs_pair(1))
+ jacobian_matrix(4:6,:) = partials_arr(1:3,:,obs_pair(2)) / &
+ cosdec0(obs_pair(2))
+ DEALLOCATE(partials_arr, stat=err)
+ jac_sph_inv = ABS(determinant(jacobian_matrix, errstr))
+ IF (LEN_TRIM(errstr) /= 0) THEN
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "Unsuccessful computation of determinant of orbital
element " // &
+ "jacobian matrix: " // TRIM(errstr), 1)
+ CALL matrix_print(jacobian_matrix, stderr, errstr)
+ errstr = ""
+ CYCLE ! main loop
+ END IF
+
+ IF (first) THEN
+ chi2_ = chi2
+ jac_sph_inv_ = jac_sph_inv
+ first = .FALSE.
+ CYCLE ! start MCMC loop
+ END IF
+
+ !!
+ !! MAKE DECISION WHETHER PROPOSED ORBIT IS ACCEPTED
+ !!
+ avalue = EXP(-0.5_bp*(chi2-chi2_))*jac_sph_inv_/jac_sph_inv
+ ! JV 16.2.2015: do we ever go here???
+ IF (avalue > HUGE(avalue)) THEN
+ IF (info_verb >= 4) THEN
+ WRITE(stdout,*) "huge avalue", MIN(chi2,chi2_)
+ END IF
+ avalue = chi2_/chi2
+ burn_in = .TRUE.
+ END IF
+
+ ! Decision making accept or reeject depending on pdf
+ IF (avalue >= 1.0_bp) THEN
+ accepted = .TRUE.
+ norb_acc = norb_acc + 1
+ IF (info_verb >= 3) THEN
+ WRITE(stdout,"(A,1X,I0,1X,A,1X,E12.4,2X,A,6(1X,F20.10))") &
+ "Orbit", this%sor_norb_cmp + 1, &
+ "accepted. pdf higher:", avalue, &
+ "State:", state(1), state(2:3)/rad_deg, &
+ state(4), state(5:6)/rad_deg
+ END IF
+ ELSE IF (.NOT. burn_in) THEN
+ CALL randomNumber(ran)
+ IF (avalue > ran) THEN
+ accepted = .TRUE.
+ norb_acc = norb_acc + 1
+ IF (info_verb >= 3) THEN
+ WRITE(stdout,"(A,1X,I0,1X,A,1X,E12.4,2X,A,6(1X,F20.10))") &
+ "Orbit", this%sor_norb_cmp + 1, &
+ "accepted. pdf random:", avalue, &
+ "State:", state(1), state(2:3)/rad_deg, &
+ state(4), state(5:6)/rad_deg
+ END IF
+ END IF
+ END IF
+
+ IF (accepted) THEN
+ !WRITE(*,*) "ANALYSIS", chi2, chi2_, jac_sph_inv_, jac_sph_inv,
jac_sph_inv_/jac_sph_inv, avalue, ran, accepted
+ !!
+ !! PROPOSED ORBIT ACCEPTED
+ !!
+ n0_ = n0
+ WHERE (n0_ == 0)
+ n0_ = 1
+ END WHERE
+ ! - rms:
+ rms =
SQRT(SUM(residuals(:,:)**2.0_bp,dim=1,mask=this%obs_masks_prm)/n0_)
+
+ ! CHECK residuals: if larger than acceptance windows, burn in is
not over yet.
+ ! Require RMS_PRM (typically, 10) subsequently accepted orbits.
+ if (burn_in .and. rms(2) < this%res_accept_prm(1,2) .and. rms(3)
< this%res_accept_prm(1,3) &
+ .and. count_rms <= rms_prm) then
+ select case (count_rms)
+ case(0)
+ rms_found = .true.
+ count_rms = 1
+! case (this%sor_burnin_prm)
+ case (rms_prm)
+ burn_in = .false.
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(A,2(L1,1X),I0,1X,2(F6.3))") 'BURN-IN
OVER: ', burn_in,rms_found,count_rms,rms(2:3)/rad_asec
+ WRITE(stdout,"(A,I0,1X,A)") ' after
checking ',norb_acc,'accepted orbits'
+ WRITE(stdout,"(A,2(F5.2))") ' against residual
acceptance windows ',this%res_accept_prm(1,2:3)/rad_asec
+ END IF
+ case default
+ if (rms_found) count_rms = count_rms + 1
+ end select
+ else
+ rms_found = .false.
+ count_rms = 0
+ end if
+
+ ! Only save orbits after burn_in is over:
+ IF (.NOT.burn_in) THEN
+ this%sor_norb_cmp = this%sor_norb_cmp + 1
+ this%rchi2_arr_cmp(this%sor_norb_cmp) = chi2 - REAL(ndof,bp)
+ this%pdf_arr_cmp(this%sor_norb_cmp) = &
+ EXP(-0.5_bp * (chi2 - REAL(ndof,bp))) / &
+ jac_sph_inv
+ this%res_arr_cmp(this%sor_norb_cmp,1:nobs,1:6) =
residuals(1:nobs,1:6)
+ this%orb_arr_cmp(this%sor_norb_cmp) = copy(orb)
+! ELSE IF (count_dchi2 > this%count_dchi2_prm) then
+! burn_in = .FALSE.
+ END IF
+ chi2_ = chi2
+ state_ = state
+ jac_sph_inv_ = jac_sph_inv
+
+ IF (chi2min > chi2) THEN
+ chi2min = chi2
+ END IF
+ accepted = .FALSE.
+
+ END IF
+
+ IF (this%sor_norb_cmp > 1) THEN
+ this%repetition_arr_cmp(this%sor_norb_cmp) =
this%repetition_arr_cmp(this%sor_norb_cmp) + 1
+ END IF
+
+ IF (info_verb >= 2 .AND. MOD(this%sor_ntrial_cmp,5000) == 0) THEN
+ WRITE(stdout,"(2X,A,3(I0,1X))") "Nr of accepted orbits and
trials: ", &
+ this%sor_norb_cmp, this%sor_ntrial_cmp
+ END IF
+
+ END DO
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(A)") "END OF MCMC RANGING"
+ END IF
+
+ this%orb_arr_cmp => reallocate(this%orb_arr_cmp, this%sor_norb_cmp)
+ this%repetition_arr_cmp => reallocate(this%repetition_arr_cmp,
this%sor_norb_cmp)
+ this%res_arr_cmp => reallocate(this%res_arr_cmp, this%sor_norb_cmp,
nobs, 6)
+ this%rchi2_arr_cmp => reallocate(this%rchi2_arr_cmp, this%sor_norb_cmp)
+ this%pdf_arr_cmp => reallocate(this%pdf_arr_cmp, this%sor_norb_cmp)
+
+ IF (this%sor_norb_cmp > 0) THEN
+ 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))
+ 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)
+ WRITE(stdout,"(A,1X,I0)") "Final number of orbits:",
this%sor_norb_cmp
+ WRITE(stdout,"(A,1X,I0)") "Final number of trials:",
this%sor_ntrial_cmp
+ WRITE(stdout,"(A,1X,E7.2)") "Acceptance rate: ", &
+ REAL(this%sor_norb_cmp)/REAL(this%sor_ntrial_cmp)
+ WRITE(stdout,"(A,1X,E8.3)") "chi2 min: ", chi2min
+ END IF
+ ELSE
+ error = .TRUE.
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "No sample orbits found!", 1)
+! RETURN
+ END IF
+
+ DO i=1,nobs
+ CALL NULLIFY(obs_scoords(i))
+ CALL NULLIFY(obsy_ccoords(i))
+ END DO
+ DEALLOCATE(obs_masks, obs_coords, comp_coords, cosdec0, &
+ residuals, obs_scoords, obsy_ccoords, &
+ information_matrix_obs, stat=err)
+ IF (err /= 0) THEN
+ error = .TRUE.
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "Could not deallocate memory (20).", 1)
+ RETURN
+ END IF
+
+ END SUBROUTINE MCMCRanging3
+
!! *Description:*
!!
@@ -6666,7 +7513,7 @@
jac_sph_inv
this%res_arr_cmp(this%sor_norb_cmp,1:nobs,1:6) =
residuals(1:nobs,1:6)
this%orb_arr_cmp(this%sor_norb_cmp) = copy(orb)
- ELSE
+ ELSE
burn_in = .FALSE.
END IF
chi2_ = chi2
@@ -6700,17 +7547,24 @@
this%rchi2_arr_cmp => reallocate(this%rchi2_arr_cmp, this%sor_norb_cmp)
this%pdf_arr_cmp => reallocate(this%pdf_arr_cmp, this%sor_norb_cmp)
- 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))
- 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)
- WRITE(stdout,"(A,1X,I0)") "Final number of orbits:",
this%sor_norb_cmp
- WRITE(stdout,"(A,1X,I0)") "Final number of trials:",
this%sor_ntrial_cmp
- WRITE(stdout,"(A,1X,E7.2)") "Acceptance rate: ", &
- REAL(this%sor_norb_cmp)/REAL(this%sor_ntrial_cmp)
- WRITE(stdout,"(A,1X,E8.3)") "chi2 min: ", chi2min
+ IF (this%sor_norb_cmp > 0) THEN
+ 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))
+ 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)
+ WRITE(stdout,"(A,1X,I0)") "Final number of orbits:",
this%sor_norb_cmp
+ WRITE(stdout,"(A,1X,I0)") "Final number of trials:",
this%sor_ntrial_cmp
+ WRITE(stdout,"(A,1X,E7.2)") "Acceptance rate: ", &
+ REAL(this%sor_norb_cmp)/REAL(this%sor_ntrial_cmp)
+ WRITE(stdout,"(A,1X,E8.3)") "chi2 min: ", chi2min
+ END IF
+ ELSE
+ error = .TRUE.
+ CALL errorMessage("StochasticOrbit / MCMCRanging2", &
+ "No sample orbits found!", 1)
+ ! RETURN
END IF
DO i=1,nobs
@@ -6798,7 +7652,7 @@
"TRACE BACK (50)", 1)
STOP
END IF
- CALL constrainRangeDistributions(storb, storb%obss)
+ CALL constrainRangeDistributions(storb, storb%obss)
IF (error) THEN
CALL errorMessage("StochasticOrbit / autoMCMCRanging", &
"TRACE BACK (55)", 1)
@@ -16258,7 +17112,7 @@
!! velocity at specified epoch in two-body approximation.
!!
!! Simplified:
- !! - No random observation selection
+ !! - No random observation selection (note: is available but cannot
be used because of the multi-propagation scheme)
!!
!! Returns error.
!!
@@ -16404,6 +17258,13 @@
RETURN
END IF
+ IF (this%sor_random_obs_prm) THEN
+ error = .TRUE.
+ CALL errorMessage("StochasticOrbit / StatisticalRanging", &
+ "Use of random observation pairs not accepted.", 1)
+ RETURN
+ END IF
+
object_id = getID(this%obss)
obs_arc = getObservationalTimespan(this%obss)
@@ -17812,7 +18673,7 @@
errstr = ""
IF (err_verb >= 1) THEN
CALL matrix_print(information_matrix_elem, stderr,
errstr)
- END IF
+ END IF
errstr = ""
CYCLE
END IF
@@ -18768,7 +19629,7 @@
"RA,Dec residuals corresponding to the", COUNT(mask), &
"orbits having residuals smaller than the acceptance window
[asec]:"
DO i=1,SIZE(mask)
- IF (mask(i)) THEN
+ IF (info_verb >= 3 .AND. mask(i)) THEN
DO j=1,SIZE(residuals,dim=1)
WRITE(stdout,"(2X,2(A,1X,I0,1X),A,2(1X,F10.3))") &
"Orbit #", i, "& observation #", j, ":",
residuals(j,i,2:3)/rad_asec
@@ -18828,7 +19689,7 @@
"RA,Dec residuals corresponding to the", 10-COUNT(mask_), &
"additional orbits included [asec]:"
DO i=1,SIZE(mask)
- IF (mask(i) .AND. .NOT.mask_(i)) THEN
+ IF (info_verb >= 3 .AND. mask(i) .AND. .NOT.mask_(i)) THEN
DO j=1,SIZE(residuals,dim=1)
WRITE(stdout,"(2X,2(A,1X,I0,1X),A,2(1X,F10.3))") &
"Orbit #", i, "& observation #", j, ":",
residuals(j,i,2:3)/rad_asec
@@ -19040,14 +19901,14 @@
DEALLOCATE(sor_rho_arr_cmp, stat=err)
RETURN
END IF
- IF (info_verb >= 2) THEN
+ IF (info_verb >= 3) THEN
WRITE(stdout,"(2X,A,1X,A)") "constrainRangeDistributions:", &
"Constrained (rho1,rho2-rho1) distribution [au]: "
END IF
j = 0
***The diff for this file has been truncated for email.***
Author: jxvir...@gmail.com
Date: Wed Feb 25 08:00:17 2015 UTC
Log: Improvements to MCMCRanging. Updates to stdout from Ranging and
ConstrainRangeDistributions.
https://code.google.com/p/oorb/source/detail?r=267
Modified:
/trunk/classes/StochasticOrbit_class.f90
=======================================
--- /trunk/classes/StochasticOrbit_class.f90 Mon Sep 15 18:11:32 2014 UTC
+++ /trunk/classes/StochasticOrbit_class.f90 Wed Feb 25 08:00:17 2015 UTC
@@ -168,6 +168,8 @@
INTEGER :: sor_niter_cmp = -1
INTEGER :: sor_niter_prm = -1
INTEGER :: sor_rho_histo_cmp = 1
+ ! Define number of burn-in orbits required in MCMC ranging
+ INTEGER :: sor_burnin_prm = 1
LOGICAL, DIMENSION(4) :: sor_iterate_bounds_prm
= .TRUE.
LOGICAL :: sor_random_obs_prm
= .FALSE.
! Allow generation of rho from Gaussian p.d.f.
@@ -768,6 +770,7 @@
this%sor_niter_cmp = -1
this%sor_niter_prm = -1
this%sor_rho_histo_cmp = 1
+ this%sor_burnin_prm = 1
this%sor_random_obs_prm = .FALSE.
this%sor_gaussian_pdf_prm = .FALSE.
this%dchi2_rejection_prm = .TRUE.
@@ -1059,6 +1062,7 @@
copy_SO%sor_niter_cmp = this%sor_niter_cmp
copy_SO%sor_niter_prm = this%sor_niter_prm
copy_SO%sor_rho_histo_cmp = this%sor_rho_histo_cmp
+ copy_SO%sor_burnin_prm = this%sor_burnin_prm
copy_SO%sor_random_obs_prm = this%sor_random_obs_prm
copy_SO%sor_gaussian_pdf_prm = this%sor_gaussian_pdf_prm
copy_SO%dchi2_rejection_prm = this%dchi2_rejection_prm
@@ -5960,8 +5964,176 @@
END SUBROUTINE includeObservations
+ SUBROUTINE autoMCMCRanging3(this)
+ IMPLICIT NONE
+ TYPE(StochasticOrbit), INTENT(inout) :: this
+ TYPE(StochasticOrbit) :: storb
+ REAL(bp), DIMENSION(6) :: stdev, stdev_
+ REAL(bp), DIMENSION(4) :: sor_rho_init
+ REAL(bp) :: dchi2, ddchi2
+ INTEGER :: i, stdev_flag, niter_, burn_in_
+
+ ! Not used (Move this to option file?)
+ !this%sor_burnin_prm = 10
+ ! First iteration
+ storb = copy(this)
+ storb%sor_norb_prm = 20
+ storb%sor_ntrial_prm = this%sor_ntrial_sw_prm
+ !burn_in_=this%sor_burnin_prm
+ !storb%sor_burnin_prm = 1
+
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(A)") "***************"
+ WRITE(stdout,"(A,2(1X,I0,1X,A))") "FIRST SAMPLING FOR", &
+ storb%sor_norb_prm, "SAMPLE ORBITS WITH MAX", &
+ storb%sor_ntrial_prm, "TRIALS..."
+ WRITE(stdout,"(A)") "***************"
+ END IF
+ CALL MCMCRanging3(storb)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging3", &
+ "TRACE BACK (10)", 1)
+ STOP
+ END IF
+ CALL constrainRangeDistributions(storb, storb%obss)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging3", &
+ "TRACE BACK (15)", 1)
+ STOP
+ END IF
+ CALL setRangeBounds(storb)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging3", &
+ "TRACE BACK (20)", 1)
+ STOP
+ END IF
+ sor_rho_init(1:4) = getRangeBounds(storb)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging3", &
+ "TRACE BACK (25)", 1)
+ STOP
+ END IF
+ CALL NULLIFY(storb)
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(A)") "Updated range bounds:"
+ WRITE(stdout,*) sor_rho_init
+ END IF
+
+ ! Second iteration
+ storb = copy(this)
+ storb%sor_norb_prm = this%sor_norb_sw_prm
+ storb%sor_ntrial_prm = this%sor_ntrial_sw_prm
+ storb%sor_rho_prm(1,1:2) = sor_rho_init(1:2)
+ storb%sor_rho_prm(2,1:2) = sor_rho_init(3:4)
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(A)") "***************"
+ WRITE(stdout,"(A,1X,I0,1X,A)") "SECOND SAMPLING FOR", &
+ storb%sor_norb_prm, "SAMPLE ORBITS..."
+ WRITE(stdout,"(A)") "***************"
+ END IF
+ CALL MCMCRanging3(storb)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging3", &
+ "TRACE BACK (30)", 1)
+ STOP
+ END IF
+ CALL constrainRangeDistributions(storb, storb%obss)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging3", &
+ "TRACE BACK (35)", 1)
+ STOP
+ END IF
+ CALL setRangeBounds(storb)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging3", &
+ "TRACE BACK (40)", 1)
+ STOP
+ END IF
+ sor_rho_init(1:4) = getRangeBounds(storb)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging3", &
+ "TRACE BACK (45)", 1)
+ STOP
+ END IF
+
+ ! Compute moments
+ stdev_=getStandardDeviations(storb,storb%element_type_prm)
+ CALL NULLIFY(storb)
+
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(A)") "Updated range bounds:"
+ WRITE(stdout,*) sor_rho_init
+ END IF
+
+ ! Final iterations
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(A)") "***************"
+ WRITE(stdout,"(A,1X,I0,1X,A)") "FINAL SAMPLING FOR",
this%sor_norb_prm, "SAMPLE ORBITS..."
+ WRITE(stdout,"(A)") "***************"
+ END IF
+
+ this%sor_niter_cmp = 2
+ niter_ = this%sor_niter_cmp
+ stdev_flag = 1
+ DO WHILE (stdev_flag > 0) !.OR. (this%dchi2_rejection_prm .AND.
ddchi2 > 2.0_bp))
+ !print*,this%sor_niter_cmp,this%sor_niter_prm + niter_
+ IF (this%sor_niter_cmp >= (this%sor_niter_prm + niter_)) THEN
+ EXIT
+ END IF
+ this%sor_niter_cmp = this%sor_niter_cmp + 1
+ !print*, this%sor_niter_cmp
+
+ 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)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging3", &
+ "TRACE BACK (50)", 1)
+ STOP
+ END IF
+
+ ! Check moments
+ stdev=getStandardDeviations(this,this%element_type_prm)
+ stdev_flag = 0
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(A)") "CHECKING orbital element standard
deviations"
+ WRITE(stdout,"(A)") " from previous and current round"
+ END IF
+ DO i=1,6
+ IF (ABS(stdev(i) - stdev_(i)) > 0.1*stdev_(i))
stdev_flag=stdev_flag + 1
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(2(F10.6,1X),I0)") stdev_(i), stdev(i),
stdev_flag
+ END IF
+ END DO
+ stdev_=stdev
+
+ dchi2 = MAXVAL(this%rchi2_arr_cmp) - MINVAL(this%rchi2_arr_cmp)
+ ddchi2 = ABS(dchi2 - this%dchi2_prm)
+ print*,'ddchi2',ddchi2
+
+ storb=copy(this)
+ CALL constrainRangeDistributions(storb, storb%obss)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging3", &
+ "TRACE BACK (55)", 1)
+ STOP
+ END IF
+ CALL setRangeBounds(storb)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / autoMCMCRanging3", &
+ "TRACE BACK (60)", 1)
+ STOP
+ END IF
+ sor_rho_init(1:4) = getRangeBounds(storb)
+ CALL NULLIFY(storb)
+
+ END DO
+
+ END SUBROUTINE autoMCMCRanging3
+
+
SUBROUTINE autoMCMCRanging2(this)
@@ -6097,7 +6269,682 @@
+ !! *Description:*
+ !!
+ !! Markov Chain Monte Carlo ranging method.
+ !! - Burn in phase: discard accepted orbits until (R.A., Dec.)
residual rms are reasonable.
+ !!
+ !! Returns error.
+ !!
+ SUBROUTINE MCMCRanging3(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) :: obs_scoord1, obs_scoord2, obs_help
+ TYPE(CartesianCoordinates), DIMENSION(:), POINTER :: obsy_ccoords
+ TYPE(CartesianCoordinates) :: obs_ccoord_helio1, &
+ obs_ccoord_helio2, obs_ccoord_topo1, obs_ccoord_topo2
+ TYPE(Time) :: tt
+ REAL(bp), DIMENSION(:,:,:), POINTER :: information_matrix_obs, &
+ partials_arr
+ REAL(bp), DIMENSION(:,:), ALLOCATABLE :: residuals, &
+ comp_coords, &
+ obs_coords, &
+ jacobians
+ REAL(bp), DIMENSION(:), ALLOCATABLE :: cosdec0
+ REAL(bp), DIMENSION(6,6):: information_matrix_elem, &
+ jacobian_matrix
+ REAL(bp), DIMENSION(6) :: elements, rans, state, state_, &
+ proposal_density
+ REAL(bp), DIMENSION(3) :: pos, rms
+ REAL(bp) :: chi2, ran, t1, t2, obs_coord, jac_sph_inv, &
+ jac_sph_inv_, chi2_, chi2min, avalue, a, q
+ INTEGER, DIMENSION(:,:), POINTER :: obs_pair_arr
+ INTEGER, DIMENSION(6) :: n0, n0_
+ INTEGER, DIMENSION(2) :: obs_pair
+ INTEGER :: ndof, i, itrial, err, nobs, j, count_rms, norb_acc
+ INTEGER, parameter :: rms_prm = 10 ! Make user-defined parameter?
this%rms_prm
+ LOGICAL, DIMENSION(:,:), POINTER :: obs_masks
+ LOGICAL :: accepted, first = .TRUE., burn_in = .TRUE., rms_found
= .FALSE.
+
+ IF (info_verb >= 2) THEN
+ WRITE(*,*)"Starting MCMC ranging"
+ END IF
+
+ IF (.NOT. this%is_initialized_prm) THEN
+ error = .TRUE.
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "Object has not yet been initialized.", 1)
+ RETURN
+ END IF
+
+ ! Allocate memory for solution storing - orbits
+ IF (ASSOCIATED(this%orb_arr_cmp)) THEN
+ DEALLOCATE(this%orb_arr_cmp, stat=err)
+ IF (err /= 0) THEN
+ error = .TRUE.
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "Could not deallocate memory (5).", 1)
+ RETURN
+ END IF
+ END IF
+ ALLOCATE(this%orb_arr_cmp(this%sor_norb_prm), stat=err)
+ IF (err /= 0) THEN
+ error = .TRUE.
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "Could not allocate memory (5).", 1)
+ RETURN
+ END IF
+ nobs = getNrOfObservations(this%obss)
+ ALLOCATE(this%res_arr_cmp(this%sor_norb_prm,nobs,6), &
+ this%pdf_arr_cmp(this%sor_norb_prm), &
+ residuals(nobs,6), obs_coords(nobs,6), &
+ comp_coords(nobs, 6), cosdec0(nobs), &
+ this%reg_apr_arr_cmp(this%sor_norb_prm), &
+ this%jac_arr_cmp(this%sor_norb_prm,3), &
+ this%rchi2_arr_cmp(this%sor_norb_prm), &
+ this%rms_arr_cmp(this%sor_norb_prm,6), &
+ this%repetition_arr_cmp(this%sor_norb_prm), &
+ stat=err)
+ IF (err /= 0) THEN
+ error = .TRUE.
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "Could not allocate memory (10).", 1)
+ RETURN
+ END IF
+ this%repetition_arr_cmp = 0
+ ! Allocate memory for solution storing - pdf
+ IF (ASSOCIATED(this%pdf_arr_cmp)) THEN
+ DEALLOCATE(this%pdf_arr_cmp, stat=err)
+ IF (err /= 0) THEN
+ error = .TRUE.
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "Could not deallocate memory (10).", 1)
+ RETURN
+ END IF
+ END IF
+ ALLOCATE(this%pdf_arr_cmp(this%sor_norb_prm), stat=err)
+ IF (err /= 0) THEN
+ error = .TRUE.
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "Could not allocate memory (15).", 1)
+ RETURN
+ END IF
+
+ ! Get observed sky positions (original observational data)
+ obs_scoords => getObservationSCoords(this%obss)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "TRACE BACK (5)", 1)
+ RETURN
+ END IF
+ obs_coords = 0.0_bp
+ DO i=1,nobs
+ CALL rotateToEquatorial(obs_scoords(i))
+ obs_coords(i,:) = getCoordinates(obs_scoords(i))
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "TRACE BACK (10)", 1)
+ CALL NULLIFY(obs_scoords(i))
+ RETURN
+ END IF
+ ! CALL NULLIFY(obs_scoords(i))
+ cosdec0(i) = COS(obs_coords(i,3))
+ END DO
+
+ obs_masks => getObservationMasks(this%obss)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "TRACE BACK (15)", 1)
+ RETURN
+ END IF
+ ! Observation number counter (observation mask must be up-to-date!),
+ ! construct cosine array:
+ DO i=1,6
+ n0(i) = COUNT(obs_masks(:,i))
+ END DO
+ ndof = COUNT(obs_masks) - 6
+
+ information_matrix_obs => getBlockDiagInformationMatrix(this%obss)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "TRACE BACK (20)", 1)
+ DEALLOCATE(obsy_ccoords, stat=err)
+ DEALLOCATE(residuals, stat=err)
+ DEALLOCATE(information_matrix_obs, stat=err)
+ DEALLOCATE(obs_coords, stat=err)
+ DEALLOCATE(comp_coords, stat=err)
+ RETURN
+ END IF
+
+ obsy_ccoords => getObservatoryCCoords(this%obss)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "TRACE BACK (25)", 1)
+ RETURN
+ END IF
+
+ ! Use predefined observation pair if not using random selection
+ IF (.NOT. this%sor_random_obs_prm) THEN
+ obs_pair = RESHAPE(this%sor_pair_arr_prm, (/ 2 /))
+ ELSE
+ error = .TRUE.
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "Use of random observation pairs not accepted.", 1)
+ RETURN
+ END IF
+ ! Copy observation pair
+ obs_scoord1 = copy(obs_scoords(obs_pair(1)))
+ obs_scoord2 = copy(obs_scoords(obs_pair(2)))
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "TRACE BACK (30)", 1)
+ RETURN
+ END IF
+
+ ! proposal density; assuming 3sigma rho bounds
+ proposal_density(1) = (this%sor_rho_prm(1,2) -
this%sor_rho_prm(1,1))/6.0_bp
+ proposal_density(2) = 1.0_bp /
SQRT(information_matrix_obs(obs_pair(1),2,2))
+ proposal_density(3) = 1.0_bp /
SQRT(information_matrix_obs(obs_pair(1),3,3))
+ proposal_density(4) = (this%sor_rho_prm(2,2) -
this%sor_rho_prm(2,1))/6.0_bp
+ proposal_density(5) = 1.0_bp /
SQRT(information_matrix_obs(obs_pair(2),2,2))
+ proposal_density(6) = 1.0_bp /
SQRT(information_matrix_obs(obs_pair(2),3,3))
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(A)") "Proposal density for the first observation
(y_l [au], RA [as], Dec [as]): "
+ WRITE(stdout,"(F15.7)") proposal_density(1)
+ WRITE(stdout,"(F15.7)") proposal_density(2)/rad_asec
+ WRITE(stdout,"(F15.7)") proposal_density(3)/rad_asec
+ WRITE(stdout,"(A)") "Proposal density for the second observation
(y_r [au], RA [as], Dec [as]): "
+ WRITE(stdout,"(F15.7)") proposal_density(4)
+ WRITE(stdout,"(F15.7)") proposal_density(5)/rad_asec
+ WRITE(stdout,"(F15.7)") proposal_density(6)/rad_asec
+ END IF
+
+ ! rho1, RA1, Dec1, drho12, RA2, Dec2 of the proposed orbit
+ state_(1:3) = getPosition(obs_scoord1)
+ state_(4:6) = getPosition(obs_scoord2)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "TRACE BACK (50)", 1)
+ DEALLOCATE(obs_masks, stat=err)
+ RETURN
+ END IF
+ 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
+ this%sor_ntrial_cmp = -1 ! first loop is setting the parameters
+ chi2min = HUGE(chi2min)
+ first = .TRUE.
+ itrial = 0
+
+ norb_acc=0
+ count_rms=0
+ burn_in = .true.
+ DO WHILE (this%sor_norb_cmp < this%sor_norb_prm .AND. &
+ this%sor_ntrial_cmp < this%sor_ntrial_prm)
+
+ ! Get next state(topocentric distances and sky positions for
+ ! the two observations) around the previous computed state
+ IF (itrial == 0 .AND. info_verb >= 3) THEN
+ WRITE(stdout,"(A)") "Generating new state..."
+ END IF
+
+ this%sor_ntrial_cmp = this%sor_ntrial_cmp + 1
+ itrial = itrial + 1
+ IF (info_verb >= 4) THEN
+ WRITE(stdout,"(A,I0)") "Trial state #", itrial
+ END IF
+
+ IF (first) THEN
+ state = state_
+ ELSE
+ CALL randomGaussian(rans)
+ state(1) = state_(1) + rans(1) * proposal_density(1) + rans(4) *
proposal_density(4)
+ state(2:3) = state_(2:3) + rans(2:3) * proposal_density(2:3)
+ state(4) = state_(4) + rans(1) * proposal_density(1) - rans(4) *
proposal_density(4)
+ state(5:6) = state_(5:6) + rans(5:6) * proposal_density(5:6)
+ END IF
+
+ ! Check if the distances are negative
+ IF (state(1) <= 0.0_bp .OR. state(4) <= 0.0_bp) THEN
+ IF (info_verb >= 5) THEN
+ WRITE(stdout,"(2X,A,F10.7,A)") &
+ "Failed (One or both topocentric distances smaller than
the Earth radius.)"
+ END IF
+ CYCLE
+ END IF
+
+ ! Create new spherical coordinates with generated states in
equatorial frame
+ CALL NULLIFY(obs_scoord1)
+ CALL NULLIFY(obs_scoord2)
+ 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", &
+ "TRACE BACK (65)", 1)
+ RETURN
+ END IF
+
+ ! Create new topocentric cartesian coordinates of the observations
in equatorial frame
+ CALL NULLIFY(obs_ccoord_topo1)
+ CALL NULLIFY(obs_ccoord_topo2)
+ CALL NEW(obs_ccoord_topo1, obs_scoord1)
+ CALL NEW(obs_ccoord_topo2, obs_scoord2)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "TRACE BACK (70)", 1)
+ RETURN
+ END IF
+
+ ! Rotate to ecliptic
+ CALL rotateToEcliptic(obs_ccoord_topo1)
+ CALL rotateToEcliptic(obs_ccoord_topo2)
+ CALL rotateToEcliptic(obsy_ccoords(obs_pair(1)))
+ CALL rotateToEcliptic(obsy_ccoords(obs_pair(2)))
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "TRACE BACK (75)", 1)
+ RETURN
+ END IF
+
+ ! Cartesian heliocentric coordinates
+ CALL NULLIFY(obs_ccoord_helio1)
+ CALL NULLIFY(obs_ccoord_helio2)
+ 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", &
+ "TRACE BACK (80)", 1)
+ RETURN
+ END IF
+
+ 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", &
+ "TRACE BACK (85)", 1)
+ RETURN
+ END IF
+
+ ! Find orbit candidate at the epoch of the first observation by
+ ! using the chosen method to solve the 2-point boundary value
+ ! problem:
+ CALL NULLIFY(orb)
+ CALL NEW(orb, obs_ccoord_helio1, obs_ccoord_helio2, &
+ this%sor_2point_method_prm, this%apriori_a_max_prm, &
+ center=this%center_prm, &
+ perturbers=this%perturbers_prm, &
+ integrator=this%integrator_prm, &
+ integration_step=this%integration_step_prm)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "TRACE BACK (86)", 1)
+ error = .FALSE.
+ CYCLE
+ END IF
+
+ CALL setParameters(orb, &
+ dyn_model=this%dyn_model_prm, &
+ perturbers=this%perturbers_prm, &
+ integrator=this%integrator_prm, &
+ integration_step=this%integration_step_prm)
+
+ ! checking if the epochs are the same
+ IF (.NOT. equal(this%t_inv_prm,getTime(orb))) THEN
+ CALL propagate(orb, this%t_inv_prm)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "TRACE BACK (90)", 1)
+ error = .FALSE.
+ CYCLE
+ END IF
+ END IF
+
+ IF (this%informative_apriori_prm .AND. .NOT. first) THEN
+ ! Semimajor axis:
+ a = -1.0_bp
+ IF (this%apriori_a_min_prm >= 0.0_bp .OR. &
+ this%apriori_a_max_prm >= 0.0_bp) THEN
+ a = getSemimajorAxis(orb)
+ IF (this%apriori_a_min_prm >= 0.0_bp .AND. &
+ a < this%apriori_a_min_prm) THEN
+ ! Semimajor axis too small
+ IF (info_verb >= 4) THEN
+ WRITE(stdout,"(2X,A,F13.7,A)") &
+ "Failed (semimajor axis too small: ", a, " au)"
+ END IF
+ CYCLE
+ END IF
+ IF (this%apriori_a_max_prm >= 0.0_bp .AND. &
+ a > this%apriori_a_max_prm) THEN
+ ! Semimajor axis too large
+ IF (info_verb >= 4) THEN
+ WRITE(stdout,"(2X,A,F10.7,A)") &
+ "Failed (semimajor axis too large: ", a, " au)"
+ END IF
+ CYCLE
+ END IF
+ END IF
+ ! Periapsis distance:
+ IF (this%apriori_periapsis_min_prm >= 0.0_bp .OR. &
+ this%apriori_periapsis_max_prm >= 0.0_bp) THEN
+ IF (a >= 0.0_bp) THEN
+ CALL getPeriapsisDistance(orb, q, a)
+ ELSE
+ CALL getPeriapsisDistance(orb, q)
+ END IF
+ IF (error) THEN
+ error = .FALSE.
+ CYCLE
+ END IF
+ ! Periapsis distance too small:
+ IF (this%apriori_periapsis_min_prm >= 0.0_bp .AND. &
+ q < this%apriori_periapsis_min_prm) THEN
+ IF (info_verb >= 4) THEN
+ WRITE(stdout,"(2X,A,F13.7,A)") &
+ "Failed (periapsis distance too small: ", a, " au)"
+ END IF
+ CYCLE
+ END IF
+ ! Periapsis distance too large:
+ IF (this%apriori_periapsis_max_prm >= 0.0_bp .AND. &
+ q > this%apriori_periapsis_max_prm) THEN
+ IF (info_verb >= 4) THEN
+ WRITE(stdout,"(2X,A,F10.7,A)") &
+ "Failed (periapsis distance too large: ", a, " au)"
+ END IF
+ CYCLE
+ END IF
+ END IF
+ ! Apoapsis distance:
+ IF (this%apriori_apoapsis_min_prm >= 0.0_bp .OR. &
+ this%apriori_apoapsis_max_prm >= 0.0_bp) THEN
+ IF (a >= 0.0_bp) THEN
+ CALL getApoapsisDistance(orb, Q, a)
+ ELSE
+ CALL getApoapsisDistance(orb, Q)
+ END IF
+ ! Apoapsis distance too small:
+ IF (this%apriori_apoapsis_min_prm >= 0.0_bp .AND. &
+ Q < this%apriori_apoapsis_min_prm) THEN
+ IF (info_verb >= 4) THEN
+ WRITE(stdout,"(2X,A,F13.7,A)") &
+ "Failed (apoapsis distance too small: ", a, " au)"
+ END IF
+ CYCLE
+ END IF
+ ! Apoapsis distance too large:
+ IF (this%apriori_apoapsis_max_prm >= 0.0_bp .AND. &
+ Q > this%apriori_apoapsis_max_prm) THEN
+ IF (info_verb >= 4) THEN
+ WRITE(stdout,"(2X,A,F10.7,A)") &
+ "Failed (apoapsis distance too large: ", a, " au)"
+ END IF
+ CYCLE
+ END IF
+ END IF
+ END IF
+
+ IF (info_verb >= 3) THEN
+ WRITE(stdout,"(A,1X,I0,1X,A)") "New state generated using",
itrial, "trials."
+ END IF
+ itrial = 0
+
+ !!
+ !! COMPUTE CHI2 AND PDF FOR PROPOSED ORBIT
+ !!
+
+ ! get computed sky positions - spherical coordinates
+ CALL getEphemerides(orb, obsy_ccoords, comp_scoords,
partials_arr=partials_arr)
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "TRACE BACK (45)", 1)
+ DEALLOCATE(comp_scoords, stat=err)
+ DEALLOCATE(partials_arr, stat=err)
+ !DEALLOCATE(obsy_ccoords, stat=err)
+ !DEALLOCATE(obs_masks, stat=err)
+ !DEALLOCATE(residuals, stat=err)
+ !DEALLOCATE(information_matrix_obs, stat=err)
+ !DEALLOCATE(obs_coords, stat=err)
+ !DEALLOCATE(comp_coords, stat=err)
+ !RETURN
+ error = .FALSE.
+ CYCLE
+ END IF
+
+ ! Residuals
+ comp_coords = 0.0_bp
+ DO i=1,nobs
+ CALL rotateToEquatorial(comp_scoords(i))
+ comp_coords(i,:) = getCoordinates(comp_scoords(i))
+ IF (error) THEN
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "TRACE BACK (55)", 1)
+ DEALLOCATE(obsy_ccoords, stat=err)
+ DEALLOCATE(obs_masks, stat=err)
+ DEALLOCATE(residuals, stat=err)
+ DEALLOCATE(information_matrix_obs, stat=err)
+ DEALLOCATE(obs_coords, stat=err)
+ DEALLOCATE(comp_coords, stat=err)
+ RETURN
+ END IF
+ CALL NULLIFY(comp_scoords(i))
+ END DO
+ DEALLOCATE(comp_scoords)
+ residuals(1:nobs,1:6) = obs_coords(1:nobs,1:6) - &
+ comp_coords(1:nobs,1:6)
+ residuals(1:nobs,2) = residuals(1:nobs,2) * &
+ COS(obs_coords(1:nobs,3))
+ DO i=1,nobs
+ IF (ABS(residuals(i,2)) > pi) THEN
+ obs_coord = obs_coords(i,2)
+ IF (obs_coord < comp_coords(i,2)) THEN
+ obs_coord = obs_coord + two_pi
+ ELSE
+ comp_coords(i,2) = comp_coords(i,2) + two_pi
+ END IF
+ residuals(i,2) = (obs_coord - &
+ comp_coords(i,2)) * COS(obs_coords(i,3))
+ END IF
+ END DO
+
+ chi2 = chi_square(residuals, information_matrix_obs, errstr=errstr)
+ IF (LEN_TRIM(errstr) /= 0) THEN
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "TRACE BACK (60)", 1)
+ DEALLOCATE(obs_masks, stat=err)
+ RETURN
+ END IF
+
+ ! Jacobians
+ ! Determinant of Jacobian between topocentric
+ ! coordinates (inverse problem coordinates)
+ ! and orbital parameters required for output
+ ! ("Topocentric Wrt Cartesian/Keplerian"):
+ jacobian_matrix(1:3,:) = partials_arr(1:3,:,obs_pair(1)) / &
+ cosdec0(obs_pair(1))
+ jacobian_matrix(4:6,:) = partials_arr(1:3,:,obs_pair(2)) / &
+ cosdec0(obs_pair(2))
+ DEALLOCATE(partials_arr, stat=err)
+ jac_sph_inv = ABS(determinant(jacobian_matrix, errstr))
+ IF (LEN_TRIM(errstr) /= 0) THEN
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "Unsuccessful computation of determinant of orbital
element " // &
+ "jacobian matrix: " // TRIM(errstr), 1)
+ CALL matrix_print(jacobian_matrix, stderr, errstr)
+ errstr = ""
+ CYCLE ! main loop
+ END IF
+
+ IF (first) THEN
+ chi2_ = chi2
+ jac_sph_inv_ = jac_sph_inv
+ first = .FALSE.
+ CYCLE ! start MCMC loop
+ END IF
+
+ !!
+ !! MAKE DECISION WHETHER PROPOSED ORBIT IS ACCEPTED
+ !!
+ avalue = EXP(-0.5_bp*(chi2-chi2_))*jac_sph_inv_/jac_sph_inv
+ ! JV 16.2.2015: do we ever go here???
+ IF (avalue > HUGE(avalue)) THEN
+ IF (info_verb >= 4) THEN
+ WRITE(stdout,*) "huge avalue", MIN(chi2,chi2_)
+ END IF
+ avalue = chi2_/chi2
+ burn_in = .TRUE.
+ END IF
+
+ ! Decision making accept or reeject depending on pdf
+ IF (avalue >= 1.0_bp) THEN
+ accepted = .TRUE.
+ norb_acc = norb_acc + 1
+ IF (info_verb >= 3) THEN
+ WRITE(stdout,"(A,1X,I0,1X,A,1X,E12.4,2X,A,6(1X,F20.10))") &
+ "Orbit", this%sor_norb_cmp + 1, &
+ "accepted. pdf higher:", avalue, &
+ "State:", state(1), state(2:3)/rad_deg, &
+ state(4), state(5:6)/rad_deg
+ END IF
+ ELSE IF (.NOT. burn_in) THEN
+ CALL randomNumber(ran)
+ IF (avalue > ran) THEN
+ accepted = .TRUE.
+ norb_acc = norb_acc + 1
+ IF (info_verb >= 3) THEN
+ WRITE(stdout,"(A,1X,I0,1X,A,1X,E12.4,2X,A,6(1X,F20.10))") &
+ "Orbit", this%sor_norb_cmp + 1, &
+ "accepted. pdf random:", avalue, &
+ "State:", state(1), state(2:3)/rad_deg, &
+ state(4), state(5:6)/rad_deg
+ END IF
+ END IF
+ END IF
+
+ IF (accepted) THEN
+ !WRITE(*,*) "ANALYSIS", chi2, chi2_, jac_sph_inv_, jac_sph_inv,
jac_sph_inv_/jac_sph_inv, avalue, ran, accepted
+ !!
+ !! PROPOSED ORBIT ACCEPTED
+ !!
+ n0_ = n0
+ WHERE (n0_ == 0)
+ n0_ = 1
+ END WHERE
+ ! - rms:
+ rms =
SQRT(SUM(residuals(:,:)**2.0_bp,dim=1,mask=this%obs_masks_prm)/n0_)
+
+ ! CHECK residuals: if larger than acceptance windows, burn in is
not over yet.
+ ! Require RMS_PRM (typically, 10) subsequently accepted orbits.
+ if (burn_in .and. rms(2) < this%res_accept_prm(1,2) .and. rms(3)
< this%res_accept_prm(1,3) &
+ .and. count_rms <= rms_prm) then
+ select case (count_rms)
+ case(0)
+ rms_found = .true.
+ count_rms = 1
+! case (this%sor_burnin_prm)
+ case (rms_prm)
+ burn_in = .false.
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(A,2(L1,1X),I0,1X,2(F6.3))") 'BURN-IN
OVER: ', burn_in,rms_found,count_rms,rms(2:3)/rad_asec
+ WRITE(stdout,"(A,I0,1X,A)") ' after
checking ',norb_acc,'accepted orbits'
+ WRITE(stdout,"(A,2(F5.2))") ' against residual
acceptance windows ',this%res_accept_prm(1,2:3)/rad_asec
+ END IF
+ case default
+ if (rms_found) count_rms = count_rms + 1
+ end select
+ else
+ rms_found = .false.
+ count_rms = 0
+ end if
+
+ ! Only save orbits after burn_in is over:
+ IF (.NOT.burn_in) THEN
+ this%sor_norb_cmp = this%sor_norb_cmp + 1
+ this%rchi2_arr_cmp(this%sor_norb_cmp) = chi2 - REAL(ndof,bp)
+ this%pdf_arr_cmp(this%sor_norb_cmp) = &
+ EXP(-0.5_bp * (chi2 - REAL(ndof,bp))) / &
+ jac_sph_inv
+ this%res_arr_cmp(this%sor_norb_cmp,1:nobs,1:6) =
residuals(1:nobs,1:6)
+ this%orb_arr_cmp(this%sor_norb_cmp) = copy(orb)
+! ELSE IF (count_dchi2 > this%count_dchi2_prm) then
+! burn_in = .FALSE.
+ END IF
+ chi2_ = chi2
+ state_ = state
+ jac_sph_inv_ = jac_sph_inv
+
+ IF (chi2min > chi2) THEN
+ chi2min = chi2
+ END IF
+ accepted = .FALSE.
+
+ END IF
+
+ IF (this%sor_norb_cmp > 1) THEN
+ this%repetition_arr_cmp(this%sor_norb_cmp) =
this%repetition_arr_cmp(this%sor_norb_cmp) + 1
+ END IF
+
+ IF (info_verb >= 2 .AND. MOD(this%sor_ntrial_cmp,5000) == 0) THEN
+ WRITE(stdout,"(2X,A,3(I0,1X))") "Nr of accepted orbits and
trials: ", &
+ this%sor_norb_cmp, this%sor_ntrial_cmp
+ END IF
+
+ END DO
+ IF (info_verb >= 2) THEN
+ WRITE(stdout,"(A)") "END OF MCMC RANGING"
+ END IF
+
+ this%orb_arr_cmp => reallocate(this%orb_arr_cmp, this%sor_norb_cmp)
+ this%repetition_arr_cmp => reallocate(this%repetition_arr_cmp,
this%sor_norb_cmp)
+ this%res_arr_cmp => reallocate(this%res_arr_cmp, this%sor_norb_cmp,
nobs, 6)
+ this%rchi2_arr_cmp => reallocate(this%rchi2_arr_cmp, this%sor_norb_cmp)
+ this%pdf_arr_cmp => reallocate(this%pdf_arr_cmp, this%sor_norb_cmp)
+
+ IF (this%sor_norb_cmp > 0) THEN
+ 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))
+ 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)
+ WRITE(stdout,"(A,1X,I0)") "Final number of orbits:",
this%sor_norb_cmp
+ WRITE(stdout,"(A,1X,I0)") "Final number of trials:",
this%sor_ntrial_cmp
+ WRITE(stdout,"(A,1X,E7.2)") "Acceptance rate: ", &
+ REAL(this%sor_norb_cmp)/REAL(this%sor_ntrial_cmp)
+ WRITE(stdout,"(A,1X,E8.3)") "chi2 min: ", chi2min
+ END IF
+ ELSE
+ error = .TRUE.
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "No sample orbits found!", 1)
+! RETURN
+ END IF
+
+ DO i=1,nobs
+ CALL NULLIFY(obs_scoords(i))
+ CALL NULLIFY(obsy_ccoords(i))
+ END DO
+ DEALLOCATE(obs_masks, obs_coords, comp_coords, cosdec0, &
+ residuals, obs_scoords, obsy_ccoords, &
+ information_matrix_obs, stat=err)
+ IF (err /= 0) THEN
+ error = .TRUE.
+ CALL errorMessage("StochasticOrbit / MCMCRanging3", &
+ "Could not deallocate memory (20).", 1)
+ RETURN
+ END IF
+
+ END SUBROUTINE MCMCRanging3
+
!! *Description:*
!!
@@ -6666,7 +7513,7 @@
jac_sph_inv
this%res_arr_cmp(this%sor_norb_cmp,1:nobs,1:6) =
residuals(1:nobs,1:6)
this%orb_arr_cmp(this%sor_norb_cmp) = copy(orb)
- ELSE
+ ELSE
burn_in = .FALSE.
END IF
chi2_ = chi2
@@ -6700,17 +7547,24 @@
this%rchi2_arr_cmp => reallocate(this%rchi2_arr_cmp, this%sor_norb_cmp)
this%pdf_arr_cmp => reallocate(this%pdf_arr_cmp, this%sor_norb_cmp)
- 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))
- 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)
- WRITE(stdout,"(A,1X,I0)") "Final number of orbits:",
this%sor_norb_cmp
- WRITE(stdout,"(A,1X,I0)") "Final number of trials:",
this%sor_ntrial_cmp
- WRITE(stdout,"(A,1X,E7.2)") "Acceptance rate: ", &
- REAL(this%sor_norb_cmp)/REAL(this%sor_ntrial_cmp)
- WRITE(stdout,"(A,1X,E8.3)") "chi2 min: ", chi2min
+ IF (this%sor_norb_cmp > 0) THEN
+ 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))
+ 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)
+ WRITE(stdout,"(A,1X,I0)") "Final number of orbits:",
this%sor_norb_cmp
+ WRITE(stdout,"(A,1X,I0)") "Final number of trials:",
this%sor_ntrial_cmp
+ WRITE(stdout,"(A,1X,E7.2)") "Acceptance rate: ", &
+ REAL(this%sor_norb_cmp)/REAL(this%sor_ntrial_cmp)
+ WRITE(stdout,"(A,1X,E8.3)") "chi2 min: ", chi2min
+ END IF
+ ELSE
+ error = .TRUE.
+ CALL errorMessage("StochasticOrbit / MCMCRanging2", &
+ "No sample orbits found!", 1)
+ ! RETURN
END IF
DO i=1,nobs
@@ -6798,7 +7652,7 @@
"TRACE BACK (50)", 1)
STOP
END IF
- CALL constrainRangeDistributions(storb, storb%obss)
+ CALL constrainRangeDistributions(storb, storb%obss)
IF (error) THEN
CALL errorMessage("StochasticOrbit / autoMCMCRanging", &
"TRACE BACK (55)", 1)
@@ -16258,7 +17112,7 @@
!! velocity at specified epoch in two-body approximation.
!!
!! Simplified:
- !! - No random observation selection
+ !! - No random observation selection (note: is available but cannot
be used because of the multi-propagation scheme)
!!
!! Returns error.
!!
@@ -16404,6 +17258,13 @@
RETURN
END IF
+ IF (this%sor_random_obs_prm) THEN
+ error = .TRUE.
+ CALL errorMessage("StochasticOrbit / StatisticalRanging", &
+ "Use of random observation pairs not accepted.", 1)
+ RETURN
+ END IF
+
object_id = getID(this%obss)
obs_arc = getObservationalTimespan(this%obss)
@@ -17812,7 +18673,7 @@
errstr = ""
IF (err_verb >= 1) THEN
CALL matrix_print(information_matrix_elem, stderr,
errstr)
- END IF
+ END IF
errstr = ""
CYCLE
END IF
@@ -18768,7 +19629,7 @@
"RA,Dec residuals corresponding to the", COUNT(mask), &
"orbits having residuals smaller than the acceptance window
[asec]:"
DO i=1,SIZE(mask)
- IF (mask(i)) THEN
+ IF (info_verb >= 3 .AND. mask(i)) THEN
DO j=1,SIZE(residuals,dim=1)
WRITE(stdout,"(2X,2(A,1X,I0,1X),A,2(1X,F10.3))") &
"Orbit #", i, "& observation #", j, ":",
residuals(j,i,2:3)/rad_asec
@@ -18828,7 +19689,7 @@
"RA,Dec residuals corresponding to the", 10-COUNT(mask_), &
"additional orbits included [asec]:"
DO i=1,SIZE(mask)
- IF (mask(i) .AND. .NOT.mask_(i)) THEN
+ IF (info_verb >= 3 .AND. mask(i) .AND. .NOT.mask_(i)) THEN
DO j=1,SIZE(residuals,dim=1)
WRITE(stdout,"(2X,2(A,1X,I0,1X),A,2(1X,F10.3))") &
"Orbit #", i, "& observation #", j, ":",
residuals(j,i,2:3)/rad_asec
@@ -19040,14 +19901,14 @@
DEALLOCATE(sor_rho_arr_cmp, stat=err)
RETURN
END IF
- IF (info_verb >= 2) THEN
+ IF (info_verb >= 3) THEN
WRITE(stdout,"(2X,A,1X,A)") "constrainRangeDistributions:", &
"Constrained (rho1,rho2-rho1) distribution [au]: "
END IF
j = 0
***The diff for this file has been truncated for email.***
Reply all
Reply to author
Forward
0 new messages