Yosuke Tsuchiya - Fab Futures - Data Science

Home About

Week 2-1: Fitting

Assignment: Fit a function to your data

Plan

I will use Print log data of 3D Bency by Prusa Core One collected by serial communication

Review: Structure of Prusa Log Data

temperture:

T:171.37/170.00 B:85.00/85.00 X:40.15/36.00 A:49.92/0.00 @:94 B@:93 C@:36.32 HBR@:255

• T: Hotend temperture current/target
• B: Bed temperture current/target
• X: heatbreak temperture (?) current/target
• @: hotend power
• @B: bed heater power
• @C: ?? (Print Fan??)
• HBR@: hotend Fan Power (by observation)

position:

X:249.00 Y:-2.50 Z:15.00 E:7.30 Count A:24650 B:25150 Z:5950

• X: current X position of the printer
• Y: current Y position of the printer
• Z: current Z position of the printer
• E: current extrusion amount or position of the extruder
• Count section A,B,Z::??? (according to RepRap G-Code documentation, the value after "Count" seems stepper motor function).

Read the logfiles

!pip install opencv-python

Requirement already satisfied: opencv-python in /Users/yosuke/.pyenv/versions/3.13.2/lib/python3.13/site-packages (4.11.0.86)
Requirement already satisfied: numpy>=1.21.2 in /Users/yosuke/.pyenv/versions/3.13.2/lib/python3.13/site-packages (from opencv-python) (2.2.4)

[notice] A new release of pip is available: 25.1.1 -> 25.3
[notice] To update, run: pip install --upgrade pip

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from datetime import datetime,timedelta,timezone

def get_printer_temp_data(path):

    with open(path,"r") as tf:
        temp_data = tf.readlines()
    
    resjson = []
    for item in temp_data:
        tempdata = item.split(' ')
        
        addjson = {}
        logtime = tempdata[0]
        addjson['logtime'] = logtime.replace('-',' ')
        
        nozzle_temps = str(tempdata[1]).replace('T:','')
        nozzle_temp_current,nozzle_temp_setting = nozzle_temps.split('/')
        addjson['hotend_temp_current'] = float(nozzle_temp_current)
        addjson['hotend_temp_setting'] = float(nozzle_temp_setting)
        
        bed_temps = str(tempdata[2]).replace('B:','')
        bed_temp_current,bed_temp_setting = bed_temps.split('/')
        addjson['bed_temp_current'] = float(bed_temp_current)
        addjson['bed_temp_setting'] = float(bed_temp_setting)
    
        heatbreak_temps = str(tempdata[3]).replace('X:','')
        heatbreak_temp_current,heatbreak_temp_setting = heatbreak_temps.split('/')
        addjson['heatbreak_temp_current'] = float(heatbreak_temp_current)
        addjson['heatbreak_temp_setting'] = float(heatbreak_temp_setting)
    
        hotend_power = str(tempdata[5])
        hotend_power = hotend_power.replace('@:','')
        #0-255 (mapping 0.0-100.0)
        hotend_power = (float(hotend_power) / 255) * 100
        addjson['hotend_power'] = hotend_power
        
        bed_heater_power = str(tempdata[6])
        bed_heater_power = bed_heater_power.replace('B@:','')
        bed_heater_power = (float(bed_heater_power) / 255) * 100
        addjson['bed_heater_power'] = bed_heater_power

        hotend_fan_power = str(tempdata[8])
        hotend_fan_power = hotend_fan_power.replace('HBR@:','')
        hotend_fan_power = (float(hotend_fan_power) / 255) * 100
        addjson['hotend_fan_power'] = hotend_fan_power

        tp = datetime.strptime(logtime,'%Y/%m/%d-%H:%M:%S.%f')

        addjson['timestamp'] = tp
        addjson['ts_nano'] = tp.timestamp()
        
        resjson.append(addjson)
        
    df_temp = pd.DataFrame(resjson)
    
    return df_temp

def get_printer_pos_data(path):
    with open(path,"r") as pf:
        position_datas = pf.readlines()

    x_last = 0.0
    y_last = 0.0
    e_last = 0.0
    
    posjson = []
    for item in position_datas:
        posdata = item.split(' ')
    
        addjson = {}
        logtime = posdata[0]
        addjson['logtime'] = logtime.replace('-',' ')
        
        x_pos = str(posdata[1]).replace('X:','')
        x = float(x_pos)
        addjson['X_pos'] = x
        
    
        y_pos = str(posdata[2]).replace('Y:','')
        y = float(y_pos)
        addjson['Y_pos'] = y

        point = np.array([x,y])
        point_last = np.array([x_last,y_last])
        travel_distance = np.linalg.norm(point - point_last)
        addjson['travel_distance'] = travel_distance
        x_last = x
        y_last = y
        
    
        z_pos = str(posdata[3]).replace('Z:','')
        addjson['Z_pos'] = float(z_pos)
    
        e_pos = str(posdata[4]).replace('E:','')
        e = float(e_pos)
        addjson['E_pos'] = e

        e_move = e - e_last
        addjson['e_move'] = e_move
        e_last = e

        count_a_pos = str(posdata[6]).replace('A:','')
        addjson['count_a_pos'] = float(count_a_pos)

        count_b_pos = str(posdata[7]).replace('B:','')
        addjson['count_b_pos'] = float(count_b_pos)

        count_z_pos = str(posdata[8]).replace('Z:','')
        addjson['count_z_pos'] = float(count_z_pos)
        
        tp = datetime.strptime(logtime,'%Y/%m/%d-%H:%M:%S.%f')
        addjson['timestamp'] = tp
        addjson['ts_nano'] = tp.timestamp()
        
        posjson.append(addjson)
    
    df_pos = pd.DataFrame(posjson)
    return df_pos

df_temp = get_printer_temp_data('./datasets/printer_data_temperature_benchy.txt')
df_pos = get_printer_pos_data('./datasets/printer_data_position_benchy.txt')

Last Week Review: Visualization

hotend_temp and hotend_power

import matplotlib.pyplot as plt

nozzle_temp = df_temp['hotend_temp_current']
hotend_power = df_temp['hotend_power']

bed_temp = df_temp['bed_temp_current']
bed_heater_power = df_temp['bed_heater_power']
hotend_fan_power = df_temp['hotend_fan_power']
logtime = df_temp['timestamp']

fig,ax1 = plt.subplots(figsize=(4,3))
ax2 = ax1.twinx()

ax1.plot(logtime,nozzle_temp,label="nozzle_temp",color="tab:red")
ax2.plot(logtime,hotend_power,label="hotend_power",color="tab:orange")
ax2.plot(logtime,hotend_fan_power,label="hotend_fan_power",color="tab:blue")

ax1.set_ylim(top=260.0)
ax1.set_ylim(bottom=0.0)
ax2.set_ylim(top=110)
ax2.set_ylim(bottom=0.0)

plt.title("nozzle_temp,hotend_power and hotend_fan_power")
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0, fontsize='x-small')
plt.show()

Benchy 3D/3D Plotting

fig = plt.figure(figsize=(5,5))
ax = fig.add_subplot(projection='3d')

x_pos = df_pos['X_pos']
y_pos = df_pos['Y_pos']
z_pos = df_pos['Z_pos']
e_pos = df_pos['E_pos']

#d_values = df_pos['timestamp'] #.astype(int) / 10**9
d_valnes = df_pos['E_pos']

ax.scatter(x_pos,y_pos,z_pos,c=e_pos,cmap='plasma',s=1)
ax.set_xlim(0.0,240.0)
ax.set_ylim(0.0,240.0)
ax.set_zlim(0.0,150.0)

ax.set_xlabel('X_pos')
ax.set_ylabel('Y_pos')
ax.set_zlabel('Z_pos')

plt.show()

fig = plt.figure(figsize=(5,4))
fig.suptitle('x,y,z position and endeffector position')

plt.subplot(2,2,1)
plt.scatter(x_pos,y_pos,c=e_pos,cmap="plasma",s=5)
plt.xlim(75,175)
plt.ylim(75,175)
plt.title("x_pos,y_pos")

plt.subplot(2,2,2)
plt.scatter(x_pos,z_pos,c=e_pos,cmap="plasma",s=5)
plt.xlim(75,180)
plt.ylim(0,75)
plt.title("x_pos,z_pos")

plt.subplot(2,2,3)
plt.scatter(y_pos,z_pos,c=e_pos,cmap="plasma",s=5)
plt.xlim(75,150)
plt.ylim(0,75)
plt.title("y_pos,z_pos")

plt.show()

Visualizing Extruder Movement

fig,ax1 = plt.subplots(figsize=(5,3))

e_pos = df_pos['E_pos']
logtime = df_pos['timestamp']
ax1.scatter(logtime,e_pos,label="e_pos",c=e_pos,cmap="plasma",s=10)

ax1.set_ylim(top=120.0)
ax1.set_ylim(bottom=-1.0)
ax1.set_xlabel("time")
ax1.set_ylabel("Extruder move")

plt.show()

Trial 1: Hotend_Temperature x Hotend_Power

polynomial fit

strategy: use the printer log that are 10 minites later after starting the printing for fitting. Then, evaluate the log in all printing time.

df_temp_filter = df_temp.query('timestamp >= "2025/11/24 21:40:00.000000" & timestamp <= "2025/11/24 21:45:00"')
temps = df_temp_filter['hotend_temp_current']
temps = temps.to_numpy()
hotend_power = df_temp_filter['hotend_power']
timestamp = df_temp_filter['ts_nano']
timestamp = timestamp.to_numpy()

plt.figure(figsize=(4,3))
plt.plot(temps,hotend_power,'o',label="real data")
plt.xlabel("hotend temperature")
plt.ylabel("hotend power")
plt.title("after 10 - 15 min since started printing",)

plt.ylim(0.0,100.0)
plt.legend()
plt.show()

n=3
coeff1 = np.polyfit(temps,hotend_power,1)
coeff2 = np.polyfit(temps,hotend_power,2)
coeffn = np.polyfit(temps,hotend_power,n)
print(coeff1,coeff2,coeffn)

plt.figure(figsize=(4,3))

pfita = np.poly1d(coeff1)
yfita = pfita(temps)
plt.plot(temps,hotend_power,'o',label="real data")
plt.plot(temps,yfita,'r-',label="coeff 1")

pfitb = np.poly1d(coeff2)
yfitb = pfitb(temps)
plt.plot(temps,yfitb,'g-',label="coeff 2")

pfitc = np.poly1d(coeffn)
yfitc = pfitc(temps)


plt.plot(temps,yfitc,'y-',label=f"coeff {n}")

plt.xlabel("hotend temperature")
plt.ylabel("hotend power")
plt.title("Polyfit after 10 - 15 min since started printing")
plt.ylim(0.0,100.0)
plt.legend()
plt.show()

[ -10.67181976 2782.7472887 ] [ 4.77662259e+00 -2.44738626e+03  3.13544790e+05] [ 4.43083274e+00 -3.38557285e+03  8.62285634e+05 -7.32052004e+07]

#続いて、プリント開始（21:30）からプリント終了（22:14ごろまで）のログを出力
df_temp_exp = df_temp.query('timestamp >= "2025/11/24 21:30:00.000000" & timestamp <= "2025/11/24 22:14:00"')

xfit = df_temp_exp['hotend_temp_current'].to_numpy()
hotend_exp_real = df_temp_exp['hotend_power'].to_numpy()

pfit1 = np.poly1d(coeff1)
yfit1 = pfit1(xfit)

fig = plt.figure(figsize=(6,5))
fig.patch.set_facecolor('white')

plt.rcParams["axes.facecolor"] = 'white'
plt.rcParams["axes.edgecolor"] = 'black'

plt.plot(xfit,hotend_exp_real,'o',label="real data")
plt.plot(xfit,yfit1,'r-',label="coeff 1")

pfit2 = np.poly1d(coeff2)
yfit2 = pfit2(xfit)
plt.plot(xfit,yfit2,'g-',label="coeff 2",)

pfitn = np.poly1d(coeffn)
yfitn = pfitn(xfit)
plt.plot(xfit,yfitn,'y-',label="coeff 3")

plt.xlabel("hotend temperature")
plt.ylabel("hotend power")
plt.title("Polyfit whole time in printing")
plt.ylim(0.0,100.0)

plt.legend()
#plt.savefig('./notupload/2-1-polyfit.png')
plt.show()

Radial Basis function fit

df_temp_filter = df_temp.query('timestamp >= "2025/11/24 21:30:00.000000" & timestamp <= "2025/11/24 22:14:00"')

temps = df_temp_filter['hotend_temp_current']
temps = temps.to_numpy()
hotend_power = df_temp_filter['hotend_power']
hotend_power.to_numpy()
timestamp = df_temp_filter['ts_nano']
timestamp = timestamp.to_numpy()

print(np.shape(temps),np.shape(timestamp))

ncenter = 300
indices = np.random.uniform(low=0,high=len(temps),size=ncenter).astype(int)
centers = temps[indices]

M = np.abs(np.outer(temps,np.ones(ncenter))
           - np.outer(np.ones(len(temps)),centers))**3
b,residuals,rank,values = np.linalg.lstsq(M,hotend_power)

xfit = np.linspace(np.min(temps),np.max(temps),50)
yfit = (np.abs(np.outer(xfit,np.ones(ncenter)) 
               - np.outer(np.ones(50),centers))**3)@b

plt.figure(figsize=(5,4))
plt.plot(temps,hotend_power,'o',label="real data")
plt.plot(xfit,yfit,'go',label="RBF Fitting")
plt.ylim(30.0,100.0)

plt.xlabel("hotend temperature")
plt.ylabel("hotend power")
plt.title("Radial Basis Function Fit of whole printing time")
plt.ylim(0.0,100.0)
plt.legend()
plt.show()

(2617,) (2617,)

Nonlinear least squares

import numpy as np
from scipy.optimize import least_squares
import matplotlib.pyplot as plt

npts = 100
x = temps
y = hotend_power
coeff = np.array([1,0.5,5,2])

def f(coeff,x):
    return (coeff[0]*(coeff[1]+np.tanh(coeff[2]*(x-coeff[3]))))
def residuals(coeff,x,y):
    return f(coeff,x) - y

result2 = least_squares(residuals,coeff,args=(x,y),max_nfev=2)
result10 = least_squares(residuals,coeff,args=(x,y),max_nfev=6)
resultend = least_squares(residuals,coeff,args=(x,y))

fig = plt.figure(figsize=(5,4))
fig.canvas.header_visible = False
plt.plot(x,y,'bo',label='data')
plt.plot(x,f(coeff,x),'b-',label='start')
plt.plot(x,f(result2.x,x),'c-',label='2 evaluations')
plt.plot(x,f(result10.x,x),'g-',label='6 evaluations')
plt.plot(x,f(resultend.x,x),'r-',label='end')
plt.xlabel("hotend temperature")
plt.ylabel("hotend power")
plt.title("Radial Basis Function Fit of whole printing time")

plt.ylim(0.0,100.0)
plt.legend()
plt.show()

Relation between "Travel Distance of X,Y" and "Extruder movement"

strategy:

use the data after 10 min ~ 15 min since printing started.

• calculate "travel distance" between (x_pos_current,y_pos_current) and (x_pos_last,y_pos_last) with euclidean distance
• calculate "e_move"(endeffector move) between e_move and e_move_last

Apply it to whole data in printing time.

Polynomial Fit

First, fit the logdata between 10-15 min after printing started.

df_pos_filter = df_pos.query('timestamp >= "2025/11/24 21:40:00" & timestamp <= "2025/11/24 21:45:00"')

df_pos_filter = df_pos_filter.query('e_move > 0.01 & e_move < 5.0 & travel_distance > 0.05')
travel_distance = df_pos_filter['travel_distance']
travel_distance = travel_distance.to_numpy()
e_move = df_pos_filter['e_move']
e_move = e_move.to_numpy()
timestamp = df_pos_filter['ts_nano']
timestamp = timestamp.to_numpy()

plt.figure(figsize=(5,4))
plt.title("10 - 15 min after printing started.")
plt.plot(travel_distance,e_move,'o',label="real data")
plt.xlabel("travel_distance")
plt.ylabel("endeffecotr move")

plt.legend()
plt.show()

n = 10

coeff1 = np.polyfit(travel_distance,e_move,1)
coeff2 = np.polyfit(travel_distance,e_move,2)
coeffn = np.polyfit(travel_distance,e_move,n)
print(coeff1,coeff2,coeffn)

pfita = np.poly1d(coeff1)
yfit = pfita(travel_distance)

pfitb = np.poly1d(coeff2)
yfitb = pfitb(travel_distance)

pfitc = np.poly1d(coeffn)
yfitc = pfitc(travel_distance)

plt.figure(figsize=(5,4))
plt.plot(travel_distance,e_move,'o',label="real data")
plt.plot(travel_distance,yfit,'ro', label="coeff 1")
plt.plot(travel_distance,yfitb,'go',label="coeff 2")
plt.plot(travel_distance,yfitc,'yo',label="coeff 3")

plt.title("Polyfit in 10 - 15 min after printing started.")
plt.xlabel("travel_distance")
plt.ylabel("endeffecotr move")

plt.legend()
plt.show()

[0.02319148 1.40998274] [-8.65444177e-04  6.16322697e-02  1.12264543e+00] [-7.42388849e-14  2.05355954e-11 -2.22931764e-09  1.19330099e-07
 -3.04071669e-06  1.48027382e-05  1.00273650e-03 -2.21734758e-02
  1.59239609e-01 -1.98864512e-01  8.55067542e-01]

Then, apply it to whole time printing log.

df_pos_exp = df_pos.query('timestamp >= "2025/11/24 21:30:00" & timestamp <= "2025/11/24 22:14:00"')
df_pos_exp = df_pos_exp.query('e_move > 0.05 & e_move < 5.00')
travel_distance_exp = df_pos_exp['travel_distance']
e_move_exp = df_pos_exp['e_move']

travel_distance_fit = df_pos_exp['travel_distance']
pfit1 = np.poly1d(coeff1)
yfit1 = pfit1(travel_distance_fit)

plt.figure(figsize=(5,4))
plt.plot(travel_distance_exp,e_move_exp,'o',label="real data")
plt.plot(travel_distance_fit,yfit1,'ro',label="coeff 1")

pfit2 = np.poly1d(coeff2)
yfit2 = pfit2(travel_distance_fit)
plt.plot(travel_distance_fit,yfit2,'go',label="coeff 2")

pfitn = np.poly1d(coeffn)
yfitn = pfitn(travel_distance_fit)
plt.plot(travel_distance_fit,yfitn,'yo',label="coeff 3")

plt.title("Polyfit in whole printing time.")
plt.xlabel("travel_distance")
plt.ylabel("endeffecotr move")

plt.legend()
plt.show()

Radial Basis Function Fit

df_pos_filter2 = df_pos.query('timestamp >= "2025/11/24 21:30:00" & timestamp <= "2025/11/24 22:14:00"')
df_pos_filter2 = df_pos_filter2.query('e_move > 0.05 & e_move < 5.00')
travel_distance = df_pos_filter2['travel_distance']
e_move = df_pos_filter2['e_move']

ncenter = 500
indices = np.random.uniform(low=0,high=len(travel_distance),size=ncenter).astype(int)
centers = temps[indices]

M = np.abs(np.outer(travel_distance,np.ones(ncenter))
           - np.outer(np.ones(len(travel_distance)),centers))**3
b,residuals,rank,values = np.linalg.lstsq(M,e_move)

xfit = np.linspace(np.min(travel_distance),np.max(travel_distance),50)
yfit = (np.abs(np.outer(xfit,np.ones(ncenter)) 
               - np.outer(np.ones(50),centers))**3)@b

plt.figure(figsize=(5,4))
plt.plot(travel_distance,e_move,'o',label="real data")
plt.plot(xfit,yfit,'co',label=f"RBF Fit (centers={ncenter})")

plt.title("Radial Basis Function Fit in whole printing time.")
plt.xlabel("travel_distance")
plt.ylabel("endeffecotr move")

plt.legend()
plt.show()

Nonliner Least Square Fit

import numpy as np
from scipy.optimize import least_squares
import matplotlib.pyplot as plt

npts = 100
x = travel_distance
y = e_move
coeff = np.array([1,0.5,5,2])

def f(coeff,x):
    return (coeff[0]*(coeff[1]+np.tanh(coeff[2]*(x-coeff[3]))))
def residuals(coeff,x,y):
    return f(coeff,x) - y

result2 = least_squares(residuals,coeff,args=(x,y),max_nfev=2)
result10 = least_squares(residuals,coeff,args=(x,y),max_nfev=10)
resultend = least_squares(residuals,coeff,args=(x,y))

fig = plt.figure(figsize=(5,4))
fig.canvas.header_visible = False
plt.plot(x,y,'o',label='data')
plt.plot(x,f(coeff,x),'bo',label='start')
plt.plot(x,f(result2.x,x),'co',label='2 evaluations')
plt.plot(x,f(result10.x,x),'go',label='10 evaluations')
plt.plot(x,f(resultend.x,x),'ro',label='end')

plt.title("Nonlinear Least Square Fit in whole printing time.")
plt.xlabel("travel_distance")
plt.ylabel("endeffecotr move")
plt.legend()
plt.savefig('./notupload/2-1-nonlinear.png')
plt.show()