code: # -*- 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 #############################
  
      
apis: numpy.sum numpy.zeros numpy.argmin numpy.abs numpy.asarray numpy.max numpy.min numpy.array numpy.linalg.pinv numpy.concatenate numpy.meshgrid numpy.linspace numpy.convolve numpy.int