<jupyter_start><jupyter_text>Simplified Ikeda sensitivity<jupyter_code>#from jupyterthemes import jtplot
#jtplot.style(theme='onedork', context='notebook', ticks=True, grid=False)
%matplotlib inline
%load_ext autoreload
%autoreload 2
import pandas as pd
pd.options.display.max_rows = 999
pd.options.display.max_columns = 999
import numpy as np
import matplotlib.pyplot as plt
from pylab import rcParams
rcParams['figure.figsize'] = 15, 5

from rolldecayestimators.simplified_ikeda import calculate_roll_damping, SimplifiedIkedaInputError
from rolldecayestimators import equations
import sympy as sp
from rolldecayestimators import symbols
from rolldecayestimators.substitute_dynamic_symbols import lambdify
from mdldb.tables import Run
from rolldecayestimators.ikeda_estimator import IkedaEstimator
from rolldecayestimators.transformers import CutTransformer, LowpassFilterDerivatorTransformer, ScaleFactorTransformer, OffsetTransformer
from sklearn.pipeline import Pipeline

from rolldecay import database

import data
import copy
from rolldecay import database
from sympy.physics.vector.printing import vpprint, vlatex
from IPython.display import display, Math, Latex
from rolldecay.paper_writing import save_fig
from rolldecay.sensitivity import calculate
db = database.get_db()
sql =  """
SELECT * from
loading_conditions
INNER JOIN models
    ON loading_conditions.model_number==models.model_number
    INNER JOIN ships
        ON ships.name == models.ship_name
"""

df = pd.read_sql_query(sql=sql, con=db.engine)
mask = df['Volume']>0
df=df.loc[mask].copy()
df.describe()
df['DRAFT']=(df.TA + df.TF)/2
df['CB'] = df.Volume/(df.lpp*df.beam*df.DRAFT)
mask = df['CB'].abs()!=np.inf
df=df.loc[mask]
df['CB'].hist(bins=30)
mask = ((0 < df['CB']) & (df['CB'] <=1))
df = df.loc[mask]
mask = df['Volume'] == df['Volume'].median()
ship = df.loc[mask].iloc[0]
ship['phi_max'] = 5
ship['omega0'] = 0.5
ship.Volume/(ship.lpp*ship.beam*(ship.TA+ship.TF)/2)
from rolldecay.sensitivity import variate_ship
N=30
changes = np.linspace(0.5,1.5,N)
df = variate_ship(ship=ship, key='lpp', changes=changes)
result = df.apply(func=calculate, axis=1)
result.plot()
from rolldecay.sensitivity import plot_variation, calculate_variation
plot_variation(ship=ship, key='lpp')
variations = ['lpp','beam','DRAFT','phi_max','BKL','BKB','omega0','kg','CB','A0']
N=30
changes = np.linspace(0.85,1.32,N)

ylims=[]
axes=[]
for variation in variations:
    
    fig,ax=plt.subplots()
    plot_variation(ship=ship, key=variation, ax=ax, changes=changes, catch_error=True)
    ax.set_title('Variation of %s: %0.2f' % (variation, ship[variation]))
    
    ylims.append(list(ax.get_ylim()))
    axes.append(ax)
    
for ax in axes:
    ax.set_ylim(np.min(ylims),np.max(ylims))
    
plot_variation(ship=ship, key='lpp')
ikeda_parameters = ['lpp','beam','DRAFT','phi_max','BKL','BKB','omega0','kg','CB','A0']
ship_=ship[ikeda_parameters]
ship_['phi_max']=8
ship_
calculate(ship_)
ship2 = ship_.copy()
scale = 0.01
ship2['lpp']/=scale
ship2['beam']/=scale
ship2['DRAFT']/=scale
ship2['BKL']/=scale
ship2['BKB']/=scale
ship2['kg']/=scale
ship2['omega0']*=np.sqrt(scale)
calculate(ship2)
ship2 = ship_.copy()
scale = ship_['lpp']
ship2['lpp']/=scale
ship2['beam']/=scale
ship2['DRAFT']/=scale
ship2['BKL']/=scale
ship2['BKB']/=scale
ship2['kg']/=scale
ship2['omega0']*=np.sqrt(scale)

#ship2['BKB']*=10
calculate(ship2)
variations = list(set(ikeda_parameters) - set(['lpp']))

N=100
changes = np.linspace(0.85,1.43,N)

ylims=[]
axes=[]
for variation in variations:
    
    fig,ax=plt.subplots()
    plot_variation(ship=ship2, key=variation, ax=ax, changes=changes, catch_error=True, 
                   plot_change_factor=False)
     
    ylims.append(list(ax.get_ylim()))
    axes.append(ax)
    ax.grid(True)
    
for ax in axes:
    ax.set_ylim(np.min(ylims),np.max(ylims))

df = variate_ship(ship=ship2, key='beam', changes=changes)
result = calculate_variation(df=df, catch_error=True)
result = pd.concat([result,df], axis=1)
result['Disp'] = result['CB']*result['lpp']*result['beam']*result['DRAFT']
equations.B44_hat_equation
B44_equation = sp.Eq(symbols.B_44,sp.solve(equations.B44_hat_equation,symbols.B_44)[0])
B44_equation
B44_lambda = lambdify(sp.solve(B44_equation,symbols.B_44)[0])
result['B44'] = B44_lambda(result['B44HAT'], result['Disp'], beam=result['beam'], g=9.81, rho=1000)
result.plot(x='beam', y = 'B44')
variations = list(set(ikeda_parameters) - set(['lpp']))

N=100
changes = np.linspace(0.7,1.32,N)

ylims=[]
axes=[]
for variation in variations:
    
    fig,ax=plt.subplots()
    df = variate_ship(ship=ship2, key=variation, changes=changes)
    result = calculate_variation(df=df, catch_error=True)
    #result.fillna(0,inplace=True)
        
    df['Disp'] = df['CB']*df['lpp']*df['beam']*df['DRAFT']
    
    dampings = ['B44','BW','BE','BBK']
    for key in dampings:
        result[key] = B44_lambda(np.array(result[key+'HAT']), np.array(df['Disp']), beam=df['beam'].astype(float),
                                 g=9.81, rho=1000)
    
    result[variation]=df[variation]
    result.plot(x=variation, y=dampings, ax=ax)
    ax.plot(ship2[variation],0,'ro', label='reference')
    
    ylims.append(list(ax.get_ylim()))
    axes.append(ax)
    
for ax in axes:
    ax.set_ylim(np.min(ylims),np.max(ylims))
from rolldecayestimators.ikeda_speed import bilge_keel
variation='BKB'
N=100
changes = np.linspace(0.5,3,N)
df = variate_ship(ship=ship2, key=variation, changes=changes)
result = calculate_variation(df=df, catch_error=True)
dddf = variate_ship(ship=ship2, key=variation, changes=changes)
ddresult = calculate_variation(df=df, catch_error=True)
R=0.3/175
df=df.astype(float)

def calculate_bilge_keel(row,R):
    OG=-row.kg + row.DRAFT
    Ho=(row['beam']/(2*row['DRAFT'])  )
    A=(row['A0']*row['beam']*row['DRAFT'])
    
    Bp44BK_N0, Bp44BK_H0, B44BK_L, B44BKW0 = bilge_keel(w=row['omega0'], 
               fi_a=np.deg2rad(row['phi_max']), V=0, B=row['beam'], 
               d=row['DRAFT'], 
               A=A, bBK=row['BKB'], R=R, g=9.81, OG=OG, Ho=Ho, ra=1000)
    
    s = pd.Series(name=row.name)
    s['Bp44BK_N0'] = Bp44BK_N0
    s['Bp44BK_H0'] = Bp44BK_H0
    s['B44BK_L'] = B44BK_L
    s['B44BKW0'] = B44BKW0
    
    LBK = row.BKL
    B44BK_N0 = Bp44BK_N0*LBK
    B44BK_H0 = Bp44BK_H0*LBK
    B44BK_L = B44BK_L
    # B44BKW0 = B44BKW0 * dim...
    B44_BK = B44BK_N0 + B44BK_H0 + B44BK_L
    s['B44_BK'] = B44_BK
    
    return s
result2 = df.apply(func=calculate_bilge_keel, R=R, axis=1)
key='BBK'
df['Disp'] = df['CB']*df['lpp']*df['beam']*df['DRAFT']
result[key] = B44_lambda(np.array(result[key+'HAT']), np.array(df['Disp']), beam=df['beam'].astype(float),
                                 g=9.81, rho=1000)
fig,ax=plt.subplots()
result.plot(y='BBK', ax=ax)
result2.plot(y='B44_BK', ax=ax)
df_rolldecay = database.load(rolldecay_table_name='rolldecay_quadratic_b', limit_score=0.90,
                             exclude_table_name='rolldecay_exclude')
db = database.get_db()
loading_conditions=df_rolldecay.copy()
scale = np.array(loading_conditions['lpp'])
loading_conditions['lpp']/=scale
loading_conditions['beam']/=scale
loading_conditions['DRAFT']=(loading_conditions['TA']+loading_conditions['TF'])/2/scale
loading_conditions['BKL']/=scale
loading_conditions['BKB']/=scale
loading_conditions['kg']/=scale
volume = loading_conditions['Volume']/scale**3
lpp=loading_conditions['lpp']
T=loading_conditions['DRAFT']
beam=loading_conditions['beam']
loading_conditions['phi_max']=np.rad2deg(loading_conditions['phi_start'].abs())

loading_conditions['CB']=volume/(lpp*beam*T)
variations = list(set(ikeda_parameters) - set(['lpp']))

N=100
ship2.omega0 = loading_conditions['omega0'].median()
ylimss=[]
axes=[]
fig,axes=plt.subplots(nrows=len(variations))

for variation,ax in zip(variations,axes):
    
    s = loading_conditions[variation]
    changes = np.linspace(s.min(),s.max(),N)/ship2[variation]
    
    
    df = variate_ship(ship=ship2, key=variation, changes=changes)
    result = calculate_variation(df=df, catch_error=True)
    #result.fillna(0,inplace=True)
        
    df['Disp'] = df['CB']*df['lpp']*df['beam']*df['DRAFT']
    
    dampings = ['B44','BW','BE','BBK']
    for key in dampings:
        result[key] = B44_lambda(np.array(result[key+'HAT']), np.array(df['Disp']), beam=df['beam'].astype(float),
                                 g=9.81, rho=1000)
    
    result[variation]=df[variation]
    result.plot(x=variation, y=dampings, ax=ax)
    ylims = ax.get_ylim()
    ax.plot([ship2[variation],ship2[variation]],ylims,'r--', label='median ship')
    ax.get_legend().remove()
    
axes[0].legend()
plt.tight_layout()

save_fig(fig=fig, name='ikeda_variation', full_page=True)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1720356408498972, 'fraction_numerical': 0.03164267763308202, 'mean_word_length': 5.249107780157031, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '296909', 'license': 'MIT', 'n_tokens_mistral': 3696, 'n_tokens_neox': 3459, 'n_words': 515, 'path': 'notebooks/rolldecay/04_simplified_ikeda/05.5_maa_simplified_ikeda_sensitivity.ipynb', 'repo_name': 'martinlarsalbert/rolldecay', 'size': '8754'}