<jupyter_start><jupyter_text>Ikeda for many ships<jupyter_code># %load imports.py
"""
These is the standard setup for the notebooks.
"""

%matplotlib inline
%load_ext autoreload
%autoreload 2

from jupyterthemes import jtplot
jtplot.style(theme='onedork', context='notebook', ticks=True, grid=False)

import pandas as pd
pd.options.display.max_rows = 999
pd.options.display.max_columns = 999
pd.set_option("display.max_columns", None)
import numpy as np
import os
import matplotlib.pyplot as plt
from collections import OrderedDict
#plt.style.use('paper')

#import data
import copy
from mdldb.run import Run

from sklearn.pipeline import Pipeline
from rolldecayestimators.transformers import CutTransformer, LowpassFilterDerivatorTransformer, ScaleFactorTransformer, OffsetTransformer
from rolldecayestimators.direct_estimator_cubic import EstimatorQuadraticB, EstimatorCubic
from rolldecayestimators.ikeda_estimator import IkedaQuadraticEstimator
import rolldecayestimators.equations as equations
import rolldecayestimators.lambdas as lambdas
from rolldecayestimators.substitute_dynamic_symbols import lambdify
import rolldecayestimators.symbols as symbols
import sympy as sp

from sympy.physics.vector.printing import vpprint, vlatex
from IPython.display import display, Math, Latex

from sklearn.metrics import r2_score
from src.data import database
from mdldb import tables
import shipflowmotionshelpers.shipflowmotionshelpers as helpers
import src.visualization.visualize as visualize

from pyscores2.indata import Indata
from pyscores2.runScores2 import Calculation
from pyscores2.output import OutputFile
from pyscores2 import TDPError
import pyscores2
from rolldecayestimators.ikeda import Ikeda, IkedaR
from rolldecayestimators.simplified_ikeda_class import SimplifiedIkeda
import subprocess
from rolldecayestimators import measure
from copy import deepcopy
df_all_sections_id = pd.read_csv('../data/external/all_sections.csv', sep=';')
df_all_sections_id.head()
section_groups=df_all_sections_id.groupby(by='loading_condition_id')
loading_condition_ids = df_all_sections_id['loading_condition_id'].unique()
mask=pd.notnull(loading_condition_ids)
loading_condition_ids=loading_condition_ids[mask]
df_rolldecay = database.load(rolldecay_table_name='rolldecay_quadratic_b', limit_score=0.99, 
                             exclude_table_name='rolldecay_exclude')

df_rolldecay['g'] = 9.81
df_rolldecay['V'] = df_rolldecay['ship_speed']*1.852/3.6     
df_rolldecay['Fn'] = df_rolldecay['V']/np.sqrt(df_rolldecay['lpp']*df_rolldecay['g'])
mask=df_rolldecay['loading_condition_id'].isin(loading_condition_ids)
df=df_rolldecay.loc[mask].copy()
df['BKB'].fillna(0, inplace=True)
def create_ikeda(row, indata, output_file, fi_a):

    w = row.omega0
    scale_factor=row.scale_factor
    V = row.ship_speed*1.852/3.6/np.sqrt(scale_factor)
    R = 0.01*row.beam/scale_factor
    lBK=row.BKL/scale_factor
    bBK=row.BKB/scale_factor
    kg = row.kg/scale_factor
    
    ikeda = Ikeda.load_scoresII(V=V, w=w, fi_a=fi_a, indata=indata, output_file=output_file, 
                                scale_factor=scale_factor, BKL=lBK, BKB=bBK, kg=kg, sigma_limit=0.99)
        
    return ikeda
def calculate_ikeda(ikeda):

    output = pd.DataFrame()
    output['B_44_hat']   = ikeda.calculate_B44()
    output['B_W0_hat']   = float(ikeda.calculate_B_W0())
    output['B_W_hat']    = float(ikeda.calculate_B_W())
    output['B_F_hat']    = ikeda.calculate_B_F()
    output['B_E_hat']    = ikeda.calculate_B_E()
    output['B_BK_hat']   = ikeda.calculate_B_BK()
    output['B_L_hat']    = float(ikeda.calculate_B_L())
    output['Bw_div_Bw0'] = float(ikeda.calculate_Bw_div_Bw0())
    return output
def get_B_1_B2(s1,s2,fi_as:np.ndarray):    
    
    # Derive linear and quadratic part for all components:
    
    if isinstance(s1,pd.Series):
        result=pd.Series(name=row.name)
        columns = s1.index
    elif isinstance(s1,pd.DataFrame):
        result=pd.DataFrame()
        columns = s1.columns
    else:
        raise ValueError('s1 must be pd.Series or pd.DataFrame')
    
    x = fi_as*8*row.omega0/(3 * np.pi)
    B_2 = (s2 - s1) / (x[1] - x[0])
    B_1 = s1 - B_2 * x[0]
    
    # Save all of the component as one linear term: _1 and a quadratic term: _2
    
    
    for key in columns:
        new_name_1 = '%s_1' % key
        result[new_name_1] = s1[key]

        new_name_2 = '%s_2' % key
        result[new_name_2] = s2[key]
        
    return result
results = pd.DataFrame()
fi_as = np.deg2rad([1,10])
ikedas = {}
for run_name, row in df.iterrows():
    loading_condition_id=row['loading_condition_id']
    
    save_name='%s.in' % row.loading_condition_id
    save_path=os.path.join('scores',save_name)
    indata = Indata()
    if os.path.exists(save_path):
        indata.open(indataPath=save_path)
    else:
        continue
    
    #calculation = Calculation(outDataDirectory='scores/result')
    output_file_name = 'loading_condition_id_%i.out' % loading_condition_id
    output_file_path = os.path.join(r'scores/result', output_file_name)
    
    coefficients_file_name = 'loading_condition_id_%i-COEFF.out' % loading_condition_id
    coefficients_file_path = os.path.join(r'scores/result', coefficients_file_name)
    
    if os.path.exists(coefficients_file_path) and (os.path.getsize(coefficients_file_path) > 0):
        output_file = OutputFile(filePath=output_file_path)
    else:
        continue
    
    ikeda = create_ikeda(row=row, indata=indata, output_file=output_file, fi_a=fi_as)
    ikedas[run_name] = ikeda
    result_datas = ikeda.calculate()
    
    # Derive linear and quadratic part for all components:
    s1=result_datas.iloc[0]
    s2=result_datas.iloc[1]
    result = get_B_1_B2(s1=s1, s2=s2, fi_as=fi_as)
    result['id'] = run_name
    
    results=results.append(result)
    
results.head()<jupyter_output><empty_output><jupyter_text>Make comparison with model tests<jupyter_code>df2 = df.copy()
scale_factor = df2['scale_factor']
df2['Volume']/=(scale_factor**3)
df2['Disp']=df2['Volume']
df2['beam']/=scale_factor
df2['lpp']/=scale_factor
df2['kg']/=scale_factor
df2['BKB']/=scale_factor
df2['BKL']/=scale_factor
df['Disp']=df2['Volume']
results['Disp']=df2.loc[results.index,'Disp']
results['omega0']=df2.loc[results.index,'omega0']

#df_compare = measure.linearized_matrix(df_rolldecay=df2, df_ikeda=results, do_hatify=False)
#mask=df_compare['B_e_hat_ikeda'].notnull()
#df_compare=df_compare.loc[mask].copy()
def linearize_model_tests(df_rolldecay, phi_as = np.deg2rad(np.linspace(1,10,10)), g=9.81, rho=1000):

    df = pd.DataFrame()
    
    for phi_a in phi_as:
        df_ = measure.linearize_model_test(phi_a=phi_a, df_rolldecay=df_rolldecay, g=g, rho=rho)
        df_['phi_a']=phi_a
        df =df.append(df_, ignore_index=True)
        
    return df

def linearize_ikedas(df_ikeda, phi_as = np.deg2rad(np.linspace(1,10,10)), g=9.81, rho=1000,
                     components = ['B_44_hat', 'B_F_hat', 'B_W_hat', 'B_E_hat', 'B_BK_hat', 'B_L_hat']):

    df = pd.DataFrame()
    
    
    for phi_a in phi_as:
        df_ = measure.linearize_si(phi_a=phi_a, df_ikeda=df_ikeda, do_hatify=False, components=components)
        df_['phi_a']=phi_a
        df =df.append(df_, ignore_index=True)
        
    return df
df_rolldecays=linearize_model_tests(df_rolldecay=df2)
df_ikedas=linearize_ikedas(df_ikeda=results)
df_ikedas.dropna(subset=['B_e_hat'], inplace=True)
df_rolldecays.describe()
df_ikedas.describe()
df_compare = pd.merge(left=df_rolldecays, right=df_ikedas, how='inner', left_on=('id','phi_a'), right_on=('id','phi_a'), suffixes=('','_ikeda'))<jupyter_output><empty_output><jupyter_text><jupyter_code>fig,ax=plt.subplots()
df_compare.plot(x='B_e_hat', y='B_e_hat_ikeda', ax=ax, style='o', alpha=0.30, ms=10, markeredgewidth=0)

xlim = ax.get_xlim()
ylim = ax.get_ylim()
lim = np.max([xlim[1],ylim[1]])
ax.set_xlim(0,lim)
ax.set_ylim(0,lim)
ax.plot([0,lim],[0,lim],'r-')

ax.grid(True)
ax.set_aspect('equal', 'box')
ax.get_legend().remove()
mask = df_compare['ship_speed'] < 1
df_compare_zero_speed = df_compare.loc[mask]
df_compare_speed = df_compare.loc[~mask]<jupyter_output><empty_output><jupyter_text><jupyter_code>fig,ax=plt.subplots()
df_compare_zero_speed.plot(x='B_e_hat', y='B_e_hat_ikeda', ax=ax, style='o', alpha=0.30, ms=10, markeredgewidth=0)

ax.set_xlabel(r'$\hat{B_{e}}^{Model}$')
ax.set_ylabel(r'$\hat{B_{e}}^{Ikeda}$')

xlim = ax.get_xlim()
ylim = ax.get_ylim()
lim = np.max([xlim[1],ylim[1]])
ax.set_xlim(0,lim)
ax.set_ylim(0,lim)
ax.plot([0,lim],[0,lim],'r-')

ax.grid(True)
ax.set_aspect('equal', 'box')
ax.get_legend().remove()
fig,ax=plt.subplots()
df_compare_speed.plot(x='B_e_hat', y='B_e_hat_ikeda', ax=ax, style='o', alpha=0.30, ms=10, markeredgewidth=0)

ax.set_xlabel(r'$\hat{B_{e}}^{Model}$')
ax.set_ylabel(r'$\hat{B_{e}}^{Ikeda}$')

xlim = ax.get_xlim()
ylim = ax.get_ylim()
lim = np.max([xlim[1],ylim[1]])
ax.set_xlim(0,lim)
ax.set_ylim(0,lim)
ax.plot([0,lim],[0,lim],'r-')

ax.grid(True)
ax.set_aspect('equal', 'box')
ax.get_legend().remove()
fig,ax=plt.subplots()

mask_bilge_keel = df_compare_zero_speed.BKL > 0

A0 = df_compare_zero_speed['A0']
bins = np.linspace(A0.min(),A0.max(),20)

df_compare_zero_speed.loc[mask_bilge_keel].hist('A0', bins=bins, ax=ax, label='bilge keel', color='green')
df_compare_zero_speed.loc[~mask_bilge_keel].hist('A0', bins=bins, ax=ax, label='no bilge keel', color='red')

ax.annotate('limit', xy=(0.99,0), xytext=(0.99,5), arrowprops=dict(facecolor='black', shrink=0.05))
ax.legend()
fig,ax=plt.subplots()
for run_id,df_ in df_compare_zero_speed.loc[mask_bilge_keel].groupby(by='id'):
    df_.plot(x='B_e_hat', y='B_e_hat_ikeda', ax=ax, style='.-', alpha=0.30)


for run_id,df_ in df_compare_zero_speed.loc[~mask_bilge_keel].groupby(by='id'):
    df_.plot(x='B_e_hat', y='B_e_hat_ikeda', ax=ax, style='ro', alpha=0.70, ms=10, markeredgewidth=0)

ax.set_xlabel(r'$\hat{B_{e}}^{Model}$')
ax.set_ylabel(r'$\hat{B_{e}}^{Ikeda}$')

xlim = ax.get_xlim()
ylim = ax.get_ylim()
lim = np.max([xlim[1],ylim[1]])
ax.set_xlim(0,lim)
ax.set_ylim(0,lim)
ax.plot([0,lim],[0,lim],'r-')

ax.grid(True)
ax.set_aspect('equal', 'box')
ax.get_legend().remove()
df_ = df_compare_zero_speed.loc[mask_bilge_keel]
id=(df_['B_e_hat'] - df_['B_e_hat_ikeda']).abs().sort_values(ascending=False).index[0]
run_id = df_compare_zero_speed.loc[id]['id']
run_ids = df_compare['id'].unique()
df_ = df_rolldecay.loc[run_ids]
df_speed_variations = df_.groupby(by = 'loading_condition_id').filter(lambda x: len(x['ship_speed'].unique()) > 3)
loading_conditions = df_speed_variations.groupby(by = 'loading_condition_id').first()
#loading_conditions
loading_condition_ids = [34,107]

ymaxs = []
axess = []
for loading_condition_id in loading_condition_ids:
    
    mask = df_rolldecay.loading_condition_id == loading_condition_id
    df_runs = df_rolldecay.loc[mask]
    
    run_groups = df_runs.groupby(by='Fn')
    
    fig,axes=plt.subplots(ncols=len(run_groups))
    fig.set_size_inches(30,8)
    
    ## Amplitude variation:
    for (Fn, df),ax in zip(run_groups,axes):
        
        axess.append(ax)
        
        run_id = df.iloc[0]['id']
        ikeda = deepcopy(ikedas[run_id])
        phi_as = np.deg2rad(np.linspace(0,10,100))
    
        result = ikeda.calculate(fi_a=phi_as)
        result['phi_a'] = phi_as
        result.set_index('phi_a', inplace=True)
        
        visualize.plot_area(result, ax=ax, interesting_ = ['B_W_hat', 'B_F_hat', 'B_E_hat', 'B_L_hat', 'B_BK_hat'])
        ax.set_title('Fn:%0.2f' % Fn)
    
        ylims = ax.get_ylim()
        ymaxs.append(ylims[1])
        
ymax=np.max(ymaxs)
for ax in axess:
    ax.set_ylim((0,ymax))

ymaxs = []
axess = []
for loading_condition_id in loading_condition_ids:
    
    mask = df_rolldecay.loading_condition_id == loading_condition_id
    run_id = df_rolldecay.loc[mask].index[0]
    
    ## Speed variation:
    ikeda = deepcopy(ikedas[run_id])
    phi_as = np.deg2rad([10])
    V = np.linspace(df_runs['V'].min(), df_runs['V'].max(), 100)                     
    V/=np.sqrt(ikeda.scale_factor)
    g = 9.81
    Fn = V/np.sqrt(ikeda.lpp*g)
    
    ikeda.V = V
    #ikeda.bBK = ikeda.bBK[0]*np.ones(len(V)).T
    #ikeda.lBK = ikeda.lBK[0]*np.ones(len(V)).T
    
    result = ikeda.calculate(fi_a=phi_as)
    result['Fn'] = Fn
    result.set_index('Fn', inplace=True)
    fig,ax=plt.subplots()
    axess.append(ax)
    visualize.plot_area(result, ax=ax, interesting_ = ['B_W_hat', 'B_F_hat', 'B_E_hat', 'B_L_hat', 'B_BK_hat'])
    
    ylims = ax.get_ylim()
    ymaxs.append(ylims[1])

ymax=np.max(ymaxs)
for ax in axess:
    ax.set_ylim((0,ymax))<jupyter_output>c:\dev\rolldecay-estimators\rolldecayestimators\ikeda_speed.py:283: RuntimeWarning: divide by zero encountered in true_divide
  Cf = 1.328*sqrt(2*pi*visc/(3.22*r_f**2*fi_a**2*w))
c:\dev\rolldecay-estimators\rolldecayestimators\ikeda_speed.py:184: RuntimeWarning: invalid value encountered in true_divide
  CD = 22.5 * bBK / (pi * l * fi_a * f) + 2.4;
c:\dev\rolldecay-estimators\rolldecayestimators\ikeda_speed.py:190: RuntimeWarning: divide by zero encountered in true_divide
  So = (0.3 * pi * l * fi_a * f / bBK + 1.95) * bBK;
c:\dev\rolldecay-estimators\rolldecayestimators\ikeda_speed.py:190: RuntimeWarning: invalid value encountered in true_divide
  So = (0.3 * pi * l * fi_a * f / bBK + 1.95) * bBK;
c:\dev\rolldecay-estimators\rolldecayestimators\ikeda_speed.py:190: RuntimeWarning: invalid value encountered in multiply
  So = (0.3 * pi * l * fi_a * f / bBK + 1.95) * bBK;
c:\dev\rolldecay-estimators\rolldecayestimators\ikeda_speed.py:210: RuntimeWarning: invalid value encountered in true[...]<jupyter_text><jupyter_code>print(df_runs.iloc[0].project_path)
df_speed_variations = df_rolldecay.groupby(by = 'loading_condition_id').filter(lambda x: len(x) > 3)
loading_conditions = df_speed_variations.groupby(by = 'loading_condition_id').first()
#print(loading_conditions.loc[95].project_path)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14717439451310996, 'fraction_numerical': 0.021433164249482033, 'mean_word_length': 4.879042419151617, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 24, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '229697', 'license': 'MIT', 'n_tokens_mistral': 5887, 'n_tokens_neox': 5674, 'n_words': 832, 'path': 'notebooks/10.2_ikeda_many.ipynb', 'repo_name': 'rddaz2013/Prediction-of-roll-motion-using-fully-nonlinear-potential-flow-and-Ikedas-method', 'size': '13997'}