text:
"""
    abstract type AbstractPEPO{S<:IndexSpace} <: AbstractArray{S,2} end

Abstract supertype of all PEPO-like objects, i.e. networks of `PEPOTensor` types.
"""
abstract type AbstractPEPO{S<:IndexSpace} <: AbstractMatrix{S} end

# TODO: Change requirement on PEPOTensor to trace = 1, check for Hermiticity.
"""
    struct PEPOTensor{S<:IndexSpace} <: AbstractTensorMap{S,3,3}

Concrete type containing an `AbstractTensorMap{S,3,3}` with unit norm. By convention,
the codomain is ordered (1,2,out) and the domain is ordered (3,4,in).
"""
struct PEPOTensor{
    S<:IndexSpace,
    I<:Sector,
    A<:Union{<:DenseMatrix,TensorKit.SectorDict{I,<:DenseMatrix}},
    F₁,
    F₂,
} <: AbstractPhysicalTensor{S,3,3, I, A, F₁, F₂}
    tensor::TensorMap{S,3,3,I,A,F₁,F₂}
    function PEPOTensor(
        tensor::TensorMap{S,3,3,I,A,F₁,F₂}
    ) where {
        S<:IndexSpace,
        I<:Sector,
        A<:Union{<:DenseMatrix,TensorKit.SectorDict{I,<:DenseMatrix}},
        F₁,
        F₂,
    }
        if norm(tensor) ≈ one(ComplexF64)
            return new{S,I,A,F₁,F₂}(tensor)
        else
            @warn "TensorMap does not have norm equal to 1"
            return new{S,I,A,F₁,F₂}(tensor)
        end
    end
end

"""
    PEPOTensor(D::Int, x::Real, y::Real, z::Real)
    PEPOTensor(D::Int, a::Vector{<:Real})

Construct a non-symmetric `PEPOTensor` with bond dimension `D` corresponding to a given
Bloch vector. Each element of the Bloch vector can be provided as a seperate arguments,
or a vector of real numbers can be provided.
"""
function PEPOTensor(D::Int, x::Real, y::Real, z::Real)
    return PEPOTensor(bloch_pepo(D, x, y, z))
end
PEPOTensor(D::Int, a::Vector{<:Real}) = PEPOTensor(D, a...)

function PEPOTensor(bond::S, phys::S, x::Real, y::Real, z::Real) where {S<:IndexSpace}
    return PEPOTensor(bloch_pepo(bond, phys, x, y, z))
end

function PEPOTensor(bond::S, phys::S, a::Vector{<:Real}) where {S<:IndexSpace}
    return PEPOTensor(bloch_pepo(bond, phys, a...))
end

"""
    struct FinitePEPO{S<:IndexSpace} <: AbstractPEPO{S}

Concrete subtype of `AbstractPEPO` for representing finite sized square lattices of
single-site PEPO tensors.
"""
struct FinitePEPO{S<:IndexSpace} <: AbstractPEPO{S}
    sites::Array{PEPOTensor{S},2}
    function FinitePEPO(sites::Array{PEPOTensor{S},2}) where {S<:IndexSpace}
        _check_bond_axis(sites)
        _check_bond_axis(sites; permute=true)
        return new{S}(sites)
    end
end

# FIXME: Refactor this function such that the array of tensor does not need to be copied.
# This will also allow the index of sites in the error message to be displayed more
# easily. Solution would be to transpose the indexing rather than transposing the array.
function _check_bond_axis(sites; permute=false)
    if permute
        permutedims(sites)
    end
    N = Base.size(sites)[1] - 1  # Number of rows in lattice
    M = Base.size(sites)[2] - 1  # Number of columns in lattice
    for i in N  # Loop down the rows (from top to bottom)
        for j in M  # Loop along the columns (from left to right)
            if !(_check_bond(sites[i, j], sites[i, j + 1]))
                throw(
                    SpaceMismatch(
                        "Domain of site ($i,$j) and codomain" *
                        " of site ($i,$(j+1)) do not match",
                    ),
                )
            end
        end
    end

    return nothing
end

function _check_bond(site1, site2)
    if domain(site1.tensor)[2] == codomain(site2.tensor)[1]
        return true
    else
        return false
    end
end

"""
    struct InfinitePEPO{S<:IndexSpace} <: AbstractPEPO{S}

Concrete subtype of `AbstractPEPO` for representing infinite sized square lattices
constructed by tessellating a single unit-cell.
"""
struct InfinitePEPO{S<:IndexSpace} <: AbstractPEPO{S}
    sites::Array{PEPOTensor{S},2}
end

function Base.getindex(A::AbstractPEPO, i1::Int, i2::Int)
    return getfield(A, :sites)[i1, i2]
end

function Base.size(A::AbstractPEPO)
    return Base.size(getfield(A, :sites))
end

function bloch_pepo_data(D::Int, a1::Real, a2::Real, a3::Real)
    if norm([a1, a2, a3]) > one(a1)
        throw(
            ArgumentError(
                "Magnitude of Bloch vector is $(norm([a1,a2,a3])). A valid Bloch vector " *
                "must have magnitude ≤ 1.",
            ),
        )
    end

    if D < 1
        throw(
            ArgumentError(
                "Bond dimension D is $(D). The bond dimension must be a positive integer " *
                " greater than or equal to 1",
            ),
        )
    end
    I = complex([1 0; 0 1])
    X = complex([0 1; 1 0])
    Y = complex([1 -im; im 0])
    Z = complex([1 0; 0 -1])
    tens = complex(zeros(D, 2, 2))
    tens[1, :, :] .= 1//2 * (I + a1 * X + a2 * Y + a3 * Z)
    @tensor begin
        A[i j d; k l dp] :=
            deepcopy(tens)[i, d, d1] *
            deepcopy(tens)[j, d1, d2] *
            deepcopy(tens)[k, d2, d3] *
            deepcopy(tens)[l, d3, dp]
    end
    A ./= maximum(abs.(A))
    return A
end

function bloch_pepo(D::Int, a1::Real, a2::Real, a3::Real)
    A = bloch_pepo_data(D, a1, a2, a3)
    return TensorMap(A, ℂ^D * ℂ^D * ℂ^2, ℂ^D * ℂ^D * ℂ^2)
end

bloch_pepo(D::Int, a::Vector{<:Real}) = bloch_pepo(D, a...)

function bloch_pepo(
    bond::S, phys::S, a1::Real, a2::Real, a3::Real
) where {S<:GradedSpace}
    A = bloch_pepo_data(dim(bond), a1, a2, a3)
    return TensorMap(A, bond * bond * phys, bond * bond * phys)
end

function vector_bloch_pepo(D, a1,a2,a3)
    I = Matrix{ComplexF64}([1 0; 0 1])
    X = Matrix{ComplexF64}([0 1; 1 0])
    Y = Matrix{ComplexF64}([1 -im; im 0])
    Z = Matrix{ComplexF64}([1 0; 0 -1])

    blochmat = (Matrix{ComplexF64}(1/2*(I + a1 * X + a2 * Y + a3 * Z)))^4
    blochvec = reshape(blochmat, 4)
    data = zeros(ComplexF64, (D,D,4,D,D))
    data[1,1,:,1,1] = blochvec
    ρ = TensorMap(data, ℂ^D*ℂ^D * ℂ^4, ℂ^D *ℂ^D )
    @assert convert(Array, ρ) == data
    return ρ
end

# TODO: structure source file.

max_stars_repo_name: jack-dunham/PEPOKit

max_issues_repo_name: jack-dunham/PEPOKit

max_forks_repo_name: jack-dunham/PEPOKit


text:
@inline function writeexp(buf, pos, v::T,
    precision=-1, plus=false, space=false, hash=false,
    expchar=UInt8('e'), decchar=UInt8('.'), trimtrailingzeros=false) where {T <: Base.IEEEFloat}
    @assert 0 < pos <= length(buf)
    x = Float64(v)
    neg = signbit(x)
    # special cases
    if x == 0
        if neg
            buf[pos] = UInt8('-')
            pos += 1
        elseif plus
            buf[pos] = UInt8('+')
            pos += 1
        elseif space
            buf[pos] = UInt8(' ')
            pos += 1
        end
        buf[pos] = UInt8('0')
        pos += 1
        if precision > 0
            buf[pos] = decchar
            pos += 1
            for _ = 1:precision
                buf[pos] = UInt8('0')
                pos += 1
            end
        elseif hash
            buf[pos] = decchar
            pos += 1
        end
        buf[pos] = expchar
        buf[pos + 1] = UInt8('+')
        buf[pos + 2] = UInt8('0')
        buf[pos + 3] = UInt8('0')
        return pos + 4
    elseif isnan(x)
        buf[pos] = UInt8('N')
        buf[pos + 1] = UInt8('a')
        buf[pos + 2] = UInt8('N')
        return pos + 3
    elseif !isfinite(x)
        if neg
            buf[pos] = UInt8('-')
            pos += 1
        elseif plus
            buf[pos] = UInt8('+')
            pos += 1
        elseif space
            buf[pos] = UInt8(' ')
            pos += 1
        end
        buf[pos] = UInt8('I')
        buf[pos + 1] = UInt8('n')
        buf[pos + 2] = UInt8('f')
        return pos + 3
    end

    bits = Core.bitcast(UInt64, x)
    mant = bits & MANTISSA_MASK
    exp = Int((bits >> 52) & EXP_MASK)

    if exp == 0
        e2 = 1 - 1023 - 52
        m2 = mant
    else
        e2 = exp - 1023 - 52
        m2 = (Int64(1) << 52) | mant
    end
    nonzero = false
    precision += 1
    if neg
        buf[pos] = UInt8('-')
        pos += 1
    elseif plus
        buf[pos] = UInt8('+')
        pos += 1
    elseif space
        buf[pos] = UInt8(' ')
        pos += 1
    end
    digits = 0
    printedDigits = 0
    availableDigits = 0
    e = 0
    if e2 >= -52
        idx = e2 < 0 ? 0 : indexforexp(e2)
        p10bits = pow10bitsforindex(idx)
        len = lengthforindex(idx)
        i = len - 1
        while i >= 0
            j = p10bits - e2
            #=@inbounds=# mula, mulb, mulc = POW10_SPLIT[POW10_OFFSET[idx + 1] + i + 1]
            digits = mulshiftmod1e9(m2 << 8, mula, mulb, mulc, j + 8)
            if printedDigits != 0
                if printedDigits + 9 > precision
                    availableDigits = 9
                    break
                end
                pos = append_nine_digits(digits, buf, pos)
                printedDigits += 9
            elseif digits != 0
                availableDigits = decimallength(digits)
                e = i * 9 + availableDigits - 1
                if availableDigits > precision
                    break
                end
                if precision > 1
                    pos = append_d_digits(availableDigits, digits, buf, pos, decchar)
                else
                    buf[pos] = UInt8('0') + digits
                    pos += 1
                    if hash
                        buf[pos] = decchar
                        pos += 1
                    end
                end
                printedDigits = availableDigits
                availableDigits = 0
            end
            i -= 1
        end
    end
    if e2 < 0 && availableDigits == 0
        idx = div(-e2, 16)
        i = MIN_BLOCK_2[idx + 1]
        while i < 200
            j = 120 + (-e2 - 16 * idx)
            p = POW10_OFFSET_2[idx + 1] + i - MIN_BLOCK_2[idx + 1]
            if p >= POW10_OFFSET_2[idx + 2]
                digits = 0
            else
                #=@inbounds=# mula, mulb, mulc = POW10_SPLIT_2[p + 1]
                digits = mulshiftmod1e9(m2 << 8, mula, mulb, mulc, j + 8)
            end
            if printedDigits != 0
                if printedDigits + 9 > precision
                    availableDigits = 9
                    break
                end
                pos = append_nine_digits(digits, buf, pos)
                printedDigits += 9
            elseif digits != 0
                availableDigits = decimallength(digits)
                e = -(i + 1) * 9 + availableDigits - 1
                if availableDigits > precision
                    break
                end
                if precision > 1
                    pos = append_d_digits(availableDigits, digits, buf, pos, decchar)
                else
                    buf[pos] = UInt8('0') + digits
                    pos += 1
                    if hash
                        buf[pos] = decchar
                        pos += 1
                    end
                end
                printedDigits = availableDigits
                availableDigits = 0
            end
            i += 1
        end
    end
    maximum = precision - printedDigits
    if availableDigits == 0
        digits = 0
    end
    lastDigit = 0
    if availableDigits > maximum
        for k = 0:(availableDigits - maximum - 1)
            lastDigit = digits % 10
            digits = div(digits, 10)
        end
    end
    roundUp = 0
    if lastDigit != 5
        roundUp = lastDigit > 5
    else
        rexp = precision - e
        requiredTwos = -e2 - rexp
        trailingZeros = requiredTwos <= 0 ||
            (requiredTwos < 60 && pow2(m2, requiredTwos))
        if rexp < 0
            requiredFives = -rexp
            trailingZeros = trailingZeros & pow5(m2, requiredFives)
        end
        roundUp = trailingZeros ? 2 : 1
    end
    if printedDigits != 0
        if digits == 0
            for _ = 1:maximum
                buf[pos] = UInt8('0')
                pos += 1
            end
        else
            pos = append_c_digits(maximum, digits, buf, pos)
        end
    else
        if precision > 1
            pos = append_d_digits(maximum, digits, buf, pos, decchar)
        else
            buf[pos] = UInt8('0') + digits
            pos += 1
            if hash
                buf[pos] = decchar
                pos += 1
            end
        end
    end
    if roundUp != 0
        roundPos = pos
        while true
            roundPos -= 1
            if roundPos == 0 || buf[roundPos] == UInt8('-')
                buf[roundPos + 1] = UInt8('1')
                e += 1
                break
            end
            c = roundPos > 0 ? buf[roundPos] : 0x00
            if c == decchar
                continue
            elseif c == UInt8('9')
                buf[roundPos] = UInt8('0')
                roundUp = 1
                continue
            else
                if roundUp == 2 && UInt8(c) % 2 == 0
                    break
                end
                buf[roundPos] = c + 1
                break
            end
        end
    end
    if trimtrailingzeros
        while buf[pos - 1] == UInt8('0')
            pos -= 1
        end
        if buf[pos - 1] == decchar && !hash
            pos -= 1
        end
    end
    buf[pos] = expchar
    pos += 1
    if e < 0
        buf[pos] = UInt8('-')
        pos += 1
        e = -e
    else
        buf[pos] = UInt8('+')
        pos += 1
    end
    if e >= 100
        c = e % 10
        unsafe_copyto!(buf, pos, DIGIT_TABLE, 2 * div(e, 10) + 1, 2)
        buf[pos + 2] = UInt8('0') + c
        pos += 3
    else
        unsafe_copyto!(buf, pos, DIGIT_TABLE, 2 * e + 1, 2)
        pos += 2
    end
    return pos
end

ext: jl

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: TimoLarson/julia

max_forks_repo_name: TimoLarson/julia


text:
import Base: one, zero, eps, issubnormal, isnan, isfinite
import Base.convert

one(T::Type{Sigmoid{N, ES, mode}}) where {N, ES, mode} = reinterpret(T, @invertbit)
zero(T::Type{Sigmoid{N, ES, mode}}) where {N, ES, mode} = reinterpret(T, zero(@UInt))

@generated function maxpos(T::Type{Sigmoid{N, ES, mode}}) where {N, ES, mode}
  v = (@signbit) - increment(Sigmoid{N, ES, mode})
  :(reinterpret(T, $v))
end
realmax(T::Type{Sigmoid{N, ES, mode}}) where {N, ES, mode} = maxpos(T)

@generated function minpos(T::Type{Sigmoid{N, ES, mode}}) where {N, ES, mode}
  v = zero(@UInt) + increment(Sigmoid{N, ES, mode})
  :(reinterpret(T, $v))
end
realmin(T::Type{Sigmoid{N, ES, mode}}) where {N, ES, mode} = minpos(T)

@generated function maxneg(T::Type{Sigmoid{N, ES, mode}}) where {N, ES, mode}
  v = zero(@UInt) - increment(Sigmoid{N, ES, mode})
  :(reinterpret(T, $v))
end

@generated function minneg(T::Type{Sigmoid{N, ES, mode}}) where {N, ES, mode}
  v = (@signbit) + increment(Sigmoid{N, ES, mode})
  :(reinterpret(T, $v))
end

@generated function eps(T::Type{Sigmoid{N, ES, mode}}) where {N, ES, mode}
  if mode == :ubit
  else
  end
end

#no subnormals, no NaNs
issubnormal(x::Sigmoid{N, ES, mode}) where {N, ES, mode} = false
isnan(x::Sigmoid{N, ES, mode}) where {N, ES, mode} = false
isfinite(x::Sigmoid{N, ES, mode}) where {N, ES, mode} = (@u(x) != @signbit)
#isfinite{N,ES}(x::Sigmoid{N, ES, :upper}) = false
#isfinite{N,ES}(x::Sigmoid{N, ES, :lower}) = false
#isfinite{N,ES}(x::Sigmoid{N, ES, :cross}) = false

#iszeroinf
iszeroinf(x::Sigmoid) = reinterpret(@UInt, x) & (~@signbit) == 0

#special constant type symbols
struct ∞; end
struct ∞n; end

Base.:-(::Type{∞}) = ∞n
convert(T::Type{Sigmoid{N, ES, mode}}, ::Type{∞}) where {N, ES, mode} = T(Inf)
convert(T::Type{Sigmoid{N, ES, mode}}, ::Type{∞n}) where {N, ES, mode} = T(Inf)

export ∞

Base.typemin(p::Type{S}) where {S <: Sigmoid} = minneg(p)

Base.typemax(p::Type{S}) where {S <: Sigmoid} = maxpos(p)

max_issues_repo_name: MohHizzani/SigmoidNumbers

max_forks_repo_name: MohHizzani/SigmoidNumbers


text:
#=
ControlFlowCamberFormula.jl
Jacob Child
Flow Lab Onboarding
Feb 3, 2022
Background: A Camber line is "a measure of the curvature of an airfoil"
    It is halfway between the upper and lower curves of the airfoil

    Equation
    if X<=p
        zbar = c*(2*p*x-x^2)/(p^2)
    else if x > p
        zbar = c*(1-2*p+2*p*x-x^2)/((1-p)^2)
    else
        println("error")

Pseudocode - ask for the c and p inputs , run an if statement and have A
for loop run through x values, 0:.1:1;
=#

println("Hello! We are going to calculate the the camber of an airfoil \n")
println("\n input c \n")
input1 = readline()
c = parse(Float64, input1)
println("\n input p \n")
input2 = readline()
p = parse(Float64, input2)
#println(p) #to check if it parsed correctly
x = [0:.1:1.1;]    #makes an array stopping at 1 and step size .1
global zbar = x[1:11]
global k = 1
while x[k] <= 1
    if x[k] <= p
        global zbar[k] = c*(2*p*x[k]-x[k]^2)/(p^2)
        global k = k + 1 #iterate

    elseif x[k] > p
        global zbar[k] = c*(1-2*p+2*p*x[k]-x[k]^2)/((1-p)^2)
        global k = k + 1 #iterate

    else
        println("error")
        break
    end

end
println("\n")
println(zbar)

pop!(x) #removes the last value of the x array to allow DataFrames to work
using DataFrames #only works in Julia REPL?
DataFrame()
DataFrame("X" => x, "zbar" => zbar)


text:
seed!(42)
L = 3
D = 10
N = 20
nSamples = 10
b = 5
x = rand(N, D)
y = rand(N)

@testset "NumericalVI" begin
    @test NumericalVI(:quad) isa QuadratureVI
    @test NumericalVI(:mc) isa MCIntegrationVI
    @test_throws ErrorException NumericalVI(:blah)
    @test NumericalSVI(b, :quad) isa QuadratureVI
    @test NumericalSVI(b, :mc) isa MCIntegrationVI
    @test_throws ErrorException NumericalSVI(b, :blah)
    @test repr(NumericalVI(:quad)) == "Numerical Inference by Quadrature"
    @test repr(NumericalVI(:mc)) ==
          "Numerical Inference by Monte Carlo Integration"
end


text:
#-------------------------------------------------------------------------------
# This function creates an instance of SNEASY+FAIR.
#-------------------------------------------------------------------------------

# Load required packages.
using Mimi
using MimiFAIR
using MimiSNEASY
using DataFrames
using CSVFiles

function create_sneasy_fairch4(;rcp_scenario::String="RCP85", start_year::Int=1765, end_year::Int=2300, etminan_ch4_forcing::Bool=true)

    # ---------------------------------------------
    # Load and clean up necessary data.
    # ---------------------------------------------

    # Load RCP emissions and concentration scenario values (RCP options = "RCP26" and "RCP85").
    rcp_emissions      = DataFrame(load(joinpath(@__DIR__, "..", "..", "data", "model_data", rcp_scenario*"_emissions.csv"), skiplines_begin=36))
    rcp_concentrations = DataFrame(load(joinpath(@__DIR__, "..", "..", "data", "model_data", rcp_scenario*"_concentrations.csv"), skiplines_begin=37))
    rcp_forcing        = DataFrame(load(joinpath(@__DIR__, "..", "..", "data", "model_data", rcp_scenario*"_midyear_radforcings.csv"), skiplines_begin=58))

    # Load RCP scenario data for fossil CH₄ fraction as a proportion of total anthropogenic CH₄ emissions.
    rcp_fossil_ch4_frac = DataFrame(load(joinpath(@__DIR__, "..", "..", "data", "model_data", "FAIR_Fossil_CH4_fraction.csv"), skiplines_begin=4))[:,Symbol(rcp_scenario)]

    # Load FAIR greenhouse species data.
    fair_ghg_data = DataFrame(load(joinpath(@__DIR__, "..", "..", "data", "model_data", "FAIR_GHG_species_data.csv"), skiplines_begin=10))

    # Find start and end year indices to crop RCP scenario data to correct model time horizon.
    rcp_indices = findall((in)(collect(start_year:end_year)), rcp_emissions.YEARS)

    # Set pre-industrial atmospheric CO₂, CH₄, and N₂O concentrations to RCP values in 1765.
    CO₂_0 = rcp_concentrations[rcp_concentrations.YEARS .== 1765, :CO2][1]
    CH₄_0 = rcp_concentrations[rcp_concentrations.YEARS .== 1765, :CH4][1]
    N₂O_0 = rcp_concentrations[rcp_concentrations.YEARS .== 1765, :N2O][1]

    # Calculate carbon dioxide emissions as well as aerosol and exogenous RCP radiative forcing scenarios.
    rcp_co2_emissions     = rcp_emissions.FossilCO2 .+ rcp_emissions.OtherCO2
    rcp_aerosol_forcing   = rcp_forcing.TOTAER_DIR_RF .+ rcp_forcing.CLOUD_TOT_RF
    rcp_exogenous_forcing = rcp_forcing.TOTAL_INCLVOLCANIC_RF .- rcp_forcing.CO2_RF .- rcp_forcing.CH4_RF .- rcp_forcing.TROPOZ_RF .- rcp_forcing.CH4OXSTRATH2O_RF .- rcp_aerosol_forcing

    # Calculate CH₄ and CO₂ emission to concentration conversion (based on atmospheric mas of 5.1352e18 kg).
    mol_wt_air     = fair_ghg_data[fair_ghg_data.gas .== "AIR", :mol_weight][1]
    mol_wt_ch4     = fair_ghg_data[fair_ghg_data.gas .== "CH4", :mol_weight][1]
    mol_wt_carbon  = fair_ghg_data[fair_ghg_data.gas .== "C", :mol_weight][1]
    CH₄_emiss2conc = (5.1352e18 / 1.0e18) * (mol_wt_ch4 / mol_wt_air)
    C_emiss2conc   = (5.1352e18 / 1.0e18) * (mol_wt_carbon / mol_wt_air)

    # ------------------------------------------------------------
    # Initialize Mimi-SNEASY and add new CH₄ components to model.
    # ------------------------------------------------------------

    # Get an instance of Mimi-SNEASY.
    m = MimiSNEASY.get_model(start_year=start_year, end_year=end_year)

    # Remove old radiative forcing components.
    delete!(m, :rfco2)
    delete!(m, :radiativeforcing)

    # Add in new components.
    add_comp!(m, rf_total, before = :doeclim)
    add_comp!(m, rf_co2_etminan, before = :rf_total)
    add_comp!(m, MimiFAIR.trop_o3_rf, before = :rf_co2_etminan)
    add_comp!(m, total_co2_emissions, before= :trop_o3_rf)
    add_comp!(m, MimiFAIR.ch4_cycle, before = :total_co2_emissions)

    # Add in user-specified CH₄ radiative forcing component.
    # Note: FAIR uses Etminan et al. CH₄ radiative forcing equations. If not using these, switch to older Myhre et al. equations.
    if etminan_ch4_forcing == true
        add_comp!(m, MimiFAIR.ch4_rf, before = :rf_co2_etminan)
    else
        add_comp!(m, rf_ch4_myhre_fair, before = :rf_co2_etminan)
    end

    # ---------------------------------------------
    # Set component parameters.
    # ---------------------------------------------

    # ---- Common parameters ----
    Mimi.set_external_param!(m, :CH₄_0, CH₄_0)
    Mimi.set_external_param!(m, :N₂O_0, N₂O_0)
    Mimi.set_external_param!(m, :N₂O, rcp_concentrations.N2O[rcp_indices], param_dims=[:time])

    # ---- Carbon Cycle ---- #
    update_param!(m, :atmco20, CO₂_0)

    # ---- Methane Cycle ---- #
    set_param!(m, :ch4_cycle, :fossil_emiss_CH₄, rcp_emissions.CH4[rcp_indices])
    set_param!(m, :ch4_cycle, :natural_emiss_CH₄, ones(length(start_year:end_year)).*190.5807)
    connect_param!(m, :ch4_cycle, :CH₄_0, :CH₄_0)
    set_param!(m, :ch4_cycle, :τ_CH₄, 9.3)
    set_param!(m, :ch4_cycle, :fossil_frac, rcp_fossil_ch4_frac[rcp_indices])
    set_param!(m, :ch4_cycle, :oxidation_frac, 0.61)
    set_param!(m, :ch4_cycle, :mol_weight_CH₄, fair_ghg_data[fair_ghg_data.gas .== "CH4", :mol_weight][1])
    set_param!(m, :ch4_cycle, :mol_weight_C, fair_ghg_data[fair_ghg_data.gas .== "C", :mol_weight][1])
    set_param!(m, :ch4_cycle, :emiss2conc_ch4, CH₄_emiss2conc)
    set_param!(m, :ch4_cycle, :gtc2ppm, C_emiss2conc)

    # ---- Total Carbon Dioxide Emissions ---- #
    set_param!(m, :total_co2_emissions, :exogenous_CO₂_emissions, rcp_co2_emissions[rcp_indices])

    # ---- Methane Radiative Forcing ---- #
    if etminan_ch4_forcing == true
        connect_param!(m, :ch4_rf, :CH₄_0, :CH₄_0)
        connect_param!(m, :ch4_rf, :N₂O_0, :N₂O_0)
        set_param!(m, :ch4_rf, :scale_CH₄, 1.0)
        set_param!(m, :ch4_rf, :a₃, -1.3e-6)
        set_param!(m, :ch4_rf, :b₃, -8.2e-6)
        connect_param!(m, :ch4_rf, :N₂O, :N₂O)
        set_param!(m, :ch4_rf, :H₂O_share, 0.12)
    else
        connect_param!(m, :rf_ch4_myhre_fair, :N₂O_0, :N₂O_0)
        connect_param!(m, :rf_ch4_myhre_fair, :CH₄_0, :CH₄_0)
        set_param!(m, :rf_ch4_myhre_fair, :scale_CH₄, 1.0)
        set_param!(m, :rf_ch4_myhre_fair, :H₂O_share, 0.15)
    end

    # ---- Tropospheric Ozone Radiative Forcing ---- #
    set_param!(m, :trop_o3_rf, :NOx_emissions, rcp_emissions.NOx[rcp_indices])
    set_param!(m, :trop_o3_rf, :CO_emissions, rcp_emissions.CO[rcp_indices])
    set_param!(m, :trop_o3_rf, :NMVOC_emissions, rcp_emissions.NMVOC[rcp_indices])
    set_param!(m, :trop_o3_rf, :mol_weight_N, fair_ghg_data[fair_ghg_data.gas .== "N", :mol_weight][1])
    set_param!(m, :trop_o3_rf, :mol_weight_NO, fair_ghg_data[fair_ghg_data.gas .== "NO", :mol_weight][1])
    connect_param!(m, :trop_o3_rf, :CH₄_0, :CH₄_0)
    set_param!(m, :trop_o3_rf, :T0, 0.0)
    set_param!(m, :trop_o3_rf, :fix_pre1850_RCP, true)

    # ---- Carbon Dioxide Radiative Forcing ---- #
    set_param!(m, :rf_co2_etminan, :a₁, -2.4e-7)
    set_param!(m, :rf_co2_etminan, :b₁, 7.2e-4)
    set_param!(m, :rf_co2_etminan, :c₁, -2.1e-4)
    set_param!(m, :rf_co2_etminan, :CO₂_0, CO₂_0)
    connect_param!(m, :rf_co2_etminan, :N₂O_0, :N₂O_0)
    connect_param!(m, :rf_co2_etminan, :N₂O, :N₂O)
    set_param!(m, :rf_co2_etminan, :rf_scale_CO₂, co2_rf_scale(3.7, CO₂_0, N₂O_0))

    # ---- Total Radiative Forcing ---- #
    set_param!(m, :rf_total, :α, 1.0)
    # TODO It would be more elegant if the `rf_aerosol` parameter in the model didn't exist
    # at this point
    connect_param!(m, :rf_total, :rf_aerosol, :rf_aerosol)
    update_param!(m, :rf_aerosol, rcp_aerosol_forcing[rcp_indices])
    set_param!(m, :rf_total, :rf_exogenous, rcp_exogenous_forcing[rcp_indices])

    # ----------------------------------------------------------
    # Create connections between Mimi SNEASY+Hector components.
    # ----------------------------------------------------------
    connect_param!(m, :doeclim,             :forcing,             :rf_total,            :total_forcing)
    connect_param!(m, :total_co2_emissions, :oxidized_CH₄_to_CO₂, :ch4_cycle,           :oxidised_CH₄_GtC)
    connect_param!(m, :ccm,                 :CO2_emissions,       :total_co2_emissions, :total_CO₂_emissions)
    connect_param!(m, :ccm,                 :temp,                :doeclim,             :temp)
    connect_param!(m, :rf_co2_etminan,      :CO₂,                 :ccm,                 :atmco2)
    connect_param!(m, :trop_o3_rf,          :CH₄,                 :ch4_cycle,           :CH₄)
    connect_param!(m, :trop_o3_rf,          :temperature,         :doeclim,             :temp)

    # Create different connections if using updated Etminan et al. CH₄ forcing equations.
    if etminan_ch4_forcing == true
        connect_param!(m, :ch4_rf,            :CH₄,        :ch4_cycle,         :CH₄)
        connect_param!(m, :rf_total,          :rf_CH₄,     :ch4_rf,            :forcing_CH₄)
        connect_param!(m, :rf_total,          :rf_CH₄_H₂O, :ch4_rf,            :forcing_CH₄_H₂O)

    else
        connect_param!(m, :rf_ch4_myhre_fair, :CH₄,        :ch4_cycle,         :CH₄)
        connect_param!(m, :rf_total,          :rf_CH₄,     :rf_ch4_myhre_fair, :forcing_CH₄)
        connect_param!(m, :rf_total,          :rf_CH₄_H₂O, :rf_ch4_myhre_fair, :forcing_CH₄_H₂O)

    end

    connect_param!(m, :rf_total, :rf_O₃,  :trop_o3_rf,     :forcing_trop_O₃)
    connect_param!(m, :rf_total, :rf_CO₂, :rf_co2_etminan, :rf_CO₂)

    # Return constructed model.
    return m
end


text:
## load in the data
ENV["OMP_NUM_THREADS"] = 16
using PowerSpectra
using Healpix
using PyCall
using DelimitedFiles
nmt = pyimport("pymaster")

nside = 256
lmax = nside2lmax(nside)

m₁ = PowerSpectra.planck256_polmap("100", "hm1")
m₂ = PowerSpectra.planck256_polmap("100", "hm2")
maskT₁ = PowerSpectra.planck256_mask("100", "hm1", :T)
maskP₁ = PowerSpectra.planck256_mask("100", "hm1", :P)
maskT₂ = PowerSpectra.planck256_mask("100", "hm2", :T)
maskP₂ = PowerSpectra.planck256_mask("100", "hm2", :P)
flat_map = HealpixMap{Float64, RingOrder}(ones(nside2npix(nside)) )
scale!(m₁, 1e6)
scale!(m₂, 1e6)

##
f₁_0 = nmt.NmtField(maskT₁.pixels, [m₁.i.pixels])
f₁_2 = nmt.NmtField(maskP₁.pixels, [m₁.q.pixels, m₁.u.pixels])
f₂_0 = nmt.NmtField(maskT₂.pixels, [m₂.i.pixels])
f₂_2 = nmt.NmtField(maskP₂.pixels, [m₂.q.pixels, m₂.u.pixels])

b = nmt.NmtBin.from_nside_linear(nside, 1);

##
# Compute MASTER estimator
# spin-0 x spin-0
cl_00 = nmt.compute_full_master(f₁_0, f₂_0, b)
# spin-0 x spin-2
cl_02 = nmt.compute_full_master(f₁_0, f₂_2, b)
# spin-0 x spin-2
cl_20 = nmt.compute_full_master(f₁_2, f₂_0, b)
# spin-2 x spin-2
cl_22 = nmt.compute_full_master(f₁_2, f₂_2, b)

##
using JLD2
@save "test/data/planck_spec.jld2" {compress=true} cl_00 cl_02 cl_20 cl_22

##
w = nmt.NmtWorkspace()
@time w.compute_coupling_matrix(f₁_0, f₂_0, b)
mcm00 = w.get_coupling_matrix()
@save "test/data/mcm00.jld2" {compress=true} mcm00

# writedlm("test/mcm_TT_diag.txt", diag(w.get_coupling_matrix()[1:lmax, 1:lmax])[3:767])

##
w = nmt.NmtWorkspace()
@time w.compute_coupling_matrix(f₁_2, f₂_2, b)
mcm22 = w.get_coupling_matrix()
@save "test/data/mcm22.jld2" {compress=true} mcm22
# writedlm("test/mcm_EE_diag.txt", diag(w.get_coupling_matrix()[1:4:4*lmax, 1:4:4*lmax])[3:767])

##
w = nmt.NmtWorkspace()
@time w.compute_coupling_matrix(f₁_0, f₂_2, b)
mcm02 = w.get_coupling_matrix()
@save "test/data/mcm02.jld2" {compress=true} mcm02

##
w = nmt.NmtWorkspace()
@time w.compute_coupling_matrix(f₁_2, f₂_0, b)
mcm20 = w.get_coupling_matrix()
@save "test/data/mcm20.jld2" {compress=true} mcm20

# writedlm("test/mcm_TE_diag.txt", diag(w.get_coupling_matrix()[1:2:2*lmax, 1:2:2*lmax])[3:767])

##
using BenchmarkTools
# @btime mcm(:TE, map2alm(maskT₁), map2alm(maskT₂))
@btime mcm((:EE_BB, :EB_BE), map2alm(maskT₁), map2alm(maskT₂))

##
@btime map2alm(maskT₁), map2alm(maskT₂)

##

# flat_beam = SpectralVector(ones(3*nside))
# flat_mask = HealpixMap{Float64, RingOrder}(ones(nside2npix(nside)) )

# mask = readMapFromFITS("test/example_mask_1.fits", 1, Float64)
# clf()
# hp.mollview(mask.pixels, title="Simple Mask")
# gcf()

max_stars_repo_name: xzackli/PlanckCov

max_issues_repo_name: xzackli/PlanckCov


text:
include("16.jl")
function angleofsquare41111(r::Robot)
    angleofsquare1111(r::Robot)
    for side in (Nord,West,Sud,Ost)
        movehole11111(r,side)
        putmarker!(r)
    end

end

function angleofsquare1111(r::Robot)
    while ((!isborder(r,Sud)) || (!isborder(r,Ost)))
        if (!isborder(r,Sud))
            move!(r,Sud)
        end
        if (!isborder(r,Ost))
            move!(r,Ost)
        end
    end
end
function movehole11111(r::Robot,side::HorizonSide)
    while (!isborder(r,side))
        move!(r,side)
    end
end

max_stars_repo_name: bubasooo/proggggggggg

max_issues_repo_name: bubasooo/proggggggggg

max_forks_repo_name: bubasooo/proggggggggg


text:
using BinaryBuilder

# Collection of sources required to build Gumbo
sources = [
    "https://github.com/google/gumbo-parser.git" =>
    "aa91b27b02c0c80c482e24348a457ed7c3c088e0",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd gumbo-parser/
./autogen.sh
./configure --prefix=$prefix --host=$target
make
make install

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    BinaryProvider.Linux(:i686, :glibc),
    BinaryProvider.Linux(:x86_64, :glibc),
    BinaryProvider.Linux(:aarch64, :glibc),
    BinaryProvider.Linux(:armv7l, :glibc),
    BinaryProvider.Linux(:powerpc64le, :glibc),
    BinaryProvider.MacOS(),
    BinaryProvider.Windows(:i686),
    BinaryProvider.Windows(:x86_64)
]

# The products that we will ensure are always built
products(prefix) = Product[
    LibraryProduct(prefix, "libgumbo", :libgumbo)
]

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

]

# Parse out some command-line arguments
BUILD_ARGS = ARGS

# This sets whether we should build verbosely or not
verbose = "--verbose" in BUILD_ARGS
BUILD_ARGS = filter!(x -> x != "--verbose", BUILD_ARGS)

# This flag skips actually building and instead attempts to reconstruct a
# build.jl from a GitHub release page.  Use this to automatically deploy a
# build.jl file even when sharding targets across multiple CI builds.
only_buildjl = "--only-buildjl" in BUILD_ARGS
BUILD_ARGS = filter!(x -> x != "--only-buildjl", BUILD_ARGS)

if !only_buildjl
    # If the user passed in a platform (or a few, comma-separated) on the
    # command-line, use that instead of our default platforms
    if length(BUILD_ARGS) > 0
        platforms = platform_key.(split(BUILD_ARGS[1], ","))
    end
    info("Building for $(join(triplet.(platforms), ", "))")

    # Build the given platforms using the given sources
    autobuild(pwd(), "Gumbo", platforms, sources, script, products;
                                      dependencies=dependencies, verbose=verbose)
else
    # If we're only reconstructing a build.jl file on Travis, grab the information and do it
    if !haskey(ENV, "TRAVIS_REPO_SLUG") || !haskey(ENV, "TRAVIS_TAG")
        error("Must provide repository name and tag through Travis-style environment variables!")
    end
    repo_name = ENV["TRAVIS_REPO_SLUG"]
    tag_name = ENV["TRAVIS_TAG"]
    product_hashes = product_hashes_from_github_release(repo_name, tag_name; verbose=verbose)
    bin_path = "https://github.com/$(repo_name)/releases/download/$(tag_name)"
    dummy_prefix = Prefix(pwd())
    print_buildjl(pwd(), products(dummy_prefix), product_hashes, bin_path)

    if verbose
        info("Writing out the following reconstructed build.jl:")
        print_buildjl(STDOUT, product_hashes; products=products(dummy_prefix), bin_path=bin_path)
    end
end


text:
module RasterDataSources
@doc let
    path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    read(path, String)
end RasterDataSources

using Dates,
      GeoData,
      HTTP,
      Requires,
      URIs,
      ZipFile

abstract type RasterDataSource end
abstract type RasterDataSet end

struct BioClim <: RasterDataSet end
struct Climate <: RasterDataSet end
struct Weather <: RasterDataSet end
struct LandCover <: RasterDataSet end
struct HabitatHeterogeneity <: RasterDataSet end

export WorldClim, CHELSA, EarthEnv, AWAP, ALWB

export BioClim, Climate, Weather, LandCover, HabitatHeterogeneity

export Values, Deciles

export getraster

export geoarray, stack, series

include("shared.jl")
include("worldclim/shared.jl")
include("worldclim/bioclim.jl")
include("worldclim/climate.jl")
include("worldclim/weather.jl")
include("chelsa/bioclim.jl")
include("earthenv/shared.jl")
include("earthenv/landcover.jl")
include("earthenv/habitatheterogeneity.jl")
include("awap/awap.jl")
include("alwb/alwb.jl")

function __init__()
    @require GeoData="9b6fcbb8-86d6-11e9-1ce7-23a6bb139a78" begin
        include("geodata.jl")
    end
end

end # module


text:
#=
mutable struct DspModel
    # Linear constraints in CSR matrix format
    start::Vector{Cint}
    index::Vector{Cint}
    value::Vector{Float64}

    rlbd::Vector{Float64} # Row lower bounds
    rubd::Vector{Float64} # Row upper bounds

    obj::Vector{Float64}       # Linear objective cofficients
    qobjrowindex::Vector{Int}  # Quadratic objective matrix row index
    qobjcolindex::Vector{Int}  # Quadratic objective matrix column index
    qobjvalue::Vector{Float64} # Quadratic objective matrix values
    qobjnum::Int               # ??

    clbd::Vector{Float64} # Column lower bounds
    cubd::Vector{Float64} # Column upper bounds
    ctype::String         # Column types
    cname::Vector{String} # Column names

    nqrows::Int            # Number of quadratic constraints
    linnzcnt::Vector{Int}  # Number of nonzero elements in linear term for each constraint
    quadnzcnt::Vector{Int} # Number of nonzero elements in quadratic terms for each constraint
    rhs
    sense
    linstart
    linind
    linval
    quadstart
    quadrow
    quadcol
    quadval
end
=#

function classifyConstrs(m::SJ.StructuredModel)

    nchildren = length(SJ.getchildren(m))
    dspenv.is_qc = Dict{Int64,Bool}()

    # first-stage quadratic constraints
    dspenv.is_qc[0] = false

    numQc = 0
    numLc = 0

    dspenv.quadConstrs[0] = Dict{Int64,AbstractConstraint}()
    dspenv.linConstrs[0] = Dict{Int64,AbstractConstraint}()
    quadConstrs_temp = Dict{Int64,AbstractConstraint}()

    # partition the constraint set into linear constraints and quadratic constraints
    ## count the number of linear and quadratic constraints and adjust the indices of linear constraints so that its keys = 1:numLc[id]
    for i = 1:length(m.constraints)
        c = m.constraints[i]
        if typeof(c.func) <: GenericAffExpr
            numLc += 1
            dspenv.linConstrs[0][numLc] = c
        elseif (typeof(c.func) <: GenericQuadExpr)
            numQc += 1
            quadConstrs_temp[numQc] = c
        else
            @error("Current version accepts constraints that are either quadratic or affine.")
        end
    end

    if numQc > 0
        dspenv.is_qc[0] = true
        # adjust the indices of quadratic constraints so that its keys = numLc+1:nrows_core
        for i = 1:numQc
            dspenv.quadConstrs[0][i + numLc] = quadConstrs_temp[i]
        end

        m.constraints = merge(dspenv.linConstrs[0], dspenv.quadConstrs[0])
        print("Core numLc: ", numLc, ", numQc: ", numQc, "\n")
    end

    # second-stage quadratic quadratic constraints
    for (id, blk) in SJ.getchildren(m)

        dspenv.is_qc[id] = false

        numQc = 0
        numLc = 0

        dspenv.quadConstrs[id] = Dict{Int64,AbstractConstraint}()
        dspenv.linConstrs[id] = Dict{Int64,AbstractConstraint}()
        quadConstrs_temp = Dict{Int64,AbstractConstraint}()

        # partition the constraint set into linear constraints and quadratic constraints
        ## count the number of linear and quadratic constraints and adjust the indices of linear constraints so that its keys = 1:numLc[id]
        for i = 1:length(blk.constraints)
            c = blk.constraints[i]
            if typeof(c.func) <: GenericAffExpr
                numLc += 1
                dspenv.linConstrs[id][numLc] = c
            elseif (typeof(c.func) <: GenericQuadExpr)
                numQc += 1
                quadConstrs_temp[numQc] = c
            else
                @error("Current version accepts constraints that are either quadratic or affine.")
            end
        end

        if numQc > 0
            dspenv.is_qc[id] = true
            # adjust the indices of quadratic constraints so that its keys = numLc+1:nrows_core
            for i = 1:numQc
                dspenv.quadConstrs[id][i + numLc] = quadConstrs_temp[i]
            end

            blk.constraints = merge(dspenv.linConstrs[id], dspenv.quadConstrs[id])
            # print("Scen ", id, " numLc: ", numLc, ", numQc: ", numQc, "\n")
        end
    end
end

function get_model_data(m::SJ.StructuredModel, id::Int = 0)

    # retrieve classified constraints
    linConstr = dspenv.linConstrs[id]
    quadConstrs = dspenv.quadConstrs[id]

    # Get a column-wise sparse matrix
    start, index, value, rlbd, rubd = get_constraint_matrix(m, linConstr)

    # column information
    clbd = Vector{Float64}(undef, num_variables(m))
    cubd = Vector{Float64}(undef, num_variables(m))
    ctype = ""
    cname = Vector{String}(undef, num_variables(m))
    for i in 1:num_variables(m)
        vref = SJ.StructuredVariableRef(m, i)
        v = m.variables[vref.idx]
        if v.info.integer
            ctype = ctype * "I"
        elseif v.info.binary
            ctype = ctype * "B"
        else
            ctype = ctype * "C"
        end
        if v.info.has_fix
            clbd[vref.idx] = v.info.fixed_value
            cubd[vref.idx] = v.info.fixed_value
        elseif v.info.binary
            clbd[vref.idx] = 0.0
            cubd[vref.idx] = 1.0
        else
            clbd[vref.idx] = v.info.has_lb ? v.info.lower_bound : -Inf
            cubd[vref.idx] = v.info.has_ub ? v.info.upper_bound : Inf
        end
        cname[vref.idx] = m.varnames[vref.idx]
    end

    # objective coefficients
    obj = zeros(num_variables(m))
    if !(objective_function_type(m) <: Real)
        for (v,coef) in objective_function(m).terms
            obj[v.idx] = coef
        end
    end

    if objective_sense(m) == MOI.MAX_SENSE
        dspenv.objective_sense = -1.
        obj .*= -1
    end

    nqrows, linnzcnt, quadnzcnt, rhs, sense, linstart, linind, linval, quadstart, quadrow, quadcol, quadval = get_qc_data(m, quadConstrs)
    return start, index, value, rlbd, rubd, obj, clbd, cubd, ctype, cname, nqrows, linnzcnt, quadnzcnt, rhs, sense, linstart, linind, linval, quadstart, quadrow, quadcol, quadval
end

function get_constraint_matrix(m::SJ.StructuredModel, linConstrs::Dict{Int64,AbstractConstraint})

    is_parent = SJ.getparent(m) == nothing ? true : false

    num_rows = 0 # need to count
    num_cols = num_variables(m)
    if !is_parent
        num_cols += num_variables(SJ.getparent(m))
    end

    # count the number of nonzero elements
    nnz = 0
    for (i,cons) in linConstrs
        nnz += length(cons.func.terms)
        num_rows += 1
    end

    rind = Vector{Int}(undef, nnz)
    cind = Vector{Int}(undef, nnz)
    value = Vector{Float64}(undef, nnz)
    rlbd = Vector{Float64}(undef, num_rows)
    rubd = Vector{Float64}(undef, num_rows)

    pos = 1
    for (i,cons) in linConstrs
        for (var,coef) in cons.func.terms
            rind[pos] = i
            if is_parent
                cind[pos] = var.idx
            elseif JuMP.owner_model(var) == SJ.getparent(m)
                cind[pos] = var.idx
            else
                cind[pos] = var.idx + num_variables(SJ.getparent(m))
            end
            value[pos] = coef
            pos += 1
        end

        # row bounds
        rlbd[i], rubd[i] = row_bounds_from_moi(cons.set)
    end
    @assert(pos-1==nnz)

    # NOTE: DSP takes CSR (row-wise sparse matrix) format.
    # So I simply switch the entries for columns and rows.
    mat = sparse(cind, rind, value, num_cols, num_rows)
    dropzeros!(mat)

    # sparse description
    start = convert(Vector{Cint}, mat.colptr .- 1)
    index = convert(Vector{Cint}, mat.rowval .- 1)
    value = mat.nzval

    return start, index, value, rlbd, rubd
end

function get_qc_data(m::SJ.StructuredModel, quadConstrs::Dict{Int64,AbstractConstraint})

    is_parent = SJ.getparent(m) == nothing ? true : false

    nqrows = length(quadConstrs)

    linnzcnt = Vector{Int}(undef, nqrows)
    quadnzcnt = Vector{Int}(undef, nqrows)

    linstart = Vector{Int}(undef, nqrows)
    quadstart = Vector{Int}(undef, nqrows)

    rhs = Vector{Float64}(undef, nqrows)
    sense = Vector{Int}(undef, nqrows)

    i = 1

    for (key,con) in quadConstrs

        # Remove terms with `0` coefficients.
        drop_zeros!(con.func)

        # affine terms
        linnzcnt[i] = length(con.func.aff.terms)

        if i > 1
            linstart[i] = linstart[i-1] + linnzcnt[i-1]
        else
            linstart[i] = 0
        end

        # print(i, "th linnzcnt: ", linnzcnt[i], "\n")

        # for (var,coef) in con.func.aff.terms
        #     print("lin var: " , var, ", coef: ", coef, "\n")
        # end

        # quad terms
        quadnzcnt[i] = length(con.func.terms);

        if i > 1
            quadstart[i] = quadstart[i-1] + quadnzcnt[i-1]
        else
            quadstart[i] = 0
        end

        # print(i, "th quadnzcnt: ", quadnzcnt[i], "\n")

        # for (var,coef) in con.func.terms
        #     print("quad col: " , var.a, "\n")
        #     print("quad row: " ,var.b, "\n")
        #     print("coef: ", coef, "\n")
        # end

        if typeof(con.set) <: MOI.LessThan
            sense[i] = 'L'
            # sense = sense * "L"
            rhs[i] = con.set.upper
        elseif typeof(con.set) <: MOI.GreaterThan
            sense[i] = 'G'
            # sense = sense * "G"
            rhs[i] = con.set.lower
        elseif typeof(con.set) <: MOI.EqualTo
            @error("quadratic constraints are either 'LessThan' or 'GreaterThan'.")
        elseif typeof(con.set) <: MOI.Interval
            @error("quadratic constraints are either 'LessThan' or 'GreaterThan'.")
        end

        i += 1
    end

    @assert(i == nqrows + 1)

    total_linnzcnt = nqrows == 0 ? 0 : linstart[nqrows] + linnzcnt[nqrows]
    total_quadnzcnt = nqrows == 0 ? 0 : quadstart[nqrows] + quadnzcnt[nqrows]

    linind = Vector{Int}(undef, total_linnzcnt)
    linval = Vector{Float64}(undef, total_linnzcnt)

    quadrow = Vector{Int}(undef, total_quadnzcnt)
    quadcol = Vector{Int}(undef, total_quadnzcnt)
    quadval = Vector{Float64}(undef, total_quadnzcnt)

    linpos = 1
    quadpos = 1
    for (key,con) in quadConstrs
        # affine terms
        for (var,coef) in con.func.aff.terms
            if is_parent
                linind[linpos] = var.idx - 1
            else
                if JuMP.owner_model(var) == SJ.getparent(m)
                    linind[linpos] = var.idx - 1
                else
                    linind[linpos] = var.idx - 1 + num_variables(SJ.getparent(m))
                end
            end
            linval[linpos] = coef
            linpos += 1
            # print("lin var: " , var, ", coef: ", coef, "\n")
        end
        # quad terms
        for (var,coef) in con.func.terms
            if is_parent
                quadrow[quadpos] = var.a.idx - 1
                quadcol[quadpos] = var.b.idx - 1
            else
                if JuMP.owner_model(var.a) == SJ.getparent(m)
                    quadrow[quadpos] = var.a.idx - 1
                else
                    quadrow[quadpos] = var.a.idx - 1 + num_variables(SJ.getparent(m))
                end
                if JuMP.owner_model(var.b) == SJ.getparent(m)
                    quadcol[quadpos] = var.b.idx - 1
                else
                    quadcol[quadpos] = var.b.idx - 1 + num_variables(SJ.getparent(m))
                end
            end
            quadval[quadpos] = coef
            quadpos += 1
            # print("quad col: " , var.a, "\n")
            # print("quad row: " ,var.b, "\n")
            # print("coef: ", coef, "\n")
        end
    end

    if nqrows > 0
        @assert(linpos-1==sum(linnzcnt[k] for k=1:nqrows))
        @assert(quadpos-1==sum(quadnzcnt[k] for k=1:nqrows))
    end

    return nqrows, linnzcnt, quadnzcnt, rhs, sense, linstart, linind, linval, quadstart, quadrow, quadcol, quadval
end

row_bounds_from_moi(set::MOI.LessThan) = (-Inf, set.upper)
row_bounds_from_moi(set::MOI.GreaterThan) = (set.lower, Inf)
row_bounds_from_moi(set::MOI.EqualTo) = (set.value, set.value)
row_bounds_from_moi(set::MOI.Interval) = @error("Interval row bounds are not supported.")


text:
# Used for objectives and solvers where the gradient and Hessian is available/exists
mutable struct TwiceDifferentiableHV{T,TDF,THv,TX} <: AbstractObjective
    f
    fdf
    hv
    F::T
    DF::TDF
    Hv::THv
    x_f::TX
    x_df::TX
    x_hv::TX
    v_hv::TX
    f_calls::Vector{Int}
    df_calls::Vector{Int}
    hv_calls::Vector{Int}
end

# compatibility with old constructor
function TwiceDifferentiableHV(f, fdf, h, x::TX, F::T, G::TG = copy(x), H::THv = copy(x)) where {T, TG, THv, TX}
    x_f, x_df, x_hv, v_hv = x_of_nans(x), x_of_nans(x), x_of_nans(x), x_of_nans(x)
    TwiceDifferentiableHV{T,TG, THv, TX}(f, fdf, h,
                                        copy(F), copy(G), copy(H),
                                        x_f, x_df, x_hv, v_hv,
                                        [0,], [0,], [0,])
end

function TwiceDifferentiableHV(f, fdf, h, x::AbstractVector{T}) where T
    return TwiceDifferentiableHV(f, fdf, h, x, real(zero(T)))
end

# assumption: fg takes (f, g, x)
function TwiceDifferentiableHV(::Nothing, fg, hv, x::AbstractVector{T}) where T
    f   =     x  -> fg(0, nothing, x)
    fg! = (g, x) -> fg(0, g, x)
    return TwiceDifferentiableHV(f, fg!, hv, x)
end

function gradient!!(obj::TwiceDifferentiableHV, x)
    obj.df_calls .+= 1
    copyto!(obj.x_df, x)
    obj.fdf(obj.DF, x)
end

function hv_product!(obj::TwiceDifferentiableHV, x, v)
    if x != obj.x_hv ||  v != obj.v_hv
        hv_product!!(obj, x, v)
    end
    obj.Hv
end
function hv_product!!(obj::TwiceDifferentiableHV, x, v)
    obj.hv_calls .+= 1
    copyto!(obj.x_hv, x)
    copyto!(obj.v_hv, v)
    obj.hv(obj.Hv, x, v)
end
# Deprecate the following?
# Using it makes code non-generic (it requires storage for the result)
hv_product(obj::TwiceDifferentiableHV) = obj.Hv


text:
module LanguageDetect

import Random
import Distributions: Truncated, Normal
import IterTools: partition, chain

# fix random seed
Random.seed!(0)

include("profiles.jl")
include("utils.jl")

global const _languages = Ref{Vector{String}}(Vector{String}())
global const _word_lang_prob = Ref{Dict{String, Vector{Float64}}}(Dict{String, Vector{Float64}}())
global const _ngrams = 3
# gauss distribution
global const _dist = Truncated(Normal(0., 1.), 0., 1.)

struct Language
    language::String
    probability::Float64
end

function init()
    #= Initializer
    =#
    profiles = load_profiles(joinpath(dirname(pathof(LanguageDetect)), "profiles"))
    for (index, profile) in enumerate(profiles)
        push!(_languages[], profile.name)
        for (word, freq) in profile.freq
            if !haskey(_word_lang_prob[], word)
                _word_lang_prob[][word] = zeros(length(profiles))
            end
            N = length(word)
            if 1 <= N <= 3
                _word_lang_prob[][word][index] = 1.0 * freq / profile.n_words[N]
            end
        end
    end
end

function extract_ngram(s:: String, n::Int)
    grams = String[]
    for gram in partition(s, n, 1)
        push!(grams, join(gram))
    end
    return grams
end

function detect(text::String; max_length::Int = 10000, trials::Int = 7,
                alpha_default::Float64 = 0.5, alpha_width::Float64 = 0.05,
                base_freq::Int = 10000, conv_threshold::Float64 = 0.99999,
                max_iters::Int = 1000, proba_threshold::Float64 = 0.1)::Vector{Language}
    #= Detect language
    Args:
      text: string, input text
      max_length: int, specify max text length
      trials: int, hyperparameter
      alpha_default: float, hyperparameter
      alpha_width: float, hyperparameter
      base_freq: int, hyperparameter
      conv_threshold: float, hyperparameter
      max_iters: int, hyperparameter
      proba_threshold: float, hyperparameter
    =#
    text = clean_text(text)
    # truncation
    text = split(text, "")
    text = join(text[begin: min(length(text), max_length)])
    # build ngrams
    ngrams = [gram for gram in chain([extract_ngram(text, i) for i in 1:_ngrams]...)
              if haskey(_word_lang_prob[], gram) && gram != " "]

    @assert length(ngrams) > 0 "No features in text"

    lang_probas = zeros(length(_languages[]))

    for t in 1:trials
        # init probability
        probas = zeros(length(_languages[]))
        fill!(probas, 1.0 / length(_languages[]))
        alpha = alpha_default + rand(_dist) * alpha_width

        for i in 0:max_iters
            weight = alpha / base_freq
            for j in eachindex(probas)
                @inbounds probas[j] *= weight + _word_lang_prob[][ngrams[rand(1:end)]][j]
            end
            if i % 5 == 0
                probas /= sum(probas)
                if maximum(probas) > conv_threshold
                    break
                end
            end
        end

        for j in eachindex(lang_probas)
            @inbounds lang_probas[j] += probas[j] / trials
        end
    end

    # sort probabilities
    result = [Language(lang, proba) for (lang, proba) in zip(_languages[], lang_probas) if proba > proba_threshold]
    return sort(result, by=x -> -x.probability)
end

# initialize
init()

end # module


text:
using BlockArrays, LinearAlgebra, FillArrays, Base64, Test

# avoid fast-paths for view
bview(a, b) = Base.invoke(view, Tuple{AbstractArray,Any}, a, b)

@testset "Array block interface" begin
    @test 1[Block()] == 1

    A = randn(5)
    @test blocksize(A) == (1,)
    @test blocksize(A,1) == 1
    @test blocksizes(A) == ([5],)
    @test blocksizes(A,1) == [5]
    @test A[Block(1)] == A
    view(A,Block(1))[1] = 2
    @test A[1] == 2
    @test_throws BlockBoundsError A[Block(2)]

    A = randn(5,5)
    @test A[Block(1,1)] == A

    @test similar(BlockArray{Float64,1}, Base.OneTo(5)) isa Vector{Float64}
    @test similar(BlockArray{Float64,1}, 5) isa Vector{Float64}
    @test similar(BlockArray{Float64,1}, (5,)) isa Vector{Float64}
end

@testset "Triangular/Symmetric/Hermitian block arrays" begin
    @testset "square blocks" begin
        A = PseudoBlockArray{ComplexF64}(undef, 1:4, 1:4)
        A .= reshape(1:length(A), size(A))
        U = UpperTriangular(A)
        S = Symmetric(A)
        H = Hermitian(A)

        @test blocksize(U) == blocksize(S) == blocksize(A)
        @test blockaxes(U) == blockaxes(S) == blockaxes(A)

        @test U[Block(2,2)] == UpperTriangular(A[2:3,2:3])
        @test U[Block(2,3)] == A[2:3,4:6]
        @test U[Block(3,2)] == zeros(3,2)
        @test S[Block(2,2)] == Symmetric(A[2:3,2:3])
        @test S[Block(2,3)] == A[2:3,4:6]
        @test S[Block(3,2)] == transpose(A[2:3,4:6])
        @test H[Block(2,2)] == Hermitian(A[2:3,2:3])
        @test H[Block(2,3)] == A[2:3,4:6]
        @test H[Block(3,2)] == A[2:3,4:6]'

        if VERSION < v"1.6-"
            @test stringmime("text/plain", UpperTriangular(A)) == "10×10 UpperTriangular{Complex{Float64},PseudoBlockArray{Complex{Float64},2,Array{Complex{Float64},2},$(typeof(axes(A)))}} with indices 1:1:10×1:1:10:\n 1.0+0.0im  │  11.0+0.0im  21.0+0.0im  │  31.0+0.0im  41.0+0.0im  51.0+0.0im  │  61.0+0.0im  71.0+0.0im  81.0+0.0im   91.0+0.0im\n ───────────┼──────────────────────────┼──────────────────────────────────────┼─────────────────────────────────────────────────\n     ⋅      │  12.0+0.0im  22.0+0.0im  │  32.0+0.0im  42.0+0.0im  52.0+0.0im  │  62.0+0.0im  72.0+0.0im  82.0+0.0im   92.0+0.0im\n     ⋅      │       ⋅      23.0+0.0im  │  33.0+0.0im  43.0+0.0im  53.0+0.0im  │  63.0+0.0im  73.0+0.0im  83.0+0.0im   93.0+0.0im\n ───────────┼──────────────────────────┼──────────────────────────────────────┼─────────────────────────────────────────────────\n     ⋅      │       ⋅           ⋅      │  34.0+0.0im  44.0+0.0im  54.0+0.0im  │  64.0+0.0im  74.0+0.0im  84.0+0.0im   94.0+0.0im\n     ⋅      │       ⋅           ⋅      │       ⋅      45.0+0.0im  55.0+0.0im  │  65.0+0.0im  75.0+0.0im  85.0+0.0im   95.0+0.0im\n     ⋅      │       ⋅           ⋅      │       ⋅           ⋅      56.0+0.0im  │  66.0+0.0im  76.0+0.0im  86.0+0.0im   96.0+0.0im\n ───────────┼──────────────────────────┼──────────────────────────────────────┼─────────────────────────────────────────────────\n     ⋅      │       ⋅           ⋅      │       ⋅           ⋅           ⋅      │  67.0+0.0im  77.0+0.0im  87.0+0.0im   97.0+0.0im\n     ⋅      │       ⋅           ⋅      │       ⋅           ⋅           ⋅      │       ⋅      78.0+0.0im  88.0+0.0im   98.0+0.0im\n     ⋅      │       ⋅           ⋅      │       ⋅           ⋅           ⋅      │       ⋅           ⋅      89.0+0.0im   99.0+0.0im\n     ⋅      │       ⋅           ⋅      │       ⋅           ⋅           ⋅      │       ⋅           ⋅           ⋅      100.0+0.0im"
        else
            @test stringmime("text/plain", UpperTriangular(A)) == "10×10 UpperTriangular{ComplexF64, PseudoBlockMatrix{ComplexF64, Matrix{ComplexF64}, $(typeof(axes(A)))}} with indices 1:1:10×1:1:10:\n 1.0+0.0im  │  11.0+0.0im  21.0+0.0im  │  31.0+0.0im  41.0+0.0im  51.0+0.0im  │  61.0+0.0im  71.0+0.0im  81.0+0.0im   91.0+0.0im\n ───────────┼──────────────────────────┼──────────────────────────────────────┼─────────────────────────────────────────────────\n     ⋅      │  12.0+0.0im  22.0+0.0im  │  32.0+0.0im  42.0+0.0im  52.0+0.0im  │  62.0+0.0im  72.0+0.0im  82.0+0.0im   92.0+0.0im\n     ⋅      │       ⋅      23.0+0.0im  │  33.0+0.0im  43.0+0.0im  53.0+0.0im  │  63.0+0.0im  73.0+0.0im  83.0+0.0im   93.0+0.0im\n ───────────┼──────────────────────────┼──────────────────────────────────────┼─────────────────────────────────────────────────\n     ⋅      │       ⋅           ⋅      │  34.0+0.0im  44.0+0.0im  54.0+0.0im  │  64.0+0.0im  74.0+0.0im  84.0+0.0im   94.0+0.0im\n     ⋅      │       ⋅           ⋅      │       ⋅      45.0+0.0im  55.0+0.0im  │  65.0+0.0im  75.0+0.0im  85.0+0.0im   95.0+0.0im\n     ⋅      │       ⋅           ⋅      │       ⋅           ⋅      56.0+0.0im  │  66.0+0.0im  76.0+0.0im  86.0+0.0im   96.0+0.0im\n ───────────┼──────────────────────────┼──────────────────────────────────────┼─────────────────────────────────────────────────\n     ⋅      │       ⋅           ⋅      │       ⋅           ⋅           ⋅      │  67.0+0.0im  77.0+0.0im  87.0+0.0im   97.0+0.0im\n     ⋅      │       ⋅           ⋅      │       ⋅           ⋅           ⋅      │       ⋅      78.0+0.0im  88.0+0.0im   98.0+0.0im\n     ⋅      │       ⋅           ⋅      │       ⋅           ⋅           ⋅      │       ⋅           ⋅      89.0+0.0im   99.0+0.0im\n     ⋅      │       ⋅           ⋅      │       ⋅           ⋅           ⋅      │       ⋅           ⋅           ⋅      100.0+0.0im"
        end

        V = view(A, Block.(1:2), Block.(1:2))
        @test blockisequal(axes(Symmetric(V)), axes(view(A,Block.(1:2), Block.(1:2))))
    end

    @testset "rect blocks" begin
        A = PseudoBlockArray{ComplexF64}(undef, 1:3, fill(3,2))
        A .= randn(6,6) .+ im * randn(6,6)
        U = UpperTriangular(A)
        S = Symmetric(A)
        H = Hermitian(A)

        @test blockisequal(axes(S), (axes(A,2),axes(A,2)))
        @test blockisequal(axes(U), axes(A))
        @test blocksize(S) == (blocksize(S,1),blocksize(S,2)) == (2,2)
        @test blocksize(U) == (blocksize(U,1),blocksize(U,2)) == blocksize(A)
        @test blockaxes(S) == (blockaxes(S,1),blockaxes(S,2)) == (Block.(1:2), Block.(1:2))
        @test blockaxes(U) == (blockaxes(U,1),blockaxes(U,2)) == blockaxes(A)

        @test U[Block(2,2)] == A[Block(2,2)]
        @test U[Block(3,2)] == triu(A[Block(3,2)])
        @test U[Block(3,1)] == zeros(3,3)
        @test S[Block(2,2)] == Symmetric(A[4:6,4:6])
        @test S[Block(1,2)] == A[1:3,4:6]
        @test S[Block(2,1)] == transpose(A[1:3,4:6])
        @test H[Block(2,2)] == Hermitian(A[4:6,4:6])
        @test H[Block(1,2)] == A[1:3,4:6]
        @test H[Block(2,1)] == A[1:3,4:6]'

        @test !BlockArrays.hasmatchingblocks(U)
        @test BlockArrays.hasmatchingblocks(S)
        @test BlockArrays.hasmatchingblocks(H)
    end
end

@testset "Adjoint/Transpose block arrays" begin
    A = PseudoBlockArray{ComplexF64}(undef, (1:4), (2:5))
    A .= randn.() .+ randn.().*im

    @test blocksize(A') == (4,4)
    @test blocksize(Transpose(A)) == (4,4)

    @test A'[Block(2,2)] == A[Block(2,2)]' == A[2:3,3:5]'
    @test transpose(A)[Block(2,2)] == transpose(A[2:3,3:5])
    @test A'[Block(2,3)] == A[Block(3,2)]'
    @test transpose(A)[Block(2,3)] == transpose(A[Block(3,2)])

    @test BlockArray(A') == A'
    @test BlockArray(transpose(A)) == transpose(A)
end

@testset "Diagonal BlockArray" begin
    A = mortar(Diagonal(fill([1 2],2)))
    @test A isa BlockMatrix{Int,Diagonal{Matrix{Int}, Vector{Matrix{Int}}}}
    @test A[Block(1,2)] == [0 0]
    @test_throws BlockBoundsError A[Block(1,3)]
    @test A == [1 2 0 0; 0 0 1 2]

    N = 3
    D = Diagonal(mortar(Fill.(-(0:N)-(0:N).^2, 1:2:2N+1)))
    @test axes(D) isa NTuple{2,BlockedUnitRange}
    @test blockisequal(axes(D,1), axes(parent(D),1))
    @test D == Diagonal(Vector(parent(D)))
    @test MemoryLayout(D) isa BlockArrays.DiagonalLayout{<:BlockArrays.BlockLayout}
end

@testset "non-standard block axes" begin
    A = BlockArray([1 2; 3 4], Fill(1,2),Fill(1,2))
    @test A isa BlockMatrix{Int,Matrix{Matrix{Int}},NTuple{2,BlockedUnitRange{StepRange{Int64,Int64}}}}
    A = BlockArray([1 2; 3 4], Fill(1,2),[1,1])
    @test A isa BlockMatrix{Int,Matrix{Matrix{Int}},Tuple{BlockedUnitRange{StepRange{Int64,Int64}},BlockedUnitRange{Vector{Int}}}}
end

@testset "block Fill" begin
    A = Fill(2,(blockedrange([1,2,2]),))
    @test A[Block(1)] ≡ Fill(2,1)
    @test A[Block.(1:2)] == [2,2,2]
    @test A[Block.(1:2)] isa Fill
    @test 2A ≡ Fill(4,axes(A))

    F = Fill(2, (blockedrange([1,2,2]),blockedrange(1:3)))
    @test F[Block(2,2)] ≡ F[Block(2),Block(2)] ≡ Fill(2, 2,2)
    @test F[Block(2),1:5] ≡ Fill(2, 2,5)
    @test F[1:5,Block(2)] ≡ Fill(2, 5,2)
    @test F[:,Block(2)] ≡ Fill(2, (axes(F,1),Base.OneTo(2)))
    @test F[Block(2),:] ≡ Fill(2, (Base.OneTo(2),axes(F,2)))
    @test F[Block.(1:2),Block.(1:2)] == Fill(2, (blockedrange([1,2]),blockedrange(1:2)))

    O = Ones{Int}((blockedrange([1,2,2]),blockedrange(1:3)))
    @test O[Block(2,2)] ≡ O[Block(2),Block(2)] ≡ Ones{Int}(2,2)
    @test O[Block(2),1:5] ≡ Ones{Int}(2,5)
    @test O[1:5,Block(2)] ≡ Ones{Int}(5,2)
    @test O[:,Block(2)] ≡ Ones{Int}((axes(O,1),Base.OneTo(2)))
    @test O[Block(2),:] ≡ Ones{Int}((Base.OneTo(2),axes(O,2)))
    @test O[Block.(1:2),Block.(1:2)] == Ones{Int}((blockedrange([1,2]),blockedrange(1:2)))

    Z = Zeros{Int}((blockedrange([1,2,2]),blockedrange(1:3)))
    @test Z[Block(2,2)] ≡ Z[Block(2),Block(2)] ≡ Zeros{Int}(2,2)
    @test Z[Block(2),1:5] ≡ Zeros{Int}(2,5)
    @test Z[1:5,Block(2)] ≡ Zeros{Int}(5,2)
    @test Z[:,Block(2)] ≡ Zeros{Int}((axes(Z,1),Base.OneTo(2)))
    @test Z[Block(2),:] ≡ Zeros{Int}((Base.OneTo(2),axes(Z,2)))
    @test Z[Block.(1:2),Block.(1:2)] == Zeros{Int}((blockedrange([1,2]),blockedrange(1:2)))

    B = Eye((blockedrange([1,2]),))
    @test B[Block(2,2)] == Matrix(I,2,2)

    C = Eye((blockedrange([1,2,2]),blockedrange([2,2])))
    @test C[Block(2,2)] == [0 0; 1.0 0]

    U = UpperTriangular(Ones((blockedrange([1,2]),blockedrange([2,1]))))

    if VERSION < v"1.6-"
        @test stringmime("text/plain", A) == "5-element Fill{Int64,1,Tuple{BlockedUnitRange{Array{Int64,1}}}} with indices 1:1:5: entries equal to 2"
        @test stringmime("text/plain", B) == "3×3 Diagonal{Float64,Ones{Float64,1,Tuple{BlockedUnitRange{Array{Int64,1}}}}} with indices 1:1:3×1:1:3"
        @test stringmime("text/plain", U) == "3×3 UpperTriangular{Float64,Ones{Float64,2,Tuple{BlockedUnitRange{Array{Int64,1}},BlockedUnitRange{Array{Int64,1}}}}} with indices 1:1:3×1:1:3:\n 1.0  1.0  │  1.0\n ──────────┼─────\n  ⋅   1.0  │  1.0\n  ⋅    ⋅   │  1.0"
    else
        @test stringmime("text/plain", A) == "5-element Fill{Int64, 1, Tuple{BlockedUnitRange{Vector{Int64}}}} with indices 1:1:5: entries equal to 2"
        @test stringmime("text/plain", B) == "3×3 Diagonal{Float64, Ones{Float64, 1, Tuple{BlockedUnitRange{Vector{Int64}}}}} with indices 1:1:3×1:1:3"
        @test stringmime("text/plain", U) == "3×3 UpperTriangular{Float64, Ones{Float64, 2, Tuple{BlockedUnitRange{Vector{Int64}}, BlockedUnitRange{Vector{Int64}}}}} with indices 1:1:3×1:1:3:\n 1.0  1.0  │  1.0\n ──────────┼─────\n  ⋅   1.0  │  1.0\n  ⋅    ⋅   │  1.0"
    end

    @testset "views" begin
        # This in theory can be dropped because `view` returns the block, but we keep in case needed
        a = BlockArray(randn(6), 1:3)
        fill!(bview(a, Block(2)), 2)
        @test view(a, Block(2)) == [2,2]
        f = mortar([Fill(1,2), Fill(2,3)])
        @test FillArrays.getindex_value(bview(f, Block(2))) == 2
    end
end


text:
module ObjectiveC

using Lazy, MacroTools

# Types & Reflection

import Base: show, convert, unsafe_convert, super, methods
export @sel_str, signature

int2bool(x::Integer) = x != 0

include("primitives.jl")
include("methods.jl")

# Calling Machinery

export @objc, @classes

include("call.jl")
include("syntax.jl")

# Class Creation

include("classes.jl")

# API wrappers
include("foundation.jl")

end

# include("cocoa/cocoa.jl")


text:
@testset "Maps and Reductions" begin
    @testset "sum" begin
        sizes = (3, 4, 7)
        @testset "dims = $dims" for dims in (:, 1)
            fkwargs = (dims=dims,)
            @testset "Array{$N, $T}" for N in eachindex(sizes), T in (Float64, ComplexF64)
                s = sizes[1:N]
                x = randn(T, s...)
                ẋ = randn(T, s...)
                x̄ = randn(T, s...)
                y = sum(x; dims=dims)
                Δy = randn(eltype(y), size(y)...)
                frule_test(sum, (x, ẋ); fkwargs=fkwargs)
                rrule_test(sum, Δy, (x, x̄); fkwargs=fkwargs)
            end
        end
    end  # sum

    @testset "sum abs2" begin
        sizes = (3, 4, 7)
        @testset "dims = $dims" for dims in (:, 1)
            fkwargs = (dims=dims,)
            @testset "Array{$N, $T}" for N in eachindex(sizes), T in (Float64, ComplexF64)
                s = sizes[1:N]
                x, ẋ, x̄ = randn(T, s...), randn(T, s...), randn(T, s...)
                y = sum(abs2, x; dims=dims)
                Δy = randn(eltype(y), size(y)...)
                @testset "frule" begin
                    # can't use frule_test here because it doesn't yet ignore nothing tangents
                    y_ad, ẏ_ad = frule((Zero(), Zero(), ẋ), sum, abs2, x; dims=dims)
                    @test y_ad == y
                    ẏ_fd = jvp(_fdm, z -> sum(abs2, z; dims=dims), (x, ẋ))
                    @test ẏ_ad ≈ ẏ_fd
                end
                @testset "rrule" begin
                    rrule_test(sum, Δy, (abs2, nothing), (x, x̄); fkwargs=fkwargs)
                end
            end
        end
    end  # sum abs2
end


text:
##  Test continuous wavelet transform
##  with time series from predator-prey chemostat experiments
##
##  This file is a documented walk-through to reproduce
##  Extended Data Fig. 1 of the publication (associated
##  to the data set)
##  Blasius et al. (2019) Long-term cyclic persistence in
##  an experimental predator-prey system.
##
##  Author Bernd Blasius (bernd.blasius@gmail.com)

module W

using CSV, Dierckx, DSP, PyPlot, StatsBase

include("../src/cwt.jl")
include("../src/tools.jl")
include("../src/wcs.jl")

function do_wavelet(t,x)
  # calculate wavelet transform and power for the chemostat data
  # x: time series; dt: sampling rate
  dt  = 1           # sampling rate: 1 per day
  pad = 1           # pad the time series with zeroes (recommended)
  nvoice = 100      # this will do 100 sub-octaves per octave
  dj = 1 ./ nvoice
  s0 = 2.0          # this says start at a scale of 2 days
  noctave = 3       # calculate scales of 3 octaves
  j1      = noctave * nvoice
  mother  = "MORLET"

  x = x .- mean(x)    # normalize
  x = x ./ std(x)

  wave,period,scale,coi = wavelet(x,dt,pad,dj,s0,j1,mother)
  power = abs2.(wave)   # wavelet power
  # Global wavelet spectrum & significance levels:
  global_ws = mean(sqrt.(power),dims=2)   # time-average over all times

  wave, period, power, global_ws, scale, coi, nvoice
end

function ExtDataFig1()
    # load-in the data
    df = CSV.read("../data/chemostat_experiments/C1.csv")  # read-in data from C1
    tt  = df[:,1]  # time (days)
    mm  = df[:,2]  # algae (10^6 cells / ml)
    bb  = df[:,3]  # rotifers (animals / ml)
    err = df[:,4]  # egg ratio
    ee  = df[:,5]  # eggs      (eggs / ml)
    dd  = df[:,6]  # dead animals (dead animals / ml)
    ss  = df[:,7]   # nitrogen concentration in external medium (mumol N/l)

    # interpolate to constant sampling rate
    dt = 1   # sampling rate  (1 day)
    t,m  = sample_regular(mm,tt,dt,true)
    t,b  = sample_regular(bb,tt,dt,true)
    t,er = sample_regular(err,tt,dt,true)
    t,e  = sample_regular(ee,tt,dt,true)
    t,d  = sample_regular(dd,tt,dt,true)

    n = length(t)
    tmax = maximum(t)

    # continuous wavelet transform and wavelet power
    # wavalet parameters are already fitted for these chemoostat data
    (wave1, period1, power1, global_ws1, scale, coi, nvoice) = do_wavelet(t,m)
    (wave2, period2, power2, global_ws2, scale, coi, nvoice) = do_wavelet(t,b)
    (wave3, period3, power3, global_ws3, scale, coi, nvoice) = do_wavelet(t,er)
    (wave4, period4, power4, global_ws4, scale, coi, nvoice) = do_wavelet(t,e)
    (wave5, period5, power5, global_ws5, scale, coi, nvoice) = do_wavelet(t,d)

    # calculate wavelet cross spectrum (CWS)
    wcs2 = wave1 .* conj.(wave2)
    wcs3 = wave1 .* conj.(wave3)
    wcs4 = wave1 .* conj.(wave4)
    wcs5 = wave1 .* conj.(wave5)

    # time scale smoothing
    spower1 = smooth(power1,dt,nvoice,period2)
    spower2 = smooth(power2,dt,nvoice,period2)
    spower3 = smooth(power3,dt,nvoice,period2)
    spower4 = smooth(power4,dt,nvoice,period2)
    spower5 = smooth(power5,dt,nvoice,period2)
    swcs2   = smooth(wcs2,dt,nvoice,period2)
    swcs3   = smooth(wcs3,dt,nvoice,period2)
    swcs4   = smooth(wcs4,dt,nvoice,period2)
    swcs5   = smooth(wcs5,dt,nvoice,period2)

    swcs_abs2 = abs.(swcs2)   # absolute value of wavelet cross spectrum
    swcs_abs3 = abs.(swcs3)   # absolute value of wavelet cross spectrum
    swcs_abs4 = abs.(swcs4)   # absolute value of wavelet cross spectrum
    swcs_abs5 = abs.(swcs5)   # absolute value of wavelet cross spectrum
    sphs2 = angle.(swcs2)     # wavelet phase
    sphs3 = angle.(swcs3)     # wavelet phase
    sphs4 = angle.(swcs4)     # wavelet phase
    sphs5 = angle.(swcs5)     # wavelet phase

    # calculate wavelet coherence
    swco2 = swcs_abs2 ./ (sqrt.(spower1) .* sqrt.(spower2))
    swco3 = swcs_abs3 ./ (sqrt.(spower1) .* sqrt.(spower3))
    swco4 = swcs_abs4 ./ (sqrt.(spower1) .* sqrt.(spower4))
    swco5 = swcs_abs5 ./ (sqrt.(spower1) .* sqrt.(spower5))

    # calculate global CWS  ...
    global_wcso = mean(swco2,dims=2)      # time-average over all times
    global_wcs = mean(swcs_abs2,dims=2)   # time-average over all times
    # ... and normalize
    gws1 = global_ws1 / maximum(global_ws1)
    gws2 = global_ws2 / maximum(global_ws2)
    gwcso = global_wcso / maximum(global_wcso)
    gwcs = global_wcs / maximum(global_wcs)

    # prefixed scale band
    # consider all scales within 0.6 octaves from peaks in global spectrum
    delta_ind = trunc(Int, 0.6*nvoice)
    ind_max   = argmax(gwcs)[1]        # peak scale in global WCS
    ind_global1 = ind_max - delta_ind  # indices  ...
    ind_global2 = ind_max + delta_ind
    period_min = period2[ind_global1]  # .. and corresponding scales
    period_max = period2[ind_global2]

    # find strongest co-oscillating component
    # i.e., maximal cross-spectrum at each time instance
    wcs_max_ind2,angle_max2,wco_max2 = get_maxInd(swco2,sphs2,n,ind_global1,ind_global2)
    wcs_max_ind3,angle_max3,wco_max3 = get_maxInd(swco3,sphs3,n,ind_global1,ind_global2)
    wcs_max_ind4,angle_max4,wco_max4 = get_maxInd(swco4,sphs4,n,ind_global1,ind_global2)
    wcs_max_ind5,angle_max5,wco_max5 = get_maxInd(swco5,sphs5,n,ind_global1,ind_global2)

    # Identifiy coherent oscillation regions
    cv_wco = 0.83    # --> 95% significance level according to Monte-Carlo simulations
    indices_wco = findall(x-> x >= cv_wco, wco_max2)

    # line of maximal wcs
    wco_scale = log2.(period2[wcs_max_ind2])

    # remove points outside coi
    relevant = wco_scale .< log2.(coi)
    indices_wco = indices_wco[findall(x-> relevant[x], indices_wco)]

    # remove points where maximal wco is at border of scale range
    indices_wco = indices_wco[findall(x-> (x != ind_global2 && x != ind_global1)
                               , wcs_max_ind2[indices_wco])]

    # find contiguous time windows of correlated oscillations
    iwco = indices_wco
    # find all breaks, ie difference of 2 indices larger than 1 ..
    inds = findall(x->x>1, diff(iwco))
    # .. and then put the indices of oscillation times together
    strt_inds = vcat(iwco[1], iwco[inds .+ 1])
    stop_inds = vcat(iwco[inds .- 1] .+ 1, iwco[end])
    delta_inds = stop_inds .- strt_inds .+ 1

    # calculate mean angles and standard deviation
    phi2, var2, std2 = circstats(pi*angle_max2[indices_wco])
    phi3, var3, std3 = circstats(pi*angle_max3[indices_wco])
    phi4, var4, std4 = circstats(pi*angle_max4[indices_wco])
    phi5, var5, std5 = circstats(pi*angle_max5[indices_wco])
    fac = 180.0/pi
    println("mean angle differences (in degrees)")
    println("angle 2 ", fac*phi2, " std ", fac*std2)
    println("angle 3 ", fac*phi3, " std ", fac*std3)
    println("angle 4 ", fac*phi4, " std ", fac*std4)
    println("angle 5 ", fac*phi5, " std ", fac*std5)

    # calculate histograms of phase differences
    ph_bins  = 0.01  # phase bins = pi/100
    ph_width = 15    # histogram smoothed over 15 bins (=pi/8 = 27 deg)
    h_angle2, h_bins2 = phasehist(angle_max2, indices_wco, ph_bins, ph_width)
    h_angle3, h_bins3 = phasehist(angle_max3, indices_wco, ph_bins, ph_width)
    h_angle4, h_bins4 = phasehist(angle_max4, indices_wco, ph_bins, ph_width)
    h_angle5, h_bins5 = phasehist(angle_max5, indices_wco, ph_bins, ph_width)

    # standardize phase histograms for plotting
    max_hist_value = maximum([maximum(h_angle2), maximum(h_angle3), maximum(h_angle4),
    maximum(h_angle5)])
    h_angle2 /= max_hist_value
    h_angle3 /= max_hist_value
    h_angle4 /= max_hist_value
    h_angle5 /= max_hist_value

    # calculate histogram of relevant period length
    period_hist = fit(Histogram, wco_scale[indices_wco], nbins=12)
    per_hists = period_hist.weights
    per_hists_bins = period_hist.edges[1][2:end]

    # repeat the last histogram value at the end to get a smooth polar diagram
    h_bins2 = collect(h_bins2)
    push!(h_bins2, h_bins2[1]+ 2*pi)
    push!(h_angle2, h_angle2[1])
    h_bins3 = collect(h_bins3)
    push!(h_bins3, h_bins3[1]+ 2*pi)
    push!(h_angle3, h_angle3[1])
    h_bins4 = collect(h_bins4)
    push!(h_bins4, h_bins4[1]+ 2*pi)
    push!(h_angle4, h_angle4[1])
    h_bins5 = collect(h_bins5)
    push!(h_bins5, h_bins5[1]+ 2*pi)
    push!(h_angle5, h_angle5[1])

    # plotting scale range
    lminp = log2(minimum(period1))
    lmaxp = log2(maximum(period1))
    lrange = (trunc(Int, lminp) :  trunc(Int, lmaxp))

    # ********* Plotting ************
    figure("ExtDataFig1",figsize=(10,10))  # 6 columns
    clf()
    figleft_width_x = 0.65
    figrigth_width_x = 0.15
    figleft_x = 0.1
    figright_x = 0.8
    fig_width_y = 0.1
    fig_sep_y   = 0.06
    fig_y6=0.06
    fig_y5=fig_y6 + fig_width_y + fig_sep_y
    fig_y4=fig_y5 + fig_width_y + fig_sep_y
    fig_y3=fig_y4 + fig_width_y + fig_sep_y
    fig_y2=fig_y3 + fig_width_y + fig_sep_y
    fig_y1=fig_y2 + fig_width_y + fig_sep_y

    # ---------------------------------
    # Time series
    ax = PyPlot.axes([figleft_x, fig_y1, figleft_width_x, fig_width_y])
    plot(t,m/maximum(m),"g",linewidth=2)
    ylabel("abundance", labelpad=8, fontsize=10)
    plot(t,b/maximum(b),"r",linewidth=2)
    xlim(0, tmax)
    ylim(0,1)
    yticks([0,0.5,1])
    ax1 = ax.transAxes
    text(-0.04, 1.10, "a", fontsize=16, fontweight="bold", transform = ax1)

    # --------------------------------------------
    # phase plane
    ax = PyPlot.axes([figright_x, fig_y1, figrigth_width_x, fig_width_y])
    time1=90; time2=130
    # do a bit of filtering for the phase-portrait plot
    mfit = preprocess1(t,m,dt)
    bfit = preprocess1(t,b,dt)
    mf = mfit[time1:time2]   #restrict data to time window
    bf = bfit[time1:time2]
    mf = mf .+ abs(minimum(mf))
    bf = bf .+ abs(minimum(bf))
    mf = mf / maximum(mf)
    bf = bf / maximum(bf)
    plot(mf, bf)
    arrow(0.57,0.9,-0.1,-0.02,color="grey", width=0.02)   # thin arrow beside the line
    yticks([0,1])
    xticks([0,1])
    xlim(0,1)
    ylim(0,1)
    xlabel("prey", labelpad = -5, fontsize=10)  # labelpad: place the label a bit closer to axis
    ylabel("predator", labelpad = -2, fontsize=10)
    ax1=ax.transAxes
    text(-0.14,1.10,"b",fontsize=16, fontweight="bold",transform = ax1)

    # --------------------------------------------
    # power spectrum algae
    ax = PyPlot.axes([figleft_x, fig_y2, figleft_width_x, fig_width_y])
    pcolormesh(t, log2.(period1),log2.(spower1), cmap=PyPlot.cm.jet)
    ylim(lmaxp,lminp)
    yticks(lrange, 2 .^ lrange)
    xlim(0,tmax)
    plot(t,ones(length(t))*log2(period_min),"k")    # lower border line
    plot(t,ones(length(t))*log2(period_max),"k")    # upper border line
    plot(t,log2.(coi),"k")
    plot(t,log2.(period2[wcs_max_ind2]),"k")    # line of maximal wcs
    ylabel("period length (d)", labelpad=12, fontsize=10)
    ax1 = ax.transAxes
    text(-0.04, 1.10, "c", fontsize=16, fontweight="bold", transform = ax1)

    # --------------------------------------------
    # global power spectrum
    ax = PyPlot.axes([figright_x, fig_y2, figrigth_width_x, fig_width_y])
    plot(gws1,log2.(period2),"g")
    plot([0.,1], [1,1]*log2.(period_min),"k")
    plot([0.,1], [1,1]*log2.(period_max),"k")
    xlabel("power", fontsize=10)
    ylim(lmaxp,lminp)
    yticks(lrange, 2 .^ lrange)
    xlim(0,1)
    xticks([0, 0.25, 0.5, 0.75, 1])
    ax1=ax.transAxes
    text(-0.14,1.10,"d",fontsize=16, fontweight="bold",transform = ax1)

    # --------------------------------------------
    # power spectrum rotifers
    ax = PyPlot.axes([figleft_x, fig_y3, figleft_width_x, fig_width_y])
    pcolormesh(t, log2.(period1),log2.(spower2), cmap=PyPlot.cm.jet)
    ylim(lmaxp,lminp)
    yticks(lrange, 2 .^ lrange)
    xlim(0,tmax)
    plot(t,ones(length(t))*log2(period_min),"k")    # lower border line
    plot(t,ones(length(t))*log2(period_max),"k")    # upper border line
    plot(t,log2.(coi),"k") #cone-of-influence, anything "below" is dubious
    plot(t,log2.(period2[wcs_max_ind2]),"k")    # line of maximal wcs
    ylabel("period length (d)", labelpad=12, fontsize=10)
    ax1 = ax.transAxes
    text(-0.04, 1.10, "e", fontsize=16, fontweight="bold", transform = ax1)

    # --------------------------------------------
    # global power spectrum
    ax = PyPlot.axes([figright_x, fig_y3, figrigth_width_x, fig_width_y])
    plot(gws2,log2.(period2),"r")
    plot([0.,1], [1,1]*log2.(period_min),"k")
    plot([0.,1], [1,1]*log2.(period_max),"k")
    xlabel("power", fontsize=10)
    ylim(lmaxp,lminp)
    yticks(lrange, 2 .^ lrange)
    xlim(0,1)
    xticks([0, 0.25, 0.5, 0.75, 1])
    ax1=ax.transAxes
    text(-0.14,1.10,"f",fontsize=16, fontweight="bold",transform = ax1)

    # --------------------------------------------
    # wavelet cross spectrum
    ax = PyPlot.axes([figleft_x, fig_y4, figleft_width_x, fig_width_y])
    label = "g"
    pcolormesh(t, log2.(period1),log2.(swcs_abs2), cmap=PyPlot.cm.jet)
    ylim(lmaxp,lminp)
    yticks(lrange, 2 .^ lrange)
    xlim(0,tmax)
    plot(t,ones(length(t))*log2(period_min),"k")    # lower border line
    plot(t,ones(length(t))*log2(period_max),"k")    # upper border line
    plot(t,log2.(coi),"k") #cone-of-influence, anything "below" is dubious
    plot(t,log2.(period2[wcs_max_ind2]),"k")    # line of maximal wcs
    ylabel("period length (d)", labelpad=12, fontsize=10)
    ax1 = ax.transAxes
    text(-0.04, 1.10, label, fontsize=16, fontweight="bold", transform = ax1)

    # --------------------------------------------
    # global power spectrum
    ax = PyPlot.axes([figright_x,fig_y4, figrigth_width_x, fig_width_y])
    label = "h"
    plot(gwcs,log2.(period2),"k")
    plot([0.,1], [1,1]*log2.(period_min),"k")
    plot([0.,1], [1,1]*log2.(period_max),"k")
    xlabel("power", fontsize=10)
    ylim(lmaxp,lminp)
    yticks(lrange, 2 .^ lrange)
    xlim(0,1)
    xticks([0, 0.25, 0.5, 0.75, 1])
    ax1=ax.transAxes
    text(-0.14,1.10,label,fontsize=16, fontweight="bold",transform = ax1)

    # --------------------------------------------
    # plot WCO
    ax = PyPlot.axes([figleft_x, fig_y5, figleft_width_x, fig_width_y])
    label = "i"
    #pcolormesh(t, log2.(period1),swco2, cmap=PyPlot.cm.jet, alpha = 0.5)
    levels = range(0, stop=1, length=200)
    contourf(t, log2.(period1),swco2,levels, cmap=PyPlot.cm.jet )
    contour(t, log2.(period1),swco2,[cv_wco], colors="black")
    plot(t,log2.(coi),"k") #cone-of-influence, anything "below" is dubious
    plot(t,ones(length(t))*log2(period_min),"k")    # lower border line
    plot(t,ones(length(t))*log2(period_max),"k")    # upper border line
    ylim(lmaxp,lminp)
    yticks(lrange, 2 .^ lrange)
    xlim(0,tmax)
    for i=1:length(strt_inds)
      i1 = strt_inds[i]; i2 = stop_inds[i]
      plot(t[i1:i2],wco_scale[i1:i2],"k-",linewidth=3)
    end
    ylabel("period length (d)", labelpad = 12, fontsize=10)
    ax1 = ax.transAxes
    text(-0.04, 1.10, label, fontsize=16, fontweight="bold", transform = ax1)

    # --------------------------------------------
    # plot histogram of relevant period length
    ax = PyPlot.axes([figright_x,fig_y5, figrigth_width_x, fig_width_y])
    label = "j"
    per_hists = per_hists / maximum(per_hists)
    mmm = maximum(per_hists)
    plot([0.,mmm], [1,1]*log2.(period_min),"k")
    plot([0.,mmm], [1,1]*log2.(period_max),"k")
    wd = 0.9*per_hists_bins[3]-per_hists_bins[2]
    barh(per_hists_bins,per_hists,height=wd,align="center",alpha=0.6)
    xlabel("histogram", fontsize=10)
    ylim(lmaxp,lminp)
    yticks(lrange, 2 .^ lrange)
    xlim(0,1)
    xticks([0, 0.25, 0.5, 0.75, 1])
    ax1=ax.transAxes
    text(-0.14,1.10,label,fontsize=16, fontweight="bold",transform = ax1)

    # --------------------------------------------
    ax = PyPlot.axes([figleft_x, fig_y6, figleft_width_x, fig_width_y])
    label = "k"
    plot([0, tmax],[0.0,0.0],color="grey", "--") #, color="gray")
    plot([0, tmax],[0.5,0.5],color="grey", "--") #, color="gray")
    plot([0, tmax],[1.0,1.0],color="grey", "--") #, color="gray")
    plot(t,angle_max2,"--", color="grey")
    plot(t,angle_max3,"--", color="grey")
    plot(t,angle_max4,"--", color="grey")
    plot(t,angle_max5,"--", color="grey")
    for i=1:length(strt_inds)
      i1 = strt_inds[i]; i2 = stop_inds[i]
      plot(t[i1:i2],angle_max2[i1:i2],"r-",linewidth=2)
      plot(t[i1:i2],angle_max3[i1:i2],"b-",linewidth=2)
      plot(t[i1:i2],angle_max4[i1:i2],"k-",linewidth=2)
      plot(t[i1:i2],angle_max5[i1:i2],"y-",linewidth=2)
    end
    xlim(0,tmax)
    yticks([-0.5, 0, 0.5, 1.0, 1.5])
    ylim(-0.5,1.5)
    xlabel("time in d", fontsize=10)
    ylabel(L"phase difference in $\pi$", labelpad = -2, fontsize=10)
    ax1 = ax.transAxes
    text(-0.04, 1.10, label, fontsize=16,  fontweight="bold",transform = ax1)

    # --------------------------------------------
    ax = PyPlot.axes([figright_x,fig_y6, figrigth_width_x, fig_width_y], polar="true") # Create a polar axis
    label = "l"
    plot(h_bins2, h_angle2,"r")
    fill_between(h_bins2,0, h_angle2,color="red",alpha=0.5)
    plot(h_bins3, h_angle3,"b")
    fill_between(h_bins3,0, h_angle3,color="blue",alpha=0.5)
    plot(h_bins4, h_angle4,"k")
    fill_between(h_bins4,0, h_angle4,color="black",alpha=0.5)
    plot(h_bins5, h_angle5,"y")
    fill_between(h_bins5,0, h_angle5,color="yellow",alpha=0.5)
    xticks([0, pi/2, pi, 3*pi/2], ["0", L"$0.5\pi$", L"$\pi$", L"$1.5 \pi$" ])
    ax.tick_params(pad=-5)  # place the tick_labels a bit closer
    yticks([0,25,0.5,0.75],[])
    plot([0,0], [0,1], "g" , linewidth=2)   # reference line for algae
    ax1=ax.transAxes
    text(-0.57,1.10,label,fontsize=16, fontweight="bold",transform = ax1)

    #savefig("ExtDataFig1.png")
    println("done")

end

ExtDataFig1()

end # modul

max_stars_repo_name: berndblasius/WaveletAnalysis

max_issues_repo_name: berndblasius/WaveletAnalysis

max_forks_repo_name: berndblasius/WaveletAnalysis


text:
function custom_generator{T<:Union{AbstractFloat,Integer}}(a::Number, b::Number, type_a::Bool, type_b::Bool,
   var_type::Type{T}, p::Float64 = 1.5)
    if a > b
        error("Generotor received lower bound that is greater then upper bound:
        $a > $b")
    end
	if a == -Inf && b == Inf
               if var_type<:AbstractFloat
                   return (size)->(rand(var_type)-.5)*var_type(p)*size
               elseif var_type<:Integer
                   return (size)->rand(Vector{var_type}(-size:size))
               else
                   error("Generator for type $var_type is not specified.")
               end
	elseif b == Inf && type_a == false  # a<<x
		if var_type<:Integer
			return (size) -> var_type(a)+var_type(100)*rand(Vector{var_type}(1:size))
		elseif var_type<:AbstractFloat
			return (size) -> var_type(a)+rand(var_type)*var_type(p)^(size-1)
		else
			error("Generator for type $var_type is not specified.")
		end
	elseif b == Inf && type_a == true  # a<x
		if var_type<:Integer
			return (size) -> var_type(a)+rand(Vector{var_type}(1:size))
		elseif var_type<:AbstractFloat
			return (size) -> var_type(a)+rand(var_type)*var_type(p)*size
		else
			error("Generator for type $var_type is not specified.")
		end
	elseif a == -Inf && type_b == false  # x<<b
		if var_type<:Integer
			return (size) -> var_type(b)-var_type(100)*rand(Vector{var_type}(1:size))
		elseif var_type<:AbstractFloat
			return (size) -> var_type(b)-rand(var_type)*var_type(p)^(size-1)
		else
			error("Generator for type $var_type is not specified.")
		end
	elseif a == -Inf && type_b == true  # x<b
		if var_type<:Integer
			return (size) -> var_type(b)-rand(Vector{var_type}(1:size))
		elseif var_type<:AbstractFloat
			return (size) -> var_type(b)-rand(var_type)*var_type(p)*size
		else
			error("Generator for type $var_type is not specified.")
		end
	elseif  abs(a)!=Inf && abs(b)!=Inf && var_type<:AbstractFloat
		if (type_a,type_b) == (true,true)  # x<a<b
			return (size)->var_type(a)+(var_type(b)-var_type(a))*rand(var_type)  # choose uniform [a,b]
		elseif (type_a,type_b) == (true,false) # a<x<<b
			return (size)->var_type(a)+(var_type(b)-var_type(a))*rand(var_type)/var_type(p)^(size-1)
		elseif (type_a,type_b) == (false,true) # a<<x<b
			return (size)->var_type(b)-(var_type(b)-var_type(a))*rand(var_type)/var_type(p)^(size-1)
		end
	else
		error("Could not specify generator for $a,$b,$var_type")
	end
end


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Tests related to the gravity models.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] http://icgem.gfz-potsdam.de/home
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

# File: ./src/earth/gravity_models/gravity_models.jl
# ==================================================

# Function compute_g
# ------------------

################################################################################
#                                 TEST RESULTS
################################################################################
#
# For each embedded model, the norm of the gravitational field obtained from
# `compute_g` is compared to the online computation available in [1]. The inputs
# for the online computation were:
#
#   * Functional selection: `gravitation_ell`.
#   * Reference system: `WGS84`.
#   * Height over Ellipsoid: 0 m.
#   * Grid step: 10 deg.
#   * Tide System: Use unmodified model.
#
################################################################################

@testset "Function compute_g" begin
    # Read the inputs.
    tests_out = (readdlm("./test_results/gravitation/EGM96_6ddd3d7d77ec891a29f7c870494a8a2e215adae4256574e5de3a233380fbecee.gdf";
                         skipstart = 34),
                 readdlm("./test_results/gravitation/JGM2_435951880302b0efcb6c8ad16ce89d1840104941a3cfcd2ad0e25831a6c0f136.gdf";
                         skipstart = 34),
                 readdlm("./test_results/gravitation/JGM3_2d444404c0f11e585c91a810332530020c8ef0a1b51697aae6b9cbedcf3d67cc.gdf";
                         skipstart = 34))

    coefs = (load_gravity_model(EGM96()),
             load_gravity_model(JGM2()),
             load_gravity_model(JGM3()))

    # Compare the results.
    for m = 1:length(coefs)
        coef = coefs[m]
        test_out = tests_out[m]

        for k = 1:size(test_out,1)
            # Get latitude and longitude.
            lat = test_out[k,2]*pi/180
            lon = test_out[k,1]*pi/180

            # Get the expected result in m/s² (the online result is in mGal).
            norm_g_online = test_out[k,3]/100000

            # Use the model to compute the gravity using all the coefficients.
            norm_g = norm(compute_g(coef, geodetic_to_ecef(lat, lon, 0)))

            # Compare the results.
            @test norm_g_online ≈ norm_g atol=5e-8 rtol=0
        end
    end
end

# Function compute_U
# ------------------

################################################################################
#                                 TEST RESULTS
################################################################################
#
# For each embedded model, the gravitational potential obtained from `compute_U`
# is compared to the online computation available in [1]. The inputs for the
# online computation were:
#
#   * Functional selection: `potential_ell`.
#   * Reference system: `WGS84`.
#   * Height over Ellipsoid: 0 m.
#   * Grid step: 10 deg.
#   * Tide System: Use unmodified model.
#
################################################################################

@testset "Function compute_U" begin
    # Read the inputs.
    tests_out = (readdlm("./test_results/potential/EGM96_08b822daf0755f2f048a0af2d1afa5dedc5afa5738e991071e7b0d12cf3337ae.gdf";
                         skipstart = 35),
                 readdlm("./test_results/potential/JGM2_7801c122125724d1f5022463d2f87728edca77158ec4af1643168e31ee3ea384.gdf";
                         skipstart = 35),
                 readdlm("./test_results/potential/JGM3_719cfffa3f04b93dce106013c2b745451b656bdd0550afa7f7d4700e03901f88.gdf";
                         skipstart = 35))

    coefs = (load_gravity_model(EGM96()),
             load_gravity_model(JGM2()),
             load_gravity_model(JGM3()))

    # Compare the results.
    for m = 1:length(coefs)
        coef = coefs[m]
        test_out = tests_out[m]

        for k = 1:size(test_out,1)
            # Get latitude and longitude.
            lat = test_out[k,2]*pi/180
            lon = test_out[k,1]*pi/180

            # Get the expected result in m²/s².
            U_online = test_out[k,3]

            # Use the model to compute the gravity using all the coefficients.
            U = compute_U(coef, geodetic_to_ecef(lat, lon, 0))

            # Compare the results.
            @test U_online ≈ U atol=1e-6 rtol=0
        end
    end
end

# Issue #22 - Gravity model maximum degree
# ----------------------------------------

@testset "Issue #22 - Gravity model maximum degree" begin
    v  = (1 .+ rand(3))*R0
    gm = load_gravity_model(EGM96())

    g_0 = compute_g(gm, v, -1)
    g_1 = compute_g(gm, v, 1)
    g_2 = compute_g(gm, v, 2)
    g_3 = compute_g(gm, v, 360)
    g_4 = compute_g(gm, v, 360, 0)

    @test g_0 == g_3
    @test g_1 != g_2
    @test g_1 != g_3
    @test g_2 != g_3
    @test g_3 != g_4

    U_0 = compute_U(gm, v, -1)
    U_1 = compute_U(gm, v, 1)
    U_2 = compute_U(gm, v, 2)
    U_3 = compute_U(gm, v, 360)
    U_4 = compute_U(gm, v, 360, 0)

    @test U_0 == U_3
    @test U_1 != U_2
    @test U_1 != U_3
    @test U_2 != U_3
    @test U_3 != U_4
end


text:
using StaticArrays

struct dealPBN
  trump::Cint
  first::Cint
  currentTrickSuit::SVector{3, Cint}
  currentTrickRank::SVector{3, Cint}
  remainCards::SVector{80, Cchar}
end

struct solvedPlay
  number::Cint
  tricks::SVector{53, Cint}
end
solvedPlay() = solvedPlay(0, zeros(Cint, 53))

struct playTracePbn
  number::Cint
  cards::SVector{106, Cint}
end

struct futureTricks
  nodes::Cint
  cards::Cint
  suit::SVector{13, Cint}
  rank::SVector{13, Cint}
  equals::SVector{13, Cint}
  score::SVector{13, Cint}
end
futureTricks() = futureTricks(0, 0, zeros(Cint, 13), zeros(Cint, 13), zeros(Cint, 13), zeros(Cint, 13))

"""Not yet tested"""
function AnalysePlay()
  remainCards = MVector{80, Cchar}(undef)
  pbn = "N:QJ6.K652.J85.T98 873.J97.AT764.Q4 K5.T83.KQ9.A7652 AT942.AQ4.32.KJ3"
  for (i, c) in enumerate(pbn)
    remainCards[i] = c
  end
  remainCards[length(pbn)+1] = 0
  spades = 0
  north = 0
  deal = dealPBN(spades, north, zeros(Cint, 3), zeros(Cint, 3), remainCards)
  playTrace = playTracePbn(0, zeros(Cint, 106))
  solved = Ref{solvedPlay}(solvedPlay())

  result = ccall((:AnalysePlayPBN, "./libdds"), Cint, (dealPBN, playTracePbn, Ref{solvedPlay}, Cint),
    deal,
    playTrace,
    solved,
    Threads.threadid() - 1
    )
  return (result, futp[])
end

# Returns (solveBoardResult, futureTricks)
function SolveBoardPBN(deal)
  futp = Ref{futureTricks}(futureTricks())

  result = ccall((:SolveBoardPBN, "./libdds"), Cint, (dealPBN, Cint, Cint, Cint, Ref{futureTricks}, Cint),
    deal,
    -1, #target
    2, #solutions
    0, #mode
    futp,
    Threads.threadid() - 1
    )
  return (result, futp[])
end

max_stars_repo_name: btheodorou/Bridge-NN-Research

max_issues_repo_name: btheodorou/Bridge-NN-Research

max_forks_repo_name: btheodorou/Bridge-NN-Research


text:
function uniform_state_distribution(pomdp::POMDP)
    d = DiscreteBelief(n_states(pomdp))
    return d
end


text:
using Krylov
using LinearOperators
using MatrixMarket
# using ProfileView

function residuals(A, b, shifts, x)
  nshifts = size(shifts, 1);
  r = [ (b - A * x[:,i] - shifts[i] * x[:,i]) for i = 1 : nshifts ];
  return r;
end

# mtx = "data/1138bus.mtx";
mtx = "data/bcsstk09.mtx";

A = MatrixMarket.mmread(mtx);
n = size(A, 1);
b = ones(n); b_norm = norm(b);

# Define a linear operator with preallocation.
Ap = zeros(n);
op = LinearOperator(n, n, true, true, p -> A_mul_B!(1.0,  A, p, 0.0, Ap))

# Solve Ax=b.
(x, stats) = cg_lanczos(op, b);
@time (x, stats) = cg_lanczos(op, b);
show(stats);
r = b - A * x;
@printf("Relative residual without shift: %8.1e\n", norm(r)/norm(b));

# Solve (A+αI)x = b sequentially.
shifts = [1, 2, 3, 4];
(x, stats) = cg_lanczos_shift_seq(op, b, shifts, verbose=false);
# @profile (x, stats) = cg_lanczos_shift_seq(op, b, shifts);
@time (x, stats) = cg_lanczos_shift_seq(op, b, shifts, verbose=false);
show(stats);
r = residuals(A, b, shifts, x);
resids = map(norm, r) / b_norm;
@printf("Relative residuals with shifts:\n");
for resid in resids
  @printf(" %8.1e", resid);
end
@printf("\n");

# ProfileView.view()