text:
using ..Atmos
using ..Atmos: recover_thermo_state
using ..Mesh.Topologies
using ..Mesh.Grids
using ..Thermodynamics
using LinearAlgebra

"""
    setup_atmos_core_diagnostics(
        ::AtmosLESConfigType,
        interval::String,
        out_prefix::String;
        writer::AbstractWriter,
        interpol = nothing,
    )

Create the "AtmosLESCore" `DiagnosticsGroup` which contains the following
diagnostic variables, all of which are density-averaged horizontal averages
conditional upon `q_liq > 0 && w > 0` except for `core_frac`:

- core_frac: cloud core fraction
- u_core: cloud core x-velocity
- v_core: cloud core y-velocity
- w_core: cloud core z-velocity
- avg_rho_core: cloud core air density (_not_ density averaged)
- rho_core: cloud core air density
- qt_core: cloud core total specific humidity
- ql_core: cloud core liquid water specific humidity
- thv_core: cloud core virtual potential temperature
- thl_core: cloud core liquid-ice potential temperature
- ei_core: cloud core specific internal energy
- var_u_core: cloud core variance of x-velocity
- var_v_core: cloud core variance of y-velocity
- var_w_core: cloud core variance of z-velocity
- var_qt_core: cloud core variance of total specific humidity
- var_thl_core: cloud core variance of liquid-ice potential temperature
- var_ei_core: cloud core variance of specific internal energy
- cov_w_rho_core: cloud core vertical eddy flux of density
- cov_w_qt_core: cloud core vertical eddy flux of specific humidity
- cov_w_thl_core: cloud core vertical eddy flux of liquid-ice potential temperature
- cov_w_ei_core: cloud core vertical eddy flux of specific internal energy
- cov_qt_thl_core: cloud core covariance of total specific humidity and liquid-ice potential temperature
- cov_qt_ei_core: cloud core covariance of total specific humidity and specific internal energy

All these variables are output with the `z` dimension (`x3id`) on the DG grid
(`interpol` may _not_ be specified) as well as a (unlimited) `time` dimension
at the specified `interval`.
"""
function setup_atmos_core_diagnostics(
    ::AtmosLESConfigType,
    interval::String,
    out_prefix::String;
    writer = NetCDFWriter(),
    interpol = nothing,
)
    # TODO: remove this
    @assert isnothing(interpol)

    return DiagnosticsGroup(
        "AtmosLESCore",
        Diagnostics.atmos_les_core_init,
        Diagnostics.atmos_les_core_fini,
        Diagnostics.atmos_les_core_collect,
        interval,
        out_prefix,
        writer,
        interpol,
    )
end

# Simple horizontal averages
function vars_atmos_les_core_simple(m::AtmosModel, FT)
    @vars begin
        u_core::FT
        v_core::FT
        w_core::FT
        avg_rho_core::FT        # ρ
        rho_core::FT            # ρρ
        qt_core::FT             # q_tot
        ql_core::FT             # q_liq
        thv_core::FT            # θ_vir
        thl_core::FT            # θ_liq
        ei_core::FT             # e_int
    end
end
num_atmos_les_core_simple_vars(m, FT) =
    varsize(vars_atmos_les_core_simple(m, FT))
atmos_les_core_simple_vars(m, array) =
    Vars{vars_atmos_les_core_simple(m, eltype(array))}(array)

function atmos_les_core_simple_sums!(atmos::AtmosModel, state, thermo, MH, sums)
    sums.u_core += MH * state.ρu[1]
    sums.v_core += MH * state.ρu[2]
    sums.w_core += MH * state.ρu[3]
    sums.avg_rho_core += MH * state.ρ
    sums.rho_core += MH * state.ρ * state.ρ
    sums.qt_core += MH * state.moisture.ρq_tot
    sums.ql_core += MH * thermo.moisture.q_liq * state.ρ
    sums.thv_core += MH * thermo.moisture.θ_vir * state.ρ
    sums.thl_core += MH * thermo.moisture.θ_liq_ice * state.ρ
    sums.ei_core += MH * thermo.e_int * state.ρ

    return nothing
end

# Variances and covariances
function vars_atmos_les_core_ho(m::AtmosModel, FT)
    @vars begin
        var_u_core::FT          # u′u′
        var_v_core::FT          # v′v′
        var_w_core::FT          # w′w′
        var_qt_core::FT         # q_tot′q_tot′
        var_thl_core::FT        # θ_liq_ice′θ_liq_ice′
        var_ei_core::FT         # e_int′e_int′

        cov_w_rho_core::FT      # w′ρ′
        cov_w_qt_core::FT       # w′q_tot′
        cov_w_thl_core::FT      # w′θ_liq_ice′
        cov_w_ei_core::FT       # w′e_int′
        cov_qt_thl_core::FT     # q_tot′θ_liq_ice′
        cov_qt_ei_core::FT      # q_tot′e_int′
    end
end
num_atmos_les_core_ho_vars(m, FT) = varsize(vars_atmos_les_core_ho(m, FT))
atmos_les_core_ho_vars(m, array) =
    Vars{vars_atmos_les_core_ho(m, eltype(array))}(array)

function atmos_les_core_ho_sums!(atmos::AtmosModel, state, thermo, MH, ha, sums)
    u = state.ρu[1] / state.ρ
    u′ = u - ha.u_core
    v = state.ρu[2] / state.ρ
    v′ = v - ha.v_core
    w = state.ρu[3] / state.ρ
    w′ = w - ha.w_core
    q_tot = state.moisture.ρq_tot / state.ρ
    q_tot′ = q_tot - ha.qt_core
    θ_liq_ice′ = thermo.moisture.θ_liq_ice - ha.thl_core
    e_int′ = thermo.e_int - ha.ei_core

    sums.var_u_core += MH * u′^2 * state.ρ
    sums.var_v_core += MH * v′^2 * state.ρ
    sums.var_w_core += MH * w′^2 * state.ρ
    sums.var_qt_core += MH * q_tot′^2 * state.ρ
    sums.var_thl_core += MH * θ_liq_ice′^2 * state.ρ
    sums.var_ei_core += MH * e_int′^2 * state.ρ

    sums.cov_w_rho_core += MH * w′ * (state.ρ - ha.avg_rho_core) * state.ρ
    sums.cov_w_qt_core += MH * w′ * q_tot′ * state.ρ
    sums.cov_w_thl_core += MH * w′ * θ_liq_ice′ * state.ρ
    sums.cov_qt_thl_core += MH * q_tot′ * θ_liq_ice′ * state.ρ
    sums.cov_qt_ei_core += MH * q_tot′ * e_int′ * state.ρ
    sums.cov_w_ei_core += MH * w′ * e_int′ * state.ρ

    return nothing
end

"""
    atmos_les_core_init(dgngrp, currtime)

Initialize the 'AtmosLESCore' diagnostics group.
"""
function atmos_les_core_init(dgngrp::DiagnosticsGroup, currtime)
    atmos = Settings.dg.balance_law
    FT = eltype(Settings.Q)
    mpicomm = Settings.mpicomm
    mpirank = MPI.Comm_rank(mpicomm)

    # FIXME properly
    if !isa(moisture_model(atmos), EquilMoist)
        @warn """
            Diagnostics $(dgngrp.name): can only be used with the `EquilMoist` moisture model
            """
        return nothing
    end

    atmos_collect_onetime(Settings.mpicomm, Settings.dg, Settings.Q)

    if mpirank == 0
        dims = OrderedDict("z" => (AtmosCollected.zvals, Dict()))

        # set up the variables we're going to be writing
        vars = OrderedDict()
        vars["core_frac"] = (("z",), FT, Dict())

        varnames = map(
            s -> startswith(s, "moisture.") ? s[10:end] : s,
            flattenednames(vars_atmos_les_core_simple(atmos, FT)),
        )
        ho_varnames = map(
            s -> startswith(s, "moisture.") ? s[10:end] : s,
            flattenednames(vars_atmos_les_core_ho(atmos, FT)),
        )
        append!(varnames, ho_varnames)
        for varname in varnames
            var = Variables[varname]
            vars[varname] = (("z",), FT, var.attrib)
        end

        # create the output file
        dprefix = @sprintf("%s_%s", dgngrp.out_prefix, dgngrp.name)
        dfilename = joinpath(Settings.output_dir, dprefix)
        noov = Settings.no_overwrite
        init_data(dgngrp.writer, dfilename, noov, dims, vars)
    end

    return nothing
end

"""
    atmos_les_core_collect(dgngrp, currtime)

Perform a global grid traversal to compute various diagnostics.
"""
function atmos_les_core_collect(dgngrp::DiagnosticsGroup, currtime)
    mpicomm = Settings.mpicomm
    dg = Settings.dg
    Q = Settings.Q
    mpirank = MPI.Comm_rank(mpicomm)
    atmos = dg.balance_law
    if !isa(moisture_model(atmos), EquilMoist)
        @warn """
            Diagnostics $(dgngrp.name): can only be used with the `EquilMoist` moisture model
            """
        return nothing
    end
    grid = dg.grid
    grid_info = basic_grid_info(dg)
    topl_info = basic_topology_info(grid.topology)
    Nqk = grid_info.Nqk
    Nqh = grid_info.Nqh
    npoints = prod(grid_info.Nq)
    nrealelem = topl_info.nrealelem
    nvertelem = topl_info.nvertelem
    nhorzelem = topl_info.nhorzrealelem

    # get needed arrays onto the CPU
    if array_device(Q) isa CPU
        state_data = Q.realdata
        aux_data = dg.state_auxiliary.realdata
        vgeo = grid.vgeo
    else
        state_data = Array(Q.realdata)
        aux_data = Array(dg.state_auxiliary.realdata)
        vgeo = Array(grid.vgeo)
    end
    FT = eltype(state_data)

    zvals = AtmosCollected.zvals

    # Visit each node of the state variables array and:
    # - generate and store the thermo variables,
    # - if core condition holds (q_liq > 0 && w > 0)
    #   - count that point in the core fraction for that z
    #   - count the point's weighting towards averaging for that z, and
    #   - accumulate the simple horizontal sums
    #
    core_MH_z = zeros(FT, Nqk * nvertelem)
    thermo_array =
        [zeros(FT, num_thermo(atmos, FT)) for _ in 1:npoints, _ in 1:nrealelem]
    simple_sums = [
        zeros(FT, num_atmos_les_core_simple_vars(atmos, FT))
        for _ in 1:(Nqk * nvertelem)
    ]
    ql_w_gt_0 = [zeros(FT, (Nqh * nhorzelem)) for _ in 1:(Nqk * nvertelem)]
    @traverse_dg_grid grid_info topl_info begin
        state = extract_state(dg, state_data, ijk, e, Prognostic())
        aux = extract_state(dg, aux_data, ijk, e, Auxiliary())
        MH = vgeo[ijk, grid.MHid, e]

        thermo = thermo_vars(atmos, thermo_array[ijk, e])
        compute_thermo!(atmos, state, aux, thermo)

        if thermo.moisture.q_liq > 0 && state.ρu[3] > 0
            idx = (Nqh * (eh - 1)) + (grid_info.Nq[2] * (j - 1)) + i
            ql_w_gt_0[evk][idx] = one(FT)
            core_MH_z[evk] += MH

            simple = atmos_les_core_simple_vars(atmos, simple_sums[evk])
            atmos_les_core_simple_sums!(atmos, state, thermo, MH, simple)
        end
    end

    # reduce horizontal sums and core fraction across ranks and compute averages
    simple_avgs = [
        zeros(FT, num_atmos_les_core_simple_vars(atmos, FT))
        for _ in 1:(Nqk * nvertelem)
    ]
    core_frac = zeros(FT, Nqk * nvertelem)
    MPI.Allreduce!(core_MH_z, +, mpicomm)
    for evk in 1:(Nqk * nvertelem)
        tot_ql_w_gt_0 = MPI.Reduce(sum(ql_w_gt_0[evk]), +, 0, mpicomm)
        tot_horz = MPI.Reduce(length(ql_w_gt_0[evk]), +, 0, mpicomm)

        MPI.Allreduce!(simple_sums[evk], simple_avgs[evk], +, mpicomm)
        simple_avgs[evk] .= simple_avgs[evk] ./ core_MH_z[evk]

        if mpirank == 0
            core_frac[evk] = tot_ql_w_gt_0 / tot_horz
        end
    end

    # complete density averaging
    simple_varnames = flattenednames(vars_atmos_les_core_simple(atmos, FT))
    for vari in 1:length(simple_varnames)
        for evk in 1:(Nqk * nvertelem)
            simple_ha = atmos_les_core_simple_vars(atmos, simple_avgs[evk])
            avg_rho = simple_ha.avg_rho_core
            if simple_varnames[vari] != "avg_rho_core"
                simple_avgs[evk][vari] /= avg_rho
            end
        end
    end

    # compute the variances and covariances
    ho_sums = [
        zeros(FT, num_atmos_les_core_ho_vars(atmos, FT))
        for _ in 1:(Nqk * nvertelem)
    ]
    @traverse_dg_grid grid_info topl_info begin
        state = extract_state(dg, state_data, ijk, e, Prognostic())
        thermo = thermo_vars(atmos, thermo_array[ijk, e])
        MH = vgeo[ijk, grid.MHid, e]

        if thermo.moisture.q_liq > 0 && state.ρu[3] > 0
            simple_ha = atmos_les_core_simple_vars(atmos, simple_avgs[evk])
            ho = atmos_les_core_ho_vars(atmos, ho_sums[evk])
            atmos_les_core_ho_sums!(atmos, state, thermo, MH, simple_ha, ho)
        end
    end

    # reduce across ranks and compute averages
    ho_avgs = [
        zeros(FT, num_atmos_les_core_ho_vars(atmos, FT))
        for _ in 1:(Nqk * nvertelem)
    ]
    for evk in 1:(Nqk * nvertelem)
        MPI.Reduce!(ho_sums[evk], ho_avgs[evk], +, 0, mpicomm)
        if mpirank == 0
            ho_avgs[evk] .= ho_avgs[evk] ./ core_MH_z[evk]
        end
    end

    # complete density averaging and prepare output
    if mpirank == 0
        varvals = OrderedDict()
        varvals["core_frac"] = core_frac

        for (vari, varname) in enumerate(simple_varnames)
            davg = zeros(FT, Nqk * nvertelem)
            for evk in 1:(Nqk * nvertelem)
                davg[evk] = simple_avgs[evk][vari]
            end
            varvals[varname] = davg
        end

        ho_varnames = flattenednames(vars_atmos_les_core_ho(atmos, FT))
        for (vari, varname) in enumerate(ho_varnames)
            davg = zeros(FT, Nqk * nvertelem)
            for evk in 1:(Nqk * nvertelem)
                simple_ha = atmos_les_core_simple_vars(atmos, simple_avgs[evk])
                avg_rho = simple_ha.avg_rho_core
                davg[evk] = ho_avgs[evk][vari] / avg_rho
            end
            varvals[varname] = davg
        end

        # write output
        append_data(dgngrp.writer, varvals, currtime)
    end

    MPI.Barrier(mpicomm)
    return nothing
end # function collect

function atmos_les_core_fini(dgngrp::DiagnosticsGroup, currtime) end


text:
using GLFW
using ModernGL

GLFW.Init()

# set up OpenGL context version
# it seems OpenGL 4.1 is the highest version supported by MacOS.
@static if Sys.isapple()
    const VERSION_MAJOR = 4
    const VERSION_MINOR = 1
end

@static if Sys.isapple()
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, VERSION_MAJOR)
    GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, VERSION_MINOR)
    GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE)
    GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE)
else
    GLFW.DefaultWindowHints()
end

window = GLFW.CreateWindow(640, 480, "Initialization")
@assert window != C_NULL "could not open window with GLFW3."

GLFW.MakeContextCurrent(window)

# dump version info
renderer = unsafe_string(glGetString(GL_RENDERER))
version = unsafe_string(glGetString(GL_VERSION))
println("Renderder: ", renderer)
println("OpenGL version supported: ", version)
@assert VersionNumber(version[1:3]) ≥ VersionNumber(3,2,0) "OpenGL version must ≥ 3.2, Please upgrade your graphic driver."

GLFW.Terminate()

ext: jl

max_stars_repo_path: OpenGL 4 Tutorials/initialization.jl

max_stars_repo_name: Gnimuc/Videre

max_issues_repo_path: OpenGL 4 Tutorials/initialization.jl

max_issues_repo_name: Gnimuc/Videre

max_forks_repo_path: OpenGL 4 Tutorials/initialization.jl

max_forks_repo_name: Gnimuc/Videre


text:
export Codecov

"""
Coverage.Codecov Module

This module provides functionality to push coverage information to the CodeCov.io
web service. It exports the `submit` and `submit_token` methods.
"""

module Codecov
    using Requests
    using Coverage
    using JSON
    using Compat

    export submit, submit_token, submit_local, submit_generic

    #=
    JSON structure for Codecov.io
    https://codecov.io/api#post-report

    {
      "coverage": {
        "path/to/file.py": [null, 1, 0, null, true, 0, 0, 1, 1],
        "path/to/other.py": [null, 0, 1, 1, "1/3", null]
      },
      "messages": {
        "path/to/other.py": {
          "1": "custom message for line 1"
        }
      }
    }
    =#

    # Turn vector of filename : coverage pairs into a dictionary
    function to_json(fcs::Vector{FileCoverage})
        cov = Dict()
        for fc in fcs
            cov[fc.filename] = vcat(nothing, fc.coverage)
        end
        return Dict("coverage" => cov)
    end

    """
    kwargs provides default values to insert into args_array, only if they are
    not already specified in args_array.
    """
    function set_defaults(args_array; kwargs...)
        defined_names = [k for (k,v) in args_array]
        for kwarg in kwargs
            if !(kwarg[1] in defined_names)
                push!(args_array, kwarg)
            end
        end
      return args_array
    end

    """
        submit(fcs::Vector{FileCoverage})

    Takes a vector of file coverage results (produced by `process_folder`),
    and submits them to Codecov.io. Assumes that this code is being run
    on TravisCI. If running locally, use `submit_local`.
    """
    function submit(fcs::Vector{FileCoverage}; kwargs...)
        kwargs = set_defaults(kwargs,
            service      = "travis-org",
            branch       = ENV["TRAVIS_BRANCH"],
            commit       = ENV["TRAVIS_COMMIT"],
            pull_request = ENV["TRAVIS_PULL_REQUEST"],
            job          = ENV["TRAVIS_JOB_ID"],
            slug         = ENV["TRAVIS_REPO_SLUG"],
            build        = ENV["TRAVIS_JOB_NUMBER"],
            )

        submit_generic(fcs; kwargs...)
    end

    import Git

    """
        submit_local(fcs::Vector{FileCoverage})

    Takes a `Vector` of file coverage results (produced by `process_folder`),
    and submits them to Codecov.io. Assumes the submission is being made from
    a local git installation.  A repository token should be specified by a
    'token' keyword argument or the CODECOV_TOKEN environment variable.
    """
    function submit_local(fcs::Vector{FileCoverage}; kwargs...)
        kwargs = set_defaults(kwargs,
            commit = Git.readchomp(`rev-parse HEAD`, dir=""),
            branch = Git.branch(dir="")
            )

        if haskey(ENV, "REPO_TOKEN")
            println("the environment variable REPO_TOKEN is deprecated, use CODECOV_TOKEN instead")
            kwargs = set_defaults(kwargs, token = ENV["REPO_TOKEN"])
        end

        submit_generic(fcs; kwargs...)
    end

    @deprecate submit_token submit_local

    """
        submit_generic(fcs::Vector{FileCoverage})

    Takes a vector of file coverage results (produced by `process_folder`),
    and submits them to a Codecov.io instance. Keyword arguments are converted
    into a generic Codecov.io API uri.  It is essential that the keywords and
    values match the Codecov upload/v2 API specification.
    The `codecov_url` keyword argument or the CODECOV_URL environment variable
    can be used to specify the base path of the uri.
    The `dry_run` keyword can be used to prevent the http request from
    being generated
    """
    function submit_generic(fcs::Vector{FileCoverage}; kwargs...)
        @assert length(kwargs) > 0

        if haskey(ENV, "CODECOV_URL")
            kwargs = set_defaults(kwargs, codecov_url = ENV["CODECOV_URL"])
        end

        if haskey(ENV, "CODECOV_TOKEN")
            kwargs = set_defaults(kwargs, token = ENV["CODECOV_TOKEN"])
        end

        codecov_url = "https://codecov.io"
        dry_run = false
        for (k,v) in kwargs
            if k == :codecov_url
                codecov_url = v
            end
            if k == :dry_run
                dry_run = true
            end
        end
        @assert codecov_url[end] != "/" "the codecov_url should not end with a /, given url $(codecov_url)"

        uri_str = "$(codecov_url)/upload/v2?"
        for (k,v) in kwargs
            if k != :codecov_url && k != :dry_run
                uri_str = "$(uri_str)&$(k)=$(v)"
            end
        end

        println("Codecov.io API URL:")
        println(uri_str)

        if !dry_run
            heads   = Dict("Content-Type" => "application/json")
            data    = to_json(fcs)
            req     = Requests.post(URI(uri_str); json = data, headers = heads)
            println("Result of submission:")
            println(Compat.UTF8String(req.data))
        end
    end

end  # module Codecov


text:
# A simple script to collect benchmark results to a CSV file

#   colwise_bench_results.txt   => colwise_bench_results.csv
#   pairwise_bench_results.txt  => pairwise_bench_results.csv
#

using Compat: String

function parse_results(lines, title1, title2)
    # just a handy function, without a lot of error handling stuff
    # assuming everything is correct

    records = Array((String, Float64, Float64, Float64), 0)

    state = 0

    for raw_line in lines
        line = strip(raw_line)
        if isempty(line) || beginswith(line, "#")
            continue
        end

        if beginswith(line, "bench")
            @assert state == 0
            m = match(r"^bench\s+(\w+)", line)
            name = m.captures[1]
            state = 1

        elseif beginswith(line, title1)
            @assert state == 1
            m = match(r"^(\w+):\s+t\s*=\s*([.\d]+)", line)
            t = m.captures[1]
            @assert t == title1
            v1 = float64(m.captures[2])
            state = 2

        elseif beginswith(line, title2)
            @assert state == 2
            m = match(r"^(\w+):\s+t\s*=\s*([.\d]+)s\s+\|\s*gain\s*=\s*([.\d]+)", line)
            t = m.captures[1]
            @assert t == title2
            v2 = float64(m.captures[2])
            gain = float64(m.captures[3])

            # store the current record
            push!(records, (name, v1, v2, gain))
            state = 0
        end
    end

    return records
end

function collect(title::String, title1, title2)
    println("Processing $title ...")
    infile = string(title, ".txt")
    outfile = string(title, ".csv")

    if !isfile(infile)
        println("File $infile not found -- ignore.")
        return
    end

    fin = open(infile, "r")
    lines = try
        readlines(fin)
    catch err
        close(fin)
        throw(err)
    end

    R = parse_results(lines, title1, title2)

    fout = open(outfile, "w")
    try
        write(fout, "| name | $title1 | $title2 | gain |\n")
        for r in R
            write(fout, "| $(r[1]) | $(r[2]) | $(r[3]) | $(r[4]) |\n")
        end
        flush(fout)
    catch err
        close(fout)
        throw(err)
    end
end

# main

collect("colwise_bench_results", "loop", "colwise")
collect("pairwise_bench_results", "loop", "pairwise")


text:
function test_numbers(::Type{T}) where T <: CasadiSymbolicObject
    @testset "$( string("Conversion of numeric ", T, " to Julia                ") )" begin
        x₁ = rand()
        x₂ = rand(3)
        x₃ = [rand(1,2); Inf -Inf]

        @test to_julia( T( 0) ) == 0.
        @test to_julia( T(-0) ) == 0.
        @test to_julia( T(-Inf) ) == -Inf
        @test to_julia( T( Inf) ) == Inf
        @test to_julia( T(x₁) ) == x₁
        @test to_julia( T(x₂) ) == x₂
        @test to_julia( T(x₃) ) == x₃
    end
end


text:
module boundarycondition

export BoundaryCondition, OpenBC, PeriodicBC

"""
Abstract base class for boundary conditions.
"""
abstract type BoundaryCondition end

"""
Open boundary conditions.
"""
abstract type OpenBC <: BoundaryCondition end

"""
Periodic boundary conditions.
"""
abstract type PeriodicBC <: BoundaryCondition end

end

max_stars_repo_name: AustenLamacraft/ExclusionProcesses

max_issues_repo_name: AustenLamacraft/ExclusionProcesses

max_forks_repo_name: AustenLamacraft/ExclusionProcesses


text:
using DistributionFits
using Test
using Random: Random

@testset "optimize error" begin
    @test_throws Exception DistributionFits.optimize(x -> x*x, DistributionFits.optimizer, -1, 1)
end
# Optim package for interactive testing
i_loadlibs = () -> begin
    push!(LOAD_PATH, expanduser("~/julia/scimltools/")) # access local package repo
end
using Optim: Optim, optimize
@testset "optimize set in __init__ after using Optim" begin
    # set in __init__
    @test DistributionFits.optimizer isa DistributionFits.OptimOptimizer
end

const tests = [
    "fitstats",
    "normal",
    "lognormal",
    "logitnormal",
    "exponential",
    "weibull",
    "gamma",
]
#tests = ["logitnormal"]

for t in tests
    @testset "Test $t" begin
        Random.seed!(345679)
        include("$t.jl")
    end
end

# test coverage of set_optimize (already called  in init)

# print method ambiguities
println("Potentially stale exports: ")
display(Test.detect_ambiguities(DistributionFits))
println()


text:
import Nemo.characteristic, Nemo.gen, Nemo.size
export gen, characteristic, size, elem_to_mat_row!, rand

function gen(R::Union{Generic.ResRing{T},Generic.ResField{T}}) where T<:PolyElem
  return R(gen(base_ring(R)))
end

function gen(R::Union{Generic.ResRing{fq_nmod_poly},Generic.ResField{fq_nmod_poly}}) ## this is not covered by above
  return R(gen(base_ring(R)))              ## and I don't know why
end

function gen(R::Union{Generic.ResRing{nmod_poly},Generic.ResField{nmod_poly}})
  return R(gen(base_ring(R)))
end

function characteristic(R::Union{Generic.ResRing{Nemo.fmpz},Generic.ResField{Nemo.fmpz}})
  return modulus(R)
end

function characteristic(R::Union{Generic.ResRing{nmod_poly},Generic.ResField{nmod_poly}})
  return characteristic(base_ring(base_ring(R)))
end

function characteristic(R::Union{Generic.ResRing{T},Generic.ResField{T}}) where T<:PolyElem
  return characteristic(base_ring(base_ring(R)))
end

# discuss: size = order? order = size?
function size(R::Nemo.Union{Generic.ResRing{Nemo.nmod_poly},Generic.ResField{Nemo.nmod_poly}})
  return characteristic(R)^degree(modulus(R))
end

function size(R::Nemo.Union{Generic.ResRing{T},Generic.ResField{T}}) where T <: ResElem
  return size(base_ring(base_ring(R)))^degree(modulus(R))
end

function size(R::Nemo.Union{Generic.ResRing{fmpz},Generic.ResField{fmpz}})
  return modulus(R)
end

function size(R::Nemo.Union{Generic.ResRing{T},Generic.ResField{T}}) where T<:PolyElem
  return size(base_ring(base_ring(R)))^degree(R.modulus)
end

function size(R::Nemo.Union{Generic.ResRing{fq_nmod_poly},Generic.ResField{fq_nmod_poly}})
  return size(base_ring(base_ring(R)))^degree(R.modulus)
end

function size(R::FqFiniteField)
  return order(R)
end

function size(R::FqNmodFiniteField)
  return order(R)
end

function size(F::GaloisField)
  return order(F)
end

function size(F::Nemo.GaloisFmpzField)
  return order(F)
end

function order(R::Nemo.NmodRing)
  return fmpz(R.n)
end

#################################################
# in triplicate.... and probably cases missing...
function elem_to_mat_row!(M::MatElem, i::Int, a::ResElem{T}) where T <: PolyElem
  z = zero(parent(M[1,1]))
  for j=0:degree(a.data)
    M[i,j+1] = coeff(a.data, j)
  end
  for j=degree(a.data)+2:ncols(M)
    M[i,j] = z
  end
end
function elem_to_mat_row!(M::MatElem, i::Int, a::ResElem{fq_poly})
  z = zero(parent(M[1,1]))
  for j=0:degree(a.data)
    M[i,j+1] = coeff(a.data, j)
  end
  for j=degree(a.data)+2:ncols(M)
    M[i,j] = z
  end
end
function elem_to_mat_row!(M::MatElem, i::Int, a::ResElem{fq_nmod_poly})
  z = zero(parent(M[1,1]))
  for j=0:degree(a.data)
    M[i,j+1] = coeff(a.data, j)
  end
  for j=degree(a.data)+2:ncols(M)
    M[i,j] = z
  end
end

function rand(R::Union{Generic.ResRing{fmpz},Generic.ResField{fmpz}})
  return R(rand(fmpz(0):(size(R)-1)))
end

function rand(R::Generic.ResField{fmpz})
  return R(rand(fmpz(0):(order(R)-1)))
end

function rand(R::Union{Generic.ResRing{T},Generic.ResField{T}}) where T<:PolyElem
  r = rand(base_ring(base_ring(R)))
  g = gen(R)
  for i=1:degree(R.modulus)
    r = r*g + rand(base_ring(base_ring(R)))
  end
  return r
end

function rand(R::Union{Generic.ResRing{fq_nmod_poly},Generic.ResField{fq_nmod_poly}})
  r = rand(base_ring(base_ring(R)))
  g = gen(R)
  for i=1:degree(R.modulus)
    r = r*g + rand(base_ring(base_ring(R)))
  end
  return r
end

function rand(R::Union{Generic.ResRing{fq_poly},Generic.ResField{fq_poly}})
  r = rand(base_ring(base_ring(R)))
  g = gen(R)
  for i=1:degree(R.modulus)
    r = r*g + rand(base_ring(base_ring(R)))
  end
  return r
end

function rand(R::Union{Generic.ResRing{nmod_poly},Generic.ResField{nmod_poly}})
  r = rand(base_ring(base_ring(R)))
  g = gen(R)
  for i=1:degree(R.modulus)
    r = r*g + rand(base_ring(base_ring(R)))
  end
  return r
end

#######################################################
##
##
##
#######################################################

function gens(R::Union{Generic.ResRing{T},Generic.ResField{T}}) where T<:PolyElem ## probably needs more cases
                                          ## as the other residue functions
  g = gen(R)
  r = Vector{typeof(g)}()
  push!(r, one(R))
  if degree(R.modulus)==1
    return r
  end
  push!(r, g)
  for i=2:degree(R.modulus)-1
    push!(r, r[end]*g)
  end
  return r
end

function gens(R::Union{Generic.ResRing{nmod_poly},Generic.ResField{nmod_poly}})
  g = gen(R)
  r = Vector{typeof(g)}()
  push!(r, one(R))
  if degree(R.modulus)==1
    return r
  end
  push!(r, g)
  for i=2:degree(R.modulus)-1
    push!(r, r[end]*g)
  end
  return r
end

function rem!(f::Nemo.nmod_poly, g::Nemo.nmod_poly, h::Nemo.nmod_poly)
  ccall((:nmod_poly_rem, libflint), Nothing, (Ref{Nemo.nmod_poly}, Ref{Nemo.nmod_poly}, Ref{Nemo.nmod_poly}), f, g, h)
  return f
end

function gcd!(f::Nemo.nmod_poly, g::Nemo.nmod_poly, h::Nemo.nmod_poly)
  ccall((:nmod_poly_gcd, libflint), Nothing, (Ref{Nemo.nmod_poly}, Ref{Nemo.nmod_poly}, Ref{Nemo.nmod_poly}), f, g, h)
  return f
end

function gcd!(f::Nemo.gfp_poly, g::Nemo.gfp_poly, h::Nemo.gfp_poly)
  ccall((:nmod_poly_gcd, libflint), Nothing, (Ref{Nemo.gfp_poly}, Ref{Nemo.gfp_poly}, Ref{Nemo.gfp_poly}), f, g, h)
  return f
end

function (R::Nemo.NmodPolyRing)(g::fmpq_poly)
  return fmpq_poly_to_nmod_poly(R, g)
end

function (R::Nemo.GFPPolyRing)(g::fmpq_poly)
  return fmpq_poly_to_gfp_poly(R, g)
end

function (R::Nemo.FmpzModPolyRing)(g::fmpq_poly)
  return fmpq_poly_to_fmpz_mod_poly(R, g)
end

function (R::Nemo.GFPFmpzPolyRing)(g::fmpq_poly)
  return fmpq_poly_to_gfp_fmpz_poly(R, g)
end


text:
#ctqual = "ColorTypes."
ctqual = ""
# fpqual = "FixedPointNumbers."
fpqual = ""

@testset "_default_fill" begin
    @test_throws UndefVarError _default_fill
    @test typeof(ImageTransformations._default_fill) <: Function

    @test @inferred(ImageTransformations._default_fill(N0f8)) === N0f8(0)
    @test @inferred(ImageTransformations._default_fill(Int)) === 0
    @test @inferred(ImageTransformations._default_fill(Float16)) === NaN16
    @test @inferred(ImageTransformations._default_fill(Float32)) === NaN32
    @test @inferred(ImageTransformations._default_fill(Float64)) === NaN

    @test @inferred(ImageTransformations._default_fill(Gray{N0f8})) === Gray{N0f8}(0)
    @test @inferred(ImageTransformations._default_fill(Gray{Float16})) === Gray{Float16}(NaN16)
    @test @inferred(ImageTransformations._default_fill(Gray{Float32})) === Gray{Float32}(NaN32)
    @test @inferred(ImageTransformations._default_fill(Gray{Float64})) === Gray{Float64}(NaN)

    @test @inferred(ImageTransformations._default_fill(GrayA{N0f8})) === GrayA{N0f8}(0,0)
    @test @inferred(ImageTransformations._default_fill(GrayA{Float16})) === GrayA{Float16}(NaN16,NaN16)
    @test @inferred(ImageTransformations._default_fill(GrayA{Float32})) === GrayA{Float32}(NaN32,NaN32)
    @test @inferred(ImageTransformations._default_fill(GrayA{Float64})) === GrayA{Float64}(NaN,NaN)

    @test @inferred(ImageTransformations._default_fill(RGB{N0f8})) === RGB{N0f8}(0,0,0)
    @test @inferred(ImageTransformations._default_fill(RGB{Float16})) === RGB{Float16}(NaN16,NaN16,NaN16)
    @test @inferred(ImageTransformations._default_fill(RGB{Float32})) === RGB{Float32}(NaN32,NaN32,NaN32)
    @test @inferred(ImageTransformations._default_fill(RGB{Float64})) === RGB{Float64}(NaN,NaN,NaN)

    @test @inferred(ImageTransformations._default_fill(RGBA{N0f8})) === RGBA{N0f8}(0,0,0,0)
    @test @inferred(ImageTransformations._default_fill(RGBA{Float16})) === RGBA{Float16}(NaN16,NaN16,NaN16,NaN16)
    @test @inferred(ImageTransformations._default_fill(RGBA{Float32})) === RGBA{Float32}(NaN32,NaN32,NaN32,NaN32)
    @test @inferred(ImageTransformations._default_fill(RGBA{Float64})) === RGBA{Float64}(NaN,NaN,NaN,NaN)

    @test @inferred(ImageTransformations._default_fill(HSV{Float16})) === HSV{Float16}(NaN16,NaN16,NaN16)
    @test @inferred(ImageTransformations._default_fill(HSV{Float32})) === HSV{Float32}(NaN32,NaN32,NaN32)
    @test @inferred(ImageTransformations._default_fill(HSV{Float64})) === HSV{Float64}(NaN,NaN,NaN)
end

@testset "box_extrapolation" begin
    @test_throws UndefVarError box_extrapolation
    @test typeof(ImageTransformations.box_extrapolation) <: Function

    img = rand(Gray{N0f8}, 2, 2)
    n0f8_str = typestring(N0f8)
    matrixf64_str = typestring(Matrix{Float64})

    etp = @inferred ImageTransformations.box_extrapolation(img)
    @test @inferred(ImageTransformations.box_extrapolation(etp)) === etp
    @test summary(etp) == "2×2 extrapolate(interpolate(::Array{Gray{N0f8},2}, BSpline(Linear())), Gray{N0f8}(0.0)) with element type $(ctqual)Gray{$(fpqual)$n0f8_str}"
    @test typeof(etp) <: Interpolations.FilledExtrapolation
    @test etp.fillvalue === Gray{N0f8}(0.0)
    @test etp.itp.coefs === img

    # to catch regressions like #60
    @test @inferred(ImageTransformations._getindex(img, @SVector([1,2]))) isa Gray{N0f8}

    @test_throws ArgumentError ImageTransformations.box_extrapolation(etp, 0)
    @test_throws ArgumentError ImageTransformations.box_extrapolation(etp, Flat())
    @test_throws ArgumentError ImageTransformations.box_extrapolation(etp, Quadratic(Flat()))
    @test_throws ArgumentError ImageTransformations.box_extrapolation(etp, Quadratic(Flat()), Flat())
    @test_throws ArgumentError ImageTransformations.box_extrapolation(etp, Constant())
    @test_throws ArgumentError ImageTransformations.box_extrapolation(etp, Constant(), Flat())
    @test_throws ArgumentError ImageTransformations.box_extrapolation(etp.itp, Constant())
    @test_throws ArgumentError ImageTransformations.box_extrapolation(etp.itp, Constant(), Flat())

    etp2 = @inferred ImageTransformations.box_extrapolation(etp.itp)
    @test summary(etp2) == "2×2 extrapolate(interpolate(::Array{Gray{N0f8},2}, BSpline(Linear())), Gray{N0f8}(0.0)) with element type $(ctqual)Gray{$(fpqual)$n0f8_str}"
    @test typeof(etp2) <: Interpolations.FilledExtrapolation
    @test etp2.fillvalue === Gray{N0f8}(0.0)
    @test etp2 !== etp
    @test etp2.itp === etp.itp

    etp2 = @inferred ImageTransformations.box_extrapolation(etp.itp, Flat())
    @test summary(etp2) == "2×2 extrapolate(interpolate(::Array{Gray{N0f8},2}, BSpline(Linear())), Flat()) with element type $(ctqual)Gray{$(fpqual)$n0f8_str}"
    @test typeof(etp2) <: Interpolations.Extrapolation
    @test etp2 !== etp
    @test etp2.itp === etp.itp

    etp = @inferred ImageTransformations.box_extrapolation(img, 1)
    @test summary(etp) == "2×2 extrapolate(interpolate(::Array{Gray{N0f8},2}, BSpline(Linear())), Gray{N0f8}(1.0)) with element type $(ctqual)Gray{$(fpqual)$n0f8_str}"
    @test typeof(etp) <: Interpolations.FilledExtrapolation
    @test etp.fillvalue === Gray{N0f8}(1.0)
    @test etp.itp.coefs === img

    etp = @inferred ImageTransformations.box_extrapolation(img, Flat())
    @test @inferred(ImageTransformations.box_extrapolation(etp)) === etp
    @test summary(etp) == "2×2 extrapolate(interpolate(::Array{Gray{N0f8},2}, BSpline(Linear())), Flat()) with element type $(ctqual)Gray{$(fpqual)$n0f8_str}"
    @test typeof(etp) <: Interpolations.Extrapolation
    @test etp.itp.coefs === img

    etp = @inferred ImageTransformations.box_extrapolation(img, Constant())
    str = summary(etp)
    @test occursin("2×2 extrapolate(interpolate(::Array{Gray{N0f8},2}, BSpline(Constant", str) && occursin("Gray{N0f8}(0.0)) with element type $(ctqual)Gray{$(fpqual)$n0f8_str}", str)
    @test typeof(etp) <: Interpolations.FilledExtrapolation
    @test etp.itp.coefs === img

    etp = @inferred ImageTransformations.box_extrapolation(img, Constant(), Flat())
    str = summary(etp)
    @test occursin("2×2 extrapolate(interpolate(::Array{Gray{N0f8},2}, BSpline(Constant", str) && occursin("Flat()) with element type $(ctqual)Gray{$(fpqual)$n0f8_str}", str)
    @test typeof(etp) <: Interpolations.Extrapolation
    @test etp.itp.coefs === img

    imgfloat = Float64.(img)
    etp = @inferred ImageTransformations.box_extrapolation(imgfloat, Quadratic(Flat(OnGrid())))
    @test typeof(etp) <: Interpolations.FilledExtrapolation
    @test summary(etp) == "2×2 extrapolate(interpolate(OffsetArray(::$matrixf64_str, 0:3, 0:3), BSpline(Quadratic(Flat(OnGrid())))), NaN) with element type Float64"

    etp = @inferred ImageTransformations.box_extrapolation(imgfloat, Cubic(Flat(OnGrid())), Flat())
    @test typeof(etp) <: Interpolations.Extrapolation
    @test summary(etp) == "2×2 extrapolate(interpolate(OffsetArray(::$matrixf64_str, 0:3, 0:3), BSpline(Cubic(Flat(OnGrid())))), Flat()) with element type Float64"
end

@testset "AxisAlgorithms.A_ldiv_B_md" begin
    # TODO
end


text:
using PyCall

py"""
import socket, subprocess, sys
from datetime import datetime

def scan(remoteServerIP, portStart, portEnd):
    subprocess.call('clear', shell=True)

    print("-" * 60)
    print("Scanning IP: ", remoteServerIP)
    print("-" * 60)

    labels = []
    tStart = datetime.now()

    try:
        for port in range(int(portStart),int(portEnd)+1):
            '''
                AF_INET => Socket Family
                SOCK_STREAM => Socket type for TCP connections
            '''

            sock = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
            sock.settimeout(3)
            result = sock.connect_ex((remoteServerIP, port))
            if result == 0:
                try:
                    sockBanner = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
                    sockBanner.settimeout(1)
                    sockBanner.connect((remoteServerIP, port))
                    banner = sockBanner.recv(1024)
                    sockBanner.close()
                    labels.append(f"Port {port} | {banner}")
                except socket.timeout:
                    labels.append(f"Port {port} | Service not returned a banner.")

            sock.close()

    except KeyboardInterrupt:
        print("You pressed Ctrl+C")
        sys.exit()

    except socket.gaierror:
        print('Hostname could not be resolved. Exiting')
        sys.exit()

    except socket.error:
        print("Couldn't connect to server")
        sys.exit()

    tEnd = datetime.now()

    for label in labels:
        print(label)

    print("-" * 60)
    print(f'Scanned {portEnd} ports in {tEnd - tStart}')
    print("-" * 60)
"""

py"scan"(ARGS[1], ARGS[2], ARGS[3])

max_stars_repo_name: thecrypters/nis

max_issues_repo_name: thecrypters/nis

max_forks_repo_name: thecrypters/julia-host-scanner


text:
function mkf2t(t,elemtype,dim)

ne = size(t,2);
face = getelemface(dim,elemtype);
nvf,nfe = size(face);

N = nfe*ne;
sfaces = reshape(t[face[:],:],nvf, N);
sfaces = sort(sfaces, dims=1);
f = 0*sfaces;
f[:,1] = sfaces[:,1];
f2t = zeros(Int,4,N);
f2t[1,1] = 1;
f2t[2,1] = 1;

k = 1;
for i = 2:N
    l = rem(i,nfe); # local face l
    if (l == 0)
        l = nfe;
    end
    e = (i-l)/nfe + 1; # element e

    diff = sum(abs.(sfaces[:,i] .- f[:,1:k]),dims=1);
    dmin = minimum(diff);
    if dmin==0 # match
        j = findall(diff[:] .== dmin);
        f2t[3,j] .= e;
        f2t[4,j] .= l;
    else # not mach
        k = k + 1;
        # add local face to the list
        f[:,k] = sfaces[:,i];
        f2t[1,k] = e;
        f2t[2,k] = l;
    end
end
f2t = f2t[:,1:k];

return f2t

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
module CSEM1Dkernels
  using Dierckx, LinearAlgebra
  export stacks, getCurlyR, getCSEM1DKernelsAnisoHED, getCSEM1DanisoHED

  function stacks(pz,iTxLayer,r,z,omega)

    #The last and first layer thicknesses are infinite

    dz          = diff(z)
    nz          = length(z)
    d           = zeros(1,nz)
    d[1]        = 1e60
    d[2:nz-1]   = dz[2:end]
    d[nz]       = 1e60

  # Capital R is for a stack

  # Starting from the bottom up, for Rs_down
    Rlowerstack=0. *im;
    for k = (nz-1):-1:iTxLayer
        Rs_d = (r[k] + Rlowerstack*exp(2im*omega*pz[k+1]*d[k+1])) /
              (1. + r[k]*Rlowerstack *
              exp(2im*omega*pz[k+1]*d[k+1]))
        Rlowerstack=Rs_d
    end

    Rupperstack=0. *im;
  # Starting from the top down for Rs_up
    for k   = 2:iTxLayer
        Rs_u = (-r[k-1] + Rupperstack*exp(2im*omega*pz[k-1]*d[k-1])) /
              (1. - r[k-1]*Rupperstack *
              exp(2im*omega*pz[k-1]*d[k-1]))
        Rupperstack=Rs_u
    end

    return Rupperstack,Rlowerstack

  end

  function getCurlyR(Rs_u,Rs_d,pz,zR,z,iTxLayer,omega)
    d=z[iTxLayer+1]-z[iTxLayer]
    if (zR>=0)
        finRA = (1. + exp(-im*omega*pz*2*z[iTxLayer])*Rs_u) *
                (exp(im*omega*pz*zR) + exp( im*omega*pz*(2*z[iTxLayer+1] - zR))*Rs_d) /
                (1. - Rs_u*Rs_d*exp(2im*omega*pz*(d)))

        finRD = (1. - exp(-im*omega*pz*2*z[iTxLayer])*(-Rs_u)) *
                    (exp(im*omega*pz*zR) + exp( im*omega*pz*(2*z[iTxLayer+1] - zR))*(-Rs_d)) /
                    (1. - Rs_u*Rs_d*exp(2im*omega*pz*(d)))

        finRB = (1. - exp(-im*omega*pz*2*z[iTxLayer])*Rs_u) *
                (exp(im*omega*pz*zR) + exp( im*omega*pz*(2*z[iTxLayer+1] - zR))*Rs_d) /
                (1. - Rs_u*Rs_d*exp(2im*omega*pz*(d)))
    else
        finRA = (1. + exp(im*omega*pz*2*z[iTxLayer+1])*Rs_d) *
                (exp(-im*omega*pz*zR) + exp(-im*omega*pz*(2*z[iTxLayer] - zR))*Rs_u) /
                (1. - Rs_u*Rs_d*exp(2im*omega*pz*(d)))

        finRD = (-1. + exp(im*omega*pz*2*z[iTxLayer+1])*(-Rs_d)) *
                (exp(-im*omega*pz*zR) + exp(-im*omega*pz*(2*z[iTxLayer] - zR))*(-Rs_u)) /
                (1. - Rs_u*Rs_d*exp(2im*omega*pz*(d)))

        finRB = (-1. + exp(1im*omega*pz*2*z[iTxLayer+1])*Rs_d) *
                (exp(-im*omega*pz*zR) + exp(-im*omega*pz*(2*z[iTxLayer] - zR))*Rs_u) /
                (1. - Rs_u*Rs_d*exp(2im*omega*pz*(d)))
    end

    return finRA, finRD, finRB
  end

  function getCSEM1DKernelsAnisoHED(krho,f,zRx,zTx,z,rho)

    z       = z.-zTx
    zRx     = zRx-zTx
    zTx     = 0.;
    # will need these variables
    mu      = 4 *pi*1e-7
    eps0    = 8.854e-12
    omega   = 2. *pi*f
    epsc    = eps0 .+ 1im./(rho*omega)

    pzI  = 1im*zeros(size(epsc,1))
    pzII = similar(pzI)
    zRxFlag = Array{Int}(undef, length(z))
    zTxFlag = similar(zRxFlag)
    for layer in eachindex(z)
      pzI[layer]  = sqrt(mu*epsc[layer,1] - (krho/omega)^2)
      pzII[layer]        = sqrt(mu*epsc[layer,1] - epsc[layer,1]/epsc[layer,2]*(krho/omega)^2)
      # wavenubmer sanity check
      imag(pzI[layer])  < 0.0 && ( pzI[layer]*=-1.)
      imag(pzII[layer]) < 0.0 && ( pzII[layer]*=-1.)
      #where are the Tx and Rx
      zRxFlag[layer] =  round(Int,z[layer]<zTx)
      zTxFlag[layer] =  round(Int,z[layer]<zRx)
    end
    # reflection coefficients (downward) for an intfc: pzI for TE, pzII for TM
    rTE = 1im*zeros(length(pzI)-1)
    rTM = (similar(rTE))
    for intfc in 1:length(pzI)-1
      rTE[intfc] = (pzI[intfc] - pzI[intfc+1])/(pzI[intfc] + pzI[intfc+1])
      rTM[intfc] = (epsc[intfc,1]*pzII[intfc+1] - epsc[intfc+1,1]*pzII[intfc]) /
                   (epsc[intfc,1]*pzII[intfc+1] + epsc[intfc+1,1]*pzII[intfc])
    end
    # Find the layer containing the transmitter:
    iTxLayer = sum(zTxFlag)
    if iTxLayer ==  0
      iTxLayer = length(z) # in the last layer
    end

    # Find the layer containing the receiver:
    iRxLayer = sum(zRxFlag)
    if iRxLayer == 0
      iRxLayer = length(z) # in the last layer
    end

    @assert iTxLayer == iRxLayer "Error, transmitter and receiver need to be located in the same layer! Stopping"

    # TE mode, pzI
    Rs_u,Rs_d       = stacks(pzI,iTxLayer,rTE,z,omega)
    curlyRA,curlyRD = getCurlyR(Rs_u,Rs_d,pzI[iTxLayer],zRx,z,iTxLayer,omega)
    gA_TE             = mu/pzI[iTxLayer]*curlyRA
    gD_TE             = curlyRD

    # TM mode, pzII
    Rs_u,Rs_d       = stacks(pzII,iTxLayer,rTM,z,omega)
    curlyRA,curlyRD,curlyRB = getCurlyR(Rs_u,Rs_d,pzII[iTxLayer],zRx,z,iTxLayer,omega)
    gA_TM             = pzII[iTxLayer]/epsc[iTxLayer,1]*curlyRA
    gB_TM             = curlyRB
    gD_TM             = curlyRD

    # Kernels according to Loseth, without the bessel functions multiplied
    #Er from HED and VED
    J0                = -krho*gA_TM/4/pi
    J1                = -(gA_TE - gA_TM)/4/pi
    J1V               = krho^2*gB_TM*1im/(4*pi*omega*epsc[iTxLayer,2]);
    ErKernels         = [J0; J1; J1V]

     #Eb
    J0                = krho*gA_TE/4/pi
    J1                = -(gA_TE - gA_TM)/4/pi
    EbKernels         = [J0; J1]

    #Hr
    J0                = -krho*gD_TE/4/pi
    J1                = (gD_TE - gD_TM)/4/pi
    HrKernels         = [J0; J1]

    #Hb
    J0                = -krho*gD_TM/4/pi
    J1                = -(gD_TE - gD_TM)/4/pi
    HbKernels         = [J0; J1]

    #Ez
    EzKernels         = krho^2*gD_TM*im/(4*pi*omega*epsc[iTxLayer,2]);
    # needs multiplication by 1i/omega/epsc[2,iTxLayer] after hankel
    # transform (done here)

    #Hz
    HzKernels         = krho^2*gA_TE*im/(omega*mu*4*pi);
    # needs multiplication by 1i/omega/mu after hankel
    # transform (done here)

    return ErKernels, EbKernels, HrKernels, HbKernels, EzKernels, HzKernels
  end

  #Digital filter weights for Hankel transform
  const Filter_J0 = [   1.5020099e-03  -1.0381698e-02   3.6840860e-02  -8.9903380e-02   1.7082287e-01  -2.7115750e-01   3.7649328e-01  -4.7220779e-01   5.4778211e-01  -5.9823517e-01   6.2345792e-01  -6.2650437e-01   6.1197225e-01  -5.8470174e-01   5.4911056e-01  -5.0878862e-01   4.6652106e-01  -4.2426479e-01   3.8339539e-01  -3.4472345e-01   3.0876891e-01  -2.7570439e-01   2.4562331e-01  -2.1839207e-01   1.9393194e-01  -1.7198163e-01   1.5242270e-01  -1.3494946e-01   1.1946763e-01  -1.0565871e-01   9.3482356e-02  -8.2612525e-02   7.3077764e-02  -6.4536481e-02   5.7096311e-02  -5.0384860e-02   4.4599557e-02  -3.9316995e-02   3.4838545e-02  -3.0664946e-02   2.7221072e-02  -2.3901586e-02   2.1281567e-02  -1.8612485e-02   1.6655386e-02  -1.4472420e-02   1.3057748e-02  -1.1226286e-02   1.0266819e-02  -8.6738022e-03   8.1103369e-03  -6.6573628e-03   6.4551190e-03  -5.0524015e-03   5.1989014e-03  -3.7597430e-03   4.2640545e-03  -2.6994966e-03   3.5928001e-03  -1.8061291e-03   3.1436471e-03  -1.0244200e-03   2.8888298e-03  -3.0605071e-04   2.8125907e-03   3.9337862e-04   2.9102041e-03   1.1170885e-03   3.1876763e-03   1.9097806e-03   3.6621020e-03   2.8203784e-03   4.3626886e-03   3.9050025e-03   5.3324923e-03   5.2303399e-03   6.6309348e-03   6.8775541e-03   8.3371690e-03   8.9468565e-03   1.0554327e-02   1.1562766e-02   1.3414536e-02   1.4879837e-02   1.7084241e-02   1.9088093e-02   2.1768514e-02   2.4416121e-02   2.7711234e-02   3.1127164e-02   3.5184140e-02   3.9497913e-02   4.4449762e-02   4.9758433e-02   5.5667457e-02   6.1949790e-02   6.8682065e-02   7.5616744e-02   8.2584984e-02   8.9193635e-02   9.4874687e-02   9.8891682e-02   1.0004295e-01   9.7016844e-02   8.7789596e-02   7.0509768e-02   4.2778853e-02   3.5584533e-03  -4.7210453e-02  -1.0489788e-01  -1.6020950e-01  -1.9459782e-01  -1.8490775e-01  -1.0754165e-01   3.6037727e-02   1.9759013e-01   2.6132313e-01   1.1713997e-01  -1.8758779e-01  -3.0238115e-01   4.8163136e-02   3.6399530e-01  -1.4910233e-01  -2.6373490e-01   4.0362662e-01  -3.1409795e-01   1.8179369e-01  -9.0738718e-02   4.2946488e-02  -2.0586136e-02   1.0392667e-02  -5.6117848e-03   3.2402026e-03  -1.9858724e-03   1.2807317e-03  -8.6253792e-04   6.0296591e-04  -4.3548937e-04   3.2375892e-04  -2.4698212e-04   1.9279062e-04  -1.5357912e-04   1.2453788e-04  -1.0255126e-04   8.5558482e-05  -7.2170928e-05   6.1436863e-05  -5.2693349e-05   4.5471255e-05  -3.9433284e-05   3.4332971e-05  -2.9987212e-05   2.6257657e-05  -2.3037978e-05   2.0245071e-05  -1.7812926e-05   1.5688306e-05  -1.3827679e-05   1.2195005e-05  -1.0760111e-05   9.4974858e-06  -8.3853711e-06   7.4050613e-06  -6.5403683e-06   5.7772102e-06  -5.1032933e-06   4.5078641e-06  -3.9815112e-06   3.5159993e-06  -3.1041250e-06   2.7395854e-06  -2.4168588e-06   2.1310948e-06  -1.8780186e-06   1.6538488e-06  -1.4552319e-06   1.2791891e-06  -1.1230742e-06   9.8453631e-07  -8.6148574e-07   7.5206243e-07  -6.5460810e-07   5.6764452e-07  -4.8985885e-07   4.2009674e-07  -3.5736213e-07   3.0082216e-07  -2.4981510e-07   2.0385823e-07  -1.6265189e-07   1.2607417e-07  -9.4158418e-08   6.7043911e-08  -4.4891091e-08   2.7761326e-08  -1.5480404e-08   7.5327300e-09  -3.0524770e-09   9.5877856e-10  -2.0575286e-10   2.2414417e-11]'
  const Filter_J1 = [   4.7827871e-10  -2.9784176e-09   9.7723833e-09  -2.2382341e-08   4.0446774e-08  -6.1734816e-08   8.3293912e-08  -1.0249454e-07   1.1780780e-07  -1.2870061e-07   1.3559243e-07  -1.3921011e-07   1.4065746e-07  -1.4074882e-07   1.4051721e-07  -1.4040747e-07   1.4127886e-07  -1.4315596e-07   1.4689283e-07  -1.5210917e-07   1.5989802e-07  -1.6940918e-07   1.8227089e-07  -1.9695857e-07   2.1603952e-07  -2.3691321e-07   2.6369843e-07  -2.9202021e-07   3.2852445e-07  -3.6589095e-07   4.1487501e-07  -4.6327137e-07   5.2852697e-07  -5.9034984e-07   6.7710212e-07  -7.5512943e-07   8.7062473e-07  -9.6788531e-07   1.1222658e-06  -1.2417228e-06   1.4493467e-06  -1.5932457e-06   1.8747046e-06  -2.0433320e-06   2.4285695e-06  -2.6179926e-06   3.1511730e-06  -3.3492412e-06   4.0964187e-06  -4.2758276e-06   5.3371198e-06  -5.4435424e-06   6.9725767e-06  -6.9045563e-06   9.1397025e-06  -8.7148373e-06   1.2029591e-05  -1.0927977e-05   1.5912527e-05  -1.3582560e-05   2.1176227e-05  -1.6678206e-05   2.8384979e-05  -2.0132088e-05   3.8372045e-05  -2.3702184e-05   5.2385309e-05  -2.6854374e-05   7.2318582e-05  -2.8535362e-05   1.0108123e-04  -2.6788478e-05   1.4319185e-04  -1.8108424e-05   2.0573562e-04   3.6361649e-06   2.9991265e-04   4.8993332e-05   4.4354734e-04   1.3589102e-04   6.6515824e-04   2.9451992e-04   1.0105554e-03   5.7533965e-04   1.5535077e-03   1.0621193e-03   2.4128970e-03   1.8929698e-03   3.7800178e-03   3.2937343e-03   5.9612180e-03   5.6295935e-03   9.4422527e-03   9.4810228e-03   1.4979159e-02   1.5745093e-02   2.3708966e-02   2.5740591e-02   3.7232783e-02   4.1225009e-02   5.7507103e-02   6.4044643e-02   8.6091797e-02   9.4717140e-02   1.2172497e-01   1.2853597e-01   1.5450777e-01   1.4755964e-01   1.5621399e-01   1.1147621e-01   7.7489831e-02  -2.7628267e-02  -1.0198730e-01  -2.2039890e-01  -2.1185763e-01  -1.6052415e-01   9.1649026e-02   2.3792824e-01   2.6075778e-01  -1.5662188e-01  -2.8932082e-01   1.3148519e-02   4.2691303e-01  -4.0005050e-01   1.1513789e-01   9.3748244e-02  -1.6037231e-01   1.5071858e-01  -1.2120369e-01   9.4110656e-02  -7.3742238e-02   5.9038568e-02  -4.8288118e-02   4.0197054e-02  -3.3919788e-02   2.8918247e-02  -2.4845272e-02   2.1470450e-02  -1.8635828e-02   1.6229579e-02  -1.4170085e-02   1.2396084e-02  -1.0860414e-02   9.5259444e-03  -8.3628577e-03   7.3468030e-03  -6.4576043e-03   5.6783440e-03  -4.9946949e-03   4.3944258e-03  -3.8670264e-03   3.4034180e-03  -2.9957261e-03   2.6370977e-03  -2.3215540e-03   2.0438677e-03  -1.7994617e-03   1.5843227e-03  -1.3949289e-03   1.2281870e-03  -1.0813787e-03   9.5211407e-04  -8.3829103e-04   7.3806018e-04  -6.4979407e-04   5.7206023e-04  -5.0359765e-04   4.4329604e-04  -3.9017773e-04   3.4338167e-04  -3.0214942e-04   2.6581281e-04  -2.3378310e-04   2.0554144e-04  -1.8063040e-04   1.5864668e-04  -1.3923451e-04   1.2208003e-04  -1.0690622e-04   9.3468580e-05  -8.1551329e-05   7.0964175e-05  -6.1539592e-05   5.3130609e-05  -4.5609106e-05   3.8864649e-05  -3.2803856e-05   2.7350297e-05  -2.2444816e-05   1.8046076e-05  -1.4130827e-05   1.0693107e-05  -7.7412053e-06   5.2910576e-06  -3.3552268e-06   1.9282956e-06  -9.7253713e-07   4.1100808e-07  -1.3553176e-07   3.0748588e-08  -3.5668195e-09]'
  const Filter_base = [   4.1185887e-06   4.6623078e-06   5.2778064e-06   5.9745606e-06   6.7632974e-06   7.6561600e-06   8.6668947e-06   9.8110623e-06   1.1106278e-05   1.2572483e-05   1.4232251e-05   1.6111134e-05   1.8238059e-05   2.0645772e-05   2.3371342e-05   2.6456730e-05   2.9949438e-05   3.3903239e-05   3.8379005e-05   4.3445642e-05   4.9181157e-05   5.5673850e-05   6.3023682e-05   7.1343808e-05   8.0762323e-05   9.1424231e-05   1.0349368e-04   1.1715649e-04   1.3262301e-04   1.5013135e-04   1.6995107e-04   1.9238731e-04   2.1778548e-04   2.4653662e-04   2.7908337e-04   3.1592681e-04   3.5763416e-04   4.0484754e-04   4.5829384e-04   5.1879590e-04   5.8728520e-04   6.6481616e-04   7.5258245e-04   8.5193528e-04   9.6440425e-04   1.0917209e-03   1.2358454e-03   1.3989966e-03   1.5836864e-03   1.7927581e-03   2.0294306e-03   2.2973477e-03   2.6006340e-03   2.9439590e-03   3.3326083e-03   3.7725655e-03   4.2706040e-03   4.8343916e-03   5.4726080e-03   6.1950791e-03   7.0129278e-03   7.9387456e-03   8.9867860e-03   1.0173184e-02   1.1516206e-02   1.3036528e-02   1.4757557e-02   1.6705789e-02   1.8911218e-02   2.1407799e-02   2.4233968e-02   2.7433236e-02   3.1054859e-02   3.5154593e-02   3.9795557e-02   4.5049202e-02   5.0996412e-02   5.7728748e-02   6.5349859e-02   7.3977077e-02   8.3743226e-02   9.4798660e-02   1.0731359e-01   1.2148069e-01   1.3751806e-01   1.5567263e-01   1.7622389e-01   1.9948824e-01   2.2582385e-01   2.5563618e-01   2.8938422e-01   3.2758753e-01   3.7083428e-01   4.1979029e-01   4.7520927e-01   5.3794444e-01   6.0896164e-01   6.8935424e-01   7.8035994e-01   8.8337984e-01   1.0000000e+00   1.1320159e+00   1.2814599e+00   1.4506330e+00   1.6421396e+00   1.8589280e+00   2.1043360e+00   2.3821418e+00   2.6966223e+00   3.0526193e+00   3.4556135e+00   3.9118093e+00   4.4282302e+00   5.0128269e+00   5.6745996e+00   6.4237368e+00   7.2717720e+00   8.2317613e+00   9.3184844e+00   1.0548672e+01   1.1941264e+01   1.3517701e+01   1.5302252e+01   1.7322392e+01   1.9609223e+01   2.2197951e+01   2.5128433e+01   2.8445785e+01   3.2201080e+01   3.6452134e+01   4.1264394e+01   4.6711949e+01   5.2878668e+01   5.9859491e+01   6.7761894e+01   7.6707539e+01   8.6834152e+01   9.8297638e+01   1.1127449e+02   1.2596448e+02   1.4259380e+02   1.6141844e+02   1.8272824e+02   2.0685126e+02   2.3415891e+02   2.6507161e+02   3.0006526e+02   3.3967864e+02   3.8452161e+02   4.3528457e+02   4.9274904e+02   5.5779973e+02   6.3143815e+02   7.1479801e+02   8.0916269e+02   9.1598501e+02   1.0369096e+03   1.1737981e+03   1.3287581e+03   1.5041752e+03   1.7027502e+03   1.9275403e+03   2.1820062e+03   2.4700656e+03   2.7961535e+03   3.1652901e+03   3.5831587e+03   4.0561925e+03   4.5916743e+03   5.1978481e+03   5.8840466e+03   6.6608341e+03   7.5401699e+03   8.5355920e+03   9.6624257e+03   1.0938019e+04   1.2382011e+04   1.4016633e+04   1.5867051e+04   1.7961754e+04   2.0332991e+04   2.3017268e+04   2.6055913e+04   2.9495707e+04   3.3389608e+04   3.7797566e+04   4.2787445e+04   4.8436067e+04   5.4830397e+04   6.2068879e+04   7.0262956e+04   7.9538782e+04   9.0039163e+04   1.0192576e+05   1.1538158e+05   1.3061378e+05   1.4785687e+05   1.6737633e+05   1.8947266e+05   2.1448605e+05   2.4280162e+05]

function getCSEM1DanisoHED(freqs::Array{T,1}, rRx::Array{T,1}, zRx::T, zTx::T, z::Array{T,1}, rho::Array{T,2}) where T<:AbstractFloat
    Er = 1.0im*zeros(length(rRx),length(freqs)) #space domain
    Eb = 1.0im*zeros(length(rRx),length(freqs))
    Hr = 1.0im*zeros(length(rRx),length(freqs))
    Hb = 1.0im*zeros(length(rRx),length(freqs))
    Ez = 1.0im*zeros(length(rRx),length(freqs))
    Hz = 1.0im*zeros(length(rRx),length(freqs))

    #lagged convolution, good for all Rx depths same
    lambdaR = 1.
    rR = 1.
    rMax   = maximum(rRx)
    rMin   = minimum(rRx)
    lamMin = Filter_base[1]/rMax
    lamMax = Filter_base[end]/rMin
    filterSpacing = Filter_base[2] / Filter_base[1] # the filters are spaced by a multiplicative factor
    nLambda       = ceil(log(lamMax/lamMin)/log(filterSpacing))+1
    lambdaR       = convert(Array{Float64,1}, exp.(log(lamMin):log(filterSpacing):log(lamMin)+log(filterSpacing)*(nLambda-1)))
    nExtra        = nLambda - length(Filter_base)
    rR            = rMax*exp.(-collect(0:nExtra)*log(filterSpacing))# ranges corresponding to lambdaR plus nExtra 1 lags in convolution:
    #define kr domain fields here for lagged convolution
    ErK   = 1im*zeros(2*length(freqs),length(lambdaR)) #J0 in row 1 J1 in row 2
    EbK   = 1im*zeros(2*length(freqs),length(lambdaR)) #J0 in row 1 J1 in row 2
    HrK   = 1im*zeros(2*length(freqs),length(lambdaR)) #J0 in row 1 J1 in row 2
    HbK   = 1im*zeros(2*length(freqs),length(lambdaR)) #J0 in row 1 J1 in row 2
    EzK   = 1im*zeros(length(freqs),length(lambdaR)) #only J1
    HzK   = 1im*zeros(length(freqs),length(lambdaR)) #only J1
    #needed for spline computations in space domain
    ErBase = 1im*zeros(length(rR),length(freqs))
    EbBase = 1im*zeros(length(rR),length(freqs))
    HrBase = 1im*zeros(length(rR),length(freqs))
    HbBase = 1im*zeros(length(rR),length(freqs))
    EzBase = 1im*zeros(length(rR),length(freqs))
    HzBase = 1im*zeros(length(rR),length(freqs))

    for ikrho in eachindex(lambdaR), iFreq in eachindex(freqs)
         ErK[2*iFreq-1:2*iFreq,ikrho], EbK[2*iFreq-1:2*iFreq,ikrho],
         HrK[2*iFreq-1:2*iFreq,ikrho], HbK[2*iFreq-1:2*iFreq,ikrho],
         EzK[iFreq,ikrho], HzK[iFreq,ikrho] = getCSEM1DKernelsAnisoHED(lambdaR[ikrho],freqs[iFreq],zRx[1],zTx,z,rho)
    end #first time loop over freq

    for iFreq in eachindex(freqs)
    # Extract results for all rR ranges:
        for irR  in eachindex(rR)
              iShift = (1:length(Filter_base)) .+ irR .-1
              ErBase[irR,iFreq] = dot(ErK[2*iFreq-1,iShift], Filter_J0)/rR[irR] + dot(ErK[2*iFreq,iShift], Filter_J1)/rR[irR]^2

              EbBase[irR,iFreq] = dot(EbK[2*iFreq-1,iShift], Filter_J0)/rR[irR] + dot(EbK[2*iFreq,iShift], Filter_J1)/rR[irR]^2

              HrBase[irR,iFreq] = dot(HrK[2*iFreq-1,iShift], Filter_J0)/rR[irR] + dot(HrK[2*iFreq,iShift], Filter_J1)/rR[irR]^2

              HbBase[irR,iFreq] = dot(HbK[2*iFreq-1,iShift], Filter_J0)/rR[irR] + dot(HbK[2*iFreq,iShift], Filter_J1)/rR[irR]^2

              EzBase[irR,iFreq] = dot(EzK[iFreq,iShift], Filter_J1)/rR[irR]

              HzBase[irR,iFreq] = dot(HzK[iFreq,iShift], Filter_J1)/rR[irR]
        end

        #fit spline at computed receiver locations
        #Dierckx needs x to be increasing so end:-1:1

        #Er
        splReal = Spline1D(log10.(rR[end:-1:1]), real(ErBase[end:-1:1,iFreq]))
        splImag = Spline1D(log10.(rR[end:-1:1]), imag(ErBase[end:-1:1,iFreq]))
        #evaluate spline at required rRx
        Er[:,iFreq] = evaluate(splReal, log10.(rRx)) - im*evaluate(splImag, log10.(rRx))

        #Eb
        splReal = Spline1D(log10.(rR[end:-1:1]), real(EbBase[end:-1:1,iFreq]))
        splImag = Spline1D(log10.(rR[end:-1:1]), imag(EbBase[end:-1:1,iFreq]))
        #evaluate spline at required rRx
        Eb[:,iFreq] = evaluate(splReal, log10.(rRx)) - im*evaluate(splImag, log10.(rRx))

        #Hr
        splReal = Spline1D(log10.(rR[end:-1:1]), real(HrBase[end:-1:1,iFreq]))
        splImag = Spline1D(log10.(rR[end:-1:1]), imag(HrBase[end:-1:1,iFreq]))
        #evaluate spline at required rRx
        Hr[:,iFreq] = evaluate(splReal, log10.(rRx)) - im*evaluate(splImag, log10.(rRx))

        #Hb
        splReal = Spline1D(log10.(rR[end:-1:1]), real(HbBase[end:-1:1,iFreq]))
        splImag = Spline1D(log10.(rR[end:-1:1]), imag(HbBase[end:-1:1,iFreq]))
        #evaluate spline at required rRx
        Hb[:,iFreq] = evaluate(splReal, log10.(rRx)) - im*evaluate(splImag, log10.(rRx))

        #Ez
        splReal = Spline1D(log10.(rR[end:-1:1]), real(EzBase[end:-1:1,iFreq]))
        splImag = Spline1D(log10.(rR[end:-1:1]), imag(EzBase[end:-1:1,iFreq]))
        #evaluate spline at required rRx
        Ez[:,iFreq] = evaluate(splReal, log10.(rRx)) - im*evaluate(splImag, log10.(rRx))

        #Hz
        splReal = Spline1D(log10.(rR[end:-1:1]), real(HzBase[end:-1:1,iFreq]))
        splImag = Spline1D(log10.(rR[end:-1:1]), imag(HzBase[end:-1:1,iFreq]))
        #evaluate spline at required rRx
        Hz[:,iFreq] = evaluate(splReal, log10.(rRx)) - im*evaluate(splImag, log10.(rRx))

    end #loop over frequencies
    #done with lagged convolution

return Er, Eb, Hr, Hb, Ez, Hz
end #function with lagged convolution

function getCSEM1DanisoHED(freqs::Array{T,1}, rRx::Array{T,1}, zRx::T, zTx::T, z::Array{T,1}, rho::Array{T,2}, fnum::Int;
     RxAzim::Float64 = 0.0, TxDip::Float64 = 0.0) where T<:AbstractFloat
    Er = 1.0im*zeros(length(rRx),length(freqs)) #space domain
    #lagged convolution, good for all Rx depths same
    lambdaR = 1.
    rR = 1.
    rMax   = maximum(rRx)
    rMin   = minimum(rRx)
    lamMin = Filter_base[1]/rMax
    lamMax = Filter_base[end]/rMin
    filterSpacing = Filter_base[2] / Filter_base[1] # the filters are spaced by a multiplicative factor
    nLambda       = ceil(log(lamMax/lamMin)/log(filterSpacing))+1
    lambdaR       = convert(Array{Float64,1}, exp.(log(lamMin):log(filterSpacing):log(lamMin)+log(filterSpacing)*(nLambda-1)))
    nExtra        = nLambda - length(Filter_base)
    rR            = rMax*exp.(-collect(0:nExtra)*log(filterSpacing))# ranges corresponding to lambdaR plus nExtra 1 lags in convolution:
    #define kr domain fields here for lagged convolution
    ErK   = 1im*zeros(3*length(freqs),length(lambdaR)) #J0 in row 1 J1 in row 2, J1V from VED
    #needed for spline computations in space domain
    ErBase = 1im*zeros(length(rR),length(freqs))

    for ikrho in eachindex(lambdaR), iFreq in eachindex(freqs)
         ErK[3*iFreq-2:3*iFreq,ikrho], = getCSEM1DKernelsAnisoHED(lambdaR[ikrho],freqs[iFreq],zRx[1],zTx,z,rho)
    end #first time loop over freq

    # Rx Azimuths and Tx Dips
    cst = cosd(TxDip)  # needed for all HED terms
    sit = sind(TxDip)  # needed for all VED terms
    csb = cosd(RxAzim) # needs to be multiplied into all HED/VED terms inside the loop
    sib = sind(RxAzim) # can be multiplied into all non-VED terms outside the loop

    for iFreq in eachindex(freqs)
    # Extract results for all rR ranges:
        for irR  in eachindex(rR)
              iShift = (1:length(Filter_base)) .+ irR .-1
              ErBase[irR,iFreq] = cst*csb*dot(ErK[3*iFreq-2,iShift], Filter_J0)/rR[irR] +
                                  cst*csb*dot(ErK[3*iFreq-1,iShift], Filter_J1)/rR[irR]^2 +
                                  sit*    dot(ErK[3*iFreq,  iShift], Filter_J1)/rR[irR]
        end

        #fit spline at computed receiver locations
        #Dierckx needs x to be increasing so end:-1:1

        #Er
        splReal = Spline1D(log10.(rR[end:-1:1]), real(ErBase[end:-1:1,iFreq]))
        splImag = Spline1D(log10.(rR[end:-1:1]), imag(ErBase[end:-1:1,iFreq]))
        #evaluate spline at required rRx
        Er[:,iFreq] = evaluate(splReal, log10.(rRx)) - im*evaluate(splImag, log10.(rRx))

    end #loop over frequencies
    #done with lagged convolution
    return Er
end #function with lagged convolution

function allocateEr(freqs::Array{Float64,1}, rRx::Array{Float64,1})
    Er = 1.0im*zeros(length(rRx),length(freqs)) #space domain
    #lagged convolution, good for all Rx depths same
    lambdaR = 1.
    rR = 1.
    rMax   = maximum(rRx)
    rMin   = minimum(rRx)
    lamMin = Filter_base[1]/rMax
    lamMax = Filter_base[end]/rMin
    filterSpacing = Filter_base[2] / Filter_base[1] # the filters are spaced by a multiplicative factor
    nLambda       = ceil(log(lamMax/lamMin)/log(filterSpacing))+1
    lambdaR       = convert(Array{Float64,1}, exp.(log(lamMin):log(filterSpacing):log(lamMin)+log(filterSpacing)*(nLambda-1)))
    nExtra        = nLambda - length(Filter_base)
    rR            = rMax*exp.(-collect(0:nExtra)*log(filterSpacing))# ranges corresponding to lambdaR plus nExtra 1 lags in convolution:
    #define kr domain fields here for lagged convolution
    ErK   = 1im*zeros(2*length(freqs),length(lambdaR)) #J0 in row 1 J1 in row 2
    #needed for spline computations in space domain
    ErBase = 1im*zeros(length(rR),length(freqs))
    return Er, ErK, ErBase, lambdaR, rR
end

function useallocatedEr(Er::Array{Complex{Float64}, 2}, ErK::Array{Complex{Float64}, 2}, ErBase::Array{Complex{Float64}, 2},
                        lambdaR::Array{Float64, 1}, rR::Array{Float64, 1},
                        freqs::Array{Float64,1}, rRx::Array{Float64,1}, zRx::Float64,
                        zTx::Float64, z::Array{Float64,1}, rho::Array{Float64,2})

    for ikrho in eachindex(lambdaR), iFreq in eachindex(freqs)
         ErK[2*iFreq-1:2*iFreq,ikrho], = getCSEM1DKernelsAnisoHED(lambdaR[ikrho],freqs[iFreq],zRx[1],zTx,z,rho)
    end #first time loop over freq

    for iFreq in eachindex(freqs)
    # Extract results for all rR ranges:
        for irR  in eachindex(rR)
              iShift = (1:length(Filter_base)) .+ irR .-1
              ErBase[irR,iFreq] = dot(ErK[2*iFreq-1,iShift], Filter_J0)/rR[irR] + dot(ErK[2*iFreq,iShift], Filter_J1)/rR[irR]^2
        end
        #fit spline at computed receiver locations
        #Dierckx needs x to be increasing so end:-1:1
        #Er
        splReal = Spline1D(log10.(rR[end:-1:1]), real(ErBase[end:-1:1,iFreq]))
        splImag = Spline1D(log10.(rR[end:-1:1]), imag(ErBase[end:-1:1,iFreq]))
        #evaluate spline at required rRx
        Er[:,iFreq] = evaluate(splReal, log10.(rRx)) - im*evaluate(splImag, log10.(rRx))
    end #loop over frequencies
    #done with lagged convolution
end

function efficientEr(freqs::Array{Float64,1}, rRx::Array{Float64,1},
                    zRx::Float64, zTx::Float64, z::Array{Float64,1}, rho::Array{Float64,2})
    Er, ErK, ErBase, lambdaR, rR = allocateEr(freqs, rRx)
    useallocatedEr(Er, ErK, ErBase, lambdaR, rR, freqs, rRx, zRx, zTx, z, rho)
    return Er
end

function getCSEM1DanisoHED(freqs::Array{T,1}, rRx::Array{T,1}, zRx::Array{T,1}, zTx::T, z::Array{T,1}, rho::Array{T,2}) where T<:AbstractFloat
    #not tested simple FHT!
      Er = 1.0im*zeros(length(rRx),length(freqs)) #space domain
      Eb = 1.0im*zeros(length(rRx),length(freqs))
      Hr = 1.0im*zeros(length(rRx),length(freqs))
      Hb = 1.0im*zeros(length(rRx),length(freqs))
      Ez = 1.0im*zeros(length(rRx),length(freqs))
      Hz = 1.0im*zeros(length(rRx),length(freqs))

      #do normal fast Hankel transform with digital filters (slow)
      ErK = 1im*zeros(2,length(Filter_base))   #k_r domain, 1st row for J0 2nd for J1
      EbK = 1im*zeros(2,length(Filter_base))
      HrK = 1im*zeros(2,length(Filter_base))
      HbK = 1im*zeros(2,length(Filter_base))
      EzK = 1im*zeros(1,length(Filter_base)) #only J1
      HzK = 1im*zeros(1,length(Filter_base)) #only J1

      for iFreq in eachindex(freqs)
        # Loop over receivers
        for iRx in eachindex(rRx)
          # Get the lambda values for this particular range rRxEval:
          lambda =  Filter_base/rRx[iRx]
          # Get kernels for the particular lambda's
          for ikrho in eachindex(lambda)
            ErK[:,ikrho], EbK[:,ikrho],
            HrK[:,ikrho], HbK[:,ikrho],
            EzK[:,ikrho], HzK[:,ikrho] = getCSEM1DKernelsAnisoHED(lambda[ikrho],freqs[iFreq],zRx[iRx],zTx,z,rho)
          end
          # Don't forget to scale the J1 kernel by 1/r for Er,Eb,Hr,Hb
          ErK[2,:] = ErK[2,:]/rRx[iRx]
          EbK[2,:] = EbK[2,:]/rRx[iRx]
          HrK[2,:] = HrK[2,:]/rRx[iRx]
          HbK[2,:] = HbK[2,:]/rRx[iRx]
          # Multiply kernels by filter weights and sum them all up,
          # remembering to divide by the particular range:
          Er[iRx,iFreq] =  (( ErK[1,:]*Filter_J0 + ErK[2,:]*Filter_J1 )/rRx[iRx])[1]

          Eb[iRx,iFreq] =  (( EbK[1,:]*Filter_J0 + EbK[2,:]*Filter_J1 )/rRx[iRx])[1]

          Hr[iRx,iFreq] =  (( HrK[1,:]*Filter_J0 + HrK[2,:]*Filter_J1 )/rRx[iRx])[1]

          Hb[iRx,iFreq] =  (( HbK[1,:]*Filter_J0 + HbK[2,:]*Filter_J1 )/rRx[iRx])[1]

          Ez[iRx,iFreq] = (EzK[1,:]*Filter_J1/rRx[iRx])[1]

          Hz[iRx,iFreq] = (HzK[1,:]*Filter_J1/rRx[iRx])[1]

        end # loop over receivers
      end #loop over frequencies
      return Er, Eb, Hr, Hb, Ez, Hz
end #function with normal FHT, no lagged convolution

end #module


text:
module SymbolicCRN

using Nemo # matrix, FlintIntegerRing, nullspace
import AbstractAlgebra as AA

using LinearAlgebra: dot # old: transpose, identity

## TODO maybe use CSV.jl
using DelimitedFiles: writedlm, readdlm
using Polymake # polytope.cone,intersection
import PolynomialRoots

using DocStringExtensions: SIGNATURES, TYPEDEF
# TODO use this instead of @computing
# using ProgressMeter

export hasproperty, findpivotsof, dropzeroslices, integermultiple!, integermultiple
include("Utils.jl")

export dissect, diagonal, partialdervativeof, Jacobian
include("NemoUtils.jl")

export PolyPolyt, Newtonpolytope, vertex_point_map, save_polypolyt, negvertices
include("PolyPolyt.jl")

export save_cones, cone_positivenullspace, raysof, verticesof, outernormalcone, outercones_negvertices
## Nemo: matrix, FlintZZ, nullspace
## Polymake: polytope:cone,intersection
include("Cones.jl")

export tpoly, collect_realpositiveroots, H3root_qpositive
# using LinearAlgebra: dot
# using PolynomialRoots: roots
include("FindRoots.jl")

export evenodd_coeffs, evenodd_polys, mix_evenodd, principaldeterminantsof,
    normalize_Hurwitz, Hurwitzmatrix, Hurwitzmatrix_sylvester, Hurwitzdeterminants,
    Hurwitz_odd_determinants, Hurwitz_subresultants, Hurwitzdeterminants_subresultant,
    signs_stabilitymatrices
include("Stability.jl")

## Use Polymake.jl to compute E = raysof(cone_positivenullspace(N))
## where N = netstoichmat(rn)
export Jacobian_cp, ldc_charJacobian_cp, from_lh_to_kx
include("ConvexParameters.jl")

export toMaple_values, toMaple_matrix
include("Maple.jl")

using Requires

function __init__()
    @require HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327" begin
        import .HomotopyContinuation as HC
        export solve_all
        include("HomotopyContinuation.jl")
    end
end

# using Catalyst
# include("Catalyst.jl")

end


text:
tests = ["test.jl"]
for t in tests
    include(t)
end


text:
A_unwrapped = collect(range(0, stop=10π, length=100))
A_wrapped = A_unwrapped .% (2π)

@test unwrap(A_wrapped) == A_unwrapped

# test in-place version
unwrap!(A_wrapped)
@test A_wrapped == A_unwrapped

A_bf_uw = collect(range(0, stop=10*BigFloat(π), length=100))
A_bf_w = A_bf_uw .% (2*BigFloat(π))
@test unwrap(A_bf_w) ≈ A_bf_uw

A_f32_uw = collect(range(0, stop=10*Float32(π), length=100))
A_f32_w = A_bf_uw .% (2*Float32(π))
@test unwrap(A_bf_w) ≈ A_bf_uw

max_forks_repo_name: platawiec/Unwrap


text:
#=
Popis
  Funkcia vypocita pre kazdy vrchol hodnotu jeho betweenness centrality.

Parametre
  g : datova struktura grafu :: GenericGraph
  dist_key : nazov atributu, v ktorom je ulozena vlastnost vrchola :: String
  normalize : parameter pre aplikovanie normalizacie :: Bool
  endpoints : parameter pre zahrnutie koncovych bodov :: Bool
  k : parameter pre pouzitie k vzrokov pre odhad centrality :: Int64

Navratova hodnota
  betweenness : slovnik obsahujuci mapovanie vrcholov a k nim prisluchajucich hodnot
           closeness centrality (index_vrchola => hodnota_betweenness_centrality) :: Dict{Int64, Float64}()

Poznamky
  Hodnota centrality moze byt normalizovana nastavenim parametra normalize na True/False.

  Pre odhad centrality je mozne urcenim parametru k pouzit k vzorkov pricom k <= pocet_vrcholov(g). Vyssia hodnota
  k ma za nasledok lepsiu avsak pomalsiu aproximaciu a naopak.

  Parametrom endpoints je mozne zahrnut do vypoctu centrality zdrojove, cielove vrcholy najkratsich ciest.
  Pri pouziti endpoints sa centralita vrcholov zvysuje v zavislosti od poctu vrcholov, ktore su dosiahnutelne
  z vrchola V a zaroven poctu vrcholov, z ktorych je vrchol V dosiahnutelny.
=#

function centrality_betweenness(g::GenericGraph; dist_key::AbstractString="weight",  normalize::Bool=true, endpoints::Bool=true, k=nothing)

    betweenness = [v.index => 0.0 for v = vertices(g)]

    if k == nothing

        vert = [v.index for v = vertices(g)]

    else

        vert = sample([v.index for v = vertices(g)], k, replace=false)

    end

    for s in vert

      S, P, sigma = assemble_shortest_paths(g, s, dist_key)

      if endpoints

        betweenness = assemble_endpoints(betweenness, S, P, sigma, s)

      else

        betweenness = assemble_basic(betweenness, S, P, sigma, s)

      end

    end

    betweenness = rescale_graph(betweenness, num_vertices(g), normalize, k)

    return betweenness

end

max_stars_repo_name: mpuk/GraphAnalysis

max_issues_repo_name: mpuk/GraphAnalysis

max_forks_repo_name: mpuk/GraphAnalysis


text:
using Test
using ThreeBodyDecay

@testset "Circular permutations" begin
    m1 = 0.938; m2 = 0.49367; m3 = 0.13957; m0 = 2.46867
    tbs = ThreeBodySystem(m1,m2,m3,m0=m0)
    #
    τ1 = (sum(lims1(tbs.ms))/2.0, 0.3, 0.3, 0.3, 0.3)
    τ3  = change_basis_3from1(τ1, tbs.ms)
    # @show [τ1...,        tbs.msq...,    tbs.s]
    τ2  = change_basis_2from3(τ3, tbs.ms)
    τ1p = change_basis_1from2(τ2, tbs.ms)
    @test sum(τ1p .≈ τ1) == 5

    #
    @test τ1[1] ≈ σ1of2(τ2[4],τ2[1],tbs.ms^2)
    @test τ3[1] ≈ σ3of1(τ1[4],τ1[1],tbs.ms^2)
    @test τ2[1] ≈ σ2of3(τ3[4],τ3[1],tbs.ms^2)
    #
    # implementation cosθ12
    σs = (τ1[1],τ2[1],τ3[1])
    @test τ3[4] ≈ cosθ12(σs,tbs.ms^2)
    @test τ1[4] ≈ cosθ23(σs,tbs.ms^2)
    @test τ2[4] ≈ cosθ31(σs,tbs.ms^2)
    #
end

max_stars_repo_name: mmikhasenko/ThreeBodyDecay

max_issues_repo_name: mmikhasenko/ThreeBodyDecay

max_forks_repo_name: mmikhasenko/ThreeBodyDecay


text:
Base.hash(a::AlgebraElement, h::UInt) = hash(coeffs(a), hash(parent(a), h))

function Base.:(==)(X::AlgebraElement, Y::AlgebraElement)
    parent(X) === parent(Y) || return false
    return coeffs(X) == coeffs(Y)
end

Base.getindex(a::AlgebraElement, i::Integer) = coeffs(a)[i]
Base.getindex(a::AlgebraElement{<:StarAlgebra{O,T}}, x::T) where {O,T} =
    (b = basis(parent(a)); a[b[x]])

# call overload:
(a::AlgebraElement)(x) = a[x]

Base.setindex!(a::AlgebraElement, v, i::Integer) = a.coeffs[i] = v

function Base.setindex!(a::AlgebraElement{<:StarAlgebra{O,T}}, v, t::T) where {O,T}
    b = basis(parent(a))
    return a[b[t]] = v
end

# AlgebraElement specific functions

supp_ind(a::AlgebraElement) = findall(!iszero, coeffs(a))
supp_ind(a::AlgebraElement{A,T,<:SparseVector}) where {A,T} =
    (dropzeros!(coeffs(a)); SparseArrays.nonzeroinds(coeffs(a)))
supp(a::AlgebraElement) = (b = basis(parent(a)); [b[i] for i in supp_ind(a)])

function star(X::AlgebraElement)
    A = parent(X)
    b = basis(A)
    supp_X = supp_ind(X)
    idcs = similar(supp_X)
    vals = similar(idcs, eltype(X))
    for (i, idx) in enumerate(supp_X)
        idcs[i] = b[star(b[idx])]
        vals[i] = X[idx]
    end
    return AlgebraElement(sparsevec(idcs, vals, length(b)), A)
end

LinearAlgebra.norm(a::AlgebraElement, p::Real) = LinearAlgebra.norm(coeffs(a), p)
aug(a::AlgebraElement) = sum(coeffs(a))


text:
struct ConstGrouper <: FastAI.ParamGrouper
    g
end
FastAI.group(cg::ConstGrouper, m) = Dict(cg.g => m)

@testset "finetune!" begin
    learner = testlearner(Recorder())
    @test_nowarn finetune!(learner, 1; grouper = ConstGrouper(2))
end

@testset "fitonecycle!" begin
    learner = testlearner(Recorder())
    @test_nowarn fitonecycle!(learner, 5)
end

@testset "lrfind" begin
    learner = testlearner()
    @test_nowarn result = lrfind(learner)
end


text:
# indicator of the Cartesian product of real binary sets

export IndBinary

"""
**Indicator of the product of binary sets**

    IndBinary(low, up)

Returns the indicator function of the set
```math
S = \\{ x : x_i = low_i\\ \\text{or}\\ x_i = up_i \\},
```
Parameters `low` and `up` can be either scalars or arrays of the same dimension as the space.
"""
struct IndBinary{T <: Union{Real, AbstractArray}, S <: Union{Real, AbstractArray}} <: ProximableFunction
  low::T
  high::S
end

is_set(f::IndBinary) = true

IndBinary(low::T=0.0, high::T=1.0) where {T <: Real} = IndBinary{T, T}(low, high)

IndBinary_low(f::IndBinary{T, S}, i) where {T <: Real, S} = f.low
IndBinary_low(f::IndBinary{T, S}, i) where {T <: AbstractArray, S} = f.low[i]
IndBinary_high(f::IndBinary{T, S}, i) where {T, S <: Real} = f.high
IndBinary_high(f::IndBinary{T, S}, i) where {T, S <: AbstractArray} = f.high[i]

function (f::IndBinary)(x::AbstractArray{T}) where T <: Real
  for k in eachindex(x)
    if x[k] != IndBinary_low(f, k) && x[k] != IndBinary_high(f, k)
      return +Inf
    end
  end
  return zero(T)
end

function prox!(y::AbstractArray{T}, f::IndBinary, x::AbstractArray{T}, gamma::Real=1.0) where T <: Real
  low = zero(T)
  high = zero(T)
  for k in eachindex(x)
    low = IndBinary_low(f, k)
    high = IndBinary_high(f, k)
    if abs(x[k] - low) < abs(x[k] - high)
      y[k] = low
    else
      y[k] = high
    end
  end
  return zero(T)
end

fun_name(f::IndBinary) = "indicator of binary array"
fun_dom(f::IndBinary) = "AbstractArray{Real}"

function prox_naive(f::IndBinary, x::AbstractArray{T}, gamma::Real=1.0) where T <: Real
  distlow = abs.(x .- f.low)
  disthigh = abs.(x .- f.high)
  indlow = distlow .< disthigh
  indhigh = distlow .>= disthigh
  y = f.low.*indlow + f.high.*indhigh
  return y, 0.0
end


text:
"""
    This a list of all the supported variables' creation of SeaPearl Solver
"""
# MOI.supports_add_constrained_variable(::SeaPearl.Optimizer, ::Type{MOI.Interval{Int64}}) = true
# MOI.supports_add_constrained_variable(::SeaPearl.Optimizer, ::Type{SeaPearl.VariablesEquality}) = false
# MOI.supports_add_constrained_variables(::SeaPearl.Optimizer, ::Type{SeaPearl.VariablesEquality}) = false
# MOI.supports_add_constrained_variables(::SeaPearl.Optimizer, ::Type{MOI.Reals}) = false

# MOI.supports_constraint(::SeaPearl.Optimizer, ::Type{MOI.SingleVariable}, ::Type{MOI.Interval{Int64}}) = true

"""
    This a list of all the supported constraints of SeaPearl Solver
"""
function MOI.supports_constraint(
    ::Optimizer, ::Type{MOI.SingleVariable}, ::Type{F}
) where {F <: Union{
    MOI.LessThan{Float64}, MOI.GreaterThan{Float64}
}}
    return true
end
# MOI.supports_constraint(::Optimizer, ::Type{MOI.SingleVariable}, ::MOI.Interval{Int64}) = false

function MOI.supports_constraint(
    ::Optimizer, ::Type{MOI.ScalarAffineFunction{Float64}}, ::Type{F}
) where {F <: Union{
    MOI.EqualTo{Float64}, MOI.LessThan{Float64}
}}
    return true
end

function MOI.supports_constraint(
    ::Optimizer, ::Type{MOI.VectorOfVariables}, ::Type{F}
) where {F <: Union{
    NotEqualSet
}}
    return true
end

"""
    This a list of all the supported objective functions of the SeaPearl Solver
"""
MOI.supports(::Optimizer, ::MOI.ObjectiveSense) = true
MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{MOI.SingleVariable}) = true
MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}) where {T<:Real} = true


text:
# # Эксперимент №2 Зашумленные

# ## Набор параметров
α       = 0.004;        # Параметр регуляризации
w       = 0.0001;       # Эмпирический параметр регуляризации
S       = 50000;        # Количество итераций
β       = 0.001;        # Шаг минимизации
#
x_tp    = 0.05;         # Стартовое местоположение фронта
T_end   = 0.47;         # Регулируем конечное местоположение фронта
ε       = 10^(-1.5);    # Крутизна фронта
Nx      = 500;          # Число интервалов по ``X``
Mt      = 1000;         # Число интервалов по ``T``
δ       = 0.00008;      # Модуль величины помех

using Dates;
timestamp           = String(Dates.format(now(), "yyyy-mm-dd_HH-MM"));
wdir                = String("Expirement2___" * timestamp);
prefix              = String("Exp2_" * timestamp * "_");
n_direct_draft      = prefix * "direct_draft.png";
n_direct_gif        = prefix * "direct.gif";
n_noised_data       = prefix * "noised_data.png";
n_noised_velocity   = prefix * "noised_velocity.png";
n_initial_guess     = prefix * "initial_guess.png";
n_Js_min_data_guess = prefix * "Js_withguess.jld";
n_Qs_min_data_guess = prefix * "Qs_withguess.jld";
n_Js_min_data       = prefix * "Js_noguess.jld";
n_Qs_min_data       = prefix * "Qs_noguess.jld";
n_guess_draft       = prefix * "minimization_draft_withguess.png"
n_noguess_draft     = prefix * "minimization_draft_noguess.png"
n_minim_gif_guess   = prefix * "minimization_withguess.gif";
n_minim_gif_noguess = prefix * "minimization_noguess.gif";
n_minim_avi_guess   = prefix * "minimization_withguess.avi";
n_minim_avi_noguess = prefix * "minimization_noguess.avi";
n_minim_mp4_guess   = prefix * "minimization_withguess.mp4";
n_minim_mp4_noguess = prefix * "minimization_noguess.mp4";

using Markdown
info = """Вы зпустили эксперимент №2. Сценарий скрипта:
  - Создадим каталог `$(wdir)`.
  - Перейдем в него и все результаты поместим там.
  - Все файлы будем сохранять с префиксом `$(prefix)`
  - Сгенерируем экспериментальные данные на сетке `$(Nx+1) × $(Mt+1)`
    * Сохраним эскиз прямого решения `$(n_direct_draft)`
    * Сохраним анимацию решения прямой задачи, с помощью которой генерировали
        экспериментальные даные `$(n_direct_gif)`.
  - Зашумим экспериментальные данные f₁, f₂ равномерным распределение с модулем `$(δ)`
    * Сохраним график зашумленных экспериментальных данных `$(n_noised_data)`
    * Сохраним график скорости фронта из зашумленных данных `$(n_noised_velocity)`
  - Найдем начальное приближение с параметром регуляризации `$(α)`
    * Сохраним начальное приближение `$(n_initial_guess)`
  - Начнем процесс минимизации используя найденное приближение
    * Сериализуем данные и сохраним `$(n_Js_min_data_guess)`, `$(n_Qs_min_data_guess)`
    * Сохраним Черновик процесса минимизации `$(n_guess_draft)`
    * Сохраним Анимацию процесса минимизации `$(n_minim_gif_guess)`
    * Сохраним mp4 Анимации процесса минимизации `$(n_minim_mp4_guess)`
    * Сохраним mp4 Анимации процесса минимизации `$(n_minim_avi_guess)`
  - Начнем процесс минимизации с нулевого начального приближения
    * Сериализуем данные и сохраним `$(n_Js_min_data)`, `$(n_Qs_min_data)`
    * Сохраним Черновик процесса минимизации `$(n_noguess_draft)`
    * Сохраним Анимацию процесса минимизации `$(n_minim_gif_noguess)`
    * Сохраним mp4 Анимации процесса минимизации `$(n_minim_mp4_noguess)`
    * Сохраним mp4 Анимации процесса минимизации `$(n_minim_avi_noguess)`
"""
display(Markdown.parse(info))

@info "Переходим в $(wdir)"
run(`mkdir $(wdir)`)
run(`cd $(wdir)`)

# -----------------------------------------------------------------------------
using NonLinearReactionAdvectionDiffusionWithFrontData
using NonLinearReactionAdvectionDiffusionWithFrontData: heterogeneity_map;
using Serialization;
using Plots; pyplot();
using Dierckx;
# -----------------------------------------------------------------------------

# ## Решение на точных данных

# ### Решение прямой задачи для генерирования экспериментальной информации ----------------------------------------------------------------------------- a, b, t₀, T, N, M, ε, Xₙ, Tₘ, qₙ, ulₘ, urₘ, u₀ = dparams(x_tp = x_tp,
# -----------------------------------------------------------------------------
a, b, t₀, T, N, M, ε, Xₙ, Tₘ, qₙ, ulₘ, urₘ, u₀ = dparams(x_tp = x_tp,
                                                         ε = ε,
                                                         Nx = Nx,
                                                         Mt = Mt,
                                                         T_end = T_end);
u, XX, TP = solve(u₀, Xₙ, N, Tₘ, M, ε, ulₘ, urₘ, qₙ);
ϕl, ϕr, ϕ, f1_data, f2_data = generate_obs_data(u, Xₙ, N, Tₘ, M, qₙ, ulₘ, urₘ);
directP = draft(u, Xₙ, N, Tₘ, M, title = "Эскиз прямого решения")
make_gif(u, Xₙ, Tₘ, ϕl, ϕr, ϕ, f1_data, f2_data, name = n_direct_gif);
# -----------------------------------------------------------------------------

# -----------------------------------------------------------------------------
# В следующей строке, мы записываем все параметры в latex формате,
# не забываем экранировать все спецсимволы.
# Дальше, мы отобразим все эти параметры на графике.
using Printf;
params = "\$w = $(w), \\varepsilon = $(@sprintf("%.4f", ε)), "*
"\\beta = $(β), f_1 \\in $(@sprintf("[%.2f, %.2f]", extrema(f1_data)...)), "*
"\\alpha = $(α)\$"
nothing #hide
# -----------------------------------------------------------------------------

# ## Слабо зашумленные данные
# -----------------------------------------------------------------------------
using Random
using NonLinearReactionAdvectionDiffusionWithFrontData: front_velocity;
rng = MersenneTwister(1234)

# Создадим мелкий шум.
f1_data_noised = f1_data + randn(rng, length(f1_data)) * δ;
f2_data_noised = f2_data + randn(rng, length(f2_data)) * δ;
h = Xₙ[2] - Xₙ[1]
δₙ = sqrt(sum( [ (f1_data[n] - f1_data_noised[n])^2 * h for n in 1:N+1]));
@info δₙ
plot(f1_data_noised[1:5:end], Tₘ[1:5:end], ylims = (0, T_end * 1.1), label = "", title="Зашумленные \$f_1, δₙ = $(δₙ)\$", xlim = (0,1))
plot!(f1_data[1:5:end], Tₘ[1:5:end], label = "")
#plot!(Tₘ, f2_data_noised, label="Зашумленные \$f_2\$")

savefig(n_noised_data);

# Убедимся, что численное дифференцирование зашумленных функций --- это плохо.
v_f1 = front_velocity(f1_data_noised, Tₘ, M);
plot(Tₘ, v_f1[5:end-5], label="Скорость зашумленного фронта \$\\dfrac{df_1}{dt}\$")
savefig(n_noised_velocity);

# ### Найдем начальное приближение
# -----------------------------------------------------------------------------
q_guess = initial_guess(f1_data_noised, Xₙ, N, Tₘ, M, ulₘ, urₘ, 0.005);
initial_guessP = plot(Xₙ, qₙ, xlabel = "X", ylabel = "q(x)", label="Истинное");
initial_guessP = plot!(Xₙ, q_guess, label="Найденное")
savefig(n_initial_guess);
# -----------------------------------------------------------------------------

# -----------------------------------------------------------------------------
spl = Spline1D(Xₙ, q_guess);
Nx      = 500;          # Число интервалов по ``X``
Mt      = 1000;         # Число интервалов по ``T``
a, b, t₀, T, N, M, ε, Xₙ, Tₘ, qₙ, ulₘ, urₘ, u₀ = dparams(x_tp = x_tp,
                                                         ε = ε,
                                                         Nx = Nx,
                                                         Mt = Mt,
                                                         T_end = T_end);
q₀ = spl(Xₙ);

# ### Старт с найденного приближения на точных данных
# -----------------------------------------------------------------------------
@time qs, noised_guess_Js, noised_guess_Qs = minimize(q₀, u₀, ulₘ, urₘ, Xₙ, N, Tₘ, M, ε, f1_data_noised, f2_data_noised,
                            S = S, β = β, w = w, showProgress = true)
serialize(n_Js_min_data_guess, noised_guess_Js);
serialize(n_Qs_min_data_guess, noised_guess_Qs);
# -----------------------------------------------------------------------------

# ### Отрисовка результатов минимизации
# -----------------------------------------------------------------------------
a, b, c = minimization_draft(qₙ, noised_guess_Qs, Xₙ, N, noised_guess_Js, zoom = true,
                              annotate_string = params, zoom_chunk = 17//18)
plot(a, size = (800, 800))
plot(b)
plot(c)
withguessP = plot(a, b, c, layout = (3,1), size = (1000, 2400));
savefig(withguessP, n_guess_draft)
nothing; #hide
# -----------------------------------------------------------------------------
frames = collect(1:div(S, 500):S);
make_minimzation_gif(
    noised_guess_Js, noised_guess_Qs, qₙ, Xₙ,
    name = n_minim_gif_guess,
    frames_to_write = frames,
    β = β
)
run(`gif2mp4 $(n_minim_gif_guess) $(n_minim_mp4_guess)`)
run(`ffmpeg -i  $(n_minim_gif_guess) $(n_minim_avi_guess) -vcodec copy -acodec copy`)
# -----------------------------------------------------------------------------

# -----------------------------------------------------------------------------
# ------------------ Старт с нулевого приближения -----------------------------
# -----------------------------------------------------------------------------
q₀ = zero(Xₙ)
# -----------------------------------------------------------------------------
@time qs, noised_noguess_Js, noised_noguess_Qs = minimize(q₀, u₀, ulₘ, urₘ, Xₙ, N, Tₘ, M, ε, f1_data, f2_data,
                            S = S, β = β, w = w, showProgress = true)
serialize(n_Js_min_data, noised_noguess_Js);
serialize(n_Qs_min_data, noised_noguess_Qs);
# -----------------------------------------------------------------------------

# ### Эскиз процесса минимизации
# -----------------------------------------------------------------------------
a, b, c = minimization_draft(qₙ, noised_noguess_Qs, Xₙ, N, noised_noguess_Js, zoom = true,
                              annotate_string = params, zoom_chunk = 17//18)
plot(a, size = (800, 800))
plot(b)
plot(c)
withguessP = plot(a, b, c, layout = (3,1), size = (1000, 2400));
savefig(withguessP, n_noguess_draft)
nothing; #hide
frames = collect(1:div(S, 500):S);
make_minimzation_gif(
    noised_noguess_Js, noised_noguess_Qs, qₙ, Xₙ,
    name = n_minim_gif_noguess,
    frames_to_write = frames,
    β = β
)
run(`gif2mp4 $(n_minim_gif_noguess) $(n_minim_mp4_noguess)`)
run(`ffmpeg -i  $(n_minim_gif_noguess) $(n_minim_avi_noguess) -vcodec copy -acodec copy`)
# -----------------------------------------------------------------------------


text:
module MOCNeutronTransport

using EzXML
using HDF5
using Pkg.Artifacts
using Printf
using Reexport

using Colors: RGBA, N0f8
export RGBA, N0f8
using Dates: now, format
using LoggingExtras: TransformerLogger, global_logger

@reexport using LinearAlgebra
@reexport using StaticArrays
@reexport using Statistics

gmsh_dir = readdir(artifact"gmsh", join=true)[1]
gmsh_jl = joinpath(gmsh_dir, "lib", "gmsh.jl")
include(gmsh_jl)
export gmsh

include("common/Common.jl")
include("quadrature/Quadrature.jl")
include("geometry/Geometry.jl")
include("mesh/Mesh.jl")
#include("MPACT/MPACT.jl")
#include("plot/Plot.jl")
#include("gmsh/Gmsh.jl")
include("raytracing/Raytracing.jl")
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
using CPLEXCP
using JavaCall
using MathOptInterface
using ConstraintProgrammingExtensions

using Test

const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIB = MOI.Bridges
const CP = ConstraintProgrammingExtensions
const COIT = CP.Test

@testset "CPLEXCP" begin
    include("api_java.jl")
    include("MOI.jl")
end


text:
# Author: Ziyi Yin, Rafael Orozco, Philipp Witte, Mathias Louboutin, Gabrio Rizzuti, Felix J. Herrmann
# Date: April, 2020

using DrWatson
@quickactivate "ExtSrcImg"

include("gen_geometry.jl")

using PyPlot, FFTW, Images
using JUDI, JOLI, JUDI4Flux
using LinearAlgebra, PyPlot, Random, Statistics
using IterativeSolvers
using JLD
using MAT
using Distributions

function weight(n, d, xs, zs, nsrc,mode)

	x = zeros(n[1],1,nsrc);
    x .= reshape((collect(1:n[1]).-1)*d[1], :, 1)
	z = zeros(1,n[2],nsrc);
    z .= reshape((collect(1:n[2]).-1)*d[2], 1, :)
	W = zeros(Float32,n[1],n[2],nsrc)
	W .= sqrt.((x.-xs).^2f0.+(z.-zs).^2f0)

    return reshape(W,n[1],n[2],1,nsrc)
end

function CG(d_obs, q0, m0, W; maxiter=5,λ=1f2)
	# use CG to solve damped least square normal equations
	r1 = d_obs-𝒢(q0, m0);
	r = 𝒢_T(r1,m0)-λ*W.^2f0.*q0;
	p = r;
	r_vec = zeros(maxiter+1);
	r_vec[1] = norm(r1)
	println("residual: $(norm(r1)) before updating q")
	for i = 1:maxiter

		if norm(r1)<1f-3
			break
		end

		temp = sum(sum(r.*r,dims=1),dims=2);
		temp1 = 𝒢(p, m0);
		temp2 = 𝒢_T(temp1,m0)+λ*W.^2f0.*p;
		alpha = temp./sum(sum(p.*temp2,dims=1),dims=2);
		q0 = q0 + alpha.* p;
		r1 = r1 - alpha.* temp1;
		r = r - alpha .* temp2;
		beta = sum(sum(r.*r,dims=1),dims=2)./temp;
		p = r + beta.*p;
		r_vec[i+1]=norm(r1)
		println("residual: $(norm(r1)) after updating q")

	end
	return q0, r1, r_vec
end

# Objective function: estimate source and compute residual + function value
function objective(d_obs, q0, m0, W; λ=1f2)

	# Data residual and source residual
  	r_data   = d_obs-𝒢(q0, m0)
	r_source = q0

	# Function value
	fval = .5f0 *sum(r_data[:,:,1,1].^2) + 0.5f0*λ*sum((W.*r_source)[:,:,1,1].^2)

	print("Function value: ", fval, "\n")
	return fval
end

function GenSimSourceMulti(xsrc_index, zsrc_index, nsrc, n)
	weights = zeros(Float32, n[1], n[2], 1, nsrc)

	num_simultaneous_sources = size(xsrc_index)[1]
	#check if Simultaneous
	if num_simultaneous_sources > 1 && nsrc == 1
		#for Simultaneous we need to put all the weights in the same slice.
		for j=1:num_simultaneous_sources

        	weights[xsrc_index[j], zsrc_index[j], 1, 1] = 1f0
    	end
	else
		for j=1:nsrc
        	weights[xsrc_index[j], zsrc_index[j], 1, j] = 1#
    	end
	end

    return weights
end

function GenSimData(q₀; model=model)   # pass q₀ directly instead of z
    # Simultaneous observed data
    d_obs = F*vec(q₀)
    d_obs = reshape(d_obs, recGeometry.nt[1], nrec, 1, nsrc)  # reshape into tensor
    return d_obs
end

function make_calcifications(x=100,z=100,spread=80, n=20)
	#make n calcifications centered around x,y with a certain spread
	normal_two_dim = MvNormal(vec([x z]), Diagonal(spread*ones(2)))
	norm_coords =convert(Array{Int64,2},round.(rand(normal_two_dim, n)))
	norm_coords_x = norm_coords[1,:]
	norm_coords_z = norm_coords[2,:]
	for i in 1:n
		#use random size and location to place "microcalcification" of high velocity
		box_size_x = rand(0:2)
		box_size_z = rand(0:2)
		v[x+norm_coords_x[i]:x+norm_coords_x[i]+box_size_x, z+norm_coords_z[i]:z+norm_coords_z[i]+box_size_z] .= 3432.03 #in micrometers/microsecond
	end
end

############################################ Models ###############################################
M = load("../data/breast_ring.jld") #load breast velocity model

v = M["v"][1:601,1:601] * 1000 #bring v to micrometers / microseconds

size_x = size(v)[1]
size_z = size(v)[2]
n = (size_x,size_z)

# define smoothed velocity model for inversion NOTE: this is before we add micro calcifications
v0 = imfilter(deepcopy(v), Kernel.gaussian(10));

#add calcifications to forward model
Random.seed!(3) #set seed to get same microcalcifications.
make_calcifications(160,130,80,20)

#setup models
dx = 25 #in micrometers
dz = 25 #in micrometers
m  = 1f0 ./v.^2;
d  = (dx,dz);
o  = (0.,0.);

#the only difference between the forward and inversion model is the velocity model
m0       = 1f0 ./v0.^2; #s^2/m^2

model  = Model(n, d, o,  m, nb=80) #extra boundary layer size (nb) to avoid reflections artifacts
model0 = Model(n, d, o, m0, nb=80) #extra boundary layer size (nb) to avoid reflections artifacts

m = reshape(model.m, n[1], n[2], 1, 1)
m0 = reshape(model0.m, n[1], n[2], 1, 1)

##setup sensor sampling time values in microseconds
dtR = 0.00698
timeR = 3.4912718*3#in microseconds.

#source will have same time as sensors
timeS = timeR
dtS = dtR

#mode = "seismic"
mode = "medical"
include("gen_geometry.jl")

if mode == "seismic"

	####----new experiment with seismic like multishots
	xrec, yrec, zrec, nrec, nsrc, xsrc_index, zsrc_index, nsrc_simult = Gen_Exp_Seismic_Calc_Breast_Blood()

	#generate tensor with sources
	q           = GenSimSourceMulti(xsrc_index, zsrc_index, nsrc, model.n);

elseif mode == "medical"

	xrec, yrec, zrec, nrec, nsrc, xsrc_index, zsrc_index, nsrc_simult = Gen_Exp_Medical_Calc_Breast_Single()

	#generate tensor with sources
	q           = GenSimSourceMulti(xsrc_index, zsrc_index, nsrc, model.n);

	#get muting region
	mute_mask = GenMuteMedical_breast();
else
	error("Please choose a correct mode.")
end

recGeometry = Geometry(xrec, yrec, zrec; dt=dtR, t=timeR, nsrc=nsrc)

# Setup wavelet
f_dom = 3.25 #in [mHz]
wavelet = ricker_wavelet(timeS,dtS,f_dom)

######################################### Set up JUDI operators ###################################

# Set up info structure for linear operators
ntComp = get_computational_nt(recGeometry, model)
info   = Info(prod(model0.n), nsrc, ntComp)

# Write shots as array
opt = Options(return_array = true, subsampling_factor=5, isic=true)

# Modeling operators
Pr  = judiProjection(info, recGeometry)
F0  = judiModeling(info, model0; options=opt)	# modeling operator true model
F  = judiModeling(info, model; options=opt)	# modeling operator initial model
Pw = judiLRWF(info, wavelet)

# Combined operators
F = Pr*F*adjoint(Pw)
F0 = Pr*F0*adjoint(Pw)

# Extended modeling CNN layers
𝒢 = ExtendedQForward(F)
𝒢_T = ExtendedQAdjoint(F)

######################################### Run Forward simulations ###################################

# Generate shots
d_obs   = GenSimData(q);  #Data in true model

######################################### Generate 3 gradient images ###################################

#### 1.) RTM
J1=judiJacobian(F0,q)
d_pred = 𝒢(q,m0)
res_rtm = d_obs-d_pred;
grad_rtm=reshape(adjoint(J1)* vec(res_rtm),n[1],n[2])

#### 2.) Extended source gradient - inversion scheme
W = 0 #no focusing
λ = 0 #no focusing

#run cg to invert for estimated source
q_iter = 5
q_inv = zeros(Float32, n[1], n[2], 1, nsrc); #initialize source guess at zero
q_inv, r_end, rvec = CG(d_obs, q_inv, m0, W; maxiter=q_iter,λ=λ)

#get jacobian using estimated source and use it to get gradient
J = judiJacobian(F0, q_inv)
grad_extend=reshape(adjoint(J)* vec(r_end),n[1],n[2])

######################################### Visualize results   ###################################
#get extents of model for plotting
pixel_size_x = 0.2 #in mm
pixel_size_z = 0.2 #in mm
extentx = n[1]*pixel_size_x  #in mm
extentz = n[2]*pixel_size_z  #in mm

#show stacked gradient image

# Mute out
r = 30
x = range(-30, 30, length=size(grad_extend, 1))
z  = range(-30, 30, length=size(grad_extend, 1))

mask = zeros(size(grad_extend))

for i =1:size(grad_extend, 1)
	for j=1:size(grad_extend, 2)
		if (x[i]^2 + z[j]^2) < r^2
			mask[i, j] = 1
		end
	end
end

mask =imfilter(mask, Kernel.gaussian(5))
scale = .4 * 1e-1

fig=figure(figsize=(6.66,5));
imshow(-mask' .* grad_extend'/norm(grad_extend,Inf), cmap="Greys", vmin=-scale, vmax=scale , extent=(0, extentx, extentz,0))
title("Extended Source Image - Inversion - Stacked", fontsize=8)
xlabel("Lateral Position [mm]", fontsize=8)
ylabel("Horizontal Position [mm]", fontsize=8)
# colorbar(fraction=0.028, pad=0.04)

tight_layout()
savefig("breast_blood_stacked.png", dpi=300, format="png")
run(`mogrify -trim breast_blood_stacked.png`)

#Show solved sources and true sources
q_inv_stack=sum(q_inv,dims=4)[:,:,1,1];
figure(8);imshow(q_inv_stack',extent=(0,extentx,extentz,0));title("Extended Source from Inversion");
xlabel("Lateral position [mm]");ylabel("Depth [mm]");
colorbar(fraction=0.028, pad=0.04);
fig = PyPlot.gcf(); fig.set_size_inches(12.0, 7.0)

q_stack=sum(q,dims=4)[:,:,1,1];
figure(8);imshow(q_stack',extent=(0,extentx,extentz,0));title("True sources");
xlabel("Lateral position [mm]");ylabel("Depth [mm]");
colorbar(fraction=0.028, pad=0.04);
fig = PyPlot.gcf(); fig.set_size_inches(12.0, 7.0)

#show models used
figure(5);imshow(m0[:,:,1,1]',extent=(0,extentx,extentz,0));title("Background");
cb = colorbar(fraction=0.038, pad=0.04);
cb[:set_label]("medium velocity [m/s]")
xlabel("Lateral position [mm]");ylabel("Depth [mm]");
fig = PyPlot.gcf(); fig.set_size_inches(12.0, 7.0)

figure(7);imshow(m[:,:,1,1]',extent=(0,extentx,extentz,0));title("True");
xlabel("Lateral position [km]");ylabel("Depth [km]");colorbar(fraction=0.028, pad=0.04);
fig = PyPlot.gcf(); fig.set_size_inches(12.0, 7.0)

#Show Experiment Setup this is the figure in the paper
PyPlot.rc("font", family="serif"); PyPlot.rc("xtick", labelsize=8); PyPlot.rc("ytick", labelsize=8)
fig=figure(figsize=(6.66,5));
imshow(reshape(v/1000, n[1], n[2])',vmin=1.47,vmax=1.70,extent=(0,extentx,extentz,0));
cb = colorbar(fraction=0.040, pad=0.04);
#cb[:set_label]("Medium Velocity [m/s]", size="large")
PyPlot.scatter((xrec/d[1])*pixel_size_x,(zrec/d[2])*pixel_size_z,marker=".",color="red",label="Receivers");
PyPlot.scatter(xsrc_index*pixel_size_x,zsrc_index*pixel_size_z,label="Sources",marker="x",color="white");
title("True Model - Experiment Setup", fontsize=8);
legend(loc=1);
xlabel("Lateral Position [mm]", fontsize=8);ylabel("Horizontal Position [mm]", fontsize=8);
xlim((0, 120)); ylim((120,0));

tight_layout()
savefig("blood_breast_setup.png", dpi=300, format="png")

######################################### Save results   ########################################
@save "breastexperiment.jld" grad_extend q_inv