-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathmove.py
119 lines (86 loc) · 3.98 KB
/
move.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
# -*- coding: utf-8 -*-
"""
Created on Sat Nov 28 14:52:59 2020
@author: emadg
"""
import numpy as np
from Log_Likelihood import Log_Likelihood
from cauchy_dist import cauchy_dist
from Identify import Identify
from Id2rho import Id2rho
import sys
def move(XnZn,globals_par,LogLc,xmc,zmc,xc,zc,rhoc,ARgc,ARTc,T,Kernel_Grv,Kernel_Mag,dg_obs,dT_obs):
Nnode=int(np.size(xc))
rho_salt_min = globals_par[1,0]
rho_salt_max = globals_par[1,1]
rho_base_min = globals_par[2,0]
rho_base_max = globals_par[2,1]
zn_min = globals_par[4,0]
dsalt = rho_salt_max-rho_salt_min
dbase = rho_base_max-rho_base_min
# possible_indx = np.arange(np.size(xmc))
for inode in np.arange(Nnode):
for ipar in np.arange(1,4): # 1 or 2 or 3
xp = xc.copy()
zp = zc.copy()
rhop = rhoc.copy()
if ipar == 1:
xp[inode] = cauchy_dist(xc[inode],0.1,0,1,xc[inode])
if np.isclose(xc[inode] , xp[inode])==1: continue
elif ipar == 2:
zp[inode] = cauchy_dist(zc[inode],0.1,zn_min,1,zc[inode])
if np.isclose(zc[inode] , zp[inode])==1: continue
else:
if rhoc[inode]<0:
rhop[inode] = cauchy_dist(rhoc[inode],0.02,rho_salt_min,rho_salt_max,rhoc[inode])
if np.isclose(rhoc[inode] , rhop[inode])==1: continue
elif rhoc[inode]>0:
rhop[inode] = cauchy_dist(rhoc[inode],0.02,rho_base_min,rho_base_max,rhoc[inode])
if np.isclose(rhoc[inode] , rhop[inode])==1: continue
if ipar<=2:
[identity_c, kcell] = Identify(xmc,zmc,xc[inode],zc[inode])
[identity_p, kcell] = Identify(xmc,zmc,xp[inode],zp[inode])
if identity_c != identity_p:
rhop[inode] = Id2rho(identity_p,globals_par)
rhop = rhop.astype('float32')
LogLp = Log_Likelihood(Kernel_Grv,Kernel_Mag,dg_obs,dT_obs,xp,zp,rhop,ARgc,ARTc,XnZn)[0]
MHP = np.exp((LogLp - LogLc)/T)
if np.random.rand()<=MHP:
LogLc = LogLp
xc = xp.copy()
zc = zp.copy()
rhoc = rhop.copy()
### Hyper Parameters (Parent Nodes)
for im in np.arange(2): # 0 for x and 1 for z
for inode in np.arange(3): # 0 or 1 or 2
rhop = rhoc.copy()
xmp = xmc.copy()
zmp = zmc.copy()
if im==0: # x
xmp[inode] = cauchy_dist(xmc[inode],0.1,0,1,xmc[inode])
if np.isclose(xmc[inode] , xmp[inode])==1: continue
elif im==1: # z
zmp[inode] = cauchy_dist(zmc[inode],0.1,zn_min,1,zmc[inode])
if np.isclose(zmc[inode] , zmp[inode])==1: continue
# else:
# indxp = np.delete(possible_indx, inode).copy()
# inodep = np.random.choice(indxp,1)
# xmp[inode] = xmc[inodep].copy()
# zmp[inode] = zmc[inodep].copy()
# xmp[inodep] = xmc[inode].copy()
# zmp[inodep] = zmc[inode].copy()
[identity_c, kcell] = Identify(xmc,zmc,xc,zc)
[identity_p, kcell] = Identify(xmp,zmp,xc,zc)
iden_diff = identity_p!=identity_c
if iden_diff.any() == True:
rhop = Id2rho(identity_p,globals_par)
rhop = (1-iden_diff)*rhoc + iden_diff*rhop
rhop = rhop.astype('float32')
LogLp = Log_Likelihood(Kernel_Grv,Kernel_Mag,dg_obs,dT_obs,xc,zc,rhop,ARgc,ARTc,XnZn)[0]
MHP = np.exp((LogLp - LogLc)/T)
if np.random.rand()<=MHP:
LogLc = LogLp
rhoc = rhop.copy()
xmc = xmp.copy()
zmc = zmp.copy()
return [LogLc,xmc,zmc,xc,zc,rhoc]