text:
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
import pmagpy.pmag as pmag
import pmagpy.ipmag as ipmag
import pmagpy.pmagplotlib as pmagplotlib
def make_plot(arch_df,edited_df,sect_depths,hole,\
              gad_inc,depth_min,depth_max,labels,spec_df=[], fignum=1,
              agemin=0,agemax=0):
    """
    Makes a downhole plot of IODP data.

    Parameters
    ___________
    arch_df : Pandas DataFrame
        dataframe of SRM archive measurements
    edited_df : Panas DataFrame
        dataframe of edited SRM archive measurements
    sect_depths : NumPy array
        array containing section depths of sections for plotting
    hole : str
        name of hole
    gad_inc : float
        inclination expected at the site from a GAD field
    depth_min : float
        minimum depth for plot
    depth_max : float
        maximum depth for plot
    labels :  Pandas Series
        series containing section names (e.g., 1H-2)
    spec_df : Pandas DataFrame
        dataframe of specimen data for plotting
    agemin : float
        if desired, minimum age for time scale plot
    agemax : float
        if non-zero an time scale plot will be generated

    """

    arch_df=arch_df[arch_df['core_depth']>depth_min]
    arch_df=arch_df[arch_df['core_depth']<=depth_max]
    edited_df=edited_df[edited_df['core_depth']>depth_min]
    edited_df=edited_df[edited_df['core_depth']<=depth_max]
    if len(spec_df)>0:
        spec_df=spec_df[spec_df['core_depth']>depth_min]
        spec_df=spec_df[spec_df['core_depth']<=depth_max]
        plot_spec=True
    else: plot_spec=False
    max_depth=arch_df.core_depth.max()
    min_depth=arch_df.core_depth.min()
    plot=1
    if agemax:
        col=5
        fig=plt.figure(fignum,(14,16))
    else:
        col=3
        fig=plt.figure(fignum,(8,20))
    ax=plt.subplot(1,col,plot)
    for d in sect_depths:
        if d<max_depth and d>min_depth:
            plt.axhline(d,color='black',linestyle='dashed',linewidth=.75)
    plot+=1
    plt.plot(np.log10(edited_df['magn_volume']*1e3),edited_df['core_depth'],\
            'co',markeredgecolor='grey')
    plt.plot(np.log10(arch_df['magn_volume']*1e3),arch_df['core_depth'],'k.',markersize=1)

    plt.ylabel('Depth (mbsf)')
    plt.xlabel('Log Intensity (mA/m)')
    plt.ylim(depth_max,depth_min)

    ax=plt.subplot(1,col,plot)
    for d in sect_depths:
        if d<max_depth and d>min_depth:
            plt.axhline(d,color='black',linestyle='dashed',linewidth=.75)
    plot+=1
    plt.plot(edited_df['dir_dec'],edited_df['core_depth'],'co',markeredgecolor='grey')
    plt.plot(arch_df['dir_dec'],arch_df['core_depth'],'k.',markersize=1)
    if plot_spec:
        plt.plot(spec_df['dir_dec'],spec_df['core_depth'],'r*',markersize=10)

    plt.axvline(180,color='red')

    plt.xlabel('Declination')
    plt.ylim(depth_max,depth_min)
    plt.ylim(depth_max,depth_min)

    plt.title(hole)
    ax=plt.subplot(1,col,plot)
    for d in sect_depths:
        if d<max_depth and d>min_depth:
            plt.axhline(d,color='black',linestyle='dashed',linewidth=.75)
    plot+=1
    plt.plot(edited_df['dir_inc'],edited_df['core_depth'],'co',markeredgecolor='grey')
    plt.plot(arch_df['dir_inc'],arch_df['core_depth'],'k.',markersize=1)
    if plot_spec:
        plt.plot(spec_df['dir_inc'],spec_df['core_depth'],'r*',markersize=10)

    plt.xlabel('Inclination')
    plt.axvline(gad_inc,color='blue',linestyle='dotted')
    plt.axvline(-gad_inc,color='blue',linestyle='dotted')
    plt.axvline(0,color='red')
    plt.xlim(-90,90)

    for k in range(len(labels.values)):
        if sect_depths[k]<max_depth and sect_depths[k]>=min_depth:
            plt.text(100,sect_depths[k],labels.values[k],verticalalignment='top')
    plt.ylim(depth_max,depth_min);
    if agemax:
        ax=plt.subplot(1,col,plot)
        ax.axis('off')
        ax=plt.subplot(1,col,plot+1)
        pmagplotlib.plot_ts(ax,agemin,agemax)
    plt.savefig('Figures/'+hole+'_'+str(fignum)+'.pdf')
    print ('Plot saved in', 'Figures/'+hole+'_'+str(fignum)+'.pdf')

def inc_hist(df,inc_key='dir_inc'):
    """
    Makes a histogram of inclination data from a data frame with 'dir_inc' as the inclination key

    Parameters
    __________
    df : Pandas DataFrame
        dataframe with  inclination data in "inc_key" column
    """
    plt.figure(figsize=(12,5))
    plt.subplot(121)
    plt.ylabel('Number of inclinations')
    sns.distplot(df[inc_key],kde=False,bins=24)
    plt.xlabel('Inclination')
    plt.xlim(-90,90)
    plt.subplot(122)
    plt.ylabel('Fraction of inclinations')
    sns.distplot(df[inc_key],bins=24)
    plt.xlabel('Inclination')
    plt.xlim(-90,90)

def demag_step(magic_dir,hole,demag_step,meas_file='srm_arch_measurements.txt',
               site_file='srm_arch_sites.txt',depth_key='core_depth',verbose=True):
    """
    Selects the desired demagnetization step, and puts the core/section/offset information
    into the returned data frame

    Parameters
    ___________
    magic_dir : str
        directory of the MagIC formatted files
    hole : str
        IODP Hole
    demag_step : float
        desired demagnetization step in tesla
    meas_file : str
        input measurement.txt format file for IODP  measurements
    site_file : str
        input sites.txt format file for  sites.
    verbose : boolean
        if True, announce return of dataframe
    Returns
    ___________
    DataFrame with selected step and additional metadata

    """
    arch_data=pd.read_csv(magic_dir+'/'+meas_file,sep='\t',header=1)
    depth_data=pd.read_csv(magic_dir+'/'+site_file,sep='\t',header=1)
    depth_data['specimen']=depth_data['site']
    depth_data=depth_data[['specimen',depth_key]]
    depth_data.sort_values(by='specimen')
    arch_data=pd.merge(arch_data,depth_data,on='specimen')
    arch_demag_step=arch_data[arch_data['treat_ac_field']==demag_step]
    pieces=arch_demag_step.specimen.str.split('-',expand=True)
    pieces.columns=['exp','hole','core','sect','A/W','offset']
    arch_demag_step['core_sects']=pieces['core'].astype('str')+'-'+pieces['sect'].astype('str')
    arch_demag_step['offset']=pieces['offset'].astype('float')
    arch_demag_step['core']=pieces['core']
    arch_demag_step['section']=pieces['sect']
    arch_demag_step['hole']=hole
    arch_demag_step.drop_duplicates(inplace=True)
    arch_demag_step.to_csv(hole+'/'+hole+'_arch_demag_step.csv',index=False)
    if verbose: print ("Here's your demagnetization step DataFrame")
    return arch_demag_step

def remove_ends(arch_demag_step,hole,core_top=80,section_ends=10):
    """
    takes an archive measurement DataFrame and peels off the section ends and core tops

    Parameters
    __________
    arch_demag_step : Pandas DataFrame
        data frame filtered by iodp_funcs.demag_step
    hole : str
        IODP hole
    core_top : float
        cm to remove from the core top
    section_ends : float
        cm to remove from the section ends

    Returns
    _______
    noends : DataFrame
        filtered data frame
    """

    noends=pd.DataFrame(columns=arch_demag_step.columns)
    core_sects=arch_demag_step.core_sects.unique()
    for core_sect in core_sects:
        cs_df=arch_demag_step[arch_demag_step['core_sects'].str.contains(core_sect)]
        if '-1' in core_sect:
            cs_df=cs_df[cs_df.offset>cs_df['offset'].min()+core_top] # top 80cm
        else:
            cs_df=cs_df[cs_df.offset>cs_df['offset'].min()+section_ends] # top 10 cm
        cs_df=cs_df[cs_df.offset<cs_df['offset'].max()-10]
        noends=pd.concat([noends,cs_df])
    noends.drop_duplicates(inplace=True)
    noends.fillna("",inplace=True)
    noends.to_csv(hole+'/'+hole+'_noends.csv',index=False)
    print ("Here's your no end DataFrame")
    return noends

def remove_disturbance(noends,hole):
    """
    takes an archive measurement DataFrame and removes disturbed intervals using DescLogic files

    Parameters
    __________
    noends : Pandas DataFrame
        data frame filtered by iodp_funcs.remove_ends
    hole : str
        IODP hole

    Returns
    _______
    nodist : DataFrame
        filtered data frame
    """
    disturbance_file=hole+'/'+hole+'_disturbances.xlsx'
    disturbance_df=pd.read_excel(disturbance_file)
    disturbance_df.dropna(subset=['Drilling disturbance intensity'],inplace=True)
    disturbance_df=disturbance_df[disturbance_df['Drilling disturbance intensity'].str.contains('high')]
    disturbance_df=disturbance_df[['Top Depth [m]','Bottom Depth [m]']]
    disturbance_df.reset_index(inplace=True)
    nodist=noends.copy(deep=True)
    for k in disturbance_df.index.tolist():
        top=disturbance_df.loc[k]['Top Depth [m]']
        bottom=disturbance_df.loc[k]['Bottom Depth [m]']
        nodist=nodist[(nodist['core_depth']<top) | (nodist['core_depth']>bottom)]
    nodist.sort_values(by='core_depth',inplace=True)
    nodist.drop_duplicates(inplace=True)
    nodist.fillna("",inplace=True)
# save for later
    nodist.to_csv(hole+'/'+hole+'_nodisturbance.csv',index=False)
    print ("Here's your no DescLogic disturbance DataFrame")
    return nodist

def no_xray_disturbance(nodist,hole):
    """
    takes an archive measurement DataFrame and removes disturbed intervals using XRAY disturbance files

    Parameters
    __________
    nodist : Pandas DataFrame
        data frame filtered by iodp_funcs.remove_disturbance
    hole : str
        IODP hole

    Returns
    _______
    no_xray_df : DataFrame
        filtered data frame
    """
    disturbance_file=hole+'/'+hole+'_disturbances.xlsx'
    disturbance_df=pd.read_excel(disturbance_file)
    disturbance_df.dropna(subset=['Drilling disturbance intensity'],inplace=True)
    disturbance_df=disturbance_df[disturbance_df['Drilling disturbance intensity'].str.contains('high')]
    xray_file=hole+'/'+hole+'_xray_disturbance.xlsx'
    xray_df=pd.read_excel(xray_file,header=2)
    no_xray_df=pd.DataFrame(columns=nodist.columns)
    xray_df=xray_df[['Core','Section','interval (offset cm)']]
    xray_df.dropna(inplace=True)
    if type(xray_df.Section)!='str':
        xray_df.Section=xray_df.Section.astype('int64')
    xray_df.reset_index(inplace=True)
    xray_df['core_sect']=xray_df['Core']+'-'+xray_df['Section'].astype('str')
    xr_core_sects=xray_df['core_sect'].tolist()
    nd_core_sects=nodist['core_sects'].tolist()
    used=[]
    # put in undisturbed cores
    for coresect in nd_core_sects:
        if coresect not in used and coresect not in xr_core_sects:
            core_df=nodist[nodist['core_sects'].str.match(coresect)]
            no_xray_df=pd.concat([no_xray_df,core_df])
            used.append(coresect)
# take out disturbed bits
    for coresect in xr_core_sects:
        core_df=nodist[(nodist['core_sects'].str.match(coresect))]
        x_core_df=xray_df[xray_df['core_sect']==coresect]
        core_sect_intervals=x_core_df['interval (offset cm)'].tolist()
        for core_sect_interval in core_sect_intervals:
            interval=core_sect_interval.split('-')
            top=int(interval[0])
            bottom=int(interval[1])
        # remove disturbed bit
            core_df=core_df[(core_df['offset']<top) | (core_df['offset']>bottom)]
        # add undisturbed bit to no_xray_df
        no_xray_df=pd.concat([no_xray_df,core_df])
    no_xray_df.sort_values(by='core_depth',inplace=True)
# save for later
    no_xray_df.drop_duplicates(inplace=True)
    no_xray_df.fillna("",inplace=True)
    no_xray_df.to_csv(hole+'/'+hole+'_noXraydisturbance.csv',index=False)
#meas_dicts = no_xray_df.to_dict('records')
#pmag.magic_write(magic_dir+'/no_xray_measurements.txt', meas_dicts, 'measurements')

    print ('Here is your DataFrame with no Xray identified disturbances')
    return no_xray_df

def adj_dec(df,hole):
    """
    takes an archive measurement DataFrame and adjusts the average declination by hole to 90 for "normal"

    Parameters
    __________
    df : Pandas DataFrame
        data frame of archive measurement data
    hole : str
        IODP hole

    Returns
    _______
    adj_dec_df : DataFrame
        adjusted declination  data frame
    core_dec_adj : dict
        dictionary of cores and the average declination
    """
    cores=df.core.unique()
    adj_dec_df=pd.DataFrame(columns=df.columns)
    core_dec_adj={}
    for core in cores:
        core_df=df[df['core']==core]
        di_block=core_df[['dir_dec','dir_inc']].values
        ppars=pmag.doprinc(di_block)
        if ppars['inc']>0: # take the antipode
            ppars['dec']=ppars['dec']-180
        core_dec_adj[core]=ppars['dec']
        core_df['adj_dec']=(core_df['dir_dec']-ppars['dec'])%360
        core_df['dir_dec']=(core_df['adj_dec']+90)%360 # set mean normal to 90 for plottingh
        adj_dec_df=pd.concat([adj_dec_df,core_df])
    adj_dec_df.fillna("",inplace=True)
    adj_dec_df.drop_duplicates(inplace=True)
    adj_dec_df.to_csv(hole+'/'+hole+'_dec_adjusted.csv',index=False)
    print ('Adjusted Declination DataFrame returned')
    return adj_dec_df,core_dec_adj

def plot_aniso(df,fignum=1,save_figs=False):
    """
    Makes a depth plot (christmas tree plot) of eigenvalues and eigenvectors for AMS data

    Parameters
    ___________
    df : Pandas DataFrame
        pandas data frame of a MagIC formatted specimens table with AMS data in it
    fignum : int
        matplotlib figure number
    save_figs : Boolean
       saves figures if True
    """
    v1_decs=df['v1_dec'].values
    v1_incs=df['v1_inc'].values
    v3_decs=df['v3_dec'].values
    v3_incs=df['v3_inc'].values
    ipmag.plot_net(fignum)
    ipmag.plot_di(dec=v1_decs,inc=v1_incs,marker='s',markersize=50,color='red')
    ipmag.plot_di(dec=v3_decs,inc=v3_incs,marker='o',markersize=50,color='black')
    plt.title('Core coordinates')
    if save_figs: plt.savefig('aniso_core.svg')
    if 'v1_dec_adj' in df.columns:
        fig2=ipmag.plot_net(fignum+1)
        v1_decs=df['v1_dec_adj'].values
        v3_decs=df['v3_dec_adj'].values

        ipmag.plot_di(dec=v1_decs,inc=v1_incs,marker='s',markersize=50,color='red')
        ipmag.plot_di(dec=v3_decs,inc=v3_incs,marker='o',markersize=50,color='black')
        plt.title('Declination Adjusted')
        if save_figs: plt.savefig('aniso_corr.svg')
    return

def convert_hole_depths(affine_file,hole_df,site,hole):
    """
    converts core depths to composite depths

    Parameters
    __________
    affine_file : str
        file name of affine downloaded from LIMS
    hole_df : Pandas DataFrame
        dataframe with data in MagIC measurement data after filtering with iodp_funcs.demag_step

    site : str
        IODP Site
    hole : str
        IODP Hole

    Returns
    _______
    hole_df : Pandas DataFrame
        hole_df with transformed composite_depth column after affine

    """

    affine=pd.read_csv(affine_file)
    affine['core']=affine['Core'].astype('str')+affine['Core type']
    affine['hole']=site+affine['Hole']
    hole_affine_df=affine[affine['hole'].str.match(hole)] # get the list for this hole
    hole_df['composite_depth']=hole_df['core_depth']
    cores=hole_affine_df.core.tolist()
    for core in cores:
        mbsf=hole_df[hole_df.core==core]['core_depth']
        if hole in affine['hole'].tolist() and  core in hole_affine_df.core.tolist():
            offset=hole_affine_df[hole_affine_df['core']==core]\
                    ['Cumulative offset (m)'].astype('float').values[0]
            hole_df.loc[hole_df.core==core,'composite_depth']=mbsf+offset
    hole_df['affine table']=affine_file
    return hole_df

def age_depth_plot(datums,paleo,size=100,depth_key='midpoint CSF-A (m)',title='Age_Model_',dmin=0,dmax=600,amin=0,amax=8,poly=3,figsize=(5,5)):
    """
    Makes an age depth plot from data frames with paleomagnetic and biostratigraphic tie points

    Parameters
    __________
    datums : Pandas DataFrame
        dataframe with paleomagnetic tie points
    paleo : Pandas DataFrame
        dataframe with biostratigraphic tie points
    size : integer
        size of symbols
    depth_key : str
        data frame column name for depths in datums dataframe
    title : str
        title for plot
    dmin :float
        minimum depth for plot
    dmax : float
        maximum depth for plot
    amin : float
        minimum age for plot
    amax : float
        maximum age for plot
    poly : int
        degree of polynomial for best fit curve
    figsize : tuple
        figure size

    Returns
    _________
    coeffs : array
        coefficients of best-fit polynomial

    """
    plt.figure(1,figsize)
# put on curve
    zero=pd.DataFrame(columns=datums.columns,index=[0])
    zero['Age (Ma)']=0
    zero[depth_key]=0
    datums=pd.concat([zero,datums])
    datums.dropna(subset=['Age (Ma)',depth_key],   inplace=True)
    datums.sort_values(by=['Age (Ma)'],inplace=True)
    coeffs=np.polyfit(datums['Age (Ma)'].values,datums[depth_key].values,poly)
    fit=np.polyval(coeffs,datums['Age (Ma)'].values)
    plt.plot(datums['Age (Ma)'].values,fit,'c-',lw=3,label='polynomial fit')

#plot the tie points
    datums.dropna(subset=['Hole'],inplace=True)
    datums['top']=datums[depth_key]-datums['range (+/-) (m)']
    datums['bot']=datums[depth_key]+datums['range (+/-) (m)']
    datums=datums[datums['bot']<dmax]
    holes=datums['Hole'].unique()
    colors=['r','b','k','g']
    pcolors=['c','m','y','r','k']
    age_key='Published Age\n(Ma)'
    mid_key='Mid depth \n(mbsf)'
    top_key='Top depth \n(mbsf)'
    bot_key='Bottom depth \n(mbsf)'
    for i in range(len(holes)):
# plot the paleo
        paleo_hole=paleo[paleo.Hole.str.match(holes[i])]
        diatoms=paleo_hole[paleo_hole.Type.str.contains('DIAT')]
        rads=paleo_hole[paleo_hole.Type.str.contains('RAD')]
        diatom_lo=diatoms[diatoms.Event.str.contains('LO')]
        diatom_lo.reset_index(inplace=True)
        diatom_fo=diatoms[diatoms.Event.str.contains('FO')]
        diatom_fo.reset_index(inplace=True)
        rad_lo=rads[rads.Event.str.contains('LO')]
        rad_lo.reset_index(inplace=True)
        rad_fo=rads[rads.Event.str.contains('FO')]
        rad_fo.reset_index(inplace=True)
        for k in rad_lo.index:
            plt.plot([rad_lo.iloc[k][age_key],rad_lo.iloc[k][age_key]],\
                 [rad_lo.iloc[k][top_key],rad_lo.iloc[k][bot_key]],'g-')
        plt.scatter(rad_lo[age_key].values,rad_lo[mid_key].values,\
                marker='v',color='w',edgecolor=pcolors[i],label=holes[i]+':Rad LO')
        for k in rad_fo.index:
            plt.plot([rad_fo.iloc[k][age_key],rad_fo.iloc[k][age_key]],\
                 [rad_fo.iloc[k][top_key],rad_fo.iloc[k][bot_key]],'g-')
        plt.scatter(rad_fo[age_key].values,rad_fo[mid_key].values,\
                marker='^',color='w',edgecolor=pcolors[i],label=holes[i]+':Rad FO')

        for k in diatom_fo.index:
            plt.plot([diatom_fo.iloc[k][age_key],diatom_fo.iloc[k][age_key]],\
                 [diatom_fo.iloc[k][top_key],diatom_fo.iloc[k][bot_key]],'b-')
        plt.scatter(diatom_fo[age_key].values,diatom_fo[mid_key].values,\
                marker='>',color='w',edgecolor=pcolors[i+1],label=holes[i]+':Diatom FO')

        for k in diatom_lo.index:
            plt.plot([diatom_lo.iloc[k][age_key],diatom_lo.iloc[k][age_key]],\
                 [diatom_lo.iloc[k][top_key],diatom_lo.iloc[k][bot_key]],'b-')
        plt.scatter(diatom_lo[age_key].values,diatom_lo[mid_key].values,\
                marker='<',color='w',edgecolor=pcolors[i+1],label=holes[i]+':Diatom LO')

# put on the pmag
        hole=datums[datums['Hole'].str.match(holes[i])]
        hole.reset_index(inplace=True)

        plt.scatter(hole['Age (Ma)'].values,hole[depth_key].values,\
            marker='*',s=size,color=colors[i],label=holes[i]+':PMAG')
        for k in hole.index:
             plt.plot([hole.iloc[k]['Age (Ma)'],hole.iloc[k]['Age (Ma)']],\
                 [hole.iloc[k]['top'],hole.iloc[k]['bot']],'k-')

        size=size-20
    plt.ylim(dmax,dmin)
    plt.xlim(amin,amax)
    plt.xlabel('Age (Ma), GTS12')
    plt.ylabel('Depth (mbsf)')
    plt.title(title)
    #plt.legend(bbox_to_anchor=(1.5,1))
    plt.legend()
    plt.savefig('Figures/'+title+'.pdf')
    return coeffs

def do_affine(affine_file,datums):
    """
    converts datum depths to affine depths

    Parameters
    __________
    affine_file : str
        file downloaded from LORE for affine transformation of depths from mbsf to CCSA
    datums : Pandas DataFrame
        datums data frame

    Returns
    _______
    datums_out : Pandas DataFrame
        datums dataframe with affined depths
    """
    cmb_top_key='Top Depth CCSF-A (m)'
    cmb_bot_key='Bottom Depth CCSF-A (m)'
    mbsf_top_key='Top Depth CSF-A (m)'
    mbsf_bot_key='Bottom Depth CSF-A (m)'
    datums.dropna(subset=['Core-Sect-Depth (top)'],inplace=True)
    holes=datums['Hole'].unique()
    affine=pd.read_csv(affine_file)
    datums['top_core_list']=datums['Core-Sect-Depth (top)'].str.split('-')
    datums['bot_core_list']=datums['Core-Sect-Depth (bottom)'].str.split('-')
    datums['top_core']=datums['top_core_list'].str.get(0)
    datums['bot_core']=datums['bot_core_list'].str.get(0)
    affine['core']=affine['Core'].astype('str')+affine['Core type']
    affine['hole']='Hole '+affine['Hole']
    aff_holes=affine['hole'].values.tolist()
    datums_out=pd.DataFrame(columns=datums.columns)
    for hole in holes:
        hole_df=datums[datums['Hole']==hole]
        hole_affine_df=affine[affine['hole']==hole]
        top_cores=hole_df['top_core'].values
        bot_cores=hole_df['bot_core'].values
        for core in top_cores:
            if hole in aff_holes and  core in hole_affine_df.core.tolist():
                offset=hole_affine_df[hole_affine_df['core']==core]\
                      ['Cumulative offset (m)'].astype('float').values[0]
                mbsf=hole_df[hole_df.top_core==core][mbsf_top_key]
                hole_df.loc[hole_df.top_core==core,cmb_top_key]=mbsf+offset
            else:
                mbsf=hole_df[hole_df.top_core==core][mbsf_top_key]
                hole_df.loc[hole_df.top_core==core,cmb_top_key]=mbsf # no offset possible

        for core in bot_cores:
            if hole in aff_holes and  core in hole_affine_df.core.tolist():
                offset=hole_affine_df[hole_affine_df['core']==core]\
                      ['Cumulative offset (m)'].astype('float').values[0]
                mbsf=hole_df[hole_df.bot_core==core][mbsf_bot_key]
                hole_df.loc[hole_df.bot_core==core,cmb_bot_key]=mbsf+offset
            else:
                mbsf=hole_df[hole_df.bot_core==core][mbsf_bot_key]
                hole_df.loc[hole_df.bot_core==core,cmb_bot_key]=mbsf # no offset possible

        datums_out=pd.concat([datums_out,hole_df])
    datums_out['affine table']=affine_file
    datums_out['midpoint CCSF-A (m)']=datums_out[cmb_bot_key]+.5*(datums_out[cmb_top_key]-datums_out[cmb_bot_key])
    return datums_out

def fix_aniso_data(aniso_df,core_dec_adj="",site_df=""):
    """
    picks apart the MagIC specimens table and
    creates columns with eigenvalues and eigenvectors for plotting

    Parameters
    __________
    aniso_df : Pandas DataFrame
        dataframe created by reading in a specimens.txt formatted file
    core_dec_adj : dict
        dictionary with {'core' : dec} for IODP cores and their average declinations
        returned by iodp_funcs.adj_dec
    site_df : Pandas dataframe
        dataframe with sites.txt information, including core_depth for enabling plotting versus depth

    Returns
    _______
    aniso_dec_adj_df : Pandas DataFrame
        dataframe with additional columns
    """

    aniso_df['aniso_v1_list']=aniso_df['aniso_v1'].str.split(":")
    aniso_df['tau1']=aniso_df['aniso_v1_list'].str.get(0).astype('float')
    aniso_df['v1_dec']=aniso_df['aniso_v1_list'].str.get(1).astype('float')
    aniso_df['v1_inc']=aniso_df['aniso_v1_list'].str.get(2).astype('float')
    aniso_df['aniso_v2_list']=aniso_df['aniso_v2'].str.split(":")
    aniso_df['tau2']=aniso_df['aniso_v2_list'].str.get(0).astype('float')
    aniso_df['v2_inc']=aniso_df['aniso_v2_list'].str.get(2).astype('float')
    aniso_df['v2_dec']=aniso_df['aniso_v2_list'].str.get(1).astype('float')
    aniso_df['aniso_v3_list']=aniso_df['aniso_v3'].str.split(":")
    aniso_df['tau3']=aniso_df['aniso_v3_list'].str.get(0).astype('float')
    aniso_df['v3_dec']=aniso_df['aniso_v3_list'].str.get(1).astype('float')
    aniso_df['v3_inc']=aniso_df['aniso_v3_list'].str.get(2).astype('float')

    pieces=aniso_df.specimen.str.split('-',expand=True)
    pieces.columns=['exp','hole','core','sect','A/W','offset']
    aniso_df['core']=pieces['core'].astype('str')
    cores=aniso_df.core.unique()
    aniso_dec_adj_df=pd.DataFrame(columns=aniso_df.columns)
    for core in cores:
        if core in core_dec_adj.keys():
            core_df=aniso_df[aniso_df.core.str.match(core)]
            core_df['v1_dec_adj']=(core_df['v1_dec']-core_dec_adj[core])%360
            core_df['v2_dec_adj']=(core_df['v2_dec']-core_dec_adj[core])%360
            core_df['v3_dec_adj']=(core_df['v3_dec']-core_dec_adj[core])%360
            aniso_dec_adj_df=pd.concat([aniso_dec_adj_df,core_df])

    site_df['specimen']=site_df['site']
    site_df=site_df[['specimen','core_depth']]
    aniso_dec_adj_df=aniso_dec_adj_df.merge(site_df,on='specimen')

    aniso_dec_adj_df.drop(columns=['aniso_v1_list','aniso_v2_list','aniso_v3_list'],inplace=True)
    return aniso_dec_adj_df

def kdes_for_incs(site_df,max_depth,interval=5,figsize=(2,12),depth_key='composite_depth',cmap='Oranges'):
    """
    creates "heat maps" of inclinations versus depth from archive half measurements

    Parameters
    __________
    site_df : Pandas DataFrame
        dataframe with archive half data from IODP cores.  Dataframe must have depth information
        either as composite_depth (returned by iodp_funcs.adj_depth) or core_depth (returned by
        iodp_funcs.demag_step or one of the editing functions (e.g., iodp_funcs.remove_ends).
    max_depth : float
        maximum depth for kde plot
    interval : float
        the interval over with the kdes will be created (e.g., 5 m)
    figsize : tuple
        size of the desired figure
    depth_key : str
        column name for the desired depth
    cmap : str
        matplotlib color map desired for the plot

    Returns
    _______
    fig,ax : matplotlib figure and axis  objects
    """

    start=0
    nplots=int(max_depth/interval)
    fig,ax=plt.subplots(nplots,1,figsize=figsize)
    fig.subplots_adjust(hspace=0,wspace=0)
    plt.xlabel('Inclination')
    ax[0].axhline(0,color='black')
    for i  in range(nplots):
        ax[i].xaxis.set_major_locator(plt.NullLocator())
        ax[i].yaxis.set_major_locator(plt.NullLocator())
        df=site_df[(site_df[depth_key]>=start)&(site_df[depth_key]<start+interval)]
        sns.kdeplot(df['dir_inc'],data2=df[depth_key],cmap=cmap,shade=True,ax=ax[i])
        ax[i].set_xlim(-90,90)
        ax[i].set_ylim(start+interval,start)
        ax[i].set_xlabel('')
        ax[i].axis('off')
        ax[i].axvline(-90,color='black')
        ax[i].axvline(90,color='black')
        start+=interval
    ax[i].axhline(interval,color='black')
    return fig,ax

def make_composite(affine_file,site,holes):
    """
    Uses an affine file downloaded from LORE to adjust depths for holes from a given IODP Site
    to adjust the depths to a composite depth.  Works on the output of iodp_func.dec_adj()

    Parameters
    __________
    affine_file : str
        affine file downloaded from LIMS for adjusting depths
    site : str
        name of site
    holes : list
        list of holes to adjust

    Returns
    _______
    site_df : Pandas DataFrame
        dataframe for site with combined data from all holes adjusted to composite depth scale
    """
    first=True
    for hole in holes: #holes:
        adj_dec_df=pd.read_csv(hole+'/'+hole+'_dec_adjusted.csv')
        adj_dec_df['hole']=hole
        if first:
            site_df=pd.DataFrame(columns=adj_dec_df.columns)
            first=False
        adj_depth=convert_hole_depths(affine_file,adj_dec_df,site,hole)
        adj_depth.to_csv(hole+'/'+hole+'_depth_adjusted.csv',index=False)
        site_df=pd.concat([site_df,adj_depth])
    return site_df

max_stars_repo_name: schwehr/PmagPy

max_issues_repo_name: schwehr/PmagPy

max_forks_repo_name: schwehr/PmagPy


text:
"""
Functions to run permutation methods
"""
import pyrft as pr
import numpy as np
from sklearn.utils import check_random_state
from scipy.stats import t

def boot_contrasts(lat_data, design, contrast_matrix, n_bootstraps = 1000, template = 'linear', replace = True, store_boots = 0, display_progress = 0):
    """ A function to compute the voxelwise t-statistics for a set of contrasts
      and their (two-sided) p-value by bootstrapping the residuals

    Parameters
    -----------------
    lat_data:  a numpy.ndarray of shape (dim, N) or an object of class field
        giving the data where dim is the spatial dimension and N is the number of subjects
        if a field then the fibersize must be 1 and the final dimension must be
        the number of subjects
    design: a numpy.ndarray of shape (N,p)
        giving the design matrix of covariates (p being the number of parameters)
    contrast_matrix: a numpy.ndarray of shape (L,p)
        corresponding to the contrast matrix, such that which each row is a
        contrast vector (where L is the number of constrasts)
    B: int,
        giving the number of bootstraps to do (default is 1000)
    t_inv: a python function
         that specifies the reference family the default is ss.t_inv_linear which
         corresponds to the linear reference family form the sansouci package
    replace:  Bool
        if True (default) then the residuals are sampled with replacement
        (i.e. a bootstrap), if False then they are sampled without replacement
        resulting in a permutation of the data
    store_boots: Bool,
          An optional input that allows the bootstrapped p-values to be stored
          if 1. Default is 0, i.e. no such storage.

    Returns
    -----------------
    minp_perm: a numpy.ndarray of shape (1, B),
          where the bth entry (for 1<=b<=B) is the minimum p-value (calculated
          over space) of the bth bootstrap. The 0th entry is the minimum p-value
          for the original data.
    orig_pvalues: an object of class field,
          of size (dim, n_constrasts) giving the p-value (calculated across subjects)
          of each of the voxels across the different contrasts and
          with the same mask as the original data
    pivotal_stats: a numpy.ndarray of shape (1,B)
          whose bth entry is the pivotal statistic of the bth bootstrap,
          i.e. min_{1 <= k <= m} t_k^-1(p_{(k:m)}(T_{n,b}^*)). These quantities
          are needed for joint error rate control. (At the moment it takes K = m.)

    Examples
    -----------------
    # 1D
    dim = 5; N = 30; categ = np.random.multinomial(2, [1/3,1/3,1/3], size = N)[:,1]
    X = pr.group_design(categ); C = np.array([[1,-1,0],[0,1,-1]]); lat_data = pr.wfield(dim,N)
    minP, orig_pvalues, pivotal_stats, _ = pr.boot_contrasts(lat_data, X, C)

    # 2D
    dim = (10,10); N = 30; categ = np.random.multinomial(2, [1/3,1/3,1/3], size = N)[:,1]
    X = pr.group_design(categ); C = np.array([[1,-1,0],[0,1,-1]]); lat_data = pr.wfield(dim,N)
    minP, orig_pvalues, pivotal_stats, _ = pr.boot_contrasts(lat_data, X, C)
    """
    #### Prep
    # Convert the data to be a field if it is not one already
    if isinstance(lat_data, np.ndarray):
        lat_data = pr.make_field(lat_data)

    # Ensure that the fibersize of the field is 1
    if isinstance(lat_data.fibersize, tuple):
        raise Exception("The fibersize of the field must be 1 dimensional")

    # Error check the inputs and obtain the size of X
    contrast_matrix, nsubj, n_params = pr.contrast_error_checking(lat_data,design,contrast_matrix)

    # Obtain the inverse template function (allowing for direct input as well!)
    if isinstance(template, str):
        _, t_inv = pr.t_ref(template)
    else:
        t_inv = template

    #### Main Function
    # Set random state
    rng = check_random_state(101)

    # Initialize the vector of pivotal statistics
    pivotal_stats = np.zeros(n_bootstraps)

    # Initialize a vector to store the minimum p-value for each permutation
    minp_perm = np.zeros(n_bootstraps)

    # Calculate the original statistic (used a the first permutation)
    orig_tstats, residuals = pr.constrast_tstats_noerrorchecking(lat_data, design, contrast_matrix)

    # Initialize the p-value field
    orig_pvalues = orig_tstats

    # Calculate the p-values
    # (using abs and multiplying by 2 to obtain the two-sided p-values)
    orig_pvalues.field = 2*(1 - t.cdf(abs(orig_tstats.field), nsubj-n_params))

    # Note need np.ravel as the size of orig_pvalues.field is (dim, L) i.e. it's not a vector!
    orig_pvalues_sorted = np.array([np.sort(np.ravel(orig_pvalues.field))])

    # Get the minimum p-value over voxels and contrasts (include the orignal in the permutation set)
    minp_perm[0] = orig_pvalues_sorted[0,0]

    # Obtain the pivotal statistic used for JER control
    pivotal_stats[0] = np.amin(t_inv(orig_pvalues_sorted))

    # Initialize the boostrap storage!
    bootstore = 0
    if store_boots:
        # Calculate the number of contrasts
        n_contrasts = contrast_matrix.shape[0]
        masksize_product = np.prod(lat_data.masksize)
        bootstore = np.zeros((n_contrasts*masksize_product, n_bootstraps))
        print(bootstore.shape)
        print(orig_pvalues_sorted.shape)
        bootstore[:,0] = orig_pvalues_sorted[0]

    # Calculate permuted stats
    # note uses the no error checking version so that the errors are not checked
    # for each bootstrap!
    lat_data_perm = lat_data
    for b in np.arange(n_bootstraps - 1):

        # Display progress
        if display_progress:
            pr.modul(b, 1)

        # Obtain a sample with replacement
        shuffle_idx = rng.choice(nsubj, nsubj, replace = replace)
        lat_data_perm.field = residuals[...,shuffle_idx]
        permuted_tstats, _ = pr.constrast_tstats_noerrorchecking(lat_data_perm, design, contrast_matrix)

        # Compute the permuted p-values
        # (using abs and multiplying by 2 to obtain the two-sided p-values)
        permuted_pvalues = 2*(1 - t.cdf(abs(permuted_tstats.field), nsubj-n_params))
        permuted_pvalues = np.array([np.sort(np.ravel(permuted_pvalues))])

        #Get the minimum p-value of the permuted data (over voxels and contrasts)
        minp_perm[b+1] = permuted_pvalues[0,0]

        #Obtain the pivotal statistic - of the permuted data - needed for JER control
        pivotal_stats[b + 1] = np.amin(t_inv(permuted_pvalues))
        # could be adjusted for K not m or in general some set A! (i.e. in the step down process)

        if store_boots:
            bootstore[:,b+1] = permuted_pvalues[0]

    return [minp_perm, orig_pvalues, pivotal_stats, bootstore]

def bootfpr(dim, nsubj, contrast_matrix, fwhm = 0, design = 0, n_bootstraps = 100, niters = 1000, pi0 = 1, alpha = 0.1, template = 'linear', replace = True, useboot = True):
    """ A function which calculates FWER and JER error rates using niters iterations

    Parameters
    -----------------
    dim: a tuple,
        giving the dimensions of the data to generate
    nsubj: int,
        giving the number of subjects to use
    C: a numpy.ndarray of shape (L,p)
          corresponding to the contrast matrix, such that which each row is a
          contrast vector (where L is the number of constrasts)
    design: a numpy.ndarray of size (N,p) or an int
        giving the covariates (p being the number of parameters), if set to be
        an integer then random category vectors are generated for each iteration
        and a corresponding design matrix selected
    fwhm: int,
        giving the fwhm with which to smooth the data (default is 0 i.e. generating
                                                white noise without smoothing)
    n_bootstraps: int,
        giving the number of bootstraps to do (default is 1000)

    niters: int,
      giving the number of iterations to use to estimate the FPR
    alpha: int,
       the alpha level at which to control (default is 0.1)
    t_inv: specifying the reference family (default is the linear reference family)
    replace:  Bool
        if True (default) then the residuals are sampled with replacement
        (i.e. a bootstrap), if False then they are sampled without replacement
        resulting in a permutation of the data
    useboot: Bool,
        determines whether to use bootstrapping to analyse the data or permutation,
        the default is True, i.e. to use bootstrapping

    Returns
    -----------------
    fpr_fwer: double,
        the false positive rate for FWER control
    fpr_jer: double,
        the false positive rate for JER control

    Examples
    -----------------
    # 1D no signal
    dim = 5; nsubj = 30; C = np.array([[1,-1,0],[0,1,-1]]);
    FWER_FPR, JER_FPR = pr.bootfpr(dim, nsubj, C)

    # 1D with signal

    # 2D
    dim = (10,10); nsubj = 30; C = np.array([[1,-1,0],[0,1,-1]]);
    FWER_FPR, JER_FPR = pr.bootfpr(dim, nsubj, C)

    # 2D with signal
    dim = (25,25); nsubj = 100; C = np.array([[1,-1,0],[0,1,-1]]);
    FWER_FPR, JER_FPR = pr.bootfpr(dim, nsubj, C, 8, 0, 100, 1000, 0.8)
    """
    # Initialize the FPR counter
    n_falsepositives_jer = 0 # jer stands for joint error rate here
    n_falsepositives_fwer = 0

    # Obtain ordered randomness
    rng = check_random_state(101)

    # If the design input is a matrix take this to be the design matrix
    # of the covariates (otherwise a random design is generated - see below)
    if not isinstance(design, int):
        design_2use = design

    # Obtain the inverse template function (allowing for direct input as well!)
    if isinstance(template, str):
        _, t_inv = pr.t_ref(template)
    else:
        # Allow the inverse function to be an input
        t_inv = template

    if len(contrast_matrix.shape) == 1:
        n_contrasts = 1
        n_groups = 1
    else:
        n_contrasts = contrast_matrix.shape[0]
        n_groups = contrast_matrix.shape[1]

    # Initialize the true signal vector
    nvox = np.prod(dim)
    m = nvox*n_contrasts
    ntrue = int(np.round(pi0 * m))
    signal_entries = np.zeros(m)
    signal_entries[ntrue:] = 1

    # Calculate the FPR
    for i in np.arange(niters):
        # Keep track of the progress.
        pr.modul(i,1)

        # Generate the data (i.e. generate stationary random fields)
        lat_data = pr.statnoise(dim,nsubj,fwhm)

        if isinstance(design, int):
            # Generate a random category vector with choices given by the design matrix
            categ = rng.choice(n_groups, nsubj, replace = True)

            # Ensure that all categories are present in the category vector
            while len(np.unique(categ)) < n_groups:
                print('had rep error')
                categ = rng.choice(n_groups, nsubj, replace = True)

            # Generate the corresponding design matrix
            design_2use = pr.group_design(categ)

        # Generate the signal by random shuffling the original signal
        # (if the proportion of signal is non-zero)
        if isinstance(dim, int):
            signal = pr.make_field(np.zeros((dim,n_contrasts)))
        else:
            signal = pr.make_field(np.zeros(dim + (n_contrasts,)))
        if pi0 < 1:
            rng = check_random_state(104)
            shuffle_idx = rng.choice(m, m, replace = False)
            shuffled_signal = signal_entries[shuffle_idx]
            spatial_signal2add = np.zeros(dim)
            for j in np.arange(n_contrasts):
                contrast_signal = shuffled_signal[j*nvox:(j+1)*nvox]
                signal.field[..., j] = contrast_signal.reshape(dim)
                spatial_signal2add += signal.field[..., j]
                subjects_with_this_contrast = np.where(categ==(j+1))[0]

                # Add the signal to the field
                for k in np.arange(len(subjects_with_this_contrast)):
                    lat_data.field[..., subjects_with_this_contrast[k]] += spatial_signal2add

        # Convert the signal to boolean
        signal.field = signal.field == 0
        if useboot:
            # Implement the bootstrap algorithm on the generated data
            minp_perm, orig_pvalues, pivotal_stats, _ = pr.boot_contrasts(lat_data, design_2use, contrast_matrix, n_bootstraps, t_inv, replace)
        else:
            perm_contrasts(lat_data, design_2use, contrast_matrix, n_bootstraps, t_inv)

        # Calculate the lambda alpha level quantile for JER control
        lambda_quant = np.quantile(pivotal_stats, alpha)

        # Calculate the null p-values
        null_pvalues = np.sort(orig_pvalues.field[signal.field])
        #null_pvalues = np.array([])
        #for j in np.arange(n_contrasts):
        #    null_pvalues = np.append(null_pvalues, orig_pvalues.field[signal.field[..., j], j])

        extended_null_pvalues = np.ones(m)
        extended_null_pvalues[0:len(null_pvalues)] = null_pvalues
        extended_null_pvalues_tinv = t_inv(np.asmatrix(extended_null_pvalues))
        null_pivotal_statistic = np.amin(extended_null_pvalues_tinv[0:len(null_pvalues)])

        # Check whether there is a JER false rejection or not
        # Use pivotal_stats[0] since it corresponds to the original pivotal statistic
        # (i.e. not bootstrapped)
        #if pivotal_stats[0] < lambda_quant:
        if null_pivotal_statistic < lambda_quant:
            n_falsepositives_jer = n_falsepositives_jer + 1

        # Calculate the alpha quantile of the permutation distribution of the minimum
        alpha_quantile = np.quantile(minp_perm, alpha)

        if np.amin(null_pvalues) < alpha_quantile:
            n_falsepositives_fwer = n_falsepositives_fwer + 1

    # Calculate the false positive rate over all iterations
    fpr_fwer = n_falsepositives_fwer/niters

    # Calculate the standard error
    std_error_fwer = 1.96*np.sqrt(fpr_fwer*(1-fpr_fwer)/niters)

    # Calculate the false positive rate over all iterations
    fpr_jer = n_falsepositives_jer/niters

    # Calculate the standard error
    std_error_jer = 1.96*np.sqrt(fpr_jer*(1-fpr_jer)/niters)

    # Print the results
    print('FWER: ', fpr_fwer, ' +/- ', round(std_error_fwer,4))
    print('JER: ', fpr_jer, ' +/- ', round(std_error_jer,4))

    return fpr_fwer, fpr_jer

def perm_contrasts(lat_data, design, contrast_vector, n_bootstraps = 100, template = 'linear'):
    """ A function to compute the voxelwise t-statistics for a set of contrasts
      and their p-value using Manly type permutation

  Parameters
  -----------------
  lat_data:  an object of class field consisting of data for N subjects
  X: a numpy.ndarray of shape (N,p)
        giving the covariates (p being the number of parameters)
  c: a numpy.ndarray of shape (1,p)
        corresponding to the contrast to use on the data
  B: int,
      giving the number of bootstraps to do (default is 1000)
  replace:  Bool
      if True (default) then the residuals are sampled with replacement
      (i.e. a bootstrap), if False then they are sampled without replacement
      resulting in a permutation of the data

  Returns
  -----------------
  tstat_field: an object of class field,
          which has spatial size the same as
                  input data and fibersize equal to the number of contrasts

  Examples
  -----------------
  dim = (10,10); N = 30; categ = np.random.multinomial(2, [1/3,1/3,1/3], size = N)[:,1]
  X = pr.group_design(categ); c = np.array([1,-1,0]); lat_data = pr.wfield(dim,N)
  minP, orig_pvalues, pivotal_stats = pr.perm_contrasts(lat_data, X, c)
    """
    # Convert the data to be a field if it is not one already
    if isinstance(lat_data, np.ndarray):
        lat_data = pr.make_field(lat_data)

    # Error check the inputs and obtain the size of X
    contrast_vector, nsubj, n_params = pr.contrast_error_checking(lat_data,design,contrast_vector)

    if contrast_vector.shape[0] > 1:
        raise Exception('c must be a row vector')

      # Obtain the inverse template function (allowing for direct input as well!)
    if isinstance(template, str):
        _, t_inv = pr.t_ref(template)
    else:
        # Allow the inverse function to be an input
        t_inv = template

    ### Main
    # Set random state
    rng = check_random_state(101)

    # Initialize a vector to store the minimum p-value for each permutation
    minp_perm = np.zeros(n_bootstraps)

    # Initialize a vector to store the pivotal statistics for each permutation
    pivotal_stats = np.zeros(n_bootstraps)

    # Calculate the original statistic (used a the first permutation)
    orig_tstats, _ = pr.constrast_tstats_noerrorchecking(lat_data, design, contrast_vector)
    orig_pvalues = orig_tstats
    orig_pvalues.field =  2*(1 - t.cdf(abs(orig_tstats.field), nsubj-n_params))

    # Note need np.ravel as the size of orig_pvalues.field is (dim, 1) i.e. it's not a vector!
    orig_pvalues_sorted = np.array([np.sort(np.ravel(orig_pvalues.field))])

    # Get the minimum p-value over voxels and contrasts (include the orignal in the permutation set)
    minp_perm[0] = orig_pvalues_sorted[0,0]

    # Obtain the pivotal statistics
    pivotal_stats[0] = np.amin(t_inv(orig_pvalues_sorted))

    # Calculate permuted stats
    # note use the no error checking version so that the errors are not checked
    # for each permutation!
    for b in np.arange(n_bootstraps - 1):
        print(b)
        shuffle_idx = rng.permutation(nsubj)
        permuted_tstats, _ = pr.constrast_tstats_noerrorchecking(lat_data, design[shuffle_idx, :], contrast_vector)
        permuted_pvalues = 2*(1 - t.cdf(abs(permuted_tstats.field), nsubj-n_params))
        permuted_pvalues = np.array([np.sort(np.ravel(permuted_pvalues))])

        #Get the minimum p-value of the permuted data (over voxels and contrasts)
        minp_perm[b+1] = permuted_pvalues[0,0]

        # Get the pivotal statistics needed for JER control
        pivotal_stats[b + 1] = np.amin(t_inv(permuted_pvalues))

    return minp_perm, orig_pvalues, pivotal_stats

ext: py

max_stars_repo_name: alexblnn/pyrft

max_issues_repo_name: alexblnn/pyrft

max_forks_repo_name: alexblnn/pyrft

reason: import numpy,from scipy


text:
import pandas as pd
import numpy as np
import re
import marcformat

class MarcExtractor(object):
    tag_marc_file = 'MARC_FILE'
    tag_filter_columns = 'FILTER_COLUMNS'
    tag_marc_output_file = 'MARC_OUTPUT_FILE'
    marcFile = ''
    marcOutFile = ''
    filteredColumns = []
    df = pd.DataFrame()
    df1 = pd.DataFrame()
    df2 = pd.DataFrame()
    chunkSize = 1000
    count = 0

    def __init__(self, config_file):
        self.__processConfigFile(config_file)
        pass

    def processDataSet(self):
        header = pd.DataFrame()
        for chunk in pd.read_csv(self.marcFile, chunksize=self.chunkSize, encoding='latin1'):
            if self.count == 0:
                header = [i.zfill(3) for i in chunk.columns]
            chunk.columns = header
            self.count += 1
            print(self.count)
            self.df2 = self.__filterColumns(chunk)
            self.__processColumns()
            self.df2 = self.__filterNone(self.df2)
            self.df2 = self.__nullToNone(self.df2)
            if self.df1.empty:
                self.df1 = self.df2
            else:
                self.df1 = pd.concat([self.df1, self.df2])

    def __processConfigFile(self, config_file):
        marc_file = re.compile('^' + self.tag_marc_file + '\s*\=\s*')
        filter_columns = re.compile('^' + self.tag_filter_columns + '\s*\=\s*')
        marc_out_file = re.compile('^' + self.tag_marc_output_file + '\s*\=\s*')
        with open(config_file, 'r') as configFile:
            for line in configFile.readlines():
                if re.search(marc_file, line):
                    self.marcFile = line.split('=')[-1].strip()
                elif re.search(filter_columns, line):
                    self.filteredColumns.extend(line.split('=')[-1].strip().replace(' ', '').split(','))
                    self.filteredColumns = [i.zfill(3) for i in self.filteredColumns]
                elif re.search(marc_out_file, line):
                    self.marcOutFile = line.split('=')[-1].strip()
        configFile.close()
        pass

    def __filterColumns(self, df):
        df = df[self.filteredColumns]
        df.columns = [marcformat.marcColumnDirectory[i] for i in self.filteredColumns]
        df = self.__filterNaN(df)
        return df

    def __processColumns(self):
        for column in self.df2:
            if column == "ISBN":
                self.df2[column] = self.df2[column].apply(lambda x: self.__clean_isbn(x))
            elif column == 'Author':
                self.df2[column] = self.df2[column].astype(str).apply(lambda x: self.__clean_author(x))
            elif column == 'Title':
                self.df2[column] = self.df2[column].apply(lambda x: self.__clean_title(x))
            elif column == 'Summary':
                self.df2[column] = self.df2[column].astype(str).apply(lambda x: self.__clean_summary(x))
            elif column == 'Topical_Term':
                self.df2['Topical_Main'] = self.df2[column].astype(str).apply(lambda x: self.__topical_main(x))
                self.df2['Topical_Geographic'] = self.df2[column].astype(str).apply(
                    lambda x: self.__topical_geographic(x))
                self.df2.drop('Topical_Term', axis=1, inplace=True)
            elif column == 'Genre':
                self.df2[column] = self.df2[column].astype(str).apply(lambda x: self.__clean_genre(x))
            elif column == 'BibID':
                self.df2[column] = self.df2[column].apply(str)
            else:
                raise Exception('ERROR: Unhandled MARC code!')
        pass

    def getResultDF(self):
        return self.df1

    def saveToOutFile(self):
        self.df1.to_csv(self.marcOutFile, index=False)

    @staticmethod
    def __lst_to_str(s):
        return "".join(str(i) for i in s)

    @staticmethod
    def __filterNaN(df):
        df.dropna(subset=['ISBN', 'Title'], inplace=True)
        return df

    @staticmethod
    def __nullToNone(df):
        df.where(pd.notnull(df), None)
        return df

    @staticmethod
    def __filterNone(df):
        df['ISBN'].replace(to_replace=[None], value=np.nan, inplace=True)
        df.dropna(subset=['ISBN', 'Title'], inplace=True)
        return df

    def __clean_isbn(self, x):
        regex = r"\ba\w+"
        if re.search(regex, x):
            match = re.search(regex, x)
            result = match.group()
            result = result.replace('a', '')
            result = str(result)
            if result:
                return result

    def __clean_author(self, x):

        regex1 = r"[$a\s]\w+"
        regex2 = r"[$]a\w+"
        if re.search(regex1, x):
            match = re.findall(regex1, x)
            match = self.__lst_to_str(match)
            match = match.replace(' ', '_')
            match = re.findall(regex2, match)
            match = self.__lst_to_str(match)
            match = match.replace('_', ' ').replace('$a', '')
            return match
        else:
            return x

    def __clean_title(self, x):
        regex1 = r"[$a+-]\w+"
        regex2 = r"[$]a\w+"
        # pattern = r"The maybe house"
        # if x.find(pattern) >= 0:
        #     print("here!")
        if re.search(regex1, x):
            x = x.replace(' ', '_').replace(',', '_').replace('-', '___').replace('\'', '____').replace('\"', '_____')
            match = re.findall(regex1, x)
            match = self.__lst_to_str(match)
            match = re.findall(regex2, match)
            match = self.__lst_to_str(match)
            match = match.replace('_____', '\"').replace('____', '\'').replace('___', '-').replace('__', ',_').replace(
                '_', ' ').replace('$a', '').strip()
            return match
        else:
            return x

    def __clean_summary(self, x):
        regex = r"\ba\w+"
        if re.search(regex, x):
            x = x.replace(' ', '_').replace(',', '_').replace('-', '___').replace('.', '____').replace('\'',
                                                                                '_____').replace(
                '\"', '______')
            match = re.search(regex, x)
            if match:
                result = match.group()
                result = result.replace('a', '', 1).replace('______', '\"').replace('_____', '\'').replace('____',
                                                                                '.').replace(
                    '___',
                    '-').replace(
                    '__', ',_').replace('_', ' ')
                return result
        else:
            return x

    def __clean_topical_term(self, x):
        x = x.replace(' ', '_')
        regex1 = r"[$]\w+"
        if re.search(regex1, x):
            match = re.findall(regex1, x)
            result = self.__lst_to_str(match).replace('_', ' ')
            return result
        else:
            return x

    def __topical_main(self, x):
        x = x.replace(' ', '_')
        regex1 = r"\b\$v\w+|\$a\w+"
        if re.search(regex1, x):
            match = re.findall(regex1, x)
            result = self.__lst_to_str(match).replace('_', ' ')
            if result.index('$a') == 0:
                result = result.replace('$a', '', 1).replace('$a', ',').replace('$v', ',')
            elif result.index('$v') == 0:
                result = result.replace('$v', '', 1).replace('$a', ',').replace('$a', ',')
            return result

    def __topical_geographic(self, x):
        x = x.replace(' ', '_')
        regex1 = r"[$]z\w+"
        if re.search(regex1, x):
            match = re.findall(regex1, x)
            result = self.__lst_to_str(match).replace('_', ' ')
            if result.index('$z') == 0:
                result = result.replace('$z', '', 1).replace('$z', ',')
            return result

    def __clean_genre(self, x):
        regex = r"\ba\w+"
        if re.search(regex, x):
            match = re.search(regex, x)
            result = match.group()
            result = result.replace('a', '', 1)
            return result
        else:
            return x

if __name__ == '__main__':
    marcExtractor = MarcExtractor('../config.ini')
    marcExtractor.processDataSet()
    marcExtractor.saveToOutFile()
    df = marcExtractor.getResultDF()
    print(df)

reason: import numpy


text:
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import numpy as np
from hypothesis import given
import hypothesis.strategies as st
from caffe2.python import core
import caffe2.python.hypothesis_test_util as hu

class TestBooleanMaskOp(hu.HypothesisTestCase):

    @given(x=hu.tensor(min_dim=1,
                       max_dim=5,
                       elements=st.floats(min_value=0.5, max_value=1.0)),
           **hu.gcs)
    def test_boolean_mask(self, x, gc, dc):
        op = core.CreateOperator("BooleanMask",
                                 ["data", "mask"],
                                 "masked_data")
        mask = np.random.choice(a=[True, False], size=x.shape[0])

        def ref(x, mask):
            return (x[mask],)

        self.assertReferenceChecks(gc, op, [x, mask], ref)
        self.assertDeviceChecks(dc, op, [x, mask], [0])

    @given(x=hu.tensor(min_dim=1,
                       max_dim=5,
                       elements=st.floats(min_value=0.5, max_value=1.0)),
           **hu.gcs)
    def test_boolean_mask_indices(self, x, gc, dc):
        op = core.CreateOperator("BooleanMask",
                                 ["data", "mask"],
                                 ["masked_data", "masked_indices"])
        mask = np.random.choice(a=[True, False], size=x.shape[0])

        def ref(x, mask):
            return (x[mask], np.where(mask)[0])

        self.assertReferenceChecks(gc, op, [x, mask], ref)
        self.assertDeviceChecks(dc, op, [x, mask], [0])


text:
import random

import numpy as np

from ...core.sampler.proposal import JumpProposal

class SkyLocationWanderJump(JumpProposal):
    """
    Jump proposal for wandering over the sky location. Does a Gaussian step in
    RA and DEC depending on the temperature.
    """

    def __call__(self, sample, **kwargs):
        temperature = 1 / kwargs.get('inverse_temperature', 1.0)
        sigma = np.sqrt(temperature) / 2 / np.pi
        sample['ra'] += random.gauss(0, sigma)
        sample['dec'] += random.gauss(0, sigma)
        return super(SkyLocationWanderJump, self).__call__(sample)

class CorrelatedPolarisationPhaseJump(JumpProposal):
    """
    Correlated polarisation/phase jump proposal. Jumps between degenerate phi/psi regions.
    """

    def __call__(self, sample, **kwargs):
        alpha = sample['psi'] + sample['phase']
        beta = sample['psi'] - sample['phase']

        draw = random.random()
        if draw < 0.5:
            alpha = 3.0 * np.pi * random.random()
        else:
            beta = 3.0 * np.pi * random.random() - 2 * np.pi
        sample['psi'] = (alpha + beta) * 0.5
        sample['phase'] = (alpha - beta) * 0.5
        return super(CorrelatedPolarisationPhaseJump, self).__call__(sample)

class PolarisationPhaseJump(JumpProposal):
    """
    Correlated polarisation/phase jump proposal. Jumps between degenerate phi/psi regions.
    """

    def __call__(self, sample, **kwargs):
        sample['phase'] += np.pi
        sample['psi'] += np.pi / 2
        return super(PolarisationPhaseJump, self).__call__(sample)


text:
import numpy as np
import torch
from rlkit.util.utils import int_to_onehot
import rlkit.torch.pytorch_util as ptu

class Preprocessor(object):
    def __init__(self, input_type, action_dim, goal_dim, obs_dim, goal_conditioned=False, network=None):
        self.input_type = input_type
        self.action_dim = action_dim
        self.goal_dim = goal_dim
        self.obs_dim = obs_dim
        self.goal_conditioned = goal_conditioned
        self.network = network

    def preprocess(self, obs):
        # if self.expert:
        #     actions = torch.tensor(int_to_onehot(self.expert.get_actions(obs), action_dim), dtype=torch.float32)
        #     return torch.cat((actions, obs[:, obs_dim:]), dim=1)
        # return x

        if self.input_type == 'action':
            network_output = self.network.get_actions(obs)
            actions = torch.tensor(int_to_onehot(network_output, self.action_dim), dtype=torch.float32)
            return torch.cat((actions, obs[:, self.obs_dim:]), dim=1)
        elif self.input_type == 'obs':
            return obs
        elif self.input_type == 'obs-action':
            # action, obs, goal
            network_output = self.network.get_actions(obs)
            actions = torch.tensor(int_to_onehot(network_output, self.action_dim), dtype=torch.float32).to(ptu.device)
            return torch.cat((actions, obs), dim=1)  # not sure if this is right
        elif self.input_type == 'latent':
            # g_enc(obs), goal
            # import ipdb; ipdb.set_trace()
            if self.goal_conditioned:
                latent, _ = self.network.encode(obs)
            else:
                latent, _ = self.network.encode(obs[:, :self.obs_dim])
            return torch.cat((latent, obs[:, self.obs_dim:]), dim=1)
        elif self.input_type == 'obs-latent':
            # g_enc(obs), obs, goal
            if self.goal_conditioned:
                latent, _ = self.network.encode(obs)
            else:
                latent, _ = self.network.encode(obs[:, :self.obs_dim])
            return torch.cat((latent, obs), dim=1)
        elif self.input_type == 'latentencgoal':
            # obs, g_enc(obs), g_enc(goal)
            if self.goal_conditioned:
                latent, _ = self.network.encode(obs)
                goal_latent, _ = self.network.encode(torch.cat((obs[:, self.obs_dim], obs[:, self.obs_dim:]), dim=1))
            else:
                latent, _ = self.network.encode(obs[:, :self.obs_dim])
                goal_latent, _ = self.network.encode(obs[:, self.obs_dim:])

            return torch.cat((latent, goal_latent, obs[:, :self.obs_dim]), dim=1)
        else:
            raise NotImplementedError

max_stars_repo_name: YashBit/accelerate-skillDiscovery

max_issues_repo_name: YashBit/accelerate-skillDiscovery

max_forks_repo_name: YashBit/accelerate-skillDiscovery


text:
import face_recognition
import cv2
import numpy as np
import speech_recognition as sr
import os

'''
# Load a sample picture and learn how to recognize it.
obama_image = face_recognition.load_image_file("obama.jpg")
obama_face_encoding = face_recognition.face_encodings(obama_image)[0]

# Load a second sample picture and learn how to recognize it.
biden_image = face_recognition.load_image_file("biden.jpg")
biden_face_encoding = face_recognition.face_encodings(biden_image)[0]

vaibhav_image = face_recognition.load_image_file("new_vaibhav.jpeg")
vaibhav_face_encoding = face_recognition.face_encodings(vaibhav_image)[0]

rugved_image = face_recognition.load_image_file("rugved.jpeg")
rugved_face_encoding = face_recognition.face_encodings(rugved_image)[0]

# Create arrays of known face encodings and their names
known_face_encodings = [
    obama_face_encoding,
    biden_face_encoding,
    vaibhav_face_encoding,
    rugved_face_encoding
]
known_face_names = [
    "Barack Obama",
    "Joe Biden",
    "vaibhav",
    "Zural"
]
'''

known_face_encodings = []
known_face_names = []
for enc in os.listdir('encoded_image'):
    known_face_names.append(enc.split('.')[0])
    known_face_encodings.append(np.fromfile('encoded_image/'+enc))

#print('Learned encoding for', len(known_face_encodings), 'images.')

# Initialize some variables
face_locations = []
face_encodings = []
face_names = []
process_this_frame = True


text:
import numpy as np
from numpy.linalg import inv
from numpy.linalg import svd
from numpy.linalg import eig
from numpy.linalg import det
from scipy.optimize import leastsq, least_squares, fmin
import pandas as pd
import numpy as np
import time
import random
import os.path
import argparse

from epipolar_geometry import EpipolarGeometry
from model import RecModel
from image_sequence import ImageSequence

class UncalibratedReconstruction:
    '''
    Class that contains high level methods to perform 3D reconstruction from a sequence of uncalibrated images.
    '''

    def __init__(self, sequence_length, width, height, triang_method=0, opt_triang=0, opt_f=1, self_foc=0):
        '''
        Constructor
        Args:
            sequence_length: number of images (views)
            width: width of the images
            height: height of the images
            triang_method: triangulation method (0: standard, 1: polynomial)
            opt_triang: optimize initial 3D point estimate
            opt_f: optimize fundamental matrix estimation
            self_foc: for self-calibration, type of focal length expected across views (0: fixed, 1: varying )

        '''
        self._eg_utils = EpipolarGeometry()

        # things that are needed throught the class
        self._sequence_length = sequence_length
        self._width = width
        self._height = height
        self._mm = (width + height) / 2

        self._triangulation_method = triang_method
        self._optimize_triangulation = opt_triang
        self._optimize_f = opt_f
        self._self_foc = self_foc

    def two_view_geometry_computation(self, view1_feat2D, view2_feat2D):
        '''
        Method to compute the fundamental matrix and epipoles for two views
        Args:
            view1_feat2D: 2D feature coordinates in view 1
            view2_feat2D: 2D feature coordinates in view 2
        Returns:
            F: the fundamental matrix
            epipole_1: view1 epipole
            epipole_2: view2 epipoles

        '''
        F = self._eg_utils.fundamental_matrix(view1_feat2D, view2_feat2D, self._optimize_f)
        epipole_1 = self._eg_utils.get_epipole(F)
        epipole_2 = self._eg_utils.get_epipole(np.transpose(F))
        return F, epipole_1, epipole_2

    def compute_reference_frame(self, epipole, F):
        '''
        Method to compute the reference frame of the reconstruction (i.e. plane at infinity in an affine or metric space).
        Args:
            epipole: the epipole
            F: the fundamental matrix
        Returns:
            p: the reference plane
            h: the homography [e]xF

        '''
        H = self._eg_utils.compute_homography(epipole, F)  # compute the homography [e]xF
        # get the reference plane
        p = np.sum(np.divide(np.eye(3) - H, np.transpose(np.asarray([epipole, epipole, epipole]))), axis=0) / 3
        # adjust reference plane to make the first two projection matrices as equal as possible
        p = fmin(self.init_plane, np.append(p, 1), xtol=1e-25, ftol=1e-25, args=(H.real, epipole.real))
        p = p[0:3]
        return p, H

    def init_plane(self, p, H, epi):
        '''
        Error function to make the difference between the first two projection matrices as small as possible
        Note: assuming that the two views used for the initial reconstruction are not too far apart (this their projection matrices are almost equal), has proven to give good results
        Args:
            p: the reference plane (i.e. plane at infinity)
            H: homography [e]x[F]
            epi: epipola
        Returns:
            error: difference between two projection matrices
        '''
        epi = np.reshape(epi, (3, 1))
        p = np.reshape(p, (1, 4))

        t = p[0, 0:3]
        t = np.reshape(t, (1, 3))

        error = sum(sum(abs(H + epi.dot(t) - p[0, 3] * np.eye(3))))

        return error

    def estimate_initial_projection_matrices(self, H, epipole_2, p):
        '''
        Method to estimate the projection matrices for the two views (i.e. P1=[I | 0], P2=[H+epi1|e])
        Args:
            H: homography [e]x[F]
            epipole_2: epipole in the 2nd view
            p: the reference plane of the reconstruction (i.e. plane at infinity)
        Returns:
            P: projection matrices for these two views

        '''
        P = np.zeros((3, 4, self._sequence_length))
        P[:, :, 0] = [[1, 0, 0, 0], [0, 1, 0, 0],
                      [0, 0, 1, 0]]  # P1=[I | 0], i.e. frist frame aligned with world frame
        epi_tmp = np.reshape(epipole_2, (3, 1))  # P2=[H+epi1|e]
        P[:, :3, 1] = H + epi_tmp.dot(np.reshape(p, (1, 3)))
        P[:, 3, 1] = epipole_2
        P[:, :, 1] = P[:, :, 1] / P[2, 2, 1]
        return P

    def get_initial_structure(self, feat_2D, P, epipole_1, epipole_2, F):
        '''
        Method to get an initial 3D structure (i.e. 3D point cloud), from the first two projection matrices through triangulation.
        Args:
            feat_2D: 2D feature coordinates for all images
            P: projection matrices for all views (only the first two views are used)
            epipole_1: view 1 epipole
            epipole_2: view 2 epipole
            F: fundamental matrix
        Returns:
            points3D: 3D point cloud
        '''
        number_of_features = feat_2D.shape[2]
        points3D = np.zeros(shape=[number_of_features, 4])
        for i in range(0, number_of_features):

            if (self._triangulation_method == 0):
                x = self._eg_utils.triangulate_points(feat_2D[0, :, i], feat_2D[1, :, i], P[:, :, 0], P[:, :, 1])
                x = x[0:3]
            elif (self._triangulation_method == 1):
                x = self._eg_utils.polynomial_triangulation(feat_2D[1, :, i], feat_2D[1, :, i], epipole_1, epipole_2, F,
                                                            P[:, :, 0], P[:, :, 1])
                x = x[0:3] / x[3]  # normalize

            if (self._optimize_triangulation == 1):
                # refine 3D point estimation (due to noise, lines of sight may not intersect perfectly). Minimizations should be carried out in the images
                # and not in the projective 3D space, thus the reprojection error is used.
                x = fmin(self.refine_3d_point, x, xtol=1e-25, ftol=1e-25, full_output=0,
                         args=(P[:, :, 0], P[:, :, 1], feat_2D[0, :, i], feat_2D[1, :, i]))

            points3D[i, :] = np.append(x, 1)

        return points3D

    def refine_3d_point(self, point3D, P1, P2, view1_feat2D, view2_feat2D):
        '''
        Method to compute the reprojection error of a 3D point in two views
        Args:
            point3D: 3D point cloud
            P1: projection matrix of view 1
            P2: projection matrix of view 2
            view1_feat2D: 2D feature coordinates in view 1
            view2_feat2D: 2D feature coordinates in view 1
        Returns:
            error: the reprojection error

        '''
        point3D = np.append(point3D, 1)
        error = sum(
            self.compute_reprojection_error_point(P1, point3D, view1_feat2D) + self.compute_reprojection_error_point(P2,
                                                                                point3D,
                                                                                view2_feat2D))
        # sdfds
        return error

    def projective_pose_estimation(self, feat_2D, P, points3D):
        '''
        Method to add views using an initial 3D structure, i.e. compute the projection matrices for all the additional views (the first two are already
        estimated in previous steps)
        Args:
            feat_2D: 2D feature coordinates for all images
            P: projection matrices
            points3d: 3D point cloud
        Returns:
            P: projection matrices for all views
        '''
        number_of_features = feat_2D.shape[2]

        AA = np.zeros(shape=[2 * number_of_features, 12])

        for i in range(2, self._sequence_length):
            for j in range(0, number_of_features):
                AA[2 * j, 0:4] = points3D[j]
                AA[2 * j, 8:12] = -feat_2D[i, 0, j] * points3D[j]
                AA[2 * j + 1, 4:8] = points3D[j]
                AA[2 * j + 1, 8:12] = -feat_2D[i, 1, j] * points3D[j]

            U, s, Vh = svd(AA)
            V = np.transpose(Vh)

            VV = V[0:12, 11]
            VV = VV / VV[10]
            VV = np.delete(VV, 10)

            # refine the estimate for the i-th projection matrix
            result = least_squares(self._eg_utils.refine_projection_matrix, VV, args=(points3D, feat_2D[i, :, :]))
            VV = result.x

            Pr = np.zeros(shape=[3, 4])
            Pr[0, :] = VV[0:4]
            Pr[1, :] = VV[4:8]
            Pr[2, :] = np.append(np.append(VV[8:10], 1), VV[10])
            P[:, :, i] = Pr

        return P

    def bundle_adjustment(self, feat_2D, P, feat3D):
        '''
        Method to refine structure and motion, i.e. refine the projection matrices and 3D points using the reprojection error
        Args:
            feat_2D: 2D feature coordinates for all images
            P: projection matrices
            points3d: 3D point cloud
        Returns:
            P: the refined projection matrices
            feat3D: the refined 3D point cloud
            error: the reprojection error
        '''
        number_of_features = feat_2D.shape[2]

        # The vector to be optimized
        X = np.reshape(P[:, :, 0], (1, 12))

        # Append the projection matrices...
        for i in range(1, self._sequence_length):
            X = np.append(X, np.reshape(P[:, :, i], (1, 12)))
        X = np.delete(X, [10, 22, (self._sequence_length - 1) * 12 + 10])

        # ...and then append the 3D points
        X = np.append(X, np.reshape(feat3D[:, 0:3], number_of_features * self._sequence_length))

        # Optimize using Levenberg-Marquardt
        result = least_squares(self._eg_utils.overall_reprojection_error, X, max_nfev=1000, method='lm',
                               args=([feat_2D]))
        X = result.x
        error = np.power(sum(self._eg_utils.overall_reprojection_error(X, feat_2D)), 2)

        # get the refined projection matrices from the optimal vector
        for i in range(0, self._sequence_length):
            P[:, :, i] = np.reshape(X[0 + i * 11:12 + i * 11], (3, 4))
            P[2, 3, i] = P[2, 2, i]
            P[2, 2, i] = 1

        # get the refined 3D coordinates from the optimal vector
        feat3D[:, 0:3] = np.reshape(
            X[self._sequence_length * 11:self._sequence_length * 11 + self._sequence_length * number_of_features * 3],
            (number_of_features, 3))

        Tp1 = np.vstack([P[:, :, 0], [0, 0, 0, 1]])
        for i in range(0, self._sequence_length):
            P[:, :, i] = P[:, :, i].dot(inv(Tp1))

        feat3D = Tp1.dot(np.transpose(feat3D))
        feat3D = np.transpose(feat3D / feat3D[3, :])

        return P, feat3D, error

    def self_calibration(self, P):
        '''
        Self calibration using the procedure described in
        M. Pollefeys, R. Koch and L. Van Gool, "Self-Calibration and Metric Reconstruction in spite of Varying and Unknown Internal Camera Parameters", Proc. International Conference on Computer Vision, Narosa Publishing House, pp.90-95, 1998.
        Args:
            P: projection matrices
        Returns:
            Tm: transformation matrix that will transform from the projective space to metric space
            K: camera intrisic parameters for each view
            error: the reprojection error

        '''
        # setup the system of equations
        AAA = np.zeros(shape=[4 * self._sequence_length - 4, 6])
        for i in range(0, self._sequence_length - 1):
            P_tmp = P[:, :, i + 1]
            AAA[0 + 4 * i, :] = [(-np.power(P_tmp[1, 1], 2) + np.power(P_tmp[0, 1], 2) - np.power(P_tmp[1, 0],
                                                                                2) + np.power(
                P_tmp[0, 0], 2)), (-2 * P_tmp[1, 0] * P_tmp[1, 3] + 2 * P_tmp[0, 0] * P_tmp[0, 3]),
                                 (-2 * P_tmp[1, 1] * P_tmp[1, 3] + 2 * P_tmp[0, 1] * P_tmp[0, 3]),
                                 (2 * P_tmp[0, 2] * P_tmp[0, 3] - 2 * P_tmp[1, 2] * P_tmp[1, 3]),
                                 (-np.power(P_tmp[1, 3], 2) + np.power(P_tmp[0, 3], 2)),
                                 (-np.power(P_tmp[1, 2], 2) + np.power(P_tmp[0, 2], 2))]
            AAA[1 + 4 * i, :] = [(P_tmp[1, 0] * P_tmp[0, 0] + P_tmp[1, 1] * P_tmp[0, 1]),
                                 (P_tmp[1, 0] * P_tmp[0, 3] + P_tmp[1, 3] * P_tmp[0, 0]),
                                 (P_tmp[1, 1] * P_tmp[0, 3] + P_tmp[1, 3] * P_tmp[0, 1]),
                                 (P_tmp[1, 2] * P_tmp[0, 3] + P_tmp[1, 3] * P_tmp[0, 2]), P_tmp[1, 3] * P_tmp[0, 3],
                                 P_tmp[1, 2] * P_tmp[0, 2]]
            AAA[2 + 4 * i, :] = [(P_tmp[2, 0] * P_tmp[0, 0] + P_tmp[2, 1] * P_tmp[0, 1]),
                                 (P_tmp[2, 0] * P_tmp[0, 3] + P_tmp[2, 3] * P_tmp[0, 0]),
                                 (P_tmp[2, 1] * P_tmp[0, 3] + P_tmp[2, 3] * P_tmp[0, 1]),
                                 (P_tmp[2, 2] * P_tmp[0, 3] + P_tmp[2, 3] * P_tmp[0, 2]), P_tmp[2, 3] * P_tmp[0, 3],
                                 P_tmp[2, 2] * P_tmp[0, 2]]
            AAA[3 + 4 * i, :] = [(P_tmp[2, 0] * P_tmp[1, 0] + P_tmp[2, 1] * P_tmp[1, 1]),
                                 (P_tmp[2, 0] * P_tmp[1, 3] + P_tmp[2, 3] * P_tmp[1, 0]),
                                 (P_tmp[2, 1] * P_tmp[1, 3] + P_tmp[2, 3] * P_tmp[1, 1]),
                                 (P_tmp[2, 2] * P_tmp[1, 3] + P_tmp[2, 3] * P_tmp[1, 2]), P_tmp[2, 3] * P_tmp[1, 3],
                                 P_tmp[2, 2] * P_tmp[1, 2]]

        U, s, Vh = svd(AAA)
        V = np.transpose(Vh)
        x = V[0:5, 5] / V[5, 5]
        jj = np.sign(x[0])
        b = x * np.sign(x[0])

        # initial estimate of the absolute conic
        W = np.asarray([[b[0], 0, 0, b[1]], [0, b[0], 0, b[2]], [0, 0, 1, b[3]], [b[1], b[2], b[3], b[4]]])

        # initial estimate of the focal lengths
        y = np.ones(shape=[self._sequence_length, 1])
        for i in range(0, self._sequence_length):
            y[i] = np.sqrt(np.abs(P[0, :, i].dot(W).dot(np.transpose(P[0, :, i]))) / (
                P[2, :, i].dot(W).dot(np.transpose(P[2, :, i]))))

        if (self._self_foc == 0):
            # optimize for fixed focal lengths
            pp2 = np.asarray([-b[1] / b[0], -b[2] / b[0], -b[3]])
            x = np.hstack((sum(y) / y.shape[0], pp2))
            x = fmin(self.fixed_f_error, x, args=(P, self._sequence_length))
            error = self.fixed_f_error(x, P, self._sequence_length)

            # fill out the camera instrisic parameters.
            K = np.zeros((3, 3, self._sequence_length))
            for i in range(0, self._sequence_length):
                K[:, :, i] = np.eye(3)
                K[0, 0, i] = x[0] * self._mm
                K[1, 1, i] = x[0] * self._mm
                K[0, 2, i] = self._width
                K[1, 2, i] = self._height

            inf_plane = x[1:4]

            # construct the transformation matrix that will take us from the projective space to to metric
            a = inv([[x[0], 0, 0], [0, x[0], 0], [0, 0, 1]])
            a = np.asarray(a) * jj
            tmp = np.asarray([0, 0, 0])
            tmp = np.reshape(tmp, (3, 1))
            Tm = np.append(a, tmp, 1)
            Tm = np.vstack((Tm, np.append(inf_plane, 1)))

        else:
            # optimize for varying focal lenghts (to do)
            print("not yet supported")

        return Tm, K, error

    def fixed_f_error(self, x, P, n):
        '''
        Error function for the self-calibration error when the focal lengths are fixed (i.e. we
        assume the same focal lengths accross the image sequence)
        '''
        K1 = np.eye(3)
        K1[0, 0] = x[0]
        K1[1, 1] = x[0]
        pp = np.asarray(x[1:4])
        pp = np.reshape(pp, (3, 1))

        W = np.asarray(np.append(K1.dot(np.transpose(K1)), -K1.dot(np.transpose(K1)).dot(pp), 1))
        tmp = np.append((np.transpose(-pp).dot(K1).dot(np.transpose(K1))),
                        np.transpose(pp).dot(K1).dot(np.transpose(K1)).dot(pp))
        W = np.vstack((W, tmp))

        a = K1.dot(np.transpose(K1))
        b = a / np.sqrt(np.trace(np.transpose(a).dot(a)))
        error = 0
        for i in range(0, n - 1):
            P_tmp = P[:, :, i + 1]
            c = P_tmp.dot(W).dot(np.transpose(P_tmp))
            d = c / np.sqrt(np.trace(np.transpose(c) * c))
            error = error + np.trace(np.transpose((b - d)).dot(b - d))

        return error

    def convert_to_metric_space(self, Tm, feat3D, P, K):
        '''
        Transform the 3D points and projective matrices to the metric space
        Args:
            Tm: transformation matrix for transformic from projective to metric space
            feat3D: 3D point cloud
            P: projection matrices (for all views)
            K: camera intrisic parameters (for all views)
        Returns:
            feat3D: 3D point cloud in metric space
            P: projectio matrices in metric space
        '''
        # transform the projective 3d coordinates to metric
        InvT = Tm.dot(np.eye(4))
        a = inv(InvT)
        InvT = InvT * a[3, 3]

        feat3D = InvT.dot(np.transpose(feat3D))
        feat3D = np.transpose(feat3D / feat3D[3, :])

        # Rescale the projection matrices to width,height (no -1, 1)
        tmp = np.eye(3)
        tmp[0, 0] = self._mm
        tmp[1, 1] = self._mm
        tmp[0, 2] = self._width
        tmp[1, 2] = self._height
        for i in range(0, self._sequence_length):
            P[:, :, i] = P[:, :, i].dot(inv(InvT))
            P[:, :, i] = tmp.dot(P[:, :, i])
            a = det(inv(K[:, :, i]).dot(P[0:3, 0:3, i]))
            P[:, :, i] = P[:, :, i] * np.sign(a) / np.power(abs(a), 0.333)

        return feat3D, P

def main(input_file, show):
    print("--------------------------------")
    print(" Uncalibrated 3D Reconstruction ")
    print("")
    print("--------------------------------")

    sequence = ImageSequence(input_file)

    start_time = time.time()
    rec_engine = UncalibratedReconstruction(sequence.length, sequence.width, sequence.height)

    # normalize coordinates
    norm_feat_2d = sequence.get_normalized_coordinates()

    print("> Estimating fundamental matrix....")
    F, epipole_1, epipole_2 = rec_engine.two_view_geometry_computation(norm_feat_2d[0], norm_feat_2d[1])

    print("> Computing reference plane....")
    # Step 2: compute the reconstruction reference plane using the epipole in the second image
    p, H = rec_engine.compute_reference_frame(epipole_2, F)

    print("> Estimating projection matrices for first two views....")
    # Step 3: Estimate projection matrices for the first two views
    P = rec_engine.estimate_initial_projection_matrices(H, epipole_2, p)

    print("> 3D point estimate triangulation....")
    # Step 4: triangulate points to get an initial estimate of the 3D point cloud
    feat3D = rec_engine.get_initial_structure(norm_feat_2d, P, epipole_1, epipole_2, F)

    print("> Estimating projection matrices for additional views....")
    # Step 5: Use the 3D point estimates to estimate the projection matrices of the remaining views
    P = rec_engine.projective_pose_estimation(norm_feat_2d, P, feat3D)

    print("> Bundle Adjustment....")
    # Step 5: Optimize 3D points and projection matrices using the reprojection error
    P, feat3D, error = rec_engine.bundle_adjustment(norm_feat_2d, P, feat3D)
    print("  - Bundle adjustment error: ", error)

    print("> Self-calibration")
    # Step 6: Self-calibration
    Tm, K, error = rec_engine.self_calibration(P)
    print("  - Self-calibration error: ", error)
    print("  - Tranformation Matrix (Projective -> Metric): ")

    print("> Converting to metric space")
    metric_feat3D, metric_P = rec_engine.convert_to_metric_space(Tm, feat3D, P, K)

    print("> Saving model...")

    recModel = RecModel()
    recModel.P = P
    recModel.points3D = metric_feat3D
    recModel.Tm = Tm

    np.savetxt('rec_model_cloud.txt', metric_feat3D, delimiter=',')
    print("	 - 3D point cloud saved in rec_model_cloud.txt ")

    recModel.export_stl_file('reconstructed_model.stl')
    print("  - STL model saved in reconstructed_model.stl")

    if (show == True):
        sequence.show()

    print("> 3D reconstruction completed in " + str(round(time.time() - start_time, 1)) + " sec!")

if __name__ == "__main__":
    parser = argparse.ArgumentParser(description='3D Reconstruction from uncalibrated images')
    parser.add_argument('--input_file', metavar='path', required=True,
                        help='Input file containing image point correspondences')
    parser.add_argument('--show', required=False, action="store_true",
                        help="Display the image sequence with the 2D features")
    args = parser.parse_args()
    main(input_file=args.input_file, show=args.show)

reason: import numpy,from numpy,from scipy


text:
import numpy as np
import tensorflow as tf

from tensorstream.tests import TestCase
from tensorstream.finance.moving_standard_deviation import MovingStandardDeviation

class MovingStandardDeviationSpec(TestCase):
  # Data from http://stockcharts.com/school/doku.php?id=chart_school:technical_indicators:standard_deviation_volatility
  def setUp(self):
    self.sheets = self.read_ods(
      self.from_test_res('moving_standard_deviation.ods', __file__))

  def test_moving_standard_deviation_single_dim(self):
    volatility = MovingStandardDeviation(10)
    values = tf.placeholder(tf.float32)
    volatility_ts, _, _ = volatility(values)

    inputs = self.sheets['Sheet1']

    with tf.Session() as sess:
      output = sess.run(volatility_ts, { values: inputs['Value'] })

    np.testing.assert_almost_equal(output,
      inputs['mstdev'].values, decimal=3)


text:
import Monsoon.HVPM as HVPM
import Monsoon.LVPM as LVPM
import Monsoon.sampleEngine as sampleEngine
import Monsoon.Operations as op
import Monsoon.pmapi as pmapi
import numpy as np

def testHVPM(serialno=None,Protocol=pmapi.USB_protocol()):
    HVMON = HVPM.Monsoon()
    HVMON.setup_usb(serialno,Protocol)
    print("HVPM Serial Number: " + repr(HVMON.getSerialNumber()))
    HVMON.fillStatusPacket()
    HVMON.setVout(3)
    HVengine = sampleEngine.SampleEngine(HVMON)
    #Output to CSV
    HVengine.enableCSVOutput("HV Main Example.csv")
    #Turning off periodic console outputs.
    HVengine.ConsoleOutput(True)

    #Setting all channels enabled
    HVengine.enableChannel(sampleEngine.channels.MainCurrent)
    HVengine.enableChannel(sampleEngine.channels.MainVoltage)
    HVengine.enableChannel(sampleEngine.channels.USBCurrent)
    HVengine.enableChannel(sampleEngine.channels.USBVoltage)
    HVengine.enableChannel(sampleEngine.channels.AuxCurrent)
    HVengine.enableChannel(sampleEngine.channels.timeStamp)

    #Setting trigger conditions
    numSamples=sampleEngine.triggers.SAMPLECOUNT_INFINITE
    HVengine.setStartTrigger(sampleEngine.triggers.GREATER_THAN,0)
    HVengine.setStopTrigger(sampleEngine.triggers.GREATER_THAN,5)
    HVengine.setTriggerChannel(sampleEngine.channels.timeStamp)

    #Actually start collecting samples
    HVengine.startSampling(numSamples)
    #startSampling() continues until the trigger conditions have been met, and then ends automatically.
    #Measurements are automatically saved to the filename passed in enableCSVOutput()

    #Disable CSV Output
    HVengine.disableCSVOutput()
    #Collect another 5 seconds worth of samples
    HVengine.startSampling(numSamples)
    #Get those samples as a Python list
    samples = HVengine.getSamples()
    #Samples has the format  [[timestamp], [mainCurrent], [usbCurrent], [auxCurrent], [mainVolts],[usbVolts]]
    #Use sampleEngine.channel to select the appropriate list index.
    timestamp = samples[sampleEngine.channels.timeStamp]
    mainCurrent = samples[sampleEngine.channels.MainCurrent]
    auxCurrent = samples[sampleEngine.channels.AuxCurrent]
    usbCurrent = samples[sampleEngine.channels.USBCurrent]
    mainVoltage = samples[sampleEngine.channels.MainVoltage]
    usbVoltage = samples[sampleEngine.channels.USBVoltage]

    #Perform analysis on the resulting data.  For example, in order to calculate, perform the following:
    #mainCurrent is given in mA.  Divide by 1000 to convert to Amps
    scaledMainCurrent = [x / 1000 for x in mainCurrent]
    #Element-wise multiply to produce Watts. Power = Current * Voltage.
    mainPower = np.multiply(scaledMainCurrent, mainVoltage)

    #When finished, close the device.
    HVMON.closeDevice();

def testLVPM(serialno=None,Protcol=pmapi.USB_protocol()):
    Mon = LVPM.Monsoon()
    Mon.setup_usb(serialno,Protcol)
    print("LVPM Serial number: " + repr(Mon.getSerialNumber()))
    Mon.fillStatusPacket()
    Mon.setVout(4.5)
    engine = sampleEngine.SampleEngine(Mon)
    engine.enableCSVOutput("Main Example.csv")
    engine.ConsoleOutput(True)
    #test main channels
    numSamples=sampleEngine.triggers.SAMPLECOUNT_INFINITE #Don't stop based on sample count, continue until the trigger conditions have been satisfied.
    engine.setStartTrigger(sampleEngine.triggers.GREATER_THAN,0) #Start when we exceed 0 s
    engine.setStopTrigger(sampleEngine.triggers.GREATER_THAN,5) #Stop when we exceed 5 s.
    engine.setTriggerChannel(sampleEngine.channels.timeStamp) #Start and stop judged by the timestamp channel.
    engine.startSampling(numSamples)

    #Disable Main channels
    engine.disableChannel(sampleEngine.channels.MainCurrent)
    engine.disableChannel(sampleEngine.channels.MainVoltage)

    engine.setStartTrigger(sampleEngine.triggers.GREATER_THAN,5)
    engine.setStopTrigger(sampleEngine.triggers.GREATER_THAN,10)
    engine.setTriggerChannel(sampleEngine.channels.timeStamp)
    #Take measurements from the USB Channel
    Mon.setVout(0)
    #Set USB Passthrough mode to 'on,' since it defaults to 'auto' and will turn off when sampling mode begins.
    Mon.setUSBPassthroughMode(op.USB_Passthrough.On)
    #Enable USB channels
    engine.enableChannel(sampleEngine.channels.USBCurrent)
    engine.enableChannel(sampleEngine.channels.USBVoltage)
    engine.enableCSVOutput("USB Test.csv")
    engine.startSampling(5000)

    #Enable every channel, take measurements
    engine.enableChannel(sampleEngine.channels.MainVoltage)
    engine.enableChannel(sampleEngine.channels.MainCurrent)
    #Enable Aux channel
    engine.enableChannel(sampleEngine.channels.AuxCurrent)
    Mon.setVout(2.5)
    engine.enableCSVOutput("All Test.csv")
    engine.startSampling(5000)

    #Enable every channel, take measurements, and retrieve them as a Python list.
    engine.disableCSVOutput()
    engine.startSampling(5000)
    samples = engine.getSamples()
    Mon.closeDevice();

def main():
    LVPMSerialNo = 60001
    HVPMSerialNo = 60000
    testLVPM(LVPMSerialNo,pmapi.USB_protocol())
    testHVPM(HVPMSerialNo,pmapi.USB_protocol())

if __name__ == "__main__":
    main()


text:
import xarray as xr
from astropy.io import fits

@xr.register_dataarray_accessor("cube")
class CubeAccessor:
    def __init__(self, xarray_obj: xr.DataArray) -> None:
        self._obj = xarray_obj

    @property
    def header(self) -> fits.header.Header:
        """Convert the header string to a header object."""
        return fits.Header.fromstring(self._obj.header.item())

max_stars_repo_name: a-lab-nagoya/xarray-cube

max_issues_repo_name: a-lab-nagoya/xarray-cube

max_forks_repo_name: a-lab-nagoya/xarray-cube

reason: from astropy


text:
# -*- coding: utf-8 -*-
"""
Created on Mon Jul 17 15:28:20 2017

@author: daniel
"""

#import Tomography as tom
#import quPy as qp
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import os
import json
import io
import scipy.constants as co
#import rydpy

c = 299792458 # m/s, speed of light CODATA 2014
a0 = 0.52917721067e-10 # m, Bohr radius
C6 = 2.3e23 * 4.36e-18 * a0**6 # Jm^6, Van-der-Waals coefficient for the 67s - 69s
hbar = 6.626070040e-34/(2 * np.pi) # Js, Planck constant, CODATA 2014
rho_peak = 2.4e12/1e-6 # peak density in cm^-3/centi^-3
d = 2.534e-29 # Cm, dipole matrix element (D. A. Steck)
Gamma_e = 2*np.pi * 6.065e6 # decay rate (D. A. Steck)
epsilon_0 = 8.854187817e-12 # dielectric constant, CODATA 2014
L = 61e-6 # medium length in m
omega_s = 2*np.pi * 384.23e12 # rad/s, transition frequency
gamma_21 = 0.0577314
chi_0 = 2*rho_peak*d**2 / (epsilon_0*hbar*Gamma_e) # prefactor of the susceptibility for the cycling transition (|R> polarization)
R_b=18e-6
n1=69
n2=100
ec=co.e
m_e=co.electron_mass

def classicalRadius(n,l=0,j=0.5):
    return (ec**2/(4*np.pi*epsilon_0*co.h*rydpy.IonEnergy(n,l,j,units="Hz")))

def orbitalRadius(n,l=0,j=0.5):
    r,rrRR,_,_=rydpy.Numerov(n,l,j)
    imax=np.argmax(rrRR)
    return r[imax]*a0

#rcl1=classicalRadius(n1)
rcl1=orbitalRadius(n1)
Vcl1=4/3.*np.pi*np.abs(rcl1)**3

rcl2=orbitalRadius(n2)
Vcl2=4/3.*np.pi*np.abs(rcl2)**3

rBend=(8*np.pi*epsilon_0*hbar**2*(n1-rydpy.QuantumDefect(n1,0,0.5))**2)/(1*ec**2*m_e)
print rBend/a0

#plt.plot(r,rrRR)
#plt.axvline(np.abs(r[imax]))

print "Average number of atoms inside Orbital for n=69"
print rho_peak*Vcl1

rhos=np.linspace(0,3e12/1e-6,10)

plot_dict={}

# plot it
f=plt.figure()

#h=pd.DataFrame(index=detuning,data=np.imag(chi))
#h2=pd.DataFrame(index=detuning,data=np.imag(chi_2L))
#plot_dict['221']={
#    'A':{'type':'plot','y':h[0].to_json(),'ylabel':u'Im($\chi / \chi_0$)','margin':(0,0.0),'num':'a'},
#    'B':{'type':'plot','y':h2[0].to_json(),'margin':(0,0.0)},
#    'C':{'type':'axv','y':0}
#}

plt.subplot(111)
plt.ylabel(u'$\Delta \phi$')
plt.plot(rhos/1e18,rhos*Vcl1)
plt.plot(rhos/1e18,rhos*Vcl2)
#plt.xscale('log')
plt.axhline(1)
#plt.yticks(np.arange(-0.4, 0.6, 0.2))
#plt.margins(0,0.1)
#plt.xlabel(u'$\varrho$ $(1e12\ cm^-3)$')
#plt.tight_layout()

h=pd.DataFrame(index=rhos/1e18,data=rhos*Vcl1)
h2=pd.DataFrame(index=rhos/1e18,data=rhos*Vcl2)
plot_dict['111']={
    'A':{'type':'plot','y':h[0].to_json(),'ylabel':r'$\left\langle N_r\right\rangle$','xlabel':r'$\varrho$ $(10^{12}\ \mathrm{cm}^{-3})$','label':'$n=69$'},
    'B':{'type':'plot','y':h2[0].to_json(),'label':'$n=100$'},
    'C':{'type':'axh','y':1}
}

#with io.open('avg_number.json', 'w+') as f:
#  f.write(unicode(json.dumps(plot_dict, ensure_ascii=False,indent=4)))

plt.show()

max_stars_repo_name: dtiarks/ThesisPlot

max_issues_repo_name: dtiarks/ThesisPlot

max_forks_repo_name: dtiarks/ThesisPlot

reason: import numpy,import scipy


text:
import numpy as np

# Problem 1
# print list of list as a np matrix
def print_matrix(M_lol):
	M = np.array(M_lol)
	print(M)

# Problem 2
# get index of first non-zero element in list row
def get_lead_ind(row):
	for i in range(len(row)):
		if row[i] != 0:
			return i
	return len(row)

# Problem 3
# get row with leading coefficient to the far left
def get_row_to_swap(M, start_i):
	val = []
	for i in range(start_i, len(M)):
		val.append(get_lead_ind(M[i]))
	return val.index(min(val)) + start_i

# Problem 4
# row operations
def add_rows_coefs(r1, c1, r2, c2):
	new_row = [0] * len(r1)
	for i in range(len(r1)):
		new_row[i] = r1[i] * c1 + r2[i] * c2
	return new_row

# Problem 5
# eliminate a coefficient with row operations
def eliminate(M, row_to_sub, best_lead_ind):
	for i in range(row_to_sub+1, len(M)):
		temp_row = add_rows_coefs(M[row_to_sub],-1*M[i][best_lead_ind],M[i],M[row_to_sub][best_lead_ind])
		M[i] = temp_row

# Problem 6
def forward_step(M):
	print("The matrix is currently:")
	print_matrix(M)
	print("Now performing the forward step")

	for row in range(len(M)):
		print("------------------------------------------------------------")
		print("Now looking at row " + str(row))
		n = get_row_to_swap(M, row)
		entry = get_lead_ind(M[n])

		print("Swapping rows %i and %i so that entry %i in the current row is non-zero" % (row,n,entry))
		M[row], M[n] = M[n], M[row]
		print("The matrix is currently:")
		print_matrix(M)

		print("Adding row %i to rows below it to eliminate coefficients in column %i" % (row,entry))
		eliminate(M, row, entry)
		print("The matrix is currently:")
		print_matrix(M)

	print("Done with the forward step")
	print("The matrix is currently:")
	print_matrix(M)

# Problem 7
def backward_step(M):
	print("The matrix is currently:")
	print_matrix(M)
	print("Now performing the backward step")

	new_M = M[::-1]

	for row in range(len(new_M)):
		print("------------------------------------------------------------")
		print("Now looking at row " + str(len(M) - 1 - row))
		eliminate(new_M, row, len(M) - 1 - row)
		print("The matrix is currently:")
		print_matrix(new_M[::-1])

	M = new_M[::-1]

	print("Now dividing each row by the leading coefficient")
	for i in range(len(M)):
		n = get_lead_ind(M[i])
		M[i] = [x/M[i][n] for x in M[i]]
	print("The matrix is currently:")
	print_matrix(M)

	# Return solutions as a list
	sol = []
	for i in M:
		sol.append(i[-1])
	return sol

# Problem 8 - Testing
# For Mx=b
# Given M
# Given b
# Find x
M = [[1, -2, 3, 22 ],
     [3, 10, 1, 314],
     [1, 5 , 3, 92 ],
     [1, 2 , 3, 4  ]]

b = [6,4,3,6]

def solve(M, b):
	for row in range(len(M)):
		M[row].append(b[row])
	forward_step(M)
	sol = backward_step(M)
	return sol
sol = solve(M, b)
print()
print("THE SOLUTION IS:")
print(sol)
print()

# Verification
# Given M
# Given x
# Check b
# Testing whether or not Mx = b
print("Verifying the input value of b:")
M = np.array([[1, -2, 3, 22],[3, 10, 1, 314],[1, 5, 3, 92],[1,2,3,4]])
x = np.array(sol)
b = np.matmul(M,x)
print(b)


text:
#VideoWriter from Alexander Mordvintsev
#https://colab.research.google.com/github/znah/notebooks/blob/master/external_colab_snippets.ipynb

import os
import numpy as np
os.environ['FFMPEG_BINARY'] = 'ffmpeg'
import moviepy.editor as mvp
from moviepy.video.io.ffmpeg_writer import FFMPEG_VideoWriter

class VideoWriter:
  def __init__(self, filename='_autoplay.mp4', fps=30.0, **kw):
    self.writer = None
    self.params = dict(filename=filename, fps=fps, **kw)

  def add(self, img):
    img = np.asarray(img)
    if self.writer is None:
      h, w = img.shape[:2]
      self.writer = FFMPEG_VideoWriter(size=(w, h), **self.params)
    if img.dtype in [np.float32, np.float64]:
      img = np.uint8(img.clip(0, 1)*255)
    if len(img.shape) == 2:
      img = np.repeat(img[..., None], 3, -1)
    self.writer.write_frame(img)

  def close(self):
    if self.writer:
      self.writer.close()

  def __enter__(self):
    return self

  def __exit__(self, *kw):
    self.close()
    if self.params['filename'] == '_autoplay.mp4':
      self.show()

  def show(self, **kw):
      self.close()
      fn = self.params['filename']
      display(mvp.ipython_display(fn, **kw))

max_stars_repo_name: PWhiddy/notebook-video-writer

max_issues_repo_name: PWhiddy/notebook-video-writer

max_forks_repo_name: PWhiddy/notebook-video-writer