model_fgg-hybrid-with-rules.py

Go to the documentation of this file.

#
# Collective Mind
#
# See cM LICENSE.txt for licensing details.
# See cM Copyright.txt for copyright details.
#
# Developer(s): (C) Grigori Fursin, started on 2011.09
#

# Should always be here
ini={}
cm_kernel=None

# Local settings
import json
import os
import csv
import math
import numpy as np
import pylab
from scipy.optimize import curve_fit

# ============================================================================
def init(i):
    return {'cm_return':0}

# ============================================================================
def is_singularity_power_of_two(x):
    """
    Separate possible singularities (for now, power of 2)

    Input:  x - float number
    Return: True if power of 2, otherwise False
    """

    y=int(x)&int(x-1)

    if y==0: r=True
    else: r=False

    return r

# ============================================================================
def lm(x, x1, a1, k1):
    """
    Linear model 

    Input:  x1,a1,k1
    Return: linear function
    """

    return a1+k1*(x-x1)

# ============================================================================
def build(i):

    """
    Build model

    Input:  {
              data                 
              ct_dimensions_input
              ct_dimensions_output
              desc
              model_name                   - (earth, svm)
              (record_data_to_file_prefix) - if !='', use this filename prefix instead of randomly generated
              additional_params            - x1 - cuts on first axis
                                             x2 - cuts on first axis
                                             x3 - cuts on first axis
            }

    Output: {
              cm_return       - return code >0 if error
              file_with_model - file with model
            }
    
    """

    # We can deal here only with one input and one output dimension for now

    data=i['data']

    dxs=i['ct_dimensions_input'][0]
    dys=i['ct_dimensions_output'][0]

    ap=i['additional_params']

    cx1=float(ap['x1'])
    cx2=float(ap['x2'])
    cx3=float(ap['x3'])

    xs=[float(q) for q in data[dxs]]
    ys=[float(q) for q in data[dys]]

    xx=np.array(xs)
    yy=np.array(ys)

    x00=[];y00=[]
    x0=[];y0=[]
    x1=[];y1=[]
    x2=[];y2=[]
    x3=[];y3=[]
    x4=[];y4=[]

    singularities={}

    for q in range(0, len(xx)):
        x=xx[q]
        y=yy[q]

        if is_singularity_power_of_two(x+1): 
           x0.append(x);y0.append(y)
           singularities[str(x)]=str(y)
        else:
          x00.append(x);y00.append(y) 
          if x<cx1:              x1.append(x);y1.append(y)
          elif x<cx2:            x2.append(x);y2.append(y)
          elif x<cx3:            x3.append(x);y3.append(y)
          else:                  x4.append(x);y4.append(y)

    xd = np.array(xx); yd = np.array(yy)
    xd00 = np.array(x00); yd00 = np.array(y00)
    xd0 = np.array(x0); yd0 = np.array(y0)
    xd1 = np.array(x1); yd1 = np.array(y1)
    xd2 = np.array(x2); yd2 = np.array(y2)
    xd3 = np.array(x3); yd3 = np.array(y3)
    xd4 = np.array(x4); yd4 = np.array(y4)

    try:
       popt1, pcov1 = curve_fit(lm, xd1, yd1)
       popt2, pcov2 = curve_fit(lm, xd2, yd2)
       popt3, pcov3 = curve_fit(lm, xd3, yd3)
       popt4, pcov4 = curve_fit(lm, xd4, yd4)
    except Exception as e:
       return {'cm_return':1, 'cm_error':'Curve fit problem ('+format(e)+')'}

    px1 = np.linspace(1, cx1, 50)
    px2 = np.linspace(cx1, cx2, 50)
    px3 = np.linspace(cx2, cx3, 50)
    px4 = np.linspace(cx3, xx.max(), 50)

    py1 = lm(px1,*popt1)
    py2 = lm(px2,*popt2)
    py3 = lm(px3,*popt3)
    py4 = lm(px4,*popt4)

    if i.get('plot','')=='yes':
       pylab.plot(xd00, yd00, 'o', label='original data')
       pylab.plot(xd0, yd0, '*', label='singularities')
       pylab.plot(px1,py1, label='lm1')
       pylab.plot(px2,py2, label='lm2')
       pylab.plot(px3,py3, label='lm3')
       pylab.plot(px4,py4, label='lm4')
       pylab.ylim(0, yy.max()+1)
       pylab.legend(loc='upper left')
       pylab.grid(True)
       pylab.show()

    # Generate output tmp file name
    fotmp=i.get('record_data_to_file_prefix','')
    if fotmp=='':
       r=cm_kernel.gen_cm_tmp_file({})
       fotmp=r['cm_path1']

    fotmp1=fotmp+'.py.model.json'

    model={'cx1':str(cx1), 'x1':str(popt1[0]), 'a1':str(popt1[1]), 'k1':str(popt1[2]),
           'cx2':str(cx2), 'x2':str(popt2[0]), 'a2':str(popt2[1]), 'k2':str(popt2[2]),
           'cx3':str(cx3), 'x3':str(popt3[0]), 'a3':str(popt3[1]), 'k3':str(popt3[2]),
           'x4':str(popt4[0]), 'a4':str(popt4[1]), 'k4':str(popt4[2]),
           'singularities':singularities}

    f=open(fotmp1,'wt')
    f.write(json.dumps(model, indent=2)+'\n')
    f.close()

    return {'cm_return':0, 'file_with_model':fotmp1, 'model':model}

# ============================================================================
def predict(i):

    """
    Predict using model

    Input:  {
              model_file 
              data                 
              ct_dimensions_input
              (ct_dimensions_output)  - for comparison
              desc                    - cM data description
              model_name              - (earth, svm)
              (max_variation_percent) - for comparison, report points where variation is more than this number (default=0.2)
            }

    Output: {
              cm_return - return code >0 if error
              (rmse)    - if comparison, root mean square error for predictions vs original
              (max_var) - list of points with variation more than max_variation_percent          
            }
    
    """

    mf=i.get('model_file','')
    if mf=='':
       return {'cm_return':1, 'cm_error':'"model_file" is not defined'}

    # Load coefficients
    f=file(i['model_file'])
    ar=json.loads(f.read())
    f.close()

    singularities=ar.get('singularities',{})
    cx1=float(ar['cx1']); x1=float(ar['x1']); a1=float(ar['a1']); k1=float(ar['k1'])
    cx2=float(ar['cx2']); x2=float(ar['x2']); a2=float(ar['a2']); k2=float(ar['k2'])
    cx3=float(ar['cx3']); x3=float(ar['x3']); a3=float(ar['a3']); k3=float(ar['k3'])
    x4=float(ar['x4']); a4=float(ar['a4']); k4=float(ar['k4'])

    data=i.get('data',{})
    desc=i.get('desc',{})
    
    dim1=i.get('ct_dimensions_input',[])
    dim2=i.get('ct_dimensions_output',[])
    
    kx=i['ct_dimensions_input'][0]
    dx=data[kx]

    ky=i['ct_dimensions_output'][0]
    d0=data[ky]

    rr={'cm_return':0}

    tp=desc.get(ky,{}).get('type','')

    mvp=i.get('max_variation_percent','0.2')
    dmvp=float(mvp)

    var=[]
    
    s=0.0
    l=range(0, len(data[ky]))
    nn=0
    for q in l:
        v0=d0[q]
        
        xx=dx[q]
        x=float(xx)

        if str(x) in singularities:
           y=singularities[str(x)]
        else:
          if x<cx1:   y=lm(x, x1, a1, k1)
          elif x<cx2: y=lm(x, x2, a2, k2)
          elif x<cx3: y=lm(x, x3, a3, k3)
          else:       y=lm(x, x4, a4, k4)
        
        v1=str(y)

        if tp=='float' or tp=='integer':
           nn+=1
           if tp=='float':
              dv0=float(v0)
              dv1=float(v1)
           else:
              dv0=int(v0)
              dv1=int(v1)
           s+=(dv0-dv1)*(dv0-dv1)
           diff=abs(dv1-dv0)/dv0
           c1=''
           if diff>dmvp: 
              c1=' ***'
              var.append(q)
           print "%7s" % data[dim1[0]][q], "%7.3f" % dv0, "%7.3f" % dv1, "%7.3f" % s, "%5.3f" % diff,c1
#        else:
#           if v0!=v1:
#              s+=1

        d0[q]=v1

    rmse=math.sqrt(s/nn)

    rr['rmse']=str(rmse)
    rr['max_var']=var

    print ''
    print 'Model RMSE='+str(rmse)

    return rr