text:
using SupplyChainOptimization

function parse_simple_data(file_name)
    lines = readlines(file_name)
    m = match(r"(?<facility_count>\d+)\s+(?<customer_count>\d+)\s+0", lines[2])
    facility_count = parse(Int, m[:facility_count])
    customer_count = parse(Int, m[:customer_count])

    sc = SupplyChain()
    product = Product("p1")

    add_product!(sc, product)

    customers = []
    for f in 1:customer_count
        customer = Customer("c$f", Location(0, 0))
        add_customer!(sc, customer)
        add_demand!(sc, customer, product; demand=[1.0])
        push!(customers, customer)
    end

    for f in 1:facility_count
        data = [parse(Float64, c.match) for c in eachmatch(r"\d+", lines[2+f])]
        storage = Storage("s$f", Location(0, 0); fixed_cost=data[2], opening_cost=0.0, closing_cost=0.0, initial_opened=false)
        add_product!(storage, product; initial_inventory=customer_count, unit_handling_cost=0.0)
        add_storage!(sc, storage)

        for c in 1:customer_count
            lane = Lane(storage, customers[c]; unit_cost=data[2+c])
            add_lane!(sc, lane)
        end
    end

    return sc
end

function parse_orlib_data_uncap(file_name)
    numbers = split(read(file_name, String))
    number_index = 1

    facility_count = parse(Int, numbers[number_index])
    number_index += 1
    customer_count = parse(Int, numbers[number_index])
    number_index += 1

    sc = SupplyChain()
    product = Product("p1")

    add_product!(sc, product)

    facilities = []
    for f in 1:facility_count
        number_index += 1
        #println(numbers[number_index])
        storage = Storage("s$f", Location(0, 0); fixed_cost=parse(Float64, numbers[number_index]), opening_cost=0.0, closing_cost=0.0, initial_opened=false)
        number_index += 1
        add_product!(storage, product; initial_inventory=customer_count, unit_handling_cost=0.0)
        add_storage!(sc, storage)
        push!(facilities, storage)
    end

    for c in 1:customer_count
        #println(numbers[number_index])
        demand = parse(Float64, numbers[number_index])
        number_index += 1
        customer = Customer("c$c", Location(0, 0))
        add_customer!(sc, customer)
        add_demand!(sc, customer, product; demand=[1.0])

        for f in 1:facility_count
            #println(numbers[number_index])
            lane = Lane(facilities[f], customer; unit_cost=parse(Float64, numbers[number_index]))
            number_index += 1
            add_lane!(sc, lane)
        end
    end

    return sc
end

function parse_orlib_data_cap(file_name, capacity=nothing)
    if isnothing(capacity)
        numbers = split(read(file_name, String))
    else
        numbers = split(replace(read(file_name, String), "capacity"=>capacity))
    end
    number_index = 1

    facility_count = parse(Int, numbers[number_index])
    number_index += 1
    customer_count = parse(Int, numbers[number_index])
    number_index += 1

    sc = SupplyChain()
    product = Product("p1")

    add_product!(sc, product)

    facilities = []
    for f in 1:facility_count
        inventory = parse(Float64, numbers[number_index])
        number_index += 1
        storage = Storage("s$f", Location(0, 0); fixed_cost=parse(Float64, numbers[number_index]), opening_cost=0.0, closing_cost=0.0, initial_opened=false)
        number_index += 1
        add_product!(storage, product; initial_inventory=inventory, unit_handling_cost=0.0)
        add_storage!(sc, storage)
        push!(facilities, storage)
    end

    for c in 1:customer_count
        #println(numbers[number_index])
        demand = parse(Float64, numbers[number_index])
        number_index += 1
        customer = Customer("c$c", Location(0, 0))
        add_customer!(sc, customer)
        add_demand!(sc, customer, product; demand=[demand])

        for f in 1:facility_count
            #println(numbers[number_index])
            lane = Lane(facilities[f], customer, parse(Float64, numbers[number_index]) / customer.demand[product][1])
            number_index += 1
            add_lane!(sc, lane)
        end
    end

    return sc
end

@test begin
    sc = parse_simple_data(raw"..\data\BildeKrarup\B\B1.1")
    SupplyChainOptimization.optimize_network!(sc)
    true
end

@test begin
    sc = parse_orlib_data_uncap(raw"..\data\ORLIB\ORLIB-cap\40\cap41.txt")
    SupplyChainOptimization.optimize_network!(sc)
    true
end


text:
using CURAND
using Base.Test

# smoke tests for high-level API
rng = create_generator()
n = 10
n = 10
m = 5.
mf = Float32(5.)
sd = .1
sdf = Float32(.1)
lambda = .5

curand(rng, Float32, n)
curand(Float32, n)
curand(rng, Float64, n)
curand(Float64, n)
curand(n)

curandn(rng, Float32, n, mf, sdf)
curandn(Float32, n, mf, sdf)
curandn(rng, Float64, n, m, sd)
curandn(Float64, n, m, sd)
curandn(n, m, sd)

curand_logn(rng, Float32, n, mf, sdf)
curand_logn(Float32, n, mf, sdf)
curand_logn(rng, Float64, n, m, sd)
curand_logn(Float64, n, m, sd)
curand_logn(n, m, sd)

curand_poisson(rng, n, lambda)
curand_poisson(n, lambda)

destroy_generator(rng)

println("ok.")

ext: jl


text:
# AUTO GENERATED FILE - DO NOT EDIT

export ''_leafletmap

"""
    ''_leafletmap(;kwargs...)
    ''_leafletmap(children::Any;kwargs...)
    ''_leafletmap(children_maker::Function;kwargs...)

A LeafletMap component.

Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional)
- `id` (String; optional)
- `_drawnShapes` (Array; optional)
- `center` (optional): . center has the following type: lists containing elements 'lng', 'lat'.
Those elements have the following types:
  - `lng` (Real; optional)
  - `lat` (Real; optional)
- `className` (String; optional)
- `closePopupOnClick` (Bool; optional)
- `doubleClickZoom` (Bool; optional)
- `dragging` (Bool; optional)
- `editToolbar` (Bool; optional)
- `editToolbarControlsOptions` (optional): . editToolbarControlsOptions has the following type: lists containing elements 'position', 'drawMarker', 'drawCircleMarker', 'drawPolyline', 'drawRectangle', 'drawPolygon', 'drawCircle', 'drawText', 'editMode', 'dragMode', 'cutPolygon', 'removalMode', 'rotateMode', 'oneBlock'.
Those elements have the following types:
  - `position` (a value equal to: 'topleft', 'topright', 'bottomleft', 'bottomright'; optional)
  - `drawMarker` (Bool; optional)
  - `drawCircleMarker` (Bool; optional)
  - `drawPolyline` (Bool; optional)
  - `drawRectangle` (Bool; optional)
  - `drawPolygon` (Bool; optional)
  - `drawCircle` (Bool; optional)
  - `drawText` (Bool; optional)
  - `editMode` (Bool; optional)
  - `dragMode` (Bool; optional)
  - `cutPolygon` (Bool; optional)
  - `removalMode` (Bool; optional)
  - `rotateMode` (Bool; optional)
  - `oneBlock` (Bool; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `maxBounds` (optional): . maxBounds has the following type: lists containing elements 'minx', 'miny', 'maxx', 'maxy'.
Those elements have the following types:
  - `minx` (Real; optional)
  - `miny` (Real; optional)
  - `maxx` (Real; optional)
  - `maxy` (Real; optional)
- `maxDrawnShapes` (Real; optional)
- `maxZoom` (Real; optional)
- `minZoom` (Real; optional)
- `scrollWheelZoom` (Bool | a value equal to: 'center'; optional)
- `showMeasurements` (Bool; optional)
- `style` (Dict; optional)
- `wheelPxPerZoomLevel` (Real; optional)
- `zoom` (Real; optional)
- `zoomControl` (Bool; optional)
- `zoomDelta` (Real; optional)
"""
function ''_leafletmap(; kwargs...)
        available_props = Symbol[:children, :id, :_drawnShapes, :center, :className, :closePopupOnClick, :doubleClickZoom, :dragging, :editToolbar, :editToolbarControlsOptions, :loading_state, :maxBounds, :maxDrawnShapes, :maxZoom, :minZoom, :scrollWheelZoom, :showMeasurements, :style, :wheelPxPerZoomLevel, :zoom, :zoomControl, :zoomDelta]
        wild_props = Symbol[]
        return Component("''_leafletmap", "LeafletMap", "feffery_leaflet_components", available_props, wild_props; kwargs...)
end

''_leafletmap(children::Any; kwargs...) = ''_leafletmap(;kwargs..., children = children)
''_leafletmap(children_maker::Function; kwargs...) = ''_leafletmap(children_maker(); kwargs...)

max_stars_repo_name: CNFeffery/feffery-leaflet-components

max_issues_repo_name: CNFeffery/feffery-leaflet-components

max_forks_repo_name: CNFeffery/feffery-leaflet-components


text:
path = dirname(@__FILE__())

# Set up arguments
m = AnSchorfheide(testing = true)
m <= Setting(:date_forecast_start, quartertodate("2015-Q4"))

forecast_args = load("$path/../reference/forecast_args.jld2")
df = forecast_args["df"]
system = forecast_args["system"]

# Read expected output
smooth_out = load("$path/../reference/smooth_out.jld2")
exp_states = smooth_out["exp_states"]
exp_shocks = smooth_out["exp_shocks"]
exp_pseudo = smooth_out["exp_pseudo"]

# Smooth without drawing states
states = Dict{Symbol, Matrix{Float64}}()
shocks = Dict{Symbol, Matrix{Float64}}()
pseudo = Dict{Symbol, Matrix{Float64}}()

@testset "Test smoother without drawing states" begin
    for smoother in [:hamilton, :koopman, :carter_kohn, :durbin_koopman]
        m <= Setting(:forecast_smoother, smoother)

        states[smoother], shocks[smoother], pseudo[smoother] =
            smooth(m, df, system; draw_states = false)

        @test @test_matrix_approx_eq exp_states states[smoother]
        @test @test_matrix_approx_eq exp_shocks shocks[smoother]
        @test @test_matrix_approx_eq exp_pseudo pseudo[smoother]
    end
end

# Smooth, drawing states
for smoother in [:carter_kohn, :durbin_koopman]
    m <= Setting(:forecast_smoother, smoother)
    smooth(m, df, system; draw_states = true)
end

nothing


text:
using Compat.Distributed

if nprocs() > 1
    error("For the Timeout package's tests to be reliable, Julia must be started with" *
          "only a single worker process.")
end

# We need 2 because one will get terminated by the ptimeout test, which means there
# wouldn't be any remaining workers for the following test
addprocs(2)

using Timeout
using Compat.Test
using Compat.Dates

#@testset "No process available" begin
#    @test_throws ArgumentError ptimeout(()->1, 20, worker=2)
#end

@testset "Misspecified options" begin
    @test_throws ArgumentError ptimeout(()->8+1, 44.4, worker="nope")
    @test_throws ArgumentError ptimeout(()->sin(π), 3, worker=myid())
    @test_throws ArgumentError ptimeout(()->cos(π), 2, worker=2, poll=-1)
end

@test !ptimeout(Second(4), worker=2) do
    while true
        sin(rand())
    end
end

@test @ptimeout 3 4.5 begin
    for i = 1:10
        cos(i)
    end
end


text:
struct Results
    ψ::Array
    ω::Array
    hist::Array
    x::Array
    y::Array
    tfinal
    steps
    Re
end

ShowStreamlines(sol::SimpleNavierStokes.Results) = contour(sol.x,sol.y,reverse(reverse(sol.ψ,dims=1),dims=2),
          aspectratio=1,framestyle=:box,
          xlims=(sol.x[1],sol.x[end]),
          ylims=(sol.y[1],sol.y[end]),
          legend=:none,grid=:none)


text:
module Move_DoublePendulum

using  Modia3D
import Modia3D.ModiaMath

vmat1 = Modia3D.Material(color="Red", transparency=0.5)
vmat2 = Modia3D.Material(color="Blue", transparency=0.5)

@signal Sine(;A=1.0, freqHz = 1.0) begin
   y = ModiaMath.RealScalar(causality=ModiaMath.Output, numericType=ModiaMath.WR)
end
function Modia3D.computeSignal(signal::Sine, sim::ModiaMath.SimulationState)
    signal.y.value  = signal.A*sin(2*pi*signal.freqHz*sim.time)
end

@assembly DoublePendulum(;length=1.0, width=0.2) begin
   world = Modia3D.Object3D(Modia3D.CoordinateSystem(0.6))

   part1 = Modia3D.Part(Modia3D.Solid(Modia3D.SolidBox(length,width,width), "Aluminium", vmat1),
                        [ [-length/2,0,0],
                          [ length/2,0,0] ])
   part2 = Modia3D.Part(Modia3D.Solid(Modia3D.SolidBox(length,width,width), "Aluminium", vmat2),
                        [ [-length/2,0,0] ])

   rev1 = Modia3D.Revolute(world          , part1.frames[1])
   rev2 = Modia3D.Revolute(part1.frames[2], part2.frames[1])

   sig1 = Sine(A=1.5,freqHz=0.5)
   sig2 = Sine(A=2.5,freqHz=1.0)
   sine1     = Modia3D.SignalToFlangeAngle(sig1.y)
   sine2     = Modia3D.SignalToFlangeAngle(sig1.y)
   Modia3D.connect(sine1,   rev1)
   Modia3D.connect(sine2,   rev2)
end

doublePendulum = DoublePendulum(sceneOptions = Modia3D.SceneOptions(visualizeFrames=true,
                                                                    defaultFrameLength=0.3,
                                                                    enableContactDetection=false))
model = Modia3D.SimulationModel( doublePendulum, analysis=ModiaMath.KinematicAnalysis )
result = ModiaMath.simulate!(model, stopTime=3.0)

ModiaMath.plot(result, ("rev1.phi", "rev2.phi") )

println("... success of Move_DoublePendulum.jl!")
end


text:
using ApproxFunOrthogonalPolynomials, ApproxFunBase, SpecialFunctions, Test
    import ApproxFunBase: testbandedoperator

@testset "Laguerre and WeightedLagueree" begin
    @testset "Evaluation" begin
        f=Fun(Laguerre(0.), [1,2,3])
        @test f(0.1) ≈ 5.215

        f = Fun(Laguerre(0.1), ones(100))
        @test f(0.2) ≈ 8.840040924281498

        x = Fun(Laguerre())
        w = (1+x) * Fun(WeightedLaguerre(),[1.0])
        @test w(0.1) ≈ exp(-0.1)*(1+0.1)
        @test last(w) == 0.0
    end

    @testset "Derivative" begin
        f = Fun(Laguerre(0.1), ones(100))
        @test (Derivative(Laguerre(0.1)) * f)(0.2) ≈ -71.44556705957386

        f = Fun(Laguerre(0.2), ones(100))
        @test (Derivative(Laguerre(0.2)) * f)(0.3) ≈ -137.05785783078218
    end

    @testset "Conversion" begin
        f = Fun(Laguerre(0.2), ones(100))
        @test (Conversion(Laguerre(0.2), Laguerre(1.2)) * f)(0.1) ≈ f(0.1)
        @test (Conversion(Laguerre(0.2), Laguerre(2.2)) * f)(0.1) ≈ f(0.1)
    end

    @testset "Derivative" begin
        f = Fun(LaguerreWeight(0.0, Laguerre(0.1)), ones(100))
        @test f'(0.2) ≈ -65.7322962859456

        B = Evaluation(LaguerreWeight(0.0, Laguerre(0.1)), false)
        @test B*f ≈ 151.53223385808576

        x = Fun(Laguerre(0.0))
        S = WeightedLaguerre(0.0)
        D = Derivative(S)
        u = [ldirichlet(); D^2 - x] \ [airyai(0.0); 0.0]
        @test u(1.0) ≈ airyai(1.0)
    end

    @testset "Multiplication" begin
        w = Fun(WeightedLaguerre(), [1.0])
        @test sum(w) == 1
        t = Fun(identity, space(w))
        @test t(10.0) == 10.0
        M = Multiplication(t^2 + 1, space(w))
        @test (M \ w)(1.0) ≈ exp(-1)/2
    end

    @testset "Multiplication rangespace" begin
        x = Fun(Laguerre())
        w = Fun(WeightedLaguerre(-0.5),[1.0])
        M = Multiplication(w,space(1+x))
        @test rangespace(M) == LaguerreWeight(-0.5,Laguerre())
    end
end


text:
# test stacking

T = Float32
fs = 100.  # sampling rate
N = 2^14   # number of points
Nwin = 10
freqmin = 10. # minimum frequency
freqmax = 20. # maximum frequency
corners = 4   # number of corners
zerophase = true # do zerophase filtering
starttime = d2u(DateTime(Date(now()))) # starttime for SeisChanel / SeisData

@testset "stack" begin
    A = rand(T,N,Nwin)
    C = CorrData()
    C.corr = A
    C.fs = fs
    # corrleation every hour
    C.t = collect(0.:3600.:(Nwin-1)*3600.) .+ starttime

    # test daystack
    Cnew = stack(C)
    @test Cnew.t == [starttime]
    @test size(Cnew.corr) == (N,1)

    # test hour stack -> this does nothing
    Cnew = stack(C,interval=Hour(1))
    @test Cnew.t == C.t
    @test Cnew.corr == C.corr

    # test 2hour stack
    Cnew = stack(C,interval=Hour(2))
    @test Cnew.t == collect(0.:7200.:(Nwin-1)*3600.) .+ starttime
    @test size(Cnew.corr) == (N,5)

    # test 2hour stack w/ robuststack
    Cnew = stack(C,interval=Hour(2),stacktype=robuststack)
    @test Cnew.t == collect(0.:7200.:(Nwin-1)*3600.) .+ starttime
    @test size(Cnew.corr) == (N,5)

    # test allstack
    Cnew = stack(C,allstack=true)
    @test Cnew.t == [starttime]
    @test size(Cnew.corr) == (N,1)

    # test allstack w/ robuststack
    Cnew = stack(C,allstack=true,stacktype=robuststack)
    @test Cnew.t == [starttime]
    @test size(Cnew.corr) == (N,1)

    # test inplace
    Cnew = stack(C)
    stack!(C)
    @test Cnew.corr == C.corr
    @test Cnew.t == C.t
end

@testset "robust stack" begin
    A = rand(T,N,Nwin) .+ T.(sin.(range(0.,10π,length=N)))
    B = robuststack(A)
    @test size(B) == (N,1)
    @test eltype(B) == eltype(A)

    # test convergence
    B = robuststack(A,ϵ=Float32(1e-3))
    @test size(B) == (N,1)
    @test eltype(B) == eltype(A)

    # test maxiter
    B = robuststack(A,maxiter=15)
    @test size(B) == (N,1)
    @test eltype(B) == eltype(A)

    # test CorrData
    C = CorrData()
    C.corr = A
    C.fs = fs
    C.t = collect(0.:3600.:(Nwin-1)*3600.) .+ starttime
    Cnew = robuststack(C)
    @test Cnew.t == [starttime]
    @test size(Cnew.corr) == (N,1)
    @test eltype(Cnew.corr) == eltype(A)

    # test in-place
    robuststack!(C)
    @test isapprox(Cnew.corr,C.corr)
    @test Cnew.t == C.t
end

@testset "median mute" begin
    A = rand(T,N,Nwin)
    high = 2.
    low = 0.5
    # set two out of 10 to have large amplute
    A[:,1] .*= 1000
    # set two out of 10 to have small amplitude
    A[:,2] ./= 1000

    # test indexes for 0 < medianamp < high
    ind = SeisNoise.medianmuteind(A, high)
    @test length(ind) == Nwin - 1

    # test indexes for low < medianamp < high
    ind = SeisNoise.medianmuteind(A, high, low)
    @test length(ind) == Nwin - 2

    # test when low > high
    @test_throws AssertionError SeisNoise.medianmuteind(A, low, high)

    # test with Nans
    # set two columns to have NaNs
    A[rand(1:N,3),3] .= NaN
    ind = SeisNoise.medianmuteind(A, high, low)
    @test length(ind) == Nwin - 3

    # test with Infs
    # set two columns to have Infs
    A[rand(1:N,3),4] .= Inf
    ind = SeisNoise.medianmuteind(A, high, low)
    @test length(ind) == Nwin - 4

    # test with CorrData
    C = CorrData()
    C.corr = A
    C.fs = fs
    C.t = collect(0.:3600.:(Nwin-1)*3600.) .+ starttime

    # test allocating version
    Cnew = medianmute(deepcopy(C), high)
    @test Cnew != C
    @test size(Cnew.corr) == (N,7)
    @test length(Cnew.t) == 7

    # test with high and low thresholds
    Cnew = medianmute(deepcopy(C), high, low)
    @test Cnew != C
    @test size(Cnew.corr) == (N,6)
    @test length(Cnew.t) == 6

    # test in-place version
    medianmute!(C, high, low)
    @test C == Cnew
end

@testset "phase weighted stack " begin
    A = rand(T,N,Nwin) .+ T.(sin.(range(0.,10π,length=N)))
    B = pws(A)
    @test size(B) == (N,1)
    @test eltype(B) == eltype(A)

    # test power
    B = pws(A,pow=1)
    @test size(B) == (N,1)
    @test eltype(B) == eltype(A)

    # test CorrData
    C = CorrData()
    C.corr = A
    C.fs = fs
    C.t = collect(0.:3600.:(Nwin-1)*3600.) .+ starttime
    Cnew = pws(C)
    @test Cnew.t == [starttime]
    @test size(Cnew.corr) == (N,1)
    @test eltype(Cnew.corr) == eltype(A)

    # test in-place
    pws!(C,pow=2.)
    @test isapprox(Cnew.corr,C.corr)
    @test Cnew.t == C.t
end

@testset "remove nan" begin
    C = CorrData()
    A = rand(T,N,Nwin)
    C.corr = A
    C.fs = fs
    C.t = collect(0.:3600.:(Nwin-1)*3600.) .+ starttime

    # add nan to 4 columns
    C.corr[unique(rand(1:N,100)),1:3:end] .= NaN

    # remove nan
    Cnew = remove_nan(C)
    @test size(Cnew.corr) == (N,6)
    @test eltype(Cnew.corr) == eltype(A)
    @test length(Cnew.t) == 6

    # test inplace
    remove_nan!(C)
    @test Cnew.corr == C.corr
    @test Cnew.t == C.t

    # test no-nan
    C = CorrData()
    A = rand(T,N,Nwin)
    C.corr = A
    C.fs = fs
    C.t = collect(0.:3600.:(Nwin-1)*3600.) .+ starttime
    Cnew = remove_nan(C)
    @test Cnew.corr == C.corr
    @test Cnew.t == C.t

    # test all NaN
    C.corr[:,:] .= NaN
    @test_throws ArgumentError remove_nan(C)
end

@testset "smooth CorrData" begin
    C = CorrData()
    A = rand(T,10001,Nwin)
    C.corr = A
    # time vector is  [T00:00:00, T01:00:00, T02:00:00, ...]
    C.t = collect(0.:3600.:(Nwin-1)*3600.) .+ starttime

    # smooth every 2 hours
    Cnew = smooth(C,Hour(2))
    # test that 2nd correlation is mean of first two hours
    @test Cnew.corr[:,2] == vec(mean(C.corr[:,1:2],dims=2))

    # add two hours to last 5 correlations
    C.t[3:10] .+= 3600 * 3
    Cnew = smooth(C,Hour(3))
    # test that 2nd correlation is mean of first two hours
    @test Cnew.corr[:,2] == vec(mean(C.corr[:,1:2],dims=2))
    # test that 3rd correlation remains unchanged
    @test Cnew.corr[:,3] == C.corr[:,3]
    # test that 5th correlation is mean of 3:5
    @test Cnew.corr[:,5] == vec(mean(C.corr[:,3:5],dims=2))

    # randomize starttimes
    C.t = collect(0.:3600.:(Nwin-1)*3600.) .+ starttime
    ind = sample(1:Nwin,Nwin,replace=false)
    Cnew = smooth(C[ind],Hour(2))

    # test correlations are ordered
    @test all(diff(Cnew.t) .== 3600)
    @test Cnew.corr[:,2] == vec(mean(C.corr[:,1:2],dims=2))

    # change time vector to Days
    C.t = collect(0.:3600.:(Nwin-1)*3600.) .* 24 .+ starttime
    Cnew = smooth(C,Day(3))
    @test Cnew.corr[:,3] == vec(mean(C.corr[:,1:3],dims=2))

    # test inplace version
    smooth!(C,Day(3))
    @test Cnew == C
end

@testset "shorten CorrData" begin
    C = CorrData()
    A = rand(T,10001,Nwin)
    C.corr = A
    C.fs = fs
    C.maxlag = 50.
    C.t = collect(0.:3600.:(Nwin-1)*3600.) .+ starttime

    # test shorten to 40s
    newlag = 40.
    Cnew = shorten(C,newlag)
    @test size(Cnew.corr) == (Int(newlag * fs * 2 + 1),Nwin)
    @test eltype(Cnew.corr) == eltype(C.corr)
    @test Cnew.maxlag == newlag

    # test with Int
    CInt = shorten(C,Int(newlag))
    @test size(CInt.corr) == (Int(newlag * fs * 2 + 1),Nwin)
    @test eltype(CInt.corr) == eltype(C.corr)
    @test CInt.maxlag == newlag

    # test newlag > maxlag
    newlag = 100.
    @test_throws AssertionError shorten(C,newlag)

    # test newlag =  0
    newlag = 0.
    @test_throws AssertionError shorten(C,newlag)

    # test newlag <=  0
    newlag = -5.
    @test_throws AssertionError shorten(C,newlag)

    # test where newlag % (1/C.fs) != 0
    newlag = π
    Cpi = shorten(C,π)
    @test Cpi.maxlag == 3.14

    # test in-place
    newlag = 40.
    Cnew = shorten(C,newlag)
    shorten!(C,newlag)
    @test Cnew.corr == C.corr
    @test Cnew.maxlag == C.maxlag
end

@testset "Adaptive Filter" begin
    # test adaptive filter
    A = rand(T,N,Nwin)
    B = deepcopy(A)
    window=10.
    adaptive_filter!(B,window,fs)
    @test size(B) == size(A)
    @test eltype(B) == eltype(B)

    # test overlap
    B = deepcopy(A)
    @test_throws ArgumentError adaptive_filter!(B,window,fs,overlap=0.)

    @test_throws ArgumentError adaptive_filter!(B,window,fs,overlap=1.)

    # test 1D vector
    v = rand(T,N)
    vcopy = deepcopy(v)
    adaptive_filter!(vcopy,window,fs)
    @test v == vcopy

    # test CorrData
    C = CorrData()
    A = rand(T,N,Nwin)
    C.corr = A
    C.fs = fs
    Cnew = adaptive_filter(C,window)
    @test size(Cnew.corr) == size(C.corr)
    @test eltype(Cnew.corr) == eltype(C.corr)

    # test CorrData in-place
    adaptive_filter!(C,window)
    @test Cnew.corr == C.corr

    # test ACF_kernel
    A = rand(T,N,Nwin)
    B = SeisNoise.ACF_kernel(A)
    @test size(B) == size(A)
    @test eltype(B) == eltype(A)

end

@testset "robustpws" begin
    A = rand(T,N,Nwin) .+ T.(sin.(range(0.,10π,length=N)))
    B = robustpws(A)
    @test size(B) == (N,1)
    @test eltype(B) == eltype(A)

    # test CorrData
    C = CorrData()
    C.corr = A
    C.fs = fs
    C.t = collect(0.:3600.:(Nwin-1)*3600.) .+ starttime
    Cnew = robustpws(C)
    @test Cnew.t == [starttime]
    @test size(Cnew.corr) == (N,1)
    @test eltype(Cnew.corr) == eltype(A)

    # test in-place
    robustpws!(C)
    @test Cnew.corr == C.corr
    @test Cnew.t == C.t
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#    Tests of table titles.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Titles" begin
    title = "This is a very very long title that will be displayed above the table."

    expected = """
This is a very very long title that will be displayed above the table.
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│      1 │  false │    1.0 │      1 │
│      2 │   true │    2.0 │      2 │
│      3 │  false │    3.0 │      3 │
│      4 │   true │    4.0 │      4 │
│      5 │  false │    5.0 │      5 │
│      6 │   true │    6.0 │      6 │
└────────┴────────┴────────┴────────┘
"""
    result = pretty_table(String, data,
                          title = title)
    @test result == expected

    expected = """
This is a very very long title that …
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│      1 │  false │    1.0 │      1 │
│      2 │   true │    2.0 │      2 │
│      3 │  false │    3.0 │      3 │
│      4 │   true │    4.0 │      4 │
│      5 │  false │    5.0 │      5 │
│      6 │   true │    6.0 │      6 │
└────────┴────────┴────────┴────────┘
"""
    result = pretty_table(String, data,
                          title = title,
                          title_same_width_as_table = true)
    @test result == expected

    expected = """
This is a very very long title that
will be displayed above the table.
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│      1 │  false │    1.0 │      1 │
│      2 │   true │    2.0 │      2 │
│      3 │  false │    3.0 │      3 │
│      4 │   true │    4.0 │      4 │
│      5 │  false │    5.0 │      5 │
│      6 │   true │    6.0 │      6 │
└────────┴────────┴────────┴────────┘
"""
    result = pretty_table(String, data,
                          title = title,
                          title_autowrap = true,
                          title_same_width_as_table = true)
    @test result == expected

    expected = """
 This is a very very long title that
 will be displayed above the table.
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│      1 │  false │    1.0 │      1 │
│      2 │   true │    2.0 │      2 │
│      3 │  false │    3.0 │      3 │
│      4 │   true │    4.0 │      4 │
│      5 │  false │    5.0 │      5 │
│      6 │   true │    6.0 │      6 │
└────────┴────────┴────────┴────────┘
"""
    result = pretty_table(String, data,
                          title = title,
                          title_alignment = :c,
                          title_autowrap = true,
                          title_same_width_as_table = true)
    @test result == expected

    expected = """
  This is a very very long title that
   will be displayed above the table.
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│      1 │  false │    1.0 │      1 │
│      2 │   true │    2.0 │      2 │
│      3 │  false │    3.0 │      3 │
│      4 │   true │    4.0 │      4 │
│      5 │  false │    5.0 │      5 │
│      6 │   true │    6.0 │      6 │
└────────┴────────┴────────┴────────┘
"""
    result = pretty_table(String, data,
                          title = title,
                          title_alignment = :r,
                          title_autowrap = true,
                          title_same_width_as_table = true)
    @test result == expected
end

@testset "Title with table cropping" begin

    # Cropping
    # --------------------------------------------------------------------------

    expected = """
This is a long long long long long…
That has two lines.
┌────────┬────────┬────────┬───────
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. ⋯
├────────┼────────┼────────┼───────
│      1 │  false │    1.0 │      ⋯
│      2 │   true │    2.0 │      ⋯
│      3 │  false │    3.0 │      ⋯
│   ⋮    │   ⋮    │   ⋮    │   ⋮  ⋱
└────────┴────────┴────────┴───────
        1 column and 3 rows omitted
"""
    result = pretty_table(String, data,
                          crop = :both,
                          display_size = (13,35),
                          title = "This is a long long long long long long title\nThat has two lines.")
    @test result == expected

    expected = """
This is a long long long long long
long title
That has two lines.
┌────────┬────────┬────────┬───────
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. ⋯
├────────┼────────┼────────┼───────
│      1 │  false │    1.0 │      ⋯
│      2 │   true │    2.0 │      ⋯
│   ⋮    │   ⋮    │   ⋮    │   ⋮  ⋱
└────────┴────────┴────────┴───────
        1 column and 4 rows omitted
"""
    result = pretty_table(String, data,
                          crop = :both,
                          display_size = (13,35),
                          title = "This is a long long long long long long title\nThat has two lines.",
                          title_autowrap = true)
    @test result == expected
end


text:
@recipe function plot(match_wave::MatchPlaneWaveMode{T}) where T <: AbstractFloat

    dx = match_wave.x_match[2] - match_wave.x_match[1]
    max_x = match_wave.x_match[end]
    x = [match_wave.x_match; (max_x:dx:(2*max_x)) .+ dx]
    @series (x, match_wave)
end

@recipe function plot(x::AbstractVector{T}, match_wave::MatchPlaneWaveMode{T};
        basis_order = match_wave.PlaneWaveModes[1].basis_order,
        hankel_indexes = 0:basis_order,
        apply = real, match_region = true) where T <: AbstractFloat

    ho = basis_order
    i1 = findmin(abs.(match_wave.discrete_wave.x .- match_wave.x_match[1]))[2]
    i2 = findmin(abs.(match_wave.discrete_wave.x .- match_wave.x_match[end]))[2]

    max_amp = maximum(apply.(match_wave.discrete_wave.amplitudes[i1:i2,(hankel_indexes) .+ (ho+1),:]))
    min_amp = minimum(apply.(match_wave.discrete_wave.amplitudes[i1:i2,(hankel_indexes) .+ (ho+1),:]))

    max_amp = (max_amp > 0) ? 1.1*max_amp : 0.0
    min_amp = (min_amp < 0) ? 1.1*min_amp : 0.0

    if match_region == true
        @series begin
            label --> "match region"
            fillalpha --> 0.3
            fill --> (0,:orange)
            line --> 0
            (x1 -> x1, y->max_amp, match_wave.x_match[1],  match_wave.x_match[end])
        end

        @series begin
            label --> ""
            fillalpha --> 0.3
            fill --> (0,:orange)
            line --> 0
            (x1 -> x1, y->min_amp, match_wave.x_match[1],  match_wave.x_match[end])
        end
    end

    @series begin
        apply --> apply
        hankel_indexes --> hankel_indexes
        # markeralpha --> 0.2
        markersize --> 4.0
        match_wave.discrete_wave
    end

    @series begin
        apply --> apply
        hankel_indexes --> hankel_indexes
        (x, match_wave.PlaneWaveModes)
    end

end

max_stars_repo_name: arturgower/EffectiveWaves


text:
# This file is a part of JuliaFEM/TetGen.
# License is MIT: see https://github.com/JuliaFEM/TetGen.jl/blob/master/LICENSE.md

using BinDeps

@BinDeps.setup
tetgen = "tetgen1.4.3"
libtet = library_dependency("libtet", aliases = ["libtet.so"])
sources = provides(Sources,URI("http://tetgen.org/files/$tetgen.tar.gz"),libtet,unpacked_dir=tetgen)
builddir = joinpath(BinDeps.builddir(libtet),libtet.name)
srcdir = joinpath(BinDeps.srcdir(libtet),tetgen)
libdir = joinpath(BinDeps.depsdir(libtet),"usr","lib")
includedir = joinpath(BinDeps.depsdir(libtet),"usr","include")

provides(SimpleBuild,
         (@build_steps begin
            GetSources(libtet)
            CreateDirectory(builddir)
            CreateDirectory(libdir)
            CreateDirectory(includedir)
            @build_steps begin
              ChangeDirectory(builddir)
              `g++ -fPIC -O0 -c $srcdir/predicates.cxx`
              `g++ -fPIC -DTETLIBRARY -c $srcdir/tetgen.cxx`
              `g++ -shared -fPIC tetgen.o predicates.o -o libtet.so`
              `cp libtet.so $libdir`
              `cp $srcdir/tetgen.h $includedir`
            end
          end),
         libtet, os = :Unix)

@BinDeps.install Dict( :libtet => :libtet)


text:
# **This example is out of date and won't work.**
using Pkg #src
Pkg.activate(@__DIR__) #src

# # ZDT3 Problem

# ## Setup
# Install the test problem suite:

#md # ```julia
#md # using Pkg
#md # Pkg.activate(tempname())
#md # Pkg.develop(url="https://github.com/manuelbb-upb/MultiObjectiveProblems.jl")
#md # using MultiObjectiveProblems
#md # ```

using MultiObjectiveProblems; #hide

# Import other dependencies:
using CairoMakie
using Morbit

# Retrieve test problem and define a `MixedMOP`
test_problem = ZDT3(2);
box = constraints(test_problem);

objectives = get_objectives(test_problem)
x₀ = get_random_point(test_problem)

mop = MixedMOP( box.lb, box.ub );
objf_cfg = ExactConfig()
for objf ∈ objectives
    add_objective!(mop, objf, objf_cfg)
end

# ## Run
# Run optimization and plot:
I = get_ideal_point(test_problem)
ac = AlgoConfig(; descent_method = :ps, reference_point = I )

x, fx, id = optimize( mop, x₀; algo_config = ac);

pset = get_pareto_set(test_problem)
PSx,PSy = get_scatter_points(pset, 100)

#~ scatter Pareto set points in grey
fig, ax, _ = scatter( PSx, PSy;
    figure = (resolution = (600, 650),),
)

#~ set axis limits to whole feasible set
xlims!(ax, (box.lb[1] .- .2, box.ub[1] .+ .2) )
ylims!(ax, (box.lb[2] .- .2, box.ub[2] .+ .2) )

#~ final iterate in red
scatter!(Tuple(x); color = :red)

fig


text:
@testset "begin and end tests" begin
    df = DataFrame([(i, j) for i in 1:3, j in 1:4])
    @test df[begin, begin] == df[1, 1]
    @test df[begin, end] == df[1, 4]
    @test df[end, begin] == df[3, 1]
    @test df[end, end] == df[3, 4]

    @test df[Not(begin), Not(begin)] == df[Not(1), Not(1)]
    @test df[Not(begin), Not(end)] == df[Not(1), Not(4)]
    @test df[Not(end), Not(begin)] == df[Not(3), Not(1)]
    @test df[Not(end), Not(end)] == df[Not(3), Not(4)]

    df[begin, begin] = (101, 101)
    @test df[begin, begin] == (101, 101)
    df[begin, end] = (101, 104)
    @test df[begin, end] == (101, 104)
    df[end, begin] = (103, 101)
    @test df[end, begin] == (103, 101)
    df[end, end] = (103, 104)
    @test df[end, end] == (103, 104)

    df[!, begin] .= [1,2,3]
    @test df[:, 1] == [1,2,3]
    df[!, end] .= [11,12,13]
    @test df[:, 4] == [11,12,13]

    @test df[begin:end, [begin, end]] == df[:, [1,4]]
    df[begin:end, [begin, end]] .= [111, 222, 333]
    @test df.x1 == df.x4 == [111, 222, 333]
    @test df[[begin, end], [begin, end]] == df[[1,3], [1,4]]
    df[[begin, end], [begin, end]] .= 1000
    @test df.x1 == df.x4 == [1000, 222, 1000]

    @test eachcol(df)[begin] == df[!, begin]
    @test eachcol(df)[end] == df[!, end]
end


text:
#=

Log-likelihood for ODEs, integrating along the tree
using flow algorithm

Ignacio Quintero Mächler

t(-_-t)

Created 23 10 2019

=#

"""
    make_Et(Et::Function,
            p0::Array{Float64,1},
            u0::Array{Float64,1},
            ts::Array{Float64,1},
            ti::Float64,

Make Extinction through time, `E(t)`.
"""
function make_Et(Et::Function,
                 p0::Array{Float64,1},
                 u0::Array{Float64,1},
                 ts::Array{Float64,1},
                 ti::Float64,
                 tf::Float64)

  prob = ODEProblem(Et, u0, (ti,tf), p0)

  int = init(prob,
             Tsit5(),
             saveat          = ts,
             save_everystep  = false,
             calck           = false,
             force_dtmin     = true,
             save_start      = true,
             initialize_save = false,
             maxiters        = 100_000_000,
             verbose         = false)

  function f(p::Array{Float64,1})
    int.p = p
    reinit!(int)
    solve!(int).u::Array{Array{Float64,1},1}
  end

  return f
end

"""
    make_Gt(At::Function,
            p0::Array{Array{Float64,1},1},
            u0::Array{Float64,2},
            ts::Array{Float64,1},
            ti::Float64,
            tf::Float64)

Make flow equation for likelihoods through time, `G(t)`.
"""
function make_Gt(At::Function,
                 p0::Array{Array{Float64,1},1},
                 u0::Array{Float64,2},
                 ts::Array{Float64,1},
                 ti::Float64,
                 tf::Float64)

  prob = ODEProblem(At, u0, (ti,tf), p0)

  int = init(prob,
             Tsit5(),
             saveat          = ts,
             save_everystep  = false,
             calck           = false,
             force_dtmin     = true,
             save_start      = true,
             initialize_save = false,
             maxiters        = 100_000_000,
             verbose         = false)

  function f(p::Array{Array{Float64,1},1})
    int.p = p
    reinit!(int)
    solve!(int).u::Array{Array{Float64,2},1}
  end

  return f
end

"""
    make_loglik(Gt     ::Function,
                Et     ::Function,
                X      ::Array{Array{Float64,1},1},
                triads ::Array{Array{Int64,1},1},
                lbts   ::Array{Int64,2},
                bts    ::Array{Float64,1},
                ns     ::Int64,
                ned    ::Int64,
                nets   ::Int64,
                ntip   ::Int64,
                λevent!::Function,
                rootll ::Function)

Make log-likelihood function using the flow algorithm.
"""
function make_loglik(Gt     ::Function,
                     Et     ::Function,
                     X      ::Array{Array{Float64,1},1},
                     triads ::Array{Array{Int64,1},1},
                     lbts   ::Array{Int64,2},
                     bts    ::Array{Float64,1},
                     ns     ::Int64,
                     ned    ::Int64,
                     nets   ::Int64,
                     λevent!::Function,
                     rootll ::Function)

  # preallocate arrays
  X0  = Array{Float64,1}(undef,ns)
  Xp1 = Array{Float64,1}(undef,ns)
  Xp2 = Array{Float64,1}(undef,ns)
  Xr  = Array{Float64,1}(undef,ns)
  wg  = Array{Float64,1}(undef,ns)
  Er  = Array{Float64,1}(undef,ns)

  rtime  = bts[end]
  netsm1 = nets - 1
  nbts   = length(bts)

  # start ll function for parameters
  function f(p::Array{Float64,1})

    @inbounds begin

      Ets = Et(p)

      # push parameters as the last vector in Ets
      push!(Ets, p)

      # estimate `Gts` according to `Ets` and `p`
      Gts = Gt(Ets)

      ll = 0.0

      if rank(Gts[nbts]) == ns

        for i in Base.OneTo(nbts)
          if rank(Gts[i]) != ns
            @show Gts[nbts]
            @show Gts[i]
          end

        end

        for trio in triads

          pr, d1, d2 = trio

          ## first daughter
          X0 = Gts[lbts[d1,2]]\X[d1]
          mul!(Xp1, Gts[lbts[d1,1]], X0)

          check_negs(Xp1, ns) && return -Inf

          ## second daughter
          X0 = Gts[lbts[d2,2]]\X[d2]
          mul!(Xp2, Gts[lbts[d2,1]], X0)

          check_negs(Xp2, ns) && return -Inf

          λtime = bts[lbts[d1,1]]

          λevent!(λtime, X, Xp1, Xp2, p, pr)

          if pr != ned
            mdlus  = sum(X[pr])
            X[pr] /= mdlus
            ll += log(mdlus)
          end
        end
      else

        for trio in triads

          pr, d1, d2 = trio

          ## first daughter
          ldiv!(X0, qr!(Gts[lbts[d1,2]], Val(true)), X[d1])
          mul!(Xp1, Gts[lbts[d1,1]], X0)

          check_negs(Xp1, ns) && return -Inf

          ## second daughter
          ldiv!(X0, qr!(Gts[lbts[d2,2]], Val(true)), X[d2])
          mul!(Xp2, Gts[lbts[d2,1]], X0)

          check_negs(Xp2, ns) && return -Inf

          λtime = bts[lbts[d1,1]]

          λevent!(λtime, X, Xp1, Xp2, p, pr)

          if pr != ned
            mdlus  = sum(X[pr])
            X[pr] /= mdlus
            ll += log(mdlus)
          end
        end

      end

      # assign root likelihoods
      @simd for i in Base.OneTo(ns)
        Xr[i] = X[ned][i]
      end

      # assign root extinction
      @simd for i in Base.OneTo(ns)
        Er[i] = Ets[netsm1][i]
      end

      # estimate weights
      normbysum!(Xr, wg, ns)

      ll += log(rootll(rtime, Xr, Er, wg, p)::Float64)
    end

    return ll
  end

  return f
end

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
export binary_quadratic_form, can_solve, discriminant,
       fundamental_discriminant, isdiscriminant, QuadBin,
       isfundamental_discriminant, prime_form, prime_power_form, cycle,
       isindefinite, ispositive_definite, isnegative_definite, isreduced

import Nemo: can_solve

###############################################################################
#
#   Constructor
#
###############################################################################

@doc Markdown.doc"""
    binary_quadratic_form(a::T, b::T, c::T) where {T <: RingElement}

Constructs the binary quadratic form $ax^2 + bxy + cy^2$.
"""
binary_quadratic_form(a::T, b::T, c::T) where {T <: RingElem} = QuadBin(parent(a), a, b, c)

binary_quadratic_form(a::Integer, b::Integer, c::Integer) = binary_quadratic_form(FlintZZ(a), FlintZZ(b), FlintZZ(c))

@doc Markdown.doc"""
    binary_quadratic_form(R::Ring, a, b, c)

Constructs the binary quadratic form $ax^2 + bxy + cy^2$ over the ring $R$.
The elements `a`, `b` and `c` are coerced into the ring $R$.
"""
binary_quadratic_form(R::Ring, a, b, c) = binary_quadratic_form(R(a), R(b), R(c))

################################################################################
#
#  Printing
#
################################################################################

function show(io::IO, ::MIME"text/plain", f::QuadBin)
  _show(io, f, get(io, :compact, false))
end

function show(io::IO, f::QuadBin)
  _show(io, f, true)
end

function _show(io::IO, f::QuadBin, compact = false)
  if !compact
    print(io, "Binary quadratic form with equation\n")
  end
  sum = Expr(:call, :+)
  a = f[1]
  if !iszero(a)
    if isone(a)
      push!(sum.args, Expr(:call, :*, Expr(:call, :^, :x, 2)))
    else
      push!(sum.args, Expr(:call, :*, AbstractAlgebra.expressify(a, context = io), Expr(:call, :^, :x, 2)))
    end
  end
  b = f[2]
  if !iszero(b)
    if isone(b)
      push!(sum.args, Expr(:call, :*, :x, :y))
    else
      push!(sum.args, Expr(:call, :*, AbstractAlgebra.expressify(b, context = io), :x, :y))
    end
  end
  c = f[3]
  if !iszero(c)
    if isone(c)
      push!(sum.args, Expr(:call, :*, Expr(:call, :^, :y, 2)))
    else
      push!(sum.args, Expr(:call, :*, AbstractAlgebra.expressify(c, context = io), Expr(:call, :^, :y, 2)))
    end
  end
  print(io, AbstractAlgebra.expr_to_string(AbstractAlgebra.canonicalize(sum)))
end

################################################################################
#
#   Base ring
#
################################################################################

base_ring(f::QuadBin{T}) where {T} = f.base_ring::parent_type(T)

###############################################################################
#
#   Coefficients
#
###############################################################################

function getindex(f::QuadBin, i::Int)
  if i == 1
    return f.a
  elseif i == 2
    return f.b
  elseif i == 3
    return f.c
  else
    throw(error("Index must be between 1 and 3"))
  end
end

################################################################################
#
#   Scalar multiplication
#
################################################################################

function Base.:(*)(c::T, f::QuadBin{T}) where T <: RingElem
  return binary_quadratic_form(c * f[1], c * f[2], c * f[3])
end

function Base.:(*)(c::fmpz, f::QuadBin{T}) where T <: RingElem
  return binary_quadratic_form(c * f[1], c * f[2], c * f[3])
end

function Base.:(*)(c::Integer, f::QuadBin)
  return binary_quadratic_form(c * f[1], c * f[2], c * f[3])
end

function divexact(f::QuadBin{T}, c::T) where T <: RingElem
  return binary_quadratic_form(divexact(f[1], c), divexact(f[2], c),
                                                  divexact(f[3], c))
end

###############################################################################
#
#  Evaluation
#
###############################################################################

(f::QuadBin)(x, y) = f[1]*x^2 + f[2] * x * y + f[3] * y^2

################################################################################
#
#  Representation
#
################################################################################

@doc Markdown.doc"""
     can_solve_with_solution(f::QuadBin, n::Int64)
For a binary quadratic form $f$ with negative discriminant and an integer $n$,
returns tuple (`true`, (`x`,`y`)) if $f$(`x`,`y`) = $n$ for integers `x`,`y`.
If no such integers exist, return (`false`,(0,0))
"""
function can_solve(f::QuadBin, n::Int64)
    if discriminant(f) >= 0
        @error("$f has non-negative discriminant")
    end
    for y in 1:Int64(floor(sqrt(Int64(4*f[1]*n)//abs(Int64(discriminant(f))))))
        #now f(x,y) quadratic in one variable -> use quadratic formula
        aq = f[1]
        bq = f[2] * y
        cq = f[3] * y^2 - n
        d = bq^2 - 4*aq*cq
        if issquare(d)
            if divides(-bq + sqrt(d), 2*aq)[1]
                return(true, (divexact(-bq + sqrt(d), 2*aq), ZZ(y)))
            end
            if divides(-bq - sqrt(d), 2*aq)[1]
                return(true, (divexact(-bq - sqrt(d), 2*aq), ZZ(y)))
            end
        end
    end
    return (false, (FlintZZ(0),ZZ(0)))
end

###############################################################################
#
#   Discriminant
#
###############################################################################

@doc Markdown.doc"""
    discriminant(f::QuadBin) -> RingElem

Return the discriminant of `f = [a, b, c]`, that is, `b^2 - 4ac`.
"""
function discriminant(f::QuadBin)
  if isdefined(f, :disc)
    return f.disc
  else
    d = f[2]^2 - 4 * f[1] * f[3]
    f.disc = d
    return d
  end
end

################################################################################
#
#  Discriminant of integral binary quadratic form
#
################################################################################

@doc Markdown.doc"""
    isdiscriminant(D)

Returns `true` if $D$ is the discriminant of an integral binary quadratic form,
otherwise returns `false`.
"""
function isdiscriminant(D::Union{Integer, fmpz})
  if D == 0
    return false
  end
  m = mod(D, 4)
  if m == 0 || m == 1
    return true
  end
  return false
end

@doc Markdown.doc"""
    isfundamental_discriminant(D)

Returns `true` if $D$ is a fundamental discriminant otherwise returns `false`.
"""
function isfundamental_discriminant(D::Union{Integer, fmpz})
  m = mod(D, 4)
  if m == 1 && issquarefree(D)
    return true
  end
  if m == 0
    h = divexact(D, 4)
    c = mod(h,4)
    if (c == 2 || c == 3) && issquarefree(h)
      return true
    end
  end
  return false
end

@doc Markdown.doc"""
     conductor(D) -> fmpz
Returns the conductor of the discriminant $D$, that is, the largest
positive integer $c$ such that $\frac{D}{c^2}$ is a discriminant.
"""
function conductor(D::Union{Integer, fmpz})
  @req isdiscriminant(D) "Value ($D) not a discriminant"
  d = divexact(D, fundamental_discriminant(D))
  return isqrt(d)
end

function fundamental_discriminant(D::Union{Integer, fmpz})
  fac = factor(D)
  sqf = one(FlintZZ)
  for (p, e) in fac
    if isodd(e)
      sqf = sqf * p
    end
  end
  # sqf = is the squarefree-part, so D = sqf * square and sqf square-free
  if mod(sign(D) * sqf, 4) == 1
    return sign(D) * sqf
  else
    return sign(D) * sqf * 4
  end
end

###############################################################################
#
#   Equality
#
###############################################################################

function ==(f1::QuadBin, f2::QuadBin)
  if base_ring(f1) != base_ring(f2)
    return false
  end
  return f1[1] == f2[1] && f1[2] == f2[2] && f1[3] == f2[3]
end

###############################################################################
#
#   Arithmetic
#
###############################################################################

conjugate(f::QuadBin) = binary_quadratic_form(f[1], -f[2], f[3])

-(f::QuadBin) = binary_quadratic_form(-f[1], -f[2], -f[3])

Generic.content(f::QuadBin{fmpz}) = gcd([f[1], f[2], f[3]])

isindefinite(f::QuadBin) = discriminant(f) > 0 ? true : false

isnegative_definite(f::QuadBin) = (discriminant(f) < 0 && f[1] < 0)

ispositive_definite(f::QuadBin) = (discriminant(f) < 0 && f[1] > 0)

Base.iszero(f::QuadBin) = (f[1] == 0 && f[2] == 0 && f[3] == 0)

isprimitive(f::QuadBin) = (isone(content(f)))

isreducible(f::QuadBin) = issquare(discriminant(f))

###############################################################################
#
#   Composition
#
###############################################################################

@doc Markdown.doc"""
    compose(f1::QuadBin, f2::QuadBin)

Returns the composition of the binary quadratic forms $f_1$ and $f_2$. The
result is not reduced, uses Dirichlet Composition.
"""
function compose(f1::QuadBin{fmpz}, f2::QuadBin{fmpz})
  #discriminants have to match
  D = discriminant(f1)
  D2 = discriminant(f2)
  if D != D2
    error("discriminants do not match")
  end
  (h, n_1, n_2) = gcdx(f1[1], f2[1])
  (e, n_3, n3) = gcdx(h, divexact(f1[2]+f2[2], 2))

  (n1, n2) = (n_3 * n_1, n_3 * n_2)
  B = n1 * divexact(f1[1] * f2[2], e) + n2 * divexact(f1[2] * f2[1], e) + n3 * divexact(f1[2] * f2[2] + D, 2*e)
  return binary_quadratic_form(divexact(f1[1]*f2[1], e^2), B, divexact(e^2*(B^2-D), 4*f1[1]*f2[1]))
end

###############################################################################
#
#   Prime Forms
#
###############################################################################

function _sqrtmod4P(d::fmpz, p::fmpz)
    if jacobi_symbol(mod(d, p), p) == -1
        error("$d is no square modulo $p")
    end
    if p == 2
        if iseven(d)
            return 2 * mod(divexact(d, 4), 2)
        else
            return 1
        end
    else
        r = sqrtmod(d, p)
        if mod(r,2) == mod(d,2)
            return r
        else
            return p-r
        end
    end
end

function _number_of_primeforms(d::fmpz, p::fmpz)
    return jacobi_symbol(mod(d, p), p) + 1
end

@doc Markdown.doc"""
     prime_form(d::fmpz, p::fmpz)
Returns an integral binary quadratic form of discriminant $d$ and leading coefficient
$p$ where $p$ is a prime number.
"""
function prime_form(d::fmpz, p::fmpz)
    if !isdiscriminant(d)
        error("$d is no discriminant")
    end
    if _number_of_primeforms(d, p) == 0
        error("prime form does not exist")
    end
    b = _sqrtmod4P(d, p)
    return binary_quadratic_form(p, b, divexact(b^2-d, 4*p))
end

function _sqrtmod4PE(d::fmpz, p::fmpz, e)
    if mod(d,p) == 0
        n = 1
        while (2*n + 2  <= e && mod(p^(2n+2), d) == 0)
            n = n + 1
        end
        if (p == 2 && !mod(2^(-2*n)*d, 4) in [0,1])
            n = n -2
        end
        d0 = divexact(d, p^(2*n))
        e0 = e - 2*n
        if e0 == 0
            return p^n*mod(d0, 2)
        else return p^n * _sqrtmod4PE(d0, p^e0, e)
        end
    else
        b = _sqrtmod4P(d, p)
        f = 1
        while f < e
            k = mod(divexact(d-b^2, b*4*p^f), p)
            b = mod(b + 2*k*p^f, 2*p^(f+1))
            f = f + 1
        end
        if b > p^e
            return b-2*p^e
        else
            return b
        end
    end
end

@doc Markdown.doc"""
     prime_power_form(d::fmpz, p::fmpz, e)
Returns an integral binary quadratic form of discriminant $d$ and leading coefficient
$p$^$e$ where $p$ is a prime number.
"""
function prime_power_form(d::fmpz, p::fmpz, e)
    if !isdiscriminant(d)
        @error("$d is no discriminant")
    end
    b = _sqrtmod4PE(d, p, e)
    c = divexact(b^2-d, 4*p^e)
    f = binary_quadratic_form(p^e, b, c)
    if discriminant(f) != d
        @error("prime power form does not exist")
    end
    return f
end

################################################################################
#
#  Equivalence
#
################################################################################

@doc Markdown.doc"""
    isequivalent(f::QuadBin{fmpz}, g::QuadBin{fmpz}; proper::Bool = false)

Return whether `f` and `g` are (properly) equivalent.
"""
function isequivalent(f::QuadBin{fmpz}, g::QuadBin{fmpz}; proper::Bool = true)
  d = discriminant(f)
  if d != discriminant(g)
    return false
  end

  if issquare(d)
    return _isequivalent_reducible(f, g, proper = proper)[1]
  end

  if isindefinite(f)
    fred = reduction(f)
    gred = reduction(g)
    if issquare(d)
      # Make sure we terminate in a form with c = 0
      while !iszero(fred[3])
        fred, = _rho(fred)
      end
      while !iszero(gred[3])
        gred, = _rho(gred)
      end
      b = fred[2]
      a = fred[1]
      a0 = gred[1]
      if proper
        return (a - a0) % (2*b) == 0
      else
        g = gcd(a, b)
        return (a * a0 - g^2 ) % (2 * b * g) == 0
      end
    end

    prop_cyc = cycle(gred, proper = true)

    is_prop = fred in prop_cyc
    if proper || is_prop
      return is_prop
    end

    # note that our definition of improper equivalence
    # differs from that of Buchmann and Vollmer
    # their action is det f * q(f(x,y))
    # ours is q(f(x,y))

    # an improper equivalence in our convention

    fred = binary_quadratic_form(fred[3], fred[2], fred[1])
    @assert isreduced(fred)
    return fred in prop_cyc
  else
    if ispositive_definite(f) && !ispositive_definite(g)
      return false
    end

    if isnegative_definite(f) && !isnegative_definite(g)
      return false
    end

    fred = reduction(f)
    gred = reduction(g)
    if fred == gred
      return true
    end

    if !proper
      f1 = reduction(binary_quadratic_form(f[3], f[2], f[1]))
      return f1 == gred
    end

    return false
  end
end

function _isequivalent_reducible(f::QuadBin{fmpz}, g::QuadBin{fmpz}; proper = true)
  if discriminant(f) != discriminant(g)
    return false
  end

  c = content(f)
  if content(g) != c
    return false
  end

  fpr = divexact(f, c)
  gpr = divexact(g, c)

  fred, Tf = _reduction(fpr)
  gred, Tg = _reduction(gpr)

  fl = fred == gred

  if proper || fl
    T = Tf * inv(Tg)
    if fl
      @assert Hecke._action(f, T) == g
    end
    return fl, T
  end

  if fred[1] == invmod(gred[1], gred[2])
    gg = binary_quadratic_form(gred[1], -gred[2], zero(fmpz))
    _, Tgg = reduction_with_transformation(gg)
    T = Tf * inv(Tg * matrix(FlintZZ, 2, 2, [1, 0, 0, -1]) * Tgg)
    @assert Hecke._action(f, T) == g
    return true, T
  end

  return false, Tf
end

################################################################################
#
#  Reduction and cycles
#
################################################################################

@doc Markdown.doc"""
    reduction(f::QuadBin{fmpz}) -> QuadBin{fmpz}

Return a reduced binary quadratic form equivalent to `f`.
"""
function reduction(f::QuadBin{fmpz})
  g, _ = _reduction(f)
  return g
end

@doc Markdown.doc"""
    reduction_with_transformation(f::QuadBin{fmpz}) -> QuadBin{fmpz}, Mat{fmpz}

Return a reduced binary quadratic form `g` equivalent to `f` and a matrix `T`
such that `f.T = g`.
"""
function reduction_with_transformation(f::QuadBin{fmpz})
  return _reduction(f)
end

function _reduction(f::QuadBin{fmpz})
  if isreducible(f)
    return _reduction_reducible(f)
  end

  if isreduced(f)
    return f, identity_matrix(FlintZZ, 2)
  end

  if isindefinite(f)
    return _reduction_indefinite(f)
  end

  throw(NotImplemented())
end

# TODO (TH):
# - Make the functions operate on (a, b, c)
# - Don't build up T, just do the operations directly on U
function _reduction_indefinite(f)
  # Compute a reduced form in the proper equivalence class of f
  RR = ArbField(53, cached = false)
  U = identity_matrix(FlintZZ, 2)
  d = sqrt(RR(discriminant(f)))
  while !isreduced(f)
    a = f[1]
    b = f[2]
    c = f[3]
    cabs = abs(c)
    # Now compute rho(f) as defined on p. 122, equation (6.12) in [BV2007]
    if !iszero(cabs)
      if cabs >= d
        s = sign(c) * round(fmpz, fmpq(cabs + b, 2 * cabs), RoundDown) # floor(cabs + b/2 * abs)
      else
        @assert d > cabs # might fail with precision too low
        e = floor(divexact(d + b, 2 * cabs))
        fl, o = unique_integer(e)
        @assert fl # might fail with precision too low
        s = sign(c) * o
      end
      T = matrix(FlintZZ, 2, 2, [0, -1, 1, s])
      U = U * T
      _f = binary_quadratic_form(c, -b + 2*s*c, c*s*s - b*s + a)
      @assert _buchmann_vollmer_action(f, T) == _f
      f = _f
    else
      if b < 0
        T = matrix(FlintZZ, 2, 2, [1, 0, 0, -1])
        U = U * T
        _f = binar_quadratic_form(a, -b, c)
        @assert _buchmann_vollmer_action(f, T) == _f
        f = _f
      else
        q, r = divrem(a, b)
        if 2*r > b
          q, r = divrem(a, -b)
          q = -q
        end
        T = matrix(FlintZZ, 2, 2, [1, 0, -q, 1])
        U = U * T
        _f = binary_quadratic_form(r, b, c)
        @assert _buchmann_vollmer_action(f, T) == _f
        f = _f
      end
    end
  end
  return f, U
end

function _reduction_reducible(f::QuadBin)
  d = discriminant(f)
  N = sqrt(d)
  @assert N^2 == d
  @assert isprimitive(f)
  if iszero(f[1])
    x = -f[3]
    y = f[2]
  else
    x = N - f[2]
    y = 2 * f[1]
  end
  @assert iszero(f(x, y))
  gg = gcd(x, y)
  x = divexact(x, gg)
  y = divexact(y, gg)
  _,w, _z = gcdx(x, y)
  z = -_z
  @assert x * w - y * z == 1
  T = matrix(FlintZZ, 2, 2, [x, z, y, w])
  g = Hecke._action(f, T)
  # Now g = [0, +/- N, g[2]]
  @assert iszero(g[1])
  @assert abs(g[2]) == N
  TT = matrix(FlintZZ, 2, 2, [0, -1, 1, 0])
  g = Hecke._action(g, TT)
  T = T * TT
  # Now g = [g[1], N, 0]
  @assert abs(g[2]) == N
  if g[2] < 0
    aa = invmod(g[1], N)
    t = divexact(a * aa' - 1)
    # a * aa - N * t == 1
    @assert a * aa - N * t == 1
    TT = matrix(FlintZZ, 2, 2, [aa, -N, -t, a])
    g = Hecke._action(g, TT)
    T = T * TT
  end
  @assert g[2] == N
  _t, r = divrem(g[1], N)
  if r < 0
    r += N
    _t -= 1
  end
  @assert r >= 0
  @assert r < N
  @assert g[1] - _t * N == r
  TT = matrix(FlintZZ, 2, 2, [1, 0, -_t, 1])
  g = Hecke._action(g, TT)
  T = T * TT
  @assert 1 <= g[1] < N && g[2] == N && iszero(g[3])
  @assert det(T) == 1
  @assert g == Hecke._action(f, T)
  return g, T
end

function _buchmann_vollmer_action(f::QuadBin, M)
  a = f[1]
  b = f[2]
  c = f[3]
  s = M[1, 1]
  t = M[1, 2]
  u = M[2, 1]
  v = M[2, 2]
  a1 = f(s, u)
  b1 = 2*(a*s*t + c*u*v) + b*(s*v + t*u)
  c1 = f(t, v)
  return det(M) * binary_quadratic_form(a1, b1, c1)
end

function _action(f::QuadBin, M)
  a = f[1]
  b = f[2]
  c = f[3]
  s = M[1, 1]
  t = M[1, 2]
  u = M[2, 1]
  v = M[2, 2]
  a1 = f(s, u)
  b1 = 2*(a*s*t + c*u*v) + b*(s*v + t*u)
  c1 = f(t, v)
  return binary_quadratic_form(a1, b1, c1)
end

@doc Markdown.doc"""
    isreduced(f::QuadBin{fmpz}) -> Bool

Return whether `f` is reduced in the following sense. Let `f = [a, b, c]`
be of discriminant `D`.

If `f` is positive definite (`D < 0` and `a > 0`), then `f` is reduced if and
only if `|b| <= a <= c`, and `b >= 0` if `a = b` or `a = c`.

If `f` is negative definite (`D < 0` and `a < 0`), then `f` is reduced if and
only if `[-a, b, -c]` is reduced.

If `f` is indefinite (`D > 0), then `f` is reduced if and only if
`|sqrt{D} - 2|a|| < b < \sqrt{D}|` or `a = 0` and `-b < 2c <= b` or `c = 0` and
`-b < 2a <= b`.
"""
function isreduced(f::QuadBin{fmpz})
  D = discriminant(f)
  a = f[1]
  b = f[2]
  c = f[3]
  if D < 0 && a > 0
    return ((-a < b <= a < c) || (0 <= b <= a == c))
  elseif D <0 && a < 0
    return ((a < b <= -a < -c) || (0 <= b <= -a == -c))
  else
    # First the two easy conditions
    if (0 == a && -b < 2*c <= b) || (0 == c && -b < 2*a <= b)
      return true
    end

    if (b^2 > D)
      return false
    end

    R = ArbField(64, cached = false)
    d = sqrt(R(D))
    z = abs(d - 2 * abs(a))
    @assert !contains(z, b)
    if z < b
      return true
    else
      return false
    end
  end
end

@doc Markdown.doc"""
    cycle(f::QuadBin{fmpz}; proper::Bool = false) -> Vector{QuadBin{fmpz}}

Return the cycle of `f` as defined by Buchmann--Vollmer (Algorithm 6.1). The
cycle consists of all reduced, equivalent forms `g`, such that first coefficient of
`f` and `g` have the same sign. The proper cycle consists of all equivalent forms,
and has either the same or twice the size of the cycle. In the latter case, the
cycle has odd length.
"""
function cycle(f::QuadBin{fmpz}; proper::Bool = false)
  @req isindefinite(f) "Quadratic form must be indefinite"
  @req isreduced(f) "Quadratic form must be reduced"

  if issquare(discriminant(f))
    throw(NotImplemented())
  end

  if proper
    # Proposition 6.10.5 in [BV2007]
    # If we decided to cache this, this must be changed
    C = cycle(f, proper = false)
    if isodd(length(C))
      append!(C, C)
    end

    for i in 1:div(length(C), 2)
      C[2*i], = _tau(C[2*i]) # tau returns also the operator
    end
    return C
  end

  return _nonproper_cycle(f)
end

function _nonproper_cycle(f::QuadBin{fmpz})
  if isdefined(f, :nonproper_cycle)
    return f.nonproper_cycle::Vector{QuadBin{fmpz}}
  end
  C = typeof(f)[f]
  Q1, T = _rhotau(f)
  while !(f == Q1)
    push!(C, Q1)
    Q1, = _rhotau(Q1)
  end
  f.nonproper_cycle = C
  return C
end

# Transform f into rho(tau(f)), as defined in equation (6.12) of
# Buchmann--Vollmer 2007.
function _rhotau(f::QuadBin{fmpz})
  RR = ArbField(64, cached = false)
  d = sqrt(RR(discriminant(f)))
  a = f[1]
  b = f[2]
  c = f[3]
  cabs = abs(c)
  if cabs >= d
    s = sign(c) * round(fmpz, fmpq(cabs + b, 2 * cabs), RoundDown) # floor(cabs + b/2 * abs)
  else
    @assert d > cabs # might fail with precision too low
    e = floor(divexact(d + b, 2 * cabs))
    fl, o = unique_integer(e)
    @assert fl # might fail with precision too low
    s = sign(c) * o
  end
  g = binary_quadratic_form(-c, -b + 2*s*c, -(a - b*s + c*s*s))
  T = matrix(FlintZZ, 2, 2, [0, 1, 1, -s])
  @assert _buchmann_vollmer_action(f, T) == g
  return (g, T)
end

# Apply the rho operator as defined by Buchmann--Vollmer
function _rho(f::QuadBin{fmpz})
  RR = ArbField(64, cached = false)
  d = sqrt(RR(discriminant(f)))
  a = f[1]
  b = f[2]
  c = f[3]
  cabs = abs(c)
  # Now compute rho(f) as defined on p. 122, equation (6.12) in [BV2007]
  if cabs >= d
    s = sign(c) * round(fmpz, fmpq(cabs + b, 2 * cabs), RoundDown) # floor(cabs + b/2 * abs)
  else
    @assert d > cabs # might fail with precision too low
    e = floor(divexact(d + b, 2 * cabs))
    fl, o = unique_integer(e)
    @assert fl # might fail with precision too low
    s = sign(c) * o
  end
  T = matrix(FlintZZ, 2, 2, [0, -1, 1, s])
  g = binary_quadratic_form(c, -b + 2*s*c, a - b*s + c*s*s)
  @assert _buchmann_vollmer_action(f, T) == g
  return g, T
end

# Apply the tau operator of Buchmann--Vollmer, which turns
# [a, b, c] into [-a, b, -c]
function _tau(f::QuadBin{fmpz})
  T = matrix(FlintZZ, 2, 2, [1, 0, 0, -1])
  g = binary_quadratic_form(-f[1], f[2], -f[3])
  @assert _buchmann_vollmer_action(f, T) == g
  return g, T
end

################################################################################
#
#  Representatives
#
################################################################################

function binary_quadratic_form_representatives(d::fmpz; proper = true, primitive = false)
  d4 = mod(d, 4)
  if d4 == 2 || d4 == 3
    error("Not a discriminant")
  end
  if d > 0
    # indefinite
    return _equivalence_classes_binary_quadratic_indefinite(d, proper = proper,
                                                            primitive = primitive)
  else
    throw(NotImplemented())
  end
end

################################################################################
#
#  Genus
#
################################################################################

function islocally_equivalent(f::QuadBin{fmpz}, g::QuadBin{fmpz})
  K, = rationals_as_number_field()
  L = _binary_quadratic_form_to_lattice(f, K)
  M = _binary_quadratic_form_to_lattice(g, K)
  return genus(L) == genus(M)
end


text:
push!(LOAD_PATH,"../src")

using Test
using TypesParticles
using TypesBasis
using BuildBasis
using BuildOverlap
using InputParser
using KineticEnergy

"""
		run_buildkineticenergy_test()

This is a function to test the user types in BuildOverlap.jl. In asserts overlap integral calculationfor two hydrogen atoms is the same as the output from Prof. James John's matlab code.

""" function run_buildkineticenergy_test()
    @testset "Kinetic Energy" begin
        system,basisset = getcalcsetup("../examples/HatreeFock.in")
        basisfunc,numelec = buildbasisfunc(system,basisset);

        Scalc = buildelecoverlap(numelec,basisfunc,basisset);
        kecalc = buildkineticenergy(numelec,basisfunc,basisset);
        ketest = [1.39568e0 0.25974e0 0.75367e0 0.22767e0;
                  0.25974e0 0.24192e0 0.22767e0 0.22327e0;
                  0.75367e0 0.22767e0 1.39568e0 0.25974e0;
                  0.22767e0 0.22327e0 0.25974e0 0.24192e0];

        tolerance = 1.0e-5;
        for i=1:4
            for j=1:4
                @test sign(kecalc[i,j]) == sign(ketest[i,j]);
                diff = abs(kecalc[i,j]-ketest[i,j]);
                @test diff <= tolerance;
	    end
        end #i

    end #testset
end #run_kineticenergy_test()


text:
### A Pluto.jl notebook ###
# v0.14.2

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 2265874e-5bac-47ed-8ca9-8ce7ee98f906
begin
	import Pkg
	Pkg.activate()
	Pkg.add(["Distributions", "Plots", "PlutoUI"])
end

# ╔═╡ b03c59d4-a28c-11eb-2c3c-052e0d3a259a
using Distributions, Plots, PlutoUI, Random

# ╔═╡ f6f8b6f1-f0f5-4b1a-8a28-5b8dbf377487
md"""
Zählexperiment mit $N$ Samples $x$ mit Erwartungswert $\lambda$:

$x_i \sim \operatorname{Poisson}(\lambda)$

Ziel: Bestimme Unsicherheit auf $\lambda$, gegeben durch Verteilung $P(\lambda | x)$

Satz von Bayes:

$P(\lambda | x) = \frac{P(x | \lambda) P(\lambda)}{\int \mathrm{d}\lambda' \, P(x | \lambda') P(\lambda')}$

Starte mit $\Gamma(k_0, \theta_0)$ als Verteilung für Prior $P(\lambda)$

Likelihood $P(x_i | \lambda)$ gegeben durch $\operatorname{Poisson}(\lambda)$

$P(x | \lambda) = P\Big(\sum_i x_i \Big| n \lambda \Big)$

Posterior nach $n$ Samples $x$ $P(\lambda | x)$ hat schöne analytische Lösung als:

$\Gamma\bigg(k_0 + \sum_i x_i, \frac{\theta_0}{n \theta_0 + 1}\bigg)$
"""

# ╔═╡ c074a4b9-1e09-44c3-b7c2-0a9edeadd8cf
md"""
λ $(@bind λ Slider(0.0:100.0; show_value=true))

n $(@bind n Slider(0:200; show_value=true))

k₀ $(@bind k_0 Slider(0.01:0.01:50; show_value=true))

θ₀ $(@bind θ_0 Slider(0.01:0.01:50; default=1, show_value=true))
"""

# ╔═╡ 424c7e22-2a8d-4049-b273-beac78f5769c
begin
	poisson = Poisson(λ)
	Random.seed!(1234)
	x = rand(poisson, n)
	prior = Gamma(k_0, θ_0)
	posterior = Gamma(k_0 + sum(x), θ_0 / (n * θ_0 + 1))
end;

# ╔═╡ 1d076946-4250-4125-afff-2b517a16fe47
plot([λ -> pdf(prior, λ) λ -> pdf(posterior, λ)], 0, 100;
	layout=2, size=(700, 300), label=["Prior" "Posterior"],
	color=[1 2], lw=1.5)

# ╔═╡ bd038557-b8da-4e4e-999f-71737b6d9cff
md"""
Im Limes für flachen Prior $k_0 \rightarrow 0$, $\theta_0 \rightarrow \infty$, wobei $\theta_0 \gg \frac{1}{k_0}$:

$P(\lambda | x) \approx \operatorname{pdf}\bigg(\Gamma\bigg(\sum_i x_i, \frac{1}{n}\bigg), \lambda\bigg) \propto \operatorname{pdf}\bigg(\operatorname{Poisson}(n\lambda), \sum_i x_i\bigg) = P(x | \lambda)$

Für Likelihood-Fits mache die Annahme, dass der Posterior (Verteilung für $\lambda$ nach Messung von $x$) der Likelihood mal Normalisierungskonstante entspricht.
"""

# ╔═╡ e6743cba-2f3c-48fb-bcc2-b319eaa446ef
plot([λ -> pdf(posterior, λ) λ -> pdf(Poisson(n*λ), sum(x))], 0, 100;
	layout=2, size=(700, 300), label=["Posterior" "Poisson Likelihood"],
	color=[1 2], lw=1.5)

# ╔═╡ Cell order:
# ╠═2265874e-5bac-47ed-8ca9-8ce7ee98f906
# ╠═b03c59d4-a28c-11eb-2c3c-052e0d3a259a
# ╟─f6f8b6f1-f0f5-4b1a-8a28-5b8dbf377487
# ╟─c074a4b9-1e09-44c3-b7c2-0a9edeadd8cf
# ╠═424c7e22-2a8d-4049-b273-beac78f5769c
# ╠═1d076946-4250-4125-afff-2b517a16fe47
# ╟─bd038557-b8da-4e4e-999f-71737b6d9cff
# ╠═e6743cba-2f3c-48fb-bcc2-b319eaa446ef

max_stars_repo_name: simeonschaub/PoissonLikelihood

max_issues_repo_name: simeonschaub/PoissonLikelihood

max_forks_repo_name: simeonschaub/PoissonLikelihood


text:
function calcbranchprops!(system::LegSystem,old=Dict("a"=>0),restart="no")
    for i in 1:length(system.branches.ID)
        # Moens-Korteweg wave speed, m/s
        push!(system.branches.c0,[sqrt(0.5*system.branches.beta[i][end]/
            system.solverparams.rho)*system.branches.A0[i]^0.25])
    end
    if restart == "no"
        push!(system.branches.termscalings,1.)
        push!(system.branches.termscalings,1.)
        push!(system.branches.termscalings,1.)
        push!(system.branches.termscalings,1.)
        push!(system.branches.termscalings,1.)
        push!(system.branches.termscalings,1.)
    elseif restart == "yes"
        append!(system.branches.termscalings,old["termscalings"])
    end
end

max_stars_repo_name: dcanuto/leg-hemodynamics

max_issues_repo_name: dcanuto/leg-hemodynamics

max_forks_repo_name: dcanuto/leg-hemodynamics


text:
module TriangularMatrices

using   LinearAlgebra,
        Base.Cartesian,
        MacroTools,
        StaticArrays # Drop, eventually?

const Sized{N,T} = Union{SVector{N,T}, NTuple{N,T}}

import LinearAlgebra: Adjoint

export  SymmetricMatrix,
        StaticSymmetricMatrix,
        LowerTriangularMatrix,
        StaticLowerTriangularMatrix,
        UpperTriangularMatrix,
        StaticUpperTriangularMatrix,
        RecursiveMatrix,
        StaticRecursiveMatrix,
        # cholesky,
        cholesky!,
        reverse_cholesky,
        reverse_cholesky!,
        xxt,
        xtx,
        # inv!,
        randmat, # Should fix that rand interface!!!
        srandmat,
        randvec,
        RecursiveVector,
        choldet!,
        invdet!

# cutoff must be even
# The indexing recursion terminates at the cutoff
const halfcutoff = 4
const cutoff = 2halfcutoff
const cutoff2 = abs2(cutoff)
# const halfcut = cld(cutoff, 2)

small_triangle(x)::Int = (x-1)*x ÷ 2
big_triangle(x)::Int = (x+1)*x ÷ 2
inv_triangle(x)::Int = (Int(sqrt(1+8x))-1) ÷ 2
@generated ValSmallT(::Val{x}) where x = Val{small_triangle(x)}()
@generated ValBigT(::Val{x}) where x = Val{big_triangle(x)}()
@generated ValInvT(::Val{x}) where x = Val{inv_triangle(x)}()

include("meta.jl")
include("recursive_indexing.jl")
include("matrix_types/triangular.jl")
include("matrix_types/symmetric.jl")
include("matrix_types/recursive_matrix.jl")
include("matrix_types/rvector.jl")
include("decomp_and_inversions/cholesky.jl")
include("decomp_and_inversions/inverse.jl")
include("arithmetic/arithmetic.jl")
include("arithmetic/gemm.jl")
include("arithmetic/symv.jl")

end # module


text:
"""
Methods for initializing the RBM and sampling from its distribution.
"""

#=
    f_logcoshsum(theta)
Inline function for calculating the log of '∏_i cosh(θ_i(NetState))'
=#
@inbounds @views @fastmath f_logcoshsum(theta) = log(prod(cosh.(theta)))

"""
    update_theta!(NetState, NetParams, flips, NetSettings)
In-place update of theta_upd for proposed spin flips.
"""
function update_theta!(NetState, NetParams, flips, NetSettings)
    @unpack m = NetSettings
    @unpack theta_upd, state = NetState
    @unpack w = NetParams

     for j=1:m
        for i in flips
             theta_upd[j] -= 2 * w[i,j] * state[i]
        end
    end
    NetState.logcoshsum_upd = f_logcoshsum(NetState.theta_upd)

end

"""
    psi_mh_revert(NetState, NetParams, flips, NetSettings)

Returns the ψ'/ψ ratio, given flips and network parameters. Does not modify
the network state. Used for energy calculation: see `find_energy`.
"""
function psi_mh_revert(NetState, NetParams, flips, NetSettings)
    @unpack theta, theta_upd, state = NetState
    @unpack a = NetParams
    update_theta!(NetState, NetParams, flips, NetSettings)
    logpsi_frac = 0.0im
    for i in flips
        logpsi_frac -= 2 * a[i] * state[i]
    end
    logpsi_frac += NetState.logcoshsum_upd - NetState.logcoshsum
    theta_upd .= theta
    return @fastmath exp(logpsi_frac)
end

"""
    psi_mh(NetState, NetParams, flips, NetSettings)

Returns the ψ'/ψ ratio, given the network state, the network parameters and the
flips. See also `flipper`  and `flipper_exchange`.
"""
function psi_mh(NetState, NetParams, flips, NetSettings)
    @unpack a = NetParams
    @unpack state = NetState

    update_theta!(NetState, NetParams, flips, NetSettings)
    logpsi_frac = 0.0im
    for i in flips
        logpsi_frac -= 2 * a[i] * state[i]
    end
    logpsi_frac += NetState.logcoshsum_upd - NetState.logcoshsum
    return @fastmath exp(logpsi_frac)
end

"""
    flipper_exchange(NetState, NetParams, NetSettings)

Performs a Metropolis-Hastings update. See also `psi_mh` and `flipper`.
Preserves magnetization.
"""
function flipper_exchange(NetState, NetParams, NetSettings)
    @unpack n = NetSettings
    @unpack state, theta, theta_upd = NetState

    flips = zeros(Int, 2)
    flips[1] = rand(1:n)
    flips[2] = flips[1] % n + 1  # ASSUMES PBC, exchange with Hamming dist 1

    while state[flips[1]] == state[flips[2]]
        flips[1] = flips[1] % n + 1
        flips[2] = flips[2] % n + 1
    end
    if abs2(psi_mh(NetState, NetParams, flips, NetSettings)) > rand()
        theta .= theta_upd
        NetState.logcoshsum = NetState.logcoshsum_upd
        state[flips] .= reverse(state[flips])
    else
        theta_upd .= theta
    end
end

"""
    flipper(NetState, NetParams, NetSettings)

Performs a Metropolis-Hastings update. See also `psi_mh`.
"""
function flipper(NetState, NetParams, NetSettings)
    @unpack nflips, n = NetSettings
    @unpack theta, state, theta_upd = NetState

    if nflips == 1
        flips  = rand(1:n)
    elseif nflips == 2
        if rand(Bool)
            flips = rand(1:n)
        else
            flips = sample(1:n, 2, replace=false)
        end
    end
    if @fastmath abs2(psi_mh(NetState, NetParams, flips, NetSettings)) > rand()
        theta .= theta_upd
        NetState.logcoshsum = NetState.logcoshsum_upd
        state[flips] *= -1
    else
        theta_upd .= theta
    end
end

"""
    find_energy(NetState, NetParams, HamiltWeights, NetSettings)

Saves the local energy in-place given Hamiltonian weight. See also
`mc_sampling()`.
"""
function find_energy(NetState, NetParams, HamiltWeights, NetSettings)
    @unpack w_x, w_y, w_z, w_xx, w_yy, w_zz = HamiltWeights
    @unpack state = NetState
    @unpack n = NetSettings

    eloc = 0.0im
    factor_tmp = 0.0im
    for i=1:n
        factor_tmp = w_x[i] + 1.0im * w_y[i] * state[i]
        if !iszero(factor_tmp)  # don't perform costly psi'/psi computation when multiplied by 0 anyway
            eloc += factor_tmp * psi_mh_revert(NetState, NetParams, i, NetSettings)
        end
        eloc += w_z[i] * state[i]
        for j=i+1:n
            state_ij_tmp = state[i] * state[j]
            factor_tmp = w_xx[i, j] - w_yy[i, j] * state_ij_tmp
            if !iszero(factor_tmp)
                eloc += factor_tmp * psi_mh_revert(NetState, NetParams, [i, j], NetSettings)
            end
            eloc += w_zz[i, j] * state_ij_tmp
        end
    end
    return eloc
end

"""
    derivs!(o_tot, NetState, NetSettings, trial, id)

In-place saving of variational derivatives `adjoint(O)` in `o_tot` for trial number `trial`.
"""
function derivs!(o_tot, NetState, NetSettings, trial, id)
    @unpack n, m = NetSettings
    @unpack state, theta = NetState

    o_tot[1:n, trial, id] .= state
    o_tot[n + 1:n + m, trial, id] .= conj.(tanh.(theta))
    @views o_tot[n + m + 1:end, trial, id] .= kron(o_tot[n + 1:n + m, trial, id], state)
end

"""
    reinit_state!(NetState, NetParams, NetSettings)

In-place reinitialization of spin state.
"""
function reinit_state!(NetState, NetParams, NetSettings)
    @unpack mag0, n = NetSettings
    @unpack state = NetState

    if mag0
        for i=1:Int(n/2)
            state[i] = -1
            state[end - i + 1] = +1
        end
        shuffle!(state)
    else
        state .= 2 * rand(0:1, n) .- 1
    end
    reset_theta!(NetState, NetParams)
end

"""
    reset_theta!(NetState, NetParams)

In-place recalculation look-up variables.
"""
function reset_theta!(NetState, NetParams)
    @unpack b, w = NetParams
    @unpack theta, theta_upd, state = NetState

    theta .= b + transpose(w) * state
    theta_upd .= theta
    NetState.logcoshsum = f_logcoshsum(theta)
end

"""
    initial_thermalize!(NetState, NetParams, NetSettings)

Initial thermalization for when the network parameters are just initialized.
"""
function initial_thermalize(NetState, NetParams, NetSettings)
    @unpack init_therm_steps, mc_steps, mag0 = NetSettings
    for i=1:init_therm_steps
        for j=1:mc_steps
            if mag0
                flipper_exchange(NetState, NetParams, NetSettings)
            else
                flipper(NetState, NetParams, NetSettings)
            end
        end
    end
end

"""
    thermalize(NetState, NetParams, NetSettings)

Thermalization after network parameter update.
"""
function thermalize(NetState, NetParams, NetSettings)
    @unpack therm_steps, mc_steps, mag0 = NetSettings
    for i=1:therm_steps
        for j=1:mc_steps
            if mag0
                flipper_exchange(NetState, NetParams, NetSettings)
            else
                flipper(NetState, NetParams, NetSettings)
            end
        end
    end
end

"""
    mc_sampling!(NetParams, o_tot, elocs, HamiltWeights, NetSettings[, warm_state])

# Arguments

- `NetParams`: the network parameters W={a,b,w}.
- `VarDervs`: a VARDERIVS-type with variational derivatives (updated in-place).
- `elocs`: a SharedArray-type, of which column `id` is updated in-place.
- `HamiltWeights`: the weights of the 1/2-local Hamiltonian. See `weights.jl`.
- `NetSettings`: The settings for optimization.
# Keywords

- `warm_state`: The `NetState` struct is dereferenced when the function is
called. However, the state could be passed to the master process without much
overhead to give the next Markov Chain a "warm start" when the network
parameters are updated.
"""
function mc_sampling!(NetParams, o_tot, elocs, HamiltWeights, NetSettings, id; warm_state=zeros(Int, 1))
    @unpack n, mc_trials, mc_steps, mag0 = NetSettings
    # Initialization
    NetState = init_state(NETSTATE, NetSettings, warm_state=warm_state)
    reset_theta!(NetState, NetParams)

    # Thermalization
    if warm_state[1] == 0  # first iter
        initial_thermalize(NetState, NetParams, NetSettings)
    else
        thermalize(NetState, NetParams, NetSettings)
    end

    # Sampling
    for t=1:mc_trials
        for i=1:mc_steps
            if mag0
                flipper_exchange(NetState, NetParams, NetSettings)
            else
                flipper(NetState, NetParams, NetSettings)
            end
        end
        derivs!(o_tot, NetState, NetSettings, t, id)
        elocs[t, id] = find_energy(NetState, NetParams, HamiltWeights, NetSettings)
    end
    return NetState.state
end

"""
    init_params(ParamStruct, NetSettings)
Constructs a `NETPARAMS` type with random network parameters. Note that these
parameters should be accessible from all processors and are thus SharedArrays.
"""
function init_params(ParamStruct, NetSettings)
    @unpack n, m, dim_rbm = NetSettings
    μ_params = 0.0
    σ_params = 0.01 / dim_rbm
    a = complex(1,0) * rand(Normal(μ_params, σ_params), n) +
        complex(0,1) * rand(Normal(μ_params, σ_params), n)
    b = complex(1,0) * rand(Normal(μ_params, σ_params), m) +
        complex(0,1) * rand(Normal(μ_params, σ_params), m)
    w = complex(1,0) * rand(Normal(μ_params, σ_params), n, m) +
        complex(0,1) * rand(Normal(μ_params, σ_params), n, m)
    netparameters = ParamStruct(SharedArray(a), SharedArray(b), SharedArray(w))
    return netparameters
end

"""
    init_state(StateStruct, NetSettings[, warm_state])
Constructs a `NETSTATE` type with a random/warm network state and the
corresponding look-up tables (θ and its cosh-product).
"""
function init_state(StateStruct, NetSettings; warm_state=zeros(Int, 1))
    @unpack n, mag0, m = NetSettings
    if warm_state[1] != 0
        state = warm_state
    elseif mag0 && warm_state[1] == 0
        state = zeros(Int, n)
        for i=1:Int(n/2)
            state[i] = -1
            state[end - i + 1] = +1
        end
        shuffle!(state)
    elseif !mag0 && warm_state[1] == 0
        state = 2 * rand(0:1, n) .- 1
    end
    theta = zeros(Complex{Float64}, m)
    theta_upd = copy(theta)
    logcoshsum = 0.0im
    logcoshsum_upd = 0.0im
    netstate = StateStruct(state, theta, theta_upd, logcoshsum, logcoshsum_upd)
end

"""
    init_vardervs(VarDervStruct, NetSettings)

Takes a Variational Derivative type as argument. Consequently initializes the
variational derivative matrix and the regularization parameter. NOTE that
adjoint(O) is saved during sampling, instead of O. This is done so that samples
can be stored in C-order in `derivs!` in `network.jl`. Generates a variational-
derivative matrix for each worker.
"""
function init_vardervs(VarDervStruct, NetSettings)       # Loads an Opt constructor and initiates
    @unpack regulator, dim_rbm, mc_trials = NetSettings
    o_tot = SharedArray(zeros(Complex{Float64}, dim_rbm, mc_trials, nworkers()))
    reg = regulator
    varderv_inst = VarDervStruct(o_tot, reg)
    return varderv_inst
end


text:
### Example set up for seismic or medical experiments
### Author: Ziyi Yin* and Rafael Orozco*
### March, 2020

function GenMuteSeismic()
	##Muting to properly visualize
	mute_region = ones(Float32,n[1],n[2]);
	mute_region[:,1:60].=0;
	mute_mask = imfilter(mute_region, Kernel.gaussian(10));
	#q_mask_smooth = imfilter(reshape(q_stack,model.n), Kernel.gaussian(30))

    return mute_mask
end

function GenMuteMedical()
	##Muting to properly visualize

	#pading for receivers
	pad_rec = 10 + 10

	#muting due to recerivers
	mute_region = ones(Float32,n[1],n[2]);
	mute_region[:,1:100].=0;
	mute_region[:,600:end].=0;
	mute_region[1:100,:].=0;
	mute_region[1100:end,:].=0;
	mute_mask = imfilter(mute_region, Kernel.gaussian(10));

	q_stack=sum(q,dims=4)[:,:,1,1];
	#muting due to sources
	q_smooth = imfilter(q_stack[:,:,1,1], Kernel.gaussian(10))
	q_smooth[q_smooth .> 0.00001] .= 1

	q_smooth_2 = imfilter(q_smooth, Kernel.gaussian(10))
	q_smooth_2_flip = 1 .- q_smooth_2
	mute_mask_full = mute_mask .* q_smooth_2_flip

    return mute_mask_full
end

function GenMuteMedical_brain()
	##Muting to properly visualize

	#muting due to recerivers

	pad = 40
	mute_region = ones(Float32,n[1],n[2]);
	mute_region[:,1:pad].=0;
	mute_region[:,n[2]-pad:end].=0;
	mute_region[1:pad,:].=0;
	mute_region[n[1]-pad:end,:].=0;
	mute_mask = imfilter(mute_region, Kernel.gaussian(10));

	q_stack=sum(q,dims=4)[:,:,1,1];
	#muting due to sources
	q_smooth = imfilter(q_stack[:,:,1,1], Kernel.gaussian(20))
	q_smooth[q_smooth .> 0.00001] .= 1

	q_smooth_2 = imfilter(q_smooth, Kernel.gaussian(10))
	q_smooth_2_flip = 1 .- q_smooth_2
	mute_mask_full = mute_mask .* q_smooth_2_flip

    return mute_mask_full
end

function GenMuteMedical_breast()
	##Muting to properly visualize

	#muting due to recerivers
	rad     = 250
	centerx = 300
	centery = 300
	mute_region = ones(Float32,n[1],n[2]);

	for i in 1:n[1]
		for j in 1:n[2]
			if (i-centerx)^2 + (j-centery)^2 > rad^2
				mute_region[i,j] = 0
			end
		end
	end

	mute_mask = imfilter(mute_region, Kernel.gaussian(10));

	q_stack=sum(q,dims=4)[:,:,1,1];
	#muting due to sources
	q_smooth = imfilter(q_stack[:,:,1,1], Kernel.gaussian(15))
	q_smooth[q_smooth .> 0.0001] .= 1

	q_smooth_2 = imfilter(q_smooth, Kernel.gaussian(5))
	q_smooth_2_flip = 1 .- q_smooth_2
	mute_mask_full = mute_mask .* q_smooth_2_flip

    return mute_mask_full
end

function GenRecGeometrySeismic()

	# Set up receiver geometry
	nrec = 500
	domain_x = (model.n[1] - 1)*model.d[1]
	xmin = 50f0*d[1]
	xmax = domain_x - xmin
	xrec = range(xmin, stop=xmax, length=nrec)
	yrec = 0f0
	zrec = range(30f0*d[2], 30f0*d[2], length=nrec)
	nsrc = 90
	# Set up receiver structure

    return xrec,yrec,zrec,nrec,nsrc
end

function GenRecGeometryMedical()

	nsrc = 1
	pad = 50f0

	domain_x = (model0.n[1] - 1)*model0.d[1]
	domain_z = (model0.n[2] - 1)*model0.d[2]

	xmin = pad*d[1]
	xmax = domain_x - xmin

	zmin = pad*d[2]
	zmax = domain_z - zmin

	#for sampling distance:
	#we put 500 receivers on top and bottom then calculate the sampling distance
	#for sides.
	nrec_top_bottom = 500
	sampling_dist = (domain_x*model0.d[1]/(nrec_top_bottom-1))
	nrec_side = Int(round(domain_z*model0.d[1] / sampling_dist))

	#lateral x distance
	xrec_top    = range(xmin, stop=xmax, length=nrec_top_bottom)
	xrec_bottom = range(xmin, stop=xmax, length=nrec_top_bottom)
	xrec_left   = range(xmin, stop=xmin, length=nrec_side)
	xrec_right  = range(xmax, stop=xmax, length=nrec_side)

	#depth z distance
	zrec_top    = range(zmin,   zmin, length=nrec_top_bottom)
	zrec_bottom = range(zmax,   zmax, length=nrec_top_bottom)
	zrec_left   = range(zmin,   zmax, length=nrec_side)
	zrec_right  = range(zmin,   zmax, length=nrec_side)

	yrec = 0f0 #2d so always 0

	xrec = vcat(xrec_top,xrec_bottom,xrec_left,xrec_right)
	zrec = vcat(zrec_top,zrec_bottom,zrec_left,zrec_right)

	#xrec = vcat(xrec_top)
	#zrec = vcat(zrec_top)

	nrec = size(xrec)[1]    # total no. of receivers

    return xrec, yrec, zrec, nrec, nsrc
end

function GenSrcGeometrySeismic(nsrc)

	# Multiple source locations (locations as indices)
	if nsrc == 1
		xsrc_index = Int.(round.((range(800,stop=800,length=nsrc))))   # indices where sim source weights are 1
		zsrc_index = Int.(round.((range(30,stop=30,length=nsrc))))
		return xsrc_index,zsrc_index
	end
	xsrc_index = Int.(round.((range(50,stop=1150,length=nsrc))))   # indices where sim source weights are 1
	zsrc_index = Int.(round.((range(30,stop=30,length=nsrc))))

    return xsrc_index,zsrc_index
end

function GenSrcGeometryMedical()
	nsrc_simult = 1

	# Multiple source locations (locations as indices)
	#xsrc_index = Int.(round.((range(50,stop=1150,length=nsrc_simult))))   # indices where sim source weights are 1
	#xsrc_index = Int.(round.((range(300,stop=1150,length=nsrc_simult))))   # indices where sim source weights are 1

	xsrc_index = Int.(round.((range(800,stop=800,length=nsrc_simult))))
	zsrc_index = Int.(round.((range(550,stop=550,length=nsrc_simult))))
	# load data and associated grid Array

	# file         = matopen(input_path*"NIC_400_full_Phantom_360.mat");
	# p0           = read(file, "p0");

	# size_rows = 501
	# size_cols = 501
	# p0_resize = imresize(p0, (size_rows, size_cols))
	# p0_resize[p0_resize .> 0] .= 1.0

	# init_x = 100
	# init_z = 200
	# in_model = zeros(Float32,n[1],n[2]);
	# in_model[init_x:init_x+size_rows-1,init_z:init_z+size_cols-1] .= p0_resize
	# #imshow(in_model')

	# # load source locations
	# inds         = findall(x->x>0.0,in_model);

	# xsrc_index = convert(Array{Int64,1},zeros(size(inds)));
	# zsrc_index = convert(Array{Int64,1},zeros(size(inds)));

	# for i=1:length(inds)
	#        xsrc_index[i] = inds[i][1];
	#        zsrc_index[i] = inds[i][2];
	# end

	# nsrc_simult = size(xsrc_index)[1]

    return xsrc_index, zsrc_index, nsrc_simult
end

function Gen_Exp_Medical_Skull_Phantom_Source()
	########----Setup receivers
	nsrc = 1

	domain_x = (model0.n[1] - 1)*model0.d[1]
	domain_z = (model0.n[2] - 1)*model0.d[2]

	file         = matopen("NIC_441_full_Phantom_60.mat");
	p0           = read(file, "p0");

	size_rows = 100
	size_cols = 100
	p0_resize = imresize(p0, (size_rows, size_cols))
	p0_resize[p0_resize .> 0] .= 1.0

	init_x = 80
	init_z = 100
	in_model = zeros(Float32,n[1],n[2]);
	in_model[init_x:init_x+size_rows-1,init_z:init_z+size_cols-1] .= p0_resize
	#imshow(in_model')

	# load source locations
	inds         = findall(x->x>0.0,in_model);

	xsrc_index = convert(Array{Int64,1},zeros(size(inds)));
	zsrc_index = convert(Array{Int64,1},zeros(size(inds)));

	for i=1:length(inds)
	       xsrc_index[i] = inds[i][1];
	       zsrc_index[i] = inds[i][2];
	end

	nsrc_simult = size(xsrc_index)[1]

	# load data and associated grid Array

	# file         = matopen(input_path*"NIC_400_full_Phantom_360.mat");
	# p0           = read(file, "p0");

	# size_rows = 501
	# size_cols = 501
	# p0_resize = imresize(p0, (size_rows, size_cols))
	# p0_resize[p0_resize .> 0] .= 1.0

	# init_x = 100
	# init_z = 200
	# in_model = zeros(Float32,n[1],n[2]);
	# in_model[init_x:init_x+size_rows-1,init_z:init_z+size_cols-1] .= p0_resize
	# #imshow(in_model')

	# # load source locations
	# inds         = findall(x->x>0.0,in_model);

	# xsrc_index = convert(Array{Int64,1},zeros(size(inds)));
	# zsrc_index = convert(Array{Int64,1},zeros(size(inds)));

	# for i=1:length(inds)
	#        xsrc_index[i] = inds[i][1];
	#        zsrc_index[i] = inds[i][2];
	# end

	# nsrc_simult = size(xsrc_index)[1]

	nrec = 512

	function circleShape(h, k ,r)
		theta = LinRange(0,2*pi, nrec)
		h .+ r*sin.(theta), k .+ r*cos.(theta)
	end

	xrec, zrec = circleShape(domain_x / 2,domain_z/2, domain_z/2-5*d[2])
	yrec = 0f0 #2d so always 0

	s=rand(1:512,100);
	xrec_rand = xrec[s];
	zrec_rand = zrec[s];
	nrec = 100

    return xrec_rand, yrec, zrec_rand, nrec, nsrc, xsrc_index, zsrc_index, nsrc_simult
#	return xrec, yrec, zrec, nrec, nsrc, xsrc_index, zsrc_index, nsrc_simult
end

function Gen_Exp_Medical_Skull_Single_Source()
	########----Setup receivers
	nsrc = 1
	nsrc_simult=1

	domain_x = (model0.n[1] - 1)*model0.d[1]
	domain_z = (model0.n[2] - 1)*model0.d[2]

	xsrc_index = Int.(round.((range(128,stop=128,length=nsrc_simult))))
	zsrc_index = Int.(round.((range(128,stop=128,length=nsrc_simult))))
	# load data and associated grid Array

	# file         = matopen(input_path*"NIC_400_full_Phantom_360.mat");
	# p0           = read(file, "p0");

	# size_rows = 501
	# size_cols = 501
	# p0_resize = imresize(p0, (size_rows, size_cols))
	# p0_resize[p0_resize .> 0] .= 1.0

	# init_x = 100
	# init_z = 200
	# in_model = zeros(Float32,n[1],n[2]);
	# in_model[init_x:init_x+size_rows-1,init_z:init_z+size_cols-1] .= p0_resize
	# #imshow(in_model')

	# # load source locations
	# inds         = findall(x->x>0.0,in_model);

	# xsrc_index = convert(Array{Int64,1},zeros(size(inds)));
	# zsrc_index = convert(Array{Int64,1},zeros(size(inds)));

	# for i=1:length(inds)
	#        xsrc_index[i] = inds[i][1];
	#        zsrc_index[i] = inds[i][2];
	# end

	# nsrc_simult = size(xsrc_index)[1]

	nrec = 512

	function circleShape(h, k ,r)
		theta = LinRange(0,2*pi, nrec)
		h .+ r*sin.(theta), k .+ r*cos.(theta)
	end

	xrec, zrec = circleShape(domain_x / 2,domain_z/2, domain_z/2-5*d[2])
	yrec = 0f0 #2d so always 0

    return xrec, yrec, zrec, nrec, nsrc, xsrc_index, zsrc_index, nsrc_simult
end

function Gen_Exp_Medical_Calc_Single_Source()

	########----Setup receivers
	nsrc = 1
	pad = 10f0

	domain_x = (model0.n[1] - 1)*model0.d[1]
	domain_z = (model0.n[2] - 1)*model0.d[2]

	xmin = pad*d[1]
	xmax = domain_x - xmin

	zmin = pad*d[2]
	zmax = domain_z - zmin

	#for sampling distance:
	#we put 500 receivers on top and bottom then calculate the sampling distance
	#for sides.
	nrec_top_bottom = 50
	sampling_dist = (domain_x*model0.d[1]/(nrec_top_bottom-1))
	nrec_side = Int(round(domain_z*model0.d[1] / sampling_dist))

	#lateral x distance
	xrec_top    = range(xmin, stop=xmax, length=nrec_top_bottom)
	xrec_bottom = range(xmin, stop=xmax, length=nrec_top_bottom)
	xrec_left   = range(xmin, stop=xmin, length=nrec_side)
	xrec_right  = range(xmax, stop=xmax, length=nrec_side)

	#depth z distance
	zrec_top    = range(zmin,   zmin, length=nrec_top_bottom)
	zrec_bottom = range(zmax,   zmax, length=nrec_top_bottom)
	zrec_left   = range(zmin,   zmax, length=nrec_side)
	zrec_right  = range(zmin,   zmax, length=nrec_side)

	yrec = 0f0 #2d so always 0

	xrec = vcat(xrec_top,xrec_bottom,xrec_left,xrec_right)
	zrec = vcat(zrec_top,zrec_bottom,zrec_left,zrec_right)

	#xrec = vcat(xrec_top)
	#zrec = vcat(zrec_top)

	nrec = size(xrec)[1]    # total no. of receivers

	########----Setup source

	#this experiment puts a single under a large wall of high velocity (bone)

	nsrc_simult = 1

	# Multiple source locations (locations as indices)
	#xsrc_index = Int.(round.((range(50,stop=1150,length=nsrc_simult))))   # indices where sim source weights are 1
	#xsrc_index = Int.(round.((range(300,stop=1150,length=nsrc_simult))))   # indices where sim source weights are 1

	xsrc_index = Int.(round.((range(300,stop=300,length=nsrc_simult))))
	zsrc_index = Int.(round.((range(300,stop=300,length=nsrc_simult))))

	# load data and associated grid Array

	# file         = matopen(input_path*"NIC_400_full_Phantom_360.mat");
	# p0           = read(file, "p0");

	# size_rows = 501
	# size_cols = 501
	# p0_resize = imresize(p0, (size_rows, size_cols))
	# p0_resize[p0_resize .> 0] .= 1.0

	# init_x = 100
	# init_z = 200
	# in_model = zeros(Float32,n[1],n[2]);
	# in_model[init_x:init_x+size_rows-1,init_z:init_z+size_cols-1] .= p0_resize
	# #imshow(in_model')

	# # load source locations
	# inds         = findall(x->x>0.0,in_model);

	# xsrc_index = convert(Array{Int64,1},zeros(size(inds)));
	# zsrc_index = convert(Array{Int64,1},zeros(size(inds)));

	# for i=1:length(inds)
	#        xsrc_index[i] = inds[i][1];
	#        zsrc_index[i] = inds[i][2];
	# end

	# nsrc_simult = size(xsrc_index)[1]

    return xrec, yrec, zrec, nrec, nsrc, xsrc_index, zsrc_index, nsrc_simult
end

function Gen_Exp_Medical_Calc_Breast_Single()
	########----Setup source
	nsrc = 1
	nsrc_simult = 1

	# Artificial source location (locations as indices)
	xsrc_index = Int.(round.((range(300,stop=300,length=nsrc_simult))))
	zsrc_index = Int.(round.((range(300,stop=300,length=nsrc_simult))))

	######----Artificial source location (locations as indices)
	#blood_inds =findall(v[x->x == 1584.0])
	#xsrc_index = blood_inds[1,:]
	#zsrc_index = blood_inds[2,:]

	######------SETUP RING TRANSDUCER ARRAY
	domain_x = (model0.n[1] - 1)*model0.d[1]
	domain_z = (model0.n[2] - 1)*model0.d[2]

	nrec = 512
	function circleShape(h, k ,r)
		theta = LinRange(0,2*pi, nrec)
		h .+ r*sin.(theta), k .+ r*cos.(theta)
	end

	xrec, zrec = circleShape(domain_x / 2,domain_z/2, domain_z/2-1600)
	yrec = 0f0 #2d so always 0

    return xrec, yrec, zrec, nrec, nsrc, xsrc_index, zsrc_index, nsrc_simult
end

function Gen_Exp_Medical_Calc_Breast_Multiple()
	########----Setup source
	nsrc = 1
	nsrc_simult = 2

	# Artificial source location (locations as indices)
	xsrc_index = Int.(round.((range(200,stop=400,length=nsrc_simult))))
	zsrc_index = Int.(round.((range(300,stop=300,length=nsrc_simult))))

	#xsrc_index = vcat(xsrc_index, Int.(round.((range(300,stop=300,length=nsrc_simult)))))
	#zsrc_index = vcat(zsrc_index, Int.(round.((range(200,stop=400,length=nsrc_simult)))))

	######----Artificial source location (locations as indices)
	#blood_inds =findall(v[x->x == 1584.0])
	#xsrc_index = blood_inds[1,:]
	#zsrc_index = blood_inds[2,:]

	######------SETUP RING TRANSDUCER ARRAY
	domain_x = (model0.n[1] - 1)*model0.d[1]
	domain_z = (model0.n[2] - 1)*model0.d[2]

	nrec = 512
	function circleShape(h, k ,r)
		theta = LinRange(0,2*pi, nrec)
		h .+ r*sin.(theta), k .+ r*cos.(theta)
	end

	xrec, zrec = circleShape(domain_x / 2,domain_z/2, domain_z/2-1600)
	yrec = 0f0 #2d so always 0

    return xrec, yrec, zrec, nrec, nsrc, xsrc_index, zsrc_index, nsrc_simult
end

function Gen_Exp_Seismic_Calc_Breast_Multiple()
	########----Setup source
	nsrc = 4
	nsrc_simult = nsrc

	# Artificial source location (locations as indices)
	#xsrc_index = []
	#zsrc_index = []

	for i in 1:nsrc
		xsrc_index = vcat(xsrc_index, rand(200:400))
		zsrc_index = vcat(zsrc_index, rand(200:400))
	end
	#xsrc_index = Int.(round.((range(200,stop=400,length=nsrc_simult))))
	#zsrc_index = Int.(round.((range(300,stop=300,length=nsrc_simult))))

	#xsrc_index = vcat(xsrc_index, Int.(round.((range(300,stop=300,length=nsrc_simult)))))
	#zsrc_index = vcat(zsrc_index, Int.(round.((range(200,stop=400,length=nsrc_simult)))))

	######----Artificial source location (locations as indices)
	#blood_inds =findall(v[x->x == 1584.0])
	#xsrc_index = blood_inds[1,:]
	#zsrc_index = blood_inds[2,:]

	######------SETUP RING TRANSDUCER ARRAY
	domain_x = (model0.n[1] - 1)*model0.d[1]
	domain_z = (model0.n[2] - 1)*model0.d[2]

	nrec = 512
	function circleShape(h, k ,r)
		theta = LinRange(0,2*pi, nrec)
		h .+ r*sin.(theta), k .+ r*cos.(theta)
	end

	xrec, zrec = circleShape(domain_x / 2,domain_z/2, domain_z/2-1600)
	yrec = 0f0 #2d so always 0

    return xrec, yrec, zrec, nrec, nsrc, xsrc_index, zsrc_index, nsrc_simult
end

function Gen_Exp_Medical_Calc_Breast_Blood()
	########----Setup source
	nsrc = 1
	nsrc_simult = 1

	# Artificial source location (locations as indices)
	#xsrc_index = Int.(round.((range(300,stop=300,length=nsrc_simult))))
	#zsrc_index = Int.(round.((range(300,stop=300,length=nsrc_simult))))

	######----Artificial source location (locations as indices)
	#blood_inds =findall(x->x == 1584.0,v)

	#xsrc_index = getindex.(blood_inds, 1)
	#zsrc_index = getindex.(blood_inds, 2)

	####MANUALLY ADD SINGLE SOURCES
	zsrc_index = [428,495,539, 309, 342, 300, 108]
	xsrc_index = [392,380,213, 403, 445, 494, 413]
	nsrc_simult = 7

	######------SETUP RING TRANSDUCER ARRAY
	domain_x = (model0.n[1] - 1)*model0.d[1]
	domain_z = (model0.n[2] - 1)*model0.d[2]

	nrec = 512
	function circleShape(h, k ,r)
		theta = LinRange(0,2*pi, nrec)
		h .+ r*sin.(theta), k .+ r*cos.(theta)
	end

	xrec, zrec = circleShape(domain_x / 2,domain_z/2, domain_z/2-1600)

	yrec = 0f0 #2d so always 0

    return xrec, yrec, zrec, nrec, nsrc, xsrc_index, zsrc_index, nsrc_simult
end

function Gen_Exp_Seismic_Calc_Breast_Blood()
	########----Setup source
	nsrc = 7
	nsrc_simult = 7

	# Artificial source location (locations as indices)
	#xsrc_index = Int.(round.((range(300,stop=300,length=nsrc_simult))))
	#zsrc_index = Int.(round.((range(300,stop=300,length=nsrc_simult))))

	######----Artificial source location (locations as indices)
	# blood_inds =findall(x->x == 1584.0,v)

	# xsrc_index = getindex.(blood_inds, 1)
	# zsrc_index = getindex.(blood_inds, 2)

	####MANUALLY ADD SINGLE SOURCES
	#zsrc_index = [428]
	#xsrc_index = [392]
	zsrc_index = [428,495,539, 309, 342, 300, 108]
	xsrc_index = [392,380,213, 403, 445, 494, 413]
	nsrc_simult = 7

	######------SETUP RING TRANSDUCER ARRAY
	domain_x = (model0.n[1] - 1)*model.d[1]
	domain_z = (model0.n[2] - 1)*model.d[2]

	nrec = 512
	function circleShape(h, k ,r)
		theta = LinRange(0,2*pi, nrec)
		h .+ r*sin.(theta), k .+ r*cos.(theta)
	end

	xrec, zrec = circleShape(domain_x / 2,domain_z/2, domain_z/2)

	yrec = 0f0 #2d so always 0

    return xrec, yrec, zrec, nrec, nsrc, xsrc_index, zsrc_index, nsrc_simult
end


text:
using Documenter, IssueReporter

makedocs()