PK!S derivative/__init__.pyfrom .loc import * from .glob import * methods = { "finitediff": finitediff, "holoborodko": holoborodko, "fft_deriv": fft_deriv, "cubic_spline": cubic_spline_deriv, "tvregdiff": tvregdiff, } def derivative(x, y, kind="finitediff", **kwargs): method = methods.get(kind) if len(y.shape) > 1: return np.apply_along_axis(method, x, y, **kwargs) else: return method(x, y, **kwargs) PK!derivative/__version__.py__version__ = "0.1.0" PK!L33derivative/glob.pyimport numpy as np from scipy import sparse from scipy.sparse import linalg as splin from scipy.interpolate import CubicSpline def cubic_spline_deriv(x, y, order=1): cs = CubicSpline(x, y) return cs.derivative(order)(x) def fft_deriv(x, y): n = len(x) dx = x[1] - x[0] dw = 2 * np.pi / (x[-1] - x[0]) w = np.fft.fftfreq(n) * n * dw return np.fft.ifft(1j * w * np.fft.fft(y)) def tvregdiff(x, y, **kwargs): return _tvregdiff(y.reshape(-1), dx=x[1] - x[0], **kwargs)[:-1] def _tvregdiff(data, itern=25, alph=0.05, u0=None, scale="small", ep=1e-6, dx=None, diagflag=0): # u = TVRegDiff( data, iter, alph, u0, scale, ep, dx, plotflag, diagflag ); # # Rick Chartrand (rickc@lanl.gov), Apr. 10, 2011 # Please cite Rick Chartrand, "Numerical differentiation of noisy, # nonsmooth data," ISRN Applied Mathematics, Vol. 2011, Article ID 164564, # 2011. # # Inputs: (First three required; omitting the final N parameters for N < 7 # or passing in [] results in default values being used.) # data Vector of data to be differentiated. # # iter Number of iterations to run the main loop. A stopping # condition based on the norm of the gradient vector g # below would be an easy modification. No default value. # # alph Regularization parameter. This is the main parameter # to fiddle with. Start by varying by orders of # magnitude until reasonable results are obtained. A # value to the nearest power of 10 is usally adequate. # No default value. Higher values increase # regularization strenght and improve conditioning. # # u0 Initialization of the iteration. Default value is the # naive derivative (without scaling), of appropriate # length (this being different for the two methods). # Although the solution is theoretically independent of # the initialization, a poor choice can exacerbate # conditioning issues when the linear system is solved. # # scale 'large' or 'small' (case insensitive). Default is # 'small'. 'small' has somewhat better boundary # behavior, but becomes unwieldly for data larger than # 1000 entries or so. 'large' has simpler numerics but # is more efficient for large-scale problems. 'large' is # more readily modified for higher-order derivatives, # since the implicit differentiation matrix is square. # # ep Parameter for avoiding division by zero. Default value # is 1e-6. Results should not be very sensitive to the # value. Larger values improve conditioning and # therefore speed, while smaller values give more # accurate results with sharper jumps. # # dx Grid spacing, used in the definition of the derivative # operators. Default is the reciprocal of the data size. # # plotflag Flag whether to display plot at each iteration. # Default is 1 (yes). Useful, but adds significant # running time. # # diagflag Flag whether to display diagnostics at each # iteration. Default is 1 (yes). Useful for diagnosing # preconditioning problems. When tolerance is not met, # an early iterate being best is more worrying than a # large relative residual. # # Output: # # u Estimate of the regularized derivative of data. Due to # different grid assumptions, length( u ) = # length( data ) + 1 if scale = 'small', otherwise # length( u ) = length( data ). # Copyright notice: # Copyright 2010. Los Alamos National Security, LLC. This material # was produced under U.S. Government contract DE-AC52-06NA25396 for # Los Alamos National Laboratory, which is operated by Los Alamos # National Security, LLC, for the U.S. Department of Energy. The # Government is granted for, itself and others acting on its # behalf, a paid-up, nonexclusive, irrevocable worldwide license in # this material to reproduce, prepare derivative works, and perform # publicly and display publicly. Beginning five (5) years after # (March 31, 2011) permission to assert copyright was obtained, # subject to additional five-year worldwide renewals, the # Government is granted for itself and others acting on its behalf # a paid-up, nonexclusive, irrevocable worldwide license in this # material to reproduce, prepare derivative works, distribute # copies to the public, perform publicly and display publicly, and # to permit others to do so. NEITHER THE UNITED STATES NOR THE # UNITED STATES DEPARTMENT OF ENERGY, NOR LOS ALAMOS NATIONAL # SECURITY, LLC, NOR ANY OF THEIR EMPLOYEES, MAKES ANY WARRANTY, # EXPRESS OR IMPLIED, OR ASSUMES ANY LEGAL LIABILITY OR # RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR USEFULNESS OF # ANY INFORMATION, APPARATUS, PRODUCT, OR PROCESS DISCLOSED, OR # REPRESENTS THAT ITS USE WOULD NOT INFRINGE PRIVATELY OWNED # RIGHTS. # BSD License notice: # Redistribution and use in source and binary forms, with or without # modification, are permitted provided that the following conditions # are met: # # Redistributions of source code must retain the above # copyright notice, this list of conditions and the following # disclaimer. # Redistributions in binary form must reproduce the above # copyright notice, this list of conditions and the following # disclaimer in the documentation and/or other materials # provided with the distribution. # Neither the name of Los Alamos National Security nor the names of its # contributors may be used to endorse or promote products # derived from this software without specific prior written # permission. # # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND # CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, # INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF # MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE # DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR # CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, # SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT # LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF # USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED # AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT # LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN # ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # POSSIBILITY OF SUCH DAMAGE. # ######################################################### # # # Python translation by Simone Sturniolo # # Rutherford Appleton Laboratory, STFC, UK (2014) # # simonesturniolo@gmail.com # # # ######################################################### def chop(v): return v[1:] # code starts here # Make sure we have a column vector data = np.array(data) if len(data.shape) != 1: print("Error - data is not a column vector") return # Get the data size. n = len(data) # Default checking. (u0 is done separately within each method.) if dx is None: dx = 1.0 / n # Different methods for small- and large-scale problems. if scale.lower() == "small": # Construct differentiation matrix. c = np.ones(n + 1) / dx D = sparse.spdiags([-c, c], [0, 1], n, n + 1) DT = D.transpose() # Construct antidifferentiation operator and its adjoint. def A(x): return chop(np.cumsum(x) - 0.5 * (x + x[0])) * dx def AT(w): return ( sum(w) * np.ones(n + 1) - np.transpose(np.concatenate(([sum(w) / 2.0], np.cumsum(w) - w / 2.0))) ) * dx # Default initialization is naive derivative if u0 is None: u0 = np.concatenate(([0], np.diff(data), [0])) u = u0 # Since Au( 0 ) = 0, we need to adjust. ofst = data[0] # Precompute. ATb = AT(ofst - data) # input: size n # Main loop. for ii in range(1, itern + 1): # Diagonal matrix of weights, for linearizing E-L equation. Q = sparse.spdiags(1.0 / (np.sqrt((D * u) ** 2 + ep)), 0, n, n) # Linearized diffusion matrix, also approximation of Hessian. L = dx * DT * Q * D # Gradient of functional. g = AT(A(u)) + ATb + alph * L * u # Prepare to solve linear equation. tol = 1e-4 maxit = 100 # Simple preconditioner. P = alph * sparse.spdiags(L.diagonal() + 1, 0, n + 1, n + 1) def linop(v): return alph * L * v + AT(A(v)) linop = splin.LinearOperator((n + 1, n + 1), linop) if diagflag: [s, info_i] = sparse.linalg.cg(linop, g, x0=None, tol=tol, maxiter=maxit, callback=None, M=P) print( "iteration {0:4d}: relative change = {1:.3e}, " "gradient norm = {2:.3e}\n".format( ii, np.linalg.norm(s[0]) / np.linalg.norm(u), np.linalg.norm(g) ) ) if info_i > 0: print("WARNING - convergence to tolerance not achieved!") elif info_i < 0: print("WARNING - illegal input or breakdown") else: [s, info_i] = sparse.linalg.cg(linop, g, x0=None, tol=tol, maxiter=maxit, callback=None, M=P) # Update solution. u = u - s elif scale.lower() == "large": # Construct antidifferentiation operator and its adjoint. def A(v): return np.cumsum(v) def AT(w): return sum(w) * np.ones(len(w)) - np.transpose(np.concatenate(([0.0], np.cumsum(w[:-1])))) # Construct differentiation matrix. c = np.ones(n) D = sparse.spdiags([-c, c], [0, 1], n, n) / dx mask = np.ones((n, n)) mask[-1, -1] = 0.0 D = sparse.dia_matrix(D.multiply(mask)) DT = D.transpose() # Since Au( 0 ) = 0, we need to adjust. data = data - data[0] # Default initialization is naive derivative. if u0 is None: u0 = np.concatenate(([0], np.diff(data))) u = u0 # Precompute. ATd = AT(data) # Main loop. for ii in range(1, itern + 1): # Diagonal matrix of weights, for linearizing E-L equation. Q = sparse.spdiags(1.0 / np.sqrt((D * u) ** 2.0 + ep), 0, n, n) # Linearized diffusion matrix, also approximation of Hessian. L = DT * Q * D # Gradient of functional. g = AT(A(u)) - ATd g = g + alph * L * u # Build preconditioner. c = np.cumsum(range(n, 0, -1)) B = alph * L + sparse.spdiags(c[::-1], 0, n, n) # droptol = 1.0e-2 R = sparse.dia_matrix(np.linalg.cholesky(B.todense())) # Prepare to solve linear equation. tol = 1.0e-4 maxit = 100 def linop(v): return alph * L * v + AT(A(v)) linop = splin.LinearOperator((n, n), linop) if diagflag: [s, info_i] = sparse.linalg.cg( linop, -g, x0=None, tol=tol, maxiter=maxit, callback=None, M=np.dot(R.transpose(), R) ) print( "iteration {0:4d}: relative change = {1:.3e}, " "gradient norm = {2:.3e}\n".format( ii, np.linalg.norm(s[0]) / np.linalg.norm(u), np.linalg.norm(g) ) ) if info_i > 0: print("WARNING - convergence to tolerance not achieved!") elif info_i < 0: print("WARNING - illegal input or breakdown") else: [s, info_i] = sparse.linalg.cg( linop, -g, x0=None, tol=tol, maxiter=maxit, callback=None, M=np.dot(R.transpose(), R) ) # Update current solution u = u + s u = u / dx return u PK!j derivative/loc.pyimport numpy as np from scipy.special import comb def finitediff(x, y, **kwargs): dx = x[1] - x[0] dy = np.zeros_like(x) dy[1:-1] = (y[2:] - y[:-2]) / (2.0 * dx) dy[0] = (-3.0 / 2 * y[0] + 2 * y[1] - y[2] / 2) / dx dy[-1] = (3.0 / 2 * y[-1] - 2 * y[-2] + y[-3] / 2) / dx return dy def holoborodko(x, y, M=2): """ https://github.com/jslavin/holoborodko_diff Implementation of Pavel Holoborodko's method of "Smooth noise-robust differentiators" see http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/ smooth-low-noise-differentiators Creates a numerical approximation to the first derivative of a function defined by data points. End point approximations are found from approximations of lower order. Greater smoothing is achieved by using a larger value for the order parameter, M. Parameters ---------- x : float array or scalar abscissa values of function or, if scalar, uniform step size y : float array ordinate values of function (same length as x if x is an array) M : int, optional (default = 2) order for the differentiator - will use surrounding 2*M + 1 points in creating the approximation to the derivative Returns ------- dydx : float array numerical derivative of the function of same size as y """ def coeffs(M): """ Generate the "Smooth noise-robust differentiators" as defined in Pavel Holoborodko's formula for c_k Parameters ---------- M : int the order of the differentiator c : float array of length M coefficents for k = 1 to M """ m = (2 * M - 2) / 2 k = np.arange(1, M + 1) c = 1.0 / 2.0 ** (2 * m + 1) * (comb(2 * m, m - k + 1) - comb(2 * m, m - k - 1)) return c if np.isscalar(x): x = x * np.arange(len(y)) N = 2 * M + 1 m = (N - 3) / 2 c = coeffs(M) df = np.zeros_like(y) nf = len(y) fk = np.zeros((M, (nf - 2 * M))) for i, cc in enumerate(c): # k runs from 1 to M k = i + 1 ill = M - k ilr = M + k iul = -M - k # this formulation is needed for the case the k = M, where the desired # index is the last one -- but range must be given as [-2*M:None] to # include that last point iur = (-M + k) or None fk[i, :] = 2.0 * k * cc * (y[ilr:iur] - y[ill:iul]) / (x[ilr:iur] - x[ill:iul]) df[M:-M] = fk.sum(axis=0) # may want to incorporate a variety of methods for getting edge values, # e.g. setting them to 0 or just using closest value with M of the ends. # For now we recursively calculate values closer to the edge with # progressively lower order approximations -- which is in some sense # ideal, though maybe not for all cases if M > 1: dflo = holoborodko(x[: 2 * M], y[: 2 * M], M=M - 1) dfhi = holoborodko(x[-2 * M :], y[-2 * M :], M=M - 1) df[:M] = dflo[:M] df[-M:] = dfhi[-M:] else: df[0] = (y[1] - y[0]) / (x[1] - x[0]) df[-1] = (y[-1] - y[-2]) / (x[-1] - x[-2]) return df # from finitediff import interpolate_by_finite_diff as ifd # # def finitediff(t, x, order=2, n_points=1, n_fit_ratio=10): # """Compute time-derivative of the data matrix X along first axis. # Args: # t : shape(n_samples) # x : array-like, shape (n_samples, n_features) # Input variables to be derived. # """ # # t_fit = np.linspace(t[0], t[-1], num=len(t) * n_fit_ratio, endpoint=True) # # t_fit[0] = t_fit[0] + (t_fit[1] - t_fit[0]) / 2 # t_fit[-1] = t_fit[-2] + (t_fit[-1] - t_fit[-2]) / 2 # return ifd(t, x, t_fit, order, n_points, n_points)[::n_fit_ratio, ..., 1:] PK!--"derivative-0.1.0.dist-info/LICENSEMIT License Copyright (c) 2019 Markus Quade Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. PK!HڽTU derivative-0.1.0.dist-info/WHEEL A н#Z;/"d&F[xzw@Zpy3Fv]\fi4WZ^EgM_-]#0(q7PK!Hڐg#derivative-0.1.0.dist-info/METADATANAEԄi0LH`!-+N?:o꺧VnCGxG,/<}m6[)g7ZbbC~/r6-*隍6`_w <6@UIDqȢKupV2Z+pvqщg_|d,w>ɿFV$U2өE 9qp_mt֒1dGumbEP_cC&