text:
import os

import numpy as np

import plotly.plotly as py
import plotly.graph_objs as go
import plotly.figure_factory as ff

from scipy.special import erf
from scipy import stats

import pandas as pd
import seaborn as sns
cp = sns.color_palette()

from sklearn.neighbors.kde import KernelDensity

#%%
# loading data and setting prelims:

rsi_df = pd.read_csv('data.csv')

r = rsi_df['0'].values

d = 2
beta = 0.05
n = 1
alpha_tcc = np.log(r)/(-n*d) - beta*d
gard_tcc_per_nd =  (alpha_tcc+beta*d)
rxdose_tcc = np.array(33/(alpha_tcc+beta*d))

t = np.arange(0,5,.1)
rsi_interval = np.round(np.arange(0.01,.81,.01),2)
dose_range = np.round(np.arange(40,82,2),0)

#%% dependency functions

def construct_kde(array, bandwidth=None):
    if bandwidth == None:
        bw = 1.2*array.std()*np.power(array.size,-1/5)
    else:
        bw = bandwidth
    kde = KernelDensity(kernel='gaussian', bandwidth=bw)
    kde.fit(array.reshape(-1,1))
    x = np.linspace(array.min(),array.max(),200)
    log_dens=kde.score_samples(x.reshape(-1,1))
    kdens=np.exp(log_dens)

    total_dens=np.sum(kdens)
    cdf_array=np.zeros(shape=len(x))
    delta=x[1]-x[0]
    for i in range(len(x)):
        cdf_array[i] = np.sum(kdens[:i])*delta

    return x,kdens, cdf_array

def cumulative_hist(hist_obj):
    counts = hist_obj[0]
    cdf_hist_obj = np.zeros(len(counts))
    total_sum = np.sum(counts)
    for i in range(len(counts)):
        cdf_hist_obj[i]  = np.sum(counts[:i])/total_sum
    return cdf_hist_obj,hist_obj[1]

def calc_cdf(array,var,bandwidth=None):
    if bandwidth == None:
        bw = 1.2*array.std()*np.power(array.size,-1/5)
    else:
        bw = bandwidth
    kde=stats.gaussian_kde(dataset=array,bw_method=bw)
    return kde.integrate_box_1d(low=0,high=var)

#%%

""" weibull fits: """

def S1(t):
    lambda_ = 0.028907
    rho_ = 0.641806
    y = np.exp(-1*np.power(lambda_*t,rho_))
    return y

def S2(t):
    lambda_ = 0.263115
    rho_ = 1.161710
    y = np.exp(-1*np.power(lambda_*t,rho_))
    return y

wb1 = dict(x=t,y=S1(t),xaxis='x2',yaxis='y2',
           line = dict(color='rgb(.5,.5,.5)'))
wb2 = dict(x=t,y=S2(t),xaxis='x2',yaxis='y2',
           line = dict(color='rgb(.5,.5,.5)'))

""" exponential fits: """

def plc_gard33(t,G33):
    """t = time in yrs, G33 = bool"""
    if G33:
        lc_gard33 = S1(t)
    else:
        lc_gard33 = S2(t)
    return lc_gard33

#%%

"""ntcp functions: """

def risk_p(dose):
    d = np.maximum(dose-18,0)
    r = d*0.068
    return r

def prob_pneumonitis(MLD):
    fx_to_lung = 8.5
    MLDh = MLD
    MLDl = 0
    b0 = -3.87
    b1 = 0.126
    prob_h = np.exp(b0+b1*MLDh)/(1+np.exp(b0+b1*MLDh))
    prob_l = np.exp(b0+b1*MLDl)/(1+np.exp(b0+b1*MLDl))
    return np.round(prob_h - prob_l,3)

def risk_e(dose):
    d = np.maximum(dose-40,0)
    r = 0.026*d
    return r

def prob_esoph(dose_h):
    if dose_h < 7.5:
        r = .01/100
    elif (dose_h >= 7.5) & (dose_h <=20):
        r = (0.08*dose_h - 0.59)/100
    elif dose_h >20:
        EUD1 = dose_h
        EUD2 = 0
        TD50 = 47
        m = 0.36
        t1 = (EUD1-TD50)/(m*TD50)
        t2 = (EUD2-TD50)/(m*TD50)
        y = (erf(t1)-erf(t2))/2
        r = np.round(y,3)
    return r

p_mce_1gy = .074

def pfs_rsi(t,rval, dose):
    lc_rsi = exp_rsi(t,rval)
    risk_e_l = (risk_p(dose)+risk_e(dose))/100
    risk_mce = 1+p_mce_1gy*dose/14
    pfs = np.power(lc_rsi,risk_mce*np.exp(risk_e_l))
    return pfs

def pfs_gard(t,G, dose):
    """t = time in yrs, G = GARD"""
    lc_gard = exp_gard(t,G)
    risk_e_l = (risk_p(dose)+risk_e(dose))/100
    risk_mce = 1+p_mce_1gy*dose/14
    pfs = np.power(lc_gard,risk_mce*np.exp(risk_e_l))
    return pfs

def pfs_gard_man(t,G, H,L,E):
    """t = time in yrs, G = GARD"""
    lc_gard = exp_gard(t,G)

    risk_e_l = prob_pneumonitis(L)+prob_esoph(E)
    risk_mce = 1+p_mce_1gy*H
    pfs = np.power(lc_gard,risk_mce*np.exp(risk_e_l))
    return pfs

def pfs_gard33(t, G33, dose):
    """plc_gard33 = S(t) local control function w/ WB model,
        dose = total dose"""

    lc_gard33 = plc_gard33(t,G33)
    risk_e_l = (risk_p(dose)+risk_e(dose))/100
    risk_mce = 1+p_mce_1gy*dose/14
    pfs = np.power(lc_gard33,risk_mce*np.exp(risk_e_l))

    return pfs

def pfs_gard33_man(t,G33, H,L,E):
    """t = time in yrs, G = GARD """
    lc_gard33 = plc_gard33(t,G33)

    risk_e_l = prob_pneumonitis(L)+prob_esoph(E)
    risk_mce = 1+p_mce_1gy*H
    pfs = np.power(lc_gard33,risk_mce*np.exp(risk_e_l))
    return pfs

#%%

""" histogram & distribution trace objects: """

bins_rsi = list(np.arange(0,.81,.03))
rsi_hist_obj = np.histogram(r,bins=bins_rsi,density=True)
bar_rsi = go.Bar(name = 'RSI',
                x=rsi_hist_obj[1][:-2], y=rsi_hist_obj[0],
                width=.02,
                marker = {'color':'rgb(.1,.7,.35)'},
                opacity=.4,
                xaxis='x2',yaxis='y2')

bins_rxdose = list(np.arange(21,141,5))+[240]
rxdose_hist_obj = np.histogram(rxdose_tcc,bins=bins_rxdose,density=True)

bar_rxdose = go.Bar(name  = 'Rx-RSI (Gy)',
                    x=rxdose_hist_obj[1][:-2],y=rxdose_hist_obj[0],
                    width = 4,
                    marker = {'color':'rgb(140,140,140)'}, # gray
                    # marker = {'color':'rgb(.1,.6,.85)'}, # blue
                    opacity=.4,
                    hoverinfo='none',
                    xaxis='x1', yaxis='y1')

hist_rsi = go.Histogram(name='RSI',
                        x=r,
                        xbins=dict(start=0, end=.81,size= .03),
                        autobinx = False,
                        histnorm='probability density', marker = {'color':'rgb(.1,.75,.3)'}, opacity=.5,
                        xaxis='x2',yaxis='y2')

hist_rxdose = go.Histogram(x=rxdose_tcc,nbinsx=60,
                           histnorm='probability density',
                           opacity=.6,
                           # marker = {'color':'rgb(.1,.6,.3)'},
                           marker = {'color':'rgb(.1,.6,.85)'}, # lighter blue
                           xaxis='x1',yaxis='y1',
                           name='Rx-RSI (Gy)',
                           hoverinfo='none')

rxdose_kde = construct_kde(rxdose_tcc)
dist_rxdose = go.Scatter(x=rxdose_kde[0],y=rxdose_kde[1],
                xaxis='x1',yaxis='y1',
                line=dict(
                    # color='rgb(.1,.6,.85)', # blue
                    color='rgb(80,80,80)',
                    width=2.5
                ),

                hoverinfo='none',
                showlegend=False)

rsi_kde = construct_kde(r)
dist_rsi = go.Scatter(x=rsi_kde[0],y=rsi_kde[1],
                xaxis='x2',yaxis='y2',
                line=dict(color='rgb(.1,.75,.3)'),
                # line=dict(color='rgb(110,30,160)'), # purple
                showlegend=False)

"""   plot layouts:   """
#%%
rsi_range_x = [0, .8]
rsi_range_y = [0,5.5]
gard_range_x=[0,150]
gard_range_y=[0,0.033]

axes_color = 'rgb(.4,.4,.4)'
axes_width = 1
frame_on = True

layout_hists = go.Layout(
    xaxis1=dict(
        title= "RxRSI (Gy)",
        domain=[0,1],
        range=gard_range_x,
        anchor='y1',
        linecolor='rgb(.4,.4,.4)',
        linewidth=axes_width,
        mirror=frame_on,
        tickmode='array',
        tickvals = list(range(0,160,10)),
        ticktext= [0,'',20,'',40,'',60,'',80,'',100,'',120,'','',''],
        ticklen=4),
    yaxis1=dict(
        domain=[0,1],
        range=gard_range_y,
        anchor='x1',
        linecolor=axes_color,
        linewidth=axes_width,
        mirror=frame_on,
        showticklabels=False,
        showgrid=False),
    xaxis2=dict(
        title='RSI',
        titlefont={'size':12},
        linecolor=axes_color,
        domain=[0.63, .98],range=rsi_range_x,
        anchor='y2',
        tickmode='array',
        tickvals = list(np.arange(0,.9,.1)),
        ticktext= [0,'','.2','','.4','','.6','','.8'],
        ticklen=3,
        tickfont={'size':10}),
    yaxis2=dict(
        showticklabels=False,
        linecolor=axes_color,
        domain=[0.65, 1],range = rsi_range_y,
        anchor='x2',
        showgrid=False),
    width=500,
    height=500,
    margin=go.layout.Margin(
        l=40,
        b=70,
        t=50,
        r=130,
        pad=0.1),
    legend=dict(x=1.05,y=.2,
        font={'size':11},
        bgcolor= 'rgb(255,255,255)',
        bordercolor=axes_color,
        borderwidth=1),
    paper_bgcolor='rgb(.97,.97,.97)'
)

layout_outcome_curves = go.Layout(
    xaxis=dict(
        title='Time (yrs)',
        range=[0, 5.2],
        showgrid = False,
        linecolor=axes_color,
        linewidth=axes_width,
        mirror =  True,
        tickmode='array',
        tickvals = list(range(6)),
        ticktext= list('012345'),
        ticklen=4),
    yaxis=dict(
        title = 'Probability',
        range=[0, 1],
        showgrid=False,
        linecolor=axes_color,
        linewidth=axes_width,
        mirror  = True,
        tickmode='array',
        tickvals = np.arange(0,1.1,.1).tolist(),
        ticktext= [0,'','.2','','.4','','.6','','.8','','1'],
        ticklen=4),
    width=400,
    height=500,
    margin=go.layout.Margin(
        l=70,
        r=40,
        b=90,
        t=70,
        pad=.1),
    legend=dict(x=.1,y=.1,
                bgcolor= 'rgb(255,255,255)',
                bordercolor=axes_color,
                borderwidth=1),
    paper_bgcolor='rgb(.97,.97,.97)'
)

#%%

import dash
import dash_core_components as dcc
import dash_html_components as html
from dash.dependencies import Input, Output, State
from dash.exceptions import PreventUpdate

#%%

external_stylesheets = ['https://codepen.io/chriddyp/pen/bWLwgP.css']

app = dash.Dash(__name__, external_stylesheets=external_stylesheets)

server = app.server

#%%

marks={.01:'0'}
for k in rsi_interval[4::5]:
    marks[str(k)]=str(k)

th_style={'font-size':12,'paddingLeft':10,'paddingTop':5,'paddingBottom':5,'backgroundColor':'rgb(245,245,245)'}
td_style={'font-size':12,'textAlign':'center','padding':5}

input_style={
    'width':'60px','height':'30px','display':'inline-block',
    'paddingRight':5,
    'paddingLeft':8,
    'marginTop':1,
    'marginBottom':1
}
radio_args=dict(
    options=[{'label': ' ', 'value': 'auto'},
    {'label': ' ', 'value': 'manual'}],
    labelStyle={'display': 'inline-block','paddingRight':'60px'},
    value='')

input_args=dict(
    placeholder='0.0',
    type='number', value='',
    debounce=True,
    min=0.0, max=100.0, step=1,
    n_blur=0, n_submit=0)

""" --------------- app layout --------------- """
#%%
app.layout = html.Div(
    [
        # top row
        html.Div([
            html.Div(
                dcc.Graph(id='histograms'),
                style={'display':'inline-block',
                       'border-style':'solid',
                       'border-width':2,
                       'border-color':'rgb(120,120,120)'}),
            html.Div(
                dcc.Graph(id='outcome-curves'),
                style={'display':'inline-block',
                       'border-style':'solid',
                       'border-width':2,
                       'border-color':'rgb(120,120,120)'}),
            html.Div(id='table-container',children=
                [
                    html.Table(
                        [
                            html.Thead([
                                html.Tr(html.Th('Parameters',colSpan=2,style={'font-size':13,
                                                                           'background-color':'rgb(220,220,220)',
                                                                            'textAlign':'center',
                                                                            'padding':5}))
                            ], style={'height':'20px'}),
                            html.Colgroup([
                                html.Col(style={'backgroundColor':'rgb(220,220,220)','width':100}),
                                html.Col(style={'width':100})
                            ]),
                            html.Tbody([
                                html.Tr([html.Th('RSI',style=th_style), html.Td(id='rsi-output',style=td_style)]),
                                html.Tr([html.Th('Dose',style=th_style), html.Td(id='dose-output',style=td_style)]),
                                html.Tr([html.Th('RxRSI',style=th_style),html.Td(id='rxdose-output',style=td_style)]),
                                html.Tr([html.Th('GARD (Tx)',style=th_style),html.Td(id='gard-output',style=td_style)]),

                                html.Tr(html.Th('Normal Tissue Doses',style={'textAlign':'center','fontSize':13,'padding':5},colSpan=2)),

                                html.Tr([html.Th('Heart',style=th_style), html.Td(id='heart-dose-output',style=td_style)]),
                                html.Tr([html.Th('Lung',style=th_style), html.Td(id='lung-dose-output',style=td_style)]),
                                html.Tr([html.Th('Esophagus',style=th_style), html.Td(id='esoph-dose-output',style=td_style)]),

                                html.Tr(html.Th('Outcomes',style={'textAlign':'center','fontSize':13,'padding':5},colSpan=2)),

                                html.Tr([html.Th('1-yr pLC',style=th_style), html.Td(id='pefs-output-1',style=td_style)]),
                                html.Tr([html.Th('2-yr pLC',style=th_style), html.Td(id='pefs-output-2',style=td_style)]),
                                html.Tr([html.Th('5-yr pLC',style=th_style), html.Td(id='pefs-output-5',style=td_style)])

                            ])
                        ], style={'margin-top':10}
                    )
                ], style={'float':'right','height':500,'margin-right':50,
                          'border-style':'hidden',
                          'border-width':2,
                          'border-color':'rgb(120,120,120)'}
            )
        ], style = {'margin-left':40, 'width':1200,'display':'block'}),

        ####### bottom row #######
        html.Div([
            # bottom left
            html.Div(children=[
                html.H6(children='RSI and Dose Selection',style={'fontSize':20}),
                html.Div([
                    html.Div(id='display-selected-rsi'),
                    html.Div(
                        dcc.Slider(
                            id = 'rsi-slider',
                            min=rsi_interval.min(),
                            max = rsi_interval.max(),
                            value = 0.20,
                            marks = marks,
                            step = 0.01,
                            updatemode='drag'
                        )
                    ),
                    html.Br(),
                    html.Br(),
                    html.Div(id='display-selected-dose'),
                    html.Div(
                        dcc.Slider(
                            id='dose-slider',
                            min=dose_range.min(),
                            max=dose_range.max(),
                            value = 60,
                            marks = {str(k):str(k) for k in dose_range},
                            step=1,
                            updatemode='drag'
                        )
                    )
                ], style = {'fontSize':14}),

            ], style={'width':450,'display':'inline-block'}),

            # bottom right
            html.Div(children=[
                    html.H6('Normal Tissue Dose',style={'fontSize':20}),

                    html.Table([
                        html.Thead([
                            html.Tr(children=[html.Th(children=['Site',html.Br(),'(mean dose)'],rowSpan=2,style={'paddingLeft':'10px','paddingBottom':'5px','paddingTop':'5px'}),
                                              html.Th('Calculation Method:',colSpan=2,style={'padding':'5px'}),
                                              html.Th(children=['Percent of',html.Br(),'Total Dose'],rowSpan=2,style={'padding':'5px'})],
                                    style={'fontSize':13,'height':20,
                                           'color':'rgb(245,245,245)','backgroundColor':'rgb(120,120,120)'}),

                            html.Tr(children=[
                                              html.Td('Use Defaults', style={'padding':'5px'}),
                                              html.Td('Manually Enter',style={'padding':'5px'})
                                              ],
                                    style={'fontSize':12,'height':20,
                                           'color':'rgb(245,245,245)','backgroundColor':'rgb(120,120,120)'})
                        ]),
                        html.Colgroup([
                            html.Col(style={'width':80,'backgroundColor':'rgb(240,240,240)'}),
                            html.Col(style={'width':90}),
                            html.Col(style={'width':100}),
                            html.Col(style={'width':115})
                        ]),
                        html.Tbody([

                            html.Tr([
                                html.Td('Apply to All',
                                        style={'font-size':11,'font-weight':'bold',
                                               'paddingLeft':'10px','paddingTop':'3px','paddingBottom':'3px',
                                               'color':'rgb(40,40,40)',
                                               'backgroundColor':'#ffffff'}),
                                html.Td([
                                    dcc.RadioItems(id = 'apply-all',
                                                  options=[{'label': '', 'value': 'auto'},
                                                  {'label': '', 'value': 'manual'}],
                                                  labelStyle={'display': 'inline-block','paddingRight':'60px'},
                                                  value='initial',
                                                  style={'display':'inline-block'})
                                    ], colSpan=2 , style={'paddingBottom':5,'paddingTop':5,'paddingLeft':'15px'}),
                                html.Td([
                                    html.Button(
                                        id='clear-button',
                                        n_clicks=0,
                                        children=html.Label('X',style={'position':'relative','top':'-10px','z-index':-1}),
                                        style={'width':20,'height':20,
                                               'padding':0,
                                               'fontSize':11,
                                               'fontWeight':'bold',
                                               # 'color':'rgb(180,20,20)',
                                               'backgroundColor':'#d8d8d8',
                                               'vertical-align':'text-bottom',
                                               'z-index':1
                                               }
                                    ),
                                    html.Div(
                                        children=['Clear'],
                                        style={'font-weight':'bold','font-size':12,'display':'inline-block','marginLeft':5}
                                    )
                                ], style={'paddingBottom':5,'paddingTop':5})
                            ], style={'height':'30px'}
                            ),

                            html.Tr([
                                html.Th(['Heart'],style={'paddingLeft':'10px'}),

                                html.Td([
                                    html.Div(id='heart_radio-table-cell',children=[
                                        dcc.RadioItems(id = 'entry-method-heart',
                                                       **radio_args,
                                                      style={'display':'inline-block'}),
                                    ]),
                                ],colSpan=2),

                                html.Td([
                                    html.Div(id='heart_input-container',children=[
                                        dcc.Input(id = 'heart-dose-input',
                                                  **input_args,
                                                  style=input_style),
                                        html.Div('%',style={'fontSize':13,'display':'inline','marginLeft':5})
                                    ],hidden=True),
                                    html.Div(id='heart-default-percent',hidden=True)
                                ], style={'paddingRight':15,'paddingTop':1,'paddingBottom':1})
                            ]),

                            html.Tr([
                                html.Th(['Lung'],style={'paddingLeft':'10px'}),
                                html.Td([
                                    html.Div(id='lung_radio-container',children=[
                                        dcc.RadioItems(id = 'entry-method-lung',
                                                       **radio_args,
                                                      style={'display':'inline-block'})
                                    ]),
                                ], colSpan=2),

                                html.Td([
                                    html.Div(id='lung_input-container',children=[
                                        dcc.Input(id = 'lung-dose-input',
                                                  **input_args,
                                                  style=input_style),
                                        html.Div('%',style={'fontSize':13,'display':'inline','marginLeft':5})
                                    ],hidden=True),
                                    html.Div(id='lung-default-percent',hidden=True)
                                ], style={'paddingRight':15,'paddingTop':1,'paddingBottom':1})
                            ]),

                            html.Tr([
                                html.Th(['Esophagus'],style={'paddingLeft':'10px'}),
                                html.Td([
                                    html.Div(id='esoph_radio-container',children=[
                                        dcc.RadioItems(id = 'entry-method-esoph',**radio_args,style={'display':'inline-block'})
                                    ]),
                                ], colSpan=2),
                                html.Td([
                                    html.Div(id='esoph_input-container',children=[
                                        dcc.Input(id = 'esoph-dose-input',**input_args,style=input_style),
                                        html.Div('%',style={'fontSize':13,'display':'inline','marginLeft':5})
                                    ],hidden=True),
                                    html.Div(id='esoph-default-percent',hidden=True)
                                ], style={'paddingRight':15,'paddingTop':1,'paddingBottom':1})
                            ])
                        ], style={'fontSize':13})
                    ]),
                    # end table

                    html.Div(id='display-store-data',style={'display':'inline-block'})

            ], style={'width':500,'margin-left':10,'margin-right':50,'float':'right'}
             ),
            dcc.Store(id='dose-output-store'),
            dcc.Store(id='memory-store'),

        ], style={'width':1050, 'margin-left':50,'margin-top':5, 'display':'block'})

    ], style={'margin-top':40}

)

@app.callback([
    Output('display-selected-rsi','children'),
    Output('display-selected-dose','children')],
    [Input('rsi-slider','value'),
    Input('dose-slider','value')])

def update_slider_text(selected_rsi,selected_dose):
    return 'RSI: %.2f' %selected_rsi, 'Treatment Dose: {} (Gy)'.format(selected_dose)

""" ********************* update hist plots ********************* """
@app.callback(
    Output('histograms','figure'),
    [Input('rsi-slider','value'),
    Input('dose-slider','value')]
)

def update_hist_figures(selected_rsi,selected_dose):
    rval = selected_rsi
    alpha_val = np.log(rval)/(-n*d)-beta*d
    rxdose_val = 33/(alpha_val+beta*d)
    dose_val = selected_dose
    gard_val = dose_val*(alpha_val+beta*d)
    G33 = True if (dose_val>=rxdose_val) else False

    gard_tcc = gard_tcc_per_nd*dose_val
    hist_gard = go.Histogram(x=gard_tcc,nbinsx=60,histnorm='probability density',opacity=.6, marker = {'color':'rgb(.1,.6,.3)'},xaxis='x1',yaxis='y1')
    traces = [bar_rsi, dist_rsi, bar_rxdose, dist_rxdose]

    ticker_color = 'rgb(.95,.6,.2)'
    traces.append(go.Scatter(
        name='Selected RSI',
        xaxis='x2',yaxis='y2',
        line = {'color':ticker_color,'width':1.5,'dash':'dash'},
        y=np.linspace(0,6,50),
        x=np.full((50),rval))
    )
    traces.append(go.Scatter(
        name = 'RxRSI',
        xaxis='x1', yaxis='y1',
        line = {'color':ticker_color,'width':2.75},
        y=np.linspace(0,gard_range_y[1],50),
        x=np.full((50),rxdose_val),
        hoverinfo='none'
        )
    )
    traces.append(go.Scatter(
        name = 'Treatment Dose',
        xaxis='x1', yaxis='y1',
        line = {'color':'rgb(.95,.2,.1)','width':2.75},
        y=np.linspace(0,gard_range_y[1],50),
        x=np.full((50),dose_val),
        hoverinfo='none')
    )

    # rxrsi hoverlabel and kde shading #
    idx = np.argwhere(rxdose_kde[0]>=rxdose_val)[0][0]
    kde_yval_rxdose=rxdose_kde[1][idx]
    percentile = np.round(calc_cdf(rxdose_tcc,rxdose_val,bandwidth=.1)*100,1)
    display_rxrsi=np.int(np.round(rxdose_val,0))
    traces.append(go.Scatter(
        name='RxRSI',
        x=[rxdose_val],
        y=[kde_yval_rxdose],
        mode='markers',
        marker=dict(color=ticker_color,size=8),
        text='<b>{} Gy</b><br>{}%'.format(display_rxrsi,percentile),
        hoverinfo='text+name',
        hoverlabel=dict(font={'size':12,'color':'white'}),
        showlegend=False
    ))

    # dose hoverlabel #
    j = np.argwhere(rxdose_kde[0]>=dose_val)[0][0]
    kde_yval_dose=rxdose_kde[1][j]
    percentile_td = np.round(calc_cdf(rxdose_tcc,dose_val,bandwidth=.1)*100,1)
    traces.append(go.Scatter(
        name='Dose',
        x=[dose_val],
        y=[kde_yval_dose],
        xaxis='x1', yaxis='y1',
        mode='markers',
        marker=dict(color='rgb(.95,.2,.1)',size=8),
        text='<b>{} Gy</b><br>{}%'.format(dose_val,percentile_td),
        hoverinfo='text+name',
        hoverlabel=dict(font={'size':12}),
        showlegend=False
        )
    )

    xmin=np.min([rxdose_val,dose_val])
    xmax=np.max([rxdose_val,dose_val])
    k = np.argwhere((rxdose_kde[0]>=xmin) & (rxdose_kde[0]<=xmax)).flatten()
    xvals_fill = rxdose_kde[0][k]
    yvals_fill=np.full(shape=(len(k)),fill_value=gard_range_y[1])
    fill_color = 'rgba(200,230,255,.2)' if rxdose_val<dose_val else 'rgba(255,180,180,.2)'
    traces.append(
        go.Scatter(
            x=xvals_fill,
            y=yvals_fill,
            xaxis='x1', yaxis='y1',
            line=dict(color='rgb(250,250,250)',width=0),
            fill='tozeroy',
            fillcolor=fill_color,
            mode='lines',
            hoverinfo='none',
            showlegend=False
        )
    )

    return {
        'data':traces,
        'layout':layout_hists
    }
"""---------------------------------------------------------------------"""

""" *************** add rsi, gard, dose to output table *************** """
"""---------------------------------------------------------------------"""
@app.callback([Output('rsi-output','children'),
               Output('dose-output','children'),
               Output('rxdose-output','children'),
               Output('gard-output','children')],
               [Input('rsi-slider','value'),
                Input('dose-slider','value')]
)
def update_output_table_rsi_gard(selected_rsi,selected_dose):
    rval = np.round(selected_rsi,2)
    alpha_val = np.log(rval)/(-n*d)-beta*d
    rxdose_val = np.round(33/(alpha_val+beta*d),0)
    dose_val = selected_dose
    gard_val = np.round(dose_val*(alpha_val+beta*d),0)
    G33 = True if (dose_val>=rxdose_val) else False

    return '%.2f' %rval, '{} Gy'.format(dose_val), '%.0f Gy' % rxdose_val, gard_val
"""------------------------------------------------------------------"""

""" ********************* plot outcome curve *********************
------------------------------------------------------------------"""
@app.callback([Output('outcome-curves','figure'),
               Output('pefs-output-1','children'),
               Output('pefs-output-2','children'),
              Output('pefs-output-5','children')],
              [Input('rsi-slider','value'),
               Input('dose-slider','value'),
               Input('dose-output-store','data'),
               # Input('heart-dose-output','children'),
               # Input('lung-dose-output','children'),
               # Input('esoph-dose-output','children')
               ]
)
def update_outcome_figure(selected_rsi,selected_dose,data):   #hdose,ldose,edose):
    rval = selected_rsi
    alpha_val = np.log(rval)/(-n*d)-beta*d
    rxdose_val = 33/(alpha_val+beta*d)
    dose_val = selected_dose
    gard_val = dose_val*(alpha_val+beta*d)
    G33 = True if (dose_val>=rxdose_val) else False

    traces = []
    tb_output = []
    normal_tissue_doses = np.array(data)

    if (np.any(normal_tissue_doses=='')) or (np.any(normal_tissue_doses==None)):
        all_doses_entered = False
    else:
        all_doses_entered = True
        hdose=data[0]
        ldose=data[1]
        edose=data[2]

    if all_doses_entered:
        pfs_array = pfs_gard33_man(t,G33,hdose,ldose,edose)
        traces.append(go.Scatter(
            name='Penalized-LC',
            line = {'color':'rgb(.8,.2,.1)'},
            x=t,
            y=pfs_array,
            visible=True,
            showlegend=True))
        tb_output.append(np.round(pfs_array[9],2))
        tb_output.append(np.round(pfs_array[19],2))
        tb_output.append(np.round(pfs_array[49],2))
        vis=True
    else:
        tb_output=['']*3
        vis=False

    traces.append(go.Scatter(
        name='Local Control',
        line = dict(color='rgb(.4,.1,.5)'),
        x=t,
        y=plc_gard33(t,G33),
        visible=False,
        showlegend=True))

    layout = layout_outcome_curves
    if all_doses_entered:
        layout['annotations']=None
    else:
        layout['annotations']=[dict(
            x=2.4,
            y=.6,
            xref='x',
            yref='y',
            text='<i>Enter Normal<br>  Tissue Doses</i>',
            showarrow=False,
            font = {'size':24,'color':'rgb(190,190,190)'},
            align='left'
        )]

    return [{'data':traces,'layout':layout}] + tb_output
"""------------------------------------------------------------------"""

""" ********************* add doses output table *********************"""
@app.callback([Output('heart-dose-output','children'),
               Output('lung-dose-output','children'),
               Output('esoph-dose-output','children'),
               Output('dose-output-store','data')],
                [Input('dose-slider','value'),
                 Input('entry-method-heart','value'),
                 Input('entry-method-lung','value'),
                 Input('entry-method-esoph','value'),
                 Input('heart-dose-input','value'),
                 Input('lung-dose-input','value'),
                 Input('esoph-dose-input','value')])
def update_output_table_doses(selected_dose,h_entry_method,l_entry_method,e_entry_method, hdose, ldose, edose):

    dose_val = selected_dose
    entry_methods = [h_entry_method,l_entry_method,e_entry_method]
    percent_total = [7.0,11.8,25.2]
    manual_percent= [hdose, ldose, edose]

    display_site_doses = []
    store_site_doses = []

    for i in range(len(entry_methods)):
        if (entry_methods[i]=='auto'):
            site_dose_i = selected_dose*percent_total[i]/100
            store_site_doses.append(site_dose_i)
            display_site_doses.append('%.1f Gy' % site_dose_i)
        elif (entry_methods[i] == 'manual') and (manual_percent[i]!=''):
            site_dose_i= selected_dose*manual_percent[i]/100
            store_site_doses.append(site_dose_i)
            display_site_doses.append('%.1f Gy' % site_dose_i)
        else:
            store_site_doses.append('')
            display_site_doses.append('')

    return display_site_doses + [store_site_doses]
"""--------------------------------------------------------"""

"""***************************** reset dose inputs ***************************"""
@app.callback([Output('heart-dose-input','value'),
               Output('lung-dose-input','value'),
               Output('esoph-dose-input','value')],
              [Input('entry-method-heart','value'),
               Input('entry-method-lung','value'),
               Input('entry-method-esoph','value')],
              [State('heart-dose-input','n_blur'),State('heart-dose-input','n_submit'),State('heart-dose-input','value'),
              State('lung-dose-input','n_blur'),State('lung-dose-input','n_submit'),State('lung-dose-input','value'),
              State('esoph-dose-input','n_blur'),State('esoph-dose-input','n_submit'),State('esoph-dose-input','value')]
)
def reset_dose_input(h_type,l_type,e_type,h_nblur,h_nsubmit,h_cur_val,l_nblur,l_nsubmit,l_cur_val,e_nblur,e_nsubmit,e_cur_val):
    entry_methods  = [h_type,l_type,e_type]
    nsubmits = [h_nsubmit,l_nsubmit,e_nsubmit]
    nblurs = [h_nblur,l_nblur,e_nblur]
    cur_vals = [h_cur_val,l_cur_val,e_cur_val]

    output_vals = []
    for i in range(3):
        if  (entry_methods[i] != 'manual') and ((nblurs[i]>0) or (nsubmits[i]>0)):
            new_val = ''
            output_vals.append(new_val)
        else:
            output_vals.append(cur_vals[i])

    return output_vals
"""----------------------------------------------------------------------"""

""" *************** update disable property of dose inputs ***************** """

@app.callback([
               Output('heart_input-container','hidden'),
               Output('lung_input-container','hidden'),
               Output('esoph_input-container','hidden'),
               Output('heart-default-percent','hidden'),
               Output('lung-default-percent','hidden'),
               Output('esoph-default-percent','hidden'),
               Output('heart-default-percent','children'),
               Output('lung-default-percent','children'),
               Output('esoph-default-percent','children')],
              [Input('entry-method-heart','value'),
               Input('entry-method-lung','value'),
               Input('entry-method-esoph','value')]
)

def update_dose_input_editable(h_entry_method,l_entry_method,e_entry_method):
    entry_methods = [h_entry_method,l_entry_method,e_entry_method]
    inputs_disabled = []
    inputs_hidden = []
    auto_pp_hidden = []
    auto_pp_value = []
    defaults = ['7.0','11.8','25.2']
    for i in range(len(entry_methods)):
        if (entry_methods[i] == 'auto'):
            inputs_hidden.append(True)
            auto_pp_hidden.append(False)
            auto_pp_value.append('{}%'.format(defaults[i]))
        elif entry_methods[i] == 'manual':
            inputs_hidden.append(False)
            auto_pp_hidden.append(True)
            auto_pp_value.append('')
        elif entry_methods[i] == '':
            inputs_hidden.append(True)
            auto_pp_hidden.append(True)
            auto_pp_value.append('')

    return  inputs_hidden + auto_pp_hidden + auto_pp_value
"""----------------------------------------------------------------------"""

""" ********************* entry method-apply to all sites ********************* """

@app.callback(
    [Output('entry-method-heart','value'),
     Output('entry-method-lung','value'),
     Output('entry-method-esoph','value')],
    [Input('apply-all','value')],
    [State('entry-method-heart','value'),
     State('entry-method-lung','value'),
     State('entry-method-esoph','value')]
)

def update_all_dose_radios(entry_method_all,hval,lval,eval):

    if (entry_method_all == 'auto'):
        val = 'auto'
    elif (entry_method_all == 'manual'):
        val = 'manual'
    elif (entry_method_all=='cleared'):
        val = ''
    elif (entry_method_all=='initial'):
        raise PreventUpdate

    return [val]*3
"""----------------------------------------------------------------------"""

@app.callback([Output('apply-all','value'),
               Output('memory-store','data')],
              [
              Input('clear-button','n_clicks'),
              Input('entry-method-heart','value'),
              Input('entry-method-lung','value'),
              Input('entry-method-esoph','value')
              ],
              [State('apply-all','value'),
               State('memory-store','data')]
)
def  update_apply_all(nclicks,hval,lval,eval,all_val,data):

    if data==None:
        data={'clicks':0}
    else:
        data=data

    button_trigger=False
    if (nclicks!=None) and (data!=None):
        if nclicks>data['clicks']:
            button_trigger=True
            data['clicks']=data['clicks']+1
            new_all_val='cleared'

            return new_all_val, data

    all_on = True if ((all_val=='auto') or (all_val=='manual')) else False

    current_vals=np.array([hval,lval,eval])
    if (button_trigger==False) and (all_on==True):
        if np.all(current_vals!='') and np.all(current_vals!=None):
            if np.any(current_vals!=all_val):
                new_all_val='initial'
                return new_all_val, data

    if button_trigger==False & all_on == False:
        raise PreventUpdate

"""----------------------------------------------------------------"""

if __name__ == '__main__':
    app.run_server(debug=True)

max_stars_repo_name: gsedor/nsclc-model

max_issues_repo_name: gsedor/nsclc-model

max_forks_repo_name: gsedor/nsclc-model


text:
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import skew, kurtosis

#generte random values
vals= np. random. normal(0,0.5, 10000)

print(np.percentile(vals, 50))
print(np.percentile(vals, 90))
print(np.percentile(vals, 20))

#first moment
print(" The first moment" , np.mean(vals))

#second moment
print( "The second moment", np.var(vals))

#3rd moment
#used to show how values have skewed to towards
#one tail if the skew is to the right then that is negative skewness
#if the skewness is to the left , then that is positive skewness
print( "The third moment", skew(vals))

#it is used to get the high of the distribution.
print("The fourth moment", kurtosis(vals))

ext: py

reason: import numpy,from scipy


text:
"""
@created_at 2017-08-13
@author Exequiel Fuentes Lettura <efulet@gmail.com>
"""

import os
import numpy as np

import matplotlib.pyplot as plt

import astropy.coordinates as coord
import astropy.units as u

from astropy.io import ascii

from archive_obs.base.util import parse_input_date
from archive_obs.conf import settings

def show(save=False):
    """
    A very simple method for reading a CSV and plot RA and Dec
    """
    path = os.path.join(settings.CSV_PATH, 'fits_headers.csv')
    data = ascii.read(path, format='csv', fast_reader=False, data_start=1, delimiter='\t')

    org = 0
    ra = coord.Angle(data['ra'].filled(np.nan) * u.degree)
    ra = ra.wrap_at(180 * u.degree)
    dec = coord.Angle(data['dec'].filled(np.nan) * u.degree)

    x = np.remainder(data['ra'] + 360 - org, 360)  # shift RA values
    ind = x > 180
    x[ind] -= 360  # scale conversion to [-180, 180]
    x = -x  # reverse the scale: East to the left
    tick_labels = np.array([150, 120, 90, 60, 30, 0, 330, 300, 270, 240, 210])
    tick_labels = np.remainder(tick_labels + 360 + org, 360)

    fig = plt.figure(figsize=(8, 6))
    ax = fig.add_subplot(111, projection='mollweide')
    ax.scatter(ra.radian, dec.radian)
    ax.set_xticklabels(tick_labels)
    ax.grid(True)

    min_date = min(data['date'])
    try:
        min_date = parse_input_date(min_date).strftime('%Y-%m-%d')
    except:
        pass
    max_date = max(data['date'])
    try:
        max_date = parse_input_date(max_date).strftime('%Y-%m-%d')
    except:
        pass

    ax.set_title('Observations made by Gemini at dates: {},{}'.format(min_date, max_date))
    ax.title.set_fontsize(15)
    ax.set_xlabel('RA')
    ax.xaxis.label.set_fontsize(12)
    ax.set_ylabel('Dec')
    ax.yaxis.label.set_fontsize(12)

    plt.show()

    if save:
        fig.savefig(os.path.join(settings.DATA_PATH, 'map.pdf'))

reason: import numpy,import astropy,from astropy


text:
import re

import numpy as np
import pandas as pd

from tfac.dataImport import form_tensor
from tfac.figures.common import getSetup
from tfac.predict import predict_known
from tensorpack import perform_CMTF

def makeFigure():
    predicted = get_predictions()
    accuracies = get_accuracy_by_dtype(predicted)
    fig = plot_results(accuracies)

    return fig

def get_accuracy_by_dtype(predictions):
    """
    Calculates model accuracy w/r to available data types.

    Parameters:
        predictions (pandas.Series): model predictions for each sample

    Returns:
        accuracies (pandas.Series): model accuracy w/r to data types
    """
    d_types = list(set(predictions['Type']))
    accuracies = pd.Series(
        index=d_types,
        dtype=float
    )

    for d_type in d_types:
        predicted = predictions.loc[predictions['Type'] == d_type]
        accuracies.loc[d_type] = np.mean(predicted['Correct'])

    accuracies = accuracies.dropna()
    accuracies = accuracies.sort_values(ascending=False)
    return accuracies

def plot_results(accuracies):
    """
    Plots model accuracy relative to different data sources.

    Parameters:
        accuracies (pandas.Series): model accuracy w/r to data types

    Returns:
        fig (matplotlib.Figure): bar plot depicting model accuracy w/r to data
            types
    """
    fig_size = (4, 4)
    layout = {
        'ncols': 1,
        'nrows': 1
    }
    axs, fig, _ = getSetup(
        fig_size,
        layout
    )

    labels = [re.sub(r'\d', '/', d_type[1:]) for d_type in accuracies.index]
    axs[0].bar(range(len(accuracies)), accuracies)
    axs[0].set_xticks(range(len(accuracies)))
    axs[0].set_xticklabels(labels, rotation=45, ha='right', fontsize=10)
    axs[0].set_ylabel('Mean Accuracy', fontsize=12)
    axs[0].set_xlabel('Datatypes Available', fontsize=12)

    return fig

def get_predictions():
    """
    Predicts samples with known outcomes via cross-validation.

    Parameters:
        None

    Returns:
        predictions (pandas.Series): model predictions for each sample
    """
    tensor, matrix, patient_data = form_tensor()
    patient_data = patient_data.loc[:, ['status', 'type']]

    components = perform_CMTF(tensor, matrix)
    components = components[1][0]

    data = pd.DataFrame(
        components,
        index=patient_data.index,
        columns=list(range(1, components.shape[1] + 1))
    )
    predicted = pd.DataFrame(
        index=patient_data.index
    )
    labels = patient_data.loc[:, 'status']

    predictions = predict_known(data, labels)
    for i in predictions.index:
        if predictions.loc[i] == labels.loc[i]:
            predicted.loc[i, 'Correct'] = 1
        else:
            predicted.loc[i, 'Correct'] = 0

    predicted.loc[:, 'Type'] = patient_data.loc[:, 'type']
    return predicted

max_stars_repo_name: meyer-lab/tfac-mrsa

max_issues_repo_name: meyer-lab/tfac-mrsa

max_forks_repo_name: meyer-lab/tfac-mrsa

reason: import numpy


text:
# Licensed under a 3-clause BSD style license - see LICENSE.rst

# TEST_UNICODE_LITERALS

from __future__ import (absolute_import, division, print_function,
                        unicode_literals)

import pytest
import numpy as np
from numpy.testing import assert_array_equal, assert_array_almost_equal

from ...nduncertainty import (StdDevUncertainty, UnknownUncertainty,
                              IncompatibleUncertaintiesException)
from ... import NDDataRef
from ...nddata import NDData

from ....units import UnitsError, Quantity
from .... import units as u

# Alias NDDataAllMixins in case this will be renamed ... :-)
NDDataArithmetic = NDDataRef

class StdDevUncertaintyUncorrelated(StdDevUncertainty):
    @property
    def supports_correlated(self):
        return False

# Test with Data covers:
# scalars, 1D, 2D and 3D
# broadcasting between them
@pytest.mark.parametrize(('data1', 'data2'), [
                         (np.array(5), np.array(10)),
                         (np.array(5), np.arange(10)),
                         (np.array(5), np.arange(10).reshape(2, 5)),
                         (np.arange(10), np.ones(10) * 2),
                         (np.arange(10), np.ones((10, 10)) * 2),
                         (np.arange(10).reshape(2, 5), np.ones((2, 5)) * 3),
                         (np.arange(1000).reshape(20, 5, 10),
                          np.ones((20, 5, 10)) * 3)
                         ])
def test_arithmetics_data(data1, data2):

    nd1 = NDDataArithmetic(data1)
    nd2 = NDDataArithmetic(data2)

    # Addition
    nd3 = nd1.add(nd2)
    assert_array_equal(data1+data2, nd3.data)
    # Subtraction
    nd4 = nd1.subtract(nd2)
    assert_array_equal(data1-data2, nd4.data)
    # Multiplication
    nd5 = nd1.multiply(nd2)
    assert_array_equal(data1*data2, nd5.data)
    # Division
    nd6 = nd1.divide(nd2)
    assert_array_equal(data1/data2, nd6.data)
    for nd in [nd3, nd4, nd5, nd6]:
        # Check that broadcasting worked as expected
        if data1.ndim > data2.ndim:
            assert data1.shape == nd.data.shape
        else:
            assert data2.shape == nd.data.shape
        # Check all other attributes are not set
        assert nd.unit is None
        assert nd.uncertainty is None
        assert nd.mask is None
        assert len(nd.meta) == 0
        assert nd.wcs is None

# Invalid arithmetic operations for data covering:
# not broadcastable data
def test_arithmetics_data_invalid():
    nd1 = NDDataArithmetic([1, 2, 3])
    nd2 = NDDataArithmetic([1, 2])
    with pytest.raises(ValueError):
        nd1.add(nd2)

# Test with Data and unit and covers:
# identical units (even dimensionless unscaled vs. no unit),
# equivalent units (such as meter and kilometer)
# equivalent composite units (such as m/s and km/h)
@pytest.mark.parametrize(('data1', 'data2'), [
    (np.array(5) * u.s, np.array(10) * u.s),
    (np.array(5) * u.s, np.arange(10) * u.h),
    (np.array(5) * u.s, np.arange(10).reshape(2, 5) * u.min),
    (np.arange(10) * u.m / u.s, np.ones(10) * 2 * u.km / u.s),
    (np.arange(10) * u.m / u.s, np.ones((10, 10)) * 2 * u.m / u.h),
    (np.arange(10).reshape(2, 5) * u.m / u.s,
     np.ones((2, 5)) * 3 * u.km / u.h),
    (np.arange(1000).reshape(20, 5, 10),
     np.ones((20, 5, 10)) * 3 * u.dimensionless_unscaled),
    (np.array(5), np.array(10) * u.s / u.h),
    ])
def test_arithmetics_data_unit_identical(data1, data2):

    nd1 = NDDataArithmetic(data1)
    nd2 = NDDataArithmetic(data2)

    # Addition
    nd3 = nd1.add(nd2)
    ref = data1 + data2
    ref_unit, ref_data = ref.unit, ref.value
    assert_array_equal(ref_data, nd3.data)
    assert nd3.unit == ref_unit
    # Subtraction
    nd4 = nd1.subtract(nd2)
    ref = data1 - data2
    ref_unit, ref_data = ref.unit, ref.value
    assert_array_equal(ref_data, nd4.data)
    assert nd4.unit == ref_unit
    # Multiplication
    nd5 = nd1.multiply(nd2)
    ref = data1 * data2
    ref_unit, ref_data = ref.unit, ref.value
    assert_array_equal(ref_data, nd5.data)
    assert nd5.unit == ref_unit
    # Division
    nd6 = nd1.divide(nd2)
    ref = data1 / data2
    ref_unit, ref_data = ref.unit, ref.value
    assert_array_equal(ref_data, nd6.data)
    assert nd6.unit == ref_unit
    for nd in [nd3, nd4, nd5, nd6]:
        # Check that broadcasting worked as expected
        if data1.ndim > data2.ndim:
            assert data1.shape == nd.data.shape
        else:
            assert data2.shape == nd.data.shape
        # Check all other attributes are not set
        assert nd.uncertainty is None
        assert nd.mask is None
        assert len(nd.meta) == 0
        assert nd.wcs is None

# Test with Data and unit and covers:
# not identical not convertible units
# one with unit (which is not dimensionless) and one without
@pytest.mark.parametrize(('data1', 'data2'), [
    (np.array(5) * u.s, np.array(10) * u.m),
    (np.array(5) * u.Mpc, np.array(10) * u.km / u.s),
    (np.array(5) * u.Mpc, np.array(10)),
    (np.array(5), np.array(10) * u.s),
    ])
def test_arithmetics_data_unit_not_identical(data1, data2):

    nd1 = NDDataArithmetic(data1)
    nd2 = NDDataArithmetic(data2)

    # Addition should not be possible
    with pytest.raises(UnitsError):
        nd1.add(nd2)
    # Subtraction should not be possible
    with pytest.raises(UnitsError):
        nd1.subtract(nd2)
    # Multiplication is possible
    nd3 = nd1.multiply(nd2)
    ref = data1 * data2
    ref_unit, ref_data = ref.unit, ref.value
    assert_array_equal(ref_data, nd3.data)
    assert nd3.unit == ref_unit
    # Division is possible
    nd4 = nd1.divide(nd2)
    ref = data1 / data2
    ref_unit, ref_data = ref.unit, ref.value
    assert_array_equal(ref_data, nd4.data)
    assert nd4.unit == ref_unit
    for nd in [nd3, nd4]:
        # Check all other attributes are not set
        assert nd.uncertainty is None
        assert nd.mask is None
        assert len(nd.meta) == 0
        assert nd.wcs is None

# Tests with wcs (not very sensible because there is no operation between them
# covering:
# both set and identical/not identical
# one set
# None set
@pytest.mark.parametrize(('wcs1', 'wcs2'), [
    (None, None),
    (None, 5),
    (5, None),
    (5, 5),
    (7, 5),
    ])
def test_arithmetics_data_wcs(wcs1, wcs2):

    nd1 = NDDataArithmetic(1, wcs=wcs1)
    nd2 = NDDataArithmetic(1, wcs=wcs2)

    if wcs1 is None and wcs2 is None:
        ref_wcs = None
    elif wcs1 is None:
        ref_wcs = wcs2
    elif wcs2 is None:
        ref_wcs = wcs1
    else:
        ref_wcs = wcs1

    # Addition
    nd3 = nd1.add(nd2)
    assert ref_wcs == nd3.wcs
    # Subtraction
    nd4 = nd1.subtract(nd2)
    assert ref_wcs == nd3.wcs
    # Multiplication
    nd5 = nd1.multiply(nd2)
    assert ref_wcs == nd3.wcs
    # Division
    nd6 = nd1.divide(nd2)
    assert ref_wcs == nd3.wcs
    for nd in [nd3, nd4, nd5, nd6]:
        # Check all other attributes are not set
        assert nd.unit is None
        assert nd.uncertainty is None
        assert len(nd.meta) == 0
        assert nd.mask is None

# Masks are completely separated in the NDArithmetics from the data so we need
# no correlated tests but covering:
# masks 1D, 2D and mixed cases with broadcasting
@pytest.mark.parametrize(('mask1', 'mask2'), [
    (None, None),
    (None, False),
    (True, None),
    (False, False),
    (True, False),
    (False, True),
    (True, True),
    (np.array(False), np.array(True)),
    (np.array(False), np.array([0, 1, 0, 1, 1], dtype=np.bool_)),
    (np.array(True),
     np.array([[0, 1, 0, 1, 1], [1, 1, 0, 1, 1]], dtype=np.bool_)),
    (np.array([0, 1, 0, 1, 1], dtype=np.bool_),
     np.array([1, 1, 0, 0, 1], dtype=np.bool_)),
    (np.array([0, 1, 0, 1, 1], dtype=np.bool_),
     np.array([[0, 1, 0, 1, 1], [1, 0, 0, 1, 1]], dtype=np.bool_)),
    (np.array([[0, 1, 0, 1, 1], [1, 0, 0, 1, 1]], dtype=np.bool_),
     np.array([[0, 1, 0, 1, 1], [1, 1, 0, 1, 1]], dtype=np.bool_)),
    ])
def test_arithmetics_data_masks(mask1, mask2):

    nd1 = NDDataArithmetic(1, mask=mask1)
    nd2 = NDDataArithmetic(1, mask=mask2)

    if mask1 is None and mask2 is None:
        ref_mask = None
    elif mask1 is None:
        ref_mask = mask2
    elif mask2 is None:
        ref_mask = mask1
    else:
        ref_mask = mask1 | mask2

    # Addition
    nd3 = nd1.add(nd2)
    assert_array_equal(ref_mask, nd3.mask)
    # Subtraction
    nd4 = nd1.subtract(nd2)
    assert_array_equal(ref_mask, nd4.mask)
    # Multiplication
    nd5 = nd1.multiply(nd2)
    assert_array_equal(ref_mask, nd5.mask)
    # Division
    nd6 = nd1.divide(nd2)
    assert_array_equal(ref_mask, nd6.mask)
    for nd in [nd3, nd4, nd5, nd6]:
        # Check all other attributes are not set
        assert nd.unit is None
        assert nd.uncertainty is None
        assert len(nd.meta) == 0
        assert nd.wcs is None

# One additional case which can not be easily incorporated in the test above
# what happens if the masks are numpy ndarrays are not broadcastable
def test_arithmetics_data_masks_invalid():

    nd1 = NDDataArithmetic(1, mask=np.array([1, 0], dtype=np.bool_))
    nd2 = NDDataArithmetic(1, mask=np.array([1, 0, 1], dtype=np.bool_))

    with pytest.raises(ValueError):
        nd1.add(nd2)
    with pytest.raises(ValueError):
        nd1.multiply(nd2)
    with pytest.raises(ValueError):
        nd1.subtract(nd2)
    with pytest.raises(ValueError):
        nd1.divide(nd2)

# Covering:
# both have uncertainties (data and uncertainty without unit)
# tested against manually determined resulting uncertainties to verify the
# implemented formulas
# this test only works as long as data1 and data2 do not contain any 0
def test_arithmetics_stddevuncertainty_basic():
    nd1 = NDDataArithmetic([1, 2, 3], uncertainty=StdDevUncertainty([1, 1, 3]))
    nd2 = NDDataArithmetic([2, 2, 2], uncertainty=StdDevUncertainty([2, 2, 2]))
    nd3 = nd1.add(nd2)
    nd4 = nd2.add(nd1)
    # Inverse operation should result in the same uncertainty
    assert_array_equal(nd3.uncertainty.array, nd4.uncertainty.array)
    # Compare it to the theoretical uncertainty
    ref_uncertainty = np.sqrt(np.array([1, 1, 3])**2 + np.array([2, 2, 2])**2)
    assert_array_equal(nd3.uncertainty.array, ref_uncertainty)

    nd3 = nd1.subtract(nd2)
    nd4 = nd2.subtract(nd1)
    # Inverse operation should result in the same uncertainty
    assert_array_equal(nd3.uncertainty.array, nd4.uncertainty.array)
    # Compare it to the theoretical uncertainty (same as for add)
    assert_array_equal(nd3.uncertainty.array, ref_uncertainty)

    # Multiplication and Division only work with almost equal array comparisons
    # since the formula implemented and the formula used as reference are
    # slightly different.
    nd3 = nd1.multiply(nd2)
    nd4 = nd2.multiply(nd1)
    # Inverse operation should result in the same uncertainty
    assert_array_almost_equal(nd3.uncertainty.array, nd4.uncertainty.array)
    # Compare it to the theoretical uncertainty
    ref_uncertainty = np.abs(np.array([2, 4, 6])) * np.sqrt(
        (np.array([1, 1, 3]) / np.array([1, 2, 3]))**2 +
        (np.array([2, 2, 2]) / np.array([2, 2, 2]))**2)
    assert_array_almost_equal(nd3.uncertainty.array, ref_uncertainty)

    nd3 = nd1.divide(nd2)
    nd4 = nd2.divide(nd1)
    # Inverse operation gives a different uncertainty!
    # Compare it to the theoretical uncertainty
    ref_uncertainty_1 = np.abs(np.array([1/2, 2/2, 3/2])) * np.sqrt(
        (np.array([1, 1, 3]) / np.array([1, 2, 3]))**2 +
        (np.array([2, 2, 2]) / np.array([2, 2, 2]))**2)
    assert_array_almost_equal(nd3.uncertainty.array, ref_uncertainty_1)
    ref_uncertainty_2 = np.abs(np.array([2, 1, 2/3])) * np.sqrt(
        (np.array([1, 1, 3]) / np.array([1, 2, 3]))**2 +
        (np.array([2, 2, 2]) / np.array([2, 2, 2]))**2)
    assert_array_almost_equal(nd4.uncertainty.array, ref_uncertainty_2)

# Tests for correlation, covering
# correlation between -1 and 1 with correlation term being positive / negative
# also with one data being once positive and once completely negative
# The point of this test is to compare the used formula to the theoretical one.
# TODO: Maybe covering units too but I think that should work because of
# the next tests. Also this may be reduced somehow.
@pytest.mark.parametrize(('cor', 'uncert1', 'data2'), [
    (-1, [1, 1, 3], [2, 2, 7]),
    (-0.5, [1, 1, 3], [2, 2, 7]),
    (-0.25, [1, 1, 3], [2, 2, 7]),
    (0, [1, 1, 3], [2, 2, 7]),
    (0.25, [1, 1, 3], [2, 2, 7]),
    (0.5, [1, 1, 3], [2, 2, 7]),
    (1, [1, 1, 3], [2, 2, 7]),
    (-1, [-1, -1, -3], [2, 2, 7]),
    (-0.5, [-1, -1, -3], [2, 2, 7]),
    (-0.25, [-1, -1, -3], [2, 2, 7]),
    (0, [-1, -1, -3], [2, 2, 7]),
    (0.25, [-1, -1, -3], [2, 2, 7]),
    (0.5, [-1, -1, -3], [2, 2, 7]),
    (1, [-1, -1, -3], [2, 2, 7]),
    (-1, [1, 1, 3], [-2, -3, -2]),
    (-0.5, [1, 1, 3], [-2, -3, -2]),
    (-0.25, [1, 1, 3], [-2, -3, -2]),
    (0, [1, 1, 3], [-2, -3, -2]),
    (0.25, [1, 1, 3], [-2, -3, -2]),
    (0.5, [1, 1, 3], [-2, -3, -2]),
    (1, [1, 1, 3], [-2, -3, -2]),
    (-1, [-1, -1, -3], [-2, -3, -2]),
    (-0.5, [-1, -1, -3], [-2, -3, -2]),
    (-0.25, [-1, -1, -3], [-2, -3, -2]),
    (0, [-1, -1, -3], [-2, -3, -2]),
    (0.25, [-1, -1, -3], [-2, -3, -2]),
    (0.5, [-1, -1, -3], [-2, -3, -2]),
    (1, [-1, -1, -3], [-2, -3, -2]),
    ])
def test_arithmetics_stddevuncertainty_basic_with_correlation(
        cor, uncert1, data2):
    data1 = np.array([1, 2, 3])
    data2 = np.array(data2)
    uncert1 = np.array(uncert1)
    uncert2 = np.array([2, 2, 2])
    nd1 = NDDataArithmetic(data1, uncertainty=StdDevUncertainty(uncert1))
    nd2 = NDDataArithmetic(data2, uncertainty=StdDevUncertainty(uncert2))
    nd3 = nd1.add(nd2, uncertainty_correlation=cor)
    nd4 = nd2.add(nd1, uncertainty_correlation=cor)
    # Inverse operation should result in the same uncertainty
    assert_array_equal(nd3.uncertainty.array, nd4.uncertainty.array)
    # Compare it to the theoretical uncertainty
    ref_uncertainty = np.sqrt(uncert1**2 + uncert2**2 +
                              2 * cor * uncert1 * uncert2)
    assert_array_equal(nd3.uncertainty.array, ref_uncertainty)

    nd3 = nd1.subtract(nd2, uncertainty_correlation=cor)
    nd4 = nd2.subtract(nd1, uncertainty_correlation=cor)
    # Inverse operation should result in the same uncertainty
    assert_array_equal(nd3.uncertainty.array, nd4.uncertainty.array)
    # Compare it to the theoretical uncertainty
    ref_uncertainty = np.sqrt(uncert1**2 + uncert2**2 -
                              2 * cor * uncert1 * uncert2)
    assert_array_equal(nd3.uncertainty.array, ref_uncertainty)

    # Multiplication and Division only work with almost equal array comparisons
    # since the formula implemented and the formula used as reference are
    # slightly different.
    nd3 = nd1.multiply(nd2, uncertainty_correlation=cor)
    nd4 = nd2.multiply(nd1, uncertainty_correlation=cor)
    # Inverse operation should result in the same uncertainty
    assert_array_almost_equal(nd3.uncertainty.array, nd4.uncertainty.array)
    # Compare it to the theoretical uncertainty
    ref_uncertainty = (np.abs(data1 * data2)) * np.sqrt(
        (uncert1 / data1)**2 + (uncert2 / data2)**2 +
        (2 * cor * uncert1 * uncert2 / (data1 * data2)))
    assert_array_almost_equal(nd3.uncertainty.array, ref_uncertainty)

    nd3 = nd1.divide(nd2, uncertainty_correlation=cor)
    nd4 = nd2.divide(nd1, uncertainty_correlation=cor)
    # Inverse operation gives a different uncertainty!
    # Compare it to the theoretical uncertainty
    ref_uncertainty_1 = (np.abs(data1 / data2)) * np.sqrt(
        (uncert1 / data1)**2 + (uncert2 / data2)**2 -
        (2 * cor * uncert1 * uncert2 / (data1 * data2)))
    assert_array_almost_equal(nd3.uncertainty.array, ref_uncertainty_1)
    ref_uncertainty_2 = (np.abs(data2 / data1)) * np.sqrt(
        (uncert1 / data1)**2 + (uncert2 / data2)**2 -
        (2 * cor * uncert1 * uncert2 / (data1 * data2)))
    assert_array_almost_equal(nd4.uncertainty.array, ref_uncertainty_2)

# Covering:
# just an example that a np.ndarray works as correlation, no checks for
# the right result since these were basically done in the function above.
def test_arithmetics_stddevuncertainty_basic_with_correlation_array():
    data1 = np.array([1, 2, 3])
    data2 = np.array([1, 1, 1])
    uncert1 = np.array([1, 1, 1])
    uncert2 = np.array([2, 2, 2])
    cor = np.array([0, 0.25, 0])
    nd1 = NDDataArithmetic(data1, uncertainty=StdDevUncertainty(uncert1))
    nd2 = NDDataArithmetic(data2, uncertainty=StdDevUncertainty(uncert2))
    nd1.add(nd2, uncertainty_correlation=cor)

# Covering:
# That propagate throws an exception when correlation is given but the
# uncertainty does not support correlation.
def test_arithmetics_with_correlation_unsupported():
    data1 = np.array([1, 2, 3])
    data2 = np.array([1, 1, 1])
    uncert1 = np.array([1, 1, 1])
    uncert2 = np.array([2, 2, 2])
    cor = 3
    nd1 = NDDataArithmetic(data1,
                           uncertainty=StdDevUncertaintyUncorrelated(uncert1))
    nd2 = NDDataArithmetic(data2,
                           uncertainty=StdDevUncertaintyUncorrelated(uncert2))

    with pytest.raises(ValueError):
        nd1.add(nd2, uncertainty_correlation=cor)

# Covering:
# only one has an uncertainty (data and uncertainty without unit)
# tested against the case where the other one has zero uncertainty. (this case
# must be correct because we tested it in the last case)
# Also verify that if the result of the data has negative values the resulting
# uncertainty has no negative values.
def test_arithmetics_stddevuncertainty_one_missing():
    nd1 = NDDataArithmetic([1, -2, 3])
    nd1_ref = NDDataArithmetic([1, -2, 3],
                               uncertainty=StdDevUncertainty([0, 0, 0]))
    nd2 = NDDataArithmetic([2, 2, -2],
                           uncertainty=StdDevUncertainty([2, 2, 2]))

    # Addition
    nd3 = nd1.add(nd2)
    nd3_ref = nd1_ref.add(nd2)
    assert_array_equal(nd3.uncertainty.array, nd3_ref.uncertainty.array)
    assert_array_equal(np.abs(nd3.uncertainty.array), nd3.uncertainty.array)

    nd3 = nd2.add(nd1)
    nd3_ref = nd2.add(nd1_ref)
    assert_array_equal(nd3.uncertainty.array, nd3_ref.uncertainty.array)
    assert_array_equal(np.abs(nd3.uncertainty.array), nd3.uncertainty.array)

    # Subtraction
    nd3 = nd1.subtract(nd2)
    nd3_ref = nd1_ref.subtract(nd2)
    assert_array_equal(nd3.uncertainty.array, nd3_ref.uncertainty.array)
    assert_array_equal(np.abs(nd3.uncertainty.array), nd3.uncertainty.array)

    nd3 = nd2.subtract(nd1)
    nd3_ref = nd2.subtract(nd1_ref)
    assert_array_equal(nd3.uncertainty.array, nd3_ref.uncertainty.array)
    assert_array_equal(np.abs(nd3.uncertainty.array), nd3.uncertainty.array)

    # Multiplication
    nd3 = nd1.multiply(nd2)
    nd3_ref = nd1_ref.multiply(nd2)
    assert_array_equal(nd3.uncertainty.array, nd3_ref.uncertainty.array)
    assert_array_equal(np.abs(nd3.uncertainty.array), nd3.uncertainty.array)

    nd3 = nd2.multiply(nd1)
    nd3_ref = nd2.multiply(nd1_ref)
    assert_array_equal(nd3.uncertainty.array, nd3_ref.uncertainty.array)
    assert_array_equal(np.abs(nd3.uncertainty.array), nd3.uncertainty.array)

    # Division
    nd3 = nd1.divide(nd2)
    nd3_ref = nd1_ref.divide(nd2)
    assert_array_equal(nd3.uncertainty.array, nd3_ref.uncertainty.array)
    assert_array_equal(np.abs(nd3.uncertainty.array), nd3.uncertainty.array)

    nd3 = nd2.divide(nd1)
    nd3_ref = nd2.divide(nd1_ref)
    assert_array_equal(nd3.uncertainty.array, nd3_ref.uncertainty.array)
    assert_array_equal(np.abs(nd3.uncertainty.array), nd3.uncertainty.array)

# Covering:
# data with unit and uncertainty with unit (but equivalent units)
# compared against correctly scaled NDDatas
@pytest.mark.parametrize(('uncert1', 'uncert2'), [
    (np.array([1, 2, 3]) * u.m, None),
    (np.array([1, 2, 3]) * u.cm, None),
    (None, np.array([1, 2, 3]) * u.m),
    (None, np.array([1, 2, 3]) * u.cm),
    (np.array([1, 2, 3]), np.array([2, 3, 4])),
    (np.array([1, 2, 3]) * u.m, np.array([2, 3, 4])),
    (np.array([1, 2, 3]), np.array([2, 3, 4])) * u.m,
    (np.array([1, 2, 3]) * u.m, np.array([2, 3, 4])) * u.m,
    (np.array([1, 2, 3]) * u.cm, np.array([2, 3, 4])),
    (np.array([1, 2, 3]), np.array([2, 3, 4])) * u.cm,
    (np.array([1, 2, 3]) * u.cm, np.array([2, 3, 4])) * u.cm,
    (np.array([1, 2, 3]) * u.km, np.array([2, 3, 4])) * u.cm,
    ])
def test_arithmetics_stddevuncertainty_with_units(uncert1, uncert2):
    # Data has same units
    data1 = np.array([1, 2, 3]) * u.m
    data2 = np.array([-4, 7, 0]) * u.m
    if uncert1 is not None:
        uncert1 = StdDevUncertainty(uncert1)
        if isinstance(uncert1, Quantity):
            uncert1_ref = uncert1.to_value(data1.unit)
        else:
            uncert1_ref = uncert1
        uncert_ref1 = StdDevUncertainty(uncert1_ref, copy=True)
    else:
        uncert1 = None
        uncert_ref1 = None

    if uncert2 is not None:
        uncert2 = StdDevUncertainty(uncert2)
        if isinstance(uncert2, Quantity):
            uncert2_ref = uncert2.to_value(data2.unit)
        else:
            uncert2_ref = uncert2
        uncert_ref2 = StdDevUncertainty(uncert2_ref, copy=True)
    else:
        uncert2 = None
        uncert_ref2 = None

    nd1 = NDDataArithmetic(data1, uncertainty=uncert1)
    nd2 = NDDataArithmetic(data2, uncertainty=uncert2)

    nd1_ref = NDDataArithmetic(data1, uncertainty=uncert_ref1)
    nd2_ref = NDDataArithmetic(data2, uncertainty=uncert_ref2)

    # Let's start the tests
    # Addition
    nd3 = nd1.add(nd2)
    nd3_ref = nd1_ref.add(nd2_ref)
    assert nd3.unit == nd3_ref.unit
    assert nd3.uncertainty.unit == nd3_ref.uncertainty.unit
    assert_array_equal(nd3.uncertainty.array, nd3.uncertainty.array)

    nd3 = nd2.add(nd1)
    nd3_ref = nd2_ref.add(nd1_ref)
    assert nd3.unit == nd3_ref.unit
    assert nd3.uncertainty.unit == nd3_ref.uncertainty.unit
    assert_array_equal(nd3.uncertainty.array, nd3.uncertainty.array)

    # Subtraction
    nd3 = nd1.subtract(nd2)
    nd3_ref = nd1_ref.subtract(nd2_ref)
    assert nd3.unit == nd3_ref.unit
    assert nd3.uncertainty.unit == nd3_ref.uncertainty.unit
    assert_array_equal(nd3.uncertainty.array, nd3.uncertainty.array)

    nd3 = nd2.subtract(nd1)
    nd3_ref = nd2_ref.subtract(nd1_ref)
    assert nd3.unit == nd3_ref.unit
    assert nd3.uncertainty.unit == nd3_ref.uncertainty.unit
    assert_array_equal(nd3.uncertainty.array, nd3.uncertainty.array)

    # Multiplication
    nd3 = nd1.multiply(nd2)
    nd3_ref = nd1_ref.multiply(nd2_ref)
    assert nd3.unit == nd3_ref.unit
    assert nd3.uncertainty.unit == nd3_ref.uncertainty.unit
    assert_array_equal(nd3.uncertainty.array, nd3.uncertainty.array)

    nd3 = nd2.multiply(nd1)
    nd3_ref = nd2_ref.multiply(nd1_ref)
    assert nd3.unit == nd3_ref.unit
    assert nd3.uncertainty.unit == nd3_ref.uncertainty.unit
    assert_array_equal(nd3.uncertainty.array, nd3.uncertainty.array)

    # Division
    nd3 = nd1.divide(nd2)
    nd3_ref = nd1_ref.divide(nd2_ref)
    assert nd3.unit == nd3_ref.unit
    assert nd3.uncertainty.unit == nd3_ref.uncertainty.unit
    assert_array_equal(nd3.uncertainty.array, nd3.uncertainty.array)

    nd3 = nd2.divide(nd1)
    nd3_ref = nd2_ref.divide(nd1_ref)
    assert nd3.unit == nd3_ref.unit
    assert nd3.uncertainty.unit == nd3_ref.uncertainty.unit
    assert_array_equal(nd3.uncertainty.array, nd3.uncertainty.array)

# Test abbreviation and long name for taking the first found meta, mask, wcs
@pytest.mark.parametrize(('use_abbreviation'), ['ff', 'first_found'])
def test_arithmetics_handle_switches(use_abbreviation):
    meta1 = {'a': 1}
    meta2 = {'b': 2}
    mask1 = True
    mask2 = False
    uncertainty1 = StdDevUncertainty([1, 2, 3])
    uncertainty2 = StdDevUncertainty([1, 2, 3])
    wcs1 = 5
    wcs2 = 100
    data1 = [1, 1, 1]
    data2 = [1, 1, 1]

    nd1 = NDDataArithmetic(data1, meta=meta1, mask=mask1, wcs=wcs1,
                           uncertainty=uncertainty1)
    nd2 = NDDataArithmetic(data2, meta=meta2, mask=mask2, wcs=wcs2,
                           uncertainty=uncertainty2)
    nd3 = NDDataArithmetic(data1)

    # Both have the attributes but option None is chosen
    nd_ = nd1.add(nd2, propagate_uncertainties=None, handle_meta=None,
                  handle_mask=None, compare_wcs=None)
    assert nd_.wcs is None
    assert len(nd_.meta) == 0
    assert nd_.mask is None
    assert nd_.uncertainty is None

    # Only second has attributes and False is chosen
    nd_ = nd3.add(nd2, propagate_uncertainties=False,
                  handle_meta=use_abbreviation, handle_mask=use_abbreviation,
                  compare_wcs=use_abbreviation)
    assert nd_.wcs == wcs2
    assert nd_.meta == meta2
    assert nd_.mask == mask2
    assert_array_equal(nd_.uncertainty.array, uncertainty2.array)

    # Only first has attributes and False is chosen
    nd_ = nd1.add(nd3, propagate_uncertainties=False,
                  handle_meta=use_abbreviation, handle_mask=use_abbreviation,
                  compare_wcs=use_abbreviation)
    assert nd_.wcs == wcs1
    assert nd_.meta == meta1
    assert nd_.mask == mask1
    assert_array_equal(nd_.uncertainty.array, uncertainty1.array)

def test_arithmetics_meta_func():
    def meta_fun_func(meta1, meta2, take='first'):
        if take == 'first':
            return meta1
        else:
            return meta2

    meta1 = {'a': 1}
    meta2 = {'a': 3, 'b': 2}
    mask1 = True
    mask2 = False
    uncertainty1 = StdDevUncertainty([1, 2, 3])
    uncertainty2 = StdDevUncertainty([1, 2, 3])
    wcs1 = 5
    wcs2 = 100
    data1 = [1, 1, 1]
    data2 = [1, 1, 1]

    nd1 = NDDataArithmetic(data1, meta=meta1, mask=mask1, wcs=wcs1,
                           uncertainty=uncertainty1)
    nd2 = NDDataArithmetic(data2, meta=meta2, mask=mask2, wcs=wcs2,
                           uncertainty=uncertainty2)

    nd3 = nd1.add(nd2, handle_meta=meta_fun_func)
    assert nd3.meta['a'] == 1
    assert 'b' not in nd3.meta

    nd4 = nd1.add(nd2, handle_meta=meta_fun_func, meta_take='second')
    assert nd4.meta['a'] == 3
    assert nd4.meta['b'] == 2

    with pytest.raises(KeyError):
        nd1.add(nd2, handle_meta=meta_fun_func, take='second')

def test_arithmetics_wcs_func():
    def wcs_comp_func(wcs1, wcs2, tolerance=0.1):
        if abs(wcs1 - wcs2) <= tolerance:
            return True
        else:
            return False

    meta1 = {'a': 1}
    meta2 = {'a': 3, 'b': 2}
    mask1 = True
    mask2 = False
    uncertainty1 = StdDevUncertainty([1, 2, 3])
    uncertainty2 = StdDevUncertainty([1, 2, 3])
    wcs1 = 99.99
    wcs2 = 100
    data1 = [1, 1, 1]
    data2 = [1, 1, 1]

    nd1 = NDDataArithmetic(data1, meta=meta1, mask=mask1, wcs=wcs1,
                           uncertainty=uncertainty1)
    nd2 = NDDataArithmetic(data2, meta=meta2, mask=mask2, wcs=wcs2,
                           uncertainty=uncertainty2)

    nd3 = nd1.add(nd2, compare_wcs=wcs_comp_func)
    assert nd3.wcs == 99.99

    with pytest.raises(ValueError):
        nd1.add(nd2, compare_wcs=wcs_comp_func, wcs_tolerance=0.00001)

    with pytest.raises(KeyError):
        nd1.add(nd2, compare_wcs=wcs_comp_func, tolerance=1)

def test_arithmetics_mask_func():
    def mask_sad_func(mask1, mask2, fun=0):
        if fun > 0.5:
            return mask2
        else:
            return mask1

    meta1 = {'a': 1}
    meta2 = {'a': 3, 'b': 2}
    mask1 = [True, False, True]
    mask2 = [True, False, False]
    uncertainty1 = StdDevUncertainty([1, 2, 3])
    uncertainty2 = StdDevUncertainty([1, 2, 3])
    wcs1 = 99.99
    wcs2 = 100
    data1 = [1, 1, 1]
    data2 = [1, 1, 1]

    nd1 = NDDataArithmetic(data1, meta=meta1, mask=mask1, wcs=wcs1,
                           uncertainty=uncertainty1)
    nd2 = NDDataArithmetic(data2, meta=meta2, mask=mask2, wcs=wcs2,
                           uncertainty=uncertainty2)

    nd3 = nd1.add(nd2, handle_mask=mask_sad_func)
    assert_array_equal(nd3.mask, nd1.mask)

    nd4 = nd1.add(nd2, handle_mask=mask_sad_func, mask_fun=1)
    assert_array_equal(nd4.mask, nd2.mask)

    with pytest.raises(KeyError):
        nd1.add(nd2, handle_mask=mask_sad_func, fun=1)

@pytest.mark.parametrize('meth', ['add', 'subtract', 'divide', 'multiply'])
def test_two_argument_useage(meth):
    ndd1 = NDDataArithmetic(np.ones((3, 3)))
    ndd2 = NDDataArithmetic(np.ones((3, 3)))

    # Call add on the class (not the instance) and compare it with already
    # tested useage:
    ndd3 = getattr(NDDataArithmetic, meth)(ndd1, ndd2)
    ndd4 = getattr(ndd1, meth)(ndd2)
    np.testing.assert_array_equal(ndd3.data, ndd4.data)

    # And the same done on an unrelated instance...
    ndd3 = getattr(NDDataArithmetic(-100), meth)(ndd1, ndd2)
    ndd4 = getattr(ndd1, meth)(ndd2)
    np.testing.assert_array_equal(ndd3.data, ndd4.data)

@pytest.mark.parametrize('meth', ['add', 'subtract', 'divide', 'multiply'])
def test_two_argument_useage_non_nddata_first_arg(meth):
    data1 = 50
    data2 = 100

    # Call add on the class (not the instance)
    ndd3 = getattr(NDDataArithmetic, meth)(data1, data2)

    # Compare it with the instance-useage and two identical NDData-like
    # classes:
    ndd1 = NDDataArithmetic(data1)
    ndd2 = NDDataArithmetic(data2)
    ndd4 = getattr(ndd1, meth)(ndd2)
    np.testing.assert_array_equal(ndd3.data, ndd4.data)

    # and check it's also working when called on an instance
    ndd3 = getattr(NDDataArithmetic(-100), meth)(data1, data2)
    ndd4 = getattr(ndd1, meth)(ndd2)
    np.testing.assert_array_equal(ndd3.data, ndd4.data)

def test_arithmetics_unknown_uncertainties():
    # Not giving any uncertainty class means it is saved as UnknownUncertainty
    ndd1 = NDDataArithmetic(np.ones((3, 3)),
                            uncertainty=UnknownUncertainty(np.ones((3, 3))))
    ndd2 = NDDataArithmetic(np.ones((3, 3)),
                            uncertainty=UnknownUncertainty(np.ones((3, 3))*2))
    # There is no way to propagate uncertainties:
    with pytest.raises(IncompatibleUncertaintiesException):
        ndd1.add(ndd2)
    # But it should be possible without propagation
    ndd3 = ndd1.add(ndd2, propagate_uncertainties=False)
    np.testing.assert_array_equal(ndd1.uncertainty.array,
                                  ndd3.uncertainty.array)

    ndd4 = ndd1.add(ndd2, propagate_uncertainties=None)
    assert ndd4.uncertainty is None

max_stars_repo_name: REMeyer/astropy

max_issues_repo_name: REMeyer/astropy

max_forks_repo_name: REMeyer/astropy

reason: import numpy,from numpy


text:
from abc import ABC
import numpy as np
from .. import pixel_formats
from ..core import Format

class _GrayBase(Format, ABC):
    """Base for all gray/yuv400 formats."""
    @staticmethod
    def chroma_subsampling():
        return 0, 0

    def unpack(self, data):
        y = data['y']
        return y, None, None

    def pack(self, yuv):
        data = np.empty(yuv[0].shape[0], dtype=self.dtype)
        data['y'][:] = yuv[0]
        return data

class Gray(_GrayBase):
    @staticmethod
    def identifier():
        return "gray"

    @property
    def dtype(self):
        return np.dtype([
            ('y', '<u1', (self._height, self._width)),
            ('u', '<u1', (0, 0)),
            ('v', '<u1', (0, 0))
        ])

class Gray10LE(_GrayBase):
    @staticmethod
    def identifier():
        return "gray10le"

    @property
    def dtype(self):
        return np.dtype([
            ('y', '<u2', (self._height, self._width)),
            ('u', '|V0', (0, 0)),
            ('v', '|V0', (0, 0))
        ])

class Gray10BE(_GrayBase):
    @staticmethod
    def identifier():
        return "gray10be"

    @property
    def dtype(self):
        return np.dtype([
            ('y', '>u2', (self._height, self._width)),
            ('u', '|V0', (0, 0)),
            ('v', '|V0', (0, 0))
        ])

class Gray16LE(_GrayBase):
    @staticmethod
    def identifier():
        return "gray16le"

    @property
    def dtype(self):
        return np.dtype([
            ('y', '<u2', (self._height, self._width)),
            ('u', '|V0', (0, 0)),
            ('v', '|V0', (0, 0))
        ])

class Gray16BE(_GrayBase):
    @staticmethod
    def identifier():
        return "gray16be"

    @property
    def dtype(self):
        return np.dtype([
            ('y', '>u2', (self._height, self._width)),
            ('u', '|V0', (0, 0)),
            ('v', '|V0', (0, 0))
        ])

class Gray9LE(_GrayBase):
    @staticmethod
    def identifier():
        return "gray9le"

    @property
    def dtype(self):
        return np.dtype([
            ('y', '<u2', (self._height, self._width)),
            ('u', '|V0', (0, 0)),
            ('v', '|V0', (0, 0))
        ])

class Gray9BE(_GrayBase):
    @staticmethod
    def identifier():
        return "gray9be"

    @property
    def dtype(self):
        return np.dtype([
            ('y', '>u2', (self._height, self._width)),
            ('u', '|V0', (0, 0)),
            ('v', '|V0', (0, 0))
        ])

class Gray12LE(_GrayBase):
    @staticmethod
    def identifier():
        return "gray12le"

    @property
    def dtype(self):
        return np.dtype([
            ('y', '<u2', (self._height, self._width)),
            ('u', '|V0', (0, 0)),
            ('v', '|V0', (0, 0))
        ])

class Gray12BE(_GrayBase):
    @staticmethod
    def identifier():
        return "gray12be"

    @property
    def dtype(self):
        return np.dtype([
            ('y', '>u2', (self._height, self._width)),
            ('u', '|V0', (0, 0)),
            ('v', '|V0', (0, 0))
        ])

class Gray14LE(_GrayBase):
    @staticmethod
    def identifier():
        return "gray14le"

    @property
    def dtype(self):
        return np.dtype([
            ('y', '<u2', (self._height, self._width)),
            ('u', '|V0', (0, 0)),
            ('v', '|V0', (0, 0))
        ])

class Gray14BE(_GrayBase):
    @staticmethod
    def identifier():
        return "gray14be"

    @property
    def dtype(self):
        return np.dtype([
            ('y', '>u2', (self._height, self._width)),
            ('u', '|V0', (0, 0)),
            ('v', '|V0', (0, 0))
        ])

pixel_formats.register(Gray)
pixel_formats.register(Gray10LE)
pixel_formats.register(Gray10BE)
pixel_formats.register(Gray16LE)
pixel_formats.register(Gray16BE)
pixel_formats.register(Gray9LE)
pixel_formats.register(Gray9BE)
pixel_formats.register(Gray12LE)
pixel_formats.register(Gray12BE)
pixel_formats.register(Gray14LE)
pixel_formats.register(Gray14BE)

max_stars_repo_name: labradon/yuvio

max_issues_repo_name: labradon/yuvio

max_forks_repo_name: labradon/yuvio


text:
"""Test `IDataLoader` data shuffling by a `DataShuffle` decorator"""

import unittest
import numpy as np

from lstm_ee.data.data_loader.dict_loader  import DictLoader
from lstm_ee.data.data_loader.data_shuffle import DataShuffle

from .tests_data_loader_base import FuncsDataLoaderBase

class TestsDataShuffle(unittest.TestCase, FuncsDataLoaderBase):
    """Test `DataShuffle` decorator"""

    def test_scalar_var(self):
        """Test shuffling of a single scalar variable"""
        seed        = 1223
        data        = { 'var' : np.array([ 1, 2, 3, 4, -1 ]) }
        data_loader = DataShuffle(DictLoader(data), seed)

        indices_original = np.arange(0, len(data['var']))
        indices_shuffled = np.array(indices_original[:])

        np.random.seed(seed)
        np.random.shuffle(indices_shuffled)
        self.assertTrue(
            np.any(~ np.isclose(indices_original, indices_shuffled))
        )

        data_shuffled = { 'var' : data['var'][indices_shuffled] }
        self._compare_scalar_vars(data_shuffled, data_loader, 'var')

    def test_varr_var(self):
        """Test shuffling of a single variable length array variable"""
        seed        = 321
        data        = { 'var' : np.array([[1, 2], [], [3], [4,5,6,7], [-1]]) }
        data_loader = DataShuffle(DictLoader(data), seed)

        indices_original = np.arange(0, len(data['var']))
        indices_shuffled = np.array(indices_original[:])

        np.random.seed(seed)
        np.random.shuffle(indices_shuffled)
        self.assertTrue(
            np.any(~ np.isclose(indices_original, indices_shuffled))
        )

        data_shuffled = { 'var' : data['var'][indices_shuffled] }
        self._compare_varr_vars(data_shuffled, data_loader, 'var')

if __name__ == '__main__':
    unittest.main()


text:
# -*- coding: utf-8 -*-
"""
Created on Mon Aug 23 18:38:23 2021

@author: ben-o_000
"""

import numpy as np

def stumpS(z):

    if z > 0:
        s = (np.sqrt(z) - np.sin(np.sqrt(z)))/(np.sqrt(z))**3
    elif z < 0:
        s = (np.sinh(np.sqrt(-z)) - np.sqrt(-z))/(np.sqrt(-z))**3
    else:
        s = 1/6

    return s

max_stars_repo_name: sgboakes/pysatellite

max_issues_repo_name: sgboakes/pysatellite

max_forks_repo_name: sgboakes/pysatellite


text:
"""
Quick and dirty script to compute shape tensor of a FITS image.

:author: Sami-Matias Niemi
:contact: smn2@mssl.ucl.ac.uk

:version: 0.1
"""
import numpy as np
import pyfits as pf
import glob as g

def computeShapeTensor(data):
    """
    Computes a shape tensor from 2D imaging array.

    :Warning: This function has been adapted from Fortran
              and thus is very slow because of the nested
              loops. Need to be rewritten.

    :param data: imaging data as a numpy array
    :type data: ndarray

    :return: half of the size of the object in x and y direction
    :rtype: dict
    """
    data = data.transpose()
    xdim, ydim = data.shape

    Qxx = 0.
    Qxy = 0.
    Qyy = 0.
    for i in range(xdim):
        for j in range(ydim):
            Qxx += data[j, i] * (i - 0.5 * (xdim - 1)) * (i - 0.5 * (xdim - 1))
            Qxy += data[j, i] * (i - 0.5 * (xdim - 1)) * (j - 0.5 * (ydim - 1))
            Qyy += data[j, i] * (j - 0.5 * (ydim - 1)) * (j - 0.5 * (ydim - 1))

    shx = (Qxx + Qyy + np.sqrt((Qxx - Qyy)**2 + 4. * Qxy * Qxy)) / 2.
    shy = (Qxx + Qyy - np.sqrt((Qxx - Qyy)**2 + 4. * Qxy * Qxy)) / 2.

    shapex = np.sqrt(shx / np.sum(data))
    shapey = np.sqrt(shy / np.sum(data))

    return dict(shapex=shapex, shapey=shapey)

if __name__ == '__main__':
    files = g.glob('*.fits')

    for file in files:
        print ('File = %s' % file)
        data = pf.getdata(file)
        shape = computeShapeTensor(data)
        print (shape)

max_stars_repo_name: dcorre/ImSimpy

max_issues_repo_name: mtourneur/ImSimpy

max_forks_repo_name: mtourneur/ImSimpy


text:
import matplotlib.pyplot as plt
import numpy as np
import torch
import torchvision
import torchvision.transforms as transforms

from model import PREPROCESS as preprocess

MNIST_train = torchvision.datasets.MNIST(
    root = '/content/drive/MyDrive/Colab Notebooks/',
    download = True,
    train = True,
    transform = transforms.ToTensor()
)

# Reorganizing Data
train_dict = {i : MNIST_train.data[MNIST_train.targets == i].unsqueeze(1) / 255 for i in range(10)}
print(*[f'{imgs.shape[0]} images of Label {label} of shape {imgs.shape[1:]}' for label, imgs in train_dict.items()], sep='\n')
# for i in range(10):
#     train_dict[i] = MNIST_train.data[MNIST.targets == i]

# One Sample Image from each class
IMAGES = {0:train_dict[0][95], 1:train_dict[1][1], 2:train_dict[2][2], 3:train_dict[3][0]}
for i in [4,6,7,8]:
    IMAGES[i] = train_dict[i][1]
IMAGES[5] = train_dict[5][63]
IMAGES[9] = train_dict[9][5]

# Visualizing Sample images
GRID = torchvision.utils.make_grid(preprocess(torch.stack([*IMAGES.values()])), nrow=5)
plt.figure(figsize=(15,5))
plt.imshow(np.transpose(GRID, (1,2,0)))

# Training only on Images of '0','1','2'
zeros_images = train_dict[0][:100]
ones_images  = train_dict[1][:100]
twos_images  = train_dict[2][:100]

train_images = torch.stack([zeros_images, ones_images, twos_images], dim=1)
# train_images.shape  # train_images[:,0,...].shape = torch.Size([100, 1, 28, 28])

max_stars_repo_name: IvLabs/Face-Unlock