Skip to content

Commit

Permalink
Changed to adaptive loss function
Browse files Browse the repository at this point in the history 
The shape value is now broken into an optimized parameter for inside experiments and between experiments
  • Loading branch information
@tsikes
tsikes committed Aug 29, 2021
1 parent 60e782f commit 4702ef3
Show file tree
Hide file tree
Showing 12 changed files with 548,014 additions and 14,452 deletions.
163 changes: 111 additions & 52 deletions Development/Random Testing/loss func integral/loss_integral.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,8 +2,15 @@
# and licensed under BSD-3-Clause. See License.txt in the top-level
# directory for license and copyright information.

import pygmo
import numpy as np
from scipy import integrate, optimize # used to integrate weights numerically
from scipy import interpolate, integrate, optimize # used to integrate weights numerically
from scipy.interpolate import CubicSpline
import matplotlib.pyplot as plt
import multiprocessing as mp

min_pos_system_value = (np.finfo(float).tiny*(1E20))**(1/2)
max_pos_system_value = (np.finfo(float).max*(1E-20))**(1/2)

def generalized_loss_fcn(x, mu=0, a=2, c=1): # defaults to L2 loss
x_c_2 = ((x-mu)/c)**2
Expand All @@ -23,76 +30,128 @@ def generalized_loss_fcn(x, mu=0, a=2, c=1): # defaults to L2 loss

return loss*c**a + mu # multiplying by c^2 is not necessary, but makes order appropriate

def triple_sigmoid(x, *var):
A1, A2, p1, p2, p3, logx01, logx02, logx03, logx04, h1, h2, h3, h4 = var
sig1 = 1/(1 + 10**((logx01 - x)*h1))
sig2 = 1/(1 + 10**((logx02 - x)*h2))
sig3 = 1/(1 + 10**((logx03 - x)*h3))
sig4 = 1/(1 + 10**((logx04 - x)*h4))
y = A1 + (A2 - A1)*(p1*sig1 + p2*sig2 + p3*sig3 + (1 - p1 - p2 - p3)*sig4)
def sig(x): # Numerically stable sigmoid function
eval = np.empty_like(x)
pos_val_f = np.exp(-x[x >= 0])
eval[x >= 0] = 1/(1 + pos_val_f)
neg_val_f = np.exp(x[x < 0])
eval[x < 0] = neg_val_f/(1 + neg_val_f)

# clip so they can be multiplied
eval[eval > 0] = np.clip(eval[eval > 0], min_pos_system_value, max_pos_system_value)
eval[eval < 0] = np.clip(eval[eval < 0], -max_pos_system_value, -min_pos_system_value)

return eval

def sig_deriv(sig_val):
return sig_val*(1-sig_val)

def eval_coef(tau, a):
return a[0] + a[1]*np.exp(a[2]*tau)

def Z_fun(tau, alpha, *vars):
x = np.log(-alpha + 2 + 1/1001)

L_x0 = eval_coef(tau, vars[0:3])
k = eval_coef(tau, vars[3:6])
A1 = [eval_coef(tau, vars[6:9]), eval_coef(tau, vars[9:12])]
A2 = [eval_coef(tau, vars[12:15]), eval_coef(tau, vars[15:18])]
p = [eval_coef(tau, vars[18:21]), eval_coef(tau, vars[21:24]), eval_coef(tau, vars[24:27])]
x0 = [eval_coef(tau, vars[27:30]), eval_coef(tau, vars[30:33]), eval_coef(tau, vars[33:36]),
eval_coef(tau, vars[36:39]), eval_coef(tau, vars[39:42])]
h = [eval_coef(tau, vars[42:45]), eval_coef(tau, vars[45:48]), eval_coef(tau, vars[48:51]),
eval_coef(tau, vars[51:54]), eval_coef(tau, vars[54:57])]

p = [p[0], 1-p[0], p[1], p[2], 1-p[1]-p[2]]

L = sig((x - L_x0)/k)
f_1 = A1[0]+(A2[0]-A1[0])*(p[0]/(1+10**((x0[0]-x)*h[0])) + p[1]/(1+10**((x0[1]-x)*h[1]))+ p[2]/(1+10**((x0[2]-x)*h[2])))
f_2 = A1[1]+(A2[1]-A1[1])*(p[3]/(1+10**((x0[3]-x)*h[3])) + p[4]/(1+10**((x0[4]-x)*h[4])))

# A1, A2, p1, p2, logx01, logx02, logx03, h1, h2, h3 = var
# sig1 = 1/(1 + 10**((logx01 - x)*h1))
# sig2 = 1/(1 + 10**((logx02 - x)*h2))
# sig3 = 1/(1 + 10**((logx03 - x)*h3))
# y = A1 + (A2 - A1)*(p1*sig1 + p2*sig2 + (1 - p1 - p2)*sig3)
return np.exp(L*f_2 + (1-L)*f_1)

return y
class pygmo_objective_fcn:
def __init__(self, bnds, tau, alpha, Z):
self.bnds = bnds
self.tau = np.unique(tau)
self.alpha = np.unique(alpha)
self.Z = Z
self.ln_Z = np.log(Z)

def fitness(self, x):
obj = 0
for i in range(len(self.tau)):
# print(self.tau[i], self.alpha[i], Z_fun(self.tau[i], self.alpha[i], *x))
try:
Z_fit = Z_fun(self.tau[i], self.alpha[i], *x)
if np.isnumeric(Z_fit):
obj += 100*(np.log(Z_fit) - self.ln_Z[i])**2
print(obj)
else:
obj += 1
except:
obj += 1

return [obj]

def get_bounds(self):
return self.bnds

def gradient(self, x):
return pygmo.estimate_gradient_h(lambda x: self.fitness(x), x)


'''
alpha_vec = np.geomspace(-1000, -1, 10000)
alpha_vec = np.concatenate([alpha_vec, np.linspace(-1, 2, int(3/(alpha_vec[-1] - alpha_vec[-2])-1))[1:]])
tau = 100
alpha_vec = np.linspace(-1000, -10, 991)
alpha_vec = np.concatenate([alpha_vec, np.linspace(-10, 2, 1201)[1:]])
tau = 250
res = []
c = 1
for alpha in alpha_vec:
int_func = lambda x: np.exp(-generalized_loss_fcn(x, a=alpha, c=1))
Z, err = integrate.quadrature(int_func, -tau, tau, tol=1E-14, maxiter=10000)
for tau in np.linspace(1, 250, 250):
for alpha in alpha_vec:
int_func = lambda x: np.exp(-generalized_loss_fcn(x, a=alpha, c=1))
Z, err = integrate.quadrature(int_func, -tau, tau, tol=1E-14, maxiter=10000)
res.append([tau, c, alpha, Z])
print(f'{alpha:0.3f} {c:0.3f} {Z:0.14e}')
res.append([tau, c, alpha, Z])
print(f'{tau:0.3f} {alpha:0.3f} {Z:0.14e}')
res = np.array(res)
np.savetxt('res.csv', res, delimiter=',')
'''

data = np.genfromtxt('res.csv', delimiter=',')

alpha = data[:,2]
Z = data[:,3]

# p0 = [2.5, 73, 0.1, 0.3, -2, -1.5, -1, -0.05, -0.2, -0.5]
# bnds = [[-100, 10, 1E-6, 1E-6, -100, -100, -100, -100, -100, -100],
# [10, 1000, 0.999999, 0.999999, 100, 100, 100, 100, 100, 100]]

# popt = -2.59255522e+01 7.41834165e+01 4.53316795e-01 2.56195379e-01 9.14556135e+00 -4.24683064e+00 -1.41210936e+00 -2.55275337e-02 -2.94402525e-01 -9.28894263e-01

# A1, A2, p1, p2, p3, logx01, logx02, logx03, logx04, h1, h2, h3, h4
p0 = [-26, 74, 0.45, 0.25, 0.1, 9, -4, -1.4, -1.0, -2.5E-2, -2.9E-1, -9.3E-1, -0.2]
data = data[:, [0,2,3]]
tau = data[:,0].reshape((250,2191))
alpha = data[:,1].reshape((250,2191))
Z = data[:,2].reshape((250,2191))

p0 = [-1.00000000e+02, 7.44435255e+01, 6.47123371e-01, 9.50095867e-02,
1.26494317e-01, 7.50612556e+01, -7.36539464e+00, -1.16617872e+00,
-2.94375758e+00, -1.31391561e-02, -1.55103228e-01, -1.18310054e+00,
-5.20960057e-01]
tau = tau[:,1:]
alpha = alpha[:,1:]
Z = Z[:,1:]

# p0 = [2.5, 73, 0.1, 0.1, 0.1, -2, -1.5, -1, -10, -0.05, -0.2, -0.5, -0.2]
bnds = [[-1E6, 10, 1E-6, 1E-6, 1E-6, -100, -100, -100, -100, -100, -100, -100, -100],
[10, 1000, 0.999999, 0.999999, 0.999999, 100, 100, 100, 100, 100, 100, 100, 100]]
kx = 3
ky = kx

popt, _ = optimize.curve_fit(triple_sigmoid, alpha, Z, p0=p0, method='dogbox', bounds=bnds, x_scale='jac', max_nfev=len(p0)*1000)
# spline_fit = interpolate.RectBivariateSpline(tau[:,0], alpha[0,:], np.log(Z), kx=kx, ky=ky, s=0.0001)

print(popt)
for i, a in enumerate(alpha):
print(a, Z[i], triple_sigmoid(a, *popt))
# with open("spline_save.csv", "ab") as f:
# for i in range(3):
# np.savetxt(f, np.matrix(spline_fit.tck[i]), delimiter=",")
# for val in [kx, ky]:
# np.savetxt(f, np.matrix(val), delimiter=",")

# t1 = get_exponent_cuv(x_y_curve, &y0, &A1, &sign);
# t1 = t2 = -1 / t1;
# t1 = t1 * 0.9;
# t2 = t2 * 1.1;
# A1 = A2 = sign * exp(A1) / 2;
tck = []
with open("spline_save.csv") as f:
for i in range(5):
tck.append(np.array(f.readline().split(','), dtype=float))

# x = C
# y = int_func
spline_fit = interpolate.RectBivariateSpline._from_tck(tck)

idx = np.argwhere(data[:,0] == 1)

# plt.plot(x_deriv, deriv(x_deriv))
# plt.plot(data[idx,1], np.log(data[idx,2]) - spline_fit([2], data[idx,1]).T)
# plt.plot(data[idx,1][1:], np.log(data[idx,2][1:]) - spline_fit([2], data[idx,1][1:]).T)
plt.plot(data[idx,1][1:], np.log(data[idx,2][1:]), data[idx,1][1:], spline_fit([1], data[idx,1][1:]).T)
plt.show()
Loading

0 comments on commit 4702ef3

Please sign in to comment.