text:
mutable struct Ocean

    id       :: Integer  # 0 = master, 1, ..., N = workers

    gi       :: Union{DisplacedPoleCoordinate.GridInfo, Nothing}
    gi_file  :: Union{AbstractString, Nothing}
    mi       :: Union{ModelMap.MapInfo, Nothing}

    Nx       :: Integer           # Number of columns in i direction
    Ny       :: Integer           # Number of columns in j direction
    Nz_bone  :: Integer           # Number of layers  in k direction

    zs_bone  :: AbstractArray{Float64, 1} # Unmasked zs bone
    topo     :: AbstractArray{Float64, 2} # Depth of the topography. Negative value if it is underwater
    zs       :: AbstractArray{Float64, 3} # Actuall zs coordinate masked by topo
    Nz       :: AbstractArray{Int64, 2} # Number of layers that is active

    K_v      :: Float64           # Diffusion coe of momentum. Used by ekman flow calculation
    Dh_T      :: Float64           # Horizontal diffusion coe of temperature
    Dv_T      :: Float64           # Vertical   diffusion coe of temperature
    Dh_S      :: Float64           # Horizontal diffusion coe of salinity
    Dv_S      :: Float64           # Vertical   diffusion coe of salinity

    fs       :: AbstractArray{Float64, 2}
    ϵs       :: AbstractArray{Float64, 2}

    mask3    :: AbstractArray{Float64, 3}
    nomotionmask3    :: AbstractArray{Float64, 3}

    mask     :: AbstractArray{Float64, 2}
    mask_idx  :: Any
    valid_idx :: AbstractArray{Int64, 2}

    b_ML     :: AbstractArray{Float64, 2}
    T_ML     :: AbstractArray{Float64, 2}
    S_ML     :: AbstractArray{Float64, 2}

    ΔT       :: AbstractArray{Float64, 2}
    ΔS       :: AbstractArray{Float64, 2}

    dΔTdt     :: AbstractArray{Float64, 2}
    dΔSdt     :: AbstractArray{Float64, 2}

    h_ML     :: AbstractArray{Float64, 2}
    h_MO     :: AbstractArray{Float64, 2}
    fric_u   :: AbstractArray{Float64, 2}
    dTdt_ent    :: AbstractArray{Float64, 2}
    dSdt_ent    :: AbstractArray{Float64, 2}

    # To calculate ocean heat transport, we need the
    # unconserved part of energy
    TSAS_clim   :: AbstractArray{Float64, 2}
    SSAS_clim   :: AbstractArray{Float64, 2}
    TFLUX_bot       :: AbstractArray{Float64, 2}
    SFLUX_bot       :: AbstractArray{Float64, 2}
    SFLUX_top       :: AbstractArray{Float64, 2}
    TFLUX_DIV_implied      :: AbstractArray{Float64, 2}
    SFLUX_DIV_implied      :: AbstractArray{Float64, 2}

    bs       :: AbstractArray{Float64, 3}

    Ts       :: AbstractArray{Float64, 3}
    Ss       :: AbstractArray{Float64, 3}

    # Ts and Ss with FLDO mixed so that it represents the temp and salt of the whole grid
    # Only master need this to output data
    Ts_mixed :: Union{Nothing, AbstractArray{Float64, 3}}
    Ss_mixed :: Union{Nothing, AbstractArray{Float64, 3}}

    FLDO           :: AbstractArray{Int64,   2}
    FLDO_ratio_top :: AbstractArray{Float64,   2}
    FLDO_ratio_bot :: AbstractArray{Float64,   2}

    qflx2atm         :: AbstractArray{Float64, 2} # The energy flux to atmosphere if freezes
    qflx2atm_pos     :: AbstractArray{Float64, 2} # The energy flux to atmosphere if freezes
    qflx2atm_neg     :: AbstractArray{Float64, 2} # The energy flux to atmosphere if freezes
    TEMP             :: AbstractArray{Float64, 2} # Total heat content
    dTEMPdt          :: AbstractArray{Float64, 2} # Total heat content change rate
    SALT             :: AbstractArray{Float64, 2} # Total salt
    dSALTdt          :: AbstractArray{Float64, 2} # Total salt change rate

    qflx_T_correction :: AbstractArray{Float64, 3}
    qflx_S_correction :: AbstractArray{Float64, 3}

    seaice_nudge_energy :: AbstractArray{Float64, 2}

    h_ML_min :: AbstractArray{Float64, 2}
    h_ML_max :: AbstractArray{Float64, 2}
    we_max   :: Float64

    # Advection related variables
    τx       :: AbstractArray{Float64, 2}
    τy       :: AbstractArray{Float64, 2}

    u        :: AbstractArray{Float64, 3}
    v        :: AbstractArray{Float64, 3}
    w        :: AbstractArray{Float64, 3}

    u_bnd    :: AbstractArray{Float64, 3}
    v_bnd    :: AbstractArray{Float64, 3}
    w_bnd    :: AbstractArray{Float64, 3}

    GRAD_bnd_x :: AbstractArray{Float64, 3}
    GRAD_bnd_y :: AbstractArray{Float64, 3}
    GRAD_bnd_z :: AbstractArray{Float64, 3}

    CURV_x     :: AbstractArray{Float64, 3}
    CURV_y     :: AbstractArray{Float64, 3}
    CURV_z     :: AbstractArray{Float64, 3}

    TFLUX_DEN_x :: AbstractArray{Float64, 3}
    TFLUX_DEN_y :: AbstractArray{Float64, 3}
    TFLUX_DEN_z :: AbstractArray{Float64, 3}

    SFLUX_DEN_x :: AbstractArray{Float64, 3}
    SFLUX_DEN_y :: AbstractArray{Float64, 3}
    SFLUX_DEN_z :: AbstractArray{Float64, 3}

    div           :: AbstractArray{Float64, 3}
    TFLUX_CONV    :: AbstractArray{Float64, 3}
    TFLUX_CONV_h  :: AbstractArray{Float64, 3}
    SFLUX_CONV    :: AbstractArray{Float64, 3}
    SFLUX_CONV_h  :: AbstractArray{Float64, 3}

    ∇∇T      :: AbstractArray{Float64, 3}
    ∇∇S      :: AbstractArray{Float64, 3}

    # Radiation Scheme
    # The parameterization is referenced to Paulson and Simpson (1977).
    # I made assumption that ζ1→0.0 as adapted in Oberhuber (1993).
    # This means `R` portion of the irradiance is going to be treated
    # as surface heat fluxes (δ like absorption).
    #
    # Reference:
    #
    # 1. Oberhuber, J. M. (1993). Simulation of the Atlantic circulation with a coupled sea
    #    ice-mixed layer-isopycnal general circulation model. Part I: Model description.
    #    Journal of Physical Oceanography, 23(5), 808-829.
    #
    # 2. Paulson, C. A., & Simpson, J. J. (1977). Irradiance measurements in the upper ocean.
    #    Journal of Physical Oceanography, 7(6), 952-956.
    #
    R               :: Float64   # Fast absorption portion of sunlight.
    ζ               :: Float64   # Light penetration depth of DO ( = ζ2 in Paulson and Simpson (1977) )

    rad_decay_coes  :: Union{AbstractArray{Float64, 3}, Nothing}
    rad_absorp_coes :: Union{AbstractArray{Float64, 3}, Nothing}

    # Climatology states
    Ts_clim_relax_time :: Union{Float64, Nothing}
    Ss_clim_relax_time :: Union{Float64, Nothing}

    Ts_clim  :: Union{AbstractArray{Float64, 3}, Nothing}
    Ss_clim  :: Union{AbstractArray{Float64, 3}, Nothing}

    # Derived quantities
    N_ocs  :: Integer           # Number of columns
    hs     :: AbstractArray{Float64, 3} # Thickness of layers
    Δzs    :: AbstractArray{Float64, 3} # Δz between layers

    in_flds :: InputFields

    lays :: Union{Nothing, NamedTuple}
    cols :: Union{Nothing, NamedTuple}

    # Accumulative variables
    acc_vars :: Union{AccumulativeVariables, Nothing}

    ASUM :: Union{AdvectionSpeedUpMatrix, Nothing}
    wksp :: Workspace

    deep_ocn_correction_start_layer :: Integer

    qflx_finding_cdm :: Union{Nothing, CyclicDataManager}

    diag :: Dict

    function Ocean(;
        id       :: Integer = 0,
        gridinfo_file :: AbstractString,
        sub_yrng :: Union{UnitRange, Nothing} = nothing,
        Nx       :: Integer,
        Ny       :: Integer,
        zs_bone  :: AbstractArray{Float64, 1},
        Ts       :: Union{AbstractArray{Float64, 3}, AbstractArray{Float64, 1}, Float64},
        Ss       :: Union{AbstractArray{Float64, 3}, AbstractArray{Float64, 1}, Float64},
        K_v      :: Float64 = 1e-2,
        Dh_T     :: Float64 = 1e3,
        Dv_T     :: Float64 = 1e-5,
        Dh_S     :: Float64 = 1e3,
        Dv_S     :: Float64 = 1e-5,
        T_ML     :: Union{AbstractArray{Float64, 2}, Float64},
        S_ML     :: Union{AbstractArray{Float64, 2}, Float64},
        h_ML     :: Union{AbstractArray{Float64, 2}, Float64, Nothing},
        h_ML_min :: Union{AbstractArray{Float64, 2}, Float64} = 30.0,
        h_ML_max :: Union{AbstractArray{Float64, 2}, Float64},
        we_max   :: Float64 =  1e-2,
        R        :: Float64 =  0.58,  # See Paulson and Simpson (1977) Type I clear water
        ζ        :: Float64 = 23.00,  # See Paulson and Simpson (1977) Type I clear water
        Ts_clim_relax_time :: Union{Float64, Nothing},
        Ss_clim_relax_time :: Union{Float64, Nothing},
        Ts_clim  :: Union{AbstractArray{Float64, 3}, AbstractArray{Float64, 1}, Nothing},
        Ss_clim  :: Union{AbstractArray{Float64, 3}, AbstractArray{Float64, 1}, Nothing},
        topo     :: Union{AbstractArray{Float64, 2}, Nothing},
        fs       :: Union{AbstractArray{Float64, 2}, Float64, Nothing} = nothing,
        ϵs       :: Union{AbstractArray{Float64, 2}, Float64, Nothing} = 1.0 / 86400.0,  # Bayesian estimation of 50m thick Ekman layer using CESM LENS POP output WVEL and τ
        in_flds  :: Union{InputFields, Nothing} = nothing,
        arrange  :: Symbol = :zxy,
        do_convective_adjustment :: Bool = false,
        qflx_finding_forcing_file :: String = "",
        verbose  :: Bool = false,
    )
        # Determine whether data should be local or shared (parallelization)
        datakind = ( id == 0 ) ? (:shared) : (:local)

        # ===== [BEG] GridInfo =====

        gridinfo = nothing
        mi = ModelMap.MapInfo{Float64}(gridinfo_file)

        _mask = allocate(datakind, Float64, Nx, Ny)

        # mask =>   lnd = 0, ocn = 1

        if id == 0

            gridinfo = DisplacedPoleCoordinate.GridInfo(
                Re,
                mi.nx,
                mi.ny,
                mi.xc,
                mi.yc,
                mi.xv,
                mi.yv,
                mi.area;
                angle_unit=:deg,
            )
            _mask .= mi.mask

        else
            if sub_yrng == nothing
                thorw(ErrorException("Init worker ocean, sub_yrng must be provided."))
            end

            gridinfo = DisplacedPoleCoordinate.GridInfo(
                Re,
                mi.nx,
                mi.ny,
                mi.xc,
                mi.yc,
                mi.xv,
                mi.yv,
                mi.area;
                angle_unit=:deg,
                sub_yrng=sub_yrng,
            )

            _mask .= view(mi.mask, :, sub_yrng)

        end

        mask_idx = (_mask .== 1.0)
        # ===== [END] GridInfo =====

        # ===== [BEGIN] topo, mask, h_ML_min, h_ML_max =====
        # Min/max of ML is tricky because it cannot be
        # deeper than the bottom boundary
        # Also, in real data topo can be 0 and not masked out

        _topo        = allocate(datakind, Float64, Nx, Ny)
        _h_ML_min    = allocate(datakind, Float64, Nx, Ny)
        _h_ML_max    = allocate(datakind, Float64, Nx, Ny)

        if topo == nothing
            _topo .= zs_bone[end]
        else
            _topo[:, :] = topo
        end

        # Arrage like (2, cnt) instead of (cnt, 2) to
        # enhance speed through memory cache
        valid_idx = allocate(datakind, Int64, 2, sum(mask_idx))

        let k = 1
            for idx in CartesianIndices((Nx, Ny))
                if _mask[idx] == 1.0
                    valid_idx[1, k] = idx[1]
                    valid_idx[2, k] = idx[2]

                    k += 1
                end
            end

            if k != size(valid_idx)[2] + 1
                throw(ErrorException("Initialization error making `valid_idx`"))
            end
        end

        #println("################ h_ML_min: ", h_ML_min)

        if typeof(h_ML_min) <: AbstractArray{Float64, 2}
            _h_ML_min[:, :] = h_ML_min
        elseif typeof(h_ML_min) <: Float64
            _h_ML_min .= h_ML_min
        end

        if typeof(h_ML_max) <: AbstractArray{Float64, 2}
            _h_ML_max[:, :] = h_ML_max
        elseif typeof(h_ML_max) <: Float64
            _h_ML_max .= h_ML_max
        end

        # Detect and fix h_ML_{max,min}
        coord_max = -zs_bone[end]
        for i=1:Nx, j=1:Ny

            if _mask[i, j] == 0
                _h_ML_max[i, j] = NaN
                _h_ML_min[i, j] = NaN
                continue
            end

            hmax = _h_ML_max[i, j]
            hmin = _h_ML_min[i, j]
            hbot = - _topo[i, j]

            if hbot < 0
                throw(ErrorException(format("Topography is negative at idx ({:d}, {:d})", i, j)))
            elseif hbot == 0
                #throw(ErrorException(format("Topography is zero at idx ({:d}, {:d})", i, j)))
                println(format("Topography is zero at idx ({:d}, {:d})", i, j))
            end

            if hmax < hmin
                throw(ErrorException(format("h_ML_max must ≥ h_ML_min. Problem happens at idx ({:d}, {:d})", i, j)))
            end

            if coord_max < hmin
                verbose && println(format("Point ({},{}) got h_min {:.2f} which is larger than coord_max {}. Tune h_ML_min to match it.", i, j, hmin, coord_max))
                hmin = coord_max
            end

            if coord_max < hmax
                verbose && println(format("Point ({},{}) got h_max {:.2f} which is larger than coord_max {}. Tune h_ML_max to match it.", i, j, hmax, coord_max))
                hmax = coord_max
            end

            if hmin > hbot
                verbose && println(format("Point ({},{}) got depth {:.2f} which is smaller than h_ML_min {}. Tune h_ML_min/max to depth.", i, j, hbot, hmin))
                hbot = hmin
            end

            if hmax > hbot
                verbose && println(format("Point ({},{}) got depth {:.2f} which is smaller than h_ML_max {}. Tune the h_ML_max to depth.", i, j, hbot, hmax))
                hmax = hbot
            end

            _h_ML_min[i, j] = hmin
            _h_ML_max[i, j] = hmax
            _topo[i, j]     = -hbot

        end

        # ===== [END] topo, mask, h_ML_min, h_ML_max =====

        # ===== [BEGIN] z coordinate =====
        zs_bone = copy(zs_bone)
        Nz_bone = length(zs_bone) - 1

        Nz   = allocate(datakind, Int64, Nx, Ny)
        zs   = allocate(datakind, Float64, Nz_bone + 1, Nx, Ny)
        hs   = allocate(datakind, Float64, Nz_bone    , Nx, Ny)
        Δzs  = allocate(datakind, Float64, Nz_bone + 1, Nx, Ny)

        zs  .= 1.0
        Nz  .= 0
        hs  .= 1.0
        Δzs .= 1.0

        for i=1:Nx, j=1:Ny

            if _mask[i, j] == 0
                continue
            end

            # Determine Nz

            # Default is that topo is deeper than
            # the bottom of zs_bone
            _Nz = Nz_bone
            for k=2:length(zs_bone)
                if zs_bone[k] <= _topo[i, j]
                    _Nz = k-1
                    #println(format("This topo gets: zs_bone[{:d}] = {:f}, _topo[{:d},{:d}]={:f}", k, zs_bone[k], i, j, _topo[i,j]))
                    break
                end
            end

            Nz[i, j] = _Nz

            # Construct vertical coordinate
            zs[1:_Nz, i, j] = zs_bone[1:_Nz]

            zs[_Nz+1, i, j] = max(_topo[i, j], zs_bone[_Nz+1])

            # Construct thickness of each layer
            hs[ 1:_Nz,   i, j] = zs[1:_Nz, i, j] - zs[2:_Nz+1, i, j]
            Δzs[2:_Nz, i, j] = (hs[1:_Nz-1, i, j] + hs[2:_Nz, i, j]) / 2.0
            Δzs[1, i, j] = Δzs[2, i, j]
            Δzs[end, i, j] = Δzs[end-1, i, j]

        end

        # ===== [END] z coordinate =====

        # ===== [BEGIN] Column information =====

        _b_ML     = allocate(datakind, Float64, Nx, Ny)
        _T_ML     = allocate(datakind, Float64, Nx, Ny)
        _S_ML     = allocate(datakind, Float64, Nx, Ny)
        _ΔT       = allocate(datakind, Float64, Nx, Ny)
        _ΔS       = allocate(datakind, Float64, Nx, Ny)
        _dΔTdt    = allocate(datakind, Float64, Nx, Ny)
        _dΔSdt    = allocate(datakind, Float64, Nx, Ny)
        _h_ML     = allocate(datakind, Float64, Nx, Ny)
        _h_MO     = allocate(datakind, Float64, Nx, Ny)
        _fric_u   = allocate(datakind, Float64, Nx, Ny)
        _dTdt_ent    = allocate(datakind, Float64, Nx, Ny)
        _dSdt_ent    = allocate(datakind, Float64, Nx, Ny)

        _TSAS_clim   = allocate(datakind, Float64, Nx, Ny)
        _SSAS_clim   = allocate(datakind, Float64, Nx, Ny)
        _TFLUX_bot       = allocate(datakind, Float64, Nx, Ny)
        _SFLUX_bot       = allocate(datakind, Float64, Nx, Ny)
        _SFLUX_top       = allocate(datakind, Float64, Nx, Ny)
        _TFLUX_DIV_implied      = allocate(datakind, Float64, Nx, Ny)
        _SFLUX_DIV_implied      = allocate(datakind, Float64, Nx, Ny)

        _bs       = allocate(datakind, Float64, Nz_bone, Nx, Ny)
        _Ts       = allocate(datakind, Float64, Nz_bone, Nx, Ny)
        _Ss       = allocate(datakind, Float64, Nz_bone, Nx, Ny)
        _FLDO     = allocate(datakind, Int64, Nx, Ny)
        _FLDO_ratio_top = allocate(datakind, Float64, Nx, Ny)
        _FLDO_ratio_bot = allocate(datakind, Float64, Nx, Ny)
        _qflx2atm  = allocate(datakind, Float64, Nx, Ny)
        _qflx2atm_pos = allocate(datakind, Float64, Nx, Ny)
        _qflx2atm_neg = allocate(datakind, Float64, Nx, Ny)
        _TEMP         = allocate(datakind, Float64, Nx, Ny)
        _dTEMPdt      = allocate(datakind, Float64, Nx, Ny)
        _SALT      = allocate(datakind, Float64, Nx, Ny)
        _dSALTdt   = allocate(datakind, Float64, Nx, Ny)
        _qflx_T_correction = allocate(datakind, Float64, Nz_bone, Nx, Ny)
        _qflx_S_correction = allocate(datakind, Float64, Nz_bone, Nx, Ny)
        _seaice_nudge_energy = allocate(datakind, Float64, Nx, Ny)

        if typeof(h_ML) <: AbstractArray{Float64, 2}
            _h_ML[:, :] = h_ML
        elseif typeof(h_ML) <: Float64
            _h_ML .= h_ML
        elseif h_ML == nothing
            _h_ML .= h_ML_min
        end

        # Need to constraint h_ML
        for i=1:Nx, j=1:Ny

            (_mask[i, j] == 0.0) && continue

            _h_ML[i, j] = boundMLD(_h_ML[i, j]; h_ML_max=_h_ML_max[i, j], h_ML_min=_h_ML_min[i, j])

        end

        if typeof(T_ML) <: AbstractArray{Float64, 2}
            _T_ML[:, :] = T_ML
        elseif typeof(T_ML) <: Float64
            _T_ML .= T_ML
        end

        if typeof(S_ML) <: AbstractArray{Float64, 2}
            _S_ML[:, :] = S_ML
        elseif typeof(S_ML) <: Float64
            _S_ML .= S_ML
        end

        if typeof(Ts) <: AbstractArray{Float64, 3}
            _Ts[:, :, :] = toZXY(Ts, arrange)
        elseif typeof(Ts) <: AbstractArray{Float64, 1}
            for i=1:Nx, j=1:Ny
                _Ts[:, i, j] = Ts
            end
        elseif typeof(Ts) <: Float64
            _Ts .= Ts
        end

        if typeof(Ss) <: AbstractArray{Float64, 3}
            _Ss[:, :, :] = toZXY(Ss, arrange)
        elseif typeof(Ss) <: AbstractArray{Float64, 1}
            for i=1:Nx, j=1:Ny
                _Ss[:, i, j] = Ss
            end
        elseif typeof(Ss) <: Float64
            _Ss .= Ss
        end

        _Ts_mixed = copy(_Ts)
        _Ss_mixed = copy(_Ts)

        # ===== [END] Column information =====

        # ===== [BEGIN] Radiation =====

        _rad_decay_coes  = allocate(datakind, Float64, Nz_bone, Nx, Ny)
        _rad_absorp_coes = allocate(datakind, Float64, Nz_bone, Nx, Ny)

        for i=1:Nx, j=1:Ny

            if _mask[i, j] == 0.0
                continue
            end

            for k=1:Nz[i, j]
                _rad_decay_coes[k, i, j]  = exp(zs[k, i, j] / ζ)         # From surface to top of the layer
                _rad_absorp_coes[k, i, j] = 1.0 - exp(- hs[k, i, j] / ζ)
            end

            # Since we assume the bottome of ocean absorbs anything
            _rad_absorp_coes[Nz[i, j], i, j] = 1.0
        end

        # ===== [END] Radiation =====

        # ===== [BEGIN] fs and ϵs =====

        _fs       = allocate(datakind, Float64, Nx, Ny)
        _ϵs       = allocate(datakind, Float64, Nx, Ny)

        if typeof(fs) <: AbstractArray{Float64, 2}
            _fs[:, :] = fs
        elseif typeof(fs) <: Float64
            _fs .= fs
        elseif fs == nothing
           _fs[:, :] = 2 * Ωe * sin.(gridinfo.c_lat)
        end

        if typeof(ϵs) <: AbstractArray{Float64, 2}
            _ϵs[:, :] = ϵs
        elseif typeof(ϵs) <: Float64
            _ϵs .= ϵs
        end

        # ===== [END] fs and ϵs =====

        # ===== [BEGIN] Advection variables =====

        _τx      = allocate(datakind, Float64, Nx, Ny)
        _τy      = allocate(datakind, Float64, Nx, Ny)
        _u       = allocate(datakind, Float64, Nz_bone, Nx, Ny)
        _v       = allocate(datakind, Float64, Nz_bone, Nx, Ny)
        _w       = allocate(datakind, Float64, Nz_bone, Nx, Ny)

        _u_bnd   = allocate(datakind, Float64, Nz_bone, Nx, Ny)
        _v_bnd   = allocate(datakind, Float64, Nz_bone, Nx, Ny+1)
        _w_bnd   = allocate(datakind, Float64, Nz_bone+1, Nx, Ny)

        _GRAD_bnd_x   = allocate(datakind, Float64, Nz_bone, Nx, Ny)
        _GRAD_bnd_y   = allocate(datakind, Float64, Nz_bone, Nx, Ny+1)
        _GRAD_bnd_z   = allocate(datakind, Float64, Nz_bone+1, Nx, Ny)

        _CURV_x       = allocate(datakind, Float64, Nz_bone, Nx, Ny)
        _CURV_y       = allocate(datakind, Float64, Nz_bone, Nx, Ny)
        _CURV_z       = allocate(datakind, Float64, Nz_bone, Nx, Ny)

        _TFLUX_DEN_x  = allocate(datakind, Float64, Nz_bone, Nx, Ny)
        _TFLUX_DEN_y  = allocate(datakind, Float64, Nz_bone, Nx, Ny+1)
        _TFLUX_DEN_z  = allocate(datakind, Float64, Nz_bone+1, Nx, Ny)

        _SFLUX_DEN_x  = allocate(datakind, Float64, Nz_bone, Nx, Ny)
        _SFLUX_DEN_y  = allocate(datakind, Float64, Nz_bone, Nx, Ny+1)
        _SFLUX_DEN_z  = allocate(datakind, Float64, Nz_bone+1, Nx, Ny)

        _div     = allocate(datakind, Float64, Nz_bone, Nx, Ny)

        _TFLUX_CONV   = allocate(datakind, Float64, Nz_bone, Nx, Ny)
        _TFLUX_CONV_h = allocate(datakind, Float64, Nz_bone, Nx, Ny)
        _SFLUX_CONV   = allocate(datakind, Float64, Nz_bone, Nx, Ny)
        _SFLUX_CONV_h = allocate(datakind, Float64, Nz_bone, Nx, Ny)

        _∇∇T     = allocate(datakind, Float64, Nz_bone, Nx, Ny)
        _∇∇S     = allocate(datakind, Float64, Nz_bone, Nx, Ny)

        # ===== [END] Advection variables =====

        # ===== [BEGIN] Climatology =====

        _Ts_clim = allocate(datakind, Float64, Nz_bone, Nx, Ny)
        if Ts_clim == nothing

        else

            if typeof(Ts_clim) <: AbstractArray{Float64, 3}

                _Ts_clim[:, :, :] = toZXY(Ts_clim, arrange)

            elseif typeof(Ts_clim) <: AbstractArray{Float64, 1}

                for i=1:Nx, j=1:Ny
                    _Ts_clim[:, i, j] = Ts_clim
                end

            end

        end

        _Ss_clim = allocate(datakind, Float64, Nz_bone, Nx, Ny)
        if Ss_clim == nothing

        else

            if typeof(Ss_clim) <: AbstractArray{Float64, 3}

                _Ss_clim[:, :, :] = toZXY(Ss_clim, arrange)

            elseif typeof(Ss_clim) <: AbstractArray{Float64, 1}

                for i=1:Nx, j=1:Ny
                    _Ss_clim[:, i, j] = Ss_clim
                end

            end

        end

        # ===== [END] Climatology =====

        # ===== [BEGIN] Mask out data =====

        _mask3        = allocate(datakind, Float64, Nz_bone, Nx, Ny)

        _mask3 .= 1.0
        # Clean up all variables
        for i=1:Nx, j=1:Ny
            _mask3[Nz[i, j] + 1:end, i, j] .= 0.0
        end

        #println("sum of _mask3: ", sum(_mask3))

        mask3_lnd_idx = (_mask3  .== 0.0)
        mask2_lnd_idx = (_mask  .== 0.0)

        for v in [_bs, _Ts, _Ss, _Ts_clim, _Ss_clim]
            if v == nothing
                continue
            end

            v[mask3_lnd_idx] .= -999
        end

        for v in [_b_ML, _T_ML, _S_ML, _h_ML, _h_ML_min, _h_ML_max]
            v[mask2_lnd_idx] .= -999
        end

        # ===== [END] Mask out data

        # ==== [BEGIN] Determine nomotionmask3 ====
        _nomotionmask3 = copy(_mask3)
        for i=1:Nx, j=1:Ny

            if _topo[i, j] > -300.0
                _nomotionmask3[:, i, j] .= 0.0
                continue
            end

            if _mask[i, j] == 1.0
                _nomotionmask3[ Nz[i, j], i, j ] = 0.0
            end

        end
        # ==== [END] Determine nomotionmask3 ====

        # ===== [BEGIN] check integrity =====

        # Check topography, h_ML_min/max and zs
        for i=1:Nx, j=1:Ny
            if _mask[i, j] == 0
                continue
            end

            if ! (-_h_ML_min[i, j] >= - _h_ML_max[i, j] >= zs[Nz[i, j] + 1, i, j] >= _topo[i, j])
                println("idx: (", i, ", ", j, ")")
                println("h_ML_min: ", _h_ML_min[i, j])
                println("h_ML_max: ", _h_ML_max[i, j])
                println("z_deepest: ", zs[Nz[i, j] + 1, i, j])
                println("topo: ", _topo[i, j])
                ErrorException("Relative relation is wrong") |> throw
            end
        end

        # Check if there is any hole in climatology

        mask3_idx = (_mask3 .== 1)
        valid_grids = sum(mask3_idx)
        total_data  = Nx * Ny * Nz_bone

        #println("Total  data count: ", total_data)
        #println("Valid  data count: ", valid_grids)
        #println("Masked data count: ", total_data - valid_grids)

        #=
        if sum(isfinite.(_Ss)) != valid_grids
            throw(ErrorException("Salinity data has holes"))
        end

        if sum(isfinite.(_Ts)) != valid_grids
            throw(ErrorException("Temperature data has holes"))
        end

        if _Ts_clim != nothing && sum(isfinite.(_Ts_clim)) != valid_grids
            throw(ErrorException("Temperature climatology has holes"))
        end

        if _Ss_clim != nothing && sum(isfinite.(_Ss_clim)) != valid_grids
            throw(ErrorException("Salinity climatology has holes"))
        end
        =#

        # Check if h_ML_min h_ML_max is negative
        if any(_h_ML_min[mask_idx] .<= 0)
            throw(ErrorException("h_ML_min should always be positive (cannot be zero or negative)"))
        end

        if any(_h_ML_max[mask_idx] .< 0)
            throw(ErrorException("h_ML_max should always be non-negative"))
        end

        # ===== [END] check integrity =====

        cols = nothing
        lays = nothing

#        if id != 0

            # ===== [BEGIN] Construct Views of Lays =====
            lays = ((
                hs      = Array{SubArray}(undef, Nz_bone),
                bs      = Array{SubArray}(undef, Nz_bone),
                Ts      = Array{SubArray}(undef, Nz_bone),
                Ss      = Array{SubArray}(undef, Nz_bone),
                u       = Array{SubArray}(undef, Nz_bone),
                v       = Array{SubArray}(undef, Nz_bone),
                div     = Array{SubArray}(undef, Nz_bone),
                ∇∇T     = Array{SubArray}(undef, Nz_bone),
                ∇∇S     = Array{SubArray}(undef, Nz_bone),
                mask3   = Array{SubArray}(undef, Nz_bone),
            ))

            for k=1:Nz_bone
                lays.hs[k]      = view(hs, k, :, :)
                lays.bs[k]      = view(_bs, k, :, :)
                lays.Ts[k]      = view(_Ts, k, :, :)
                lays.Ss[k]      = view(_Ss, k, :, :)
                lays.u[k]       = view(_u, k, :, :)
                lays.v[k]       = view(_v, k, :, :)
                lays.div[k]     = view(_div, k, :, :)
                lays.∇∇T[k]     = view(_∇∇T, k, :, :)
                lays.∇∇S[k]     = view(_∇∇S, k, :, :)
                lays.mask3[k]   = view(_mask3, k, :, :)
            end

            # ===== [END] Construct Views of Lays =====

            # ===== [BEGIN] Construct Views of Cols =====
            cols = ((
                zs  = Array{SubArray}(undef, Nx, Ny),
                Δzs = Array{SubArray}(undef, Nx, Ny),
                hs = Array{SubArray}(undef, Nx, Ny),
                w  = Array{SubArray}(undef, Nx, Ny),
                w_bnd = Array{SubArray}(undef, Nx, Ny),
                bs = Array{SubArray}(undef, Nx, Ny),
                Ts = Array{SubArray}(undef, Nx, Ny),
                Ss = Array{SubArray}(undef, Nx, Ny),
                Ts_mixed = Array{SubArray}(undef, Nx, Ny),
                Ss_mixed = Array{SubArray}(undef, Nx, Ny),
                rad_decay_coes  = Array{SubArray}(undef, Nx, Ny),
                rad_absorp_coes = Array{SubArray}(undef, Nx, Ny),
                Ts_clim = (Ts_clim == nothing) ? nothing : Array{SubArray}(undef, Nx, Ny),
                Ss_clim = (Ss_clim == nothing) ? nothing : Array{SubArray}(undef, Nx, Ny),
                qflx_T_correction = Array{SubArray}(undef, Nx, Ny),
                qflx_S_correction = Array{SubArray}(undef, Nx, Ny),
            ))

            for i=1:Nx, j=1:Ny
                cols.zs[i, j]              = view(zs,  :, i, j)
                cols.Δzs[i, j]             = view(Δzs, :, i, j)
                cols.hs[i, j]              = view(hs,  :, i, j)
                cols.w[i, j]               = view(_w,  :, i, j)
                cols.w_bnd[i, j]           = view(_w_bnd,  :, i, j)
                cols.bs[i, j]              = view(_bs, :, i, j)
                cols.Ts[i, j]              = view(_Ts, :, i, j)
                cols.Ss[i, j]              = view(_Ss, :, i, j)
                cols.Ts_mixed[i, j]        = view(_Ts_mixed, :, i, j)
                cols.Ss_mixed[i, j]        = view(_Ss_mixed, :, i, j)
                cols.rad_decay_coes[i, j]  = view(_rad_decay_coes,  :, i, j)
                cols.rad_absorp_coes[i, j] = view(_rad_absorp_coes, :, i, j)
                cols.qflx_T_correction[i, j] = view(_qflx_T_correction, :, i, j)
                cols.qflx_S_correction[i, j] = view(_qflx_S_correction, :, i, j)
            end

            if Ts_clim != nothing
                for i=1:Nx, j=1:Ny
                    cols.Ts_clim[i, j] = view(_Ts_clim, :, i, j)
                end
            end

            if Ss_clim != nothing
                for i=1:Nx, j=1:Ny
                    cols.Ss_clim[i, j] = view(_Ss_clim, :, i, j)
                end
            end

#        end
        # ===== [END] Construct Views of Cols =====

        # ===== [BEGIN] Making speed-up matrix

        if true || id != 0
            ASUM = AdvectionSpeedUpMatrix(;
                gi    = gridinfo,
                Nz    = Nz_bone,
                Nz_av = Nz,
                mask3 = _mask3,
                Δz_W  = Δzs,
                Δz_T  = hs,
                nomotionmask3 = _nomotionmask3,
            )
        else
            ASUM = nothing
        end
        # ===== [END] Making speed-up matrix

        deep_ocn_correction_start_layer = getLayerFromDepth(;
            zs = zs_bone,
            z  = -1.0,
            Nz = length(zs_bone) - 1
        )

        diag = Dict(
            :total_heat => zeros(Float64, 1),
            :total_salt => zeros(Float64, 1),
            :total_heat_budget_residue => zeros(Float64, 1),
            :total_salt_budget_residue => zeros(Float64, 1),
        )

        # ===== [BEGIN] Q-flux finding forcing file =====
        qflx_finding_cdm = nothing
        if id == 0

            # do nothing

        else
            if sub_yrng == nothing
                thorw(ErrorException("Init worker ocean, sub_yrng must be provided."))
            end

            if qflx_finding_forcing_file != ""

                println("Load qflx file: ", qflx_finding_forcing_file)
                qflx_finding_cdm = CyclicDataManager(;
                    filename     = qflx_finding_forcing_file,
                    varname_time = "time",
                    varnames     = ["TEMP", "SALT"],
                    beg_time     = 0.0,
                    cyc_time     = 365.0,
                    spatial_rng  = (:, sub_yrng, :),  # in the shape of the file
                    xyz2zxy      = true,
                )
            end

        end

        # ===== [END] Q-flux finding forcing file =====

        ocn = new(
            id,
            gridinfo,
            gridinfo_file,
            mi,
            Nx, Ny, Nz_bone,
            zs_bone, _topo, zs, Nz,
            K_v, Dh_T, Dv_T, Dh_S, Dv_S,
            _fs, _ϵs,
            _mask3, _nomotionmask3,
            _mask, mask_idx, valid_idx,
            _b_ML, _T_ML, _S_ML, _ΔT, _ΔS, _dΔTdt, _dΔSdt,
            _h_ML, _h_MO, _fric_u, _dTdt_ent, _dSdt_ent,
            _TSAS_clim, _SSAS_clim,
            _TFLUX_bot, _SFLUX_bot, _SFLUX_top,
            _TFLUX_DIV_implied, _SFLUX_DIV_implied,
            _bs,   _Ts,   _Ss,
            _Ts_mixed,   _Ss_mixed,
            _FLDO, _FLDO_ratio_top, _FLDO_ratio_bot,
            _qflx2atm, _qflx2atm_pos, _qflx2atm_neg,
            _TEMP, _dTEMPdt, _SALT, _dSALTdt,
            _qflx_T_correction, _qflx_S_correction,
            _seaice_nudge_energy,
            _h_ML_min, _h_ML_max, we_max,
            _τx, _τy,
            _u, _v, _w,
            _u_bnd, _v_bnd, _w_bnd,
            _GRAD_bnd_x, _GRAD_bnd_y, _GRAD_bnd_z,
            _CURV_x, _CURV_y, _CURV_z,
            _TFLUX_DEN_x, _TFLUX_DEN_y, _TFLUX_DEN_z,
            _SFLUX_DEN_x, _SFLUX_DEN_y, _SFLUX_DEN_z,
            _div,
            _TFLUX_CONV, _TFLUX_CONV_h,
            _SFLUX_CONV, _SFLUX_CONV_h,
            _∇∇T, _∇∇S,
            R, ζ,
            _rad_decay_coes, _rad_absorp_coes,
            Ts_clim_relax_time, Ss_clim_relax_time,
            _Ts_clim, _Ss_clim,
            Nx * Ny, hs, Δzs,
            ( in_flds == nothing ) ? InputFields(datakind, Nx, Ny) : in_flds,
            lays,
            cols,
            ( id == 0 ) ? nothing : AccumulativeVariables(Nx, Ny, Nz_bone),
            ASUM,
            Workspace(Nx=Nx, Ny=Ny, Nz=Nz_bone, shape=:zxy),
            deep_ocn_correction_start_layer,
            qflx_finding_cdm,
            diag,
        )

        updateB!(ocn)
        updateFLDO!(ocn)

        if do_convective_adjustment
            @loop_hor ocn i j let
                OC_doConvectiveAdjustment!(ocn, i, j)
            end
        end

        calTEMP!(ocn)

        return ocn
    end

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
println("Calculating summary statistics")
for m in [boon,lti,volterra,kmm]
    err = StatsBase.msd(TDM.validate(quantity(Q2),1:N2,m.M),TDM.evalmodel(m,p2,Q2,1:N2))
    ns = nash_sutcliffe(m,p2,Q2)
    sf = flatness(m,p2,Q2)
    println(typeof(m),":: ",err,"\t",ns,"\t",sf)
end

ext: jl


text:
__precompile__(true)
""""
Enhanced string literals

String literals with Swift-like format, extendable at run-time

Copyright 2016-2020 Gandalf Software, Inc., Scott P. Jones
Licensed under MIT License, see LICENSE.md
"""
module StrLiterals

using ModuleInterfaceTools

const is_empty = isempty
const is_valid = isvalid
const is_printable = isprint
const TypeOrFunc = Union{DataType,Function}

outhex(v, p=1)   = string(v, base=16, pad=p)
_sprint(f, s)    = sprint(f, s; sizehint=lastindex(s))
_sprint(f, s, c) = sprint(f, s, c; sizehint=lastindex(s))

@api develop! interpolated_parse, interpolated_parse_vec, s_parse_unicode, s_parse_legacy,
              s_print_unescaped_legacy, s_print_unescaped, s_print_escaped, s_print,
              s_escape_string, s_unescape_string, s_unescape_str, s_unescape_legacy

@api public "@f_str", "@pr_str", "@F_str", "@PR_str", "@sym_str"
export @f_str, @pr_str, @F_str, @PR_str, @sym_str

const parse_chr   = Dict{Char, Function}()
const interpolate = Dict{Char, Function}()
const string_type = Ref{TypeOrFunc}(String)

const SymStr = Union{Symbol, AbstractString}

@api develop throw_arg_err, hexerr, parse_error, check_expr, check_done

str_done(str::AbstractString, pos::Integer) = (pos > ncodeunits(str))

parse_error(s) = throw(Base.Meta.ParseError(s))
incomplete_expr_error() = parse_error("Incomplete expression")
check_expr(ex) = isa(ex, Expr) && (ex.head === :continue) && incomplete_expr_error()
check_done(str, pos, msg) = str_done(str, pos) && parse_error(string(msg, " in ", repr(str)))

"""
Create a symbol from a string (that allows for interpolation and escape sequences)
"""
macro sym_str(str) ; QuoteNode(interpolated_parse(str, Symbol)) ; end

"""
String macro with more Swift-like syntax, plus support for emojis and LaTeX names
"""
macro f_str(str) ; interpolated_parse(str, string_type[]) ; end
macro f_str(str, args...)
    for v in args ; dump(v) end
    interpolated_parse(str, string_type[])
end

"""
String macro with more Swift-like syntax, plus support for emojis and LaTeX names, also legacy
"""
macro F_str(str) ; interpolated_parse(str, string_type[], true) ; end

"""
String macros that calls print directly
"""
macro pr_str(str) ; s_print(str, false) ; end
macro PR_str(str) ; s_print(str, true) ; end

throw_arg_err(msg)      = parse_error(msg)
throw_arg_err(msg, val) = parse_error(string(msg, repr(val)))

"""
Handle Unicode character constant, of form \\u{<hexdigits>}
"""
function s_parse_unicode(io, str,  pos)
    check_done(str, pos, "Incomplete \\u{...}")
    chr, pos = iterate(str, pos)
    chr != '{' && throw_arg_err("\\u missing opening { in ", str)
    check_done(str, pos, "Incomplete \\u{...}")
    beg = pos
    chr, pos = iterate(str, pos)
    num::UInt32 = 0
    cnt = 0
    while chr != '}'
        check_done(str, pos, "\\u{ missing closing }")
        (cnt += 1) > 6 && throw_arg_err("Unicode constant too long in ", str)
        num = num<<4 + chr - ('0' <= chr <= '9' ? '0' :
                              'a' <= chr <= 'f' ? 'a' - 10 :
                              'A' <= chr <= 'F' ? 'A' - 10 :
                              throw_arg_err("\\u missing closing } in ", str))
        chr, pos = iterate(str, pos)
    end
    cnt == 0 && throw_arg_err("\\u{} has no hex digits in ", str)
    ((0x0d800 <= num <= 0x0dfff) || num > 0x10ffff) &&
        throw_arg_err("Invalid Unicode character constant ", str[beg-3:pos-1])
    print(io, Char(num))
    pos
end

"""
String interpolation parsing, allow legacy \$, \\xHH, \\uHHHH, \\UHHHHHHHH
"""
s_print_unescaped_legacy(io, str::AbstractString) = s_print_unescaped(io, str, true)

"""
String interpolation parsing
Based on code resurrected from Julia base:
https://github.com/JuliaLang/julia/blob/deab8eabd7089e2699a8f3a9598177b62cbb1733/base/string.jl
"""
function s_print_unescaped(io, str::AbstractString, flg::Bool=false)
    pos = 1
    while !str_done(str, pos)
        chr, pos = iterate(str, pos)
        if !str_done(str, pos) && chr == '\\'
            chr, pos = iterate(str, pos)
            if (chr == 'u' ||  chr == 'U' || chr == 'x')
                if flg
                    pos = s_parse_legacy(io, str, pos, chr)
                elseif chr == 'u'
                    pos = s_parse_unicode(io, str, pos)
                else
                    throw_arg_err(string("\\", chr, " only supported in legacy mode (i.e. ",
                                         "F\"...\" or PR\"...\""))
                end
            elseif haskey(parse_chr, chr)
                pos = parse_chr[chr](io, str, pos, chr)
            else
                chr = (chr == '0' ? '\0' :
                       chr == '$' ? '$'  :
                       chr == '"' ? '"'  :
                       chr == '\'' ? '\'' :
                       chr == '\\' ? '\\' :
                       chr == 'a' ? '\a' :
                       chr == 'b' ? '\b' :
                       chr == 't' ? '\t' :
                       chr == 'n' ? '\n' :
                       chr == 'v' ? '\v' :
                       chr == 'f' ? '\f' :
                       chr == 'r' ? '\r' :
                       chr == 'e' ? '\e' :
                       throw_arg_err(string("Invalid \\", chr, " sequence in "), str))
                write(io, UInt8(chr))
            end
        else
            print(io, chr)
        end
    end
end

hexerr(chr) = throw_arg_err("\\$chr used with no following hex digits")

function s_parse_legacy(io, str, pos, chr)
    str_done(str, pos) && hexerr(chr)
    beg = pos
    max = chr == 'x' ? 2 : chr == 'u' ? 4 : 8
    if str[pos] == '{'
        max == 4 || throw_arg_err("{ only allowed with \\u")
        return s_parse_unicode(io, str, pos)
    end
    num = cnt = 0
    while (cnt += 1) <= max && !str_done(str, pos)
        chr, nxt = iterate(str, pos)
        num = '0' <= chr <= '9' ? num << 4 + chr - '0' :
              'a' <= chr <= 'f' ? num << 4 + chr - 'a' + 10 :
              'A' <= chr <= 'F' ? num << 4 + chr - 'A' + 10 : break
        pos = nxt
    end
    cnt == 1 && hexerr(chr)
    if max == 2
        write(io, UInt8(num))
    elseif is_valid(Char, num)
        print(io, Char(num))
    else
        throw_arg_err("Invalid Unicode character constant ", str[beg-2:pos-1])
    end
    pos
end

s_unescape_string(str::AbstractString) = _sprint(s_print_unescaped, str)

function s_print_escaped(io, str::AbstractString, esc::Union{AbstractString, AbstractChar})
    pos = 1
    while !str_done(str, pos)
        chr, pos = iterate(str, pos)
        chr == '\0'         ? print(io, "\\0") :
        chr == '\e'         ? print(io, "\\e") :
        chr == '\\'         ? print(io, "\\\\") :
        chr in esc          ? print(io, '\\', chr) :
        '\a' <= chr <= '\r' ? print(io, '\\', "abtnvfr"[Int(chr)-6]) :
        is_printable(chr)   ? print(io, chr) : print(io, "\\u{", outhex(chr%UInt32), "}")
    end
end

s_escape_string(str::AbstractString) = _sprint(s_print_escaped, str, '\"')

s_print(str::AbstractString, flg::Bool=false) =
    s_print(str, flg, flg ? s_unescape_str : s_unescape_legacy)

function s_print(str::AbstractString, flg::Bool, unescape::Function)
    sx = interpolated_parse_vec(str, unescape, flg)
    (length(sx) == 1 && isa(sx[1], String)
     ? Expr(:call, :print, sx[1])
     : Expr(:call, :print, sx...))
end

function interpolated_parse(str::AbstractString, strfun::TypeOrFunc, flg::Bool, unescape::Function,
                            p::Function)
    sx = interpolated_parse_vec(str, unescape, flg)
    ((length(sx) == 1 && isa(sx[1], String)) ? strfun(sx[1])
     : Expr(:call, strfun, Expr(:call, :sprint, p, sx...)))
end

function interpolated_parse_vec(s::AbstractString, unescape::Function, flg::Bool=false)
    sx = []
    i = j = 1
    while !str_done(s, j)
        c, k = iterate(s, j)
        if c == '\\' && !str_done(s, k)
            c = s[k]
            if c == '('
                # Handle interpolation
                is_empty(s[i:prevind(s, j)]) ||
                    push!(sx, unescape(s[i:prevind(s, j)]))
                ex, j = parse(Expr, s, k, greedy=false)
                check_expr(ex)
                push!(sx, esc(ex))
                i = j
            elseif haskey(interpolate, c)
                i = j = interpolate[c](sx, s, unescape, i, j, k)
            elseif c == '\\'
                # Move past second backslash
                c, j = iterate(s, k)
            elseif flg && c == '$'
                is_empty(s[i:prevind(s, j)]) ||
                    push!(sx, unescape(s[i:prevind(s, j)]))
                i = k
                # Move past \\, c should point to '$'
                c, j = iterate(s, k)
            else
                j = k
            end
        elseif flg && c == '$'
            is_empty(s[i:prevind(s, j)]) ||
                push!(sx, unescape(s[i:prevind(s, j)]))
            ex, j = parse(Expr, s, k, greedy=false)
            check_expr(ex)
            push!(sx, esc(ex))
            i = j
        else
            j = k
        end
    end
    is_empty(s[i:end]) ||
        push!(sx, unescape(s[i:end]))
    sx
end

function s_unescape_str(str)
    str = s_unescape_string(str)
    is_valid(String, str) ? str : throw_arg_err("Invalid UTF-8 sequence")
end
function s_unescape_legacy(str)
    str = _sprint(s_print_unescaped_legacy, str)
    is_valid(String, str) ? str : throw_arg_err("Invalid UTF-8 sequence")
end

interpolated_parse(str::AbstractString, strfun::TypeOrFunc, flg::Bool, u::Function) =
    interpolated_parse(str, strfun, flg, u, print)

interpolated_parse(str::AbstractString, strfun::TypeOrFunc, flg::Bool=false) =
    interpolated_parse(str, strfun, flg, flg ? s_unescape_legacy : s_unescape_str)

@api freeze

end # module StrLiterals


text:
module FlippingBitsGame

using Printf, Random
import Base.size, Base.show, Base.==

struct Configuration
    M::BitMatrix
end

Base.size(c::Configuration) = size(c.M)
function Base.show(io::IO, conf::Configuration)
    M = conf.M
    nrow, ncol = size(M)
    print(io, " " ^ 3)
    for c in 1:ncol
        @printf(io, "%3i", c)
    end
    println(io, "\n", " " ^ 4, "-" ^ 3ncol)
    for r in 1:nrow
        @printf(io, "%2i |", r)
        for c in 1:ncol
            @printf(io, "%2c ", ifelse(M[r, c], 'T', 'F'))
        end
        println(io)
    end
    return nothing
end
Base.:(==)(a::Configuration, b::Configuration) = a.M == b.M

struct Index{D}
    i::Int
end
const ColIndex = Index{:C}
const RowIndex = Index{:R}

function flipbits!(conf::Configuration, c::ColIndex)
    col = @view conf.M[:, c.i]
    @. col = !col
    return conf
end
function flipbits!(conf::Configuration, r::RowIndex)
    row = @view conf.M[r.i, :]
    @. row = !row
    return conf
end

randomconfig(nrow::Integer, ncol::Integer) = Configuration(bitrand(nrow, ncol))

function randommoves!(conf::Configuration, nflips::Integer)
    nrow, ncol = size(conf)
    for _ in Base.OneTo(nflips)
        if rand() < 0.5
            flipbits!(conf, ColIndex(rand(1:ncol)))
        else
            flipbits!(conf, RowIndex(rand(1:nrow)))
        end
    end
    return conf
end

function play()
    nrow::Int, ncol::Int = 0, 0
    while nrow < 2 || ncol < 2
        print("Insert the size of the matrix (nrow [> 1] *space* ncol [> 1]):")
        nrow, ncol = parse.(Int, split(readline()))
    end
    mat = randomconfig(nrow, ncol)
    obj = deepcopy(mat)
    randommoves!(obj, 100)
    nflips = 0
    while mat != obj
        println("\n", nflips, " flips until now.")
        println("Current configuration:")
        println(mat)
        println("Objective configuration:")
        println(obj)
        print("Insert R[ind] to flip row, C[ind] to flip a column, Q to quit: ")
        line  = readline()
        input = match(r"([qrc])(\d+)"i, line)
        if input ≢ nothing && all(input.captures .≢ nothing)
            dim = Symbol(uppercase(input.captures[1]))
            ind = Index{dim}(parse(Int, input.captures[2]))
            flipbits!(mat, ind)
            nflips += 1
        elseif occursin("q", line)
            println("\nSEE YOU SOON!")
            return
        else
            println("\nINPUT NOT VALID, RETRY!\n")
        end
    end
    println("\nSUCCED! In ", nflips, " flips.")
    println(mat)
    return
end

end  # module FlippingBitsGame

using .FlippingBitsGame

FlippingBitsGame.play()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# Usage:
#     julia benchmark/run_benchmarks.jl
using BenchmarkTools
using MIRT
using FFTW: fft
using Random
using LinearMapsAA: LinearMapAA

include("benchmark_utils.jl")

on_CI = haskey(ENV, "GITHUB_ACTIONS")

const SUITE = BenchmarkGroup()

for benchmark_file in [
    "mri.jl",
    "nufft.jl",
    "regularize.jl",
    "utility.jl",
]
    include(benchmark_file)
end


text:
using Distributions

#This calculates the power for an upper-tailed z-test
#Need to add support for choose different hypotheses: lower-tailed, upper-tailed, two-sided
function power_norm(n, alpha, mu_0, mu_a, sigma)
    power_calc = 1 - cdf(Normal(0,1),quantile(Normal(0,1),1-alpha)+(mu_0-mu_a)*sqrt(n)/sigma)
end

max_stars_repo_name: dshehtanian/julia-learn

max_issues_repo_name: dshehtanian/julia-learn

max_forks_repo_name: dshehtanian/julia-learn


text:
"""
$(TYPEDEF)

A spatial inertia, or inertia matrix, represents the mass distribution of a
rigid body.

A spatial inertia expressed in frame ``i`` is defined as:
```math
I^i =
\\int_{B}\\rho\\left(x\\right)\\left[\\begin{array}{cc}
\\hat{p}^{T}\\left(x\\right)\\hat{p}\\left(x\\right) & \\hat{p}\\left(x\\right)\\\\
\\hat{p}^{T}\\left(x\\right) & I
\\end{array}\\right]dx=\\left[\\begin{array}{cc}
J & \\hat{c}\\\\
\\hat{c}^{T} & mI
\\end{array}\\right]
```
where ``\\rho(x)`` is the density of point ``x``, and ``p(x)`` are the coordinates
of point ``x`` expressed in frame ``i``.
``J`` is the mass moment of inertia, ``m`` is the total mass, and ``c`` is the
'cross part', center of mass position scaled by ``m``.

!!! warning

    The `moment` field of a `SpatialInertia` is the moment of inertia **about the origin of its `frame`**,
    not about the center of mass.
"""
struct SpatialInertia{T}
    frame::CartesianFrame3D
    moment::SMatrix{3, 3, T, 9}
    cross_part::SVector{3, T} # mass times center of mass
    mass::T

    @inline function SpatialInertia{T}(frame::CartesianFrame3D, moment::AbstractMatrix, cross_part::AbstractVector, mass) where T
        new{T}(frame, moment, cross_part, mass)
    end
end

"""
$(SIGNATURES)

Construct a `SpatialInertia` by specifying:

* `frame`: the frame in which the spatial inertia is expressed.
* `moment`: the moment of inertia expressed in `frame` (i.e., in `frame`'s axes and about the origin of `frame`).
* `cross_part`: the center of mass expressed in `frame`, multiplied by the mass.
* `mass`: the total mass.

For more convenient construction of `SpatialInertia`s, consider using the keyword argument constructor instead.
"""
@inline function SpatialInertia(frame::CartesianFrame3D, moment::AbstractMatrix{T1}, cross_part::AbstractVector{T2}, mass::T3) where {T1, T2, T3}
    SpatialInertia{promote_type(T1, T2, T3)}(frame, moment, cross_part, mass)
end

"""
$(SIGNATURES)

Construct a `SpatialInertia` by specifying:

* `frame`: the frame in which the spatial inertia is expressed.
* one of:
  * `moment`: the moment of inertia expressed in `frame` (i.e., about the origin of `frame` and in `frame`'s axes).
  * `moment_about_com`: the moment of inertia about the center of mass, in `frame`'s axes.
* `com`: the center of mass expressed in `frame`.
* `mass`: the total mass.

The `com` and `mass` keyword arguments are required, as well as exactly one of `moment` and `moment_about_com`
"""
@inline function SpatialInertia(frame::CartesianFrame3D;
        moment::Union{AbstractMatrix, Nothing}=nothing,
        moment_about_com::Union{AbstractMatrix, Nothing}=nothing,
        com::AbstractVector,
        mass)
    if !((moment isa AbstractMatrix) ⊻ (moment_about_com isa AbstractMatrix))
        throw(ArgumentError("Exactly one of `moment` or `moment_about_com` is required."))
    end
    _com = SVector{3}(com)
    if moment !== nothing
        _moment = moment
    else
        _moment = SMatrix{3, 3}(moment_about_com)
        _moment -= mass * hat_squared(_com) # parallel axis theorem
    end
    SpatialInertia(frame, _moment, mass * _com, mass)
end

# SpatialInertia-specific functions
Base.eltype(::Type{SpatialInertia{T}}) where {T} = T

# Construct/convert given another SpatialInertia
SpatialInertia{T}(inertia::SpatialInertia{T}) where {T} = inertia
@inline function SpatialInertia{T}(inertia::SpatialInertia) where {T}
    SpatialInertia{T}(inertia.frame,
        SMatrix{3, 3, T}(inertia.moment),
        SVector{3, T}(inertia.cross_part),
        T(inertia.mass))
end
@inline Base.convert(::Type{S}, inertia::SpatialInertia) where {S<:SpatialInertia} = S(inertia)

# Construct/convert to SMatrix
function StaticArrays.SMatrix{6, 6, T, 36}(inertia::SpatialInertia) where {T}
    J = SMatrix{3, 3, T}(inertia.moment)
    C = hat(SVector{3, T}(inertia.cross_part))
    m = T(inertia.mass)
    [[J C]; [C' SMatrix{3, 3}(m * I)]]
end
StaticArrays.SMatrix{6, 6, T}(inertia::SpatialInertia) where {T} = SMatrix{6, 6, T, 36}(inertia)
StaticArrays.SMatrix{6, 6}(inertia::SpatialInertia{T}) where {T} = SMatrix{6, 6, T}(inertia)
StaticArrays.SMatrix(inertia::SpatialInertia) = SMatrix{6, 6}(inertia)
StaticArrays.SArray(inertia::SpatialInertia) = SMatrix(inertia)
Base.convert(::Type{A}, inertia::SpatialInertia) where {A<:SArray} = A(inertia)

function Base.convert(::Type{Matrix}, inertia::SpatialInertia) where {T<:Matrix}
    Base.depwarn("This convert method is deprecated. Please use `$T(SMatrix(inertia))` instead or
    reconsider whether conversion to Matrix is necessary.", :convert)
    T(SMatrix(inertia))
end

function Base.Array(inertia::SpatialInertia)
    Base.depwarn("This Array constructor is deprecated. Please use `Matrix(SMatrix(inertia))` instead or
    reconsider whether conversion to Matrix is necessary.", :Array)
    Matrix(SMatrix(inertia))
end

"""
$(SIGNATURES)

Return the center of mass of the `SpatialInertia` as a [`Point3D`](@ref).
"""
center_of_mass(inertia::SpatialInertia) = Point3D(inertia.frame, inertia.cross_part / inertia.mass)

function Base.show(io::IO, inertia::SpatialInertia)
    println(io, "SpatialInertia expressed in $(name_and_id(inertia.frame)):")
    println(io, "mass: $(inertia.mass)")
    println(io, "center of mass: $(center_of_mass(inertia))")
    print(io, "moment of inertia (about origin of $(name_and_id(inertia.frame)):\n$(inertia.moment)")
end

Base.zero(::Type{SpatialInertia{T}}, frame::CartesianFrame3D) where {T} = SpatialInertia(frame, zero(SMatrix{3, 3, T}), zero(SVector{3, T}), zero(T))
Base.zero(inertia::SpatialInertia) = zero(typeof(inertia), inertia.frame)

function Base.isapprox(x::SpatialInertia, y::SpatialInertia; atol = 1e-12)
    x.frame == y.frame && isapprox(x.moment, y.moment; atol = atol) && isapprox(x.cross_part, y.cross_part; atol = atol) && isapprox(x.mass, y.mass; atol = atol)
end

@inline function Base.:+(inertia1::SpatialInertia, inertia2::SpatialInertia)
    @framecheck(inertia1.frame, inertia2.frame)
    SpatialInertia(inertia1.frame,
        inertia1.moment + inertia2.moment,
        inertia1.cross_part + inertia2.cross_part,
        inertia1.mass + inertia2.mass)
end

"""
$(SIGNATURES)

Transform the `SpatialInertia` to a different frame.
"""
@inline function transform(inertia::SpatialInertia, t::Transform3D)
    @framecheck(t.from, inertia.frame)
    J = inertia.moment
    m = inertia.mass
    c = inertia.cross_part
    R = rotation(t)
    p = translation(t)
    cnew = R * c
    Jnew = hat_squared(cnew)
    cnew += m * p
    Jnew -= hat_squared(cnew)
    minv = m > 0 ? inv(m) : zero(m)
    Jnew *= minv
    Jnew += R * J * R'
    SpatialInertia(t.to, Jnew, cnew, convert(eltype(Jnew), m))
end

function Random.rand(::Type{<:SpatialInertia{T}}, frame::CartesianFrame3D) where {T}
    # Try to generate a random but physical moment of inertia by constructing it from its eigendecomposition.

    # Create random principal moments of inertia (about the center of mass).
    # Scale the inertias to make the length scale of the equivalent inertial ellipsoid roughly ~1 unit
    ixx = rand(T) / 10.
    iyy = rand(T) / 10.

    # Ensure that the principal moments of inertia obey the triangle inequalities:
    # http://www.mathworks.com/help/physmod/sm/mech/vis/about-body-color-and-geometry.html
    lb = abs(ixx - iyy)
    ub = ixx + iyy
    izz = rand(T) * (ub - lb) + lb

    # Randomly rotate the principal axes.
    R = rand(RotMatrix3{T})
    moment_about_com = R * Diagonal(SVector(ixx, iyy, izz)) * R'

    SpatialInertia(frame, moment_about_com=moment_about_com, com=rand(SVector{3, T}) .- T(0.5), mass=rand(T))
end

"""
$(SIGNATURES)

Compute the mechanical power associated with a pairing of a wrench and a twist.
"""
LinearAlgebra.dot(w::Wrench, t::Twist) = begin @framecheck(w.frame, t.frame); dot(angular(w), angular(t)) + dot(linear(w), linear(t)) end
LinearAlgebra.dot(t::Twist, w::Wrench) = dot(w, t)

for (ForceSpaceMatrix, ForceSpaceElement) in (:MomentumMatrix => :Momentum, :MomentumMatrix => :Wrench, :WrenchMatrix => :Wrench)
    # MomentumMatrix * velocity vector --> Momentum
    # MomentumMatrix * acceleration vector --> Wrench
    # WrenchMatrix * dimensionless multipliers --> Wrench
    @eval function $ForceSpaceElement(mat::$ForceSpaceMatrix, x::AbstractVector)
        $ForceSpaceElement(mat.frame, SVector{3}(angular(mat) * x), SVector{3}(linear(mat) * x))
    end
end

@inline function Base.:*(inertia::SpatialInertia, twist::Twist)
    @framecheck(inertia.frame, twist.frame)
    ang, lin = mul_inertia(inertia.moment, inertia.cross_part, inertia.mass, angular(twist), linear(twist))
    Momentum(inertia.frame, ang, lin)
end

@inline function Base.:*(inertia::SpatialInertia, jac::GeometricJacobian)
    @framecheck(inertia.frame, jac.frame)
    Jω = angular(jac)
    Jv = linear(jac)
    J = inertia.moment
    m = inertia.mass
    c = inertia.cross_part
    ang = J * Jω + colwise(×, c, Jv)
    lin = m * Jv - colwise(×, c, Jω)
    MomentumMatrix(inertia.frame, ang, lin)
end

"""
$(SIGNATURES)

Apply the Newton-Euler equations to find the external wrench required to
make a body with spatial inertia ``I``, which has twist ``T`` with respect
to an inertial frame, achieve spatial acceleration ``\\dot{T}``.

This wrench is also equal to the rate of change of momentum of the body.
"""
@inline function newton_euler(inertia::SpatialInertia, spatial_accel::SpatialAcceleration, twist::Twist)
    I = inertia
    T = twist
    Ṫ = spatial_accel

    body = Ṫ.body
    base = Ṫ.base # TODO: should assert that this is an inertial frame somehow
    frame = Ṫ.frame

    @framecheck(I.frame, frame)
    @framecheck(T.body, body)
    @framecheck(T.base, base)
    @framecheck(T.frame, frame)

    ang, lin = mul_inertia(I.moment, I.cross_part, I.mass, angular(Ṫ), linear(Ṫ))
    angular_momentum, linear_momentum = mul_inertia(I.moment, I.cross_part, I.mass, angular(T), linear(T))
    ang += angular(T) × angular_momentum + linear(T) × linear_momentum
    lin += angular(T) × linear_momentum
    Wrench(frame, ang, lin)
end

@inline torque!(τ::AbstractVector, jac::GeometricJacobian, wrench::Wrench) = mul!(τ, transpose(jac), wrench)

@inline function torque(jac::GeometricJacobian, wrench::Wrench)
    τ = Vector{promote_type(eltype(jac), eltype(wrench))}(undef, size(jac, 2))
    torque!(τ, jac, wrench)
    τ
end

for (MatrixType, VectorType) in (:WrenchMatrix => :(Union{Twist, SpatialAcceleration}), :GeometricJacobian => :(Union{Momentum, Wrench}))
    @eval begin
        @inline function LinearAlgebra.mul!(
                dest::AbstractVector{T},
                transposed_mat::LinearAlgebra.Transpose{<:Any, <:$MatrixType},
                vec::$VectorType) where T
            mat = parent(transposed_mat)
            @boundscheck length(dest) == size(mat, 2) || throw(DimensionMismatch())
            @framecheck mat.frame vec.frame

            mat_angular = angular(mat)
            mat_linear = linear(mat)
            vec_angular = angular(vec)
            vec_linear = linear(vec)

            @inbounds begin
                @simd for row in eachindex(dest)
                    dest[row] =
                        mat_angular[1, row] * vec_angular[1] +
                        mat_angular[2, row] * vec_angular[2] +
                        mat_angular[3, row] * vec_angular[3]
                end
                @simd for row in eachindex(dest)
                    dest[row] +=
                        mat_linear[1, row] * vec_linear[1] +
                        mat_linear[2, row] * vec_linear[2] +
                        mat_linear[3, row] * vec_linear[3]
                end
            end
            dest
        end

        function Base.:*(transposed_mat::LinearAlgebra.Transpose{<:Any, <:$MatrixType}, vec::$VectorType)
            mat = parent(transposed_mat)
            @framecheck mat.frame vec.frame
            transpose(angular(mat)) * angular(vec) + transpose(linear(mat)) * linear(vec)
        end
    end
end

for ForceSpaceMatrix in (:MomentumMatrix, :WrenchMatrix)
    for (A, B) in ((ForceSpaceMatrix, :GeometricJacobian), (:GeometricJacobian, ForceSpaceMatrix))
        @eval begin
            function Base.:*(at::LinearAlgebra.Transpose{<:Any, <:$A}, b::$B)
                a = parent(at)
                @framecheck a.frame b.frame
                transpose(angular(a)) * angular(b) + transpose(linear(a)) * linear(b)
            end
            function Base.:*(a::$A, bt::LinearAlgebra.Transpose{<:Any, <:$B})
                b = parent(bt)
                @framecheck a.frame b.frame
                angular(a) * transpose(angular(b)) + linear(a) * transpose(linear(b))
            end
        end
    end
end

"""
$(SIGNATURES)

Compute the kinetic energy of a body with spatial inertia ``I``, which has
twist ``T`` with respect to an inertial frame.
"""
@inline function kinetic_energy(inertia::SpatialInertia, twist::Twist)
    @framecheck(inertia.frame, twist.frame)
    # TODO: should assert that twist.base is an inertial frame somehow
    ω = angular(twist)
    v = linear(twist)
    J = inertia.moment
    c = inertia.cross_part
    m = inertia.mass
    (ω ⋅ (J * ω) + v ⋅ (m * v + 2 * (ω × c))) / 2
end


text:
using ..PyPreferences: PyPreferences, use_system

import Conda
import Libdl
using CompilePreferences
using Pkg.Artifacts
using VersionParsing

struct PythonPreferences
    python::Union{Nothing,String}
    inprocess::Bool
    conda::Bool
    # jll::Bool
end

# Fix the environment for running `python`, and setts IO encoding to UTF-8.
# If cmd is the Conda python, then additionally removes all PYTHON* and
# CONDA* environment variables.
function pythonenv(cmd::Cmd)
    @assert cmd.env === nothing  # TODO: handle non-nothing case
    env = copy(ENV)
    if dirname(cmd.exec[1]) == abspath(Conda.PYTHONDIR)
        pythonvars = String[]
        for var in keys(env)
            if startswith(var, "CONDA") || startswith(var, "PYTHON")
                push!(pythonvars, var)
            end
        end
        for var in pythonvars
            pop!(env, var)
        end
    end
    # set PYTHONIOENCODING when running python executable, so that
    # we get UTF-8 encoded text as output (this is not the default on Windows).
    env["PYTHONIOENCODING"] = "UTF-8"
    setenv(cmd, env)
end

pyvar(python::AbstractString, mod::AbstractString, var::AbstractString) =
    chomp(read(pythonenv(`$python -c "import $mod; print($mod.$(var))"`), String))

pyconfigvar(python::AbstractString, var::AbstractString) =
    pyvar(python, "distutils.sysconfig", "get_config_var('$(var)')")
pyconfigvar(python, var, default) =
    let v = pyconfigvar(python, var)
        v == "None" ? default : v
    end

function pythonhome_of(pyprogramname::AbstractString)
    if Sys.iswindows()
        # PYTHONHOME tells python where to look for both pure python
        # and binary modules.  When it is set, it replaces both
        # `prefix` and `exec_prefix` and we thus need to set it to
        # both in case they differ. This is also what the
        # documentation recommends.  However, they are documented
        # to always be the same on Windows, where it causes
        # problems if we try to include both.
        script = """
        import sys
        if hasattr(sys, "base_exec_prefix"):
            sys.stdout.write(sys.base_exec_prefix)
        else:
            sys.stdout.write(sys.exec_prefix)
        """
    else
        script = """
        import sys
        if hasattr(sys, "base_exec_prefix"):
            sys.stdout.write(sys.base_prefix)
            sys.stdout.write(":")
            sys.stdout.write(sys.base_exec_prefix)
        else:
            sys.stdout.write(sys.prefix)
            sys.stdout.write(":")
            sys.stdout.write(sys.exec_prefix)
        """
        # https://docs.python.org/3/using/cmdline.html#envvar-PYTHONHOME
    end
    return read(pythonenv(`$pyprogramname -c $script`), String)
end
# To support `venv` standard library (as well as `virtualenv`), we
# need to use `sys.base_prefix` and `sys.base_exec_prefix` here.
# Otherwise, initializing Python in `__init__` below fails with
# unrecoverable error:
#
#   Fatal Python error: initfsencoding: unable to load the file system codec
#   ModuleNotFoundError: No module named 'encodings'
#
# This is because `venv` does not symlink standard libraries like
# `virtualenv`.  For example, `lib/python3.X/encodings` does not
# exist.  Rather, `venv` relies on the behavior of Python runtime:
#
#   If a file named "pyvenv.cfg" exists one directory above
#   sys.executable, sys.prefix and sys.exec_prefix are set to that
#   directory and it is also checked for site-packages
#   --- https://docs.python.org/3/library/venv.html
#
# Thus, we need point `PYTHONHOME` to `sys.base_prefix` and
# `sys.base_exec_prefix`.  If the virtual environment is created by
# `virtualenv`, those `sys.base_*` paths point to the virtual
# environment.  Thus, above code supports both use cases.
#
# See also:
# * https://docs.python.org/3/library/venv.html
# * https://docs.python.org/3/library/site.html
# * https://docs.python.org/3/library/sys.html#sys.base_exec_prefix
# * https://github.com/JuliaPy/PyCall.jl/issues/410

python_version_of(python) = vparse(pyvar(python, "platform", "python_version()"))

function find_libpython_py_path()
    dir = artifact"pyjulia_src_julia"
    return joinpath(dir, only(readdir(dir)), "find_libpython.py")
end

function exec_find_libpython(python::AbstractString, options, verbose::Bool)
    # Do not inline `@__DIR__` into the backticks to expand correctly.
    # See: https://github.com/JuliaLang/julia/issues/26323
    script = find_libpython_py_path()
    cmd = `$python $script $options`
    if verbose
        cmd = `$cmd --verbose`
    end
    return readlines(pythonenv(cmd))
end

# return libpython path, libpython pointer
function find_libpython(
    python::AbstractString;
    _dlopen = Libdl.dlopen,
    verbose::Bool = false,
)
    dlopen_flags = Libdl.RTLD_LAZY | Libdl.RTLD_DEEPBIND | Libdl.RTLD_GLOBAL

    libpaths = exec_find_libpython(python, `--list-all`, verbose)
    for lib in libpaths
        try
            return (lib, _dlopen(lib, dlopen_flags))
        catch e
            @warn "Failed to `dlopen` $lib" exception = (e, catch_backtrace())
        end
    end
    @warn """
    Python (`find_libpython.py`) failed to find `libpython`.
    Falling back to `Libdl`-based discovery.
    """

    # Try all candidate libpython names and let Libdl find the path.
    # We do this *last* because the libpython in the system
    # library path might be the wrong one if multiple python
    # versions are installed (we prefer the one in LIBDIR):
    libs = exec_find_libpython(python, `--candidate-names`, verbose)
    for lib in libs
        lib = splitext(lib)[1]
        try
            libpython = _dlopen(lib, dlopen_flags)
            return (Libdl.dlpath(libpython), libpython)
        catch e
            @debug "Failed to `dlopen` $lib" exception = (e, catch_backtrace())
        end
    end

    return nothing, nothing
end

function PyPreferences.use_system(python::AbstractString = "python3")
    return set(python = python)
end

function PyPreferences.use_conda()
    Conda.add("numpy")
    return set(conda = true)
end

function PyPreferences.use_inprocess()
    return set(inprocess = true)
end

conda_python_fullpath() =
    abspath(Conda.PYTHONDIR, "python" * (Sys.iswindows() ? ".exe" : ""))

#=
function use_jll()
end
=#

set(; python = nothing, inprocess = false, conda = false) =
    set(PythonPreferences(python, inprocess, conda))

function set(prefs::PythonPreferences)
    @save_preferences!(Dict(prefs))
    recompile()
    return prefs
end

function Base.Dict(prefs::PythonPreferences)
    rawprefs = Dict{String,Any}()
    if prefs.python !== nothing
        rawprefs["python"] = prefs.python
    end
    if prefs.inprocess
        rawprefs["inprocess"] = true
    end
    if prefs.conda
        rawprefs["conda"] = true
    end
    return rawprefs
end

PythonPreferences(rawprefs::AbstractDict) = PythonPreferences(
    get(rawprefs, "python", nothing),
    get(rawprefs, "inprocess", false),
    get(rawprefs, "conda", false),
)

function _load_python_preferences()
    # TODO: lookup v#.#?
    rawprefs = @load_preferences()
    isempty(rawprefs) && return nothing
    return PythonPreferences(rawprefs)
end

function load_pypreferences_code()
    return """
     $(Base.load_path_setup_code())
     PyPreferences = Base.require(Base.PkgId(
         Base.UUID("cc9521c6-0242-4dda-8d66-c47a9d9eec02"),
         "PyPreferences",
     ))
     """
end

function include_stdin_cmd()
    return ```
    $(Base.julia_cmd())
    --startup-file=no
    -e "include_string(Main, read(stdin, String))"
    ```
end

function recompile()
    code = """
    $(load_pypreferences_code())
    PyPreferences.assert_configured()
    """
    cmd = include_stdin_cmd()
    open(cmd; write = true) do io
        write(io, code)
    end
    return
end

function PyPreferences.instruction_message()
    return """
    PyPreferences.jl is not configured properly. Run:
        using Pkg
        Pkg.add("PyPreferences")
        using PyPreferences
        @doc PyPreferences
    for usage.
    """
end

function PyPreferences.assert_configured()
    if (
        PyPreferences.python === nothing ||
        PyPreferences.python_fullpath === nothing ||
        PyPreferences.libpython === nothing ||
        PyPreferences.python_version === nothing ||
        PyPreferences.PYTHONHOME === nothing
    )
        error(PyPreferences.instruction_message())
    end
end

function setup_non_failing()
    python = nothing
    inprocess = false
    conda = false
    python_fullpath = nothing
    libpython = nothing
    python_version = nothing
    PYTHONHOME = nothing

    prefs = _load_python_preferences()
    if prefs !== nothing
        python = prefs.python
        inprocess = prefs.inprocess
        conda = prefs.conda
    end
    if !inprocess
        if conda
            python = python_fullpath = conda_python_fullpath()
        elseif python === nothing
            # TODO: mimic PyCall's deps/build.jl
            python = "python3"
        end

        try
            if python !== nothing
                python_fullpath = Sys.which(python)
            end
            if python_fullpath !== nothing
                libpython, = find_libpython(python_fullpath)
                python_version = python_version_of(python_fullpath)
                PYTHONHOME = pythonhome_of(python_fullpath)
            end
        catch err
            @error(
                "Failed to configure for `$python`",
                exception = (err, catch_backtrace())
            )
        end
    end
    if python === nothing
        python = python_fullpath
    end

    return (
        python = python,
        inprocess = inprocess,
        conda = conda,
        python_fullpath = python_fullpath,
        libpython = libpython,
        python_version = python_version,
        PYTHONHOME = PYTHONHOME,
    )
end

function status_inprocess()
    print("""
    python         : $(PyPreferences.python)
    inprocess      : $(PyPreferences.inprocess)
    conda          : $(PyPreferences.conda)
    python_fullpath: $(PyPreferences.python_fullpath)
    libpython      : $(PyPreferences.libpython)
    python_version : $(PyPreferences.python_version)
    PYTHONHOME     : $(PyPreferences.PYTHONHOME)
    """)
end

function PyPreferences.status()
    # TODO: compare with in-process values
    code = """
    $(load_pypreferences_code())
    PyPreferences.Implementations.status_inprocess()
    """
    cmd = include_stdin_cmd()
    open(pipeline(cmd; stdout = stdout, stderr = stderr); write = true) do io
        write(io, code)
    end
    return
end


text:
using SparseTransforms
using StatsBase, Random, Distributions
using Test

Random.seed!(1234)

function test_random()
    n = 8
    k = 2
    σ = 1e-5
    minpower, maxpower = 0.5, 5.0
    signal = get_random_sparse_signal(n, 2^k, σ, minpower, maxpower)

    println("True locations: ", signal.locs)
    println("True strengths: ", signal.strengths)

    b = get_b(signal; method=:simple)
    Ms = get_Ms(n, b; method=:simple)
    contents = Dict()
    for loc in signal.locs
        println(loc)
        for (i, M) in enumerate(Ms)
            j = expected_bin(dec_to_bin(loc, n), M)
            if haskey(contents, (i,j))
                append!(contents[(i,j)], loc)
            else
                contents[(i,j)] = [loc]
            end
            println("\t$i: $j")
        end
        println()
    end
    println("true contents:")
    for (k,v) in contents
        println("\t$k: $v")
    end

    methods_list = [:simple, :nso, :nso, :none]
    # methods_list = [:simple, :random, :mle, :none]
    spright_wht, used_size = spright(signal, methods_list; verbose=true, report=true)

    println("SPRIGHT result: ", spright_wht)
    println("Used $(used_size / 2^n) of all time samples")
    @test length(signal.locs) == length(spright_wht)
    for (l, s) in zip(signal.locs, signal.strengths)
        @test isapprox(spright_wht[l], s, atol=5*σ)
    end
end


text:
using StructuralCausalModels

ProjDir = @__DIR__
cd(ProjDir) #do

df = DataFrame!(CSV.File(scm_path("..", "data", "WaffleDivorce.csv"), delim=';'));
df = DataFrame(
  :s => df[:, :South],
  :a => df[:, :MedianAgeMarriage],
  :m => df[:, :Marriage],
  :w => df[:, :WaffleHouses],
  :d => df[:, :Divorce]
);

fname = scm_path("..", "examples", "SR", "SR6.4.3", "sr6.4.3.dot")
Sys.isapple() && run(`open -a GraphViz.app $(fname)`)

d = OrderedDict(
  :s => [:a, :m, :w],
  :a => [:m, :d],
  :m => [:d],
  :w => [:d]
);

dag = DAG("sr6.4.3", d; df=df);
show(dag)

adjustmentsets = adjustment_sets(dag, :w, :d)
println("Adjustment sets:\n")
adjustmentsets |> display
println()

function show_differences(res)
   if length(res) > 0
    println()
    res |> display
    println()
  end
end

debug = true

a = copy(dag.a)
ar = copy(a);
vars = names(a, 1)
zna = NamedArray(zeros(Int, size(a)), (vars, vars), ("Rows", "Cols"));

c = [:s]
m = Symbol[]

println()
s = StructuralCausalModels.update_s(a, c, debug)
#@assert s == [:s, :a, :m, :w, :d]

count = 0
while true
  local res
  global count += 1
  println("\ncount = $count\n")
  println()

  global ar
  at = copy(ar)

  println("\nupdate_21\n")
  ar21 = StructuralCausalModels.update_21(ar, m, debug);
  res = findall(ar21 .!== ar)
  show_differences(res)

  ar = copy(ar21);
  println("\nupdate_22\n")
  ar22 = StructuralCausalModels.update_22(ar, s, debug)
  res = findall(ar22 .!== zna)
  show_differences(res)

  ar .+= ar22

  println("\nupdate_23\n")
  ar23 = StructuralCausalModels.update_23(ar, m, debug)
  res = findall(ar23 .!== zna)
  show_differences(res)

  ar .+= ar23

  println("\nupdate_24\n")
  ar24 = StructuralCausalModels.update_24(ar, m, debug)
  res = findall(ar24 .!== zna)
  show_differences(res)

  ar .+= ar24

  println("\nupdate_25\n")
  ar25 = StructuralCausalModels.update_25(ar, m, debug)
  res = findall(ar25 .!== zna)
  show_differences(res)

  ar .+= ar25

  println("\nupdate_26\n")
  ar26 = StructuralCausalModels.update_26(ar, m, debug)
  res = findall(ar26 .!== zna)
  show_differences(res)

  ar .+= ar26

  println("\nupdate_27\n")
  ar27 = StructuralCausalModels.update_27(ar, s, debug)
  res = findall(ar27 .!== zna)
  show_differences(res)

  ar .+= ar27

  println("\nupdate_28\n")
  ar28 = StructuralCausalModels.update_28(ar, s, debug)
  res = findall(ar28 .!== zna)
  show_differences(res)

  ar .+= ar28

  println("\nupdate_29\n")
  ar29 = StructuralCausalModels.update_29(ar, m, debug)
  res = findall(ar29 .!== zna)
  show_differences(res)

  ar .+= ar29

  println("\nupdate_30\n")
  ar30 = StructuralCausalModels.update_30(ar, m, debug)
  res = findall(ar30 .!== zna)
  show_differences(res)

  ar .+= ar30

  ar == at && break

end

println("\nEnd of 1st while loop ($count iterations).\n")

ar_old = copy(ar)
println("\nupdate_a\n")
ar = StructuralCausalModels.update_a(ar, true)
res = findall(ar .!== ar_old)
show_differences(res)

ar_old = copy(ar)
println("\nupdate_b\n")
ar = StructuralCausalModels.update_b(ar, s, true)
res = findall(ar .!== ar_old)
show_differences(res)

an = copy(ar)
count = 0
while true
  global count += 1
  println("\ncount = $count\n")
  println()

  global ar
  global an
  global res
  at = copy(ar)

  println("\nupdate_27n\n")
  a27n = StructuralCausalModels.update_27n(ar, an)
  res = findall(a27n .!== zna)
  show_differences(res)

  an .+= a27n
  ar .+= a27n

  println("\nupdate_22n\n")
  a22n = StructuralCausalModels.update_22n(ar, an)
  res = findall(a22n .!== zna)
  show_differences(res)

  an .+= a22n
  ar .+= a22n

  at == ar && break

end

println("\nEnd of 2nd while loop ($count iterations).\n")

println("\nupdate_i\n")
ai = StructuralCausalModels.update_i(ar)
res = findall(ar .!== ai)
show_differences(res)

ar = copy(ai)

select = setdiff(vars, vcat(c, m))
println(select)
fr = ar[select, select]
println()

display(fr)
println()


text:
using WAV
using ArgParse

push!(LOAD_PATH, pwd())

using Plotting

# Main entry
function main(wav_file)

	# Read the WAV file
	s, fs = wavread(wav_file)

	plot_wave(s, fs)
	plot_frec(s, fs)
	plot_spectrogram(s, fs)
end

if PROGRAM_FILE == "SignalAnalysis.jl"

	# Set backend renderer
	set_plot()

	# Parse the argument
	s = ArgParseSettings("Audio dialing analisys")

	@add_arg_table s begin
		"--audio"
		help = "The audio file (WAV)"
		required = true
	end

	parsed_args = parse_args(ARGS, s)

	main(parsed_args["audio"])
end

max_stars_repo_name: PatricioIribarneCatella/TP-Telefono

max_issues_repo_name: PatricioIribarneCatella/TP-Telefono

max_forks_repo_name: PatricioIribarneCatella/TP-Telefono


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "VerizonEctoken"
version = v"0.1.0"

# Collection of sources required to complete build
sources = [
    "https://github.com/VerizonDigital/ectoken.git" =>
    "7b8812d476f5be5b290fe2832859b9b7636f43ae",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
export OPENSSL_INCLUDE="-I $prefix/include"
export OPENSSL_LIBS="-lssl -lcrypto"
cd ectoken/c-ectoken/ecencrypt/
gcc -m64 -O2 -Wall -Werror -std=gnu99 ec_encrypt.c ectoken_v3.c base64.c -o 64/ectoken3 $OPENSSL_LIBS $OPENSSL_INCLUDE -lm -lpthread -ldl
cp 64/ectoken3 $prefix/bin/
exit

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64, libc=:glibc)
]

# The products that we will ensure are always built
products = [
    ExecutableProduct("ectoken3", :ectoken3)
]

# Dependencies that must be installed before this package can be built
dependencies = [
  "OpenSSL_jll"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
sbmlfiles = [
    # a test model from BIGG
    (
        joinpath(@__DIR__, "data", "e_coli_core.xml"),
        "http://bigg.ucsd.edu/static/models/e_coli_core.xml",
        "b4db506aeed0e434c1f5f1fdd35feda0dfe5d82badcfda0e9d1342335ab31116",
        72,
        95,
        fill(1000.0, 3),
    ),
    # a relatively new non-curated model from biomodels
    (
        joinpath(@__DIR__, "data", "T1M1133.xml"),
        "https://www.ebi.ac.uk/biomodels/model/download/MODEL1909260004.4?filename=T1M1133.xml",
        "2b1e615558b6190c649d71052ac9e0dc1635e3ad281e541bc7d4fdf2892a5967",
        2517,
        3956,
        fill(1000.0, 3),
    ),
    # a curated model from biomodels
    (
        joinpath(@__DIR__, "data", "Dasgupta2020.xml"),
        "https://www.ebi.ac.uk/biomodels/model/download/BIOMD0000000973.3?filename=Dasgupta2020.xml",
        "958b131d4df2f215dae68255433542f228601db0326d26a54efd08ddcf823489",
        2,
        6,
        fill(Inf, 3),
    ),
    # a cool model with `time` from SBML testsuite
    (
        joinpath(@__DIR__, "data", "sbml00852.xml"),
        "https://raw.githubusercontent.com/sbmlteam/sbml-test-suite/master/cases/semantic/00852/00852-sbml-l3v2.xml",
        "d013765aa358d265941420c2e3d81fcbc24b0aa4e9f39a8dc8852debd1addb60",
        4,
        3,
        fill(Inf, 3),
    ),
    # another model from SBML suite, with initial concentrations
    (
        joinpath(@__DIR__, "data", "sbml00374.xml"),
        "https://raw.githubusercontent.com/sbmlteam/sbml-test-suite/master/cases/semantic/00374/00374-sbml-l3v2.xml",
        "424683eea6bbb577aad855d95f2de5183a36e296b06ba18b338572cd7dba6183",
        4,
        2,
        fill(Inf, 2),
    ),
    # this contains some special math
    (
        joinpath(@__DIR__, "data", "sbml01565.xml"),
        "https://raw.githubusercontent.com/sbmlteam/sbml-test-suite/master/cases/semantic/01565/01565-sbml-l3v1.xml",
        "14a80fbce316eea2adb566f67b4668ad151db8954e487309852ece7f730c8c99",
        104,
        52,
        fill(Inf, 3),
    ),
    # this contains l3v1-incompatible contents
    (
        joinpath(@__DIR__, "data", "sbml01289.xml"),
        "https://raw.githubusercontent.com/sbmlteam/sbml-test-suite/master/cases/semantic/01289/01289-sbml-l3v2.xml",
        "35ffa072052970b92fa358ee0f5750394ad74958e889cb85c98ed238642de4d0",
        0,
        0,
        Float64[],
    ),
    # this contains a relational operator
    (
        joinpath(@__DIR__, "data", "sbml00191.xml"),
        "https://raw.githubusercontent.com/sbmlteam/sbml-test-suite/master/cases/semantic/00191/00191-sbml-l3v2.xml",
        "c474e94888767d70f9e9e03b32778f18069641563953de60dabac7daa7f481ce",
        4,
        2,
        fill(Inf, 2),
    ),
    # expandInitialAssignments converter gives some warning
    (
        joinpath(@__DIR__, "data", "01234-sbml-l3v2.xml"),
        "https://raw.githubusercontent.com/sbmlteam/sbml-test-suite/52d94baf97a005b6e1fdbdb6116f5c7b4a8a100c/cases/semantic/01234/01234-sbml-l3v2.xml",
        "9610ef29f2d767af627042a15bde505b068ab75bbf00b8983823800ea8ef67c8",
        0,
        0,
        Float64[],
    ),
    (
        joinpath(@__DIR__, "data", "00489-sbml-l3v2.xml"),
        "https://raw.githubusercontent.com/sbmlteam/sbml-test-suite/master/cases/semantic/00489/00489-sbml-l3v2.xml",
        "dab2bce4e5036fa47ad8137055ca5f6dec6dfcb183542ce38573ca2e5a615813",
        3,
        2,
        fill(Inf, 2),
    ),
    # has listOfEvents
    (
        joinpath(@__DIR__, "data", "00026-sbml-l3v2.xml"),
        "https://raw.githubusercontent.com/sbmlteam/sbml-test-suite/master/cases/semantic/00026/00026-sbml-l3v2.xml",
        "991381015d9408164bf00206848ba5796d0c86dc055be91968cd7f0f68daa903",
        2,
        1,
        fill(Inf, 1),
    ),
    # has all rules types
    (
        joinpath(@__DIR__, "data", "00983-sbml-l3v2.xml"),
        "https://raw.githubusercontent.com/sbmlteam/sbml-test-suite/master/cases/semantic/00983/00983-sbml-l3v2.xml",
        "b84e53cc48edd5afc314e17f05c6a30a509aadb9486b3d788c7cf8df82a7527f",
        0,
        0,
        Float64[],
    ),
    # has all kinds of default model units
    (
        joinpath(@__DIR__, "data", "00054-sbml-l3v2.xml"),
        "https://raw.githubusercontent.com/sbmlteam/sbml-test-suite/master/cases/semantic/00054/00054-sbml-l3v2.xml",
        "987038ec9bb847123c41136928462d7ed05ad697cc414cab09fcce9f5bbc8e73",
        3,
        2,
        fill(Inf,2),
    ),
]

@testset "Loading of models from various sources - $(reader)" for reader in (
    readSBML,
    readSBMLFromString,
)
    for (sbmlfile, url, hash, expected_mets, expected_rxns, expected_3_ubs) in sbmlfiles
        if !isfile(sbmlfile)
            Downloads.download(url, sbmlfile)
        end

        cksum = bytes2hex(sha256(open(sbmlfile)))
        if cksum != hash
            @warn "The downloaded model `$sbmlfile' seems to be different from the expected one. Tests will likely fail." cksum
        end

        @testset "Loading of $sbmlfile" begin
            mdl = if reader === readSBML
                readSBML(sbmlfile)
            else
                readSBMLFromString(readchomp(sbmlfile))
            end

            @test typeof(mdl) == Model

            mets, rxns, _ = stoichiometry_matrix(mdl)

            @test length(mets) == expected_mets
            @test length(rxns) == expected_rxns

            lbs, ubs = flux_bounds(mdl)
            @test length(lbs) == expected_rxns
            @test length(ubs) == expected_rxns
            @test first.(ubs)[1:min(3, length(ubs))] == expected_3_ubs

            ocs = flux_objective(mdl)
            @test length(ocs) == expected_rxns
        end
    end
end

@testset "readSBMLFromString" begin
    @test_logs (:error, r"^SBML reported error") @test_throws AssertionError readSBMLFromString(
        "",
    )
end

@testset "Time variables in math" begin
    # this test is here mainly for keeping a magical constant that we need for
    # parsing time synced with libsbml source
    contains_time(x::SBML.MathTime) = true
    contains_time(x::SBML.MathApply) = any(contains_time.(x.args))
    contains_time(_) = false

    m = readSBML(joinpath(@__DIR__, "data", "sbml00852.xml"))
    @test all(contains_time.(r.kinetic_math for (_, r) in m.reactions))
end

@testset "Units" begin
    m = readSBML(joinpath(@__DIR__, "data", "sbml00852.xml"))
    @test SBML.unitful(m.units["volume"]) == 1 * u"L"
    @test SBML.unitful(m.units["time"]) == 1 * u"s"
    @test SBML.unitful(m.units["substance"]) == 1 * u"mol"

    m = readSBML(joinpath(@__DIR__, "data", "custom.xml"))
    @test SBML.unitful(m.units["non_existent"]) == 0.00314
    @test SBML.unitful(m.units["no_dimensions"]) == 20.0
end

@testset "Initial amounts and concentrations" begin
    m = readSBML(joinpath(@__DIR__, "data", "sbml00852.xml"))

    @test all(isnothing(ic) for (k, ic) in SBML.initial_concentrations(m))
    @test length(SBML.initial_amounts(m)) == 4
    @test isapprox(sum(ia for (sp, ia) in SBML.initial_amounts(m)), 0.001)
    @test isapprox(
        sum(ic for (sp, ic) in SBML.initial_concentrations(m, convert_amounts = true)),
        0.001,
    )

    m = readSBML(joinpath(@__DIR__, "data", "sbml00374.xml"))

    @test all(isnothing(ic) for (k, ic) in SBML.initial_amounts(m))
    @test length(SBML.initial_concentrations(m)) == 4
    @test isapprox(sum(ic for (sp, ic) in SBML.initial_concentrations(m)), 0.00208)
    @test isapprox(
        sum(ia for (sp, ia) in SBML.initial_amounts(m, convert_concentrations = true)),
        0.25 * 0.00208,
    )
end

@testset "Initial assignments" begin
    m = readSBML(joinpath(@__DIR__, "data", "00489-sbml-l3v2.xml"))
    @test m.initial_assignments ==
          Dict("S1" => SBML.MathApply("*", [SBML.MathVal{Int32}(2), SBML.MathIdent("p1")]))

    m = readSBML(joinpath(@__DIR__, "data", "sbml01289.xml"))
    @test m.initial_assignments == Dict(
        "p2" => SBML.MathApply("gt5", [SBML.MathVal{Int32}(8)]),
        "p1" => SBML.MathApply("gt5", [SBML.MathVal{Int32}(3)]),
    )
end

@testset "Rules" begin
    m = readSBML(joinpath(@__DIR__, "data", "Dasgupta2020.xml"))
    @test m.rules == [
        SBML.AssignmentRule(
            "s",
            SBML.MathApply(
                "/",
                [
                    SBML.MathApply(
                        "-",
                        [SBML.MathIdent("ModelValue_6"), SBML.MathIdent("P")],
                    ),
                    SBML.MathIdent("N"),
                ],
            ),
        ),
    ]

    m = readSBML(joinpath(@__DIR__, "data", "00983-sbml-l3v2.xml"))
    @test m.rules == [
        SBML.RateRule("x", SBML.MathVal{Int32}(1)),
        SBML.AlgebraicRule(
            SBML.MathApply(
                "+",
                [
                    SBML.MathApply("*", [SBML.MathVal{Int32}(-1), SBML.MathIdent("temp")]),
                    SBML.MathApply("/", [SBML.MathTime("time"), SBML.MathVal{Int32}(2)]),
                ],
            ),
        ),
        SBML.AssignmentRule(
            "y",
            SBML.MathApply("delay", [SBML.MathIdent("x"), SBML.MathIdent("temp")]),
        ),
    ]
end

@testset "Constraints" begin
    m = readSBML(joinpath(@__DIR__, "data", "custom.xml"))
    @test only(m.constraints) == SBML.Constraint(
        SBML.MathApply(
            "and",
            [
                SBML.MathApply("lt", [SBML.MathVal{Float64}(1.0), SBML.MathIdent("S1")]),
                SBML.MathApply("lt", [SBML.MathIdent("S1"), SBML.MathVal{Float64}(100.0)]),
            ],
        ),
        " Species S1 is out of range. ",
    )
end

@testset "Extensive kinetic math" begin
    m = readSBML(joinpath(@__DIR__, "data", "sbml00852.xml"))

    subterm =
        SBML.extensive_kinetic_math(m, m.reactions["reaction1"].kinetic_math).args[1].args[2]
    @test subterm.fn == "/"
    @test subterm.args[1] == SBML.MathIdent("S1")
    @test isapprox(subterm.args[2].val, 1.0)
end

@testset "logBase and root math functions" begin
    m = readSBML(joinpath(@__DIR__, "data", "sbml01565.xml"))

    @test convert(Num, m.reactions["J23"].kinetic_math) == 0.0

    @variables S29 S29b
    @test isequal(convert(Num, m.reactions["J29"].kinetic_math), 2.0 * S29 * S29b)
end

@testset "rationals in math" begin
    m = readSBML(joinpath(@__DIR__, "data", "custom.xml"))
    @test m.reactions["rationalTest"].kinetic_math.val == 1 // 5000
end

@testset "converters work and fail gracefully" begin
    @test_logs (:error, r"^SBML reported error:") (:error, r"^SBML reported error:") @test_throws ErrorException readSBML(
        joinpath(@__DIR__, "data", "sbml01289.xml"),
        doc -> begin
            set_level_and_version(3, 1)(doc)
            convert_simplify_math(doc)
        end,
    )

    test_math =
        readSBML(
            joinpath(@__DIR__, "data", "sbml01565.xml"),
            libsbml_convert("expandInitialAssignments"),
        ).reactions["J31"].kinetic_math

    @test test_math.args[2].fn == "sin"
    @test test_math.args[2].args[1].val == 2.1

    @test_logs (:warn,) (:warn,) (:warn,) (:warn,) readSBML(
        joinpath(@__DIR__, "data", "01234-sbml-l3v2.xml"),
        doc ->
            libsbml_convert("expandInitialAssignments", ["Fatal", "Error", "Warning"])(doc),
    )
end

@testset "relational operators are decoded correctly" begin
    test_math =
        readSBML(joinpath(@__DIR__, "data", "sbml00191.xml")).reactions["reaction2"].kinetic_math

    @test test_math.args[2].fn == "geq"
end

@testset "custom show" begin
    m = readSBML(joinpath(@__DIR__, "data", "custom.xml"))
    @test repr(MIME("text/plain"), m) ==
          "Model with 1 reactions, 0 species, and 0 parameters."
    @test eval(Meta.parse(repr(m))) isa SBML.Model
end

@testset "events" begin
    m = readSBML(joinpath(@__DIR__, "data", "00026-sbml-l3v2.xml"))
    @test length(m.events) == 1
end

@testset "model units" begin
    m = readSBML(joinpath(@__DIR__, "data", "00054-sbml-l3v2.xml"))
    @test m.area_units == "area"
    @test m.extent_units == "substance"
    @test m.length_units == "metre"
    @test m.substance_units == "substance"
    @test m.time_units == "second"
    @test m.volume_units == "volume"
end


text:
using DeviationsLH
using Test

include("deviation_test_setup.jl");

@testset "All" begin
    include("deviation_test.jl");
    include("dev_vector_test.jl");
    include("change_table_test.jl");
end

# -------------