model_fgg-hybrid-sigm-sing.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
#

# NOT FINISHED!

# 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 sigm(x, x0, k, a, c):
    """
    Signmoidal model 

    Input:  x0,k,a,c
    Return: linear function
    """
    return a / (1 + np.exp(-k*(x-x0))) + c

# ============================================================================
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 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
            }

    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]

    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=[]

    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) 

    xd = np.array(xx); yd = np.array(yy)
    xd00 = np.array(x00); yd00 = np.array(y00)
    xd0 = np.array(x0); yd0 = np.array(y0)

    p0=[0,0,0,0]
    p0[0]=np.min(xx)+(np.max(xx)-np.min(xx))/2
    p0[1]=-0.1
    p0[2]=1
    p0[3]=1

    try:
       popt, pcov = curve_fit(sigm, xd00, yd00, p0)
    except Exception as e:
       return {'cm_return':1, 'cm_error':'Curve fit problem ('+format(e)+')'}

    pxd00 = np.linspace(np.min(xx), np.max(xx), 50)
    pyd00 = sigm(pxd00,*popt)

    if i.get('plot','')=='yes':
       pylab.plot(xd00, yd00, 'o', label='original data')
       pylab.plot(xd0, yd0, '*', label='singularities')
       pylab.plot(pxd00,pyd00, label='sigm')
       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={'x0':str(popt[0]), 'k':str(popt[1]), 'a':str(popt[2]), 'c':str(popt[3]),
           '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',{})
    x0=float(ar['x0']); k=float(ar['k']); a=float(ar['a']); c=float(ar['c'])

    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]))
    for q in l:
        v0=d0[q]
        
        xx=dx[q]
        x=float(xx)

        if str(x) in singularities:
           y=singularities[str(x)]
        else:
           y = sigm(x,x0,k,a,c)
        
        v1=str(y)

        if tp=='float' or tp=='integer':
           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/len(l))

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

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

    return rr