text:
from pyomo.environ import Block, Expression, units as pyunits
from watertap3.utils import financials
from watertap3.wt_units.wt_unit import WT3UnitProcess

## REFERENCE: EPA model

module_name = 'gac'
basis_year = 2012
tpec_or_tic = 'TPEC'

class UnitProcess(WT3UnitProcess):

    def autosize(system_type, flow, num_tanks, num_lines, ebct, ebct_tank, geom='vertical'):

        def autosize_constraints(system_type, flow=flow, num_tanks=num_tanks, num_lines=num_lines, ebct=ebct, geom=geom):

            if system_type == 'pressure':
                # vessel surface area
                comp_min_sa_vessel = comp_sa_min / num_lines  # based on max load, design flow, and target number of lines
                comp_max_sa_vessel = comp_sa_max / num_lines  # based on min load, design flow, and target number of lines

                if geom == 'vertical':
                    # vertical vessel height
                    comp_min_vert_height = 3  # from reasonable dimensions table
                    comp_max_vert_height = 14  # from reasonable dimensions table

                    # vertical bed depth
                    min_depth_a = comp_vol_required / num_tanks / sa_all_vessels  # ensures that bed depth is sufficient to pass diameter check on vessel design sheet
                    min_depth_b = min_bed_depth  # from guidance on input sheet
                    comp_min_vert_bed_depth = max(min_depth_a, min_depth_b)  # strictest of the above
                    max_depth_a = (max_height - freeboard) / (1 + bed_expansion)  # max bed depth for max_height, given freeboard and bed expansion
                    max_depth_b = max_bed_depth  # from guidance on input sheet
                    comp_max_vert_bed_depth = min(max_depth_a, max_depth_b)  # strictest max of the above, rounded

                    # vertical vessel diameter
                    min_diam_a = 2 * ((
                                              comp_vol_required / num_tanks / 3.14159 / comp_max_vert_bed_depth / num_lines) ** 0.5 + vessel_thickness)  # min diam to avoid bed depth becoming
                    # too deep, given required volume
                    min_diam_b = 2 * ((comp_min_sa_vessel / 3.14159) ** 0.5 + vessel_thickness)  # from min surface area
                    min_diam_c = 1.5  # from reasonable dimensions table
                    comp_min_vert_diam = 2 * max(min_diam_a, min_diam_b, min_diam_c), 0 / 2
                    max_diam_a = 14  # max diameter CDA
                    max_diam_b = 2 * ((comp_max_sa_vessel / 3.14159) ** 0.5 + vessel_thickness)  # from max surface area
                    comp_max_vert_diam = max(0.5, 2 * min(max_diam_a, max_diam_b) / 2)  # from reasonable dimensions table

                    return comp_min_vert_height, comp_max_vert_height, comp_min_vert_diam, comp_max_vert_diam, comp_min_vert_bed_depth, comp_max_vert_bed_depth

                if geom == 'horizontal':
                    # horizontal vessel diameter
                    comp_min_horiz_diam = 8  # from guidance provided by Bob Dvorin, 5/16/05
                    comp_max_horiz_diam = min(max_height, max_diam)  # strictest of max_height CDA and max_diam CDA

                    # horizontal vessel length
                    min_length_a = comp_min_sa_vessel / (comp_horiz_diam - 2 * vessel_thickness)  # from max load and actual diameter
                    min_length_b = 20
                    comp_min_horiz_length = max(min_length_a, min_length_b)
                    max_length_a = max_length  # max length CDA
                    max_length_b = comp_max_vol / (comp_horiz_diam * comp_horiz_diam * 3.14159 / 4) / 7.48  # max volume combined with actual diameter
                    comp_max_horiz_length = min(max_length_a, max_length_b)

                    # horizontal bed depth
                    min_depth_a = comp_vol_required / num_tanks / sa_all_vessels  # ensures that bed depth is sufficient for required volume of medium
                    min_depth_b = min_bed_depth  # from media constraint table
                    comp_min_horiz_bed_depth = max(min_depth_a, min_depth_b)
                    max_depth_a = (min(max_height, max_diam) - freeboard) / (1 + bed_expansion)  # max bed depth for max_height or max_diam, given freeboard and bed expansion
                    max_depth_b = max_bed_depth  # from media constraint table
                    comp_max_horiz_bed_depth = min(max_depth_a, max_depth_b)

                    return comp_min_horiz_diam, comp_max_horiz_diam, comp_min_horiz_length, comp_max_horiz_length, comp_min_horiz_bed_depth, comp_max_horiz_bed_depth

            if system_type == 'gravity':
                # bed depth
                min_depth_a = min_depth  # from min_depth CDA
                min_depth_b = flow * ebct / (num_lines * 30 * max_width) / 7.48  # assuming max-size contactor and number of contactors above, min depth for sufficient media volume
                comp_min_bed_depth_a = max(min_depth_a, min_depth_b)
                comp_max_bed_depth_a = max_depth  # from max_depth CDA

                # surface area
                comp_min_sa_a = flow / num_lines / load_max  # from loading: based on design flow per target line, max loading
                comp_max_sa_a = flow / num_lines / load_min  # from loading: based on design flow per target line, min loading

                # contactor length
                comp_min_length_a = 6  # from min_length CDA
                comp_max_length_a = 30  # from max_length CDA

                # contactor width
                comp_min_width_a = min_width  # from min_width CDA
                comp_max_width_a = max_width  # from max_width CDA

                return comp_min_bed_depth_a, comp_max_bed_depth_a, comp_min_sa_a, comp_max_sa_a, comp_min_length_a, comp_max_length_a, comp_min_width_a, comp_max_width_a

            ## PRESSURE ##

        if system_type == 'pressure':

            comp_vol_required_stg1 = flow * ebct_tank / 7.48  # volume of media in first-stage vessels
            comp_sa_min = flow / load_max  # to reach maximum loading rate
            comp_sa_max = flow / load_min  # to reach minimum loading rate
            comp_vol_required = flow * ebct / 7.481  # volume of media in all vessels

            if geom == 'vertical':
                comp_target_bed_depth_vert = target_bed_depth_over  # based on target load CDA
                comp_sa_required_vert = comp_vol_required_stg1 / comp_target_bed_depth_vert  # based on media volume and bed depth
                comp_vert_diam = 2 * ((comp_sa_required_vert / num_lines / 3.14159) ** 0.5 + vessel_thickness)  # target diameter
                sa_one_vessel = (comp_vert_diam - 2 * vessel_thickness) * (comp_vert_diam - 2 * vessel_thickness) * 3.14159 / 4  # surface area, one vessel
                sa_all_vessels = sa_one_vessel * num_lines  # surface area, all vessels

                comp_min_vert_height, comp_max_vert_height, comp_min_vert_diam, comp_max_vert_diam, comp_min_vert_bed_depth, comp_max_vert_bed_depth = autosize_constraints(system_type, geom=geom)

                # vertical design search
                comp_vert_max_media = (comp_max_vert_diam / 2 - vessel_thickness) * (
                        comp_max_vert_diam / 2 - vessel_thickness) * 3.14159 * comp_max_vert_bed_depth  # based on guidance and CDAs for max dimensions
                comp_vert_max_sa = (max_diam / 2 - vessel_thickness) * (max_diam / 2 - vessel_thickness) * 3.14159  # based on max_diam CDA

                num_lines_vert_a = comp_sa_min / comp_vert_max_sa  # based on surface area from max load
                num_lines_vert_b = comp_sa_required_vert / comp_vert_max_sa  # based on surface area from bed depth
                num_lines_vert_c = comp_vol_required_stg1 / comp_vert_max_media  # based on media volume
                comp_vert_min_number = max(num_lines_vert_a, num_lines_vert_b, num_lines_vert_c)

                # inputs

                comp_vert_bed_depth = comp_vol_required / num_tanks / sa_all_vessels  # bed depth
                comp_vert_height = (1 + bed_expansion) * comp_vert_bed_depth + freeboard  # target height
                comp_bed_depth = comp_vert_bed_depth  # bed depth for chosen geometry

                if comp_vert_diam < comp_min_vert_diam:
                    comp_vert_diam = comp_min_vert_diam
                elif comp_vert_diam > comp_max_vert_diam:
                    comp_vert_diam = comp_max_vert_diam
                else:
                    comp_vert_diam = 2 * comp_vert_diam / 2

                if comp_min_vert_bed_depth <= comp_max_vert_bed_depth:
                    if comp_bed_depth < comp_min_vert_bed_depth:
                        comp_bed_depth = comp_min_vert_bed_depth
                    elif comp_bed_depth > comp_max_vert_bed_depth:
                        comp_bed_depth = comp_max_vert_bed_depth
                    else:
                        comp_bed_depth = comp_vert_bed_depth

                if comp_min_vert_height > comp_max_vert_height:
                    comp_vert_height = -1
                    print(f'Min vert height greater than max vert height.\n\tcomp_min_vert_height = {comp_min_vert_height}\n\tcomp_max_vert_height = {comp_max_vert_height}')
                    print('I am not sure what this means yet.')
                elif comp_vert_height > comp_max_vert_height:
                    comp_vert_height = comp_max_vert_height
                elif comp_vert_height < comp_min_vert_height:
                    comp_vert_height = 2 * comp_vert_height / 2  ## NEED TO HAVE A ROUNDUP FUNCTION HERE

                return comp_vert_min_number, comp_vert_height, comp_bed_depth, comp_vert_diam

            if geom == 'horizontal':

                comp_min_horiz_diam, comp_max_horiz_diam, comp_min_horiz_length, comp_max_horiz_length, comp_min_horiz_bed_depth, comp_max_horiz_bed_depth = autosize_constraints(system_type,
                                                                                geom=geom)
                # inputs
                comp_max_vol = 27100  # max pressure vessel volume LOOKED UP
                comp_max_vol_length = comp_max_vol / 7.48 / (3.14159 * comp_max_horiz_diam * comp_max_horiz_diam / 4)  # rounded to reflect the length we will actually use on the input sheet
                comp_target_bed_depth_horiz = target_bed_depth_horiz  # based on target load CDA
                comp_sa_required_horiz = comp_vol_required_stg1 / comp_target_bed_depth_horiz  # based on media volume and bed depth

                fb = comp_max_vol_length * (comp_max_horiz_diam * comp_max_horiz_diam / 2 * np.arcsin((freeboard / comp_max_horiz_diam) ** 0.5) - (comp_max_horiz_diam / 2 - freeboard) * (
                        (comp_max_horiz_diam - freeboard) * freeboard) ** 0.5)  # freeboard - volume of a horizontal slice of a cylinder
                ud = comp_max_vol_length * (
                        comp_max_horiz_diam * comp_max_horiz_diam / 2 * np.arcsin((horiz_underdrain / comp_max_horiz_diam) ** 0.5) - (comp_max_horiz_diam / 2 - horiz_underdrain) * (
                        (comp_max_horiz_diam - horiz_underdrain) * horiz_underdrain) ** 0.5)  # underdrain
                comp_horiz_void_required = fb + ud  # total of the two void volumes above
                comp_horiz_void_area = comp_horiz_void_required / comp_max_vol_length  # assumed void volume per length
                comp_horiz_expanded_bd = comp_target_bed_depth_horiz * (1 + bed_expansion)  # target bed depth after expansion
                comp_horiz_diam = (comp_horiz_expanded_bd + (comp_horiz_expanded_bd * comp_horiz_expanded_bd + 3.14159 * comp_horiz_void_area) ** 0.5) / (3.14159 / 2)  # target diameter
                # comp_horiz_diam_old = (comp_sa_required_horiz / num_lines / length_diam_ratio) ** 0.5
                comp_horiz_length = comp_sa_required_horiz / (comp_horiz_diam - 2 * vessel_thickness) / num_lines  # target length
                sa_all_vessels = num_lines * (comp_horiz_diam - 2 * vessel_thickness) * comp_horiz_length  # surface area for all vessels
                comp_horiz_bed_depth = comp_vol_required / num_tanks / sa_all_vessels
                comp_bed_depth = comp_horiz_bed_depth  # bed depth for chosen geometry

                # horizontal design search

                comp_max_horiz_sa = (comp_max_horiz_diam - 2 * vessel_thickness) * comp_max_vol_length  # max vessel surface area
                comp_target_bed_depth_horiz = target_bed_depth_horiz  # based on target load CDA -- NEEDS TO BE ADJUSTABLE FOR FLOW AND INCLUSION OF UV AOP
                comp_sa_required_horiz = comp_vol_required_stg1 / comp_target_bed_depth_horiz  # total vessel surface area, based upon media volume and ved depth

                comp_horiz_max_media = (comp_max_vol / 7.481 - comp_horiz_void_required) / (1 + bed_expansion)  # max media volume
                num_lines_horiz_a = comp_sa_min / comp_max_horiz_sa  # based on surface area from max load
                num_lines_horiz_b = comp_sa_required_horiz / comp_max_horiz_sa  # based on surface area from bed depth
                num_lines_horiz_c = comp_vol_required_stg1 / comp_horiz_max_media  # based on media volume
                comp_horiz_min_number = max(num_lines_horiz_a, num_lines_horiz_b, num_lines_horiz_c)  # min number of lines, horizontal vessels

                if comp_min_horiz_diam <= comp_max_horiz_diam:
                    if comp_horiz_diam < comp_min_horiz_diam:
                        comp_horiz_diam = comp_min_horiz_diam
                    elif comp_horiz_diam > comp_max_horiz_diam:
                        comp_horiz_diam = comp_max_horiz_diam
                    else:
                        comp_horiz_diam = 2 * comp_horiz_diam / 2  ##  NEED ROUND UP FUNCtiON HERE

                if comp_min_horiz_length <= comp_max_horiz_length:
                    if comp_horiz_length < comp_min_horiz_length:
                        comp_horiz_length = comp_min_horiz_length
                    elif comp_horiz_length > comp_max_horiz_length:
                        comp_horiz_length = comp_max_horiz_length
                    else:
                        comp_horiz_length = comp_horiz_length

                return comp_max_horiz_sa, comp_horiz_length, comp_bed_depth, comp_horiz_min_number

        if system_type == 'gravity':
            comp_min_bed_depth_a, comp_max_bed_depth_a, comp_min_sa_a, comp_max_sa_a, comp_min_length_a, comp_max_length_a, comp_min_width_a, comp_max_width_a = autosize_constraints(system_type,
                                                                                geom=geom)

            comp_vol_required_a = flow * ebct / 7.481  # volume of media in all vessels
            # design search
            comp_vert_max_media_a = comp_max_bed_depth_a * comp_max_length_a * comp_max_width_a  # based on guidance and CDAs for max dimensions
            comp_sa_min_a = flow / load_max  # to reach maximum loading rate
            comp_sa_max_a = flow / load_min  # to reach minimum loading rate
            comp_target_bed_depth_vert_a = target_bed_depth_over  ## based on Target Bed Depth CDA -- NEEDS TO BE ADJUSTABLE FOR SMALLER FLOWS AND INCLUSION (OR NOT) OF UV AOP... SEE
            # DOCUMENTATION AND EXCEL SHEET
            comp_sa_required_vert_a = comp_vol_required_a / comp_target_bed_depth_vert_a  # based on media volume and bed depth
            num_lines_grav_a = comp_min_sa_a / comp_sa_max_a
            num_lines_grav_b = comp_sa_required_vert_a / comp_sa_max_a
            num_lines_grav_c = comp_vol_required_a / comp_vert_max_media_a
            comp_vert_min_number_a = min(num_lines_grav_a, num_lines_grav_b, num_lines_grav_c)  # below this number of lines, we are constrained to horizontal vessels

            comp_length_width_a = (comp_sa_required_vert_a / num_lines) ** 0.5  # target side length
            sa_one_vessel = comp_length_width_a ** 2  # surface area one vessel
            sa_all_vessels = num_lines * sa_one_vessel  # surface area all vessels
            comp_bed_depth_a = comp_vol_required_a / num_tanks / sa_all_vessels  # bed depth

            if comp_length_width_a < comp_min_length_a:
                comp_length_width_a = comp_min_length_a
            elif comp_length_width_a > comp_max_length_a:
                comp_length_width_a = 2 * comp_length_width_a / 2  ## NEED ROUND UP FUNCTION HERE=

            return comp_length_width_a, comp_bed_depth_a, comp_vert_min_number_a

    def get_costing(self, unit_params=None, year=None):
        financials.create_costing_block(self, basis_year, tpec_or_tic)
        t = self.flowsheet().config.time
        time = self.flowsheet().config.time.first()
        self.ebct = Var(t, initialize=unit_params['ebct'],
                        domain=NonNegativeReals,
                        bounds=(1, 30),  # units=pyunits.dimensionless,
                        doc='ebct')
        self.flow_in = pyunits.convert(self.flow_vol_in[time], to_units=pyunits.Mgallons / pyunits.day)

        self.system_type = unit_params['system_type']
        # regen = unit_params['regen']
        self.regen = 'onsite'
        try:
            self.freund1 = unit_params['freund1']  # Kf (ug/g)(L/ug)1/n
            self.freund2 = unit_params['freund2']  # dimensionless
            self.conc_in = unit_params['conc_in']
            self.conc_breakthru = unit_params['conc_breakthru']
        except:
            # defaults to values for Atrazine and 99% reduction (1 mg/L --> 0.01 mg/L)
            self.freund1 = 38700
            self.freund2 = 0.291
            self.conc_in = 1
            self.conc_breakthru = 0.01

        # if system_type == 'pressure':
        # geom should be determined by autosize function
        # or should be able to be determined by autosize function
        try:
            self.geom = unit_params['geom']
            self.pv_material = unit_params['pv_material']
            self.bw_tank_type = unit_params['bw_tank_type']
        except:
            self.geom = 'vertical'
            self.pv_material = 'stainless'
            self.bw_tank_type = 'stainless'

        ############## SYSTEM INPUTS ##############
        #     flow_in = 1
        self.bp_pct = 0  # * pyunits.dimensionless # bypass percent of flow
        self.comm_diam = 8.5  # * pyunits.ft # vessel diameter ??
        self.vessel_thickness = 0  # * pyunits.ft # thickness of vessel walls ?
        self.comm_height_length = 11  # * pyunits.ft # 'straight' height of vessel ??
        self.bed_depth = 6.7  # * pyunits.ft # vessel bed depth
        self.backwash_load_rate = 3  # * (pyunits.gallon / pyunits.minute) / pyunits.ft**2 # backwash loading rate
        self.backwash_time = 12  # * pyunits.minutes # backwashing time
        self.backwash_interval = 168  # * pyunits.hours # interval between backwash occurrences, depends on many things # 48 hours for surface water if loading rate 4-6 gpm/sf; 7+ days for ground water
        # (Richard E. Hubel, Peter Keenan, American Water Works Service Company, Inc).  2 weeks or loss of head based on operational experience (Verna Arnette, Greater Cincinnati Water Works).
        self.regen_list = ['onsite', 'offsite', 'throw_nonhaz', 'offsite_haz', 'throw_haz', 'throw_rad', 'throw_hazrad']
        self.pv_material_list = ['stainless', 'carbon with stainless', 'carbon with plastic', 'fiberglass']
        #     regen = 'onsite'

        ########################################################
        ########################################################
        ############## CRITICAL DESIGN ASSUMPTIONS ##############
        ########################################################
        ########################################################

        ############## GAC CONTACTOR DESIGN ##############
        self.bed_expansion = 0.5  # * pyunits.dimensionless # resin expansion during backwash
        self.load_min = 0.5  # gpm / ft2, # Hydraulic surface loading rate determines headloss and is usually kept between 2 -10 gpm/sf. Higher loading rates OK if hydraulics and tests allow. Typical
        # hydraulic design criteria are a surface loading rate of 2-10 gpm/ft2 (Montgomery, James M.  Water Treatment Principles and Design.  New York: John Wiley and Sons, 1985. Pp. 553.)  Most
        # regulatory agencies will not approve rates in excess of 2 gpm/ft2 without successful pilot-scale testing (American Water Works Association and American Society of Civil Engineers. Water
        # Treatment Plant Design. Third Edition. New York: McGraw Hill, 1998. Pp. 158.).
        self.load_max = 10  # gpm / ft2, # Many systems can operate to 0.5 gpm/sf (Richard E. Hubel, Peter Keenan, American Water Works Service Company, Inc.).

        ############## PRESSURE VESSEL DESIGN ##############
        if self.system_type == 'pressure':
            self.max_height = 14  # * pyunits.ft
            self.max_length = 53  # * pyunits.ft
            self.freeboard = 0.5  # * pyunits.ft # free board above media at full expansion
            self.horiz_underdrain = 1.5  # * pyunits.ft # underdrain space for horizontal vessels
            self.length_diam_ratio = 2.5
            self.min_bed_depth = 2  # * pyunits.ft # Typical pressure GAC bed depths are 2-8.5 ft
            self.max_bed_depth = 8.5  # * pyunits.ft
            self.target_bed_depth_under = 4  # * pyunits.ft # pressure designs < 1MGD
            self.target_bed_depth_over = 7  # * pyunits.ft # pressure designs > 1MGD
            self.target_bed_depth_horiz = 8  # * pyunits.ft # horizontal designs

        ############## GRAVITY BASIN DESIGN ##############
        if self.system_type == 'gravity':
            self.max_length = 30  # * pyunits.ft
            self.min_length = 6  # * pyunits.ft
            self.max_width = 30  # * pyunits.ft
            self.min_width = 6  # * pyunits.ft
            self.max_depth = 10  # * pyunits.ft
            self.min_depth = 3  # * pyunits.ft # Typical gravity GAC bed depths are 3-10 ft
            self.freeboard = 2  # * pyunits.ft # free board above media at full expansion
            self.target_bed_depth_under = 6  # * pyunits.ft # gravity designs < 1MGD
            self.target_bed_depth_uv_aop = 5.4  # * pyunits.ft # gravity with UV AOP no matter flow
            self.target_bed_depth_over = 8  # * pyunits.ft # pressure designs > 1MGD

        ############## GENERAL FOR ALL BASINS ##############
        self.num_stairs = 15  # only if concrete used
        self.overexcavate = 4  # * pyunits.ft # 2 feet each side in each dimension. This assumption used only if concrete basins are used.
        self.conc_thick = 1  # * pyunits.ft # This assumption used only if concrete basins are used

        ############## PUMPING ##############
        self.pump_safety_factor = 0.25  # * pyunits.dimensionless # safety factor for pumps
        self.num_redund_pumps = 1  # number of redundant pumps
        self.num_redund_back_pumps = 1
        self.num_pumps = 1

        ############## BACKWASH ##############
        self.backwash_rate = 12  # * pyunits.gallon / pyunits.minute / pyunits.ft ** 2 # Fluidization occurs at 2-4 gpm/sf; 12-15 gpm/sf is needed to 'expand' GAC and wash out particulates [Calgon
        # Carbon] (Richard E. Hubel, Peter Keenan, American Water Works Service Company, Inc)
        self.backwash_time = 10  # * pyunits.minute # Backwash duration is time allowed for backwash of filter. Use highest expected duration. Typically 5-10 minutes. Default is 10 minutes.
        self.back_multiplier = 1  # volume # number of backwash ctcles to provide a storage volume for clean water
        self.back_ratio = 0.6  # only if concrete basins are used
        self.back_freeboard = 3  # * pyunits.ft # only if concrete basins are  used

        ############## GAC REGENERATION AND TRANSFER ##############
        self.gac_density = 30  # * pyunits.lb / pyunits.ft ** 3 # Bulk density, not particle density. GAC density range of 22-30 lb/cubic ft; lignite GAC approximately 22 and bituminous GAC
        # approximately 28 lb/cubic ft (Randal J. Braker, General Manager, Duck River Utility Commission)
        self.regen_runtime = 0.85  # Percent of time that regeneration facility operates. Less than full time (100%) to allow for routine maintenance.
        self.regen_redund = 1  # 100% means no redundancy (i.e., no additional capacity beyond that needed to account for runtime). Use a value greater than 100% to add capacity and provide redundancy.
        self.makeup_rate = 0.1  # Makeup carbon 10% for on-site regeneration (Richard E. Hubel, Peter Keenan, American Water Works Service Company).
        self.makeup_rate_off = 0.3  # Makeup carbon needs are higher for transport and regeneration off-site -- 30% (Richard E. Hubel, Peter Keenan, American Water Works Service Company).
        self.transfer_time = 4  # * pyunits.hours # Time to empty GAC bed; for mechanical transfer only
        self.transfer_time_hold = 4  # * pyunits.hours # Time to empty holding tank/basin soilds; for mechanical transfer only
        self.virgin_gac_storage = 1.33  # bed volumes
        self.spent_gac_storage = 1.33
        self.regen_gac_storage = 1.33
        self.storage_basins = 1  # Default: one each for virgin, spent, and regenerated GAC. This assumption applies if concrete basins are used. If steel tanks are used, number calculated on
        # regeneration sheet.
        self.storage_ratio = 0.6  # only if concrete used
        self.storage_freeboard = 3  # * pyunits.ft # only if concrete used

        ########################################################
        ########################################################
        ############## CONTACTOR CONSTRAINTS ##############
        ########################################################
        ########################################################

        ############## FROM INPUTS ##############
        self.design_flow = self.flow_in  # MGD
        self.average_flow = (1 - self.bp_pct) * self.design_flow  # MGD
        # design_flow_rate = design_flow * 1E6 / (24 * 60)  ## convert to gpm  .convert(design_flow, to_units=* (pyunits.gallon / pyunits.min)
        self.design_flow_rate = pyunits.convert(self.design_flow, to_units=(pyunits.gallon / pyunits.min))
        # average_flow_rate = average_flow * 1E6 / (24 * 60)  ## convert to gpm .convert(design_flow, to_units=* (pyunits.gallon / pyunits.min)
        self.average_flow_rate = pyunits.convert(self.average_flow, to_units=(pyunits.gallon / pyunits.min))

        self.bp_flow_rate = self.bp_pct / (1 - self.bp_pct) * self.design_flow_rate  # gpm - bypass flowrate
        self.target_contaminant = 'test'  #
        #     conc_in = 10  # * (pyunits.mg / pyunits.L) # influent concentration of target contaminant
        #     conc_out = 1  # * (pyunits.mg / pyunits.L) # target effluent concentration of target contaminant
        #     conc_breakthru = 5  # * (pyunits.mg / pyunits.L) # breakthrough concentration
        #     kss = 0.01  # * pyunits.hour ** -1 # steady state system constant
        #     ebct = (np.log(conc_in / conc_out) / kss) * 60  # * pyunits.min # calculated EBCT based on Kss - conc_out/conc_in = exp(-kss * ebct)
        #     ebct = 4  # * pyunits.min # empty bed contact time - INPUT

        if self.system_type == 'pressure':
            ############## PRESSURE VESSELS ##############
            self.num_vessels = 1

            ############## FROM CENTRAL COST DATABASE ##############
            self.max_diam = 14  # * pyunits.ft # max diameter for vessel
            self.min_diam = 1.5  # * pyunits.ft # min diameter for vessel

            ############## VESSEL REQUIREMENTS BASED ON BED DEPTH AND CONTACT TIME ##############
            self.contact_material = self.design_flow_rate * self.ebct[time] / 7.481  # * pyunits.ft ** 3 # minimum contact material needed

            ############## QUANTITIES BASED ON USER-SPECIFIED VESSEL DIMENSIONS ##############

            self.comm_vol = 3.14159 * comm_diam ** 2 / 4 * comm_height_length  # ft3 - total vessel volume

            ############## CHECKS ON USER-SPECIFIED DIMENSIONS ##############
            #### ::: DESIGN DIMENSIONS FOR VERTICAL VESSELS ::: ####
            self.height = self.bed_depth * (1 + self.bed_expansion) + self.freeboard  # * pyunits.ft # check on required height of vessel
            self.min_comm_height = 6  # * pyunits.ft # minimum commercial height given diameter - NEED TO DETERMINE DYNAMICALLY
            self.max_comm_height = 20  # * pyunits.ft # maximum commercial height given diameter - NEED TO DETERMINE DYNAMICALLY

            #### ::: VOID SPACE REQUIRED FOR HORIZONTAL VESSELS ::: ####
            self.freeboard_horiz_void = self.comm_height_length * (
                    self.comm_diam * self.comm_diam / 2 * np.arcsin((self.freeboard / self.comm_diam) ** 0.5) - (self.comm_diam / 2 - self.freeboard) * ((self.comm_diam - self.freeboard) * self.freeboard) ** 0.5)
            self.underdrain_horiz_void = self.comm_height_length * (self.comm_diam * self.comm_diam / 2 * np.arcsin((self.horiz_underdrain / self.comm_diam) ** 0.5) - (self.comm_diam / 2 - self.horiz_underdrain) * (
                    (self.comm_diam - self.horiz_underdrain) * self.horiz_underdrain) ** 0.5)  # * pyunits.ft**3
            self.horiz_void = self.freeboard_horiz_void + self.underdrain_horiz_void  # * pyunits.ft**3

            #### ::: VOLUME AVAILABLE FOR MEDIA AND EXPANSION ::: ####
            self.vert_vol_avail = (self.comm_height_length - self.freeboard) * 3.14159 * self.comm_diam ** 2 / 4  # * pyunits.ft**3 - volume available for vertical vessel
            self.horiz_vol_avail = self.comm_vol - self.horiz_void  # * pyunits.ft**3 - volume available for horizontal vessel
            self.vol_tol = 0.40  # * pyunits.dimensionless # tolerance for oversized vessels; possibly could determine dynamically based on tank geometry and flow
            self.ebct_tank = self.ebct[time] / self.num_vessels
            self.distance_above = 0  # pyunits.ft # height available above media bed for horizontal vessels
            self.void_above = 0  # * pyunits.ft**3 # void space above media bed for horizontal vessels
            if self.geom == 'vertical':
                self.max_sa = ((self.comm_diam - (2 * self.vessel_thickness)) ** 2 * 3.14159 / 4)  # * pynits.ft**2 # surface area of vessel- depends on tank_geom
                self.media_vol = (self.max_sa * self.bed_depth)  # * pyunits.ft ** 3 # GAC volume per contactor
                self.flow_per_vessel = (self.media_vol * 7.48) / self.ebct_tank  # flow per vessel # depends on tank_geom
            if self.geom == 'horizontal':
                self.max_sa = ((self.comm_diam - (2 * vessel_thickness)) * self.comm_height_length)  # * pynits.ft**2 # surface area of vessel- depends on tank_geom
                self.media_vol = (self.comm_vol - self.underdrain_horiz_void - self.void_above)  # * pyunits.ft ** 3 # GAC volume per contactor
                self.flow_per_vessel = (self.media_vol * 7.48) / self.ebct_tank
            # num_treat_lines = math.ceil(design_flow_rate / flow_per_vessel)  # * pyunits.dimensionless # number of treatment trains
            self.num_treat_lines = self.design_flow_rate / self.flow_per_vessel  # * pyunits.dimensionless # number of treatment trains

            #     num_treat_lines = design_flow_rate / flow_per_vessel  # need to determine dynamically; see above line
            # self.op_num_tanks = math.ceil(num_treat_lines * num_vessels)  # * pyunits.dimensionless # number of vessels excluding redundancy
            self.op_num_tanks = self.num_treat_lines * self.num_vessels  # * pyunits.dimensionless # number of vessels excluding redundancy

            self.num_redund_tanks = 1  # number redundant vessles # NEED TO FIGURE OUT HOW TO CALCULATE DYNAMICALLY
            self.tot_num_tanks = self.op_num_tanks + num_redund_tanks  # * pyunits.dimensionless # total number of vessels needed
            self.contact_mat_per_vessel = self.contact_material / self.op_num_tanks
            self.min_sa = self.contact_mat_per_vessel / self.bed_depth
            self.diameter = (self.min_sa * 4 / 3.14159) ** 0.5 + (2 * self.vessel_thickness)  # pyunits.ft # check on required diameter of vessel
            self.gac_each = self.media_vol * self.gac_density  # * pyunits.lb # GAC quantity each contactor
            self.gac_total = self.gac_each * self.tot_num_tanks
            self.gac_expand = self.media_vol * (1 + self.bed_expansion)

        if self.system_type == 'gravity':
            ############## GRAVITY CONTACT BASINS ##############
            self.min_basin_vol = design_flow_rate * self.ebct[time] / 7.48  # * pyunits.ft ** 3
            self.basin_freeboard = 2  # * pyunits.ft # Basin freeboard above filter media at full expansion
            self.basin_width = 8  # * pyunits.ft
            self.basin_length = self.basin_width
            self.basin_op_depth = 7.9  # * pyunits.ft
            self.basin_depth = self.basin_op_depth * (1 + self.bed_expansion) + self.basin_freeboard
            self.num_redund_basins = 1  # number redundant basins # NEED TO FIGURE OUT HOW TO CALCULATE DYNAMICALLY
            # self.op_num_basins = math.ceil(self.min_basin_vol / (self.basin_width * self.basin_length * basin_op_depth))
            self.op_num_basins = self.min_basin_vol / (self.basin_width * self.basin_length * self.basin_op_depth)

            self.tot_num_basins = self.num_redund_basins + self.op_num_basins
            self.basin_conc_vol = (2 * (self.basin_width * self.tot_num_basins + self.conc_thick * (self.tot_num_basins + 1)) * self.basin_depth * self.conc_thick + (
                    self.tot_num_basins + 1) * self.basin_length * self.basin_depth * self.conc_thick + (self.basin_length + 2 * self.conc_thick) * (
                                      self.basin_width * self.tot_num_basins + self.conc_thick * (self.tot_num_basins + 1)) * self.conc_thick) / 27  # * pyunits.yard ** 3 # concrete vol for contact basins
            self.excavation_vol = ((self.tot_num_basins * self.basin_width + (self.tot_num_basins + 1) * conc_thick + overexcavate) * (self.basin_length + 2 * self.conc_thick + self.overexcavate) * (
                    self.basin_depth + self.conc_thick + self.overexcavate) + 0.5 * (3 * (self.basin_depth + self.conc_thick + self.overexcavate) * (self.basin_depth + self.conc_thick + self.overexcavate) * (
                    self.basin_length + 2 * conc_thick + overexcavate))) / 27  # * pyunits.yard ** 3 # concrete excavation vol for contact basins
            self.backfill_vol = excavation_vol - (self.tot_num_basins * self.basin_width + (self.tot_num_basins + 1) * self.conc_thick) * (self.basin_length + 2 * self.conc_thick) * (self.basin_depth + self.conc_thick) / 27
            self.railing = 2 * (self.tot_num_basins * self.basin_width + 2 * self.conc_thick / 2 + (self.tot_num_basins - 1) * self.conc_thick) + 2 * (
                    self.basin_length + 2 * self.conc_thick / 2)  # * pyunits.ft # railing length for contact basins
            self.basin_area = self.basin_width * self.basin_length  # * pyunits.ft ** 2 # surface area for each contactor
            self.loading_rate = self.design_flow_rate / (self.op_num_basins * (self.basin_width * self.basin_length))  # * pyunits.gallon / pyunits.minute / pyunits.ft ** 2 # basin loading rate
            self.media_vol = self.basin_width * self.basin_op_depth * self.basin_length
            self.gac_each = self.media_vol * self.gac_density  # * pyunits.lb # GAC quantity each contactor
            self.gac_total = self.gac_each * self.tot_num_basins
            self.gac_expand = self.media_vol * (1 + self.bed_expansion)

        ############## MEDIA QUANTITIES ##############

        self.media_vol_expanded = self.media_vol * (1 + self.bed_expansion)

        ########################################################
        ########################################################
        ############## BACKWASH AND REGENERATION ##############
        ########################################################
        ########################################################

        ############## FROM INPUTS ##############
        self.bw_pumps = 1  # variable to determine method of backwashing
        self.bw_tanks = 1  # variable to determine method of backwashing
        self.storage_basins = 1  # how many storage basins do you want for each of virgin, spent, regen GAC

        ############## BACKWASH ##############
        if self.system_type == 'pressure':
            self.water_flush = self.backwash_rate * self.max_sa  # * (pyunits.gallons / pyunits.minute) # backwash flow rate
        if self.system_type == 'gravity':
            self.water_flush = self.backwash_rate * self.basin_width * self.basin_length  # * (pyunits.gallons / pyunits.minute) # backwash flow rate
        self.total_backwash = water_flush * backwash_time  # * pyunits.gallons # total backwash volume
        self.backwash_storage = total_backwash * back_multiplier / 7.48  # * pyunits.ft ** 3 # backwash storage volume
        self.num_back_basins = 1  # default value, number of backwash basins
        self.back_basin_op_depth = (backwash_storage / num_back_basins) ** (1 / 3) * back_ratio  # * pyunits.ft # operational depth
        self.back_basin_width = ((total_backwash / num_back_basins) / back_basin_op_depth) ** 0.5
        self.back_basin_length = back_basin_width
        self.back_basin_depth = back_basin_op_depth + back_freeboard  # * pyunits.ft # total depth including freeboard
        self.back_basin_vol = back_basin_width * back_basin_length * back_basin_depth * 7.48  # pyunits.gallons # backwash basin volume
        self.back_conc_vol = (2 * (back_basin_width * num_back_basins + conc_thick * (num_back_basins + 1)) * back_basin_depth * conc_thick + (
                num_back_basins + 1) * back_basin_length * back_basin_depth * conc_thick + (back_basin_length + 2 * conc_thick) * (
                                 back_basin_width * num_back_basins + conc_thick * (num_back_basins + 1)) * conc_thick) / 27  # * pyunits.yards ** 3 # Concrete volume
        back_excavation_vol = ((num_back_basins * back_basin_width + (num_back_basins + 1) * conc_thick + overexcavate) * (back_basin_length + 2 * conc_thick + overexcavate) * (
                back_basin_depth + conc_thick + overexcavate) + 0.5 * (3 * (back_basin_depth + conc_thick + overexcavate) * (back_basin_depth + conc_thick + overexcavate) * (
                back_basin_length + 2 * conc_thick + overexcavate))) / 27  # * pyunits.yards ** 3 # Excavation volume
        back_backfill_vol = back_excavation_vol - (num_back_basins * back_basin_width + (num_back_basins + 1) * conc_thick) * (back_basin_length + 2 * conc_thick) * (
                back_basin_depth + conc_thick) / 27  # * pyunits.yards ** 3 # Backfill volume
        back_railing = 2 * (num_back_basins * back_basin_width + 2 * conc_thick / 2 + (num_back_basins - 1) * conc_thick) + 2 * (back_basin_length + 2 * conc_thick / 2)  # pyunits.ft
        max_back_tank_size = 282294  # * pyunits.gallons # maximum backwash tank size - NEED TO DETERMINE DYNAMICALLY
        num_back_tanks = total_backwash * back_multiplier / max_back_tank_size  # Number of backwash tanks

        back_tank_vol = total_backwash * back_multiplier / num_back_tanks  # * pyunits.gallons # Backwash tank volume

        ############## REGENERATION ##############
        regen_list = ['onsite', 'offsite', 'throw_nonhaz', 'offsite_haz', 'throw_haz', 'throw_rad', 'throw_hazrad']
        flow_regen = average_flow_rate * 3.785  # * pyunits.L / pyunits.min # flow rate for Freundlich calculation
        n0e = self.freund1 * gac_density * 453.59 / 0.02832  # * pyunits.mg / pyunits.m ** 3 # for BDST calculation

        if system_type == 'pressure':
            mass_regen = (self.gac_each * num_treat_lines) * 453.59  # * pyunits.gram # mass of GAC for Freundlich calculation
            velocity_regen = flow_regen * 0.001 / (self.max_sa * num_treat_lines * 0.093)  # m / min # for BDST calculation
            bed_depth_regen = bed_depth * 0.3048  # * pyunits.m
            # __ Calculates bed life based on Freundlich isotherm
            self.bed_life = (mass_regen * (self.freund1 / 1000) * (self.conc_in * 1000) ** self.freund2 / (flow_regen * (self.conc_in - self.conc_breakthru))) / (30 * 24 * 60)  # * pyunits.months #
            #         self.bed_life = self.freund1
            #         self.bed_life = (n0e / (1000 * velocity_regen * self.conc_in) * bed_depth_regen - np.log(self.conc_in / self.conc_breakthru - 1) / (self.freund2 * self.conc_in)) / (30*24*60)
            self.avg_regen_rate = self.gac_each * self.op_num_tanks / (self.bed_life / 12) / 365  # * pyunits.lb / pyunits.day

        if system_type == 'gravity':
            mass_regen = self.gac_total * 453.49  # * pyunits.gram # mass of GAC for Freundlich calculation
            velocity_regen = flow_regen * 0.001 / (self.basin_width * self.basin_length * self.tot_num_basins * 0.093)  # m / min # for BDST calculation
            bed_depth_regen = basin_op_depth * 0.3048  # * pyunits.m
            # __ Calculates bed life based on Freundlich isotherm
            self.bed_life = (mass_regen * (self.freund1 / 1000) * (self.conc_in * 1000) ** self.freund2 / (flow_regen * (self.conc_in - self.conc_breakthru))) / (30 * 24 * 60)  # * pyunits.months #
            self.avg_regen_rate = self.gac_each * self.op_num_basins / (self.bed_life / 12) / 365  # * pyunits.lb / pyunits.day

        self.gac_annual = self.avg_regen_rate * 365  # *pyunits.lb / pyunits.year # When regenerating, this quantity is made up of regenerated GAC + makeup GAC

        if regen == 'onsite':
            regen_capacity = (self.avg_regen_rate / regen_runtime) * regen_redund  # * pyunits.lb / pyunits.day # Regeneration capacity needed (on-site only)
            max_furnace_size = 100000  # * pyunits.lb / pyunits.day # maximum furnace size, looked up from table
            num_furnace = regen_capacity / max_furnace_size
            furnace_size = regen_capacity / num_furnace
            regen_annual = (1 - makeup_rate) * self.gac_annual  # * pyunits.lbs / pyunits.year # rate if regenerated on-site
            gac_makeup = makeup_rate * self.gac_annual  # * pyunits.lbs / pyunits.year # makeup rate for onsite regeneration
        if regen == 'offsite':
            regen_annual = (1 - makeup_rate_off) * self.gac_annual  # * pyunits.lbs / pyunits.year # rate if regenerated off-site
            gac_makeup = makeup_rate_off * self.gac_annual  # * pyunits.lbs / pyunits.year # makeup rate for offsite regeneration
        if regen == 'throwaway':
            gac_makeup = self.gac_annual  # * pyunits.lbs / pyunits.year # makeup rate for GAC throwaway

        ############## GAC STORAGE ##############
        # _____ USING THE SAME STORAGE TANK FOR VIRGIN, SPENT, AND/OR REGENERATED GAC _____
        gac_storage_vol = (self.gac_each / gac_density) * virgin_gac_storage  # * pyunits.ft ** 3 # GAC storage on-site only
        storage_basin_depth = (gac_storage_vol / storage_basins) ** (1 / 3) * storage_ratio  # * pyunits.ft # stage basin depth
        storage_basin_width = ((gac_storage_vol / storage_basins) / storage_basin_depth) ** 0.5  # * pyunits.ft
        storage_basin_length = storage_basin_width  # * pyunits.ft
        storage_basin_vol = storage_basin_width * storage_basin_length * storage_basin_depth * 7.48  # * pyunits.gallons
        storage_conc_vol = ((
                                    2 * storage_basin_depth * storage_basin_length + 2 * storage_basin_depth * storage_basin_width + storage_basin_width * storage_basin_length) * 1.2 * conc_thick *
                            storage_basins) / 27  # * pyunits.yard ** 3
        storage_excavation_vol = ((storage_basin_width + overexcavate) * (storage_basin_length + overexcavate) * (storage_basin_depth + overexcavate) + 0.5 * (
                3 * (storage_basin_depth + overexcavate) * (storage_basin_width + overexcavate) * (storage_basin_depth + overexcavate))) * storage_basins / 27  # * pyunits.yard ** 3
        storage_backfill_vol = storage_excavation_vol - (storage_basin_width * storage_basin_length * storage_basin_depth * storage_basins) / 27  # * pyunits.yard ** 3
        storage_railing = (storage_basin_width * 2 + storage_basin_length * 2) * storage_basins * 2  # * pyunits.ft
        num_storage_tanks = gac_storage_vol * 7.48 / max_back_tank_size
        storage_tank_vol = gac_storage_vol * 7.48 / num_storage_tanks

        ########################################################
        ########################################################
        ############## PUMPS PIPES STRUCTURE ##############
        ########################################################
        ########################################################

        ############## PUMP AND TRANSFER EQUIPMENT  ##############
        transfer_method = 'eductors'
        incl_backwash_pumps = True
        num_redund_eductors = 1

        max_pump_size = 10000  # * pyunits.gallons / pyunits.minute # largest pump size for which costs are available

        pump_rating = (design_flow_rate + bp_flow_rate) * (1 + pump_safety_factor) / num_pumps  # * pyunits.gallons / pyunits.minute # booster pump rating
        num_pumps = None

        if not num_pumps:
            if system_type == 'pressure':
                pump_rating = 0  # pyunits.gallons / pyunits.minute
            elif system_type == 'gravity':
                pump_rating_a = average_flow_rate + bp_flow_rate * average_flow_rate / design_flow_rate
                pump_rating_b = (design_flow_rate + bp_flow_rate) / 2
                pump_rating_c = max(pump_rating_a, pump_rating_b) * (1 + pump_safety_factor)
                pump_rating = min(pump_rating_c, max_pump_size)  # pyunits.gallons / pyunits.minute
        elif num_pumps == 0:
            pump_rating = 0  # pyunits.gallons / pyunits.minute
        elif num_pumps > 0:
            pump_rating = (design_flow_rate + bp_flow_rate) * (1 + pump_safety_factor) / num_pumps  # pyunits.gallons / pyunits.minute
        max_cm_pump_size = 360  # * (pyunits.gallons / pyunits.hour) # The largest size for which costs are available.
        if pump_rating == 0:
            booster_pumps = 0
            total_booster_pumps = 0
        else:
            booster_pumps = (design_flow_rate + bp_flow_rate) * (1 + pump_safety_factor) / pump_rating
            total_booster_pumps = booster_pumps + num_redund_pumps

        back_pump_flow_total = water_flush * (1 + pump_safety_factor)  # * pyunits.gallon / pyunits.minute # backwash pumping flow needed
        if not incl_backwash_pumps:
            num_back_pumps = 0
            back_pump_rating = 0
        else:
            num_back_pumps = back_pump_flow_total / max_pump_size + num_redund_back_pumps
            back_pump_rating = back_pump_flow_total / (num_back_pumps - num_redund_back_pumps)  # * pyunits.gallons / pyunits.minute # backwash pump rating

        if transfer_method == 'eductors':
            transfer_rate = self.gac_each / transfer_time
            max_eductor_size = 54 * gac_density / 7.48 * 60  # * (pyunits.lbs / pyunits.hour) # The largest size for which costs are available. Assumption used only if tanks are used.
            eductors = (transfer_rate / max_eductor_size) + num_redund_eductors

            # The regression for eductors is a linear regression made by fitting the following data to a linear trendline in Excel.
            # Taken from 'Lookup Tables' sheet in the EPA GAC model 'wbs-gac-020818.xls'
            #
            # x = [0, 1673.459893, 3153.406417, 5343.245989, 7942.176471, 9434.15508, 12875.33155] # Transfer Rate (lbs/hr)
            # y = [1.5, 2, 2.5, 3, 4, 6, 8] # Diameter (inch)
            #
            # y = m * x + b
            # diameter = 0.0005 * transfer_rate + 0.9764
            # R_squared = 0.9381

            eductor_size = 0.0005 * transfer_rate + 0.9764  # * pyunits.inches # eductor size NEED TO LOOK UP IN eductor_size_table but don't know where that is

        ########################################################
        ########################################################
        ############## RESIDUALS MANAGEMENT ##############
        ########################################################
        ########################################################

        # discharge = ['surface_water', 'potw', 'recycle', 'septic', 'evaporation_pond']
        discharge = 'surface_water'
        holding_tank = True
        res_cap_factor = 2  # Holding tank/basin capacity safety factor
        dredging_freq = 1  # * pyunits.year * -1 # Evaporation pond solids removal frequency -- Range from 1 to 3 yrs, depending on the final solid concentration required and local climate (
        # Dewatering Municipal Wastewater Sludge,EPA design Manual, Dewatering Municipal Wastewater Sludges.EPA/625/1-87/014, September 1987)
        potw_tss_limit = 250  # * (pyunits.mg / pyunits.liter) # POTW TSS discharge limit (over which fee would be imposed) -- Most common limit for cities with a limit

        ############## SPENT BACKWASH ##############
        res_vol = total_backwash  # * pyunits.gallons # Single event volume -- backwashing one vessel

        if system_type == 'pressure':
            res_stagger_time = backwash_interval * 60 / self.tot_num_tanks  # * pyunits.minute # Time between events
            res_flow_annual = res_vol * self.tot_num_tanks * (365 * 24 / backwash_interval)  # * pyunits.gallons # annual residuals discharge
        elif system_type == 'gravity':
            res_stagger_time = backwash_interval * 60 / self.tot_num_basins  # * pyunits.minute # Time between events
            res_flow_annual = res_vol * self.tot_num_basins * (365 * 24 / backwash_interval)  # * pyunits.gallons # annual residuals discharge

        if holding_tank or discharge == 'recycle':
            res_flow_actual = res_vol / res_stagger_time  # * (pyunits.gallons / pyunits.minute) # Actual residuals discharge flow rate -- Accounts for flow equalization if holding tanks used
            res_flow = res_cap_factor * res_vol / res_stagger_time
        else:
            res_flow_actual = res_vol / backwash_time  # * (pyunits.gallons / pyunits.minute) # Actual residuals discharge flow rate -- Accounts for flow equalization if holding tanks used
            res_flow = res_flow_actual

        ############## EVAPORATION POND OR HOLDING/SEPTIC TANK SOLIDS ##############
        if discharge in ['septic', 'evaporation_pond']:
            tss_in = pyunits.convert(self.conc_mass_in[time, 'tss'], to_units=(pyunits.mg / pyunits.liter))
            tss_out = pyunits.convert(self.conc_mass_out[time, 'tss'], to_units=(pyunits.mg / pyunits.liter))
            if system_type == 'pressure':
                tss_residuals = (tss_in - tss_out) * average_flow * 1E6 * (backwash_interval / 24) / (res_vol * self.tot_num_tanks)  # * (pyunits.mg / pyunits.liter) # TSS in residuals
            elif system_type == 'gravity':
                tss_residuals = (tss_in - tss_out) * average_flow * 1E6 * (backwash_interval / 24) / (res_vol * self.tot_num_basins)  # * (pyunits.mg / pyunits.liter) # TSS in residuals
            if discharge == 'evaporation_pond':
                ep_solids = tss_residuals * res_flow_annual * 3.785 / 1000 / 453.59 / dredging_freq  # * pyunits.lbs # Evaporation pond TSS accumulation

        ############## SPENT MEDIA ##############
        disp_freq = 12 / self.bed_life  # times per year, replacement/disposal frequency
        spent_media = self.gac_total / 2000  # * pyunits.tons # Spent media quantity (total per event)
        spent_media_yr = spent_media * disp_freq  # * pyunits.tons # Spent media quantity (average per year)

        ############## ::: CAPITAL ITEMS ::: ##############

        ############## HOLDING TANKS (FOR FLOW EQUALIZATION) ##############
        max_htank_size = 282094  # * pyunits.gallons # Maximum holding tank size -- The largest size for which costs are available.
        # htanks = # number of holding tanks
        # htank_vol = # * pyunits.gallons # Holding tank volume

        ############## HOLDING BASINS (FOR FLOW EQUALIZATION) ##############

        # holding_ratio = 0.6  # Holding basin depth to volume ratio --
        # holding_freeboard = 3 # * pyunits.ft # holding basin freeboard
        # hmix_hp_per = 0.25 # * (pyunits.hp / (1000 * pyunits.gallons)) # Holding tank/basin mixer horsepower (if coagulant is used)

        # hbasin_redund = # number redundant holding basins
        # hbasins = # number holding basins
        # total_hvol_ft = res_cap_factor * res_vol * hbasin_redund / 7.48 # * pyunits.ft ** 3 # Total residuals storage volume in cubic feet
        # hvol_ft = total_hvol_ft / hbasins # * pyunits.ft ** 3 # Total residuals storage volume in cubic feet
        # hbasin_op_depth = hvol_ft ** (1 / 3) * holding_ratio # * pyunits.ft ** 3 # holding basin operating depth
        # hbasin_width = (hvol_ft / hbasin_op_depth) ** 0.5 # * pyunits.ft # Holding basin width  ## NEEDS TO BE ROUNDED
        # hbasin_length = hbasin_width # * pyunits.ft # Holding basin length
        # hbasin_depth = hbasin_op_depth + holding_freeboard # * pyunits.ft # Holding basin depth including freeboard
        # hbasin_vol = hbasin_width * hbasin_length * hbasin_depth * 7.48 # * pyunits.ft # Holding basin volume  ## NEEDS TO BE ROUNDED
        # hconc_vol = (2 * (hbasin_width * hbasins + conc_thick * (hbasins + 1)) * hbasin_depth * conc_thick + (hbasins + 1) * hbasin_length * hbasin_depth * conc_thick + (
        #             hbasin_length + 2 * conc_thick) * (hbasin_width * hbasins + conc_thick * (hbasins + 1)) * conc_thick) / 27 # * pyunits.yard ** 3 # Concrete volume
        # hexcavation_vol = ((hbasins * hbasin_width + (hbasins + 1) * conc_thick + overexcavate) * (hbasin_length + 2 * conc_thick + overexcavate) * (hbasin_depth + conc_thick + overexcavate) +
        # 0.5 * (
        #             3 * (hbasin_depth + conc_thick + overexcavate) * (hbasin_depth + conc_thick + overexcavate) * (hbasin_length + 2 * conc_thick + overexcavate))) / 27 # * pyunits.yard ** 3 #
        #             Excavation volume
        # hbackfill_vol = hexcavation_vol - (hbasins * hbasin_width + (hbasins + 1) * conc_thick) * (hbasin_length + 2 * conc_thick) * (hbasin_depth + conc_thick) / 27 # * pyunits.yard ** 3 #
        # Backfill volume
        # hrailing = 2 * (hbasins * hbasin_width + 2 * conc_thick / 2 + (hbasins - 1) * conc_thick) + 2 * (hbasin_length + 2 * conc_thick / 2) # * pyunits.ft # Railing length

        ############## FERRIC CHLORIDE ADDITION ##############

        ############## POLYMER ADDITION ##############

        ############## SOLIDS TRANSFER FOR HOLDING TANKS ##############
        # res_transfer_rate = res_solids * 2000 / htransfer_time # * (pyunits.lbs / pyunits.hr) # transfer rate required
        # res_transfer_rate_gal = res_transfer_rate / (ht_solids_density / 7.48 * 60) # * (pyunits.gallons / pyunits.minute) # transfer rate in gpm
        # res_slurry_pumps = 1 # number of slurry pump systems
        # res_eductors = res_transfer_rate_gal / max_eductor_size # number of eductors
        # res_eff_transfer_rate = res_transfer_rate_gal / res_eductors # * (pyunits.gallons / pyunits.minute) # Effective transfer rate
        # res_eductor_size =

        ############## SOLIDS DRYING PAD ##############
        # solids_pad_area =
        # solids_pad_concrete =

        ############## SEPTIC SYSTEM ##############
        # max_sept_tank_vol =
        # num_sept_tanks =
        # sept_tank_vol =
        # sept_design_flow =
        # sept_infil_area =
        # sept_trench_bot_area =
        # sept_trenches =
        # sept_trench_I =
        # sept_gravel_vol =
        # sept_trench_excav_vol =
        # sept_tank_side =
        # sept_tank_excav_vol =
        # sept_dist_boxes =
        # sept_pipe_length =
        # sept_fp =

        ############## EVAPORATION POND ##############
        # ep_total_sa =
        # ep_solids_depth =
        # ep_water_depth =
        # ep_depth =
        # ep_cells =
        # ep_cell_side =
        # ep_excavation_vol =
        # ep_backfill_vol =
        # dike_volume =
        # liner_area =

        ############## ::: O & M ITEMS ::: ##############

        ############## PUMPS ##############

        ############## COAGULANT USE ##############

        ############## COAGULANT PURCHASING ##############

        ############## POLYMER PURCHASING ##############

        ############## POTW DISCHARGE FEES ##############
        tss_residuals_dis = tss_residuals  # THERE IS MORE TO DETERMINING THIS VARIABLE BUT NOT SURE IF IT APPLIES
        tss_test = ((tss_residuals_dis - potw_tss_limit) * res_flow_annual * 3.785 / 1000 / 453.59)
        if tss_test > 0:
            ex_tss_total = tss_test  # * (pyunits.lbs / pyunits.year) # Total TSS in excess of POTW limit (if applicable)
        else:
            ex_tss_total = 0  # * (pyunits.lbs / pyunits.year) # Total TSS in excess of POTW limit (if applicable)
        # POTW_base_cost_avg_cl = 15.33 $/month
        # POTW_base_cost_typ_cl = 17.25 $/month
        # POTW_TSS_cost_avg_cl = 0.24 $/lb over limit
        # POTW_TSS_cost_typ_cl = 0.41 $/lb over limit
        # POTW_vol_cost_avg_cl = 4.33 $/ 1000 gallons
        # POTW_vol_cost_typ_cl = 4.46 $/ 1000 gallons

        potw_base = 17.25 * 12  # Annual POTW base fee $ / year
        potw_vol = 0.41 * res_flow_annual / 1000  # Annnual POTW volume fee $ / year
        potw_tss = 4.46 * ex_tss_total  # Annual POTW TSS fee $ / year
        potw_fee = potw_base + potw_vol + potw_tss  # Total Annual POTW discharge fee $ / year

        ############## SPENT MEDIA DISPOSAL FEES ##############
        # Constants from cost lookup tables and inputs from input sheet or design assumptions
        max_waste_ship = 18  # * pyunits.tons # Maximum waste shipment size
        non_haz_miles = 10  # * pyunits.miles # Miles to non-hazardous solid waste disposal site
        haz_miles = 200  # * pyunits.miles # Miles to hazardous solid waste disposal site
        rad_miles = 700  # * pyunits.miles # Miles to radioactive non-hazardous solid waste disposal site
        hazrad_miles = 700  # * pyunits.miles # Miles to radioactive hazardous solid waste disposal site

        nonhaz_disposal_cost = 68.20  # $/ton Non-hazardous waste disposal
        nonhaz_transportation_cost = 0.51  # $/ton*mile Non-hazardous waste transportation
        haz_disposal_cost = 357.60  # $/ton Hazardous waste disposal
        haz_trans_min_cost = 3067.09  # $/shipment Hazardous waste minimum charge per shipment
        haz_transportation_cost = 0.11  # $/ton*mile Hazardous waste transportation
        hazrad_annual_fee = 1585.06  # $/year Radioactive hazardous waste annual fee
        hazrad_disposal_cost = 11043.39  # $/ton Radiactive hazardous waste disposal
        hazrad_trans_min_cost = haz_trans_min_cost
        hazrad_transportation_cost = 0.26  # $/ton*mile radiactive hazardous waste transportation
        rad_annual_fee = hazrad_annual_fee  # $/year radiactive waste annual fee
        rad_disposal_cost = 5145.16  # $/ton radiactive waste disposal
        rad_transportation_cost = hazrad_transportation_cost  # $/ton*mile radioactive waste transportation
        rad_trans_min_cost = hazrad_trans_min_cost

        regen_list = ['onsite', 'offsite', 'throw_nonhaz', 'offsite_haz', 'throw_haz', 'throw_rad', 'throw_hazrad']

        # 1 for regeneration onsite, 2 for regeneration off-site (non-hazardous), 3 for throwaway (non-hazardous), 4 for regeneration off-site (hazardous), 5 for throwaway (hazardous),
        # 6 for throwaway (radioactive), 7 for throwaway (hazardous & radioactive)
        if spent_media > max_waste_ship:
            disposal_shipments = spent_media / max_waste_ship * disp_freq  # Shipments per year ##  NEEDS TO BE ROUNDED UP SOMEHOW
        else:
            disposal_shipments = disp_freq

        if regen in ['onsite', 'offsite', 'offsite_haz']:
            disposal = 0  # $/year -- cost for disposal
            transport_unit = 0  # $/ton*mile -- transportation unit cost
            transport_min = 0  # $/shipment -- transportation minimum charge per shipment
            disposal_miles = 0  # miles -- transportation mileage
        if regen == 'throw_nonhaz':
            disposal = spent_media * disp_freq * nonhaz_disposal_cost  # $/year -- cost for disposal
            transport_unit = nonhaz_transportation_cost  # $/ton*mile -- transportation unit cost
            transport_min = 0  # $/shipment -- transportation minimum charge per shipment
            disposal_miles = non_haz_miles  # miles -- transportation mileage
        if regen == 'throw_haz':
            disposal = spent_media * disp_freq * haz_disposal_cost  # $/year -- cost for disposal
            transport_unit = haz_transportation_cost  # $/ton*mile -- transportation unit cost
            transport_min = haz_trans_min_cost  # $/shipment -- transportation minimum charge per shipment
            disposal_miles = haz_miles  # miles -- transportation mileage
        if regen == 'throw_rad':
            disposal = spent_media * disp_freq * rad_disposal_cost  # $/year -- cost for disposal
            transport_unit = rad_transportation_cost  # $/ton*mile -- transportation unit cost
            transport_min = rad_trans_min_cost  # $/shipment -- transportation minimum charge per shipment
            disposal_miles = rad_miles  # miles -- transportation mileage
        if regen == 'throw_hazrad':
            disposal = spent_media * disp_freq * hazrad_disposal_cost  # $/year -- cost for disposal
            transport_unit = hazrad_transportation_cost  # $/ton*mile -- transportation unit cost
            transport_min = hazrad_trans_min_cost  # $/shipment -- transportation minimum charge per shipment
            disposal_miles = hazrad_miles  # miles -- transportation mileage

        transportation = max(disposal_shipments * transport_min, spent_media * disp_freq * disposal_miles * transport_unit)  # $ maximum transportation fee

        if regen in ['onsite', 'offsite', 'throw_nonhaz', 'offsite_haz', 'throw_haz']:
            disposal_fee = disposal + transportation  # $/year -- total annual spent media disposal fee
        else:
            disposal_fee = disposal + transportation + hazrad_annual_fee  # $/year -- total annual spent media disposal fee

        ########################################################
        ########################################################
        ############## O & M ##############
        ########################################################
        ########################################################

        ############## ENERGY ##############
        # pump_hp =
        # pump_energy =
        #
        # back_pump_psi =
        # back_pump_hp =
        # back_hrs_yr =
        # back_pump_energy =

        ############## REGEN ENERGY ##############
        # regen_energy =
        # regen_ng =
        #

        pump_head = 44  # * pyunits.ft # pump head for booster pumps, can be adjusted by user input
        back_pump_head = 33  # * pyunits.ft # backwash pump head
        res_pump_psi = 25  # * pyunits.psi # residual pump head
        pump_efficiency = 0.75  # * pyunits.dimensionless # pump efficiency
        if bp_pct > 0:
            pump_hp = average_flow_rate / bp_pct * pump_head * 0.43351482 / booster_pumps / pump_efficiency / 14.7 * 33.9 * 8.34 / 550 / 60

        transfer_rate = self.gac_each / transfer_time  # * pyunits.lb / pyunits.hour # transfer rate required
        max_eductor_size = 12784  # * pyunits.lb / pyunits.hour # maximum eductor size
        num_redund_eductors = 1
        eductors = transfer_rate / max_eductor_size + num_redund_eductors

        # eductor_size = transfer_rate / eductors - num_redund_eductors) ---> look up this size in size table

        # _______________________________________________________  # _______________________________________________________  ############## ACTUALLY CALCULATING COST ##############  #
        # _______________________________________________________  # _______________________________________________________

        # CONTACTORS  ## PRESSURE VESSELS
        def fixed_cap():
            def gac_cost_curves(eqn, x):
                cost_df = pd.read_csv('data/gac_cost_eqns.csv', index_col='eqn')
                cost_df.drop(columns=['pct_deviation', 'date_modified', 'r_squared', 'max_size', 'min_size'], inplace=True)
                coeffs = dict(cost_df.loc[eqn].items())
                cost = coeffs['C1'] * x ** coeffs['C2'] + coeffs['C3'] * log(x) + coeffs['C4'] + coeffs['C5'] * exp(coeffs['C6'] * x) + coeffs['C7'] * x ** 3 + coeffs['C8'] * x ** 2 + coeffs[
                    'C9'] * x + coeffs['C10']
                return cost

            pv_material_list = ['stainless', 'carbon with stainless', 'carbon with plastic', 'fiberglass']

            if system_type == 'pressure':
                #         comm_vol = comm_vol * 7.48  # gal - converting from ft3 to gal
                pv_ss_cost = gac_cost_curves('ss_pv_eq', comm_vol * 7.48)  # cost of stainless steel pressure vessel
                pv_cs_cost = gac_cost_curves('cs_pv_eq', comm_vol * 7.48)  # cost of carbon steel pressure vessels with stainless internals
                pv_csp_cost = gac_cost_curves('csp_pv_eq', comm_vol * 7.48)  # cost of carbon steel pressure vessels with plastic internals
                pv_fg_cost = gac_cost_curves('fg_pv_eq', comm_vol * 7.48)  # cost of fiberglass pressure vessels
                if pv_material == 'stainless':
                    pv_cost = pv_ss_cost * self.tot_num_tanks
                if pv_material == 'carbon with stainless':
                    pv_cost = pv_cs_cost * self.tot_num_tanks
                if pv_material == 'carbon with plastic':
                    pv_cost = pv_csp_cost * self.tot_num_tanks
                if pv_material == 'fiberglass':
                    pv_cost = pv_fg_cost * self.tot_num_tanks

            ## GAC CONTACT BASINS
            if system_type == 'gravity':
                concrete_contact_cost = basin_conc_vol * 582.09  # 'conc_basin_cost_cl' cost per cubic yard of concrete
                first_underdrain_cost = gac_cost_curves('cont_bot_eq', self.basin_area)
                add_underdrain_cost = gac_cost_curves('cont_bot_add_eq', self.basin_area)
                back_dist_cost = gac_cost_curves('cont_top_eq', self.basin_length)
                internals_cost = first_underdrain_cost + (self.tot_num_basins - 1) * add_underdrain_cost + self.tot_num_basins * back_dist_cost
                railing_cost = (railing * 36)  # 'metal_cost_cl' cost per linear foot
                stairs_cost = (num_stairs * 541.43)  # 'metal_cost_cl' cost per riser
                excavation_cost = (excavation_vol * 30.08)  # 'excavate_cost_cl' cost per cubic yard of excavation
                backfill_compact_cost = (backfill_vol * 13.95)  # 'backfill_cost_cl' cost per cubic yard of backfill

                basin_cost = concrete_contact_cost + internals_cost + railing_cost + stairs_cost + excavation_cost + backfill_compact_cost

            # TANKS
            ## BACKWASH BASINS/TANKS

            #     bw_tank_type_list = ['concrete', 'stainless', 'fiberglass', 'hdpe']
            #     bw_tank_type = 'stainless'

            concrete_back_cost = num_back_basins * (back_conc_vol * 582.09 + back_excavation_vol * 30.08 + back_backfill_vol * 13.95 + back_railing * 36)
            bw_ss_cost = gac_cost_curves('st_bwt_eq', back_tank_vol)
            bw_fg_cost = gac_cost_curves('fg_bwt_eq', back_tank_vol)
            bw_hdpe_cost = gac_cost_curves('hdpe_bwt_eq', back_tank_vol)

            if bw_tank_type == 'concrete':
                bw_tank_cost = concrete_back_cost * num_back_basins
            if bw_tank_type == 'stainless':
                bw_tank_cost = bw_ss_cost * num_back_tanks
            if bw_tank_type == 'fiberglass':
                bw_tank_cost = bw_fg_cost * num_back_tanks
            if bw_tank_type == 'hdpe':
                bw_tank_cost = bw_hdpe_cost * num_back_tanks

            ## HOLDING TANKS
            # ht_ss_cost = gac_cost_curves('st_bwt_eq', back_tank_vol)
            # ht_fb_cost = gac_cost_curves('fg_bwt_eq', back_tank_vol)
            # ht_hdpe_cost =

            # PUMPS
            booster_pump_cost = gac_cost_curves('pump_eq', pump_rating)
            bw_pump_cost = gac_cost_curves('back_pump_eq', back_pump_rating)
            # res_pump_cost = gac_cost_curves('pump_eq', res_pump_rating)

            #  MEDIA AND CHEMICALS
            ## INITIAL GAC CHARGE
            gac_cost = gac_cost_curves('gac_eq', self.gac_total)
            # fecl3_cost =

            # ONSITE GAC REGENERATION
            ## MULTIPLE HEARTH FURNACE SYSTEM
            # furnace_cost = gac_cost_curves('on_regen_eq', furnace_size)

            ## STORAGE
            storage_ss_cost = gac_cost_curves('st_bwt_eq', storage_tank_vol) * 3  # multiplied by 3 to account for virigin, spent, and regen GAC storage
            # storage_conc_cost =

            ##### TOTAL SYSTEM COST #####
            if system_type == 'pressure':
                total_system_cost = pv_cost + bw_tank_cost + booster_pump_cost + bw_pump_cost + gac_cost + storage_ss_cost

            if system_type == 'gravity':
                total_system_cost = basin_cost + bw_tank_cost + booster_pump_cost + bw_pump_cost + gac_cost + storage_ss_cost

            return total_system_cost * 1E-6

        self.chem_dict = {}

        def electricity(flow_in):  # m3/hr
            flow_in = pyunits.convert(flow_in, to_units=(pyunits.m ** 3 / pyunits.year))
            electricity = (750 * 1E3) / flow_in  # kWh/m3
            return electricity

        self.costing.fixed_cap_inv_unadjusted = Expression(expr=fixed_cap(), doc='Unadjusted fixed capital investment')  # $M

        self.electricity = electricity(flow_in)  # kwh/m3

        financials.get_complete_costing(self.costing)

reason: from pyomo


text:
import pandas as pd
import quandl, math
import numpy as np
from sklearn import preprocessing, cross_validation, svm
from sklearn.linear_model import LinearRegression

# get Google stock data set
df = quandl.get('WIKI/GOOGL')

df = df[['Adj. Open', 'Adj. High', 'Adj. Low', 'Adj. Close', 'Adj. Volume',]]
df['HL_PCT'] = (df['Adj. High'] - df['Adj. Close']) / df['Adj. Close'] * 100
df['PCT_change'] = (df['Adj. Close'] - df['Adj. Open']) / df['Adj. Open'] * 100

df = df[['Adj. Close', 'HL_PCT', 'PCT_change', 'Adj. Volume']]

forecast_col = 'Adj. Close'
# Because ML can not work with na data
df.fillna(-99999, inplace=True)
# math.ceil(x) = smallest integer value greater than or equal to x.
# We are try to predict out 10 percent of the dataframe and you'll see that actually
# when will go out and do this.
# number of predict
forecast_out = int(math.ceil(0.01*len(df)))
print (forecast_out)
df['label'] = df[forecast_col].shift(-forecast_out)
# print (df['label'])
# print (df[forecast_col])
# delect the the column which value is NA
df.dropna(inplace=True)
# print (df['label'])

# print (df.head())

X = np.array(df.drop(['label'], 1))
y = np.array(df['label'])
# scale data to "Gaussian with zero mean and unit variance"
X = preprocessing.scale(X)

# 20% of data we want to actually use as testing data.
# Generate the training and test data
X_train, X_test, y_train, y_test = cross_validation.train_test_split(X, y, test_size=0.2)
# classifier
clf = LinearRegression(n_jobs=-1)
# clf = svm.SVR()
# clf = svm.SVR(kernel='poly')
clf.fit(X_train, y_train)
accuracy = clf.score(X_test, y_test)
print (accuracy)

ext: py

reason: import numpy


text:
import os
import time
import datetime
import pathlib
import json
import random

import torch
import numpy as np
import cv2
import carla
from PIL import Image
import weakref
import pickle

from team_code.map_agent import MapAgent
from team_code.pid_controller import PIDController
from contrastiveFC.data import scale_and_crop_image, transform_2d_points, lidar_to_histogram_features
from contrastiveFC.config import GlobalConfig
from contrastiveFC.models import ContrastiveLearningModel_ImageOnly, ControlsModel_linear_ImageOnly, ImitationLearningModel_ImageOnly
#from contrastiveFC.models3 import ImitationLearningModel_ImageOnly
from leaderboard.autoagents.agent_wrapper import AgentWrapper
import leaderboard.autoagents.agent_wrapper

SAVE_PATH = os.environ.get('SAVE_PATH', None)

WEATHERS = {
        'ClearNoon': carla.WeatherParameters.ClearNoon,
}
WEATHERS_IDS = list(WEATHERS)

def get_entry_point():
    return 'ContrastiveFC_Agent'

def _numpy(carla_vector, normalize=False):
    result = np.float32([carla_vector.x, carla_vector.y])

    if normalize:
        return result / (np.linalg.norm(result) + 1e-4)

    return result

def _location(x, y, z):
    return carla.Location(x=float(x), y=float(y), z=float(z))

def _orientation(yaw):
    return np.float32([np.cos(np.radians(yaw)), np.sin(np.radians(yaw))])

def get_collision(p1, v1, p2, v2):
    A = np.stack([v1, -v2], 1)
    b = p2 - p1

    if abs(np.linalg.det(A)) < 1e-3:
        return False, None

    x = np.linalg.solve(A, b)
    collides = all(x >= 0) and all(x <= 1) # how many seconds until collision

    return collides, p1 + x[0] * v1

def check_episode_has_noise(lat_noise_percent, long_noise_percent):
    lat_noise = False
    long_noise = False
    if random.randint(0, 101) < lat_noise_percent:
        lat_noise = True

    if random.randint(0, 101) < long_noise_percent:
        long_noise = True

    return lat_noise, long_noise

class ContrastiveFC_Agent(MapAgent):

    # for stop signs
    PROXIMITY_THRESHOLD = 30.0  # meters
    SPEED_THRESHOLD = 0.1
    WAYPOINT_STEP = 1.0  # meters

    def setup(self, path_to_conf_file):
        super().setup(path_to_conf_file)
        self.weather_id = None

        # Config file path for contrastiveFC (specs for data transforms)
        self.config_path = path_to_conf_file
        self.config = GlobalConfig()

        # Contrastive+FC Model paths
        cvm_path = '/home/yewon/new-data-dl/models/34_100to0split_nonlinearity3layers/bestContrastive_valLoss.pt'
        clm_path = '/home/yewon/new-data-dl/models/34_100to0split_nonlinearity3layers/bestControl_valLoss.pt'

        # Load trained models
        self.net_contrastive = ContrastiveLearningModel_ImageOnly().cuda()
        ###self.net_contrastive = torch.nn.DataParallel(self.net_contrastive)
        self.net_contrastive.load_state_dict(torch.load(cvm_path))
        #self.net_contrastive = torch.nn.DataParallel(self.net_contrastive)
        self.net_contrastive.eval()
        ##self.net_contrastive = torch.nn.DataParallel(self.net_contrastive)
        self.net_controls = ControlsModel_linear_ImageOnly(self.net_contrastive).cuda()
        #self.net_controls = ImitationLearningModel_ImageOnly().cuda()
        #self.net_controls = torch.nn.DataParallel(self.net_controls)
        self.net_controls.load_state_dict(torch.load(clm_path))
        #self.net_controls = torch.nn.DataParallel(self.net_controls)
        self.net_controls.eval()
        ##self.net_controls = torch.nn.DataParallel(self.net_controls)
        self.prev_lane_intervention = 0

        if SAVE_PATH is not None:
            now = datetime.datetime.now()
            string = pathlib.Path(os.environ['ROUTES']).stem + '_'
            string += '_'.join(map(lambda x: '%02d' % x, (now.month, now.day, now.hour, now.minute, now.second)))

            print(string)

            # Create the save directories
            self.save_path = pathlib.Path(os.environ['SAVE_PATH']) / string
            self.save_path.mkdir(parents=True, exist_ok=True)
            (self.save_path / 'rgb_front').mkdir(parents=True, exist_ok=True)
            (self.save_path / 'lidar').mkdir(parents=True, exist_ok=True)
            (self.save_path / 'measurements').mkdir(parents=True, exist_ok=True)

        self.start_time = time.time()

    def _init(self):
        super()._init()

        # Toggle this to set autopilot/contrastiveFC to control
        self.autopilot = False
        self.reversed_time = None

        # Initialize metrics
        self.cte = 0
        self.autonomy = 0
        self.prev_cte = 0 # needed for when vehicle has no waypoint measurements
        self.collided = False
        self.autopilot_next = False # directs run_step to use autopilot in the next run step if previous was a collision/reverse drive

        # Initialize history of times, ctes, and interventions
        self.history = {"time":[],
        				"cte":[],
        				"intervention":[]}

        # Initialize variables for intervention
        self.prev_intervention_time = 0
        self.takeovertime = 6 #3 # what we used in original code base
        self.prev_wp = None
        self.interventions = 0
        self.episode_time_len = 0

        self._turn_controller = PIDController(K_P=1.25, K_I=0.75, K_D=0.3, n=40)
        self._speed_controller = PIDController(K_P=5.0, K_I=0.5, K_D=1.0, n=40)

        # For stop signs
        self._target_stop_sign = None # the stop sign affecting the ego vehicle
        self._stop_completed = False # if the ego vehicle has completed the stop sign
        self._affected_by_stop = False # if the ego vehicle is influenced by a stop sign

        # Lane invasion sensor
        bp = self._world.get_blueprint_library().find('sensor.other.lane_invasion')
        lane_invasion_sensor = self._world.spawn_actor(bp, carla.Transform(), self._vehicle)
        weak_self = weakref.ref(self)
        lane_invasion_sensor.listen(lambda event: ContrastiveFC_Agent._on_invasion(weak_self, event))

    def _on_invasion(weak_self, event):
        self = weak_self()
        if not self:
            return
        print("Lane invasion occured")
        #text = ['%r' % str(x).split()[-1] for x in set(event.crossed_lane_markings)]
        #self._hud.notification('Crossed line %s' % ' and '.join(text))

    def _get_angle_to(self, pos, theta, target):
        R = np.array([
            [np.cos(theta), -np.sin(theta)],
            [np.sin(theta),  np.cos(theta)],
            ])

        aim = R.T.dot(target - pos)
        angle = -np.degrees(np.arctan2(-aim[1], aim[0]))
        angle = 0.0 if np.isnan(angle) else angle

        return angle

    def _get_control(self, target, far_target, tick_data):
        pos = self._get_position(tick_data)
        theta = tick_data['compass']
        speed = tick_data['speed']

        # Steering.
        angle_unnorm = self._get_angle_to(pos, theta, target)
        angle = angle_unnorm / 90

        steer = self._turn_controller.step(angle)
        steer = np.clip(steer, -1.0, 1.0)
        steer = round(steer, 3)

        # Acceleration.
        angle_far_unnorm = self._get_angle_to(pos, theta, far_target)
        should_slow = abs(angle_far_unnorm) > 45.0 or abs(angle_unnorm) > 5.0
        target_speed = 4.0 if should_slow else 7.0
        brake = self._should_brake()
        target_speed = target_speed if not brake else 0.0

        delta = np.clip(target_speed - speed, 0.0, 0.25)
        throttle = self._speed_controller.step(delta)
        throttle = np.clip(throttle, 0.0, 0.75)

        if brake:
            steer *= 0.5
            throttle = 0.0

        return steer, throttle, brake, target_speed

    def _get_control_autopilot(self, target, far_target, tick_data):

        """Same as _get_control() but want to import this into ConrolsFC_agent"""

        pos = self._get_position(tick_data)
        theta = tick_data['compass']
        speed = tick_data['speed']

        # Steering.
        angle_unnorm = self._get_angle_to(pos, theta, target)
        angle = angle_unnorm / 90

        steer = self._turn_controller.step(angle)
        steer = np.clip(steer, -1.0, 1.0)
        steer = round(steer, 3)

        # Acceleration.
        angle_far_unnorm = self._get_angle_to(pos, theta, far_target)
        should_slow = abs(angle_far_unnorm) > 45.0 or abs(angle_unnorm) > 5.0
        target_speed = 4.0 if should_slow else 7.0
        brake = self._should_brake()
        target_speed = target_speed if not brake else 0.0

        delta = np.clip(target_speed - speed, 0.0, 0.25)
        throttle = self._speed_controller.step(delta)
        throttle = np.clip(throttle, 0.0, 0.75)

        if brake:
            steer *= 0.5
            throttle = 0.0

        return steer, throttle, brake, target_speed

    def crossTrackErrorNew(self, vehicle_wp, ref_wp):
    	min_distance = 10000
    	distance = np.linalg.norm(ref_wp - vehicle_wp)
    	min_distance = min(distance, min_distance)
    	return min_distance

    def has_collided(self, speed):
    	if speed < 0.1 and speed > 0.0 and self.step > 100:
    	    return True
    	else:
    	    return False

    def is_on_left_lane(self, vehicle_wp, left_wp, right_wp):
    	vehicle_wp = np.array((vehicle_wp.transform.location.x, vehicle_wp.transform.location.y, vehicle_wp.transform.location.z))
    	left_wp = np.array((left_wp.transform.location.x, left_wp.transform.location.y, left_wp.transform.location.z))
    	right_wp = np.array((right_wp.transform.location.x, right_wp.transform.location.y, right_wp.transform.location.z))
    	print("vehicle_wp {}".format(vehicle_wp))
    	print("left_wp {}".format(left_wp))
    	print("right_wp {}".format(right_wp))
    	dist_left = np.linalg.norm(left_wp - vehicle_wp)
    	dist_right = np.linalg.norm(right_wp - vehicle_wp)
    	print("right {} left {}".format(dist_right, dist_left))
    	if dist_right > dist_left:
    	    return True
    	else:
    	    return False

    def run_step(self, input_data, timestamp):

        # Lane invasion check
        if self.step == 0 or self.step == 1:
    	    self.prev_lane_intervention = 0
        if self.prev_lane_intervention < leaderboard.autoagents.agent_wrapper.invaded_lane:
            is_off_lane = True
            print("Lane intervention occured")
            intervened_lane = True
            self.prev_lane_intervention = leaderboard.autoagents.agent_wrapper.invaded_lane
        else:
            intervened_lane = False

        # Collision check
        #if self.step == 0 or self.step == 1:
    	#    self.prev_lane_intervention = 0
        #if self.prev_lane_intervention < leaderboard.autoagents.agent_wrapper.invaded_lane:
        #    is_off_lane = True
        #    print("Lane intervention occured")
        #    intervened_lane = True
        #    self.prev_lane_intervention = leaderboard.autoagents.agent_wrapper.invaded_lane
        #else:
        #    intervened_lane = False

        if not self.initialized:
            self._init()

        data = self.tick(input_data)
        gps = self._get_position(data)

        near_node, near_command = self._waypoint_planner.run_step(gps)
        far_node, far_command = self._command_planner.run_step(gps)

        # Autopilot steering
        steer_AP, throttle, brake, target_speed = self._get_control(near_node, far_node, data)

        # Waypoints of vehicle & reference (center of lane closest to vehicle)
        wp = self._world.get_map().get_waypoint(self._vehicle.get_location(), project_to_road=True) #lane_type=(carla.LaneType.Driving | carla.LaneType.Shoulder | carla.LaneType.Sidewalk))
        #wp_vehicle = self._world.get_map().get_waypoint(self._vehicle.get_location(), project_to_road=False)
        vehicle_pos = self._vehicle.get_location()
        #print("wp:",wp.transform.location.x, wp.transform.location.y, wp.transform.location.z)
        #print("wp vehicle:", wp_vehicle.transform.location.x, wp_vehicle.transform.location.y, wp_vehicle.transform.location.z)
        #print(vehicle_pos)
        #pdb.set_trace()

        # Convert reference wp to a numpy array
        ref_wp_array = np.array((wp.transform.location.x, wp.transform.location.y))

        # Convert vehicle position to numpy array
        if vehicle_pos is not None:
            vehicle_wp_array = np.array((vehicle_pos.x, vehicle_pos.y))
            width = wp.lane_width
            offset = np.linalg.norm(ref_wp_array - vehicle_wp_array)
            if (offset > width/2):
            	junction_offset = True
            else:
            	junction_offset = False
        else:
            junction_offset = False

        # Determine if autopilot must take over
        if time.time() - self.prev_intervention_time < self.takeovertime:
            self.autopilot = True
            steer, throttle, brake, target_speed = self._get_control(near_node, far_node, data)
		    #elif (self.has_collided(data['speed']) == True) or (self.prev_wp == wp):
		    #    self.collided = True
		    #    self.interventions += 1
		    #    self.autopilot = True
		    #    steer, throttle, brake, target_speed = self._get_control(near_node, far_node, data)
		    #    self.prev_intervention_time = time.time()
        elif (intervened_lane==True) or (vehicle_pos is None): #or (str(wp.lane_type) != 'Driving') or (junction_offset == True): #or (self.autopilot_next == True): #or (is_on_left==True):
            self.interventions += 1
            self.autopilot = True
            steer, throttle, brake, target_speed = self._get_control(near_node, far_node, data)
            self.prev_intervention_time = time.time()
        else:
            # retrieve rgb and lidar data
            rgb = torch.from_numpy(data['rgb_front'])
            rgb = rgb.view(1, rgb.size(2), rgb.size(0), rgb.size(1))
            lidar_point_cloud = data['lidar'][:,:3]
            lidar_transformed = torch.from_numpy(lidar_to_histogram_features(lidar_point_cloud, crop=256)).unsqueeze(0)
            # predict controls
            self.autopilot = False
            throttle = 0
            if data['speed'] < 20:
                throttle = 0.5
            #steer = self.net_controls(rgb.to('cuda',dtype=torch.float32), lidar_transformed.to('cuda',dtype=torch.float32)).item()
            steer = self.net_controls(rgb.to('cuda',dtype=torch.float32)).item()
        self.prev_wp = wp

        # Assign values to control
        control = carla.VehicleControl()

        #if self.reversed_time is not None:
        #    reverse_time_diff = self.step - self.reversed_time
            #print("reverse time diff {} reversed time {} step {}".format(reverse_time_diff, self.reversed_time, self.step))
        #else:
        #    reverse_time_diff = 100

        if self.autopilot == True:
            #if reverse_time_diff < 7:
            #    control.steer = 0.0
            #    control.reverse = True
            #    self.autopilot_next = True
            #    print("still reversing")
            if self.has_collided(data['speed']) == True:
                control.steer = 0.2
                control.reverse = True
                self.autopilot_next = True
                print("starting to reverse")
            else:
                steer = steer_AP
                control.steer = steer
                control.reverse = False
                self.autopilot_next = False
            # Store intervention by autopilot
            self.history['intervention'] += [1]
        else:
            control.steer = steer
            control.reverse = False
            self.autopilot_next = False
            # Store not intervened
            self.history['intervention'] += [0]
        control.throttle = throttle
        control.brake = float(brake)

        #print("step number", self.step)

        # Compute cte
        if vehicle_pos is None:
            self.cte += self.prev_cte # add previous cte when car goes out of bounds
        else:
            self.prev_cte = self.crossTrackErrorNew(vehicle_wp_array, ref_wp_array)
            self.cte += self.prev_cte # cross track error

        # Store current timestamp, cte, and interventions
        self.history['time'] += [time.time()-self.start_time]
        self.history['cte'] += [self.prev_cte]

        # Save current drive statistics
        if self.save_path is not None:
            self.episode_time_len = time.time() - self.start_time # total cumulated simulation runtime
            self.save(far_node, near_command, steer, throttle, brake, target_speed, data)

        return control

    def _should_brake(self):
        actors = self._world.get_actors()

        vehicle = self._is_vehicle_hazard(actors.filter('*vehicle*'))
        light = self._is_light_red(actors.filter('*traffic_light*'))
        walker = self._is_walker_hazard(actors.filter('*walker*'))
        stop_sign = self._is_stop_sign_hazard(actors.filter('*stop*'))

        return any(x is not None for x in [vehicle, light, walker, stop_sign])

    def _point_inside_boundingbox(self, point, bb_center, bb_extent):
        A = carla.Vector2D(bb_center.x - bb_extent.x, bb_center.y - bb_extent.y)
        B = carla.Vector2D(bb_center.x + bb_extent.x, bb_center.y - bb_extent.y)
        D = carla.Vector2D(bb_center.x - bb_extent.x, bb_center.y + bb_extent.y)
        M = carla.Vector2D(point.x, point.y)

        AB = B - A
        AD = D - A
        AM = M - A
        am_ab = AM.x * AB.x + AM.y * AB.y
        ab_ab = AB.x * AB.x + AB.y * AB.y
        am_ad = AM.x * AD.x + AM.y * AD.y
        ad_ad = AD.x * AD.x + AD.y * AD.y

        return am_ab > 0 and am_ab < ab_ab and am_ad > 0 and am_ad < ad_ad

    def _get_forward_speed(self, transform=None, velocity=None):
        """ Convert the vehicle transform directly to forward speed """
        if not velocity:
            velocity = self._vehicle.get_velocity()
        if not transform:
            transform = self._vehicle.get_transform()

        vel_np = np.array([velocity.x, velocity.y, velocity.z])
        pitch = np.deg2rad(transform.rotation.pitch)
        yaw = np.deg2rad(transform.rotation.yaw)
        orientation = np.array([np.cos(pitch) * np.cos(yaw), np.cos(pitch) * np.sin(yaw), np.sin(pitch)])
        speed = np.dot(vel_np, orientation)
        return speed

    def _is_actor_affected_by_stop(self, actor, stop, multi_step=20):
        """
        Check if the given actor is affected by the stop
        """
        affected = False
        # first we run a fast coarse test
        current_location = actor.get_location()
        stop_location = stop.get_transform().location
        if stop_location.distance(current_location) > self.PROXIMITY_THRESHOLD:
            return affected

        stop_t = stop.get_transform()
        transformed_tv = stop_t.transform(stop.trigger_volume.location)

        # slower and accurate test based on waypoint's horizon and geometric test
        list_locations = [current_location]
        waypoint = self._world.get_map().get_waypoint(current_location)
        for _ in range(multi_step):
            if waypoint:
                waypoint = waypoint.next(self.WAYPOINT_STEP)[0]
                if not waypoint:
                    break
                list_locations.append(waypoint.transform.location)

        for actor_location in list_locations:
            if self._point_inside_boundingbox(actor_location, transformed_tv, stop.trigger_volume.extent):
                affected = True

        return affected

    def _is_stop_sign_hazard(self, stop_sign_list):
        if self._affected_by_stop:
            if not self._stop_completed:
                current_speed = self._get_forward_speed()
                if current_speed < self.SPEED_THRESHOLD:
                    self._stop_completed = True
                    return None
                else:
                    return self._target_stop_sign
            else:
                # reset if the ego vehicle is outside the influence of the current stop sign
                if not self._is_actor_affected_by_stop(self._vehicle, self._target_stop_sign):
                    self._affected_by_stop = False
                    self._stop_completed = False
                    self._target_stop_sign = None
                return None

        ve_tra = self._vehicle.get_transform()
        ve_dir = ve_tra.get_forward_vector()

        wp = self._world.get_map().get_waypoint(ve_tra.location)
        wp_dir = wp.transform.get_forward_vector()

        dot_ve_wp = ve_dir.x * wp_dir.x + ve_dir.y * wp_dir.y + ve_dir.z * wp_dir.z

        if dot_ve_wp > 0:  # Ignore all when going in a wrong lane
            for stop_sign in stop_sign_list:
                if self._is_actor_affected_by_stop(self._vehicle, stop_sign):
                    # this stop sign is affecting the vehicle
                    self._affected_by_stop = True
                    self._target_stop_sign = stop_sign
                    return self._target_stop_sign

        return None

    def _is_light_red(self, lights_list):
        if self._vehicle.get_traffic_light_state() != carla.libcarla.TrafficLightState.Green:
            affecting = self._vehicle.get_traffic_light()

            for light in self._traffic_lights:
                if light.id == affecting.id:
                    light.set_state(carla.libcarla.TrafficLightState.Green)

        return None

    def _is_walker_hazard(self, walkers_list):
        z = self._vehicle.get_location().z
        p1 = _numpy(self._vehicle.get_location())
        v1 = 10.0 * _orientation(self._vehicle.get_transform().rotation.yaw)

        for walker in walkers_list:
            v2_hat = _orientation(walker.get_transform().rotation.yaw)
            s2 = np.linalg.norm(_numpy(walker.get_velocity()))

            if s2 < 0.05:
                v2_hat *= s2

            p2 = -3.0 * v2_hat + _numpy(walker.get_location())
            v2 = 8.0 * v2_hat

            collides, collision_point = get_collision(p1, v1, p2, v2)

            if collides:
                return walker

        return None

    def _is_vehicle_hazard(self, vehicle_list):
        z = self._vehicle.get_location().z

        o1 = _orientation(self._vehicle.get_transform().rotation.yaw)
        p1 = _numpy(self._vehicle.get_location())
        s1 = max(10, 3.0 * np.linalg.norm(_numpy(self._vehicle.get_velocity()))) # increases the threshold distance
        v1_hat = o1
        v1 = s1 * v1_hat

        for target_vehicle in vehicle_list:
            if target_vehicle.id == self._vehicle.id:
                continue

            o2 = _orientation(target_vehicle.get_transform().rotation.yaw)
            p2 = _numpy(target_vehicle.get_location())
            s2 = max(5.0, 2.0 * np.linalg.norm(_numpy(target_vehicle.get_velocity())))
            v2_hat = o2
            v2 = s2 * v2_hat

            p2_p1 = p2 - p1
            distance = np.linalg.norm(p2_p1)
            p2_p1_hat = p2_p1 / (distance + 1e-4)

            angle_to_car = np.degrees(np.arccos(v1_hat.dot(p2_p1_hat)))
            angle_between_heading = np.degrees(np.arccos(o1.dot(o2)))

            # to consider -ve angles too
            angle_to_car = min(angle_to_car, 360.0 - angle_to_car)
            angle_between_heading = min(angle_between_heading, 360.0 - angle_between_heading)

            if angle_between_heading > 60.0 and not (angle_to_car < 15 and distance < s1):
                continue
            elif angle_to_car > 30.0:
                continue
            elif distance > s1:
                continue

            return target_vehicle

        return None

    def save(self, far_node, near_command, steer, throttle, brake, target_speed, tick_data):
        frame = self.step #// 10

        pos = self._get_position(tick_data)
        theta = tick_data['compass']
        speed = tick_data['speed']

        data = {
                'x': pos[0],
                'y': pos[1],
                'theta': theta,
                'speed': speed,
                'target_speed': target_speed,
                'x_command': far_node[0],
                'y_command': far_node[1],
                'command': near_command.value,
                'steer': steer,
                'throttle': throttle,
                'brake': brake,
                'weather': self.weather_id,
                }

        for sensor in self.sensors():
            if 'camera' in sensor['type'] and 'map' not in sensor['id']:
                if sensor['id']=='rgb_front':
                    Image.fromarray(tick_data[sensor['id']]).save(self.save_path / sensor['id'] / ('%04d.png' % frame))
            elif 'lidar' in sensor['type']:
                np.save(self.save_path / 'lidar' / ('%04d.npy' % frame), tick_data['lidar'], allow_pickle=True)

        #Image.fromarray(tick_data['topdown']).save(self.save_path / 'topdown' / ('%04d.png' % frame))
        measurements_file = self.save_path / 'measurements' / ('%04d.json' % frame)
        f = open(measurements_file, 'w')
        json.dump(data, f, indent=4)
        f.close()

        # Save metrics (cte, autonomy, interventions)
        autonomy = (1 - (self.interventions * self.takeovertime) / self.episode_time_len) * 100
        results = "Number of Interventions: {}\n \
        	Autonomy: {}\n \
        	CTE: {}\n".format(self.interventions, autonomy, self.cte/self.step)
        result_filepath = self.save_path /'drive_results.txt'
        #print("Saving results to", result_filepath)
        with open(result_filepath,'w') as fp:
        	fp.write(results)
        fp.close()

        # Save cte & interventions vs time for plotting
        history_path = self.save_path / 'history.pkl'
        with open(history_path,'wb') as f:
        	pickle.dump(self.history,f)

max_stars_repo_name: yewon-lee/transfuser

max_issues_repo_name: yewon-lee/transfuser

max_forks_repo_name: yewon-lee/transfuser