text:
import numpy as np
import matplotlib.pyplot as plt

from ccblade import CCAirfoil, CCBlade
from openmdao.api import Component, Group, Problem
from scipy.optimize import optimize

class CCBladeEvaluator(Component):

    def __init__(self, config, sdim, nalpha, nv, naero_coeffs):
        super(CCBladeEvaluator, self).__init__()

        self.sdim = sdim
        self.nsec = sdim[1]

        self.nalpha = nalpha

        self.nv = nv

        # standard values for cp lambda curves
        self.tsr_cp_start = 1
        self.tsr_cp_end = 18
        self.ntsr_cp = 100
        self.pitches_cp_start = 0
        self.pitches_cp_end = 30
        self.npitches_cp = 10

        self.use_pitch_oper = False
        self.pitch_oper = np.array([])

        #self.blend_var = config['CCBlade']['blend_var']

        self.af_name_base = 'cs_'
        self.af_name_suffix = '_aerodyn'

        self.naero_coeffs = naero_coeffs

        for k, w in config['CCBlade'].iteritems():
            try:
                setattr(self, k, w)
            except:
                pass

        self.nre = len(config['CCBlade']['res'])
        self.nmet = len(config['CCBlade']['analysis_methods'])

        self._init_params()
        self._init_unknowns()

    def _init_params(self):
        self.add_param('blade_length', 0.)

        self.add_param('rotor_diameter', 0.)

        # add planform
        self.add_param('s_st', np.zeros(self.nsec))
        self.add_param('x_st', np.zeros(self.nsec))
        self.add_param('y_st', np.zeros(self.nsec))
        self.add_param('z_st', np.zeros(self.nsec))
        self.add_param('rot_x_st', np.zeros(self.nsec))
        self.add_param('rot_y_st', np.zeros(self.nsec))
        self.add_param('rot_z_st', np.zeros(self.nsec))
        self.add_param('chord_st', np.zeros(self.nsec))
        self.add_param('rthick_st', np.zeros(self.nsec))
        self.add_param('p_le_st', np.zeros(self.nsec))

        self.add_param('cs_polars', np.zeros(  # nsec, alpha, cl ,cd, cm
            (self.nalpha, 4, self.nsec, self.nre, self.nmet)))
        self.add_param(  # amount of polar entries per table
            'n_cs_alpha', np.zeros((self.nsec, self.nre, self.nmet)).astype(int))
        self.add_param('cs_polars_tc', np.zeros((self.nsec)))

        self.add_param('airfoildata:blend_var', np.zeros(self.naero_coeffs))

    def _init_unknowns(self):
        self.add_output('pitch_v', np.zeros((self.nv)))
        self.add_output('aeloads_v', np.zeros((self.nsec, 2, self.nv)))
        self.add_output('omega_v', np.zeros((self.nv)))
        self.add_output('cp_v', np.zeros((self.nv)))
        self.add_output('tsr_v', np.zeros((self.nv)))
        self.add_output('v', np.zeros((self.nv)))
        self.add_output('P_v', np.zeros((self.nv)))

        self.add_output('load_cases_aero', np.zeros((self.nv, self.nsec, 15)))
        self.add_output('load_cases_gravity', np.zeros((self.nv, 1, 15)))
        self.add_output('load_cases_align', np.zeros((self.nv, 5)))

        self.add_output('r_stat', np.zeros((self.nsec)))

        #self.add_output('load_cases_aero_ult', np.zeros((2, self.nsec, 15)))
        #self.add_output('load_cases_align_ult', np.zeros((2, 5)))

        self.add_output('tsr_cp', np.zeros((self.ntsr_cp)))
        self.add_output('CP_cp', np.zeros((self.ntsr_cp, self.npitches_cp)))
        self.add_output('CT_cp', np.zeros((self.ntsr_cp, self.npitches_cp)))
        self.add_output('CQ_cp', np.zeros((self.ntsr_cp, self.npitches_cp)))
        self.add_output('pitches_cp', np.zeros((self.npitches_cp)))

    def solve_nonlinear(self, params, unknowns, resids):

        self.Rtip = 0.5 * params['rotor_diameter']
        self.Rhub = self.Rtip - params['blade_length']
        z = params['z_st'] * params['blade_length']
        self.r = z + self.Rhub

        # First and last section need to be shifted if they coincide with the
        # root/tip
        dr01 = self.r[1] - self.r[0]
        rnode1 = self.r[0] + 0.5 * dr01
        rnode0 = self.r[0]
        self.r[0] = 0.5 * (rnode1 + rnode0)  # override root value

        dr_2_1 = self.r[-1] - self.r[-2]
        rnode_2 = self.r[-2] + 0.5 * dr_2_1
        rnode_1 = self.r[-1]
        self.r[-1] = 0.5 * (rnode_1 + rnode_2)  # override tip value

        unknowns['r_stat'] = self.r

        # Workaround to fix nan calculated by CCblade
        # obsolete
        # self.r[0] = self.Rhub + 1.0E-11
        # self.r[-1] = self.Rtip - 1.0E-11

        rot_z = params['rot_z_st']
        # twist angle
        self.theta = - rot_z

        #file_planform = '/home/becmax/iwes-gitlab/bdtprojects/smartblades/btc_design_model/data/turbine/aero/SmartBlades-Ref12 Aerodynamic Design 1.txt'
        # planform_data = np.loadtxt(
        #    file_planform, skiprows=13, usecols=(0, 1, 3, 4))
        self.chord = params['chord_st'] * params['blade_length']

        self.blend_var = params['airfoildata:blend_var']
        self.af = []
        # init from aerodyn files
        for i, tc in enumerate(self.blend_var):

            af_name = self.af_name_base + \
                '%03d_%04d' % (i, tc * 1000) + self.af_name_suffix

            self.af.append(CCAirfoil.initFromAerodynFile(af_name + '.dat'))

        # alternative input
        '''
        # load all airfoils
        pol = params['cs_polars']  # prob.root.unknowns['af_data']

        naf = len(np.trim_zeros(params['cs_polars_tc']))  # len(af_list)

        #nsec = 0
        #nre = 0

        nmet = 0
        print 'CCBLade uses polars generated with analysis method: %s ' % self.analysis_methods[nmet]

        #cs_polars_met = p['cs_polars_met']
        cs_polars_tc = params['cs_polars_tc']
        #cs_polars_re = p['cs_polars_re']
        n_cs_alpha = params['n_cs_alpha']

        # convert to CCAirfoil

        for nsec in range(naf):  # [0, 1, 2, 3, 4, 5, 6, 7]:  #
            re = self.res  # Reynolds numbers where airfoil data are defined
            mins = []
            maxs = []
            for rey in range(self.nre):
                mins.append(np.min(pol[:, 0, nsec, rey, nmet]))
                maxs.append(np.max(pol[:, 0, nsec, rey, nmet]))

            alpha_min = max(mins)
            alpha_max = min(maxs)

            alpha_incr = pol[1, 0, nsec, 0, nmet] - pol[0, 0, nsec, 0, nmet]
            alpha = np.arange(alpha_min, alpha_max + alpha_incr, alpha_incr)

            cl = np.zeros((len(alpha), self.nre))
            cd = np.zeros((len(alpha), self.nre))
            for rey in range(self.nre):
                # cl[i, j] is the lift coefficient at alpha[i] and Re[j]
                cl[:, rey] = np.interp(alpha, pol[:n_cs_alpha[nsec, rey, nmet], 0, nsec, rey, nmet], pol[
                    :n_cs_alpha[nsec, rey, nmet], 1, nsec, rey, nmet])
                # cd[i, j] is the lift coefficient at alpha[i] and Re[j]
                cd[:, rey] = np.interp(alpha, pol[:n_cs_alpha[nsec, rey, nmet], 0, nsec, rey, nmet], pol[
                    :n_cs_alpha[nsec, rey, nmet], 2, nsec, rey, nmet])
            self.af.append(CCAirfoil(alpha, re, cl, cd))
        '''
        # create CCBlade object
        self.rotor = CCBlade(self.r, self.chord, self.theta, self.af, self.Rhub, self.Rtip, self.B, self.rho, self.mu, self.precone, self.tilt, self.yaw, self.shearExp, self.hubHt, self.nSector,
                             presweepTip=0.0, wakerotation=self.wakerotation, usecd=self.usecd, tiploss=self.tiploss, hubloss=self.hubloss, iterRe=self.iterRe)

        # find tsr opt
        cp_opt, tsr_opt, tsr_cp, CP_cp, CT_cp, CQ_cp, pitches_cp = self.evaluate_cp_lambda()
        lambda_opt = tsr_opt[0]  # at pitch zero

        # create wind speed vector
        #U_oper = np.linspace(self.Uin, self.Uout, self.nv)
        U_oper = np.linspace(self.Ustart_eval, self.Ustop_eval, self.nv)

        # create operating conditions
        pitch_angles, cp, aeloads, tsr, omega, P = self.evaluate_steady_states(
            P_elec=self.P_elec, eta_elec=self.eta_elec, lambda_opt=lambda_opt, OmegaMin=self.OmegaMin, OmegaMax=self.OmegaMax, U_oper=U_oper, use_pitch_oper=self.use_pitch_oper, pitch_oper=self.pitch_oper)
        # self.pp.close()

        # get u_rated
        self.diameter = params['rotor_diameter']

        cp_max_idx = cp.argsort()[-1]
        Urated = U_oper[cp_max_idx]
        omega_rated = omega[cp_max_idx]
        pitch_rated = pitch_angles[cp_max_idx]

        # run ultimate load cases
        # DLC 2.3 EOG v_hub = v_r+-2m/s and v_out Ultimate Abnormal
        # acc. to IEC64100-1 ed3 Eq. 17
        #V_hub = Urated
        #V_e1 = 0.8 * V_e50
        #V_g1 = 1.35 * (V_e1 - V_hub)
        #D = self.diameter

        # IEC64100-1 ed Tab.1
        #Iref = self.i_ref

        # IEC64100-1 ed3 Eq. 11
        #b = 5.6
        #sigma1 = Iref * (0.75 * V_hub + b)

        # IEC64100-1 ed3 Eq. 5
        # lambda1 = 0.7*h_hub # z <60m
        # lambda1 = 42  # z <=60m
        #V_g2 = 3.3 * (sigma1 / (1 + 0.1 * (D / lambda1)))
        #V_gust = np.min(np.array([V_g1, V_g2]))

        # DLC 6.1 EWM 50y recurrence period +- 8deg yaw
        shearExp_ult = 0.11

        self.rotor_ult = CCBlade(self.r, self.chord, self.theta, self.af, self.Rhub, self.Rtip, self.B, self.rho, self.mu, self.precone, self.tilt, self.yaw, shearExp_ult, self.hubHt, self.nSector,
                                 presweepTip=0.0, wakerotation=self.wakerotation, usecd=self.usecd, tiploss=self.tiploss, hubloss=self.hubloss, iterRe=self.iterRe)

        V_ref = self.v_ref
        V_e50 = 1.4 * V_ref
        pitch_dlc61 = 90
        yaw_errors = [+8, -8]

        aeloads_ult = self.evaluate_ultimate_loads(
            V_e50, pitch_dlc61, yaw_errors)

        # shift data to unknowns
        unknowns['pitch_v'] = - pitch_angles  # convert positive along z axis
        unknowns['aeloads_v'] = aeloads
        unknowns['omega_v'] = omega
        unknowns['cp_v'] = cp
        unknowns['tsr_v'] = tsr
        unknowns['v'] = U_oper
        unknowns['P_v'] = P

        # convert to FEPROC load case which need to be applied on t/4 point
        # x,y,z,rotx,roty,rotz,Fx,Fy,Fz,Mx,My,Mz,ax,ay,az
        # 'load_cases', np.zeros((self.nv, self.nsec, 15)))

        # tilt, pitch, x_orig, y_orig, z_orig
        #'load_cases_align', np.zeros((self.nv, 5)))

        for i, _ in enumerate(U_oper):
            unknowns['load_cases_aero'][i, :, 2] = params[
                'z_st'] * params['blade_length']
            unknowns['load_cases_aero'][
                i, :, 7] = aeloads[:, 0, i]  # flap-wise
            unknowns['load_cases_aero'][i, :, 6] = aeloads[:, 1, i]  # lead-lag

            unknowns['load_cases_align'][
                i, 1] = - pitch_angles[i]  # positive along z

            # load cases for lead-lag gravity (x-direction)
            unknowns['load_cases_gravity'][i, 0, :] = np.array(
                    [[0, 0, 0,
                      0, 0, 0,
                      9.81, 0.0, 0.0,  # g_x, g_y, gz
                      0.0, 0.0, 0.0,
                      0.0, 0.0, 0.0]])

        # for i, _ in enumerate(yaw_errors):
        #    unknowns['load_cases_aero_ult'][i, :, 2] = params[
         #       'z_st'] * params['blade_length']
         #   unknowns['load_cases_aero_ult'][
         #       i, :, 7] = aeloads_ult[:, 0, i]  # flap-wise
        #    unknowns['load_cases_aero_ult'][
        #        i, :, 6] = aeloads_ult[:, 1, i]  # lead-lag

        #    unknowns['load_cases_align_ult'][
        #        i, 1] = - pitch_dlc61  # positive along z
        unknowns['tsr_cp'] = tsr_cp
        unknowns['CP_cp'] = CP_cp
        unknowns['CT_cp'] = CT_cp
        unknowns['CQ_cp'] = CQ_cp
        unknowns['pitches_cp'] = - pitches_cp  # convert positive along z axis

    def evaluate_cp_lambda(self):
        tsr = np.linspace(self.tsr_cp_start, self.tsr_cp_end, self.ntsr_cp)
        Omega = 10.0 * np.ones_like(tsr)
        Uinf = Omega * np.pi / 30.0 * self.Rtip / tsr

        # pitch = np.zeros_like(tsr)

        pitches = np.linspace(self.pitches_cp_start,
                              self.pitches_cp_end, self.npitches_cp)

        cp_opt = np.zeros_like(pitches)
        tsr_opt = np.zeros_like(pitches)
        CP = np.zeros((len(tsr), len(pitches)))
        CT = np.zeros((len(tsr), len(pitches)))
        CQ = np.zeros((len(tsr), len(pitches)))
        for i, pit in enumerate(pitches):
            pitch = pit * np.ones_like(tsr)
            CP[:, i], CT[:, i], CQ[:, i] = self.rotor.evaluate(
                Uinf, Omega, pitch, coefficient=True)

            cp_opt[i] = np.max(CP[:, i])
            tsr_opt[i] = tsr[CP[:, i].argmax()]

        #filename_excelcoeff = '/home/becmax/iwes-gitlab/bdtprojects/smartblades/btc_design_model/reference_excel_coeff.txt'
        #excel_coeff = np.loadtxt(filename_excelcoeff)

        #TSR_excel = excel_coeff[:, 0]
        #CP_excel = excel_coeff[:, 1]

        return cp_opt, tsr_opt, tsr, CP, CT, CQ, pitches

    def evaluate_partial_load(self, lambda_opt, Urated, OmegaMin, OmegaMax, eta_elec):

        U_vector = np.arange(0, Urated, 0.1)
        Omega_vector = np.zeros_like(U_vector)

        #lambda_opt = 8.4

        for i in range(len(U_vector)):

            Omega = U_vector[i] * 30.0 / np.pi * lambda_opt / self.Rtip
            if Omega < OmegaMin:
                Omega = OmegaMin
            elif Omega > OmegaMax:
                Omega = OmegaMax
            Omega_vector[i] = Omega

        pitch_vector = np.zeros_like(U_vector)
        P, T, Q = self.rotor.evaluate(U_vector, Omega_vector, pitch_vector)

        #filename_excel = '/home/becmax/iwes-gitlab/bdtprojects/smartblades/btc_design_model/reference_excel.txt'
        #excel_data = np.loadtxt(filename_excel)

        #P_excel = excel_data[:, 2]
        #U_excel = excel_data[:, 0]

        plt.figure()
        plt.plot(U_vector, P / eta_elec, label='CCblade')
        #plt.plot(U_excel, P_excel, label='Reference')
        plt.xlabel('Wind speed (m/s)')
        plt.ylabel('Electric power (kW)')
        plt.legend(loc='best')
        plt.grid()
        # self.pp.savefig()

    def evaluate_ultimate_loads(self, V_e50, pitch_dlc61, yaw_errors):

        aeloads_ult = np.zeros((len(self.r), 2, len(yaw_errors)))
        m_tan = []
        m_norm = []
        m_res = []
        for i, yaw_error in enumerate(yaw_errors):
            aeloads_ult[:, 0, i], aeloads_ult[:, 1, i] = self.rotor_ult.distributedAeroLoads(
                V_e50, Omega=0.000, pitch=pitch_dlc61 + yaw_error, azimuth=0)

            #dr = self.r[1] - self.r[0]
            #mn = np.sum(Np * dr * self.r)
            #mt = np.sum(Tp * dr * self.r)
            # gf = 1.35  # normal
            #m_tan.append(-mt * gf / 1E06)
            #m_norm.append(mn * gf / 1E06)
            #mres = np.sqrt(mn**2 + mt**2) * gf / 1E06
            # m_res.append(mres)

        # print m_tan
        # print m_norm
        # print m_res

        return aeloads_ult

    def evaluate_steady_states(self, P_elec, eta_elec, lambda_opt, OmegaMin, OmegaMax, U_oper, use_pitch_oper, pitch_oper):

        Uin = self.Uin  # U_oper[0]
        Uout = self.Uout  # U_oper[-1]
        Omega_oper = np.zeros_like(U_oper)
        for i, U in enumerate(U_oper):
            Omega = U_oper[i] * 30.0 / np.pi * lambda_opt / self.Rtip
            if Omega < OmegaMin:
                Omega = OmegaMin
            elif Omega > OmegaMax:
                Omega = OmegaMax

            if U < Uin or U > Uout:
                Omega = 0

            Omega_oper[i] = Omega

        #Omega = 10.0
        #U = 15.0
        P_rated = P_elec / eta_elec
        # start values
        pitch_init = 1.0  # closest

        pitch_angles = np.zeros_like(U_oper)
        cp = np.zeros_like(U_oper)
        P = np.zeros_like(U_oper)
        aeloads = np.zeros((len(self.r), 2, len(U_oper)))
        #Np = np.zeros_like(U_oper)
        #Tp = np.zeros_like(U_oper)
        tsr = np.zeros_like(U_oper)

        for i, (U, Omega) in enumerate(zip(U_oper, Omega_oper)):

            P_goal = P_rated
            if U < Uin or U > Uout:
                P_goal = 0

            if not use_pitch_oper:
                # pitch_init2 = pitch_init - 1.0  # second closest
                pitch = pitch_init

                def min_func(pitch):
                    ''' function to be minimized
                    '''
                    """
                    U_p_tsr_opt = Omega * np.pi / 30.0 * self.Rtip / lambda_opt
                    P_tsr_opt, _, _ = self.rotor.evaluate(
                        [U_p_tsr_opt], [Omega], [pitch])
                    if P_tsr_opt < P_rated:
                        P_goal = P_tsr_opt
                    """
                    P, _, _ = self.rotor.evaluate([U], [Omega], [pitch])
                    P_diff = P - P_goal
                    return abs(P_diff)

                # do not pitch when omega is controlled
                '''
                if Omega < OmegaMax:
                    pitch = 0.
                else:
                '''
                pitch = optimize.brent(
                    min_func, brack=(pitch_init, pitch_init + 1))
            else:
                pitch = pitch_oper[i]

            cp[i], _, _ = self.rotor.evaluate(
                [U], [Omega], [pitch], coefficient=True)
            P[i], _, _ = self.rotor.evaluate(
                [U], [Omega], [pitch], coefficient=False)
            aeloads[:, 0, i], aeloads[:, 1, i] = self.rotor.distributedAeroLoads(
                U, Omega, pitch, azimuth=90.0)
            #aeloads[:, 0, i] = Np[i]
            #aeloads[:, 1, i] = Tp[i]
            tsr[i] = Omega / 60. * 2 * np.pi * self.Rtip / U

            pitch_init = abs(pitch)
            pitch_angles[i] = pitch

        plot = False
        if plot:
            plt.figure()
            plt.plot(U_oper, pitch_angles, label='CCblade')
            plt.plot(U_oper, Omega_oper, label='CCblade')
            plt.plot(U_oper, cp, label='CCblade')
            #plt.plot(TSR_excel, CP_excel, label='Reference simulated')
            plt.xlabel('$U$')
            plt.ylabel('$\beta$')
            plt.legend(loc='best')
            plt.grid()

        return pitch_angles, cp, aeloads, tsr, Omega_oper, P

    def evaluate_omega(self):
        print 'check'
        # first calculate correct CP by a lot of ccblade runs
        tsr = np.linspace(8.7, 8.8, 1000)
        Uinf = 6.0 * np.ones_like(tsr)
        Omega = Uinf * 30.0 / np.pi * tsr / self.Rtip
        pitch = np.zeros_like(tsr)
        CP, CT, CQ = self.rotor.evaluate(Uinf, Omega, pitch, coefficient=True)
        correctCP = max(CP)
        print 'check'
        # estimate CP using Newton-Raphson iteration using only a few ccblade
        # runs
        tsr = np.array([8.5, 8.75, 9.0])
        Uinf = 6.0 * np.ones_like(tsr)
        Omega = Uinf * 30.0 / np.pi * tsr / self.Rtip
        pitch = np.zeros_like(tsr)

        CP, CT, CQ = self.rotor.evaluate(Uinf, Omega, pitch, coefficient=True)
        CP_init = CP[1]
        CP_init2 = CP[0]
        CP_init3 = CP[2]
        tsr_init = tsr[1]
        tsr_init2 = tsr[0]
        tsr_init3 = tsr[2]
        tsr = tsr_init
        CP = CP_init

        CP_deriv = (CP_init - CP_init2) / (tsr_init - tsr_init2)
        CP_deriv_old = (CP_init2 - CP_init3) / (tsr_init2 - tsr_init3)
        CP_deriv_two = (CP_deriv - CP_deriv_old) / (tsr_init - tsr_init2)

        step = 0

        # Newton-Raphson iteration method according to
        # http://www.cup.uni-muenchen.de/ch/compchem/geom/nr.html
        # this parameter needs to be fixed to some proper value
        while abs(CP_deriv) >= 0.00005:
            tsr_old = tsr
            CP_old = CP
            CP_deriv_old = CP_deriv
            CP_deriv_two_old = CP_deriv_two

            # 0.05 is to reduce maximum stepsize (under relaxation factor)
            # script would diverge without
            tsr = tsr_old - 0.1 * (CP_deriv_old / CP_deriv_two_old)

            Omega = Uinf[0] * 30.0 / np.pi * tsr / self.Rtip

            CP, CT, CQ = self.rotor.evaluate(
                [Uinf[0]], [Omega], [0.0], coefficient=True)
            CP = CP[0]

            CP_deriv = (CP - CP_old) / (tsr - tsr_old)  # backward difference
            # backward difference
            CP_deriv_two = (CP_deriv - CP_deriv_old) / (tsr - tsr_old)

            CP_diff = CP - CP_old

            print step
            print tsr
            print Omega
            print CP_diff
            print CP
            step = step + 1

        if Omega < 5.0:
            Omega = 5.0
        if Omega > 10.0:
            Omega = 10.0

        print Omega
        print 'CP correct is:'
        print correctCP
        print 'CP difference is:'
        print (CP - correctCP)
        print 'CP relative difference is:'
        print 100 * (CP - correctCP) / correctCP

if __name__ == '__main__':

    cfg = {}
    cfg['path_airfoils'] = '/home/becmax/iwes-gitlab/bdtprojects/smartblades/btc_design_model/data/turbine/aero/geo_BlendedAirfoils_export'

    prob = Problem()
    root = prob.root = Group()
    #root.add('comp1', ReadAirfoils(cfg), promotes=['af_data'])
    root.add('comp2', CCBladeEvaluator(), promotes=['af_data'])
    prob.setup()

    print("num connections:", len(prob.root.connections))
    print("num unknowns:", len(prob.root._unknowns_dict),
          "size:", prob.root.unknowns.vec.size)
    print("num params:", len(prob.root._params_dict),
          "size:", prob.root.params.vec.size)

    prob.run()

max_stars_repo_name: mrosemeier/CCBlade

max_issues_repo_name: mrosemeier/CCBlade

max_forks_repo_name: mrosemeier/CCBlade

reason: import numpy,from scipy


text:
#!/usr/bin/env python
import os
import numpy as np
from tqdm import trange
import librenderman as rm
import scipy.io.wavfile
from deepdiva.utils.h2p_utils import H2P
from deepdiva.utils.patch_utils import split_train_override_patch

class DatasetGenerator():

    def __init__(self,
                 data_path,
                 vst_path,
                 folder_name,
                 base_preset,
                 random_parameters,
                 save_audio,
                 sample_rate,
                 midi_note_pitch,
                 midi_note_velocity,
                 note_length_seconds,
                 render_length_seconds):

        """
        Constructor
        """

        self.data_path = data_path
        self.vst_path = vst_path
        self.folder_name = folder_name
        self.base_preset = base_preset
        self.random_parameters = random_parameters
        self.save_audio = save_audio
        self.sample_rate = sample_rate
        self.midi_note_pitch = midi_note_pitch
        self.midi_note_velocity = midi_note_velocity
        self.note_length_seconds = note_length_seconds
        self.render_length_seconds = render_length_seconds

    def synth_configuration(self):
        engine = rm.RenderEngine(self.sample_rate, 512, 512)
        # Supply full path to the DIVA plugin to load the vst
        engine.load_plugin(self.vst_path)

        return engine

    def generate(self, sample_size, file_prefix):
        # Create relevant folders if not existing
        if not os.path.exists(self.data_path):
            os.makedirs(self.data_path)

        if not os.path.exists(os.path.join(self.data_path, self.folder_name)):
            os.makedirs(os.path.join(self.data_path, self.folder_name))

        if self.save_audio:
            if not os.path.exists(os.path.join(self.data_path, self.folder_name, "audio")):
                os.makedirs(os.path.join(self.data_path, self.folder_name, "audio"))

        # Set generator with synthesizer configuration
        engine = self.synth_configuration()
        generator = rm.PatchGenerator(engine)

        # Create one data example to infer size of arrays
        example_patch = generator.get_random_patch()
        engine.set_patch(example_patch)
        engine.render_patch(self.midi_note_pitch,
                            self.midi_note_velocity,
                            self.note_length_seconds,
                            self.render_length_seconds)
        example_audio = engine.get_audio_frames()

        # Get overridden parameters
        overridden_parameters = self.__get_overridden()
        overridden = [p[0] for p in overridden_parameters]

        # Assert number of data samples is specified
        assert sample_size is not None, "Please set the number of data samples required"

        # Initiate empty arrays to save data samples
        audio_set = np.zeros((sample_size, len(example_audio)), dtype=np.float32)
        patch_set = np.zeros((sample_size, len(self.random_parameters)), dtype=np.float32)

        for i in trange(sample_size):
            # Generate random patch
            random_patch = generator.get_random_patch()
            engine.set_patch(random_patch)

            # Overwrite fixed parameters and get full patch
            for param in overridden_parameters:
                engine.override_plugin_parameter(param[0], param[1])
            patch = engine.get_patch()

            # Render the patch
            engine.render_patch(self.midi_note_pitch,
                                self.midi_note_velocity,
                                self.note_length_seconds,
                                self.render_length_seconds)

            # Get decoded audio
            audio = engine.get_audio_frames()
            audio = np.array(audio, copy=True, dtype=np.float32)
            audio_set[i] = audio

            # Get synthesizer patch
            for element in sorted(overridden, reverse=True):
                patch.pop(element)
            patch_set[i] = [p[1] for p in patch]

            if self.save_audio:
                scipy.io.wavfile.write(
                    os.path.join(self.data_path, self.folder_name, "audio", f"{file_prefix}audio_{i}.wav"),
                    self.sample_rate,
                    audio)

        # Cut size audio array to exact length: sample_rate * render_length_seconds
        audio_set = audio_set[:, :int(self.sample_rate * self.render_length_seconds)]

        # Save dataset
        np.save(os.path.join(self.data_path, self.folder_name, f"{file_prefix}audio.npy"), audio_set)
        np.save(os.path.join(self.data_path, self.folder_name, f"{file_prefix}patches.npy"), patch_set)

    def __get_overridden(self):

        if self.base_preset == "MS-REV1_deepdiva.h2p":
            print("The default preset 'REV-MS1_deepdiva.h2p' is used for fixed parameters")

        assert os.path.exists(os.path.join(self.data_path, self.base_preset)), \
            f"The chosen preset {os.path.join(self.data_path, self.base_preset)} does not exist"

        used_base_preset = os.path.join(self.data_path, self.base_preset)
        h2p = H2P()
        base_patch = h2p.preset_to_patch(h2p_filename=used_base_preset)

        if self.random_parameters == []:
            print("None of the parameters is set to vary, all parameters will be fixed")

        overridden_parameters, _ = split_train_override_patch(base_patch, self.random_parameters)

        return overridden_parameters

ext: py

max_stars_repo_name: maudgrol/deepdiva

max_issues_repo_name: maudgrol/deepdiva

max_forks_repo_name: maudgrol/deepdiva

reason: import numpy,import scipy


text:
from collections import defaultdict
from logging import Logger
from typing import Dict, List
from re import VERBOSE

import icontract
import numpy as np
import pandas as pd
import prefect
from prefect import task

from secs.tasks.utilities import (
    dataframe_contains_invalid_references,
    replace_header_with_row,
    rename_columns_to_unique_names,
)

def _select_numeric_columns(df: pd.DataFrame, logger: Logger = None) -> List[str]:

    column_names_numeric = []
    for column_name in df.columns:

        column = df[column_name].copy()

        numeric_rows = """
        ^                   # beginning of string
        (?:[^A-Za-z]+ )?    # (optional) not preceded by a word
        (?:[*,])?           # (optional) preceded by * or *
        (\d+)               # capture the digits
        (?:[%])?            # (optional) followed by %
        (?: [^A-Za-z]+)?    # (optional) not followed by a word
        $                   # end of string
        """
        any_row_contains_a_valid_number = (
            column.astype(str).str.contains(numeric_rows, flags=VERBOSE).any()
        )
        if any_row_contains_a_valid_number:
            column_names_numeric.append(column_name)

    if logger:

        logger.debug(f"\n\nNumeric column names:\n{column_names_numeric}")
        column_names_non_numeric = np.setdiff1d(
            df.columns.to_list(), column_names_numeric
        )
        logger.debug(f"\nNon-numeric column names:\n{column_names_non_numeric}")

    return column_names_numeric

def _clean_numeric_columns(df: pd.DataFrame, logger: Logger = None) -> pd.DataFrame:

    column_names_numeric = _select_numeric_columns(df)
    df.loc[:, column_names_numeric] = (
        df[column_names_numeric]
        .copy()
        .replace(
            r"\s+", np.nan, regex=True
        )  # remove whitespace so can convert str to float
        .replace(
            to_replace=r"[^1-9.]", value="", regex=True
        )  # remove non-numeric characters
        .replace(r"", np.nan, regex=True)
        .fillna(0)
        .astype(np.number)  # convert string columns to numbers
        .convert_dtypes()  # infer ints
    )

    return df

def _drop_rows_where_first_column_empty(df: pd.DataFrame) -> pd.DataFrame:

    df = df.copy()

    first_column = df.columns[0]
    df = df.dropna(subset=[first_column])

    return df

@task
@icontract.ensure(lambda result: not result.empty, "Output cannot be empty!")
# @icontract.ensure(
#     lambda result: dataframe_contains_invalid_references(result),
#     "Output cannot contain invalid references!",
# )
def transform_sheet(
    excel_sheets_raw: List[pd.DataFrame], header_row: int,
) -> pd.DataFrame:

    logger = prefect.context.get("logger")

    excel_sheets_clean = [
        df.copy()
        .pipe(replace_header_with_row, header_row)
        .pipe(rename_columns_to_unique_names)
        .replace("?", np.nan)
        .replace(0, np.nan)
        .pipe(_clean_numeric_columns, logger)
        .pipe(_drop_rows_where_first_column_empty)
        for df in excel_sheets_raw
    ]

    df = pd.concat(excel_sheets_clean).reset_index(drop=True)

    return df

max_stars_repo_name: rdmolony/merge-sec-mentor-excels

max_issues_repo_name: rdmolony/merge-sec-mentor-excels

max_forks_repo_name: rdmolony/merge-sec-mentor-excels

reason: import numpy


text:
"""Collection of routines for RFI excision."""
import numpy as np
from scipy.ndimage import convolve1d

def sumthreshold_py(
    data,
    max_m=16,
    start_flag=None,
    threshold1=None,
    remove_median=True,
    correct_for_missing=True,
):
    """SumThreshold outlier detection algorithm.

    See http://www.astro.rug.nl/~offringa/SumThreshold.pdf for description of
    the algorithm.

    Parameters
    ----------
    data : np.ndarray[:, :]
        The data to flag.
    max_m : int, optional
        Maximum size to expand to.
    start_flag : np.ndarray[:, :], optional
        A boolean array of the initially flagged data.
    threshold1 : float, optional
        Initial threshold. By default use the 95 percentile.
    remove_median : bool, optional
        Subtract the median of the full 2D dataset. Default is True.
    correct_for_missing : bool, optional

    Returns
    -------
    mask : np.ndarray[:, :]
        Boolean array, with `True` entries marking outlier data.
    """
    data = np.copy(data)
    (ny, nx) = data.shape

    if start_flag is None:
        start_flag = np.isnan(data)
    flag = np.copy(start_flag)

    if remove_median:
        data -= np.median(data[~flag])

    if threshold1 is None:
        threshold1 = np.percentile(data[~flag], 95.0)

    m = 1
    while m <= max_m:

        if m == 1:
            threshold = threshold1
        else:
            threshold = threshold1 / 1.5 ** (np.log2(m))

        # The centre of the window for even windows is the bin right to the left of
        # centre. I want the origin at the leftmost bin
        if m == 1:
            centre = 0
        else:
            centre = m // 2 - 1

        ## X-axis

        data[flag] = 0.0
        count = (~flag).astype(np.float)

        # Convolution of the data
        dconv = convolve1d(
            data, weights=np.ones(m, dtype=float), origin=-centre, axis=1
        )[:, : (nx - m + 1)]

        # Convolution of the counts
        cconv = convolve1d(
            count, weights=np.ones(m, dtype=float), origin=-centre, axis=1
        )[:, : (nx - m + 1)]
        if correct_for_missing:
            cconv = m ** 0.5 * cconv ** 0.5
        flag_temp = dconv > cconv * threshold
        flag_temp += dconv < -cconv * threshold
        for ii in range(flag_temp.shape[1]):
            flag[:, ii : (ii + m)] += flag_temp[:, ii][:, np.newaxis]

        ## Y-axis

        data[flag] = 0.0
        count = (~flag).astype(np.float)
        # Convolution of the data
        dconv = convolve1d(
            data, weights=np.ones(m, dtype=float), origin=-centre, axis=0
        )[: (ny - m + 1), :]
        # Convolution of the counts
        cconv = convolve1d(
            count, weights=np.ones(m, dtype=float), origin=-centre, axis=0
        )[: (ny - m + 1), :]
        if correct_for_missing:
            cconv = m ** 0.5 * cconv ** 0.5
        flag_temp = dconv > cconv * threshold
        flag_temp += dconv < -cconv * threshold

        for ii in range(flag_temp.shape[0]):
            flag[ii : ii + m, :] += flag_temp[ii, :][np.newaxis, :]

        m *= 2

    return flag

# This routine might be substituted by a faster one later
sumthreshold = sumthreshold_py

# Scale-invariant rank (SIR) functions
def sir1d(basemask, eta=0.2):
    """Numpy implementation of the scale-invariant rank (SIR) operator.

    For more information, see arXiv:1201.3364v2.

    Parameters
    ----------
    basemask : numpy 1D array of boolean type
        Array with the threshold mask previously generated.
        1 (True) for flagged points, 0 (False) otherwise.
    eta : float
        Aggressiveness of the method: with eta=0, no additional samples are
        flagged and the function returns basemask. With eta=1, all samples
        will be flagged. The authors in arXiv:1201.3364v2 seem to be convinced
        that 0.2 is a mostly universally optimal value, but no optimization
        has been done on CHIME data.

    Returns
    -------
    mask : numpy 1D array of boolean type
        The mask after the application of the (SIR) operator. Same shape and
        type as basemask.
    """
    n = basemask.size
    psi = basemask.astype(np.float) - 1.0 + eta

    M = np.zeros(n + 1, dtype=np.float)
    M[1:] = np.cumsum(psi)

    MP = np.minimum.accumulate(M)[:-1]
    MQ = np.concatenate((np.maximum.accumulate(M[-2::-1])[-2::-1], M[-1, np.newaxis]))

    return (MQ - MP) >= 0.0

def sir(basemask, eta=0.2, only_freq=False, only_time=False):
    """Apply the SIR operator over the frequency and time axes for each product.

    This is a wrapper for `sir1d`.  It loops over times, applying `sir1d`
    across the frequency axis.  It then loops over frequencies, applying `sir1d`
    across the time axis.  It returns the logical OR of these two masks.

    Parameters
    ----------
    basemask : np.ndarray[nfreq, nprod, ntime] of boolean type
        The previously generated threshold mask.
        1 (True) for masked points, 0 (False) otherwise.
    eta : float
        Aggressiveness of the method: with eta=0, no additional samples are
        flagged and the function returns basemask. With eta=1, all samples
        will be flagged.
    only_freq : bool
        Only apply the SIR operator across the frequency axis.
    only_time : bool
        Only apply the SIR operator across the time axis.

    Returns
    -------
    mask : np.ndarray[nfreq, nprod, ntime] of boolean type
        The mask after the application of the SIR operator.
    """
    if only_freq and only_time:
        raise ValueError("Only one of only_freq and only_time can be True.")

    nfreq, nprod, ntime = basemask.shape

    newmask = basemask.astype(np.bool).copy()

    for pp in range(nprod):

        if not only_time:
            for tt in range(ntime):
                newmask[:, pp, tt] |= sir1d(basemask[:, pp, tt], eta=eta)

        if not only_freq:
            for ff in range(nfreq):
                newmask[ff, pp, :] |= sir1d(basemask[ff, pp, :], eta=eta)

    return newmask

max_stars_repo_name: Arnab-half-blood-prince/draco

max_issues_repo_name: Arnab-half-blood-prince/draco

max_forks_repo_name: Arnab-half-blood-prince/draco


text:
import torch
import numpy as np
import torch.nn as nn
import torch.nn.functional as F

import matplotlib.pyplot as plt
from fourier_neural_operator.utilities3 import *

import operator
from functools import reduce
from functools import partial

from timeit import default_timer
import scipy.io

torch.manual_seed(0)
np.random.seed(0)

from fourier_neural_operator.lowrank_2d_time import *

activation = F.relu

if __name__ == "__main__":
    ################################################################
    # configs
    ################################################################
    # TRAIN_PATH = 'data/ns_data_V10000_N1200_T20.mat'
    # TEST_PATH = 'data/ns_data_V10000_N1200_T20.mat'
    # TRAIN_PATH = 'data/ns_data_V1000_N1000_train.mat'
    # TEST_PATH = 'data/ns_data_V1000_N1000_train_2.mat'
    # TRAIN_PATH = 'data/ns_data_V1000_N5000.mat'
    # TEST_PATH = 'data/ns_data_V1000_N5000.mat'
    TRAIN_PATH = 'data/ns_data_V100_N1000_T50_1.mat'
    TEST_PATH = 'data/ns_data_V100_N1000_T50_2.mat'

    ntrain = 1000
    ntest = 200

    batch_size = 5
    batch_size2 = batch_size

    epochs = 500
    learning_rate = 0.0025
    scheduler_step = 100
    scheduler_gamma = 0.5

    print(epochs, learning_rate, scheduler_step, scheduler_gamma)

    path = 'ns_lowrank_rnn_V100_T40_N'+str(ntrain)+'_ep' + str(epochs) + '_m'
    path_model = 'model/'+path
    path_train_err = 'results/'+path+'train.txt'
    path_test_err = 'results/'+path+'test.txt'
    path_image = 'image/'+path

    runtime = np.zeros(2, )
    t1 = default_timer()

    sub = 1
    S = 64
    T_in = 10
    T = 40
    step = 1

    ################################################################
    # load dataset
    ################################################################
    reader = MatReader(TRAIN_PATH)
    train_a = reader.read_field('u')[:ntrain,::sub,::sub,:T_in]
    train_u = reader.read_field('u')[:ntrain,::sub,::sub,T_in:T+T_in]

    reader = MatReader(TEST_PATH)
    test_a = reader.read_field('u')[-ntest:,::sub,::sub,:T_in]
    test_u = reader.read_field('u')[-ntest:,::sub,::sub,T_in:T+T_in]

    print(train_u.shape)
    print(test_u.shape)
    assert (S == train_u.shape[-2])
    assert (T == train_u.shape[-1])

    train_a = train_a.reshape(ntrain,S,S,T_in)
    test_a = test_a.reshape(ntest,S,S,T_in)

    gridx = torch.tensor(np.linspace(0, 1, S), dtype=torch.float)
    gridx = gridx.reshape(1, S, 1, 1).repeat([1, 1, S, 1])
    gridy = torch.tensor(np.linspace(0, 1, S), dtype=torch.float)
    gridy = gridy.reshape(1, 1, S, 1).repeat([1, S, 1, 1])

    train_a = torch.cat((train_a, gridx.repeat([ntrain,1,1,1]), gridy.repeat([ntrain,1,1,1])), dim=-1)
    test_a = torch.cat((test_a, gridx.repeat([ntest,1,1,1]), gridy.repeat([ntest,1,1,1])), dim=-1)

    train_loader = torch.utils.data.DataLoader(torch.utils.data.TensorDataset(train_a, train_u), batch_size=batch_size, shuffle=True)
    test_loader = torch.utils.data.DataLoader(torch.utils.data.TensorDataset(test_a, test_u), batch_size=batch_size, shuffle=False)

    t2 = default_timer()

    print('preprocessing finished, time used:', t2-t1)
    device = torch.device('cuda')

    ################################################################
    # training and evaluation
    ################################################################
    model = Net2d().cuda()
    # model = torch.load('model/ns_fourier_V100_N1000_ep100_m8_w20')

    print(model.count_params())
    optimizer = torch.optim.Adam(model.parameters(), lr=learning_rate, weight_decay=1e-4)
    scheduler = torch.optim.lr_scheduler.StepLR(optimizer, step_size=scheduler_step, gamma=scheduler_gamma)

    myloss = LpLoss(size_average=False)
    gridx = gridx.to(device)
    gridy = gridy.to(device)

    for ep in range(epochs):
        model.train()
        t1 = default_timer()
        train_l2_step = 0
        train_l2_full = 0
        for xx, yy in train_loader:
            loss = 0
            xx = xx.to(device)
            yy = yy.to(device)

            for t in range(0, T, step):
                y = yy[..., t:t + step]
                im = model(xx)
                loss += myloss(im.reshape(batch_size, -1), y.reshape(batch_size, -1))

                if t == 0:
                    pred = im
                else:
                    pred = torch.cat((pred, im), -1)

                xx = torch.cat((xx[..., step:-2], im,
                                gridx.repeat([batch_size, 1, 1, 1]), gridy.repeat([batch_size, 1, 1, 1])), dim=-1)

            train_l2_step += loss.item()
            l2_full = myloss(pred.reshape(batch_size, -1), yy.reshape(batch_size, -1))
            train_l2_full += l2_full.item()

            optimizer.zero_grad()
            loss.backward()
            # l2_full.backward()
            optimizer.step()

        test_l2_step = 0
        test_l2_full = 0
        with torch.no_grad():
            for xx, yy in test_loader:
                loss = 0
                xx = xx.to(device)
                yy = yy.to(device)

                for t in range(0, T, step):
                    y = yy[..., t:t + step]
                    im = model(xx)
                    loss += myloss(im.reshape(batch_size, -1), y.reshape(batch_size, -1))

                    if t == 0:
                        pred = im
                    else:
                        pred = torch.cat((pred, im), -1)

                    xx = torch.cat((xx[..., step:-2], im,
                                    gridx.repeat([batch_size, 1, 1, 1]), gridy.repeat([batch_size, 1, 1, 1])), dim=-1)

                test_l2_step += loss.item()
                test_l2_full += myloss(pred.reshape(batch_size, -1), yy.reshape(batch_size, -1)).item()

        t2 = default_timer()
        scheduler.step()
        print(ep, t2 - t1, train_l2_step / ntrain / (T / step), train_l2_full / ntrain, test_l2_step / ntest / (T / step),
              test_l2_full / ntest)
    # torch.save(model, path_model)

    # pred = torch.zeros(test_u.shape)
    # index = 0
    # test_loader = torch.utils.data.DataLoader(torch.utils.data.TensorDataset(test_a, test_u), batch_size=1, shuffle=False)
    # with torch.no_grad():
    #     for x, y in test_loader:
    #         test_l2 = 0;
    #         x, y = x.cuda(), y.cuda()
    #
    #         out = model(x)
    #         out = y_normalizer.decode(out)
    #         pred[index] = out
    #
    #         test_l2 += myloss(out.view(1, -1), y.view(1, -1)).item()
    #         print(index, test_l2)
    #         index = index + 1

    # scipy.io.savemat('pred/'+path+'.mat', mdict={'pred': pred.cpu().numpy()})


text:
import numpy as np
import pandas as pd
from sklearn.decomposition import PCA
from matplotlib import pyplot as plt

print("Download datasets ....")
digits_train = pd.read_csv(
    'https://archive.ics.uci.edu/ml/machine-learning-databases/optdigits/optdigits.tra',
    header=None)
digits_test = pd.read_csv(
    'https://archive.ics.uci.edu/ml/machine-learning-databases/optdigits/optdigits.tes',
    header=None)
print("Datasets downloaded!")

X_digits = digits_train[np.arange(64)]
y_digits = digits_train[64]

estimator = PCA(n_components=2)
X_pca = estimator.fit_transform(X_digits)

def plot_pca_scatter():
    colors = ['black', 'blue', 'purple', 'yellow', 'white', 'red', 'lime', 'cyan', 'orange', 'gray']
    for i in range(len(colors)):
        px = X_pca[:, 0][y_digits.to_numpy() == i]
        py = X_pca[:, 1][y_digits.to_numpy() == i]
        plt.scatter(px, py, c=colors[i])
    plt.legend(np.arange(0, 10).astype(str))
    plt.xlabel('First Principal Component!')
    plt.ylabel('Second Principal Component!')
    plt.show()

plot_pca_scatter()


text:
import torch
import torch.nn as nn
import torch.optim as optim
import torch.nn.functional as F
from torchvision.datasets.mnist import MNIST
import torchvision.transforms as transforms
from torch.utils.data import DataLoader
from torch.autograd import Variable

import numpy as np
import argparse
from matplotlib import pyplot as plt
from PIL import Image

from helpers import get_device, rotate_img, one_hot_embedding
from data import dataloaders, digit_one
from train import train_model
from test import rotating_image_classification, test_single_image
from losses import edl_mse_loss, edl_digamma_loss, edl_log_loss, relu_evidence
from lenet import LeNet

def main():

    parser = argparse.ArgumentParser()
    mode_group = parser.add_mutually_exclusive_group(required=True)
    mode_group.add_argument(
        "--train", action="store_true", help="To train the network."
    )
    mode_group.add_argument("--test", action="store_true", help="To test the network.")
    mode_group.add_argument(
        "--examples", action="store_true", help="To example MNIST data."
    )
    parser.add_argument(
        "--epochs", default=10, type=int, help="Desired number of epochs."
    )
    parser.add_argument(
        "--dropout", action="store_true", help="Whether to use dropout or not."
    )
    parser.add_argument(
        "--uncertainty", action="store_true", help="Use uncertainty or not."
    )
    uncertainty_type_group = parser.add_mutually_exclusive_group()
    uncertainty_type_group.add_argument(
        "--mse",
        action="store_true",
        help="Set this argument when using uncertainty. Sets loss function to Expected Mean Square Error.",
    )
    uncertainty_type_group.add_argument(
        "--digamma",
        action="store_true",
        help="Set this argument when using uncertainty. Sets loss function to Expected Cross Entropy.",
    )
    uncertainty_type_group.add_argument(
        "--log",
        action="store_true",
        help="Set this argument when using uncertainty. Sets loss function to Negative Log of the Expected Likelihood.",
    )
    args = parser.parse_args()

    if args.examples:
        examples = enumerate(dataloaders["val"])
        batch_idx, (example_data, example_targets) = next(examples)
        fig = plt.figure()
        for i in range(6):
            plt.subplot(2, 3, i + 1)
            plt.tight_layout()
            plt.imshow(example_data[i][0], cmap="gray", interpolation="none")
            plt.title("Ground Truth: {}".format(example_targets[i]))
            plt.xticks([])
            plt.yticks([])
        plt.savefig("./images/examples.jpg")

    elif args.train:
        num_epochs = args.epochs
        use_uncertainty = args.uncertainty
        num_classes = 10

        model = LeNet(dropout=args.dropout)

        if use_uncertainty:
            if args.digamma:
                criterion = edl_digamma_loss
            elif args.log:
                criterion = edl_log_loss
            elif args.mse:
                criterion = edl_mse_loss
            else:
                parser.error("--uncertainty requires --mse, --log or --digamma.")
        else:
            criterion = nn.CrossEntropyLoss()

        optimizer = optim.Adam(model.parameters(), lr=1e-3, weight_decay=0.005)

        exp_lr_scheduler = optim.lr_scheduler.StepLR(optimizer, step_size=7, gamma=0.1)

        device = get_device()
        model = model.to(device)

        model, metrics = train_model(
            model,
            dataloaders,
            num_classes,
            criterion,
            optimizer,
            scheduler=exp_lr_scheduler,
            num_epochs=num_epochs,
            device=device,
            uncertainty=use_uncertainty,
        )

        state = {
            "epoch": num_epochs,
            "model_state_dict": model.state_dict(),
            "optimizer_state_dict": optimizer.state_dict(),
        }

        if use_uncertainty:
            if args.digamma:
                torch.save(state, "./results/model_uncertainty_digamma.pt")
                print("Saved: ./results/model_uncertainty_digamma.pt")
            if args.log:
                torch.save(state, "./results/model_uncertainty_log.pt")
                print("Saved: ./results/model_uncertainty_log.pt")
            if args.mse:
                torch.save(state, "./results/model_uncertainty_mse.pt")
                print("Saved: ./results/model_uncertainty_mse.pt")

        else:
            torch.save(state, "./results/model.pt")
            print("Saved: ./results/model.pt")

    elif args.test:

        use_uncertainty = args.uncertainty
        device = get_device()
        model = LeNet()
        model = model.to(device)
        optimizer = optim.Adam(model.parameters())

        if use_uncertainty:
            if args.digamma:
                checkpoint = torch.load("./results/model_uncertainty_digamma.pt")
                filename = "./results/rotate_uncertainty_digamma.jpg"
            if args.log:
                checkpoint = torch.load("./results/model_uncertainty_log.pt")
                filename = "./results/rotate_uncertainty_log.jpg"
            if args.mse:
                checkpoint = torch.load("./results/model_uncertainty_mse.pt")
                filename = "./results/rotate_uncertainty_mse.jpg"

        else:
            checkpoint = torch.load("./results/model.pt")
            filename = "./results/rotate.jpg"

        model.load_state_dict(checkpoint["model_state_dict"])
        optimizer.load_state_dict(checkpoint["optimizer_state_dict"])

        model.eval()

        rotating_image_classification(
            model, digit_one, filename, uncertainty=use_uncertainty
        )

        test_single_image(model, "./data/one.jpg", uncertainty=use_uncertainty)
        test_single_image(model, "./data/yoda.jpg", uncertainty=use_uncertainty)

if __name__ == "__main__":
    main()


text:
# -*- coding: utf-8 -*-
"""
/***************************************************************************
 ZomatoNightlifeMapper
                                 A QGIS plugin
 The plugin generates nightlife popularity grid for the given area of interest using Zomato's Nightlife Index.
                              -------------------
        begin                : 2016-06-16
        git sha              : $Format:%H$
        copyright            : (C) 2016 by Serhat YILMAZTURK
        email                : serhat@yilmazturk.info
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/
"""
from PyQt4.QtCore import QSettings, QTranslator, qVersion, QCoreApplication, QUrl, pyqtSlot, QVariant
from PyQt4.QtGui import QAction, QIcon, QMessageBox, QFileDialog
from PyQt4 import uic
# Initialize Qt resources from file resources.py
import resources
# Import the code for the dialog
from zomato_mapper_dialog import ZomatoNightlifeMapperDialog
import os.path
from qgis.utils import iface
from qgis.gui import QgsMessageBar
from qgis.core import QgsVectorLayer, QgsGeometry, QgsFeature, QgsPoint, QgsMapLayerRegistry, QgsField, \
    QgsVectorFileWriter
from numpy import arange
import itertools
import urllib2
import json
import time
import datetime

js_ne_lat = "$(ne_lat)"
js_ne_lng = "$(ne_lng)"
js_sw_lat = "$(sw_lat)"
js_sw_lng = "$(sw_lng)"

class ZomatoNightlifeMapper:
    """QGIS Plugin Implementation."""
    def __init__(self, iface):
        """Constructor.

        :param iface: An interface instance that will be passed to this class
            which provides the hook by which you can manipulate the QGIS
            application at run time.
        :type iface: QgsInterface
        """
        # Save reference to the QGIS interface
        self.iface = iface
        # initialize plugin directory
        self.plugin_dir = os.path.dirname(__file__)
        # initialize locale
        locale = QSettings().value('locale/userLocale')[0:2]
        locale_path = os.path.join(
            self.plugin_dir,
            'i18n',
            'ZomatoNightlifeMapper_{}.qm'.format(locale))

        if os.path.exists(locale_path):
            self.translator = QTranslator()
            self.translator.load(locale_path)

            if qVersion() > '4.3.3':
                QCoreApplication.installTranslator(self.translator)

        # Create the dialog (after translation) and keep reference
        self.dlg = ZomatoNightlifeMapperDialog()
        self.dlg.setFixedSize(382, 232)

        # Declare instance attributes
        self.actions = []
        self.menu = self.tr(u'&Zomato Nightlife Mapper')
        # TODO: We are going to let the user set this up in a future iteration
        self.toolbar = self.iface.addToolBar(u'ZomatoNightlifeMapper')
        self.toolbar.setObjectName(u'ZomatoNightlifeMapper')

        self.dlg.pushButton_aoi.clicked.connect(self.popup_show)
        self.dlg.pushButton_grid.clicked.connect(self.show_grid)

    # noinspection PyMethodMayBeStatic
    def tr(self, message):
        """Get the translation for a string using Qt translation API.

        We implement this ourselves since we do not inherit QObject.

        :param message: String for translation.
        :type message: str, QString

        :returns: Translated version of message.
        :rtype: QString
        """
        # noinspection PyTypeChecker,PyArgumentList,PyCallByClass
        return QCoreApplication.translate('ZomatoNightlifeMapper', message)

    def add_action(
        self,
        icon_path,
        text,
        callback,
        enabled_flag=True,
        add_to_menu=True,
        add_to_toolbar=True,
        status_tip=None,
        whats_this=None,
        parent=None):
        """Add a toolbar icon to the toolbar.

        :param icon_path: Path to the icon for this action. Can be a resource
            path (e.g. ':/plugins/foo/bar.png') or a normal file system path.
        :type icon_path: str

        :param text: Text that should be shown in menu items for this action.
        :type text: str

        :param callback: Function to be called when the action is triggered.
        :type callback: function

        :param enabled_flag: A flag indicating if the action should be enabled
            by default. Defaults to True.
        :type enabled_flag: bool

        :param add_to_menu: Flag indicating whether the action should also
            be added to the menu. Defaults to True.
        :type add_to_menu: bool

        :param add_to_toolbar: Flag indicating whether the action should also
            be added to the toolbar. Defaults to True.
        :type add_to_toolbar: bool

        :param status_tip: Optional text to show in a popup when mouse pointer
            hovers over the action.
        :type status_tip: str

        :param parent: Parent widget for the new action. Defaults None.
        :type parent: QWidget

        :param whats_this: Optional text to show in the status bar when the
            mouse pointer hovers over the action.

        :returns: The action that was created. Note that the action is also
            added to self.actions list.
        :rtype: QAction
        """

        icon = QIcon(icon_path)
        action = QAction(icon, text, parent)
        action.triggered.connect(callback)
        action.setEnabled(enabled_flag)

        if status_tip is not None:
            action.setStatusTip(status_tip)

        if whats_this is not None:
            action.setWhatsThis(whats_this)

        if add_to_toolbar:
            self.toolbar.addAction(action)

        if add_to_menu:
            self.iface.addPluginToMenu(
                self.menu,
                action)

        self.actions.append(action)

        return action

    def initGui(self):
        """Create the menu entries and toolbar icons inside the QGIS GUI."""

        icon_path = ':/plugins/ZomatoNightlifeMapper/icon.png'
        self.add_action(
            icon_path,
            text=self.tr(u'Zomato Nightlife Mapper'),
            callback=self.run,
            parent=self.iface.mainWindow())

    def unload(self):
        """Removes the plugin menu item and icon from QGIS GUI."""
        for action in self.actions:
            self.iface.removePluginMenu(
                self.tr(u'&Zomato Nightlife Mapper'),
                action)
            self.iface.removeToolBarIcon(action)
        # remove the toolbar
        del self.toolbar

    def popup_show(self):
        @pyqtSlot()
        def close_popup():
            global ne_lat, ne_lng, sw_lat, sw_lng
            ne_lat = self.ui.webView_aoi.page().mainFrame().evaluateJavaScript(js_ne_lat)
            ne_lng = self.ui.webView_aoi.page().mainFrame().evaluateJavaScript(js_ne_lng)
            sw_lat = self.ui.webView_aoi.page().mainFrame().evaluateJavaScript(js_sw_lat)
            sw_lng = self.ui.webView_aoi.page().mainFrame().evaluateJavaScript(js_sw_lng)
            self.ui.close()

        self.ui = uic.loadUi(os.path.join(os.path.dirname(__file__), 'webView_gmap.ui'))
        self.ui.webView_aoi.load(QUrl('http://yilmazturk.info/zomato/index.html'))
        self.ui.setWindowTitle('Define area of interest')
        self.ui.pushButton.clicked.connect(close_popup)
        self.ui.show()

    def show_grid(self):
        try:
            point_interval = self.dlg.lineEdit_interval.text()
            if ',' in point_interval:
                iface.messageBar().pushMessage(u"Warning: ", "Please use point (.) as separator",
                                               level=QgsMessageBar.WARNING, duration=3)
            else:
                start_time = time.time()

                latlist = list(arange(sw_lat['0'], ne_lat['0'], float(point_interval)))
                lnglist = list(arange(sw_lng['0'], ne_lng['0'], float(point_interval)))

                msgBox = QMessageBox()
                msgBox.setWindowTitle('Info')
                message = "%s points will be generated. Continue?" % str(len(list(itertools.product(latlist, lnglist))))
                msgBox.setText(message)
                msgBox.setInformativeText("Zomato API daily limit is up to 1000 calls.")
                msgBox.setStandardButtons(QMessageBox.Yes | QMessageBox.No)
                pop = msgBox.exec_()

                if pop == QMessageBox.Yes:
                    directory = QFileDialog.getExistingDirectory(None, "Select folder to save shapefile...")
                    shp_dir = directory.replace("\\", "/")

                    if len(directory) > 1:
                        current_date = str(datetime.datetime.now().date()).replace("-", "")
                        layer_name = 'grid_' + current_date
                        coord_pairs = []
                        memory_layer = QgsVectorLayer("Point?crs=epsg:4326", layer_name, "memory")
                        memory_layer.startEditing()
                        provider = memory_layer.dataProvider()
                        provider.addAttributes([QgsField("NLINDEX", QVariant.Double)])

                        for i in itertools.product(latlist, lnglist):
                            api_key = self.dlg.lineEdit_apikey.text()
                            url = "https://developers.zomato.com/api/v2.1/geocode?lat=%f&lon=%f" % (i[0], i[1])
                            headers = {"User-agent": "curl/7.43.0", "user_key": api_key, "Accept": "application/json"}
                            req = urllib2.Request(url, None, headers)
                            response = urllib2.urlopen(req)
                            data = json.load(response)
                            nightlife_index = data['popularity']['nightlife_index']

                            geometry = QgsGeometry.fromPoint(QgsPoint(i[1], i[0]))
                            feature = QgsFeature()
                            feature.setGeometry(geometry)
                            feature.setAttributes([nightlife_index])
                            coord_pairs.append(feature)

                        memory_layer.dataProvider().addFeatures(coord_pairs)
                        memory_layer.updateExtents()
                        memory_layer.commitChanges()
                        QgsMapLayerRegistry.instance().addMapLayer(memory_layer)

                        shp_path = shp_dir + "/" + layer_name
                        QgsVectorFileWriter.writeAsVectorFormat(memory_layer, shp_path,
                                                                "utf-8", None, "ESRI Shapefile")

                        end_time = time.time()
                        time_difference = end_time - start_time
                        time_difference_rounded = "%.2f" % time_difference

                        iface.messageBar().pushMessage(u"Info: ",
                                                       str(time_difference_rounded) + " seconds have elapsed for " + str(
                                                           len(list(itertools.product(latlist, lnglist)))) + " points.",
                                                       level=QgsMessageBar.SUCCESS, duration=5)
                    else:
                        iface.messageBar().pushMessage(u"Warning: ", u"Folder is not selected...",
                                                       level=QgsMessageBar.WARNING, duration=5)
                else:
                    pass
        except:
            iface.messageBar().pushMessage(u"Error: ", u"Check your API key and/or point interval value!",
                                           level=QgsMessageBar.CRITICAL, duration=5)

    def run(self):
        """Run method that performs all the real work"""
        # show the dialog
        self.dlg.show()
        # Run the dialog event loop
        result = self.dlg.exec_()
        # See if OK was pressed
        if result:
            # Do something useful here - delete the line containing pass and
            # substitute with your code.
            pass

max_stars_repo_name: syilmazturk/ZomatoNightlifeMapper

max_issues_repo_name: syilmazturk/ZomatoNightlifeMapper

max_forks_repo_name: syilmazturk/ZomatoNightlifeMapper

reason: from numpy


text:
import glob
import os
import operator
import sys

import numpy as np
import matplotlib.pyplot as plt

import yt

yt.funcs.mylog.setLevel(50)

class Profile:
    """read a plotfile using yt and store the 1d profile for T and enuc"""

    def __init__(self, plotfile):
        ds = yt.load(plotfile)

        time = float(ds.current_time)
        ad = ds.all_data()

        # Sort the ray values by 'x' so there are no discontinuities
        # in the line plot
        srt = np.argsort(ad['x'])
        x_coord = np.array(ad['x'][srt])
        temp = np.array(ad['Temp'][srt])
        enuc = np.array(ad['enuc'][srt])

        self.time = time
        self.x = x_coord
        self.T = temp
        self.enuc = enuc

    def find_x_for_T(self, T_0=1.e9):
        """ given a profile x(T), find the x_0 that corresponds to T_0 """

        # our strategy here assumes that the hot ash is in the early
        # part of the profile.  We then find the index of the first
        # point where T drops below T_0
        idx = np.where(self.T < T_0)[0][0]

        T1 = self.T[idx-1]
        x1 = self.x[idx-1]

        T2 = self.T[idx]
        x2 = self.x[idx]

        slope = (x2 - x1)/(T2 - T1)

        return x1 + slope*(T_0 - T1)

def get_velocity(p1, p2):
        """look at the last 2 plotfiles and estimate the velocity by
        finite-differencing"""

        # we'll do this by looking at 3 different temperature
        # thresholds and averaging
        T_ref = [2.e9, 3.e9, 4.e9]

        vs = []
        for T0 in T_ref:
            x1 = p1.find_x_for_T(T0)
            x2 = p2.find_x_for_T(T0)
            vs.append((x1 - x2)/(p1.time - p2.time))

        vs = np.array(vs)
        v = np.mean(vs)
        v_sigma = np.std(vs)
        return v, v_sigma

if __name__ == "__main__":

    run_dir = sys.argv[1]

    print(run_dir)

    # get all the plotfiles
    pfiles = glob.glob("{}/*plt?????".format(run_dir))
    pfiles.sort()

    profiles = []
    for pf in pfiles:
        profiles.append(Profile(pf))

    for n in range(1, len(profiles)):
        p1 = profiles[n-1]
        p2 = profiles[n]
        print(0.5*(p1.time + p2.time), get_velocity(p1, p2))

max_stars_repo_name: MargotF/Castro

max_issues_repo_name: MargotF/Castro

max_forks_repo_name: MargotF/Castro


text:
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#  Copyright 2021-  QuOCS Team
#
#  Licensed under the Apache License, Version 2.0 (the "License");
#  you may not use this file except in compliance with the License.
#  You may obtain a copy of the License at
#
#      http://www.apache.org/licenses/LICENSE-2.0
#
#  Unless required by applicable law or agreed to in writing, software
#  distributed under the License is distributed on an "AS IS" BASIS,
#  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#  See the License for the specific language governing permissions and
#  limitations under the License.
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

import numpy as np

from quocslib.pulses.BasePulse import BasePulse
from quocslib.pulses.basis.ChoppedBasis import ChoppedBasis
from quocslib.tools.randomgenerator import RandomNumberGenerator

class Walsh(ChoppedBasis):
    # TODO Complete this basis
    amplitude_variation: float
    optimized_control_parameters: np.ndarray
    optimized_super_parameters: np.ndarray
    time_grid: np.ndarray

    def __init__(self, map_index: int, pulse_dictionary: dict, rng: RandomNumberGenerator = None):
        """

        :param int map_index: Index number to use to get the control parameters for the Fourier basis
        :param dict pulse_dictionary: The dictionary of the pulse defined here. Only the basis dictionary is used btw
        """
        #################
        # Basis dependent settings
        #################
        basis_dict = pulse_dictionary["basis"]
        # Super Parameter number i.e. the basis vector number in the pulse parametrization
        self.super_parameter_number = basis_dict.setdefault("basis_vector_number", 1)
        # Number of control parameters to be optimized
        self.control_parameters_number = 2 * self.super_parameter_number
        #################
        # Standard Basis Settings: amplitude limits, amplitude variation for the simplex,
        # distribution of super parameters, etc ...
        ################
        # Constructor of the parent classes, i.e. Base Pulse and Chopped Basis
        super().__init__(map_index=map_index, **pulse_dictionary)
        #################
        # Basis dependent settings
        #################
        # Scale coefficients: average distance of the points in the intial simplex
        self.scale_coefficients = (
            self.amplitude_variation
            / np.sqrt(2)
            * np.ones((self.control_parameters_number,))
        )
        # Initial value of the parameters in the pulse parametrization
        self.offset_coefficients = np.zeros((self.control_parameters_number,))

    def _get_shaped_pulse(self) -> np.array:
        """Definition of the pulse parametrization. It is called at every function evaluation to build the pulse"""
        #################
        # Standard Basis Settings: amplitude limits, amplitude variation for the simplex,
        # distribution of super parameters, etc ...
        ################
        # Pulse initialization
        pulse = np.zeros(self.bins_number)
        # Final time definition
        final_time = self.final_time
        # Pulse creation
        xx = self.optimized_control_parameters
        w = self.super_parameter_distribution_obj.w
        t = self.time_grid
        #################
        # Basis dependent settings
        #################
        for ii in range(self.super_parameter_number):
            pulse += xx[2 * ii] * np.sin(2 * np.pi * w[ii] * t / final_time) + xx[
                2 * ii + 1
            ] * np.cos(2 * np.pi * w[ii] * t / final_time)
        return self._limit_constraints(pulse)

    def _limit_constraints(self, pulse: np.ndarray) -> np.ndarray:
        """Limit the dCRAB pulse accordingly to the system constraints"""
        # Add the initial guess
        optimal_pulse = self.add_initial_guess(pulse)
        # Constraint the pulse
        optimal_pulse = self.limit_pulse(optimal_pulse)
        # Substract the initial guess pulse
        optimal_pulse = optimal_pulse - self._get_initial_guess()
        return optimal_pulse

max_stars_repo_name: marcorossignolo/QuOCS

max_issues_repo_name: marcorossignolo/QuOCS

max_forks_repo_name: marcorossignolo/QuOCS


text:
import sys
import requests
import numpy as np

from flask import Flask, request, jsonify
from PIL import Image, ImageOps
sys.path.append('src')

from utils.utils import calc_mean_score
from evaluater.predict import main
from handlers.model_builder import Nima
from handlers.data_generator import TestDataGenerator

base_model_name = 'MobileNet'

weights_file_aesthetic = 'models/MobileNet/weights_mobilenet_aesthetic_0.07.hdf5'
weights_file_technical = 'models/MobileNet/weights_mobilenet_technical_0.11.hdf5'

def build_model(base_model_name, weights_file):
    nima = Nima(base_model_name, weights=None)
    nima.build()
    nima.nima_model.load_weights(weights_file)
    return nima

model_aesthetic = build_model(base_model_name, weights_file_aesthetic)
model_technical = build_model(base_model_name, weights_file_technical)

app = Flask(__name__)

def resize_image(img):
    target_size = 224
    img.thumbnail((target_size, target_size))

    delta_w = target_size - img.width
    delta_h = target_size - img.height
    padding = (delta_w//2, delta_h//2, delta_w-(delta_w//2), delta_h-(delta_h//2))
    return ImageOps.expand(img, padding)

@app.route('/', methods=['POST'])
def predict():
    req = request.json
    result = []
    for image in req['images']:
        stream = requests.get(image, stream=True).raw
        img = Image.open(stream)

        preprocess_fn = model_aesthetic.preprocessing_function()
        img = resize_image(img)
        img = preprocess_fn(np.array(img))

        img = img[np.newaxis,:,:,:]

        aesthetic_pred = model_aesthetic.nima_model.predict(img)
        technical_pred = model_technical.nima_model.predict(img)

        aesthetic_pred = calc_mean_score(aesthetic_pred[0])
        technical_pred = calc_mean_score(technical_pred[0])

        result.append({
            'image': image,
            'aesthetic': aesthetic_pred,
            'technical': technical_pred
        })

    return jsonify(result)

if __name__ == "__main__":
    app.run(host="0.0.0.0")

max_stars_repo_name: samuelbr/image-quality-assessment

max_issues_repo_name: samuelbr/image-quality-assessment

max_forks_repo_name: samuelbr/image-quality-assessment


text:
"""
A debug module for the AMI430 instrument driver. We cannot rely on the physical instrument to be present
which is why we need to mock it.
"""
import re
from datetime import datetime

import numpy as np
import pytest
from hypothesis import given, settings
from hypothesis.strategies import floats
from hypothesis.strategies import tuples

from qcodes.instrument_drivers.american_magnetics.AMI430 import AMI430, AMI430_3D
from qcodes.math.field_vector import FieldVector

field_limit = [  # If any of the field limit functions are satisfied we are in the safe zone.
    lambda x, y, z: x == 0 and y == 0 and z < 3,  # We can have higher field along the z-axis if x and y are zero.
    lambda x, y, z: np.linalg.norm([x, y, z]) < 2
]

@pytest.fixture(scope='module')
def current_driver(request):
    """
    Start three mock instruments representing current drivers for the x, y and z directions.
    """

    driver = AMI430_3D(
        "AMI430-3D",
        AMI430("x", testing=True),
        AMI430("y", testing=True),
        AMI430("z", testing=True),
        field_limit
    )

    return driver

def get_instruments_ramp_messages(current_driver):
    """
    Listen to the mock instruments and parse the messages to extract the ramping targets. This is useful in determining
    if the set targets arrive at the individual instruments correctly. The time stamps are useful to test that the
    messages arrive in the correct order.
    """

    search_string = "\[(.*)\] ([x, y, z]): Ramping to (.*)"
    # We expect the log messages to be in the format "[<time stamp>] <name>: <message>", where name is either x, y or z
    reported_ramp_targets = {}
    messages = current_driver.get_mocker_messages()

    for msg in messages:
        result = re.search(search_string, msg)
        if result is not None:
            time_string, device_name, value = result.groups()
            time_value = datetime.strptime(time_string, "%d:%m:%Y-%H:%M:%S.%f")
            reported_ramp_targets[device_name] = {"value": float(value), "time": time_value}

    return reported_ramp_targets

random_coordinates = {
    "cartesian": tuples(
        floats(min_value=0, max_value=1),  # x
        floats(min_value=0, max_value=1),  # y
        floats(min_value=0, max_value=1)  # z
    ),
    "spherical": tuples(
        floats(min_value=0, max_value=1),  # r
        floats(min_value=0, max_value=180),  # theta
        floats(min_value=0, max_value=180)  # phi
    ),
    "cylindrical": tuples(
        floats(min_value=0, max_value=1),  # rho
        floats(min_value=0, max_value=180),  # phi
        floats(min_value=0, max_value=1)  # z
    )
}

@given(random_coordinates["cartesian"])
@settings(max_examples=10)
def test_cartesian_sanity(current_driver, set_target):
    """
    A sanity check to see if the driver remember vectors in any random configuration in cartesian coordinates
    """
    current_driver.cartesian(set_target)
    get_target = current_driver.cartesian()

    assert np.allclose(set_target, get_target)

    # test that we can get the individual coordinates
    x = current_driver.x()
    y = current_driver.y()
    z = current_driver.z()
    assert np.allclose(set_target, [x, y, z])

@given(random_coordinates["spherical"])
@settings(max_examples=10)
def test_spherical_sanity(current_driver, set_target):
    """
    A sanity check to see if the driver remember vectors in any random configuration in spherical coordinates
    """
    current_driver.spherical(set_target)
    get_target = current_driver.spherical()

    assert np.allclose(set_target, get_target)

    # test that we can get the individual coordinates
    r = current_driver.field()
    theta = current_driver.theta()
    phi = current_driver.phi()
    assert np.allclose(set_target, [r, theta, phi])

@given(random_coordinates["cylindrical"])
@settings(max_examples=10)
def test_cylindrical_sanity(current_driver, set_target):
    """
    A sanity check to see if the driver remember vectors in any random configuration in cylindrical coordinates
    """
    current_driver.cylindrical(set_target)
    get_target = current_driver.cylindrical()

    assert np.allclose(set_target, get_target)

    # test that we can get the individual coordinates
    rho = current_driver.rho()
    z = current_driver.z()
    phi = current_driver.phi()
    assert np.allclose(set_target, [rho, phi, z])

@given(random_coordinates["cartesian"])
@settings(max_examples=10)
def test_cartesian_setpoints(current_driver, set_target):
    """
    Check that the individual x, y, z instruments are getting the set points as intended. We can do this because the
    instruments are printing log messages to an IO stream. We intercept these messages and extract the log lines which
    mention that the instrument is ramping to certain values. These values should match the values of the input. In
    this test we are verifying this for cartesian coordinates.
    """
    current_driver.cartesian(set_target)

    reported_ramp_targets = get_instruments_ramp_messages(current_driver)
    get_target = {k: v["value"] for k, v in reported_ramp_targets.items()}

    set_vector = FieldVector(*set_target)
    get_vector = FieldVector(**get_target)
    assert set_vector.is_equal(get_vector)

@given(random_coordinates["spherical"])
@settings(max_examples=10)
def test_spherical_setpoints(current_driver, set_target):
    """
    Check that the individual x, y, z instruments are getting the set points as intended. We can do this because the
    instruments are printing log messages to an IO stream. We intercept these messages and extract the log lines which
    mention that the instrument is ramping to certain values. These values should match the values of the input. In
    this test we are verifying this for spherical coordinates.
    """
    current_driver.spherical(set_target)

    reported_ramp_targets = get_instruments_ramp_messages(current_driver)
    get_target = {k: v["value"] for k, v in reported_ramp_targets.items()}

    set_vector = FieldVector(**dict(zip(["r", "theta", "phi"], set_target)))
    get_vector = FieldVector(**get_target)
    assert set_vector.is_equal(get_vector)

@given(random_coordinates["cylindrical"])
@settings(max_examples=10)
def test_cylindrical_setpoints(current_driver, set_target):
    """
    Check that the individual x, y, z instruments are getting the set points as intended. We can do this because the
    instruments are printing log messages to an IO stream. We intercept these messages and extract the log lines which
    mention that the instrument is ramping to certain values. These values should match the values of the input. In
    this test we are verifying this for cylindrical coordinates.
    """
    current_driver.cylindrical(set_target)

    reported_ramp_targets = get_instruments_ramp_messages(current_driver)
    get_target = {k: v["value"] for k, v in reported_ramp_targets.items()}

    set_vector = FieldVector(**dict(zip(["rho", "phi", "z"], set_target)))
    get_vector = FieldVector(**get_target)
    assert set_vector.is_equal(get_vector)

@given(random_coordinates["cartesian"])
@settings(max_examples=10)
def test_measured(current_driver, set_target):
    """
    Simply call the measurement methods and verify that no exceptions are raised.
    """
    current_driver.cartesian(set_target)

    cartesian = current_driver.cartesian_measured()
    cartesian_x = current_driver.x_measured()
    cartesian_y = current_driver.y_measured()
    cartesian_z = current_driver.z_measured()

    assert FieldVector(*set_target).is_equal(FieldVector(x=cartesian_x, y=cartesian_y, z=cartesian_z))
    assert np.allclose(cartesian, [cartesian_x, cartesian_y, cartesian_z])

    spherical = current_driver.spherical_measured()
    spherical_field = current_driver.field_measured()
    spherical_theta = current_driver.theta_measured()
    spherical_phi = current_driver.phi_measured()

    assert FieldVector(*set_target).is_equal(FieldVector(
        r=spherical_field,
        theta=spherical_theta,
        phi=spherical_phi)
    )
    assert np.allclose(spherical, [spherical_field, spherical_theta, spherical_phi])

    cylindrical = current_driver.cylindrical_measured()
    cylindrical_rho = current_driver.rho_measured()

    assert FieldVector(*set_target).is_equal(FieldVector(rho=cylindrical_rho, phi=spherical_phi, z=cartesian_z))
    assert np.allclose(cylindrical, [cylindrical_rho, spherical_phi, cartesian_z])

def test_ramp_down_first(current_driver):
    """
    To prevent quenching of the magnets, we need the driver to always be within the field limits. Part of the strategy
    of making sure that this is the case includes always performing ramp downs first. For example, if magnets x and z
    need to be ramped up, while magnet y needs to be ramped down, ramp down magnet y first and then ramp up x and z.
    This is tested here.
    """
    names = ["x", "y", "z"]
    set_point = np.array([0.3, 0.4, 0.5])
    current_driver.cartesian(set_point)  # Put the magnets in a well defined state
    delta = np.array([-0.1, 0.1, 0.1])  # We begin with ramping down x first while ramping up y and z

    for count, ramp_down_name in enumerate(names):
        set_point += np.roll(delta, count)  # The second iteration will ramp down y while ramping up x an z.
        # Check if y is adjusted first. We will perform the same test with z in the third iteration

        current_driver.cartesian(set_point)
        reported_ramp_targets = get_instruments_ramp_messages(
            current_driver)  # get the logging outputs from the instruments.
        times = {k: v["time"] for k, v in reported_ramp_targets.items()}  # extract the time stamps

        for ramp_up_name in np.delete(names, count):
            assert times[ramp_down_name] < times[ramp_up_name]  # ramp down occurs before ramp up.

def test_field_limit_exception(current_driver):
    """
    Test that an exception is raised if we intentionally set the field beyond the limits. Together with the
    no_test_ramp_down_first test this should prevent us from ever exceeding set point limits.
    In this test we generate a regular grid in three-D and assert that the driver can be set to a set point if
    any of of the requirements given by field_limit is satisfied. An error *must* be raised if none of the
    safety limits are satisfied.
    """
    x = np.linspace(-3, 3, 11)
    y = np.copy(x)
    z = np.copy(x)
    set_points = zip(*[i.flatten() for i in np.meshgrid(x, y, z)])

    for set_point in set_points:
        should_not_raise = any([is_safe(*set_point) for is_safe in field_limit])

        if should_not_raise:
            current_driver.cartesian(set_point)
        else:
            with pytest.raises(Exception) as excinfo:
                current_driver.cartesian(set_point)

            assert "field would exceed limit" in excinfo.value.args[0]

def test_cylindrical_poles(current_driver):
    """
    Test that the phi coordinate is remembered even if the resulting vector is equivalent to the null vector
    """
    rho, phi, z = 0.4, 30.0, 0.5
    current_driver.cylindrical((0.0, phi, 0.0))  # This is equivalent to the null vector
    current_driver.rho(rho)
    current_driver.z(z)

    rho_m, phi_m, z_m = current_driver.cylindrical()  # after setting the rho and z values we should have the vector as
    # originally intended.
    assert np.allclose([rho_m, phi_m, z_m], [rho, phi, z])

def test_spherical_poles(current_driver):
    """
    Test that the theta and phi coordinates are remembered even if the resulting vector is equivalent to the null vector
    """

    field, theta, phi = 0.5, 30.0, 50.0
    current_driver.spherical((0.0, theta, phi))  # If field=0, then this is equivalent to the null vector
    # x, y, z = (0, 0, 0)
    current_driver.field(field)

    field_m, theta_m, phi_m = current_driver.spherical()  # after setting the field we should have the vector as
    # originally intended.
    assert np.allclose([field_m, theta_m, phi_m], [field, theta, phi])

def test_warning_increased_max_ramp_rate():
    """
    Test that a warning is raised if we increase the maximum current ramp rate. We want the user to be really sure
    what he or she is doing, as this could risk quenching the magnet
    """
    max_ramp_rate = AMI430.default_current_ramp_limit
    target_ramp_rate = max_ramp_rate + 0.01  # Increasing the maximum ramp rate should raise a warning

    with pytest.raises(Warning) as excinfo:
        AMI430("testing_increased_max_ramp_rate", testing=True, current_ramp_limit=target_ramp_rate)

    assert "Increasing maximum ramp rate" in excinfo.value.args[0]

def test_ramp_rate_exception(current_driver):
    """Test that an exception is raised if we try to set the ramp rate to a higher value then is allowed"""
    max_ramp_rate = AMI430.default_current_ramp_limit
    target_ramp_rate = max_ramp_rate + 0.01
    ix = current_driver._instrument_x

    with pytest.raises(Exception) as excinfo:
        ix.ramp_rate(target_ramp_rate)

    assert "must be between 0 and {} inclusive".format(max_ramp_rate) in excinfo.value.args[0]

max_stars_repo_name: aaroncslau/Qcodes

max_issues_repo_name: qdev-dk/Qcodes

max_forks_repo_name: aaroncslau/Qcodes


text:
"""Code here is only about plotting. There should be no non-trivial
pandas code. There should be no Flask code."""

import io

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import seaborn as sns; sns.set_style('whitegrid')

from matplotlib import cm # color

from .data import my_data

def plot_demo(data_frame):
    """Demonstrates a particular plot for a data frame. Returns tuple of
    image_data and image_type"""
    img = io.BytesIO()

    # matplot/ seaborn style setting
    mpl.rcParams['patch.force_edgecolor'] = True
    sns.set(rc={'figure.figsize': (5.0, 5.0)})

    # seaborn countplot
    p1 = sns.countplot(x='Group', data=data_frame, palette='coolwarm')

    # rotate xticklabels to prevent the labels being overlapped
    p1.set_xticklabels(p1.get_xticklabels(), rotation=90)

    # adjust subplot params
    # to avoid axis labels,titles or ticklabels being clipped
    p1.figure.tight_layout()
    p1.figure.savefig(img, format='png')

    # results
    plt.close()
    image_data = img.getvalue()
    img.close()
    image_type = "image/png"
    return image_data, image_type

def home_bar_plot(data_frame):
    img = io.BytesIO()

    # seaborn style setting
    sns.set_style('whitegrid')
    sns.set(rc={'figure.figsize': (7.0, 4.5)})

    # parameters for grp
    grp = data_frame['grp']
    grp_sizes = data_frame['grp_sizes']

    # seaborn barplot
    p1 = sns.barplot(x='grp', y='grp_sizes',
                     data=data_frame, palette='coolwarm')
    p1.set_xlabel('Group')
    p1.set_ylabel('Total TB')
    p1.set_xticklabels(p1.get_xticklabels(), rotation=90)

    # save plot
    p1.figure.tight_layout()
    p1.figure.savefig(img, format='png')

    # results
    plt.close()
    image_data = img.getvalue()
    img.close()
    image_type = "image/png"
    return image_data, image_type

def home_pie_plot(data_frame):
    img = io.BytesIO()

    # matplot/ seaborn style setting
    sns.set_style('whitegrid')
    mpl.rcParams['patch.force_edgecolor'] = True
    sns.set(rc={'figure.figsize': (2.5, 2.5)})

    # returns a list of values
    appl = data_frame['appl'].tolist()
    appl_sizes = data_frame['appl_sizes'].tolist()

    # pie chart, where the slices will be ordered and plotted counter-clockwise
    labels = appl
    sizes = appl_sizes

    # only "explode" one slice (eg 'Whole EXOME')
    explode = (0, 0.1, 0, 0)

    fig1, ax1 = plt.subplots()
    ax1.pie(sizes, explode=explode, labels=labels, autopct='%1.1f%%',
            labeldistance=1.2, wedgeprops={'linewidth': 0.8, 'edgecolor': 'w'},
            shadow=True, startangle=90)
    ax1.axis('equal')

    # savefig
    fig1.tight_layout()
    plt.savefig(img, format='png', bbox_inches='tight')

    # results
    plt.close()
    image_data = img.getvalue()
    img.close()
    image_type = "image/png"
    return image_data, image_type

def bar_plot(data_frame):
    """input is query data_frame and output is seaborn countplot"""
    img = io.BytesIO()

    # matplot/ seaborn setting
    mpl.rcParams['patch.force_edgecolor'] = True
    sns.set(rc={'figure.figsize': (10, 6)}) # figure size in inches

    # order = data_frane['Group'].value_counts().index # query data
    at = my_data.at # all time data
    order = at['Group'].value_counts().index

    # seaborn countplot
    p1 = sns.countplot(x='Group',
                      data=data_frame,
                      # hue='Group',
                      palette='coolwarm',
                      order = order)

    # rotation xlabels to prevent the labels being overlapped
    p1.set_xticklabels(p1.get_xticklabels(), rotation=-45)

    # adjust subplot params
    # to avoid axis labels, titles or ticklabels being clipped
    p1.figure.tight_layout()
    p1.figure.savefig(img, format='png')

    # put legend out of the figure when hue is turned on
    plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0)

    # results
    plt.close()
    image_data = img.getvalue()
    img.close()
    image_type = "image/png"
    return image_data, image_type

def grp_pie_plot(data_frame):
    img = io.BytesIO()

    # the slices will be ordered and plotted counter-clockwise
    grp = data_frame.groupby('Group')
    grp2 = (data_frame.groupby('Group')['Numeric Total MB'].sum()).reset_index()
    grp2['Total TB'] = (grp2['Numeric Total MB'] / 1_000_000)
    grp3 = grp2
    # grp3 = grp2[grp2['Total TB'] > 100] # filter out negligence

    grp = grp3['Group']
    grp_sizes = grp3['Total TB']

    # parameters for make_pie
    title = "Pie Chart by Group"
    labels = grp
    sizes = grp_sizes
    num_rows = len(grp3)
    explode_index = 0
    angle = 100

    make_pie(title, labels, sizes, num_rows, explode_index, angle)

    # resize the figure box
    # similar to calling figure.tight_layout()
    plt.savefig(img, format='png', bbox_inches='tight')

    # results
    plt.close()
    image_data = img.getvalue()
    img.close()
    image_type = "image/png"
    return image_data, image_type

def appl_pie_plot(data_frame):
    img = io.BytesIO()

    # the slices will be ordered and plotted counter-clockwise:
    appl = data_frame.groupby('Application')
    appl2 = (data_frame.groupby('Application')['Numeric Total MB'].sum()).reset_index()
    appl2['Total TB'] = (appl2['Numeric Total MB'] / 1_000_000)
    appl3 = appl2
    # appl3 = appl2[appl2['Total TB'] > 10] # filter out negligence

    appl = appl3['Application']
    appl_sizes = appl3['Total TB']

    # parameters for make_pie
    title = "Pie Chart by Application"
    labels = appl
    sizes = appl_sizes
    num_rows = len(appl3)
    explode_index = 0
    angle = 110

    make_pie(title, labels, sizes, num_rows, explode_index, angle)

    # resize the figure box
    # similar to calling figure.tight_layout()
    plt.savefig(img, format='png', bbox_inches='tight')

    # results
    plt.close()
    image_data = img.getvalue()
    img.close()
    image_type = "image/png"
    return image_data, image_type

def make_pie(title, labels, sizes, num_rows, explode_index, angle):
    # explode index of the pie slice
    explode = make_explode(num_rows, explode_index)

    fig1, ax1 = plt.subplots()

    # set explode=None to turn off
    ax1.pie(sizes, labels=labels, labeldistance=1.2,
            explode=explode, autopct='%1.1f%%',
            shadow=True, startangle=angle)

    # matplot/ seaborn style setting
    mpl.rcParams['patch.force_edgecolor'] = True
    plt.rcParams["figure.figsize"] = (4, 4) # set fig size in inches

    # fig title
    fig1.suptitle(title, fontsize=13)

    # legend options
    # bbox_to_anchor(x0, y0, width, height)
    plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)

    # set color
    cs=cm.Set1(np.arange(40)/40.)

    # equal aspect ratio ensures that pie is drawn as a circle
    ax1.axis('equal')

    # TODO vertical call out labels

    # TODO correct overlap labels

def make_explode(num_rows, explode_index):
    base = make_base(num_rows)
    base[explode_index] = 0.05
    return tuple(base)

def make_base(length):
    return [0,]*length

max_stars_repo_name: BCM-HGSC/qcweb

max_issues_repo_name: BCM-HGSC/qcweb

max_forks_repo_name: BCM-HGSC/qcweb


text:
#!/usr/bin/env python3

import re

import numpy as np
import pandas as pd

from bs4 import BeautifulSoup

# read in the guidebook-exported sessions file
df_schedule = pd.read_csv('guidebook-exports/schedule.csv')

# read in the guidebook-exported papers file
df_posters = pd.read_csv('guidebook-exports/posters.csv')

# read in the session link template
df_session_links = pd.read_csv('Sessions_Link_Template.csv')
rows = []

# get the HTML from the program page we want to parse
with open('../_pages/testing/program.md', 'r') as htmlf:
    schedule_html = re.search(r'
.*
', htmlf.read(), re.DOTALL).group()

# parse this HTML using BeautifulSoup
soup = BeautifulSoup(schedule_html, 'html.parser')

# go over each paper session, find each paper and create the link from the session
# to its papers
for session_id in ['poster-1', 'poster-2']:
    session_html = soup.find(id="session-{}".format(session_id))
    for poster_html in session_html.find_all("tr", id="poster"):
        poster_title_str, author_str = list(poster_html.find_all("td")[0].stripped_strings)
        poster_title_str = poster_title_str.strip().rstrip('.')

        # get the guidebook custom list ID for this poster title
        guidebook_poster_id = df_posters[df_posters['Name'] == poster_title_str]['Item ID (Optional)'].values[0]
        assert isinstance(guidebook_poster_id, np.int64)

        # now get the guidebook ID for this session
        lookup_session_id = 'Posters & Dinner' if session_id == 'poster-1' else "Posters, Demos & Dinner"
        guidebook_session_id = df_schedule[df_schedule['Session Title'] == lookup_session_id]['Session ID'].values[0]
        assert isinstance(guidebook_session_id, np.int64)

        d = {'Session ID (Optional)': guidebook_session_id,
             'Session Name (Optional)': '',
             'Link To Session ID (Optional)': '',
             'Link To Session Name (Optional)':'',
             'Link To Custom List Item ID (Optional)': guidebook_poster_id,
             'Link To Custom List Item Name (Optional)':'',
             'Link To Form Name (Optional)':'',
             'Link To URLs (Optional)':'',
             'URL Names (Optional)':''}
        rows.append(d)

df_session_links = df_session_links.append(rows)
df_session_links['Session ID (Optional)'] = df_session_links['Session ID (Optional)'].astype(int)
df_session_links['Link To Custom List Item ID (Optional)'] = df_session_links['Link To Custom List Item ID (Optional)'].astype(int)
df_session_links.to_csv('acl2017-session-to-poster-links.csv', index=False)


text:
from numpy.lib.financial import ipmt
from tool.utils import *
from tool.torch_utils import *
from tool.darknet2pytorch import Darknet
import cv2
import argparse

"""hyper parameters"""
use_cuda = False

def detect_BEV(cfgfile, weightfile, imgfile):
    """Detect elements in BEV map with yolov4_BEV

    Args:
        cfgfile (str): Path to .cfg file
        weightfile (str): Path to .weights file
        imgfile (str): Path to image on which we want to run BEV detection
    """

    # load model
    m = Darknet(cfgfile, model_type="BEV_grid")
    m.print_network()
    m.load_weights(weightfile, cut_off=128)
    print("Loading weights from %s... Done!" % (weightfile))

    # push to GPU
    if use_cuda:
        m.cuda()

    # load names
    namesfile = "names/BEV.names"
    class_names = load_class_names(namesfile)

    # read sample image
    img = cv2.imread(imgfile)
    sized = cv2.resize(img, (m.width, m.height))
    sized = cv2.cvtColor(sized, cv2.COLOR_BGR2RGB)

    # create batch
    sized = np.expand_dims(sized, 0)
    sized = np.concatenate((sized, sized), 0)

    # run inference
    start = time.time()
    boxes = do_detect(m, sized, 0.4, 0.6, use_cuda)
    finish = time.time()
    print("%s: Predicted in %f seconds." % (imgfile, (finish - start)))

    # TODO: plot boxes in BEV

def get_args():
    parser = argparse.ArgumentParser("Test your image by trained model.")
    parser.add_argument(
        "-cfgfile",
        type=str,
        default="./cfg/model/yolov4_BEV_grid.cfg",
        help="path of cfg file",
        dest="cfgfile",
    )
    parser.add_argument(
        "-weightfile",
        type=str,
        default="./checkpoints/yolov4.weights",
        help="path of trained model.",
        dest="weightfile",
    )
    parser.add_argument(
        "-imgfile",
        type=str,
        default="../data/KITTI/train/images/000001.png",
        help="path of your image file.",
        dest="imgfile",
    )
    args = parser.parse_args()

    return args

if __name__ == "__main__":
    args = get_args()
    detect_BEV(args.cfgfile, args.weightfile, args.imgfile)