PrivacyScanData/truncate/new2/pyrod-master/tool/tools.py

# -*- coding: utf-8 -*-
"""
Created on Wed Aug 15 21:50:58 2018

@author: USER
"""

# Codes are free to use. Do whatever you want

from __future__ import absolute_import

"""Read raw data"""

####################### LIBRARY #############################

# exceptions library
from exceptions import (Data_Format_Exception,
                        Data_Match_Exception)

# Python stdlib imports
import datetime
from math import factorial

# data processing library
import numpy as np

# pyrod library

####################### CONSTANT ############################

# constant 

####################### FUNCTIONS ###########################

'.......................optimise.........................'

# f - fitting data
# y - experiment data
# mask - mask data

def R_square(f, y, mask):
    
    if not len(f) == len(y) == len(mask):
        raise Data_Match_Exception('Please input equal length')
    
    def nplist(data):
        
        # check and transform data
        try:
            
            # check np array
            if isinstance(data, np.ndarray):
                pass
            # check list
            elif isinstance(data, list):
                rl = np.array(data)
            # check np mat
            elif isinstance(data, np.matrix):
                rl = np.asarray(data).reshape(-1)
            # for other unpoackable datatype
            else:
                # init a list first
                l = []
                # unpack raw data with for
                for e in data:
                    l.append(e)
                # trans to np array
                rl = np.array(l)
                
        # unknown type
        except Data_Format_Exception:
            
            print('unknown data type')
            
        return rl

    # tranform to np array; apply mask 
    rf, ry = nplist(f)*nplist(mask), nplist(y)*nplist(mask)

    # calculate r square
    ss_tot = np.sum((ry - np.sum(ry)/len(ry))**2)
    ss_res = np.sum((ry - rf)**2)
    
    r2 = 1 - ss_res/ss_tot
    
    return r2


def opt_step_brute(func,x0_range,grid_size = 10,step = 2):
    
    """
       Brute method is much too slow and big.
       However, its usefull and simple. To improve it, we try to step it
       
       x0_range: range of variable, [x1-,x1+],[x2-,x2+]
       currently,only two axes are avaialble
    """
    # current step is 3
    step = 3
    
    # grid_size and step have to be integer
    try:
        grid_size = int(grid_size)
        step = int(step)
        
    except ValueError:
        raise ValueError("grid_size and step have to be of type int")
    
    # one dimensional step brute method
    if len(x0_range) == 1:
        
        # store func(grid_data) result
        grid_list0 = []
        x0 = np.linspace(x0_range[0][0],x0_range[0][1],grid_size)
        
        # func(grid_data)
        for px in range(grid_size):
            grid_list0.append(func(x0[px]))
        # store min in step1
        min_idx = np.argmin(grid_list0)
        
        # continue step2
        grid_list1 = []
        x1 = x0[min_idx]
        delta = (abs(x0_range[0][1] - x0_range[0][0]))/grid_size
        
        x2 = np.linspace(x1-delta,x1+delta,grid_size)
        for sx in range(grid_size):
            grid_list1.append(func(x2[sx]))
            
        min_step2 = x2[np.argmin(grid_list1)]
        
    elif len(x0_range) == 2:
    
        # step1: grid the x0_range
        min_step1 = []
        au = np.linspace(x0_range[0][0],x0_range[0][1],grid_size)
        av = np.linspace(x0_range[1][0],x0_range[1][1],grid_size)
        
        # find minimum in xu and xv grid
        def grid_min(xu,xv):
            
            x0_grid = np.meshgrid(xu, xv)
            
            #grid list
            grid_list = np.mat(np.zeros([grid_size**2,3]))
            idx = 0
            
            # pu-- for postion in u axes
            for pu in range(grid_size):
                # pv--for postion in v axes
                for pv in range(grid_size):
                    
                    grid_list[idx,0] = x0_grid[0][pu,pv]
                    grid_list[idx,1] = x0_grid[1][pu,pv]
                    grid_list[idx,2] = func([x0_grid[0][pu,pv],
                                             x0_grid[1][pu,pv]])
                    idx = idx + 1
            # find the minimum in step1
            min_idx = np.argmin(grid_list[:,2])
            
            return grid_list[min_idx,:]
        
        # append the firt minimum before rocking
        min_step1.append(grid_min(au,av))
        
        # start rocking, try to avoid local minmum
        bu = au - (au[1]-au[0])/2
        bv = av - (av[1]-av[0])/2
        
        min_step1.append(grid_min(bu,bv))
        
        # step 2
        # step 2 new x range
        u_min = np.min([min_step1[0][0,0],
                        min_step1[1][0,0]])
        u_max = np.max([min_step1[0][0,0],
                        min_step1[1][0,0]])
        deta_u = u_max - u_min
        v_min = np.min([min_step1[0][0,1],
                        min_step1[1][0,1]])
        v_max = np.max([min_step1[0][0,1],
                        min_step1[1][0,1]])
        deta_v = v_max - v_min
        # new u and v
        cu = np.linspace(u_min-deta_u, u_min+deta_u, grid_size)
        cv = np.linspace(v_min-deta_v, v_min+deta_v, grid_size)
        
        min_step2 = grid_min(cu,cv).tolist()
    
    return min_step2
    
            
'......................smooth.........................'

def savitzky_golay(y, window_size, order, deriv=0, rate=1):
    
    """    
       Smooth (and optionally differentiate) data with a Savitzky-Golay filter.
       The Savitzky-Golay filter removes high frequency noise from data.
       It has the advantage of preserving the original shape and
       features of the signal better than other types of filtering
       approaches, such as moving averages techniques.

      ----------
      .. [1] A. Savitzky, M. J. E. Golay, Smoothing and Differentiation of
         Data by Simplified Least Squares Procedures. Analytical
         Chemistry, 1964, 36 (8), pp 1627-1639.
      .. [2] Numerical Recipes 3rd Edition: The Art of Scientific Computing
         W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery
         Cambridge University Press ISBN-13: 9780521880688
    """
    
    # integer value
    try:
        window_size = np.abs(np.int(window_size))
        order = np.abs(np.int(order))
    except ValueError:
        raise ValueError("window_size and order have to be of type int")
     
    if window_size % 2 != 1 or window_size < 1:
        raise TypeError("window_size size must be a positive odd number")
    if window_size < order + 2:
        raise TypeError("window_size is too small for the polynomials order")
        
    order_range = range(order+1)
    half_window = (window_size -1) // 2
    
    # precompute coefficients
    
    b = np.mat([[k**i for i in order_range] for k in range(-half_window, half_window+1)])
    m = np.linalg.pinv(b).A[deriv] * rate**deriv * factorial(deriv)
    
    # pad the signal at the extremes with
    # values taken from the signal itself
    
    firstvals = y[0] - np.abs( y[1:half_window+1][::-1] - y[0] )
    lastvals = y[-1] + np.abs(y[-half_window-1:-1][::-1] - y[-1])
    y = np.concatenate((firstvals, y, lastvals))
     
    return np.convolve( m[::-1], y, mode='valid')

######################## CLASSS #############################