Source code for pygimli.frameworks.harmfit

# -*- coding: utf-8 -*-
import numpy as np

import pygimli as pg


[docs]class HarmFunctor(object): """"""
[docs] def __init__(self, A, coeff, xmin, xSpan): self.A_ = A self.coeff_ = coeff self.xmin_ = xmin self.xSpan_ = xSpan
def __call__(self, x): nc = len(self.coeff_) / 2 A = np.ones(nc * 2) * 0.5 A[1] = 3. * x A[2::2] = np.sin(2.0 * np.pi * np.arange(1, nc) * (x - self.xmin_) / self.xSpan_) A[3::2] = np.cos(2.0 * np.pi * np.arange(1, nc) * (x - self.xmin_) / self.xSpan_) return sum(A * self.coeff_)
[docs]def harmfitNative(y, x=None, nc=None, xc=None, err=None): """ python based curve-fit by harmonic functions yc = harmfitNativ(x,y[,nc,xc,err]) y .. values of a curve to be fitted x .. abscissa, if none [0..len(y)) nc .. number of coefficients xc .. abscissa to fit on (otherwise equal to x) err .. data error """ y = np.asarray(y) if x is None: x = list(range(len(y))) x = np.asarray(x) if xc is not None: xc = np.asarray(xc) else: xc = x if err is None: err = np.ones((1, len(x))) # ones(size(x)) end if nc is None or nc == 0: # nc=round(length(x)/30) % number of coefficients nc = int(len(x) / 30) xspan = max(x) - min(x) xmi = min(x) # A=ones(length(x),nc*2+2)/2 %nc*(sin+cos)+offset+drift A = np.ones((len(x), nc * 2 + 2)) / 2.0 # %nc*(sin+cos)+offset+drift # B=ones(length(xc),nc*2+2)/2 B = np.ones((len(xc), nc * 2 + 2)) / 2.0 A[:, 1] = x[:] * 3.0 B[:, 1] = xc[:] * 3.0 for i in range(1, nc + 1): # A(:,i*2+1)=sin(2*i*pi*(x-xmi)/xspan) # A(:,i*2+2)=cos(2*i*pi*(x-xmi)/xspan) # B(:,i*2+1)=sin(2*i*pi*(xc(:)-xmi)/xspan) # B(:,i*2+2)=cos(2*i*pi*(xc(:)-xmi)/xspan) A[:, i * 2 + 0] = np.sin(2.0 * i * np.pi * (x - xmi) / xspan) A[:, i * 2 + 1] = np.cos(2.0 * i * np.pi * (x - xmi) / xspan) B[:, i * 2 + 0] = np.sin(2.0 * i * np.pi * (xc - xmi) / xspan) B[:, i * 2 + 1] = np.cos(2.0 * i * np.pi * (xc - xmi) / xspan) # w = 1.0 / err w(~isfinite(w))=0 w = 1.0 / err # w[(~isfinite(w)]=0 # coeff=(spdiags(w,0,length(w),length(w))*A)\(y.*w) coeff, res, rank, s = np.linalg.lstsq(np.diag(w, 0) * A, (y * w)[0, :], rcond=None) return sum((B * coeff).T), HarmFunctor(A, coeff, xmi, xspan)
[docs]def harmfit(y, x=None, error=None, nc=42, resample=None, lam=0.1, window=None, verbose=False, dosave=False, lineSearch=True, robust=False, maxiter=20): """HARMFIT - GIMLi based curve-fit by harmonic functions Parameters ---------- y : 1d-array - values to be fitted x : 1d-array(len(y)) - data abscissa data. default: [0 .. len(y)) error : 1d-array(len(y)) error of y. default (absolute error = 0.01) nc : int - Number of harmonic coefficients resample : 1d-array - resample y to x using fitting coeffients window : int - just fit data inside window bounds Returns ------- response : 1d-array(len(resample) or len(x)) - smoothed values coefficients : 1d-array - fitting coefficients """ if x is None: x = np.arange(len(y)) xToFit = None yToFit = None if window is not None: idx = pg.find((x >= window[0]) & (x < window[1])) # idx = getIndex(x , lambda v: v > window[0] and v < window[1]) xToFit = x(idx) yToFit = y(idx) if error is not None: error = error(idx) else: xToFit = x yToFit = y fop = pg.core.HarmonicModelling(nc, xToFit, verbose) inv = pg.Inversion(yToFit, fop, verbose, dosave) if error is not None: inv.setAbsoluteError(error) else: inv.setAbsoluteError(0.01) inv.setMarquardtScheme(0.8) if error is not None: inv.stopAtChi1(True) inv.setLambda(lam) inv.setMaxIter(maxiter) inv.setLineSearch(lineSearch) inv.setRobustData(robust) # inv.setConstraintType(0) coeff = inv.run() if resample is not None: ret = fop.response(coeff, resample) if window is not None: # print pg.find((resample < window[0]) | (resample >= window[1])) ret.setVal(0.0, pg.find((resample < window[0]) | (resample >= window[1]))) # idx = getIndex(resample, # lambda v: v <= window[0] or v >= window[1]) # for i in idx: ret[i] = 0.0 return ret, coeff else: return inv.response(), coeff
# def harmfit(...)