From 897a516f8c2f9105b82abced32978e2ff37c15ca Mon Sep 17 00:00:00 2001 From: BANG XIANG YONG Date: Tue, 27 Jul 2021 10:35:27 +0100 Subject: [PATCH 01/43] placeholder for noise jitter example --- agentMET4FOF_tutorials/noise_jitter_removal/__init__.py | 0 .../noise_jitter_removal/noise_jitter_removal.py | 0 2 files changed, 0 insertions(+), 0 deletions(-) create mode 100644 agentMET4FOF_tutorials/noise_jitter_removal/__init__.py create mode 100644 agentMET4FOF_tutorials/noise_jitter_removal/noise_jitter_removal.py diff --git a/agentMET4FOF_tutorials/noise_jitter_removal/__init__.py b/agentMET4FOF_tutorials/noise_jitter_removal/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/agentMET4FOF_tutorials/noise_jitter_removal/noise_jitter_removal.py b/agentMET4FOF_tutorials/noise_jitter_removal/noise_jitter_removal.py new file mode 100644 index 00000000..e69de29b From 045698fe5b317ba8805971d69b59c7c0403224dc Mon Sep 17 00:00:00 2001 From: BANG XIANG YONG Date: Tue, 27 Jul 2021 10:42:34 +0100 Subject: [PATCH 02/43] (wip) njremoval : first example on launching NJRemoval Agent --- .../noise_jitter_removal.py | 85 +++++++++++++++++++ 1 file changed, 85 insertions(+) diff --git a/agentMET4FOF_tutorials/noise_jitter_removal/noise_jitter_removal.py b/agentMET4FOF_tutorials/noise_jitter_removal/noise_jitter_removal.py index e69de29b..3f412ebc 100644 --- a/agentMET4FOF_tutorials/noise_jitter_removal/noise_jitter_removal.py +++ b/agentMET4FOF_tutorials/noise_jitter_removal/noise_jitter_removal.py @@ -0,0 +1,85 @@ +from agentMET4FOF.agents import AgentMET4FOF, AgentNetwork, MonitorAgent +# from agentMET4FOF.streams import SineGeneratorJitter +from Sinegen import SineGeneratorJitter + +import numpy as np +from NJRemove.NJRemoval_class_withmcmc import MCMCMH_NJ + +from agentMET4FOF.agents.signal_agents import StaticSineGeneratorWithJitterAgent, NoiseAgent + + +def njr(fs, ydata, N, niter, tol, m0w, s0w, m0t, s0t, Mc, M0, Nc, Q): + analyse_fun = MCMCMH_NJ(fs, ydata, N, niter, tol, m0w, s0w, m0t, s0t, Mc, M0, Nc, Q) + yhat1= analyse_fun.AnalyseSignalN() + return yhat1 + + +######################################## +class NJRemoved(AgentMET4FOF): + def init_parameters(self, fs=100, ydata = np.array([]), N=15, niter=100, tol=1e-9, m0w = 10, s0w = 0.0005, m0t = 10, s0t = 0.0002*100/8, Mc=5000, M0=100, Nc=100, Q=50 ): + self.fs = fs + self.ydata = ydata + self.N = N + self.niter = niter + self.tol = tol + self.m0w = m0w + self.s0w = s0w + self.m0t = m0t + self.s0t = s0t + self.Mc = Mc + self.M0 = M0 + self.Nc = Nc + self.Q = Q + + + def on_received_message(self, message): + ddata = message['data'] + self.ydata = np.append(self.ydata, ddata) + if np.size(self.ydata) == self.N: + t = njr(self.fs, self.ydata, self.N, self.niter, self.tol, self.m0w, self.s0w, self.m0t, self.s0t, self.Mc, self.M0,self.Nc, self.Q) + self.send_output(self.ydata[7] - t) + self.ydata = self.ydata[1:self.N] + + +class SineGeneratorAgent(AgentMET4FOF): + def init_parameters(self): + self.stream = SineGeneratorJitter(jittersd=0.0005, noisesd=0.0002) + + def agent_loop(self): + if self.current_state == "Running": + sine_data = self.stream.next_sample() # dictionary + self.send_output(sine_data['quantities']) + + +def main(): + # start agent network server + agentNetwork = AgentNetwork(backend="mesa") + # init agents + + jitter_gen_agent = agentNetwork.add_agent(agentType=StaticSineGeneratorWithJitterAgent) + + noise_agent = agentNetwork.add_agent(agentType=NoiseAgent) + + njremove_agent = agentNetwork.add_agent(agentType=NJRemoved) + + monitor_agent = agentNetwork.add_agent(agentType=MonitorAgent) + monitor_agent2 = agentNetwork.add_agent(agentType=MonitorAgent) + + + # connect agents : jitter generator -> noise -> njremoval agent + agentNetwork.bind_agents(jitter_gen_agent, noise_agent) + agentNetwork.bind_agents(noise_agent, njremove_agent) + + # connect monitor agents + agentNetwork.bind_agents(jitter_gen_agent, monitor_agent) + agentNetwork.bind_agents(njremove_agent, monitor_agent2) + + # set all agents states to "Running" + agentNetwork.set_running_state() + + # allow for shutting down the network after execution + return agentNetwork + + +if __name__ == '__main__': + main() \ No newline at end of file From 769aa146b903bd636e733a985af8c6df0d59d160 Mon Sep 17 00:00:00 2001 From: BANG XIANG YONG Date: Tue, 27 Jul 2021 10:48:00 +0100 Subject: [PATCH 03/43] wip (njremoval) : clean up unnecessary imports and classes --- .../noise_jitter_removal/noise_jitter_removal.py | 14 -------------- 1 file changed, 14 deletions(-) diff --git a/agentMET4FOF_tutorials/noise_jitter_removal/noise_jitter_removal.py b/agentMET4FOF_tutorials/noise_jitter_removal/noise_jitter_removal.py index 3f412ebc..a8e84aad 100644 --- a/agentMET4FOF_tutorials/noise_jitter_removal/noise_jitter_removal.py +++ b/agentMET4FOF_tutorials/noise_jitter_removal/noise_jitter_removal.py @@ -1,7 +1,4 @@ from agentMET4FOF.agents import AgentMET4FOF, AgentNetwork, MonitorAgent -# from agentMET4FOF.streams import SineGeneratorJitter -from Sinegen import SineGeneratorJitter - import numpy as np from NJRemove.NJRemoval_class_withmcmc import MCMCMH_NJ @@ -40,17 +37,6 @@ def on_received_message(self, message): self.send_output(self.ydata[7] - t) self.ydata = self.ydata[1:self.N] - -class SineGeneratorAgent(AgentMET4FOF): - def init_parameters(self): - self.stream = SineGeneratorJitter(jittersd=0.0005, noisesd=0.0002) - - def agent_loop(self): - if self.current_state == "Running": - sine_data = self.stream.next_sample() # dictionary - self.send_output(sine_data['quantities']) - - def main(): # start agent network server agentNetwork = AgentNetwork(backend="mesa") From e195d192a47e2503155e2fd0eab0c75e98f2b7c8 Mon Sep 17 00:00:00 2001 From: Anupam Prasad Vedurmudi Date: Tue, 27 Jul 2021 12:20:06 +0200 Subject: [PATCH 04/43] wip(noise_jitter_removal_agent): methods and classes copied over from npl repo --- .../agents/noise_jtter_removal_agents.py | 868 ++++++++++++++++++ 1 file changed, 868 insertions(+) create mode 100644 agentMET4FOF/agents/noise_jtter_removal_agents.py diff --git a/agentMET4FOF/agents/noise_jtter_removal_agents.py b/agentMET4FOF/agents/noise_jtter_removal_agents.py new file mode 100644 index 00000000..67cbc9d0 --- /dev/null +++ b/agentMET4FOF/agents/noise_jtter_removal_agents.py @@ -0,0 +1,868 @@ +import math + +import numpy as np +from scipy.optimize import minimize + +from agentMET4FOF.agents import AgentMET4FOF, AgentNetwork, MonitorAgent +from agentMET4FOF.streams.signal_streams import StaticSineGeneratorWithJitter + +"""### mcmcci: MCMC convergence indices for multiple chains.""" + +def mcmcci(A,M0): + ''' + mcmcci: MCMC convergence indices for multiple chains. + ---------------------------------------------------------------------------- + KJ, LRW, PMH + + Version 2020-04-22 + + ---------------------------------------------------------------------------- + Inputs: + A(M, N): Chain samples, N chains of length M. + + M0: Length of burn-in, M > M0 >= 0. + + Outputs: + Rhat: Convergence index. Rhat is expected to be greater than 1. The closer + Rhat is to 1, the better the convergence. + + Neff: Estimate of the effective number of independent draws. Neff is + expected to be less than (M-M0)*N. + ---------------------------------------------------------------------------- + Note: If the calculated value of Rhat is < 1, then Rhat is set to 1 and Neff + set to (M-M0)*N, their limit values. + + Note: If N = 1 or M0 > M-2, Rhat = 0; Neff = 0. + ''' + A = np.array(A) + + M, N = A.shape + + # Initialisation + Rhat = 0 + Neff = 0 + + # The convergence statistics can only be evaluated if there are multiple chains + # and the chain length is greater than the burn in length + + if N > 1 and M > M0 + 1: + Mt = M - M0 + + # Chain mean and mean of means + asub = A[M0:,:] + ad = np.mean(asub,axis = 0) + add = asub.mean() + + # Within group standard deviation + ss = np.std(asub,axis = 0) + + # Between groups variance. + dd = np.square(ad - add) + B = (Mt*np.sum(dd))/(N-1) + + # Within groups variance. + W = np.sum(np.square(ss))/N + + # V plus + Vp = (Mt-1)*W/Mt + B/Mt + + # Convergence statistic, effective number of independent samples + Rhat = np.sqrt(Vp/W) + Neff = Mt*N*Vp/B + + Rhat = np.maximum(Rhat,1) + Neff = np.minimum(Neff,Mt*N) + + return Rhat, Neff + +"""### mcsums: Summary information from MC samples.""" + +def mcsums(A,M0,Q): + ''' + mcsums: Summary information from MC samples. + ------------------------------------------------------------------------- + KJ, LRW, PMH + Version 2020-04-22 + + ------------------------------------------------------------------------- + Inputs: + A(M,N): An array that stores samples of size M x N. + + M0: Burn-in period with M > M0 >= 0. + + Q(nQ,1): Percentiles specifications, 0 <= Q(l) <= 100. + + Outputs: + abar(n,1): Mean for each sample. + + s(n,1): Standard deviation for sample. + + aQ(nQ,n): Percentiles corresponding to Q. + + ''' + # Size of samples after burn-in + A = np.array(A) + + M, N = A.shape + + m = (M - M0)*N + + # Initialise percentile vector + nQ = Q.size + aQ = np.zeros(nQ) + + # Samples from N parallel chains after burn-in period + aaj = A[M0:, :] + aaj = aaj.flatten() + + # Mean and standard deviation of samples + abar = np.mean(aaj) + s = np.std(aaj) + + # Percentiles of samples + aQ = np.percentile(aaj,Q) + + return abar, s, aQ + +"""### jumprwg: Jumping distribution for the Metropolis Hastings Gaussian random walk algorithm""" + +def jumprwg(A, L): + ''' + jumprwg: Jumping distribution for the Metropolis Hastings Gaussian random + walk algorithm + ------------------------------------------------------------------------- + KJ, LRW, PMH + Version 2020-04-22 + ------------------------------------------------------------------------- + Inputs: + A(n,N): Samples at the current iteration + + L(n,n): Cholesky factor of variance of parameter vector. + + Outputs: + As(n,N): Proposed parameter array which is randomly sampled from the + jumping distribution + + dp0: The difference between the logarithm of the jumping distribution + associated with moving from A(:,j) to As(:,j) and that associated with + moving from As(:,j) to A(:,j), up to an additive constant. + log P0(a|as) - log P0(as|a) + + ''' + # Number of parameters and parallel chains + n, N = A.shape + + # random draw from a Gaussian distribution + e = np.random.normal(0, 1, size=(n,N)) + + # proposed draw from a Gaussian distribution with mean A and variance LL' + As = A + np.matmul(L,e) + + # For a Gaussian random walk, since log P0(a|as) = log P0(as|a), dp0 will always be zero + dp0 = np.zeros(N) + + return As, dp0 + +"""### Cubic function and its first and second derivative""" + +def fgh_cubic(alpha,t): + ''' + ------------------------------------------------------------------------- + Cubic function and its first and second derivative + ------------------------------------------------------------------------- + KJ, LRW, PMH + Version 2020-04-22 + ------------------------------------------------------------------------- + Inputs: + alpha(4,N): Alpha parameters + + t(m,1): Times + + Outputs: + f(m,N): Cubic function + + f1(m,N): Derivative of cubic + + f2(m,N): Second derivative of cubic + ''' + + # length of data and number of paramaters + m = t.size + + # design matrix + C = np.array([np.ones(m), t, t**2, t**3]) + + # derivate info + C1 = np.array([np.ones(m), 2*t, 3*t**2]) + C2 = np.array([2*np.ones(m), 6*t]) + + # cubic and derivatives + f = np.matmul(C.T,alpha) + f1 = np.matmul(C1.T,alpha[1:]) + f2 = np.matmul(C2.T,alpha[2:]) + + return f, f1, f2 + +"""### Log of the gaussian pdf""" + +def ln_gauss_pdf_v(x,mu,sigma): + ''' + ------------------------------------------------------------------------- + Log of the Gaussian pdf + ------------------------------------------------------------------------- + KJ, LRW, PMH + Version 2020-03-12 + -------------------------------------------------------------------------- + Inputs: + x(m,1): Points at which pdf is to be evaluated + + mu: Mean of distribution + + sigma: Standard deviation of the distribution + + Output: + logf: Log of the Gaussian pdf at x with mean mu and + std sigma + ''' + + + try: + # When inputs are high dimensional arrays/matrices + xx = np.matlib.repmat(x,mu.shape[1],1) + xx = xx.T + logk = - np.log(2*math.pi)/2 - np.log(sigma) + logp = -((xx - mu)**2)/(2*sigma**2) + # Log of the Gaussian PDF + logf = logk + logp + + except IndexError: + # When inputs are vectors + logk = - np.log(2*math.pi)/2 - np.log(sigma) + logp = -((x - mu)**2)/(2*sigma**2) + # Log of the Gaussian PDF + logf = logk + logp + + return logf + +"""### Target dist for noise and jitter posterior dist""" + +def tar_at(at, y, x, m0w, s0w, m0t, s0t): + ''' + ------------------------------------------------------------------------- + Target dist for noise and jitter posterior dist + ------------------------------------------------------------------------- + KJ, LRW, PMH + Version 2020-04-22 + -------------------------------------------------------------------------- + Inputs: + at(n+2,N): Parameters alpha, log(1/tau^2) and log(1/w^2) + + y(m,1): Signal + + x(m,1): time at which signal was recorded + + s0w and s0t: prior estimates of tau and omega + + m0w and m0t: degree of belief in prior estimates for tau and omega + + Output: + T: Log of the posterior distribution + ''' + + # Size of parameter vector + at = np.array(at) + p = at.shape[0] + + # Number of alphas + n = p - 2 + + # Extract parameters + alpha = at[0:n] + phi1 = np.exp(at[-2]) + phi2 = np.exp(at[-1]) + taus = np.ones(phi1.shape)/np.sqrt(phi1) + omegas = np.ones(phi2.shape)/np.sqrt(phi2) + + # Gamma priors for phis + prior_phi1 = (m0t/2)*np.log(phi1) - phi1*m0t*s0t**2/2 + prior_phi2 = (m0w/2)*np.log(phi2) - phi2*m0w*s0w**2/2 + + + # function that evaluates the cubic function with user specified cubic parameters + fun = lambda aa: fgh_cubic(aa, x) + + # cubic, expectation and variance + [st,st1,st2] = fun(alpha) + expect = st + 0.5*(taus**2)*st2 + vari = (taus**2)*(st1**2) + omegas**2 + + # Likelihood + lik = sum(ln_gauss_pdf_v(y,expect,np.sqrt(vari))) + + # Posterior + T = lik + prior_phi1 + prior_phi2 + + return T + +"""###mcmcmh: Metrolopolis-Hasting MCMC algorithm generating N chains of length M for a parameter vector A of length n.""" + +def mcmcmh(M, N, M0, Q, A0, tar, jump): + ''' + mcmcmh: Metrolopolis-Hasting MCMC algorithm generating N chains of length + M for a parameter vector A of length n. + + For details about the algorithm please refer to: + Gelman A, Carlin JB, Stern HS, Dunson DB, Vehtari A, Rubin DB. + Bayesian data analysis. CRC press; 2013 Nov 1. + ------------------------------------------------------------------------- + KJ, LRW, PMH + Version 2020-04-22 + ------------------------------------------------------------------------- + Inputs: + M: Length of the chains. + + N: Number of chains. + + M0: Burn in period. + + Q(nQ,1): Percentiles 0 <= Q(k) <= 100. + + A0(n,N): Array of feasible starting points: the target distribution + evaluated at A0(:,j) is strictly positive. + + Outputs: + + S(2+nQ,n): Summary of A - mean, standard deviation and percentile limits, + where the percentile limits are given by Q. + + aP(N,1): Acceptance percentages for AA calculated for each chain. + + Rh(n,1): Estimate of convergence. Theoretically Rh >= 1, and the closer + to 1, the more evidence of convergence. + + Ne(n,1): Estimate of the number of effective number of independent draws. + + AA(M,N,n): Array storing the chains: A(i,j,k) is the kth element of the + parameter vector stored as the ith member of the jth chain. + AA(1,j,:) = A0(:,j). + + IAA(M,N): Acceptance indices. IAA(i,j) = 1 means that the proposal + as(n,1) generated at the ith step of the jth chain was accepted so + that AA(i,j,:) = as. IAA(i,j) = 0 means that the proposal as(n,1) + generated at the ith step of the jth chain was rejected so that + AA(i,j,:) = AA(i-1,j,:), i > 1. The first set of proposal coincide with + A0 are all accepted, so IAA(1,j) = 1. + ''' + A0 = np.array(A0) + Q = np.array(Q) + # number of parameters for which samples are to be drawn + n = A0.shape[0] + + # number of percentiles to be evaluated + nQ = Q.size + + + # Initialising output arrays + AA = np.empty((M, N, n)) + IAA = np.zeros((M, N)) + Rh = np.empty((n)) + Ne = np.empty((n)) + S = np.empty((2 + nQ, n)) + + + # starting values of the sample and associated log of target density + aq = A0 + lq = tar(aq) + + # Starting values must be feasible for each chain + Id = lq > -np.Inf + + if sum(Id) < N: + print("Initial values must be feasible for all chains") + return None + + + # Run the chains in parallel + for q in range(M): + # draw from the jumping distribution and calculate + # d = log P0(aq|as) - log P0(as|aq) + + asam, d = jump(aq) + + # log of the target density for the new draw as + ls = tar(asam) + + # Metropolis-Hastings acceptance ratio + rq = np.exp(ls - lq + d) + + # draws from the uniform distribution for acceptance rule + uq = np.random.rand(N) + # index of samples that have been accepted + ind = uq < rq + + # updating the sample and evaluating acceptance indices + aq[:, ind] = asam[:, ind] + lq[ind] = ls[ind] + IAA[q, ind] = 1 + + # Store Metropolis Hastings sample + AA[q, :, :] = np.transpose(aq) + + + # acceptance probabilities for each chain + aP = 100*np.sum(IAA,0)/M + + # Convergence and summary statistics for each of the n parameters + for j in range(n): + # test convergence + RN = mcmcci(np.squeeze(AA[:,:,j]), M0) + Rh[j] = RN[0] + Ne[j] = RN[1] + + # provide summary information + asq = np.squeeze(AA[:,:,j]) + SS = mcsums(asq,M0,Q) + S[0,j] = SS[0] + S[1,j] = SS[1] + S[2:,j] = SS[2] + + return S, aP, Rh, Ne, AA, IAA + +def mcmcm_main(datay, datax, m0w, s0w, m0t, s0t, Mc, M0, Nc, Q): + + at0 = np.array((1,1,1,1, np.log(1/s0w**2), np.log(1/s0t**2))) + # function that evaluates the log of the target distribution at given parameter values + tar = lambda at: tar_at(at, datay, datax, m0w, s0w, m0t, s0t) + # function that evaluates the negative log of the target distribution to evaluate MAP estimates + mapp = lambda at: -tar(at) + + res = minimize(mapp, at0) + pars = res.x + V = res.hess_inv + L = np.linalg.cholesky(V) + + + # Function that draws sample from a Gaussian random walk jumping distribution + jump = lambda A: jumprwg(A, L) + + + rr = np.random.normal(0,1,size=(6,Nc)) + + A0 = np.matlib.repmat(pars.T,Nc,1).T + np.matmul(L,rr) + + + sam = mcmcmh(Mc,Nc,M0,Q,A0,tar,jump) + return 1/np.sqrt(np.exp(sam[0][0,-2:])) + +# Decaying exponential +def DecayExpFunction(a, b, f, x): + ''' + decaying exponential function evaluated at x with parameters a, b and f + ''' + return a * np.exp(-b * x) * np.sin(2 * np.pi * f * x) + + +# Decaying exponential - first derivative +def DecayExpFunction1der(a, b, f, x): + ''' + first derviative of the decaying exponential function evaluated at x with parameters a, b and f + ''' + return a * np.exp(-b * x) * ((2 * np.pi * f) * np.cos(2 * np.pi * f * x) - b * np.sin(2 * np.pi * f * x)) + + +# Decaying exponential - second derivative +def DecayExpFunction2der(a, b, f, x): + ''' + second derviative of the decaying exponential function evaluated at x with parameters a, b and f + ''' + return a * np.exp(-b * x) * ((-np.power((2 * np.pi * f), 2)) * np.sin(2 * np.pi * f * x) + - ((4 * b * np.pi * f) * np.cos(2 * np.pi * f * x)) + b * np.sin(2 * np.pi * f * x)) + + +class MCMCMH_NJ(): + ''' + Bayesian Noise and jitter reduction algorithm. MCMC used to determine the noise and jitter variances. + Noise and jitter variances are then used in an iterative algorithm to remove the noise and jitter from the signal + ''' + + def __init__(self, fs, ydata, N, niter, tol, m0w, s0w, m0t, s0t, Mc, M0, Nc, Q): + 'Setting initial variables' + + # variables for AnalyseSignalN and NJAlgorithm + self.fs = fs + self.ydata = ydata + # self.xdata = xdata + self.N = N + self.niter = niter + self.tol = tol + self.m0w = m0w + self.s0w = s0w + self.m0t = m0t + self.s0t = s0t + self.Mc = Mc + self.M0 = M0 + self.Nc = Nc + self.Q = Q + # outs = MCMCMH.mcmcm_decayexp.main(ydata, xdata, m0w, s0w, m0t, s0t, Mc, M0, Nc, Q) + # self.jitterSD = outs[0] + # self.noiseSD = outs[1] + + def AnalyseSignalN(self): + ''' + Analyse signal to remove noise and jitter providing signal estimates with associated + uncertainty. Uses normalised independent variable + ''' + + # Initialisation + self.N = np.int64(self.N) # Converting window length integer to int64 format + m = np.size(self.ydata) # Signal data length + m = np.int64(m) # Converting signal length to int64 format + + # Setting initial variables + n = (self.N - 1) // 2 + # Covariance matric for signal values + Vmat = np.zeros((np.multiply(2, self.N) - 1, np.multiply(2, self.N) - 1)) + # Sensitivtity vecotrs + cvec = np.zeros(n, self.N) + # Sensitivtity matrix + Cmat = np.zeros((self.N, np.multiply(2, self.N) - 1)) + # Initial signal estimate + yhat0 = np.full((m, 1), np.nan) + # Final signal estimate + yhat = np.full((m, 1), np.nan) + # Uncertainty infromation for estimates + vyhat = np.full((m, self.N), np.nan) + # Consistency matrix + R = np.full((m, 1), np.nan) + + # Values of normalised independent variables + datax = np.divide(np.arange(-n, n + 1), self.fs) + + outs = mcmcm_main(self.ydata, datax, self.m0w, self.s0w, self.m0t, self.s0t, self.Mc, + self.M0, self.Nc, self.Q) + self.jitterSD = outs[0] + self.noiseSD = outs[1] + # Loop through indices L of window + L = 0 + if np.size(self.ydata) > 10: + # for L in range(1, m-self.N+1): + # while + # Index k of indices L of windows + k = L + n + # print(k) + # Extract data in window + datay = self.ydata[L:L + self.N] + # Inital polynomial approximation + p = np.polyfit(datax, datay, 3) + pval = np.polyval(p, datax) + yhat0[k] = pval[n] + + # Applying algortithm to remove noise and jitter + [yhat[k], ck, vark, R[k]] = MCMCMH_NJ.NJAlgorithm(self, datax, datay, p, pval) + print(yhat[k]) + # First n windows, start building the covariance matrix Vmat for the data + if L < n + 1: + Vmat[L - 1, L - 1] = vark + + # for windows n+1 to 2n, continue building the covariance matrix Vmat and start stroing the + # sensitivtity vectors ck in cvec + elif L > n and L < np.multiply(2, n) + 1: + Vmat[L - 1, L - 1] = vark + cvec[L - n - 1, :] = ck + + + # For windows between 2n+1 and 4n, continue to build Vmat and cvec, and start building the sensitivtity + # matrix Cmat from the sensitivtity vecotrs. Also, evaluate uncertainties for pervious estimates. + elif L > np.multiply(2, n) and L < np.multiply(4, n) + 2: + + Vmat[L - 1, L - 1] = vark + # Count for building sensitivtity matrix + iC = L - np.multiply(2, n) + # Start building sensitivtity matrix from cvec + Cmat[iC - 1, :] = np.concatenate((np.zeros((1, iC - 1)), cvec[0, :], np.zeros((1, self.N - iC))), + axis=None) + + # Removing the first row of cvec and shift every row up one - creating + # empty last row + cvec = np.roll(cvec, -1, axis=0) + cvec[-1, :] = 0 + # Update empty last row + cvec[n - 1, :] = ck + Cmatk = Cmat[0:iC, 1:self.N - 1 + iC] + + # Continue building Vmat + Vmatk = Vmat[1:self.N - 1 + iC, 1:self.N - 1 + iC] + V = np.matmul(np.matmul(Cmatk, Vmatk), np.transpose(Cmatk)) + vhempty = np.empty((1, self.N - iC)) + vhempty[:] = np.nan + # Begin building vyhat + vyhat[L, :] = np.concatenate((vhempty, V[iC - 1, :]), axis=None) + + + # For the remaining windows, update Vmat, Cmat and cvec. Continue to + # evaluate the uncertainties for previous estimates. + elif L > np.multiply(4, n) + 1: + # Update Vmat + Vmat = np.delete(Vmat, 0, axis=0) + Vmat = np.delete(Vmat, 0, axis=1) + Vmatzeros = np.zeros([Vmat.shape[0] + 1, Vmat.shape[1] + 1]) + Vmatzeros[:Vmat.shape[0], :Vmat.shape[1]] = Vmat + Vmat = Vmatzeros + Vmat[2 * self.N - 2, 2 * self.N - 2] = vark + + # Building updated Cmat matrix + Cmat_old = np.concatenate((Cmat[1:self.N, 1:2 * self.N - 1], np.zeros([self.N - 1, 1])), axis=1) + Cmat_new = np.concatenate((np.zeros([1, self.N - 1]), cvec[0, :]), axis=None) + Cmat = np.concatenate((Cmat_old, Cmat_new[:, None].T), axis=0) + + # Update cvec + cvec = np.roll(cvec, -1, axis=0) + cvec[-1, :] = 0 + + # Continue building vyhat + V = np.matmul(np.matmul(Cmat, Vmat), np.transpose(Cmat)) + vyhat[L, :] = V[self.N - 2, :] + + L += 1 + + return (yhat[k]) + + def NJAlgorithm(self, datax, datay, p0, p0x): + ''' + Noise and Jitter Removal Algorithm- Iterative scheme that preprocesses data to reduce the effects of + noise and jitter, resulting in an estimate of the true signal along with its associated uncertainty. + + Refer paper for details: https://ieeexplore.ieee.org/document/9138266 + ''' + + # Initialisatio + iter_ = 0 + delta = np.multiply(2, self.tol) + + # Values of basis function at central point in window + N = np.size(datax) + n = np.divide(self.N - 2, 2) + k = np.int64(n + 1) + t = np.array([np.power(datax[k], 3), np.power(datax[k], 2), datax[k], 1]) + # Deisgn Matrix + X = np.array([np.power(datax, 3) + 3 * np.multiply(np.power(self.jitterSD, 2), datax), + np.power(datax, 2) + np.power(self.jitterSD, 2), datax, np.ones(np.size(datax))]) + X = X.T + + # Iterative algortithm + while delta >= self.tol: + # Increment number of iterations + iter_ = iter_ + 1 + + # Step 2 - Polynomial fitting over window + pd = np.polyder(p0) + pdx = np.polyval(pd, datax) + + # Step 4 + # Weight calculation + w = np.divide(1, [np.sqrt(np.power(self.jitterSD, 2) * np.power(pdx, 2) + + np.power(self.noiseSD, 2))]) + w = w.T + + # Calculating polynomial coeffs + Xt = np.matmul(np.diagflat(w), X) + C = np.matmul(np.linalg.pinv(Xt), np.diagflat(w)) + datay = datay.T + p1 = np.matmul(C, datay) + p1x = np.polyval(p1, datax) + + # Step 5 - stablise process + delta = np.max(np.abs(p1x - p0x)) + p0 = p1 + p0x = p1x + + if iter_ == self.niter: + print('Maximum number of iterations reached') + break + + # Evaluate outputs + c = np.matmul(t, C) + yhat = np.matmul(c, datay) + pd = np.polyder(p0) + pdx = np.polyval(pd, datax[k]) + vk = np.power(self.jitterSD, 2) * np.power(pdx, 2) + np.power(self.noiseSD, 2) + R = np.power(np.linalg.norm(np.matmul(np.diagflat(w), (datay - np.matmul(X, p0)))), 2) + + return yhat, c, vk, R + + +def random_gaussian_whrand(M, mu, sigma, istate1, istate2): + ''' + Generates random numbers from a Gaussian Distribution using the Wichmann–Hill random number generator + + Inputs: + M: Number of random numbers required + mu: Mean of the Gaussian + sigma: Std of the Gaussian + istate1: vector of 4 integers + istate2: vector of 4 different integers + + Output: + xdist: Random Gaussian vector of size M + ''' + mn = 0 + ndist = np.zeros(M) + while mn < M: + rr1 = whrand(istate1, 1) + rr2 = whrand(istate2, 1) + + v1, istate1 = rr1 + v2, istate2 = rr2 + + if mn < M: + ndist[mn] = math.sqrt(-2 * math.log(v1)) * math.cos(2 * math.pi * v2) + mn = mn + 1 + + if mn < M: + ndist[mn] = math.sqrt(-2 * math.log(v1)) * math.sin(2 * math.pi * v2) + mn = mn + 1 + + xdist = mu + sigma * ndist + + return xdist, istate1, istate2 + + +def whrand(istate, N): + ''' + Generates uniform random numbers between 0 and 1 using the Wichmann–Hill random number generator + + Inputs: + istate: vector of 4 integers + N: Number of random numbers required + + Outputs: + r: Vector uniform random numbers of size M + istate: Output vector of 4 integers + ''' + + # Constants + a = np.array([11600, 47003, 23000, 33000]) + b = np.array([185127, 45688, 93368, 65075]) + c = np.array([10379, 10479, 19423, 8123]) + d = np.array([456, 420, 300, 0]) + 2147483123 + + r = np.zeros(N) + for i in range(N): + # Update states + for j in range(4): + istate[j] = a[j] * np.mod(istate[j], b[j]) - c[j] * np.fix(istate[j] / b[j]) + if istate[j] < 0: + istate[j] = istate[j] + d[j] + + # Evaluate random number + w = np.sum(np.divide(istate, d)) + r[i] = np.remainder(w, 1) + + return r, istate + + +def njr(fs, ydata, N, niter, tol, m0w, s0w, m0t, s0t, Mc, M0, Nc, Q): + analyse_fun = MCMCMH_NJ(fs, ydata, N, niter, tol, m0w, s0w, m0t, s0t, Mc, M0, Nc, Q) + yhat1= analyse_fun.AnalyseSignalN() + return yhat1 + +class NJRemoved(AgentMET4FOF): + def init_parameters(self, fs=100, ydata = np.array([]), N=15, niter=100, tol=1e-9, m0w = 10, s0w = 0.0005, m0t = 10, s0t = 0.0002*100/8, Mc=5000, M0=100, Nc=100, Q=50 ): + self.fs = fs + self.ydata = ydata + self.N = N + self.niter = niter + self.tol = tol + self.m0w = m0w + self.s0w = s0w + self.m0t = m0t + self.s0t = s0t + self.Mc = Mc + self.M0 = M0 + self.Nc = Nc + self.Q = Q + + + def on_received_message(self, message): + ddata = message['data'] + self.ydata = np.append(self.ydata, ddata) + if np.size(self.ydata) == self.N: + t = njr(self.fs, self.ydata, self.N, self.niter, self.tol, self.m0w, self.s0w, self.m0t, self.s0t, self.Mc, self.M0,self.Nc, self.Q) + self.send_output(self.ydata[7] - t) + self.ydata = self.ydata[1:self.N] + +def njr(fs, ydata, N, niter, tol, m0w, s0w, m0t, s0t, Mc, M0, Nc, Q): + analyse_fun = MCMCMH_NJ(fs, ydata, N, niter, tol, m0w, s0w, m0t, s0t, Mc, M0, Nc, Q) + yhat1= analyse_fun.AnalyseSignalN() + return yhat1 + + +######################################## +class NJRemoved(AgentMET4FOF): + def init_parameters(self, fs=100, ydata = np.array([]), N=15, niter=100, tol=1e-9, m0w = 10, s0w = 0.0005, m0t = 10, s0t = 0.0002*100/8, Mc=5000, M0=100, Nc=100, Q=50 ): + self.fs = fs + self.ydata = ydata + self.N = N + self.niter = niter + self.tol = tol + self.m0w = m0w + self.s0w = s0w + self.m0t = m0t + self.s0t = s0t + self.Mc = Mc + self.M0 = M0 + self.Nc = Nc + self.Q = Q + + + def on_received_message(self, message): + ddata = message['data'] + self.ydata = np.append(self.ydata, ddata) + if np.size(self.ydata) == self.N: + t = njr(self.fs, self.ydata, self.N, self.niter, self.tol, self.m0w, self.s0w, self.m0t, self.s0t, self.Mc, self.M0,self.Nc, self.Q) + self.send_output(self.ydata[7] - t) + self.ydata = self.ydata[1:self.N] + + +class SineGeneratorAgent(AgentMET4FOF): + def init_parameters(self): + self.stream = StaticSineGeneratorWithJitter(jittersd=0.0005, noisesd=0.0002) + + def agent_loop(self): + if self.current_state == "Running": + sine_data = self.stream.next_sample() # dictionary + self.send_output(sine_data['quantities']) + +def main(): + # start agent network server + agentNetwork = AgentNetwork(backend="mesa") + # init agents + + gen_agent = agentNetwork.add_agent(agentType=SineGeneratorAgent) + + njremove_agent = agentNetwork.add_agent(agentType=NJRemoved) + + monitor_agent = agentNetwork.add_agent(agentType=MonitorAgent) + monitor_agent2 = agentNetwork.add_agent(agentType=MonitorAgent) + + + # connect agents : We can connect multiple agents to any particular agent + + agentNetwork.bind_agents(gen_agent, njremove_agent) + + # connect + agentNetwork.bind_agents(gen_agent, monitor_agent) + + # agentNetwork.bind_agents(njremove_agent, moitor_agent2) + + agentNetwork.bind_agents(njremove_agent, monitor_agent2) + + # set all agents states to "Running" + agentNetwork.set_running_state() + + # allow for shutting down the network after execution + return agentNetwork + + +if __name__ == '__main__': + main() \ No newline at end of file From 5bed0c6d41338e6030179f3edd18d6d768c77d62 Mon Sep 17 00:00:00 2001 From: Anupam Prasad Vedurmudi Date: Tue, 27 Jul 2021 14:33:05 +0200 Subject: [PATCH 05/43] wip(noise_jitter_removal_agent): re-added matlib import --- agentMET4FOF/agents/noise_jtter_removal_agents.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/agentMET4FOF/agents/noise_jtter_removal_agents.py b/agentMET4FOF/agents/noise_jtter_removal_agents.py index 67cbc9d0..97dc5e78 100644 --- a/agentMET4FOF/agents/noise_jtter_removal_agents.py +++ b/agentMET4FOF/agents/noise_jtter_removal_agents.py @@ -1,6 +1,7 @@ import math import numpy as np +import numpy.matlib from scipy.optimize import minimize from agentMET4FOF.agents import AgentMET4FOF, AgentNetwork, MonitorAgent @@ -826,7 +827,7 @@ def on_received_message(self, message): class SineGeneratorAgent(AgentMET4FOF): def init_parameters(self): - self.stream = StaticSineGeneratorWithJitter(jittersd=0.0005, noisesd=0.0002) + self.stream = StaticSineGeneratorWithJitter() def agent_loop(self): if self.current_state == "Running": From 81ca2466f626e5b5bdfa0cdcfac77d68c2216ba0 Mon Sep 17 00:00:00 2001 From: Anupam Prasad Vedurmudi Date: Tue, 27 Jul 2021 12:20:06 +0200 Subject: [PATCH 06/43] wip(noise_jitter_removal_agent): methods and classes copied over from npl repo --- .../agents/noise_jtter_removal_agents.py | 868 ++++++++++++++++++ 1 file changed, 868 insertions(+) create mode 100644 agentMET4FOF/agents/noise_jtter_removal_agents.py diff --git a/agentMET4FOF/agents/noise_jtter_removal_agents.py b/agentMET4FOF/agents/noise_jtter_removal_agents.py new file mode 100644 index 00000000..67cbc9d0 --- /dev/null +++ b/agentMET4FOF/agents/noise_jtter_removal_agents.py @@ -0,0 +1,868 @@ +import math + +import numpy as np +from scipy.optimize import minimize + +from agentMET4FOF.agents import AgentMET4FOF, AgentNetwork, MonitorAgent +from agentMET4FOF.streams.signal_streams import StaticSineGeneratorWithJitter + +"""### mcmcci: MCMC convergence indices for multiple chains.""" + +def mcmcci(A,M0): + ''' + mcmcci: MCMC convergence indices for multiple chains. + ---------------------------------------------------------------------------- + KJ, LRW, PMH + + Version 2020-04-22 + + ---------------------------------------------------------------------------- + Inputs: + A(M, N): Chain samples, N chains of length M. + + M0: Length of burn-in, M > M0 >= 0. + + Outputs: + Rhat: Convergence index. Rhat is expected to be greater than 1. The closer + Rhat is to 1, the better the convergence. + + Neff: Estimate of the effective number of independent draws. Neff is + expected to be less than (M-M0)*N. + ---------------------------------------------------------------------------- + Note: If the calculated value of Rhat is < 1, then Rhat is set to 1 and Neff + set to (M-M0)*N, their limit values. + + Note: If N = 1 or M0 > M-2, Rhat = 0; Neff = 0. + ''' + A = np.array(A) + + M, N = A.shape + + # Initialisation + Rhat = 0 + Neff = 0 + + # The convergence statistics can only be evaluated if there are multiple chains + # and the chain length is greater than the burn in length + + if N > 1 and M > M0 + 1: + Mt = M - M0 + + # Chain mean and mean of means + asub = A[M0:,:] + ad = np.mean(asub,axis = 0) + add = asub.mean() + + # Within group standard deviation + ss = np.std(asub,axis = 0) + + # Between groups variance. + dd = np.square(ad - add) + B = (Mt*np.sum(dd))/(N-1) + + # Within groups variance. + W = np.sum(np.square(ss))/N + + # V plus + Vp = (Mt-1)*W/Mt + B/Mt + + # Convergence statistic, effective number of independent samples + Rhat = np.sqrt(Vp/W) + Neff = Mt*N*Vp/B + + Rhat = np.maximum(Rhat,1) + Neff = np.minimum(Neff,Mt*N) + + return Rhat, Neff + +"""### mcsums: Summary information from MC samples.""" + +def mcsums(A,M0,Q): + ''' + mcsums: Summary information from MC samples. + ------------------------------------------------------------------------- + KJ, LRW, PMH + Version 2020-04-22 + + ------------------------------------------------------------------------- + Inputs: + A(M,N): An array that stores samples of size M x N. + + M0: Burn-in period with M > M0 >= 0. + + Q(nQ,1): Percentiles specifications, 0 <= Q(l) <= 100. + + Outputs: + abar(n,1): Mean for each sample. + + s(n,1): Standard deviation for sample. + + aQ(nQ,n): Percentiles corresponding to Q. + + ''' + # Size of samples after burn-in + A = np.array(A) + + M, N = A.shape + + m = (M - M0)*N + + # Initialise percentile vector + nQ = Q.size + aQ = np.zeros(nQ) + + # Samples from N parallel chains after burn-in period + aaj = A[M0:, :] + aaj = aaj.flatten() + + # Mean and standard deviation of samples + abar = np.mean(aaj) + s = np.std(aaj) + + # Percentiles of samples + aQ = np.percentile(aaj,Q) + + return abar, s, aQ + +"""### jumprwg: Jumping distribution for the Metropolis Hastings Gaussian random walk algorithm""" + +def jumprwg(A, L): + ''' + jumprwg: Jumping distribution for the Metropolis Hastings Gaussian random + walk algorithm + ------------------------------------------------------------------------- + KJ, LRW, PMH + Version 2020-04-22 + ------------------------------------------------------------------------- + Inputs: + A(n,N): Samples at the current iteration + + L(n,n): Cholesky factor of variance of parameter vector. + + Outputs: + As(n,N): Proposed parameter array which is randomly sampled from the + jumping distribution + + dp0: The difference between the logarithm of the jumping distribution + associated with moving from A(:,j) to As(:,j) and that associated with + moving from As(:,j) to A(:,j), up to an additive constant. + log P0(a|as) - log P0(as|a) + + ''' + # Number of parameters and parallel chains + n, N = A.shape + + # random draw from a Gaussian distribution + e = np.random.normal(0, 1, size=(n,N)) + + # proposed draw from a Gaussian distribution with mean A and variance LL' + As = A + np.matmul(L,e) + + # For a Gaussian random walk, since log P0(a|as) = log P0(as|a), dp0 will always be zero + dp0 = np.zeros(N) + + return As, dp0 + +"""### Cubic function and its first and second derivative""" + +def fgh_cubic(alpha,t): + ''' + ------------------------------------------------------------------------- + Cubic function and its first and second derivative + ------------------------------------------------------------------------- + KJ, LRW, PMH + Version 2020-04-22 + ------------------------------------------------------------------------- + Inputs: + alpha(4,N): Alpha parameters + + t(m,1): Times + + Outputs: + f(m,N): Cubic function + + f1(m,N): Derivative of cubic + + f2(m,N): Second derivative of cubic + ''' + + # length of data and number of paramaters + m = t.size + + # design matrix + C = np.array([np.ones(m), t, t**2, t**3]) + + # derivate info + C1 = np.array([np.ones(m), 2*t, 3*t**2]) + C2 = np.array([2*np.ones(m), 6*t]) + + # cubic and derivatives + f = np.matmul(C.T,alpha) + f1 = np.matmul(C1.T,alpha[1:]) + f2 = np.matmul(C2.T,alpha[2:]) + + return f, f1, f2 + +"""### Log of the gaussian pdf""" + +def ln_gauss_pdf_v(x,mu,sigma): + ''' + ------------------------------------------------------------------------- + Log of the Gaussian pdf + ------------------------------------------------------------------------- + KJ, LRW, PMH + Version 2020-03-12 + -------------------------------------------------------------------------- + Inputs: + x(m,1): Points at which pdf is to be evaluated + + mu: Mean of distribution + + sigma: Standard deviation of the distribution + + Output: + logf: Log of the Gaussian pdf at x with mean mu and + std sigma + ''' + + + try: + # When inputs are high dimensional arrays/matrices + xx = np.matlib.repmat(x,mu.shape[1],1) + xx = xx.T + logk = - np.log(2*math.pi)/2 - np.log(sigma) + logp = -((xx - mu)**2)/(2*sigma**2) + # Log of the Gaussian PDF + logf = logk + logp + + except IndexError: + # When inputs are vectors + logk = - np.log(2*math.pi)/2 - np.log(sigma) + logp = -((x - mu)**2)/(2*sigma**2) + # Log of the Gaussian PDF + logf = logk + logp + + return logf + +"""### Target dist for noise and jitter posterior dist""" + +def tar_at(at, y, x, m0w, s0w, m0t, s0t): + ''' + ------------------------------------------------------------------------- + Target dist for noise and jitter posterior dist + ------------------------------------------------------------------------- + KJ, LRW, PMH + Version 2020-04-22 + -------------------------------------------------------------------------- + Inputs: + at(n+2,N): Parameters alpha, log(1/tau^2) and log(1/w^2) + + y(m,1): Signal + + x(m,1): time at which signal was recorded + + s0w and s0t: prior estimates of tau and omega + + m0w and m0t: degree of belief in prior estimates for tau and omega + + Output: + T: Log of the posterior distribution + ''' + + # Size of parameter vector + at = np.array(at) + p = at.shape[0] + + # Number of alphas + n = p - 2 + + # Extract parameters + alpha = at[0:n] + phi1 = np.exp(at[-2]) + phi2 = np.exp(at[-1]) + taus = np.ones(phi1.shape)/np.sqrt(phi1) + omegas = np.ones(phi2.shape)/np.sqrt(phi2) + + # Gamma priors for phis + prior_phi1 = (m0t/2)*np.log(phi1) - phi1*m0t*s0t**2/2 + prior_phi2 = (m0w/2)*np.log(phi2) - phi2*m0w*s0w**2/2 + + + # function that evaluates the cubic function with user specified cubic parameters + fun = lambda aa: fgh_cubic(aa, x) + + # cubic, expectation and variance + [st,st1,st2] = fun(alpha) + expect = st + 0.5*(taus**2)*st2 + vari = (taus**2)*(st1**2) + omegas**2 + + # Likelihood + lik = sum(ln_gauss_pdf_v(y,expect,np.sqrt(vari))) + + # Posterior + T = lik + prior_phi1 + prior_phi2 + + return T + +"""###mcmcmh: Metrolopolis-Hasting MCMC algorithm generating N chains of length M for a parameter vector A of length n.""" + +def mcmcmh(M, N, M0, Q, A0, tar, jump): + ''' + mcmcmh: Metrolopolis-Hasting MCMC algorithm generating N chains of length + M for a parameter vector A of length n. + + For details about the algorithm please refer to: + Gelman A, Carlin JB, Stern HS, Dunson DB, Vehtari A, Rubin DB. + Bayesian data analysis. CRC press; 2013 Nov 1. + ------------------------------------------------------------------------- + KJ, LRW, PMH + Version 2020-04-22 + ------------------------------------------------------------------------- + Inputs: + M: Length of the chains. + + N: Number of chains. + + M0: Burn in period. + + Q(nQ,1): Percentiles 0 <= Q(k) <= 100. + + A0(n,N): Array of feasible starting points: the target distribution + evaluated at A0(:,j) is strictly positive. + + Outputs: + + S(2+nQ,n): Summary of A - mean, standard deviation and percentile limits, + where the percentile limits are given by Q. + + aP(N,1): Acceptance percentages for AA calculated for each chain. + + Rh(n,1): Estimate of convergence. Theoretically Rh >= 1, and the closer + to 1, the more evidence of convergence. + + Ne(n,1): Estimate of the number of effective number of independent draws. + + AA(M,N,n): Array storing the chains: A(i,j,k) is the kth element of the + parameter vector stored as the ith member of the jth chain. + AA(1,j,:) = A0(:,j). + + IAA(M,N): Acceptance indices. IAA(i,j) = 1 means that the proposal + as(n,1) generated at the ith step of the jth chain was accepted so + that AA(i,j,:) = as. IAA(i,j) = 0 means that the proposal as(n,1) + generated at the ith step of the jth chain was rejected so that + AA(i,j,:) = AA(i-1,j,:), i > 1. The first set of proposal coincide with + A0 are all accepted, so IAA(1,j) = 1. + ''' + A0 = np.array(A0) + Q = np.array(Q) + # number of parameters for which samples are to be drawn + n = A0.shape[0] + + # number of percentiles to be evaluated + nQ = Q.size + + + # Initialising output arrays + AA = np.empty((M, N, n)) + IAA = np.zeros((M, N)) + Rh = np.empty((n)) + Ne = np.empty((n)) + S = np.empty((2 + nQ, n)) + + + # starting values of the sample and associated log of target density + aq = A0 + lq = tar(aq) + + # Starting values must be feasible for each chain + Id = lq > -np.Inf + + if sum(Id) < N: + print("Initial values must be feasible for all chains") + return None + + + # Run the chains in parallel + for q in range(M): + # draw from the jumping distribution and calculate + # d = log P0(aq|as) - log P0(as|aq) + + asam, d = jump(aq) + + # log of the target density for the new draw as + ls = tar(asam) + + # Metropolis-Hastings acceptance ratio + rq = np.exp(ls - lq + d) + + # draws from the uniform distribution for acceptance rule + uq = np.random.rand(N) + # index of samples that have been accepted + ind = uq < rq + + # updating the sample and evaluating acceptance indices + aq[:, ind] = asam[:, ind] + lq[ind] = ls[ind] + IAA[q, ind] = 1 + + # Store Metropolis Hastings sample + AA[q, :, :] = np.transpose(aq) + + + # acceptance probabilities for each chain + aP = 100*np.sum(IAA,0)/M + + # Convergence and summary statistics for each of the n parameters + for j in range(n): + # test convergence + RN = mcmcci(np.squeeze(AA[:,:,j]), M0) + Rh[j] = RN[0] + Ne[j] = RN[1] + + # provide summary information + asq = np.squeeze(AA[:,:,j]) + SS = mcsums(asq,M0,Q) + S[0,j] = SS[0] + S[1,j] = SS[1] + S[2:,j] = SS[2] + + return S, aP, Rh, Ne, AA, IAA + +def mcmcm_main(datay, datax, m0w, s0w, m0t, s0t, Mc, M0, Nc, Q): + + at0 = np.array((1,1,1,1, np.log(1/s0w**2), np.log(1/s0t**2))) + # function that evaluates the log of the target distribution at given parameter values + tar = lambda at: tar_at(at, datay, datax, m0w, s0w, m0t, s0t) + # function that evaluates the negative log of the target distribution to evaluate MAP estimates + mapp = lambda at: -tar(at) + + res = minimize(mapp, at0) + pars = res.x + V = res.hess_inv + L = np.linalg.cholesky(V) + + + # Function that draws sample from a Gaussian random walk jumping distribution + jump = lambda A: jumprwg(A, L) + + + rr = np.random.normal(0,1,size=(6,Nc)) + + A0 = np.matlib.repmat(pars.T,Nc,1).T + np.matmul(L,rr) + + + sam = mcmcmh(Mc,Nc,M0,Q,A0,tar,jump) + return 1/np.sqrt(np.exp(sam[0][0,-2:])) + +# Decaying exponential +def DecayExpFunction(a, b, f, x): + ''' + decaying exponential function evaluated at x with parameters a, b and f + ''' + return a * np.exp(-b * x) * np.sin(2 * np.pi * f * x) + + +# Decaying exponential - first derivative +def DecayExpFunction1der(a, b, f, x): + ''' + first derviative of the decaying exponential function evaluated at x with parameters a, b and f + ''' + return a * np.exp(-b * x) * ((2 * np.pi * f) * np.cos(2 * np.pi * f * x) - b * np.sin(2 * np.pi * f * x)) + + +# Decaying exponential - second derivative +def DecayExpFunction2der(a, b, f, x): + ''' + second derviative of the decaying exponential function evaluated at x with parameters a, b and f + ''' + return a * np.exp(-b * x) * ((-np.power((2 * np.pi * f), 2)) * np.sin(2 * np.pi * f * x) + - ((4 * b * np.pi * f) * np.cos(2 * np.pi * f * x)) + b * np.sin(2 * np.pi * f * x)) + + +class MCMCMH_NJ(): + ''' + Bayesian Noise and jitter reduction algorithm. MCMC used to determine the noise and jitter variances. + Noise and jitter variances are then used in an iterative algorithm to remove the noise and jitter from the signal + ''' + + def __init__(self, fs, ydata, N, niter, tol, m0w, s0w, m0t, s0t, Mc, M0, Nc, Q): + 'Setting initial variables' + + # variables for AnalyseSignalN and NJAlgorithm + self.fs = fs + self.ydata = ydata + # self.xdata = xdata + self.N = N + self.niter = niter + self.tol = tol + self.m0w = m0w + self.s0w = s0w + self.m0t = m0t + self.s0t = s0t + self.Mc = Mc + self.M0 = M0 + self.Nc = Nc + self.Q = Q + # outs = MCMCMH.mcmcm_decayexp.main(ydata, xdata, m0w, s0w, m0t, s0t, Mc, M0, Nc, Q) + # self.jitterSD = outs[0] + # self.noiseSD = outs[1] + + def AnalyseSignalN(self): + ''' + Analyse signal to remove noise and jitter providing signal estimates with associated + uncertainty. Uses normalised independent variable + ''' + + # Initialisation + self.N = np.int64(self.N) # Converting window length integer to int64 format + m = np.size(self.ydata) # Signal data length + m = np.int64(m) # Converting signal length to int64 format + + # Setting initial variables + n = (self.N - 1) // 2 + # Covariance matric for signal values + Vmat = np.zeros((np.multiply(2, self.N) - 1, np.multiply(2, self.N) - 1)) + # Sensitivtity vecotrs + cvec = np.zeros(n, self.N) + # Sensitivtity matrix + Cmat = np.zeros((self.N, np.multiply(2, self.N) - 1)) + # Initial signal estimate + yhat0 = np.full((m, 1), np.nan) + # Final signal estimate + yhat = np.full((m, 1), np.nan) + # Uncertainty infromation for estimates + vyhat = np.full((m, self.N), np.nan) + # Consistency matrix + R = np.full((m, 1), np.nan) + + # Values of normalised independent variables + datax = np.divide(np.arange(-n, n + 1), self.fs) + + outs = mcmcm_main(self.ydata, datax, self.m0w, self.s0w, self.m0t, self.s0t, self.Mc, + self.M0, self.Nc, self.Q) + self.jitterSD = outs[0] + self.noiseSD = outs[1] + # Loop through indices L of window + L = 0 + if np.size(self.ydata) > 10: + # for L in range(1, m-self.N+1): + # while + # Index k of indices L of windows + k = L + n + # print(k) + # Extract data in window + datay = self.ydata[L:L + self.N] + # Inital polynomial approximation + p = np.polyfit(datax, datay, 3) + pval = np.polyval(p, datax) + yhat0[k] = pval[n] + + # Applying algortithm to remove noise and jitter + [yhat[k], ck, vark, R[k]] = MCMCMH_NJ.NJAlgorithm(self, datax, datay, p, pval) + print(yhat[k]) + # First n windows, start building the covariance matrix Vmat for the data + if L < n + 1: + Vmat[L - 1, L - 1] = vark + + # for windows n+1 to 2n, continue building the covariance matrix Vmat and start stroing the + # sensitivtity vectors ck in cvec + elif L > n and L < np.multiply(2, n) + 1: + Vmat[L - 1, L - 1] = vark + cvec[L - n - 1, :] = ck + + + # For windows between 2n+1 and 4n, continue to build Vmat and cvec, and start building the sensitivtity + # matrix Cmat from the sensitivtity vecotrs. Also, evaluate uncertainties for pervious estimates. + elif L > np.multiply(2, n) and L < np.multiply(4, n) + 2: + + Vmat[L - 1, L - 1] = vark + # Count for building sensitivtity matrix + iC = L - np.multiply(2, n) + # Start building sensitivtity matrix from cvec + Cmat[iC - 1, :] = np.concatenate((np.zeros((1, iC - 1)), cvec[0, :], np.zeros((1, self.N - iC))), + axis=None) + + # Removing the first row of cvec and shift every row up one - creating + # empty last row + cvec = np.roll(cvec, -1, axis=0) + cvec[-1, :] = 0 + # Update empty last row + cvec[n - 1, :] = ck + Cmatk = Cmat[0:iC, 1:self.N - 1 + iC] + + # Continue building Vmat + Vmatk = Vmat[1:self.N - 1 + iC, 1:self.N - 1 + iC] + V = np.matmul(np.matmul(Cmatk, Vmatk), np.transpose(Cmatk)) + vhempty = np.empty((1, self.N - iC)) + vhempty[:] = np.nan + # Begin building vyhat + vyhat[L, :] = np.concatenate((vhempty, V[iC - 1, :]), axis=None) + + + # For the remaining windows, update Vmat, Cmat and cvec. Continue to + # evaluate the uncertainties for previous estimates. + elif L > np.multiply(4, n) + 1: + # Update Vmat + Vmat = np.delete(Vmat, 0, axis=0) + Vmat = np.delete(Vmat, 0, axis=1) + Vmatzeros = np.zeros([Vmat.shape[0] + 1, Vmat.shape[1] + 1]) + Vmatzeros[:Vmat.shape[0], :Vmat.shape[1]] = Vmat + Vmat = Vmatzeros + Vmat[2 * self.N - 2, 2 * self.N - 2] = vark + + # Building updated Cmat matrix + Cmat_old = np.concatenate((Cmat[1:self.N, 1:2 * self.N - 1], np.zeros([self.N - 1, 1])), axis=1) + Cmat_new = np.concatenate((np.zeros([1, self.N - 1]), cvec[0, :]), axis=None) + Cmat = np.concatenate((Cmat_old, Cmat_new[:, None].T), axis=0) + + # Update cvec + cvec = np.roll(cvec, -1, axis=0) + cvec[-1, :] = 0 + + # Continue building vyhat + V = np.matmul(np.matmul(Cmat, Vmat), np.transpose(Cmat)) + vyhat[L, :] = V[self.N - 2, :] + + L += 1 + + return (yhat[k]) + + def NJAlgorithm(self, datax, datay, p0, p0x): + ''' + Noise and Jitter Removal Algorithm- Iterative scheme that preprocesses data to reduce the effects of + noise and jitter, resulting in an estimate of the true signal along with its associated uncertainty. + + Refer paper for details: https://ieeexplore.ieee.org/document/9138266 + ''' + + # Initialisatio + iter_ = 0 + delta = np.multiply(2, self.tol) + + # Values of basis function at central point in window + N = np.size(datax) + n = np.divide(self.N - 2, 2) + k = np.int64(n + 1) + t = np.array([np.power(datax[k], 3), np.power(datax[k], 2), datax[k], 1]) + # Deisgn Matrix + X = np.array([np.power(datax, 3) + 3 * np.multiply(np.power(self.jitterSD, 2), datax), + np.power(datax, 2) + np.power(self.jitterSD, 2), datax, np.ones(np.size(datax))]) + X = X.T + + # Iterative algortithm + while delta >= self.tol: + # Increment number of iterations + iter_ = iter_ + 1 + + # Step 2 - Polynomial fitting over window + pd = np.polyder(p0) + pdx = np.polyval(pd, datax) + + # Step 4 + # Weight calculation + w = np.divide(1, [np.sqrt(np.power(self.jitterSD, 2) * np.power(pdx, 2) + + np.power(self.noiseSD, 2))]) + w = w.T + + # Calculating polynomial coeffs + Xt = np.matmul(np.diagflat(w), X) + C = np.matmul(np.linalg.pinv(Xt), np.diagflat(w)) + datay = datay.T + p1 = np.matmul(C, datay) + p1x = np.polyval(p1, datax) + + # Step 5 - stablise process + delta = np.max(np.abs(p1x - p0x)) + p0 = p1 + p0x = p1x + + if iter_ == self.niter: + print('Maximum number of iterations reached') + break + + # Evaluate outputs + c = np.matmul(t, C) + yhat = np.matmul(c, datay) + pd = np.polyder(p0) + pdx = np.polyval(pd, datax[k]) + vk = np.power(self.jitterSD, 2) * np.power(pdx, 2) + np.power(self.noiseSD, 2) + R = np.power(np.linalg.norm(np.matmul(np.diagflat(w), (datay - np.matmul(X, p0)))), 2) + + return yhat, c, vk, R + + +def random_gaussian_whrand(M, mu, sigma, istate1, istate2): + ''' + Generates random numbers from a Gaussian Distribution using the Wichmann–Hill random number generator + + Inputs: + M: Number of random numbers required + mu: Mean of the Gaussian + sigma: Std of the Gaussian + istate1: vector of 4 integers + istate2: vector of 4 different integers + + Output: + xdist: Random Gaussian vector of size M + ''' + mn = 0 + ndist = np.zeros(M) + while mn < M: + rr1 = whrand(istate1, 1) + rr2 = whrand(istate2, 1) + + v1, istate1 = rr1 + v2, istate2 = rr2 + + if mn < M: + ndist[mn] = math.sqrt(-2 * math.log(v1)) * math.cos(2 * math.pi * v2) + mn = mn + 1 + + if mn < M: + ndist[mn] = math.sqrt(-2 * math.log(v1)) * math.sin(2 * math.pi * v2) + mn = mn + 1 + + xdist = mu + sigma * ndist + + return xdist, istate1, istate2 + + +def whrand(istate, N): + ''' + Generates uniform random numbers between 0 and 1 using the Wichmann–Hill random number generator + + Inputs: + istate: vector of 4 integers + N: Number of random numbers required + + Outputs: + r: Vector uniform random numbers of size M + istate: Output vector of 4 integers + ''' + + # Constants + a = np.array([11600, 47003, 23000, 33000]) + b = np.array([185127, 45688, 93368, 65075]) + c = np.array([10379, 10479, 19423, 8123]) + d = np.array([456, 420, 300, 0]) + 2147483123 + + r = np.zeros(N) + for i in range(N): + # Update states + for j in range(4): + istate[j] = a[j] * np.mod(istate[j], b[j]) - c[j] * np.fix(istate[j] / b[j]) + if istate[j] < 0: + istate[j] = istate[j] + d[j] + + # Evaluate random number + w = np.sum(np.divide(istate, d)) + r[i] = np.remainder(w, 1) + + return r, istate + + +def njr(fs, ydata, N, niter, tol, m0w, s0w, m0t, s0t, Mc, M0, Nc, Q): + analyse_fun = MCMCMH_NJ(fs, ydata, N, niter, tol, m0w, s0w, m0t, s0t, Mc, M0, Nc, Q) + yhat1= analyse_fun.AnalyseSignalN() + return yhat1 + +class NJRemoved(AgentMET4FOF): + def init_parameters(self, fs=100, ydata = np.array([]), N=15, niter=100, tol=1e-9, m0w = 10, s0w = 0.0005, m0t = 10, s0t = 0.0002*100/8, Mc=5000, M0=100, Nc=100, Q=50 ): + self.fs = fs + self.ydata = ydata + self.N = N + self.niter = niter + self.tol = tol + self.m0w = m0w + self.s0w = s0w + self.m0t = m0t + self.s0t = s0t + self.Mc = Mc + self.M0 = M0 + self.Nc = Nc + self.Q = Q + + + def on_received_message(self, message): + ddata = message['data'] + self.ydata = np.append(self.ydata, ddata) + if np.size(self.ydata) == self.N: + t = njr(self.fs, self.ydata, self.N, self.niter, self.tol, self.m0w, self.s0w, self.m0t, self.s0t, self.Mc, self.M0,self.Nc, self.Q) + self.send_output(self.ydata[7] - t) + self.ydata = self.ydata[1:self.N] + +def njr(fs, ydata, N, niter, tol, m0w, s0w, m0t, s0t, Mc, M0, Nc, Q): + analyse_fun = MCMCMH_NJ(fs, ydata, N, niter, tol, m0w, s0w, m0t, s0t, Mc, M0, Nc, Q) + yhat1= analyse_fun.AnalyseSignalN() + return yhat1 + + +######################################## +class NJRemoved(AgentMET4FOF): + def init_parameters(self, fs=100, ydata = np.array([]), N=15, niter=100, tol=1e-9, m0w = 10, s0w = 0.0005, m0t = 10, s0t = 0.0002*100/8, Mc=5000, M0=100, Nc=100, Q=50 ): + self.fs = fs + self.ydata = ydata + self.N = N + self.niter = niter + self.tol = tol + self.m0w = m0w + self.s0w = s0w + self.m0t = m0t + self.s0t = s0t + self.Mc = Mc + self.M0 = M0 + self.Nc = Nc + self.Q = Q + + + def on_received_message(self, message): + ddata = message['data'] + self.ydata = np.append(self.ydata, ddata) + if np.size(self.ydata) == self.N: + t = njr(self.fs, self.ydata, self.N, self.niter, self.tol, self.m0w, self.s0w, self.m0t, self.s0t, self.Mc, self.M0,self.Nc, self.Q) + self.send_output(self.ydata[7] - t) + self.ydata = self.ydata[1:self.N] + + +class SineGeneratorAgent(AgentMET4FOF): + def init_parameters(self): + self.stream = StaticSineGeneratorWithJitter(jittersd=0.0005, noisesd=0.0002) + + def agent_loop(self): + if self.current_state == "Running": + sine_data = self.stream.next_sample() # dictionary + self.send_output(sine_data['quantities']) + +def main(): + # start agent network server + agentNetwork = AgentNetwork(backend="mesa") + # init agents + + gen_agent = agentNetwork.add_agent(agentType=SineGeneratorAgent) + + njremove_agent = agentNetwork.add_agent(agentType=NJRemoved) + + monitor_agent = agentNetwork.add_agent(agentType=MonitorAgent) + monitor_agent2 = agentNetwork.add_agent(agentType=MonitorAgent) + + + # connect agents : We can connect multiple agents to any particular agent + + agentNetwork.bind_agents(gen_agent, njremove_agent) + + # connect + agentNetwork.bind_agents(gen_agent, monitor_agent) + + # agentNetwork.bind_agents(njremove_agent, moitor_agent2) + + agentNetwork.bind_agents(njremove_agent, monitor_agent2) + + # set all agents states to "Running" + agentNetwork.set_running_state() + + # allow for shutting down the network after execution + return agentNetwork + + +if __name__ == '__main__': + main() \ No newline at end of file From 844cf922068ef12566b4153e2c0ede89abf3bd7d Mon Sep 17 00:00:00 2001 From: Anupam Prasad Vedurmudi Date: Tue, 27 Jul 2021 14:33:05 +0200 Subject: [PATCH 07/43] wip(noise_jitter_removal_agent): re-added matlib import --- agentMET4FOF/agents/noise_jtter_removal_agents.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/agentMET4FOF/agents/noise_jtter_removal_agents.py b/agentMET4FOF/agents/noise_jtter_removal_agents.py index 67cbc9d0..97dc5e78 100644 --- a/agentMET4FOF/agents/noise_jtter_removal_agents.py +++ b/agentMET4FOF/agents/noise_jtter_removal_agents.py @@ -1,6 +1,7 @@ import math import numpy as np +import numpy.matlib from scipy.optimize import minimize from agentMET4FOF.agents import AgentMET4FOF, AgentNetwork, MonitorAgent @@ -826,7 +827,7 @@ def on_received_message(self, message): class SineGeneratorAgent(AgentMET4FOF): def init_parameters(self): - self.stream = StaticSineGeneratorWithJitter(jittersd=0.0005, noisesd=0.0002) + self.stream = StaticSineGeneratorWithJitter() def agent_loop(self): if self.current_state == "Running": From f949255e3a78f48dba86f95fae349220470e05fb Mon Sep 17 00:00:00 2001 From: Bjoern Ludwig Date: Fri, 30 Jul 2021 12:39:47 +0200 Subject: [PATCH 08/43] wip(noise_jitter): combine noise and jitter generation and removal into one tutorial package --- agentMET4FOF_tutorials/noise/__init__.py | 0 agentMET4FOF_tutorials/{jitter => noise_jitter}/__init__.py | 0 .../add_noise_to_existing_signal.py} | 0 .../generate_sine_with_jitter.py} | 0 .../remove_noise_and_jitter.py} | 0 agentMET4FOF_tutorials/noise_jitter_removal/__init__.py | 0 tests/test_tutorials.py | 4 ++-- 7 files changed, 2 insertions(+), 2 deletions(-) delete mode 100644 agentMET4FOF_tutorials/noise/__init__.py rename agentMET4FOF_tutorials/{jitter => noise_jitter}/__init__.py (100%) rename agentMET4FOF_tutorials/{noise/sine_with_noise.py => noise_jitter/add_noise_to_existing_signal.py} (100%) rename agentMET4FOF_tutorials/{jitter/sine_with_jitter.py => noise_jitter/generate_sine_with_jitter.py} (100%) rename agentMET4FOF_tutorials/{noise_jitter_removal/noise_jitter_removal.py => noise_jitter/remove_noise_and_jitter.py} (100%) delete mode 100644 agentMET4FOF_tutorials/noise_jitter_removal/__init__.py diff --git a/agentMET4FOF_tutorials/noise/__init__.py b/agentMET4FOF_tutorials/noise/__init__.py deleted file mode 100644 index e69de29b..00000000 diff --git a/agentMET4FOF_tutorials/jitter/__init__.py b/agentMET4FOF_tutorials/noise_jitter/__init__.py similarity index 100% rename from agentMET4FOF_tutorials/jitter/__init__.py rename to agentMET4FOF_tutorials/noise_jitter/__init__.py diff --git a/agentMET4FOF_tutorials/noise/sine_with_noise.py b/agentMET4FOF_tutorials/noise_jitter/add_noise_to_existing_signal.py similarity index 100% rename from agentMET4FOF_tutorials/noise/sine_with_noise.py rename to agentMET4FOF_tutorials/noise_jitter/add_noise_to_existing_signal.py diff --git a/agentMET4FOF_tutorials/jitter/sine_with_jitter.py b/agentMET4FOF_tutorials/noise_jitter/generate_sine_with_jitter.py similarity index 100% rename from agentMET4FOF_tutorials/jitter/sine_with_jitter.py rename to agentMET4FOF_tutorials/noise_jitter/generate_sine_with_jitter.py diff --git a/agentMET4FOF_tutorials/noise_jitter_removal/noise_jitter_removal.py b/agentMET4FOF_tutorials/noise_jitter/remove_noise_and_jitter.py similarity index 100% rename from agentMET4FOF_tutorials/noise_jitter_removal/noise_jitter_removal.py rename to agentMET4FOF_tutorials/noise_jitter/remove_noise_and_jitter.py diff --git a/agentMET4FOF_tutorials/noise_jitter_removal/__init__.py b/agentMET4FOF_tutorials/noise_jitter_removal/__init__.py deleted file mode 100644 index e69de29b..00000000 diff --git a/tests/test_tutorials.py b/tests/test_tutorials.py index 54758be5..e1ff7d0a 100644 --- a/tests/test_tutorials.py +++ b/tests/test_tutorials.py @@ -1,7 +1,7 @@ -from agentMET4FOF_tutorials.jitter.sine_with_jitter import ( +from agentMET4FOF_tutorials.noise_jitter.generate_sine_with_jitter import ( demonstrate_sine_with_jitter_agent_use, ) -from agentMET4FOF_tutorials.noise.sine_with_noise import ( +from agentMET4FOF_tutorials.noise_jitter.add_noise_to_existing_signal import ( demonstrate_noise_agent_use, ) from agentMET4FOF_tutorials.redundancy.redundancy_agent_four_signals import ( From e64b5d4e255f2ef00d36c76bf1089bb2031ac6ca Mon Sep 17 00:00:00 2001 From: Bjoern Ludwig Date: Fri, 30 Jul 2021 12:46:22 +0200 Subject: [PATCH 09/43] wip(noise_jitter_removal_agents): fix typo in module name --- ...se_jtter_removal_agents.py => noise_jitter_removal_agents.py} | 1 - 1 file changed, 1 deletion(-) rename agentMET4FOF/agents/{noise_jtter_removal_agents.py => noise_jitter_removal_agents.py} (99%) diff --git a/agentMET4FOF/agents/noise_jtter_removal_agents.py b/agentMET4FOF/agents/noise_jitter_removal_agents.py similarity index 99% rename from agentMET4FOF/agents/noise_jtter_removal_agents.py rename to agentMET4FOF/agents/noise_jitter_removal_agents.py index 97dc5e78..b853fdb2 100644 --- a/agentMET4FOF/agents/noise_jtter_removal_agents.py +++ b/agentMET4FOF/agents/noise_jitter_removal_agents.py @@ -1,5 +1,4 @@ import math - import numpy as np import numpy.matlib from scipy.optimize import minimize From 59e8da1c75058f2284ee08a9ab5a6f036f5a49a9 Mon Sep 17 00:00:00 2001 From: Bjoern Ludwig Date: Fri, 30 Jul 2021 12:47:22 +0200 Subject: [PATCH 10/43] wip(noise_jitter_removal_agents): remove unneeded imports, classes and functions --- .../agents/noise_jitter_removal_agents.py | 47 +------------------ 1 file changed, 1 insertion(+), 46 deletions(-) diff --git a/agentMET4FOF/agents/noise_jitter_removal_agents.py b/agentMET4FOF/agents/noise_jitter_removal_agents.py index b853fdb2..b8cfcc30 100644 --- a/agentMET4FOF/agents/noise_jitter_removal_agents.py +++ b/agentMET4FOF/agents/noise_jitter_removal_agents.py @@ -3,8 +3,7 @@ import numpy.matlib from scipy.optimize import minimize -from agentMET4FOF.agents import AgentMET4FOF, AgentNetwork, MonitorAgent -from agentMET4FOF.streams.signal_streams import StaticSineGeneratorWithJitter +from agentMET4FOF.agents import AgentMET4FOF """### mcmcci: MCMC convergence indices for multiple chains.""" @@ -822,47 +821,3 @@ def on_received_message(self, message): t = njr(self.fs, self.ydata, self.N, self.niter, self.tol, self.m0w, self.s0w, self.m0t, self.s0t, self.Mc, self.M0,self.Nc, self.Q) self.send_output(self.ydata[7] - t) self.ydata = self.ydata[1:self.N] - - -class SineGeneratorAgent(AgentMET4FOF): - def init_parameters(self): - self.stream = StaticSineGeneratorWithJitter() - - def agent_loop(self): - if self.current_state == "Running": - sine_data = self.stream.next_sample() # dictionary - self.send_output(sine_data['quantities']) - -def main(): - # start agent network server - agentNetwork = AgentNetwork(backend="mesa") - # init agents - - gen_agent = agentNetwork.add_agent(agentType=SineGeneratorAgent) - - njremove_agent = agentNetwork.add_agent(agentType=NJRemoved) - - monitor_agent = agentNetwork.add_agent(agentType=MonitorAgent) - monitor_agent2 = agentNetwork.add_agent(agentType=MonitorAgent) - - - # connect agents : We can connect multiple agents to any particular agent - - agentNetwork.bind_agents(gen_agent, njremove_agent) - - # connect - agentNetwork.bind_agents(gen_agent, monitor_agent) - - # agentNetwork.bind_agents(njremove_agent, moitor_agent2) - - agentNetwork.bind_agents(njremove_agent, monitor_agent2) - - # set all agents states to "Running" - agentNetwork.set_running_state() - - # allow for shutting down the network after execution - return agentNetwork - - -if __name__ == '__main__': - main() \ No newline at end of file From 009ed8e8f14447692f1cdbbe7ccaa94b1ac12c53 Mon Sep 17 00:00:00 2001 From: Anupam Prasad Vedurmudi Date: Fri, 30 Jul 2021 13:36:56 +0200 Subject: [PATCH 11/43] wip(noise_jitter_removal_agents): brings module up-to-date and removes duplicated classes --- .../agents/noise_jitter_removal_agents.py | 31 ------------------- 1 file changed, 31 deletions(-) diff --git a/agentMET4FOF/agents/noise_jitter_removal_agents.py b/agentMET4FOF/agents/noise_jitter_removal_agents.py index b8cfcc30..aeb1a938 100644 --- a/agentMET4FOF/agents/noise_jitter_removal_agents.py +++ b/agentMET4FOF/agents/noise_jitter_removal_agents.py @@ -759,37 +759,6 @@ def whrand(istate, N): return r, istate - -def njr(fs, ydata, N, niter, tol, m0w, s0w, m0t, s0t, Mc, M0, Nc, Q): - analyse_fun = MCMCMH_NJ(fs, ydata, N, niter, tol, m0w, s0w, m0t, s0t, Mc, M0, Nc, Q) - yhat1= analyse_fun.AnalyseSignalN() - return yhat1 - -class NJRemoved(AgentMET4FOF): - def init_parameters(self, fs=100, ydata = np.array([]), N=15, niter=100, tol=1e-9, m0w = 10, s0w = 0.0005, m0t = 10, s0t = 0.0002*100/8, Mc=5000, M0=100, Nc=100, Q=50 ): - self.fs = fs - self.ydata = ydata - self.N = N - self.niter = niter - self.tol = tol - self.m0w = m0w - self.s0w = s0w - self.m0t = m0t - self.s0t = s0t - self.Mc = Mc - self.M0 = M0 - self.Nc = Nc - self.Q = Q - - - def on_received_message(self, message): - ddata = message['data'] - self.ydata = np.append(self.ydata, ddata) - if np.size(self.ydata) == self.N: - t = njr(self.fs, self.ydata, self.N, self.niter, self.tol, self.m0w, self.s0w, self.m0t, self.s0t, self.Mc, self.M0,self.Nc, self.Q) - self.send_output(self.ydata[7] - t) - self.ydata = self.ydata[1:self.N] - def njr(fs, ydata, N, niter, tol, m0w, s0w, m0t, s0t, Mc, M0, Nc, Q): analyse_fun = MCMCMH_NJ(fs, ydata, N, niter, tol, m0w, s0w, m0t, s0t, Mc, M0, Nc, Q) yhat1= analyse_fun.AnalyseSignalN() From c37a47092259c0fe2d75dd3c345adcde5ce4805d Mon Sep 17 00:00:00 2001 From: Anupam Prasad Vedurmudi Date: Fri, 30 Jul 2021 13:53:16 +0200 Subject: [PATCH 12/43] wip(noise_jitter_removal_agents): cleans up noise jitter removal tutorial --- .../agents/noise_jitter_removal_agents.py | 2 +- .../noise_jitter/remove_noise_and_jitter.py | 58 ++++--------------- 2 files changed, 13 insertions(+), 47 deletions(-) diff --git a/agentMET4FOF/agents/noise_jitter_removal_agents.py b/agentMET4FOF/agents/noise_jitter_removal_agents.py index aeb1a938..c2d1b629 100644 --- a/agentMET4FOF/agents/noise_jitter_removal_agents.py +++ b/agentMET4FOF/agents/noise_jitter_removal_agents.py @@ -766,7 +766,7 @@ def njr(fs, ydata, N, niter, tol, m0w, s0w, m0t, s0t, Mc, M0, Nc, Q): ######################################## -class NJRemoved(AgentMET4FOF): +class NoiseJitterRemovalAgent(AgentMET4FOF): def init_parameters(self, fs=100, ydata = np.array([]), N=15, niter=100, tol=1e-9, m0w = 10, s0w = 0.0005, m0t = 10, s0t = 0.0002*100/8, Mc=5000, M0=100, Nc=100, Q=50 ): self.fs = fs self.ydata = ydata diff --git a/agentMET4FOF_tutorials/noise_jitter/remove_noise_and_jitter.py b/agentMET4FOF_tutorials/noise_jitter/remove_noise_and_jitter.py index a8e84aad..6ab7ab36 100644 --- a/agentMET4FOF_tutorials/noise_jitter/remove_noise_and_jitter.py +++ b/agentMET4FOF_tutorials/noise_jitter/remove_noise_and_jitter.py @@ -1,63 +1,29 @@ -from agentMET4FOF.agents import AgentMET4FOF, AgentNetwork, MonitorAgent -import numpy as np -from NJRemove.NJRemoval_class_withmcmc import MCMCMH_NJ +from agentMET4FOF.agents.base_agents import MonitorAgent +from agentMET4FOF.agents.noise_jitter_removal_agents import NoiseJitterRemovalAgent +from agentMET4FOF.agents.signal_agents import StaticSineWithJitterGeneratorAgent, NoiseAgent +from agentMET4FOF.network import AgentNetwork -from agentMET4FOF.agents.signal_agents import StaticSineGeneratorWithJitterAgent, NoiseAgent - -def njr(fs, ydata, N, niter, tol, m0w, s0w, m0t, s0t, Mc, M0, Nc, Q): - analyse_fun = MCMCMH_NJ(fs, ydata, N, niter, tol, m0w, s0w, m0t, s0t, Mc, M0, Nc, Q) - yhat1= analyse_fun.AnalyseSignalN() - return yhat1 - - -######################################## -class NJRemoved(AgentMET4FOF): - def init_parameters(self, fs=100, ydata = np.array([]), N=15, niter=100, tol=1e-9, m0w = 10, s0w = 0.0005, m0t = 10, s0t = 0.0002*100/8, Mc=5000, M0=100, Nc=100, Q=50 ): - self.fs = fs - self.ydata = ydata - self.N = N - self.niter = niter - self.tol = tol - self.m0w = m0w - self.s0w = s0w - self.m0t = m0t - self.s0t = s0t - self.Mc = Mc - self.M0 = M0 - self.Nc = Nc - self.Q = Q - - - def on_received_message(self, message): - ddata = message['data'] - self.ydata = np.append(self.ydata, ddata) - if np.size(self.ydata) == self.N: - t = njr(self.fs, self.ydata, self.N, self.niter, self.tol, self.m0w, self.s0w, self.m0t, self.s0t, self.Mc, self.M0,self.Nc, self.Q) - self.send_output(self.ydata[7] - t) - self.ydata = self.ydata[1:self.N] - -def main(): +def demonstrate_noise_jitter_removal_agent(): # start agent network server agentNetwork = AgentNetwork(backend="mesa") # init agents - jitter_gen_agent = agentNetwork.add_agent(agentType=StaticSineGeneratorWithJitterAgent) + sine_with_jitter_agent = agentNetwork.add_agent(agentType=StaticSineWithJitterGeneratorAgent) noise_agent = agentNetwork.add_agent(agentType=NoiseAgent) - njremove_agent = agentNetwork.add_agent(agentType=NJRemoved) - - monitor_agent = agentNetwork.add_agent(agentType=MonitorAgent) - monitor_agent2 = agentNetwork.add_agent(agentType=MonitorAgent) + njremove_agent = agentNetwork.add_agent(agentType=NoiseJitterRemovalAgent) + monitor_agent = agentNetwork.add_agent(agentType=MonitorAgent, name="Sine with Noise and Jitter") + monitor_agent2 = agentNetwork.add_agent(agentType=MonitorAgent, name="Output of Noise-Jitter Removal Agent") # connect agents : jitter generator -> noise -> njremoval agent - agentNetwork.bind_agents(jitter_gen_agent, noise_agent) + agentNetwork.bind_agents(sine_with_jitter_agent, noise_agent) agentNetwork.bind_agents(noise_agent, njremove_agent) # connect monitor agents - agentNetwork.bind_agents(jitter_gen_agent, monitor_agent) + agentNetwork.bind_agents(noise_agent, monitor_agent) agentNetwork.bind_agents(njremove_agent, monitor_agent2) # set all agents states to "Running" @@ -68,4 +34,4 @@ def main(): if __name__ == '__main__': - main() \ No newline at end of file + demonstrate_noise_jitter_removal_agent() From 91950684fe6dbd90d325d679ea458337ce6700fa Mon Sep 17 00:00:00 2001 From: Anupam Prasad Vedurmudi Date: Fri, 30 Jul 2021 14:29:07 +0200 Subject: [PATCH 13/43] wip(noise_jitter_removal_agents): cleans up noise_jitter_removal_agents.py partially --- .../agents/noise_jitter_removal_agents.py | 320 ++++++------------ 1 file changed, 111 insertions(+), 209 deletions(-) diff --git a/agentMET4FOF/agents/noise_jitter_removal_agents.py b/agentMET4FOF/agents/noise_jitter_removal_agents.py index c2d1b629..c217bc1a 100644 --- a/agentMET4FOF/agents/noise_jitter_removal_agents.py +++ b/agentMET4FOF/agents/noise_jitter_removal_agents.py @@ -75,7 +75,6 @@ def mcmcci(A,M0): return Rhat, Neff """### mcsums: Summary information from MC samples.""" - def mcsums(A,M0,Q): ''' mcsums: Summary information from MC samples. @@ -123,47 +122,7 @@ def mcsums(A,M0,Q): return abar, s, aQ -"""### jumprwg: Jumping distribution for the Metropolis Hastings Gaussian random walk algorithm""" - -def jumprwg(A, L): - ''' - jumprwg: Jumping distribution for the Metropolis Hastings Gaussian random - walk algorithm - ------------------------------------------------------------------------- - KJ, LRW, PMH - Version 2020-04-22 - ------------------------------------------------------------------------- - Inputs: - A(n,N): Samples at the current iteration - - L(n,n): Cholesky factor of variance of parameter vector. - - Outputs: - As(n,N): Proposed parameter array which is randomly sampled from the - jumping distribution - - dp0: The difference between the logarithm of the jumping distribution - associated with moving from A(:,j) to As(:,j) and that associated with - moving from As(:,j) to A(:,j), up to an additive constant. - log P0(a|as) - log P0(as|a) - - ''' - # Number of parameters and parallel chains - n, N = A.shape - - # random draw from a Gaussian distribution - e = np.random.normal(0, 1, size=(n,N)) - - # proposed draw from a Gaussian distribution with mean A and variance LL' - As = A + np.matmul(L,e) - - # For a Gaussian random walk, since log P0(a|as) = log P0(as|a), dp0 will always be zero - dp0 = np.zeros(N) - - return As, dp0 - """### Cubic function and its first and second derivative""" - def fgh_cubic(alpha,t): ''' ------------------------------------------------------------------------- @@ -243,65 +202,7 @@ def ln_gauss_pdf_v(x,mu,sigma): return logf -"""### Target dist for noise and jitter posterior dist""" - -def tar_at(at, y, x, m0w, s0w, m0t, s0t): - ''' - ------------------------------------------------------------------------- - Target dist for noise and jitter posterior dist - ------------------------------------------------------------------------- - KJ, LRW, PMH - Version 2020-04-22 - -------------------------------------------------------------------------- - Inputs: - at(n+2,N): Parameters alpha, log(1/tau^2) and log(1/w^2) - - y(m,1): Signal - - x(m,1): time at which signal was recorded - - s0w and s0t: prior estimates of tau and omega - - m0w and m0t: degree of belief in prior estimates for tau and omega - - Output: - T: Log of the posterior distribution - ''' - - # Size of parameter vector - at = np.array(at) - p = at.shape[0] - - # Number of alphas - n = p - 2 - - # Extract parameters - alpha = at[0:n] - phi1 = np.exp(at[-2]) - phi2 = np.exp(at[-1]) - taus = np.ones(phi1.shape)/np.sqrt(phi1) - omegas = np.ones(phi2.shape)/np.sqrt(phi2) - - # Gamma priors for phis - prior_phi1 = (m0t/2)*np.log(phi1) - phi1*m0t*s0t**2/2 - prior_phi2 = (m0w/2)*np.log(phi2) - phi2*m0w*s0w**2/2 - - - # function that evaluates the cubic function with user specified cubic parameters - fun = lambda aa: fgh_cubic(aa, x) - - # cubic, expectation and variance - [st,st1,st2] = fun(alpha) - expect = st + 0.5*(taus**2)*st2 - vari = (taus**2)*(st1**2) + omegas**2 - - # Likelihood - lik = sum(ln_gauss_pdf_v(y,expect,np.sqrt(vari))) - - # Posterior - T = lik + prior_phi1 + prior_phi2 - return T """###mcmcmh: Metrolopolis-Hasting MCMC algorithm generating N chains of length M for a parameter vector A of length n.""" @@ -427,57 +328,6 @@ def mcmcmh(M, N, M0, Q, A0, tar, jump): return S, aP, Rh, Ne, AA, IAA -def mcmcm_main(datay, datax, m0w, s0w, m0t, s0t, Mc, M0, Nc, Q): - - at0 = np.array((1,1,1,1, np.log(1/s0w**2), np.log(1/s0t**2))) - # function that evaluates the log of the target distribution at given parameter values - tar = lambda at: tar_at(at, datay, datax, m0w, s0w, m0t, s0t) - # function that evaluates the negative log of the target distribution to evaluate MAP estimates - mapp = lambda at: -tar(at) - - res = minimize(mapp, at0) - pars = res.x - V = res.hess_inv - L = np.linalg.cholesky(V) - - - # Function that draws sample from a Gaussian random walk jumping distribution - jump = lambda A: jumprwg(A, L) - - - rr = np.random.normal(0,1,size=(6,Nc)) - - A0 = np.matlib.repmat(pars.T,Nc,1).T + np.matmul(L,rr) - - - sam = mcmcmh(Mc,Nc,M0,Q,A0,tar,jump) - return 1/np.sqrt(np.exp(sam[0][0,-2:])) - -# Decaying exponential -def DecayExpFunction(a, b, f, x): - ''' - decaying exponential function evaluated at x with parameters a, b and f - ''' - return a * np.exp(-b * x) * np.sin(2 * np.pi * f * x) - - -# Decaying exponential - first derivative -def DecayExpFunction1der(a, b, f, x): - ''' - first derviative of the decaying exponential function evaluated at x with parameters a, b and f - ''' - return a * np.exp(-b * x) * ((2 * np.pi * f) * np.cos(2 * np.pi * f * x) - b * np.sin(2 * np.pi * f * x)) - - -# Decaying exponential - second derivative -def DecayExpFunction2der(a, b, f, x): - ''' - second derviative of the decaying exponential function evaluated at x with parameters a, b and f - ''' - return a * np.exp(-b * x) * ((-np.power((2 * np.pi * f), 2)) * np.sin(2 * np.pi * f * x) - - ((4 * b * np.pi * f) * np.cos(2 * np.pi * f * x)) + b * np.sin(2 * np.pi * f * x)) - - class MCMCMH_NJ(): ''' Bayesian Noise and jitter reduction algorithm. MCMC used to determine the noise and jitter variances. @@ -537,7 +387,7 @@ def AnalyseSignalN(self): # Values of normalised independent variables datax = np.divide(np.arange(-n, n + 1), self.fs) - outs = mcmcm_main(self.ydata, datax, self.m0w, self.s0w, self.m0t, self.s0t, self.Mc, + outs = self.mcmcm_main(self.ydata, datax, self.m0w, self.s0w, self.m0t, self.s0t, self.Mc, self.M0, self.Nc, self.Q) self.jitterSD = outs[0] self.noiseSD = outs[1] @@ -689,81 +539,128 @@ def NJAlgorithm(self, datax, datay, p0, p0x): return yhat, c, vk, R + @staticmethod + def mcmcm_main(datay, datax, m0w, s0w, m0t, s0t, Mc, M0, Nc, Q): -def random_gaussian_whrand(M, mu, sigma, istate1, istate2): - ''' - Generates random numbers from a Gaussian Distribution using the Wichmann–Hill random number generator + at0 = np.array((1, 1, 1, 1, np.log(1 / s0w ** 2), np.log(1 / s0t ** 2))) + # function that evaluates the log of the target distribution at given parameter values + tar = lambda at: MCMCMH_NJ.tar_at(at, datay, datax, m0w, s0w, m0t, s0t) + # function that evaluates the negative log of the target distribution to evaluate MAP estimates + mapp = lambda at: -tar(at) - Inputs: - M: Number of random numbers required - mu: Mean of the Gaussian - sigma: Std of the Gaussian - istate1: vector of 4 integers - istate2: vector of 4 different integers + res = minimize(mapp, at0) + pars = res.x + V = res.hess_inv + L = np.linalg.cholesky(V) - Output: - xdist: Random Gaussian vector of size M - ''' - mn = 0 - ndist = np.zeros(M) - while mn < M: - rr1 = whrand(istate1, 1) - rr2 = whrand(istate2, 1) + # Function that draws sample from a Gaussian random walk jumping distribution + jump = lambda A: MCMCMH_NJ.jumprwg(A, L) - v1, istate1 = rr1 - v2, istate2 = rr2 + rr = np.random.normal(0, 1, size=(6, Nc)) - if mn < M: - ndist[mn] = math.sqrt(-2 * math.log(v1)) * math.cos(2 * math.pi * v2) - mn = mn + 1 + A0 = np.matlib.repmat(pars.T, Nc, 1).T + np.matmul(L, rr) - if mn < M: - ndist[mn] = math.sqrt(-2 * math.log(v1)) * math.sin(2 * math.pi * v2) - mn = mn + 1 + sam = mcmcmh(Mc, Nc, M0, Q, A0, tar, jump) + return 1 / np.sqrt(np.exp(sam[0][0, -2:])) - xdist = mu + sigma * ndist + """### Target dist for noise and jitter posterior dist""" - return xdist, istate1, istate2 + @staticmethod + def tar_at(at, y, x, m0w, s0w, m0t, s0t): + ''' + ------------------------------------------------------------------------- + Target dist for noise and jitter posterior dist + ------------------------------------------------------------------------- + KJ, LRW, PMH + Version 2020-04-22 + -------------------------------------------------------------------------- + Inputs: + at(n+2,N): Parameters alpha, log(1/tau^2) and log(1/w^2) + y(m,1): Signal -def whrand(istate, N): - ''' - Generates uniform random numbers between 0 and 1 using the Wichmann–Hill random number generator + x(m,1): time at which signal was recorded - Inputs: - istate: vector of 4 integers - N: Number of random numbers required + s0w and s0t: prior estimates of tau and omega - Outputs: - r: Vector uniform random numbers of size M - istate: Output vector of 4 integers - ''' + m0w and m0t: degree of belief in prior estimates for tau and omega + + Output: + T: Log of the posterior distribution + ''' - # Constants - a = np.array([11600, 47003, 23000, 33000]) - b = np.array([185127, 45688, 93368, 65075]) - c = np.array([10379, 10479, 19423, 8123]) - d = np.array([456, 420, 300, 0]) + 2147483123 + # Size of parameter vector + at = np.array(at) + p = at.shape[0] - r = np.zeros(N) - for i in range(N): - # Update states - for j in range(4): - istate[j] = a[j] * np.mod(istate[j], b[j]) - c[j] * np.fix(istate[j] / b[j]) - if istate[j] < 0: - istate[j] = istate[j] + d[j] + # Number of alphas + n = p - 2 + + # Extract parameters + alpha = at[0:n] + phi1 = np.exp(at[-2]) + phi2 = np.exp(at[-1]) + taus = np.ones(phi1.shape) / np.sqrt(phi1) + omegas = np.ones(phi2.shape) / np.sqrt(phi2) + + # Gamma priors for phis + prior_phi1 = (m0t / 2) * np.log(phi1) - phi1 * m0t * s0t ** 2 / 2 + prior_phi2 = (m0w / 2) * np.log(phi2) - phi2 * m0w * s0w ** 2 / 2 + + # function that evaluates the cubic function with user specified cubic parameters + fun = lambda aa: fgh_cubic(aa, x) + + # cubic, expectation and variance + [st, st1, st2] = fun(alpha) + expect = st + 0.5 * (taus ** 2) * st2 + vari = (taus ** 2) * (st1 ** 2) + omegas ** 2 + + # Likelihood + lik = sum(ln_gauss_pdf_v(y, expect, np.sqrt(vari))) + + # Posterior + T = lik + prior_phi1 + prior_phi2 + + return T + + """### jumprwg: Jumping distribution for the Metropolis Hastings Gaussian random walk algorithm""" + @staticmethod + def jumprwg(A, L): + ''' + jumprwg: Jumping distribution for the Metropolis Hastings Gaussian random + walk algorithm + ------------------------------------------------------------------------- + KJ, LRW, PMH + Version 2020-04-22 + ------------------------------------------------------------------------- + Inputs: + A(n,N): Samples at the current iteration + + L(n,n): Cholesky factor of variance of parameter vector. + + Outputs: + As(n,N): Proposed parameter array which is randomly sampled from the + jumping distribution + + dp0: The difference between the logarithm of the jumping distribution + associated with moving from A(:,j) to As(:,j) and that associated with + moving from As(:,j) to A(:,j), up to an additive constant. + log P0(a|as) - log P0(as|a) + + ''' + # Number of parameters and parallel chains + n, N = A.shape - # Evaluate random number - w = np.sum(np.divide(istate, d)) - r[i] = np.remainder(w, 1) + # random draw from a Gaussian distribution + e = np.random.normal(0, 1, size=(n, N)) - return r, istate + # proposed draw from a Gaussian distribution with mean A and variance LL' + As = A + np.matmul(L, e) -def njr(fs, ydata, N, niter, tol, m0w, s0w, m0t, s0t, Mc, M0, Nc, Q): - analyse_fun = MCMCMH_NJ(fs, ydata, N, niter, tol, m0w, s0w, m0t, s0t, Mc, M0, Nc, Q) - yhat1= analyse_fun.AnalyseSignalN() - return yhat1 + # For a Gaussian random walk, since log P0(a|as) = log P0(as|a), dp0 will always be zero + dp0 = np.zeros(N) + return As, dp0 ######################################## class NoiseJitterRemovalAgent(AgentMET4FOF): @@ -782,11 +679,16 @@ def init_parameters(self, fs=100, ydata = np.array([]), N=15, niter=100, tol=1e self.Nc = Nc self.Q = Q + @staticmethod + def njr(fs, ydata, N, niter, tol, m0w, s0w, m0t, s0t, Mc, M0, Nc, Q) -> object: + analyse_fun = MCMCMH_NJ(fs, ydata, N, niter, tol, m0w, s0w, m0t, s0t, Mc, M0, Nc, Q) + yhat1 = analyse_fun.AnalyseSignalN() + return yhat1 def on_received_message(self, message): ddata = message['data'] self.ydata = np.append(self.ydata, ddata) if np.size(self.ydata) == self.N: - t = njr(self.fs, self.ydata, self.N, self.niter, self.tol, self.m0w, self.s0w, self.m0t, self.s0t, self.Mc, self.M0,self.Nc, self.Q) + t = self.njr(self.fs, self.ydata, self.N, self.niter, self.tol, self.m0w, self.s0w, self.m0t, self.s0t, self.Mc, self.M0,self.Nc, self.Q) self.send_output(self.ydata[7] - t) self.ydata = self.ydata[1:self.N] From 35ebe670b7d1f46fc8ef71a2d4ce3f94316d2d5b Mon Sep 17 00:00:00 2001 From: Bjoern Ludwig Date: Fri, 30 Jul 2021 14:30:09 +0200 Subject: [PATCH 14/43] wip(NoiseJitterRemovalAgent): when working with signals that contain time stamps or targets we only use the quantities --- agentMET4FOF/agents/noise_jitter_removal_agents.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/agentMET4FOF/agents/noise_jitter_removal_agents.py b/agentMET4FOF/agents/noise_jitter_removal_agents.py index c217bc1a..32a337e1 100644 --- a/agentMET4FOF/agents/noise_jitter_removal_agents.py +++ b/agentMET4FOF/agents/noise_jitter_removal_agents.py @@ -686,7 +686,10 @@ def njr(fs, ydata, N, niter, tol, m0w, s0w, m0t, s0t, Mc, M0, Nc, Q) -> object: return yhat1 def on_received_message(self, message): - ddata = message['data'] + if isinstance(message["data"], dict): + ddata = message["data"]["quantities"] + else: + ddata = message["data"] self.ydata = np.append(self.ydata, ddata) if np.size(self.ydata) == self.N: t = self.njr(self.fs, self.ydata, self.N, self.niter, self.tol, self.m0w, self.s0w, self.m0t, self.s0t, self.Mc, self.M0,self.Nc, self.Q) From 54dd135078e21264a971b464a9aa1880399da406 Mon Sep 17 00:00:00 2001 From: Bjoern Ludwig Date: Fri, 30 Jul 2021 14:31:10 +0200 Subject: [PATCH 15/43] wip(noise_and_jitter): blacken the two corresponding files --- .../agents/noise_jitter_removal_agents.py | 256 +++++++++++------- .../noise_jitter/remove_noise_and_jitter.py | 19 +- 2 files changed, 177 insertions(+), 98 deletions(-) diff --git a/agentMET4FOF/agents/noise_jitter_removal_agents.py b/agentMET4FOF/agents/noise_jitter_removal_agents.py index 32a337e1..b2509b75 100644 --- a/agentMET4FOF/agents/noise_jitter_removal_agents.py +++ b/agentMET4FOF/agents/noise_jitter_removal_agents.py @@ -7,8 +7,9 @@ """### mcmcci: MCMC convergence indices for multiple chains.""" -def mcmcci(A,M0): - ''' + +def mcmcci(A, M0): + """ mcmcci: MCMC convergence indices for multiple chains. ---------------------------------------------------------------------------- KJ, LRW, PMH @@ -32,7 +33,7 @@ def mcmcci(A,M0): set to (M-M0)*N, their limit values. Note: If N = 1 or M0 > M-2, Rhat = 0; Neff = 0. - ''' + """ A = np.array(A) M, N = A.shape @@ -48,35 +49,38 @@ def mcmcci(A,M0): Mt = M - M0 # Chain mean and mean of means - asub = A[M0:,:] - ad = np.mean(asub,axis = 0) + asub = A[M0:, :] + ad = np.mean(asub, axis=0) add = asub.mean() # Within group standard deviation - ss = np.std(asub,axis = 0) + ss = np.std(asub, axis=0) # Between groups variance. dd = np.square(ad - add) - B = (Mt*np.sum(dd))/(N-1) + B = (Mt * np.sum(dd)) / (N - 1) # Within groups variance. - W = np.sum(np.square(ss))/N + W = np.sum(np.square(ss)) / N # V plus - Vp = (Mt-1)*W/Mt + B/Mt + Vp = (Mt - 1) * W / Mt + B / Mt # Convergence statistic, effective number of independent samples - Rhat = np.sqrt(Vp/W) - Neff = Mt*N*Vp/B + Rhat = np.sqrt(Vp / W) + Neff = Mt * N * Vp / B - Rhat = np.maximum(Rhat,1) - Neff = np.minimum(Neff,Mt*N) + Rhat = np.maximum(Rhat, 1) + Neff = np.minimum(Neff, Mt * N) return Rhat, Neff + """### mcsums: Summary information from MC samples.""" -def mcsums(A,M0,Q): - ''' + + +def mcsums(A, M0, Q): + """ mcsums: Summary information from MC samples. ------------------------------------------------------------------------- KJ, LRW, PMH @@ -97,13 +101,13 @@ def mcsums(A,M0,Q): aQ(nQ,n): Percentiles corresponding to Q. - ''' + """ # Size of samples after burn-in A = np.array(A) M, N = A.shape - m = (M - M0)*N + m = (M - M0) * N # Initialise percentile vector nQ = Q.size @@ -118,13 +122,16 @@ def mcsums(A,M0,Q): s = np.std(aaj) # Percentiles of samples - aQ = np.percentile(aaj,Q) + aQ = np.percentile(aaj, Q) return abar, s, aQ + """### Cubic function and its first and second derivative""" -def fgh_cubic(alpha,t): - ''' + + +def fgh_cubic(alpha, t): + """ ------------------------------------------------------------------------- Cubic function and its first and second derivative ------------------------------------------------------------------------- @@ -142,29 +149,31 @@ def fgh_cubic(alpha,t): f1(m,N): Derivative of cubic f2(m,N): Second derivative of cubic - ''' + """ # length of data and number of paramaters m = t.size # design matrix - C = np.array([np.ones(m), t, t**2, t**3]) + C = np.array([np.ones(m), t, t ** 2, t ** 3]) # derivate info - C1 = np.array([np.ones(m), 2*t, 3*t**2]) - C2 = np.array([2*np.ones(m), 6*t]) + C1 = np.array([np.ones(m), 2 * t, 3 * t ** 2]) + C2 = np.array([2 * np.ones(m), 6 * t]) # cubic and derivatives - f = np.matmul(C.T,alpha) - f1 = np.matmul(C1.T,alpha[1:]) - f2 = np.matmul(C2.T,alpha[2:]) + f = np.matmul(C.T, alpha) + f1 = np.matmul(C1.T, alpha[1:]) + f2 = np.matmul(C2.T, alpha[2:]) return f, f1, f2 + """### Log of the gaussian pdf""" -def ln_gauss_pdf_v(x,mu,sigma): - ''' + +def ln_gauss_pdf_v(x, mu, sigma): + """ ------------------------------------------------------------------------- Log of the Gaussian pdf ------------------------------------------------------------------------- @@ -181,33 +190,32 @@ def ln_gauss_pdf_v(x,mu,sigma): Output: logf: Log of the Gaussian pdf at x with mean mu and std sigma - ''' - + """ try: - # When inputs are high dimensional arrays/matrices - xx = np.matlib.repmat(x,mu.shape[1],1) - xx = xx.T - logk = - np.log(2*math.pi)/2 - np.log(sigma) - logp = -((xx - mu)**2)/(2*sigma**2) - # Log of the Gaussian PDF - logf = logk + logp + # When inputs are high dimensional arrays/matrices + xx = np.matlib.repmat(x, mu.shape[1], 1) + xx = xx.T + logk = -np.log(2 * math.pi) / 2 - np.log(sigma) + logp = -((xx - mu) ** 2) / (2 * sigma ** 2) + # Log of the Gaussian PDF + logf = logk + logp except IndexError: - # When inputs are vectors - logk = - np.log(2*math.pi)/2 - np.log(sigma) - logp = -((x - mu)**2)/(2*sigma**2) - # Log of the Gaussian PDF - logf = logk + logp + # When inputs are vectors + logk = -np.log(2 * math.pi) / 2 - np.log(sigma) + logp = -((x - mu) ** 2) / (2 * sigma ** 2) + # Log of the Gaussian PDF + logf = logk + logp return logf - """###mcmcmh: Metrolopolis-Hasting MCMC algorithm generating N chains of length M for a parameter vector A of length n.""" + def mcmcmh(M, N, M0, Q, A0, tar, jump): - ''' + """ mcmcmh: Metrolopolis-Hasting MCMC algorithm generating N chains of length M for a parameter vector A of length n. @@ -252,7 +260,7 @@ def mcmcmh(M, N, M0, Q, A0, tar, jump): generated at the ith step of the jth chain was rejected so that AA(i,j,:) = AA(i-1,j,:), i > 1. The first set of proposal coincide with A0 are all accepted, so IAA(1,j) = 1. - ''' + """ A0 = np.array(A0) Q = np.array(Q) # number of parameters for which samples are to be drawn @@ -261,7 +269,6 @@ def mcmcmh(M, N, M0, Q, A0, tar, jump): # number of percentiles to be evaluated nQ = Q.size - # Initialising output arrays AA = np.empty((M, N, n)) IAA = np.zeros((M, N)) @@ -269,7 +276,6 @@ def mcmcmh(M, N, M0, Q, A0, tar, jump): Ne = np.empty((n)) S = np.empty((2 + nQ, n)) - # starting values of the sample and associated log of target density aq = A0 lq = tar(aq) @@ -281,7 +287,6 @@ def mcmcmh(M, N, M0, Q, A0, tar, jump): print("Initial values must be feasible for all chains") return None - # Run the chains in parallel for q in range(M): # draw from the jumping distribution and calculate @@ -308,34 +313,34 @@ def mcmcmh(M, N, M0, Q, A0, tar, jump): # Store Metropolis Hastings sample AA[q, :, :] = np.transpose(aq) - # acceptance probabilities for each chain - aP = 100*np.sum(IAA,0)/M + aP = 100 * np.sum(IAA, 0) / M # Convergence and summary statistics for each of the n parameters for j in range(n): # test convergence - RN = mcmcci(np.squeeze(AA[:,:,j]), M0) + RN = mcmcci(np.squeeze(AA[:, :, j]), M0) Rh[j] = RN[0] Ne[j] = RN[1] # provide summary information - asq = np.squeeze(AA[:,:,j]) - SS = mcsums(asq,M0,Q) - S[0,j] = SS[0] - S[1,j] = SS[1] - S[2:,j] = SS[2] + asq = np.squeeze(AA[:, :, j]) + SS = mcsums(asq, M0, Q) + S[0, j] = SS[0] + S[1, j] = SS[1] + S[2:, j] = SS[2] return S, aP, Rh, Ne, AA, IAA -class MCMCMH_NJ(): - ''' + +class MCMCMH_NJ: + """ Bayesian Noise and jitter reduction algorithm. MCMC used to determine the noise and jitter variances. Noise and jitter variances are then used in an iterative algorithm to remove the noise and jitter from the signal - ''' + """ def __init__(self, fs, ydata, N, niter, tol, m0w, s0w, m0t, s0t, Mc, M0, Nc, Q): - 'Setting initial variables' + "Setting initial variables" # variables for AnalyseSignalN and NJAlgorithm self.fs = fs @@ -357,10 +362,10 @@ def __init__(self, fs, ydata, N, niter, tol, m0w, s0w, m0t, s0t, Mc, M0, Nc, Q): # self.noiseSD = outs[1] def AnalyseSignalN(self): - ''' + """ Analyse signal to remove noise and jitter providing signal estimates with associated uncertainty. Uses normalised independent variable - ''' + """ # Initialisation self.N = np.int64(self.N) # Converting window length integer to int64 format @@ -387,8 +392,18 @@ def AnalyseSignalN(self): # Values of normalised independent variables datax = np.divide(np.arange(-n, n + 1), self.fs) - outs = self.mcmcm_main(self.ydata, datax, self.m0w, self.s0w, self.m0t, self.s0t, self.Mc, - self.M0, self.Nc, self.Q) + outs = self.mcmcm_main( + self.ydata, + datax, + self.m0w, + self.s0w, + self.m0t, + self.s0t, + self.Mc, + self.M0, + self.Nc, + self.Q, + ) self.jitterSD = outs[0] self.noiseSD = outs[1] # Loop through indices L of window @@ -400,14 +415,16 @@ def AnalyseSignalN(self): k = L + n # print(k) # Extract data in window - datay = self.ydata[L:L + self.N] + datay = self.ydata[L : L + self.N] # Inital polynomial approximation p = np.polyfit(datax, datay, 3) pval = np.polyval(p, datax) yhat0[k] = pval[n] # Applying algortithm to remove noise and jitter - [yhat[k], ck, vark, R[k]] = MCMCMH_NJ.NJAlgorithm(self, datax, datay, p, pval) + [yhat[k], ck, vark, R[k]] = MCMCMH_NJ.NJAlgorithm( + self, datax, datay, p, pval + ) print(yhat[k]) # First n windows, start building the covariance matrix Vmat for the data if L < n + 1: @@ -419,7 +436,6 @@ def AnalyseSignalN(self): Vmat[L - 1, L - 1] = vark cvec[L - n - 1, :] = ck - # For windows between 2n+1 and 4n, continue to build Vmat and cvec, and start building the sensitivtity # matrix Cmat from the sensitivtity vecotrs. Also, evaluate uncertainties for pervious estimates. elif L > np.multiply(2, n) and L < np.multiply(4, n) + 2: @@ -428,8 +444,10 @@ def AnalyseSignalN(self): # Count for building sensitivtity matrix iC = L - np.multiply(2, n) # Start building sensitivtity matrix from cvec - Cmat[iC - 1, :] = np.concatenate((np.zeros((1, iC - 1)), cvec[0, :], np.zeros((1, self.N - iC))), - axis=None) + Cmat[iC - 1, :] = np.concatenate( + (np.zeros((1, iC - 1)), cvec[0, :], np.zeros((1, self.N - iC))), + axis=None, + ) # Removing the first row of cvec and shift every row up one - creating # empty last row @@ -437,17 +455,16 @@ def AnalyseSignalN(self): cvec[-1, :] = 0 # Update empty last row cvec[n - 1, :] = ck - Cmatk = Cmat[0:iC, 1:self.N - 1 + iC] + Cmatk = Cmat[0:iC, 1 : self.N - 1 + iC] # Continue building Vmat - Vmatk = Vmat[1:self.N - 1 + iC, 1:self.N - 1 + iC] + Vmatk = Vmat[1 : self.N - 1 + iC, 1 : self.N - 1 + iC] V = np.matmul(np.matmul(Cmatk, Vmatk), np.transpose(Cmatk)) vhempty = np.empty((1, self.N - iC)) vhempty[:] = np.nan # Begin building vyhat vyhat[L, :] = np.concatenate((vhempty, V[iC - 1, :]), axis=None) - # For the remaining windows, update Vmat, Cmat and cvec. Continue to # evaluate the uncertainties for previous estimates. elif L > np.multiply(4, n) + 1: @@ -455,13 +472,18 @@ def AnalyseSignalN(self): Vmat = np.delete(Vmat, 0, axis=0) Vmat = np.delete(Vmat, 0, axis=1) Vmatzeros = np.zeros([Vmat.shape[0] + 1, Vmat.shape[1] + 1]) - Vmatzeros[:Vmat.shape[0], :Vmat.shape[1]] = Vmat + Vmatzeros[: Vmat.shape[0], : Vmat.shape[1]] = Vmat Vmat = Vmatzeros Vmat[2 * self.N - 2, 2 * self.N - 2] = vark # Building updated Cmat matrix - Cmat_old = np.concatenate((Cmat[1:self.N, 1:2 * self.N - 1], np.zeros([self.N - 1, 1])), axis=1) - Cmat_new = np.concatenate((np.zeros([1, self.N - 1]), cvec[0, :]), axis=None) + Cmat_old = np.concatenate( + (Cmat[1 : self.N, 1 : 2 * self.N - 1], np.zeros([self.N - 1, 1])), + axis=1, + ) + Cmat_new = np.concatenate( + (np.zeros([1, self.N - 1]), cvec[0, :]), axis=None + ) Cmat = np.concatenate((Cmat_old, Cmat_new[:, None].T), axis=0) # Update cvec @@ -474,15 +496,15 @@ def AnalyseSignalN(self): L += 1 - return (yhat[k]) + return yhat[k] def NJAlgorithm(self, datax, datay, p0, p0x): - ''' + """ Noise and Jitter Removal Algorithm- Iterative scheme that preprocesses data to reduce the effects of noise and jitter, resulting in an estimate of the true signal along with its associated uncertainty. Refer paper for details: https://ieeexplore.ieee.org/document/9138266 - ''' + """ # Initialisatio iter_ = 0 @@ -494,8 +516,14 @@ def NJAlgorithm(self, datax, datay, p0, p0x): k = np.int64(n + 1) t = np.array([np.power(datax[k], 3), np.power(datax[k], 2), datax[k], 1]) # Deisgn Matrix - X = np.array([np.power(datax, 3) + 3 * np.multiply(np.power(self.jitterSD, 2), datax), - np.power(datax, 2) + np.power(self.jitterSD, 2), datax, np.ones(np.size(datax))]) + X = np.array( + [ + np.power(datax, 3) + 3 * np.multiply(np.power(self.jitterSD, 2), datax), + np.power(datax, 2) + np.power(self.jitterSD, 2), + datax, + np.ones(np.size(datax)), + ] + ) X = X.T # Iterative algortithm @@ -509,8 +537,15 @@ def NJAlgorithm(self, datax, datay, p0, p0x): # Step 4 # Weight calculation - w = np.divide(1, [np.sqrt(np.power(self.jitterSD, 2) * np.power(pdx, 2) - + np.power(self.noiseSD, 2))]) + w = np.divide( + 1, + [ + np.sqrt( + np.power(self.jitterSD, 2) * np.power(pdx, 2) + + np.power(self.noiseSD, 2) + ) + ], + ) w = w.T # Calculating polynomial coeffs @@ -526,7 +561,7 @@ def NJAlgorithm(self, datax, datay, p0, p0x): p0x = p1x if iter_ == self.niter: - print('Maximum number of iterations reached') + print("Maximum number of iterations reached") break # Evaluate outputs @@ -535,7 +570,9 @@ def NJAlgorithm(self, datax, datay, p0, p0x): pd = np.polyder(p0) pdx = np.polyval(pd, datax[k]) vk = np.power(self.jitterSD, 2) * np.power(pdx, 2) + np.power(self.noiseSD, 2) - R = np.power(np.linalg.norm(np.matmul(np.diagflat(w), (datay - np.matmul(X, p0)))), 2) + R = np.power( + np.linalg.norm(np.matmul(np.diagflat(w), (datay - np.matmul(X, p0)))), 2 + ) return yhat, c, vk, R @@ -567,7 +604,7 @@ def mcmcm_main(datay, datax, m0w, s0w, m0t, s0t, Mc, M0, Nc, Q): @staticmethod def tar_at(at, y, x, m0w, s0w, m0t, s0t): - ''' + """ ------------------------------------------------------------------------- Target dist for noise and jitter posterior dist ------------------------------------------------------------------------- @@ -587,7 +624,7 @@ def tar_at(at, y, x, m0w, s0w, m0t, s0t): Output: T: Log of the posterior distribution - ''' + """ # Size of parameter vector at = np.array(at) @@ -624,9 +661,10 @@ def tar_at(at, y, x, m0w, s0w, m0t, s0t): return T """### jumprwg: Jumping distribution for the Metropolis Hastings Gaussian random walk algorithm""" + @staticmethod def jumprwg(A, L): - ''' + """ jumprwg: Jumping distribution for the Metropolis Hastings Gaussian random walk algorithm ------------------------------------------------------------------------- @@ -647,7 +685,7 @@ def jumprwg(A, L): moving from As(:,j) to A(:,j), up to an additive constant. log P0(a|as) - log P0(as|a) - ''' + """ # Number of parameters and parallel chains n, N = A.shape @@ -662,9 +700,25 @@ def jumprwg(A, L): return As, dp0 + ######################################## class NoiseJitterRemovalAgent(AgentMET4FOF): - def init_parameters(self, fs=100, ydata = np.array([]), N=15, niter=100, tol=1e-9, m0w = 10, s0w = 0.0005, m0t = 10, s0t = 0.0002*100/8, Mc=5000, M0=100, Nc=100, Q=50 ): + def init_parameters( + self, + fs=100, + ydata=np.array([]), + N=15, + niter=100, + tol=1e-9, + m0w=10, + s0w=0.0005, + m0t=10, + s0t=0.0002 * 100 / 8, + Mc=5000, + M0=100, + Nc=100, + Q=50, + ): self.fs = fs self.ydata = ydata self.N = N @@ -681,7 +735,9 @@ def init_parameters(self, fs=100, ydata = np.array([]), N=15, niter=100, tol=1e @staticmethod def njr(fs, ydata, N, niter, tol, m0w, s0w, m0t, s0t, Mc, M0, Nc, Q) -> object: - analyse_fun = MCMCMH_NJ(fs, ydata, N, niter, tol, m0w, s0w, m0t, s0t, Mc, M0, Nc, Q) + analyse_fun = MCMCMH_NJ( + fs, ydata, N, niter, tol, m0w, s0w, m0t, s0t, Mc, M0, Nc, Q + ) yhat1 = analyse_fun.AnalyseSignalN() return yhat1 @@ -692,6 +748,20 @@ def on_received_message(self, message): ddata = message["data"] self.ydata = np.append(self.ydata, ddata) if np.size(self.ydata) == self.N: - t = self.njr(self.fs, self.ydata, self.N, self.niter, self.tol, self.m0w, self.s0w, self.m0t, self.s0t, self.Mc, self.M0,self.Nc, self.Q) + t = self.njr( + self.fs, + self.ydata, + self.N, + self.niter, + self.tol, + self.m0w, + self.s0w, + self.m0t, + self.s0t, + self.Mc, + self.M0, + self.Nc, + self.Q, + ) self.send_output(self.ydata[7] - t) - self.ydata = self.ydata[1:self.N] + self.ydata = self.ydata[1 : self.N] diff --git a/agentMET4FOF_tutorials/noise_jitter/remove_noise_and_jitter.py b/agentMET4FOF_tutorials/noise_jitter/remove_noise_and_jitter.py index 6ab7ab36..8f15f594 100644 --- a/agentMET4FOF_tutorials/noise_jitter/remove_noise_and_jitter.py +++ b/agentMET4FOF_tutorials/noise_jitter/remove_noise_and_jitter.py @@ -1,6 +1,9 @@ from agentMET4FOF.agents.base_agents import MonitorAgent from agentMET4FOF.agents.noise_jitter_removal_agents import NoiseJitterRemovalAgent -from agentMET4FOF.agents.signal_agents import StaticSineWithJitterGeneratorAgent, NoiseAgent +from agentMET4FOF.agents.signal_agents import ( + StaticSineWithJitterGeneratorAgent, + NoiseAgent, +) from agentMET4FOF.network import AgentNetwork @@ -9,14 +12,20 @@ def demonstrate_noise_jitter_removal_agent(): agentNetwork = AgentNetwork(backend="mesa") # init agents - sine_with_jitter_agent = agentNetwork.add_agent(agentType=StaticSineWithJitterGeneratorAgent) + sine_with_jitter_agent = agentNetwork.add_agent( + agentType=StaticSineWithJitterGeneratorAgent + ) noise_agent = agentNetwork.add_agent(agentType=NoiseAgent) njremove_agent = agentNetwork.add_agent(agentType=NoiseJitterRemovalAgent) - monitor_agent = agentNetwork.add_agent(agentType=MonitorAgent, name="Sine with Noise and Jitter") - monitor_agent2 = agentNetwork.add_agent(agentType=MonitorAgent, name="Output of Noise-Jitter Removal Agent") + monitor_agent = agentNetwork.add_agent( + agentType=MonitorAgent, name="Sine with Noise and Jitter" + ) + monitor_agent2 = agentNetwork.add_agent( + agentType=MonitorAgent, name="Output of Noise-Jitter Removal Agent" + ) # connect agents : jitter generator -> noise -> njremoval agent agentNetwork.bind_agents(sine_with_jitter_agent, noise_agent) @@ -33,5 +42,5 @@ def demonstrate_noise_jitter_removal_agent(): return agentNetwork -if __name__ == '__main__': +if __name__ == "__main__": demonstrate_noise_jitter_removal_agent() From 045866df0d42b09c2a825d932ef2d4ca35c2440b Mon Sep 17 00:00:00 2001 From: Bjoern Ludwig Date: Fri, 30 Jul 2021 14:35:27 +0200 Subject: [PATCH 16/43] wip(remove_noise_and_jitter): optimize imports --- agentMET4FOF_tutorials/noise_jitter/remove_noise_and_jitter.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/agentMET4FOF_tutorials/noise_jitter/remove_noise_and_jitter.py b/agentMET4FOF_tutorials/noise_jitter/remove_noise_and_jitter.py index 8f15f594..394f1fa2 100644 --- a/agentMET4FOF_tutorials/noise_jitter/remove_noise_and_jitter.py +++ b/agentMET4FOF_tutorials/noise_jitter/remove_noise_and_jitter.py @@ -1,8 +1,8 @@ from agentMET4FOF.agents.base_agents import MonitorAgent from agentMET4FOF.agents.noise_jitter_removal_agents import NoiseJitterRemovalAgent from agentMET4FOF.agents.signal_agents import ( - StaticSineWithJitterGeneratorAgent, NoiseAgent, + StaticSineWithJitterGeneratorAgent, ) from agentMET4FOF.network import AgentNetwork From ef407fc1a6ecb23318e9e20411574c92f8542b5e Mon Sep 17 00:00:00 2001 From: anupam-prasad Date: Fri, 30 Jul 2021 15:00:51 +0200 Subject: [PATCH 17/43] wip(noise_jitter_removal_agents): incorporates all methods into MCMCMH_JN class --- .../agents/noise_jitter_removal_agents.py | 654 +++++++++--------- 1 file changed, 318 insertions(+), 336 deletions(-) diff --git a/agentMET4FOF/agents/noise_jitter_removal_agents.py b/agentMET4FOF/agents/noise_jitter_removal_agents.py index b2509b75..8cbebd51 100644 --- a/agentMET4FOF/agents/noise_jitter_removal_agents.py +++ b/agentMET4FOF/agents/noise_jitter_removal_agents.py @@ -5,333 +5,6 @@ from agentMET4FOF.agents import AgentMET4FOF -"""### mcmcci: MCMC convergence indices for multiple chains.""" - - -def mcmcci(A, M0): - """ - mcmcci: MCMC convergence indices for multiple chains. - ---------------------------------------------------------------------------- - KJ, LRW, PMH - - Version 2020-04-22 - - ---------------------------------------------------------------------------- - Inputs: - A(M, N): Chain samples, N chains of length M. - - M0: Length of burn-in, M > M0 >= 0. - - Outputs: - Rhat: Convergence index. Rhat is expected to be greater than 1. The closer - Rhat is to 1, the better the convergence. - - Neff: Estimate of the effective number of independent draws. Neff is - expected to be less than (M-M0)*N. - ---------------------------------------------------------------------------- - Note: If the calculated value of Rhat is < 1, then Rhat is set to 1 and Neff - set to (M-M0)*N, their limit values. - - Note: If N = 1 or M0 > M-2, Rhat = 0; Neff = 0. - """ - A = np.array(A) - - M, N = A.shape - - # Initialisation - Rhat = 0 - Neff = 0 - - # The convergence statistics can only be evaluated if there are multiple chains - # and the chain length is greater than the burn in length - - if N > 1 and M > M0 + 1: - Mt = M - M0 - - # Chain mean and mean of means - asub = A[M0:, :] - ad = np.mean(asub, axis=0) - add = asub.mean() - - # Within group standard deviation - ss = np.std(asub, axis=0) - - # Between groups variance. - dd = np.square(ad - add) - B = (Mt * np.sum(dd)) / (N - 1) - - # Within groups variance. - W = np.sum(np.square(ss)) / N - - # V plus - Vp = (Mt - 1) * W / Mt + B / Mt - - # Convergence statistic, effective number of independent samples - Rhat = np.sqrt(Vp / W) - Neff = Mt * N * Vp / B - - Rhat = np.maximum(Rhat, 1) - Neff = np.minimum(Neff, Mt * N) - - return Rhat, Neff - - -"""### mcsums: Summary information from MC samples.""" - - -def mcsums(A, M0, Q): - """ - mcsums: Summary information from MC samples. - ------------------------------------------------------------------------- - KJ, LRW, PMH - Version 2020-04-22 - - ------------------------------------------------------------------------- - Inputs: - A(M,N): An array that stores samples of size M x N. - - M0: Burn-in period with M > M0 >= 0. - - Q(nQ,1): Percentiles specifications, 0 <= Q(l) <= 100. - - Outputs: - abar(n,1): Mean for each sample. - - s(n,1): Standard deviation for sample. - - aQ(nQ,n): Percentiles corresponding to Q. - - """ - # Size of samples after burn-in - A = np.array(A) - - M, N = A.shape - - m = (M - M0) * N - - # Initialise percentile vector - nQ = Q.size - aQ = np.zeros(nQ) - - # Samples from N parallel chains after burn-in period - aaj = A[M0:, :] - aaj = aaj.flatten() - - # Mean and standard deviation of samples - abar = np.mean(aaj) - s = np.std(aaj) - - # Percentiles of samples - aQ = np.percentile(aaj, Q) - - return abar, s, aQ - - -"""### Cubic function and its first and second derivative""" - - -def fgh_cubic(alpha, t): - """ - ------------------------------------------------------------------------- - Cubic function and its first and second derivative - ------------------------------------------------------------------------- - KJ, LRW, PMH - Version 2020-04-22 - ------------------------------------------------------------------------- - Inputs: - alpha(4,N): Alpha parameters - - t(m,1): Times - - Outputs: - f(m,N): Cubic function - - f1(m,N): Derivative of cubic - - f2(m,N): Second derivative of cubic - """ - - # length of data and number of paramaters - m = t.size - - # design matrix - C = np.array([np.ones(m), t, t ** 2, t ** 3]) - - # derivate info - C1 = np.array([np.ones(m), 2 * t, 3 * t ** 2]) - C2 = np.array([2 * np.ones(m), 6 * t]) - - # cubic and derivatives - f = np.matmul(C.T, alpha) - f1 = np.matmul(C1.T, alpha[1:]) - f2 = np.matmul(C2.T, alpha[2:]) - - return f, f1, f2 - - -"""### Log of the gaussian pdf""" - - -def ln_gauss_pdf_v(x, mu, sigma): - """ - ------------------------------------------------------------------------- - Log of the Gaussian pdf - ------------------------------------------------------------------------- - KJ, LRW, PMH - Version 2020-03-12 - -------------------------------------------------------------------------- - Inputs: - x(m,1): Points at which pdf is to be evaluated - - mu: Mean of distribution - - sigma: Standard deviation of the distribution - - Output: - logf: Log of the Gaussian pdf at x with mean mu and - std sigma - """ - - try: - # When inputs are high dimensional arrays/matrices - xx = np.matlib.repmat(x, mu.shape[1], 1) - xx = xx.T - logk = -np.log(2 * math.pi) / 2 - np.log(sigma) - logp = -((xx - mu) ** 2) / (2 * sigma ** 2) - # Log of the Gaussian PDF - logf = logk + logp - - except IndexError: - # When inputs are vectors - logk = -np.log(2 * math.pi) / 2 - np.log(sigma) - logp = -((x - mu) ** 2) / (2 * sigma ** 2) - # Log of the Gaussian PDF - logf = logk + logp - - return logf - - -"""###mcmcmh: Metrolopolis-Hasting MCMC algorithm generating N chains of length M for a parameter vector A of length n.""" - - -def mcmcmh(M, N, M0, Q, A0, tar, jump): - """ - mcmcmh: Metrolopolis-Hasting MCMC algorithm generating N chains of length - M for a parameter vector A of length n. - - For details about the algorithm please refer to: - Gelman A, Carlin JB, Stern HS, Dunson DB, Vehtari A, Rubin DB. - Bayesian data analysis. CRC press; 2013 Nov 1. - ------------------------------------------------------------------------- - KJ, LRW, PMH - Version 2020-04-22 - ------------------------------------------------------------------------- - Inputs: - M: Length of the chains. - - N: Number of chains. - - M0: Burn in period. - - Q(nQ,1): Percentiles 0 <= Q(k) <= 100. - - A0(n,N): Array of feasible starting points: the target distribution - evaluated at A0(:,j) is strictly positive. - - Outputs: - - S(2+nQ,n): Summary of A - mean, standard deviation and percentile limits, - where the percentile limits are given by Q. - - aP(N,1): Acceptance percentages for AA calculated for each chain. - - Rh(n,1): Estimate of convergence. Theoretically Rh >= 1, and the closer - to 1, the more evidence of convergence. - - Ne(n,1): Estimate of the number of effective number of independent draws. - - AA(M,N,n): Array storing the chains: A(i,j,k) is the kth element of the - parameter vector stored as the ith member of the jth chain. - AA(1,j,:) = A0(:,j). - - IAA(M,N): Acceptance indices. IAA(i,j) = 1 means that the proposal - as(n,1) generated at the ith step of the jth chain was accepted so - that AA(i,j,:) = as. IAA(i,j) = 0 means that the proposal as(n,1) - generated at the ith step of the jth chain was rejected so that - AA(i,j,:) = AA(i-1,j,:), i > 1. The first set of proposal coincide with - A0 are all accepted, so IAA(1,j) = 1. - """ - A0 = np.array(A0) - Q = np.array(Q) - # number of parameters for which samples are to be drawn - n = A0.shape[0] - - # number of percentiles to be evaluated - nQ = Q.size - - # Initialising output arrays - AA = np.empty((M, N, n)) - IAA = np.zeros((M, N)) - Rh = np.empty((n)) - Ne = np.empty((n)) - S = np.empty((2 + nQ, n)) - - # starting values of the sample and associated log of target density - aq = A0 - lq = tar(aq) - - # Starting values must be feasible for each chain - Id = lq > -np.Inf - - if sum(Id) < N: - print("Initial values must be feasible for all chains") - return None - - # Run the chains in parallel - for q in range(M): - # draw from the jumping distribution and calculate - # d = log P0(aq|as) - log P0(as|aq) - - asam, d = jump(aq) - - # log of the target density for the new draw as - ls = tar(asam) - - # Metropolis-Hastings acceptance ratio - rq = np.exp(ls - lq + d) - - # draws from the uniform distribution for acceptance rule - uq = np.random.rand(N) - # index of samples that have been accepted - ind = uq < rq - - # updating the sample and evaluating acceptance indices - aq[:, ind] = asam[:, ind] - lq[ind] = ls[ind] - IAA[q, ind] = 1 - - # Store Metropolis Hastings sample - AA[q, :, :] = np.transpose(aq) - - # acceptance probabilities for each chain - aP = 100 * np.sum(IAA, 0) / M - - # Convergence and summary statistics for each of the n parameters - for j in range(n): - # test convergence - RN = mcmcci(np.squeeze(AA[:, :, j]), M0) - Rh[j] = RN[0] - Ne[j] = RN[1] - - # provide summary information - asq = np.squeeze(AA[:, :, j]) - SS = mcsums(asq, M0, Q) - S[0, j] = SS[0] - S[1, j] = SS[1] - S[2:, j] = SS[2] - - return S, aP, Rh, Ne, AA, IAA - class MCMCMH_NJ: """ @@ -415,7 +88,7 @@ def AnalyseSignalN(self): k = L + n # print(k) # Extract data in window - datay = self.ydata[L : L + self.N] + datay = self.ydata[L: L + self.N] # Inital polynomial approximation p = np.polyfit(datax, datay, 3) pval = np.polyval(p, datax) @@ -432,13 +105,13 @@ def AnalyseSignalN(self): # for windows n+1 to 2n, continue building the covariance matrix Vmat and start stroing the # sensitivtity vectors ck in cvec - elif L > n and L < np.multiply(2, n) + 1: + elif n < L < np.multiply(2, n) + 1: Vmat[L - 1, L - 1] = vark cvec[L - n - 1, :] = ck # For windows between 2n+1 and 4n, continue to build Vmat and cvec, and start building the sensitivtity # matrix Cmat from the sensitivtity vecotrs. Also, evaluate uncertainties for pervious estimates. - elif L > np.multiply(2, n) and L < np.multiply(4, n) + 2: + elif np.multiply(2, n) < L < np.multiply(4, n) + 2: Vmat[L - 1, L - 1] = vark # Count for building sensitivtity matrix @@ -597,11 +270,9 @@ def mcmcm_main(datay, datax, m0w, s0w, m0t, s0t, Mc, M0, Nc, Q): A0 = np.matlib.repmat(pars.T, Nc, 1).T + np.matmul(L, rr) - sam = mcmcmh(Mc, Nc, M0, Q, A0, tar, jump) + sam = MCMCMH_NJ.mcmcmh(Mc, Nc, M0, Q, A0, tar, jump) return 1 / np.sqrt(np.exp(sam[0][0, -2:])) - """### Target dist for noise and jitter posterior dist""" - @staticmethod def tar_at(at, y, x, m0w, s0w, m0t, s0t): """ @@ -645,7 +316,7 @@ def tar_at(at, y, x, m0w, s0w, m0t, s0t): prior_phi2 = (m0w / 2) * np.log(phi2) - phi2 * m0w * s0w ** 2 / 2 # function that evaluates the cubic function with user specified cubic parameters - fun = lambda aa: fgh_cubic(aa, x) + fun = lambda aa: MCMCMH_NJ.fgh_cubic(aa, x) # cubic, expectation and variance [st, st1, st2] = fun(alpha) @@ -653,14 +324,90 @@ def tar_at(at, y, x, m0w, s0w, m0t, s0t): vari = (taus ** 2) * (st1 ** 2) + omegas ** 2 # Likelihood - lik = sum(ln_gauss_pdf_v(y, expect, np.sqrt(vari))) + lik = sum(MCMCMH_NJ.ln_gauss_pdf_v(y, expect, np.sqrt(vari))) # Posterior T = lik + prior_phi1 + prior_phi2 return T - """### jumprwg: Jumping distribution for the Metropolis Hastings Gaussian random walk algorithm""" + @staticmethod + def fgh_cubic(alpha, t): + """ + ------------------------------------------------------------------------- + Cubic function and its first and second derivative + ------------------------------------------------------------------------- + KJ, LRW, PMH + Version 2020-04-22 + ------------------------------------------------------------------------- + Inputs: + alpha(4,N): Alpha parameters + + t(m,1): Times + + Outputs: + f(m,N): Cubic function + + f1(m,N): Derivative of cubic + + f2(m,N): Second derivative of cubic + """ + + # length of data and number of paramaters + m = t.size + + # design matrix + C = np.array([np.ones(m), t, t ** 2, t ** 3]) + + # derivate info + C1 = np.array([np.ones(m), 2 * t, 3 * t ** 2]) + C2 = np.array([2 * np.ones(m), 6 * t]) + + # cubic and derivatives + f = np.matmul(C.T, alpha) + f1 = np.matmul(C1.T, alpha[1:]) + f2 = np.matmul(C2.T, alpha[2:]) + + return f, f1, f2 + + @staticmethod + def ln_gauss_pdf_v(x, mu, sigma): + """ + ------------------------------------------------------------------------- + Log of the Gaussian pdf + ------------------------------------------------------------------------- + KJ, LRW, PMH + Version 2020-03-12 + -------------------------------------------------------------------------- + Inputs: + x(m,1): Points at which pdf is to be evaluated + + mu: Mean of distribution + + sigma: Standard deviation of the distribution + + Output: + logf: Log of the Gaussian pdf at x with mean mu and + std sigma + """ + + try: + # When inputs are high dimensional arrays/matrices + xx = np.matlib.repmat(x, mu.shape[1], 1) + xx = xx.T + logk = -np.log(2 * math.pi) / 2 - np.log(sigma) + logp = -((xx - mu) ** 2) / (2 * sigma ** 2) + # Log of the Gaussian PDF + logf = logk + logp + + except IndexError: + # When inputs are vectors + logk = -np.log(2 * math.pi) / 2 - np.log(sigma) + logp = -((x - mu) ** 2) / (2 * sigma ** 2) + # Log of the Gaussian PDF + logf = logk + logp + + return logf @staticmethod def jumprwg(A, L): @@ -700,6 +447,241 @@ def jumprwg(A, L): return As, dp0 + @staticmethod + def mcmcmh(M, N, M0, Q, A0, tar, jump): + """ + mcmcmh: Metrolopolis-Hasting MCMC algorithm generating N chains of length + M for a parameter vector A of length n. + + For details about the algorithm please refer to: + Gelman A, Carlin JB, Stern HS, Dunson DB, Vehtari A, Rubin DB. + Bayesian data analysis. CRC press; 2013 Nov 1. + ------------------------------------------------------------------------- + KJ, LRW, PMH + Version 2020-04-22 + ------------------------------------------------------------------------- + Inputs: + M: Length of the chains. + + N: Number of chains. + + M0: Burn in period. + + Q(nQ,1): Percentiles 0 <= Q(k) <= 100. + + A0(n,N): Array of feasible starting points: the target distribution + evaluated at A0(:,j) is strictly positive. + + Outputs: + + S(2+nQ,n): Summary of A - mean, standard deviation and percentile limits, + where the percentile limits are given by Q. + + aP(N,1): Acceptance percentages for AA calculated for each chain. + + Rh(n,1): Estimate of convergence. Theoretically Rh >= 1, and the closer + to 1, the more evidence of convergence. + + Ne(n,1): Estimate of the number of effective number of independent draws. + + AA(M,N,n): Array storing the chains: A(i,j,k) is the kth element of the + parameter vector stored as the ith member of the jth chain. + AA(1,j,:) = A0(:,j). + + IAA(M,N): Acceptance indices. IAA(i,j) = 1 means that the proposal + as(n,1) generated at the ith step of the jth chain was accepted so + that AA(i,j,:) = as. IAA(i,j) = 0 means that the proposal as(n,1) + generated at the ith step of the jth chain was rejected so that + AA(i,j,:) = AA(i-1,j,:), i > 1. The first set of proposal coincide with + A0 are all accepted, so IAA(1,j) = 1. + """ + A0 = np.array(A0) + Q = np.array(Q) + # number of parameters for which samples are to be drawn + n = A0.shape[0] + + # number of percentiles to be evaluated + nQ = Q.size + + # Initialising output arrays + AA = np.empty((M, N, n)) + IAA = np.zeros((M, N)) + Rh = np.empty((n)) + Ne = np.empty((n)) + S = np.empty((2 + nQ, n)) + + # starting values of the sample and associated log of target density + aq = A0 + lq = tar(aq) + + # Starting values must be feasible for each chain + Id = lq > -np.Inf + + if sum(Id) < N: + print("Initial values must be feasible for all chains") + return None + + # Run the chains in parallel + for q in range(M): + # draw from the jumping distribution and calculate + # d = log P0(aq|as) - log P0(as|aq) + + asam, d = jump(aq) + + # log of the target density for the new draw as + ls = tar(asam) + + # Metropolis-Hastings acceptance ratio + rq = np.exp(ls - lq + d) + + # draws from the uniform distribution for acceptance rule + uq = np.random.rand(N) + # index of samples that have been accepted + ind = uq < rq + + # updating the sample and evaluating acceptance indices + aq[:, ind] = asam[:, ind] + lq[ind] = ls[ind] + IAA[q, ind] = 1 + + # Store Metropolis Hastings sample + AA[q, :, :] = np.transpose(aq) + + # acceptance probabilities for each chain + aP = 100 * np.sum(IAA, 0) / M + + # Convergence and summary statistics for each of the n parameters + for j in range(n): + # test convergence + RN = MCMCMH_NJ.mcmcci(np.squeeze(AA[:, :, j]), M0) + Rh[j] = RN[0] + Ne[j] = RN[1] + + # provide summary information + asq = np.squeeze(AA[:, :, j]) + SS = MCMCMH_NJ.mcsums(asq, M0, Q) + S[0, j] = SS[0] + S[1, j] = SS[1] + S[2:, j] = SS[2] + + return S, aP, Rh, Ne, AA, IAA + + @staticmethod + def mcmcci(A, M0): + """ + mcmcci: MCMC convergence indices for multiple chains. + ---------------------------------------------------------------------------- + KJ, LRW, PMH + + Version 2020-04-22 + + ---------------------------------------------------------------------------- + Inputs: + A(M, N): Chain samples, N chains of length M. + + M0: Length of burn-in, M > M0 >= 0. + + Outputs: + Rhat: Convergence index. Rhat is expected to be greater than 1. The closer + Rhat is to 1, the better the convergence. + + Neff: Estimate of the effective number of independent draws. Neff is + expected to be less than (M-M0)*N. + ---------------------------------------------------------------------------- + Note: If the calculated value of Rhat is < 1, then Rhat is set to 1 and Neff + set to (M-M0)*N, their limit values. + + Note: If N = 1 or M0 > M-2, Rhat = 0; Neff = 0. + """ + A = np.array(A) + + M, N = A.shape + + # Initialisation + Rhat = 0 + Neff = 0 + + # The convergence statistics can only be evaluated if there are multiple chains + # and the chain length is greater than the burn in length + + if N > 1 and M > M0 + 1: + Mt = M - M0 + + # Chain mean and mean of means + asub = A[M0:, :] + ad = np.mean(asub, axis=0) + add = asub.mean() + + # Within group standard deviation + ss = np.std(asub, axis=0) + + # Between groups variance. + dd = np.square(ad - add) + B = (Mt * np.sum(dd)) / (N - 1) + + # Within groups variance. + W = np.sum(np.square(ss)) / N + + # V plus + Vp = (Mt - 1) * W / Mt + B / Mt + + # Convergence statistic, effective number of independent samples + Rhat = np.sqrt(Vp / W) + Neff = Mt * N * Vp / B + + Rhat = np.maximum(Rhat, 1) + Neff = np.minimum(Neff, Mt * N) + + return Rhat, Neff + + @staticmethod + def mcsums(A, M0, Q): + """ + mcsums: Summary information from MC samples. + ------------------------------------------------------------------------- + KJ, LRW, PMH + Version 2020-04-22 + + ------------------------------------------------------------------------- + Inputs: + A(M,N): An array that stores samples of size M x N. + + M0: Burn-in period with M > M0 >= 0. + + Q(nQ,1): Percentiles specifications, 0 <= Q(l) <= 100. + + Outputs: + abar(n,1): Mean for each sample. + + s(n,1): Standard deviation for sample. + + aQ(nQ,n): Percentiles corresponding to Q. + + """ + # Size of samples after burn-in + A = np.array(A) + + M, N = A.shape + + m = (M - M0) * N + + # Initialise percentile vector + nQ = Q.size + aQ = np.zeros(nQ) + + # Samples from N parallel chains after burn-in period + aaj = A[M0:, :] + aaj = aaj.flatten() + + # Mean and standard deviation of samples + abar = np.mean(aaj) + s = np.std(aaj) + + # Percentiles of samples + aQ = np.percentile(aaj, Q) + + return abar, s, aQ + ######################################## class NoiseJitterRemovalAgent(AgentMET4FOF): From b6ca636fc7eeae693819bd82893328d3e838fdd0 Mon Sep 17 00:00:00 2001 From: Bjoern Ludwig Date: Fri, 30 Jul 2021 15:07:34 +0200 Subject: [PATCH 18/43] feat(SineWithJitter): introduce a streaming data stream of a sine signal with jitter and the corresponding agent --- agentMET4FOF/agents/signal_agents.py | 61 +++++++++++++- agentMET4FOF/streams/signal_streams.py | 80 ++++++++++++++++++- .../noise_jitter/generate_sine_with_jitter.py | 10 ++- 3 files changed, 147 insertions(+), 4 deletions(-) diff --git a/agentMET4FOF/agents/signal_agents.py b/agentMET4FOF/agents/signal_agents.py index 060b5f64..809e263c 100644 --- a/agentMET4FOF/agents/signal_agents.py +++ b/agentMET4FOF/agents/signal_agents.py @@ -4,7 +4,11 @@ import pandas as pd from .base_agents import AgentMET4FOF -from ..streams.signal_streams import SineGenerator, StaticSineWithJitterGenerator +from ..streams.signal_streams import ( + SineGenerator, + SineWithJitterGenerator, + StaticSineWithJitterGenerator, +) __all__ = ["SineGeneratorAgent", "StaticSineWithJitterGeneratorAgent", "NoiseAgent"] @@ -92,6 +96,7 @@ def agent_loop(self): class NoiseAgent(AgentMET4FOF): """An agent adding white noise to the incoming signal""" + _noise_std: float @property @@ -155,3 +160,57 @@ def _compute_noisy_signal_from_clean_signal( data_in_message["quantities"] ) self.send_output(fully_assembled_resulting_data) + + +class SineWithJitterGeneratorAgent(SineGeneratorAgent): + """An agent streaming a sine signal + + Takes samples from the :py:mod:`SineWithJitterGenerator` and pushes them sample by + sample to connected agents via its output channel. + """ + + def init_parameters( + self, + sfreq: Optional[int] = 10, + sine_freq: Optional[float] = np.reciprocal(2 * np.pi), + amplitude: Optional[float] = 1.0, + initial_phase: Optional[float] = 0.0, + jitter_std: Optional[float] = 0.02, + ): + """Initialize the input data + + Initialize the input data stream as an instance of the + :class:`SineWithJitterGenerator` class. + + Parameters + ---------- + sfreq : int, optional + sampling frequency which determines the time step when :meth:`.next_sample` + is called, defaults to 10 + sine_freq : float, optional + frequency of the generated sine wave, defaults to :math:`\frac{1}{2 \pi}` + amplitude : float, optional + amplitude of the generated sine wave, defaults to 1.0 + initial_phase : float, optional + initial phase (at t=0) of the generated sine wave, defaults to 0.0 + jitter_std : float, optional + the standard deviation of the distribution to randomly draw jitter from, + defaults to 0.02 + """ + self._sine_stream = SineWithJitterGenerator( + sfreq=sfreq, + sine_freq=sine_freq, + amplitude=amplitude, + initial_phase=initial_phase, + jitter_std=jitter_std, + ) + + def agent_loop(self): + """Model the agent's behaviour + + On state *Running* the agent will extract sample by sample the input data + streams content and push it via invoking :meth:`AgentMET4FOF.send_output`. + """ + if self.current_state == "Running": + sine_data = self._sine_stream.next_sample() + self.send_output(sine_data) diff --git a/agentMET4FOF/streams/signal_streams.py b/agentMET4FOF/streams/signal_streams.py index 20192f8b..5caeb348 100644 --- a/agentMET4FOF/streams/signal_streams.py +++ b/agentMET4FOF/streams/signal_streams.py @@ -1,10 +1,15 @@ -from typing import Optional +from typing import Dict, Optional import numpy as np from .base_streams import DataStreamMET4FOF -__all__ = ["SineGenerator", "CosineGenerator", "StaticSineWithJitterGenerator"] +__all__ = [ + "SineGenerator", + "CosineGenerator", + "SineWithJitterGenerator", + "StaticSineWithJitterGenerator", +] class SineGenerator(DataStreamMET4FOF): @@ -151,3 +156,74 @@ def __init__(self, num_cycles=1000, jitter_std=0.02): timestamps_with_jitter = np.random.normal(loc=timestamps, scale=jitter_std) signal_values_at_timestamps = np.sin(timestamps_with_jitter) self.set_data_source(quantities=signal_values_at_timestamps, time=timestamps) + + +class SineWithJitterGenerator(SineGenerator): + r"""Represents a streamed sine signal with jitter + + Parameters + ---------- + sfreq : int, optional + sampling frequency which determines the time step when :meth:`.next_sample` + is called, defaults to 10 + sine_freq : float, optional + frequency of wave function, defaults to :math:`\frac{1}{2 \pi}` + amplitude : float, optional + amplitude of the wave function, defaults to 1.0 + initial_phase : float, optional + initial phase of the wave function, defaults to 0.0 + jitter_std : float, optional + the standard deviation of the distribution to randomly draw jitter from, + defaults to 0.02 + """ + + _jitter_std: float + + @property + def jitter_std(self): + """The standard deviation of the distribution to randomly draw jitter from""" + return self._jitter_std + + def __init__( + self, + sfreq: Optional[int] = 10, + sine_freq: Optional[float] = np.reciprocal(2 * np.pi), + amplitude: Optional[float] = 1.0, + initial_phase: Optional[float] = 0.0, + jitter_std: Optional[float] = 0.02, + ): + self._jitter_std = jitter_std + super().__init__( + sfreq=sfreq, + sine_freq=sine_freq, + amplitude=amplitude, + initial_phase=initial_phase, + ) + + def _next_sample_generator( + self, batch_size: Optional[int] = 1 + ) -> Dict[str, np.ndarray]: + """Generate the next batch of samples from the sine function with jitter + + Parameters + ---------- + batch_size : int, optional + number of batches to get from data stream, defaults to 1 + + Returns + ------- + Dict[str, np.ndarray] + latest samples of the sine signal with jitter in the form:: + + dict like { + "quantities":