Skip to content

Commit

Permalink
[BUG] Fix bugs
Browse files Browse the repository at this point in the history
  • Loading branch information
@suzuki-shm
suzuki-shm committed Jul 23, 2019
1 parent a90d085 commit 22a6b7f
Show file tree
Hide file tree
Showing 3 changed files with 137 additions and 71 deletions.
156 changes: 88 additions & 68 deletions sphere/sphere_mcplot.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -203,44 +203,54 @@ def normalized_constraint_se(kappa, psi, lambda_):

# inverse Batschelet transformation
# (or inverse symmetric extended transformation)
def inv_batschelet_trans_se_newton(theta, loc, lambda_):
t = theta
count = 0
err = 1.0
while(err > np.finfo(float).eps):
f = t - (1.0+lambda_) * np.sin(t-loc) / 2.0 - theta
fd = 1.0 - (1.0+lambda_) * np.cos(t-loc) / 2.0
tp = calc_tp(t, f, fd)
t = tp
err = np.abs(f).max()
count += 1
if count == 30:
break
return t


def inv_batschelet_trans_se_bisection(theta, loc, lambda_):
dsize = (lambda_.shape[0], theta.size)
t1 = np.ones(dsize) * -2.0 * np.pi
t2 = np.ones(dsize) * 2.0 * np.pi
count = 0
err = 1.0
while(err > np.finfo(float).eps):
t = (t1 + t2) / 2.0
f = t - (1.0+lambda_) * np.sin(t-loc) / 2.0 - theta
t1[f < 0] = t[f < 0]
t2[f >= 0] = t[f >= 0]
err = np.abs(f).max()
count += 1
if count == 100:
break
return t


def inv_trans_se(theta, loc, lambda_):
def inv_batschelet_trans_se_newton(theta, loc, lambda_):
t = theta
count = 0
err = 1.0
while(err > np.finfo(float).eps):
f = t - (1.0+lambda_) * np.sin(t-loc) / 2.0 - theta
fd = 1.0 - (1.0+lambda_) * np.cos(t-loc) / 2.0
tp = calc_tp(t, f, fd)
t = tp
err = np.abs(f).max()
count += 1
if count == 30:
break
return t
if type(lambda_) is list:
la = np.array(lambda_)
elif type(lambda_) is int or type(lambda_) is float:
la = np.array([lambda_])
else:
la = lambda_

def inv_batschelet_trans_se_bisection(theta, loc, lambda_):
dsize = (lambda_.shape[0], theta.size)
t1 = np.ones(dsize) * -2.0 * np.pi
t2 = np.ones(dsize) * 2.0 * np.pi
count = 0
err = 1.0
while(err > np.finfo(float).eps):
t = (t1 + t2) / 2.0
f = t - (1.0+lambda_) * np.sin(t-loc) / 2.0 - theta
t1[f < 0] = t[f < 0]
t2[f >= 0] = t[f >= 0]
err = np.abs(f).max()
count += 1
if count == 100:
break
return t
if np.abs(lambda_).max() < 0.8:
t = inv_batschelet_trans_se_newton(theta, loc, lambda_)
t = inv_batschelet_trans_se_newton(theta, loc, la)
else:
t = inv_batschelet_trans_se_bisection(theta, loc, lambda_)
return ((1.0 - lambda_) / (1.0 + lambda_) * theta +
2.0 * lambda_ / (1.0 + lambda_) * t)
t = inv_batschelet_trans_se_bisection(theta, loc, la)
return ((1.0 - la) / (1.0 + la) * theta +
2.0 * la / (1.0 + la) * t)


def invsevonmises_pdf(theta, loc, kappa, lambda_):
Expand Down Expand Up @@ -297,48 +307,58 @@ def miaevonmises_pdf(theta, loc, kappa, nu):


# inverse mode invariance asymmetry extended transformation
def inv_trans_sin2(theta, loc, nu):
def inv_trans_sin2_newton(theta, loc, nu):
t = theta
f = t + nu * np.sin(t-loc)**2.0 - theta
def inv_trans_sin2_newton(theta, loc, nu):
t = theta
f = t + nu * np.sin(t-loc)**2.0 - theta
fd = 1.0 + 2.0 * nu * np.sin(t-loc) * np.cos(t-loc)
tp = calc_tp(t, f, fd)
t = tp
err = np.abs(f).max()
count = 0
while(err > np.finfo(float).eps):
f = t + nu * np.sin(t-loc)**2 - theta
fd = 1.0 + 2.0 * nu * np.sin(t-loc) * np.cos(t-loc)
tp = calc_tp(t, f, fd)
t = tp
err = np.abs(f).max()
count = 0
while(err > np.finfo(float).eps):
f = t + nu * np.sin(t-loc)**2 - theta
fd = 1.0 + 2.0 * nu * np.sin(t-loc) * np.cos(t-loc)
tp = calc_tp(t, f, fd)
t = tp
err = np.abs(f).max()
count += 1
if count == 30:
break
return t

def inv_trans_sin2_bisection(theta, loc, nu):
dsize = (nu.shape[0], theta.size)
t1 = np.ones(dsize) * -2.0 * np.pi
t2 = np.ones(dsize) * 2.0 * np.pi
count += 1
if count == 30:
break
return t


def inv_trans_sin2_bisection(theta, loc, nu):
dsize = (nu.shape[0], theta.size)
t1 = np.ones(dsize) * -2.0 * np.pi
t2 = np.ones(dsize) * 2.0 * np.pi
t = (t1 + t2) / 2.0
f = t + nu * np.sin(t-loc)**2 - theta
err = np.abs(f).max()
count = 0
while(err > np.finfo(float).eps):
t = (t1 + t2) / 2.0
f = t + nu * np.sin(t-loc)**2 - theta
f = t + nu * np.sin(t-loc)**2.0 - theta
t1[f < 0] = t[f < 0]
t2[f >= 0] = t[f >= 0]
err = np.abs(f).max()
count = 0
while(err > np.finfo(float).eps):
t = (t1 + t2) / 2.0
f = t + nu * np.sin(t-loc)**2.0 - theta
t1[f < 0] = t[f < 0]
t2[f >= 0] = t[f >= 0]
err = np.abs(f).max()
count += 1
if count == 100:
break
return t
count += 1
if count == 100:
break
return t


def inv_trans_sin2(theta, loc, nu):
if type(nu) is list:
n = np.array(nu)
elif type(nu) is int or type(nu) is float:
n = np.array([nu])
else:
n = nu

if np.abs(nu).max() < 0.8:
return inv_trans_sin2_newton(theta, loc, nu)
return inv_trans_sin2_newton(theta, loc, n)
else:
return inv_trans_sin2_bisection(theta, loc, nu)
return inv_trans_sin2_bisection(theta, loc, n)


def invmiaecardioid_pdf(theta, loc, rho, nu):
Expand Down
19 changes: 16 additions & 3 deletions sphere/sphere_waic.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@
from sphere.stan_utils import load_log_lik
from sphere.stan_utils import get_waic
from sphere.sphere_utils import get_logger
from sphere.sphere_utils import load_multiple_depth_file


def argument_parse(argv=None):
Expand All @@ -24,25 +25,37 @@ def argument_parse(argv=None):
default="both",
choices=["original", "bda3", "both"],
type=str)
parser.add_argument("-d", "--depth",
dest="d",
nargs="*",
default=None,
type=str)
args = parser.parse_args(argv)
return vars(args)


def main(args, logger):
results = []
if args["d"] is not None:
ddf = load_multiple_depth_file(args["d"])
ddf = ddf[ddf["depth"] != 0]
depth = ddf["depth"].values
else:
depth = None

for f in args["log_lik_files"]:
name = os.path.basename(f)
log_lik = load_log_lik(f)
if args["t"] != "both":
waic = get_waic(log_lik, args["t"])
waic = get_waic(log_lik, depth, args["t"])
results.append({
"file_name": name,
"file_path": f,
"waic_{0}".format(args["t"]): waic
})
else:
waic_original = get_waic(log_lik, "original")
waic_bda3 = get_waic(log_lik, "bda3")
waic_original = get_waic(log_lik, depth, "original")
waic_bda3 = get_waic(log_lik, depth, "bda3")
results.append({
"file_name": name,
"file_path": f,
Expand Down
33 changes: 33 additions & 0 deletions sphere/tests/test_sphere_mcplot.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@

import unittest
import os
import numpy as np
from sphere import sphere_mcplot
from sphere import sphere_estimate
from sphere.sphere_utils import get_logger
Expand Down Expand Up @@ -625,6 +626,38 @@ def test_sphere_mcplot_command_invmiaevm_opt(self):
args = sphere_mcplot.argument_parse(argv)
sphere_mcplot.main(args, SphereMcplotTest.logger)

def test_inv_trans_sin2_newton(self):
t = np.linspace(-np.pi, np.pi, 10)
sphere_mcplot.inv_trans_sin2(t, loc=0, nu=0)

def test_inv_trans_sin2_bisection(self):
t = np.linspace(-np.pi, np.pi, 10)
sphere_mcplot.inv_trans_sin2(t, loc=0, nu=1.0)

def test_inv_trans_sin2_newton_list(self):
t = np.linspace(-np.pi, np.pi, 10)
sphere_mcplot.inv_trans_sin2(t, loc=0, nu=[0])

def test_inv_trans_sin2_bisection_list(self):
t = np.linspace(-np.pi, np.pi, 10)
sphere_mcplot.inv_trans_sin2(t, loc=0, nu=[1.0])

def test_inv_trans_sin2_newton_1darray(self):
t = np.linspace(-np.pi, np.pi, 10)
sphere_mcplot.inv_trans_sin2(t, loc=0, nu=np.array([0]))

def test_inv_trans_sin2_bisection_1darray(self):
t = np.linspace(-np.pi, np.pi, 10)
sphere_mcplot.inv_trans_sin2(t, loc=0, nu=np.array([1.0]))

def test_inv_trans_sin2_newton_2darray(self):
t = np.linspace(-np.pi, np.pi, 10)
sphere_mcplot.inv_trans_sin2(t, loc=0, nu=np.array([[0]]))

def test_inv_trans_sin2_bisection_2darray(self):
t = np.linspace(-np.pi, np.pi, 10)
sphere_mcplot.inv_trans_sin2(t, loc=0, nu=np.array([[1.0]]))


if __name__ == '__main__':
unittest.main()

0 comments on commit 22a6b7f

Please sign in to comment.