text:
#1D Matrix operator utility functions.
spI(n) = spdiagm(n,n, 0 => ones(n))
∂1d(n,dx) = spdiagm(n,n+1,0 => -ones(n), 1 => ones(n))/dx
c(n) = spdiagm(n,n+1,0 => ones(n), 1 => ones(n))/2
χ(n) = spdiagm(n,n+2, 1 => ones(n))

"""
    opvec(model::AbstractModel,vec::AbstractVector)

Function to multiply a vector by the momentum operator.
"""
function opvec(model::AbstractModel,vec::AbstractVector)
    @unpack gh,gu,gv,gc=model.fields
    @assert length(vec)==(gu.n+gv.n)
    uspread=gu.spread*vec[1:gu.n]
    vspread=gv.spread*vec[(gu.n+1):(gu.n+gv.n)]
    extra = gh.dimplicit[]*(gu.∂x*(gu.crop*(gu.h[:].*uspread)).+gv.∂y*(gv.crop*(gv.h[:].*vspread)))
    opvecprod=
    [
     #x-component
     gu.samp*(gu.∂x'*(2gh.dneghηav[]*(2gu.∂x*uspread .+ gv.∂y*vspread)) .+
              gu.∂y'*(gc.crop*(gc.cent'*(gh.dneghηav[]*(gc.cent*(gc.crop*( gu.∂y*uspread .+ gv.∂x*vspread )))))) .+
              gu.dnegβeff[]*uspread .+ (gu.h[:].*(gu.∂x'*(extra))) )
        ;
     #y-component
     gv.samp*(gv.∂y'*(2gh.dneghηav[]*(2gv.∂y*vspread .+ gu.∂x*uspread)) .+
              gv.∂x'*(gc.crop*(gc.cent'*(gh.dneghηav[]*(gc.cent*(gc.crop*( gv.∂x*vspread .+ gu.∂y*uspread )))))) .+
              gv.dnegβeff[]*vspread .+ (gv.h[:].*(gv.∂y'*(extra))) )
    ]
    return opvecprod
end
"""
    restrictvec(model::AbstractModel,vec::AbstractVector)

Function to restrict a vector from the fine grid to the coarse grid, used in multigrid preconditioner.
"""
function restrictvec(model::AbstractModel,vec::AbstractVector)
    @unpack wu,wv,gu,gv=model.fields
    @assert length(vec)==(gu.n+gv.n)
    vecx=vec[1:gu.n]
    vecy=vec[(gu.n+1):(gu.n+gv.n)]
    restrictvec=
    [
     #x-component
     wu.samp[]*(wu.idwt'*(gu.spread*vecx))
        ;
     #y-component
     wv.samp[]*(wv.idwt'*(gv.spread*vecy))
    ]
    return restrictvec
end
"""
    prolongvec(model::AbstractModel,waveletvec::AbstractVector)

Function to prolong a vector from the coarse grid to the fine grid, used in multigrid preconditioner.
"""
function prolongvec(model::AbstractModel,waveletvec::AbstractVector)
    @unpack wu,wv,gu,gv=model.fields
    @assert length(waveletvec)==(wu.n[]+wv.n[])
    waveletvecx = waveletvec[1:wu.n[]]
    waveletvecy = waveletvec[(wu.n[]+1):(wu.n[]+wv.n[])]
    prolongvec =
    [
     #x-component
     gu.samp*(wu.idwt*(wu.spread[]*waveletvecx))
        ;
     #y-component
     gv.samp*(wv.idwt*(wv.spread[]*waveletvecy))
    ]
    return prolongvec
end

"""
pos_fraction(z1;mask=mask) -> area_fraction, area_fraction_u, area_fraction_v

Return fraction of each grid cell with function z1 above zero. Uses bilinear
interpolation of values at three nearest cell centers to represent the function.
In:
   z1:             m x n array of gridded function values.
   mask:           m x n mask   1 = valid data, 0= invalid data.
Out:
   area_fraction:   m x n array showing area fraction of interpolated z1>0 on h-grid.
   area_fraction_u: m+1 x n array showing area fraction of interpolated z1>0 on u-grid.
   area_fraction_v: m x n+1 array showing area fraction of interpolated z1>0 on v-grid.
"""
function pos_fraction(z1::AbstractArray{T,2};mask=trues(size(z1))) where {T};

    m,n=size(z1);

    area_fraction=zeros(T,m,n);
    area_fraction_u=zeros(T,m+1,n);
    area_fraction_v=zeros(T,m,n+1);

    for quadrant=1:4

        area_fraction_quadrant=zeros(T,m,n);

        #Define (x_i,y_i,z_i) triples i=1,2,3 for each point
        #and its two nearest neighbours in each of two orthogonal directions.
        #Distance unit is half a grid cell so each grid cell
        #is made up of four unit-square quadrants.
        #Quadrant 1: defined by neighbours in grid East & grid North direction
        #In local coordinate system.
        #(x_1,y_1,z_1)=(0,0,z_1);
        #(x_2,y_2,z_2)=(2,0,z_2); Grid East
        #(x_3,y_3,z_3)=(0,2,z_3); Grid North
        #Triples for quadrants 2-4 are defined by neighbours in the other
        #directions.
        #Quadrant 2: Neighbours from grid North and grid West
        #Quadrant 3: Neighbours from grid West and grid South
        #Quadrant 4: Neighbours from grid South and grid East
        #If there is no valid neighbour in any direction the central value is
        #reused.
        if quadrant == 1
              z2=z1[[2:m;m],1:n];
              z2[.!mask[[2:m;m],1:n]]=z1[.!mask[[2:m;m],1:n]];
              z3=z1[1:m,[2:n;n]];
              z3[.!mask[1:m,[2:n;n]]]=z1[.!mask[1:m,[2:n;n]]];
        elseif quadrant == 2
              z2=z1[1:m,[2:n;n]];
              z2[.!mask[1:m,[2:n;n]]]=z1[.!mask[1:m,[2:n;n]]];
              z3=z1[[1;1:(m-1)],1:n];
              z3[.!mask[[1;1:(m-1)],1:n]]=z1[.!mask[[1;1:(m-1)],1:n]];
        elseif quadrant == 3
              z2=z1[[1;1:(m-1)],1:n];
              z2[.!mask[[1;1:(m-1)],1:n]]=z1[.!mask[[1;1:(m-1)],1:n]];
              z3=z1[1:m,[1;1:(n-1)]];
              z3[.!mask[1:m,[1;1:(n-1)]]]=z1[.!mask[1:m,[1;1:(n-1)]]];
        elseif quadrant == 4
              z2=z1[1:m,[1;1:(n-1)]];
              z2[.!mask[1:m,[1;1:(n-1)]]]=z1[.!mask[1:m,[1;1:(n-1)]]];
              z3=z1[[2:m;m],1:n];
              z3[.!mask[[2:m;m],1:n]]=z1[.!mask[[2:m;m],1:n]];
        else
                error("Quadrant not defined");
        end

        #Find triples where sign changes.
        sc=findall(.!((sign.(z1) .== sign.(z2)) .& (sign.(z1) .== sign.(z3))));
        nsc=length(sc);

        #Initialise array.
        areasc=zeros(T,nsc);

        #Equation (y=ax+b) of zero contour for planar elements;
        #Planar elements are defined by (x_i,y_i,z_i) triples i=1,2,3.
        #N.B. distance unit for x and y is half a grid cell.
        #Flip x and y if needed by exchanging z2 and z3, to avoid infinite a;
        flip=(abs.(z3[sc] .- z1[sc]) .< abs.(z2[sc] .- z1[sc]));
        a=zeros(T,nsc);
        b=zeros(T,nsc);
        a[.!flip].=-(z2[sc[.!flip]].-z1[sc[.!flip]])./(z3[sc[.!flip]].-z1[sc[.!flip]]);
        b[.!flip].=-2*z1[sc[.!flip]]./(z3[sc[.!flip]].-z1[sc[.!flip]]);
        a[flip].=-(z3[sc[flip]].-z1[sc[flip]])./(z2[sc[flip]].-z1[sc[flip]]);
        b[flip].=-2*z1[sc[flip]]./(z2[sc[flip]].-z1[sc[flip]]);

        #Areas of various useful triangles and trapezoids.
        a1::Vector{T}=0.5*(b.^2)./a;
        a2::Vector{T}=0.5*((1.0 .- b).^2)./a;
        a3::Vector{T}=0.5*a.+b;

        #Test which edges of unit square are intercepted by zero contour.
        test1::Vector{Int64} = (b .> 0.0) .+ (b .> 1.0);
        test2::Vector{Int64} = ((a+b) .> 0.0) .+ ((a+b) .> 1.0);

        #There are nine possible combinations ix=1:9.
        ix= 1 .+ test1 .+ 3*test2;

        #Work out area fraction of unit square for each different combination.
        areasc[ix.==1].= 0.0;
        areasc[ix.==2].=-a1[ix.==2];
        areasc[ix.==3].=a2[ix.==3].-a1[ix.==3];
        areasc[ix.==4].=a1[ix.==4].+a3[ix.==4];
        areasc[ix.==5].=a3[ix.==5];
        areasc[ix.==6].=a2[ix.==6].+a3[ix.==6];
        areasc[ix.==7].=1.0.-a2[ix.==7].+a1[ix.==7];
        areasc[ix.==8].=1.0.-a2[ix.==8];
        areasc[ix.==9].=1.0;

        #Checks to see whether this is the area above zero or below zero.
        areasc[b.<0.0].=1.0.-areasc[b.<0.0];
        areasc[z1[sc].<0.0].=1.0.-areasc[z1[sc].<0.0];

        #Compile areas for sign changes with trivial cases.
        area_fraction_quadrant[sc].=areasc;
        area_fraction_quadrant[(z1.>0.0) .& (z2.>0.0) .& (z3.>0.0)].=1.0;
        area_fraction_quadrant[(z1.<=0.0) .& (z2.<=0.0) .& (z3.<=0.0)].=0.0;

        #Add areas for each quadrant together.
        #N.B. distance unit for x and y is half a grid cell.
        #Summing over quadrants gives the fraction of each grid cell above zero.
        area_fraction.=area_fraction.+0.25*area_fraction_quadrant;

        if quadrant == 1
            area_fraction_u.=area_fraction_u.+0.25*[zeros(T,1,n);area_fraction_quadrant];
        elseif quadrant == 2
            area_fraction_u.=area_fraction_u.+0.25*[area_fraction_quadrant;zeros(T,1,n)];
        elseif quadrant == 3
            area_fraction_u.=area_fraction_u.+0.25*[area_fraction_quadrant;zeros(T,1,n)];
        elseif quadrant == 4
            area_fraction_u.=area_fraction_u.+0.25*[zeros(T,1,n);area_fraction_quadrant];
        else
            error("Quadrant not defined");
        end

        if quadrant == 1
            area_fraction_v.=area_fraction_v.+0.25*[zeros(T,m,1) area_fraction_quadrant];
        elseif quadrant == 2
            area_fraction_v.=area_fraction_v.+0.25*[zeros(T,m,1) area_fraction_quadrant];
        elseif quadrant == 3
            area_fraction_v.=area_fraction_v.+0.25*[area_fraction_quadrant zeros(T,m,1)];
        elseif quadrant == 4
            area_fraction_v.=area_fraction_v.+0.25*[area_fraction_quadrant zeros(T,m,1)];
        else
            error("Quadrant not defined");
        end

    end

    return area_fraction, area_fraction_u, area_fraction_v

end

#MISMIP+ bed elevation
function mismip_plus_bed(x,y)
    xbar = 300000.0
    b0 = -150.0; b2 = -728.8; b4 = 343.91; b6 = -50.75
    wc = 24000.0; fc = 4000.0; dc = 500.0
    bx(x)=b0+b2*(x/xbar)^2+b4*(x/xbar)^4+b6*(x/xbar)^6
    by(y)= dc*( (1+exp(-2(y-wc)/fc))^(-1) + (1+exp(2(y+wc)/fc))^(-1) )
    b = max(bx(x) + by(y), -720.0)
    return b
end

"""
     get_glx(model)

Return the grounding line in the form x = x(y). Assumes each y-row has at least
one grid point where ice grounded and one where ice floating.
"""
function get_glx(model)
    @unpack fields, grid = model
    glmask=diff(sign.(fields.gh.haf),dims=1).==-2 #calculate where sign of height above floating passes thru zero
    glx1=grid.xxh[1:end-1,:][glmask] #x co-ordiates upstream of grounding line
    glx2=grid.xxh[2:end,:][glmask] #x co-ordinates immediately downstream
    haf1=fields.gh.haf[1:end-1,:][glmask] #Height above floatation immediately upstream of gl
    haf2=fields.gh.haf[2:end,:][glmask]
    glx=glx1+(glx2-glx1).*(zero(haf1)-haf1)./(haf2-haf1) #interpolate between grid points to find grounding line
    return glx
end

"""
    icedraft(s,h,sea_level_wrt_geoid)

Compute ice draft.
"""
icedraft(s,h,sea_level_wrt_geoid)=max(h-(s-sea_level_wrt_geoid),zero(typeof(h-(s-sea_level_wrt_geoid))))

"""
    height_above_floatation(h,b,params)

Compute height above floatation.

"""
height_above_floatation(h,b,params) = h - (params.density_ocean/params.density_ice)*(params.sea_level_wrt_geoid - b)

"""
    volume_above_floatation(h,b,params)

Compute the volume above floatation: integrated height above floatation for cells with positive height above floatation
"""
volume_above_floatation(h,b,params,grid) = sum(sum(height_above_floatation.(h,b,params)[height_above_floatation.(h,b, params) .> 0])) .* grid.dx .* grid.dy

"""
    glen_b(temperature,damage,params)

Compute stiffness parameter B in Glen flow law.

"""

function glen_b(temperature,damage,glen_a_ref, glen_n, glen_a_activation_energy, glen_temperature_ref, gas_const)
    glen_a0 = glen_a_ref*exp(+glen_a_activation_energy/(glen_temperature_ref*gas_const) )
    glen_b = (1-damage)*( glen_a0*exp(-glen_a_activation_energy/(temperature*gas_const)) )^(-1.0/glen_n)
    return glen_b
end

"""
    get_u_mask(h_mask)

Find mask of valid grid points on u-grid corresponding to a mask defined on h-grid.

"""
function get_u_mask(h_mask)
    #include all u faces next to a selected center
    (nx,ny)=size(h_mask)
    u_mask=falses(nx+1,ny)
    u_mask[1:end-1,1:end]=u_mask[1:end-1,1:end].|h_mask
    u_mask[2:end,1:end]=u_mask[2:end,1:end].|h_mask
    return u_mask
end
"""
    get_v_mask(h_mask)

Find mask of valid grid points on v-grid corresponding to a mask defined on h-grid.

"""
function get_v_mask(h_mask)
    #include all v faces next to a selected center
    (nx,ny)=size(h_mask)
    v_mask=falses(nx,ny+1)
    v_mask[1:end,1:end-1]=v_mask[1:end,1:end-1].|h_mask
    v_mask[1:end,2:end]=v_mask[1:end,2:end].|h_mask
    return v_mask
end
"""
    get_c_mask(h_mask)

Find mask of valid grid points on c-grid corresponding to a mask defined on h-grid.

"""
function get_c_mask(h_mask)
    #select cell corners with four neighbouring cell centers in h_mask
    c_mask=h_mask[1:end-1,1:end-1] .& h_mask[1:end-1,2:end] .& h_mask[2:end,1:end-1] .& h_mask[2:end,2:end]
    return c_mask
end
"""
    clip(trial_mask)

Find mask of valid grid points on h-grid corresponding to a trial mask, also defined on h-grid.
Clip any grid points from the trial mask that cannot be used in the model.

"""
function clip(trial_mask)

    good_corners=get_c_mask(trial_mask)

    #include all centers next to a selected corner
    (nx,ny)=size(trial_mask)
    mask=falses(nx,ny)
    mask[1:end-1,1:end-1]=mask[1:end-1,1:end-1].|good_corners
    mask[1:end-1,2:end]=mask[1:end-1,2:end].|good_corners
    mask[2:end,1:end-1]=mask[2:end,1:end-1].|good_corners
    mask[2:end,2:end]=mask[2:end,2:end].|good_corners

    return mask
end


text:
"""
    Bisect_Abs(S::BnBSolver,B::BnBModel,N::Vector{Interval{Float64}})

Returns two interval boxes 'X1,X2' created by bisecting 'N' in the highest width dimension.
"""

function bisect(X::MCInterval)
    m = mid(X)
    return (MCInterval(X.lo, m), MCInterval(m, X.hi))
end

function Bisect_Abs(S::BnBSolver,B::BnBModel{T},N::Vector{T}) where {T}
  i::Int64 = indmax(diam.(N))
  N1::T,N2::T = bisect(N[i])
  X1::Vector{T} = deepcopy(N)
  X2::Vector{T} = deepcopy(N)
  X1[i] = N1
  X2[i] = N2
  return X1,X2
end

"""
    Bisect_Rel(S::BnBSolver,B::BnBModel,N::Vector{Interval{Float64}})

Returns two interval boxes 'X1,X2' created by bisecting 'N' in the highest width dimension after scaling by initial box size.
"""
function Bisect_Rel(S::BnBSolver,B::BnBModel{T},N::Vector{T}) where {T}
  i::Int64 = indmax(diam.(N)./diam.(B.Init_Box))
  N1::T,N2::T = bisect(N[i])
  X1::Vector{T} = deepcopy(N)
  X2::Vector{T} = deepcopy(N)
  X1[i] = N1
  X2[i] = N2
  return X1,X2
end

"""
    Bisect_Abs_Imp(S::BnBSolver,B::BnBModel,N::Vector{Interval{Float64}})

Returns two interval boxes 'X1,X2' created by bisecting 'N' in the highest width
dimension greater than 'nx'.
"""
function Bisect_Abs_Imp(S::BnBSolver,B::BnBModel{T},N::Vector{T},nx::Int64) where {T}
  i::Int64 = indmax(diam.(N[(nx+1):end]))
  N1::T,N2::T = bisect(N[nx+i])
  X1::Vector{T} = deepcopy(N)
  X2::Vector{T} = deepcopy(N)
  X1[nx+i] = N1
  X2[nx+i] = N2
  return X1,X2
end

"""
    Bisect_Rel_Imp(S::BnBSolver,B::BnBModel,N::Vector{Interval{Float64}})

Returns two interval boxes 'X1,X2' created by bisecting 'N' in the highest width
dimension greater than 'nx' after scaling by initial box size.
"""
function Bisect_Rel_Imp(S::BnBSolver,B::BnBModel{T},N::Vector{T},nx::Int64) where {T}
  i::Int64 = indmax(diam.(N[(nx+1):end])./diam.(B.Init_Box[(nx+1):end]))
  N1::T,N2::T = bisect(N[nx+i])
  X1::Vector{T} = deepcopy(N)
  X2::Vector{T} = deepcopy(N)
  X1[nx+i] = N1
  X2[nx+i] = N2
  return X1,X2
end

ext: jl


text:
include("graphs.jl")
using DataFrames
using Plots
using Statistics
using HypothesisTests
using CSV
using Weave
using IJulia

function benchmark(exp::Experiment)
    name = exp.name
    workspace = exp.workspace

    if !ispath(workspace * "/notebooks")
        mkpath(workspace * "/notebooks")
    end

    convert_doc(
        "./src/benchmark/notebooks/overview.jmd",
        workspace * "/notebooks/overview.ipynb",
    )

    IJulia.jupyterlab(dir = workspace * "/notebooks")
end

function benchmark(path::String = "_results")
    df = DataFrames.DataFrame()
    for item in readdir(path)
        println(item)
        maindf = DataFrame(CSV.File(path * '/' * item))
        rename!(maindf, :iteration => :runtime, :size => :problem_size)
        df = vcat(df, maindf[:, [:algorithm, :runtime, :problem_size]])
    end
    resultpvalue(df)
    return benchmark(df)
end

function benchmark(traces::Array{Trace,1})
    df = tracestodataframe(traces)
    resultpvalue(df)
    return benchmark(df)
end

function tracestodataframe(traces::Array{Trace,1})
    result = DataFrames.DataFrame()
    for trace in traces
        push!(
            result,
            (
                algorithm = string(typeof(trace.algorithm)),
                runtime = trace.optimum[2],
                problem_size = length(trace.individual.problem),
            ),
        )
    end
    result
end

function jobstodataframe(jobs::Array{Job,1}) end

function benchmark(data::DataFrames.DataFrame)
    algorithms = keys(groupby(data, "algorithm"))
    plts = []

    append!(plts, [plot()])
    for algo in algorithms
        classical_0 = filter(row -> row["algorithm"] == algo[1], data)
        classical_0 = combine(groupby(classical_0, :problem_size), :runtime => mean)
        y = classical_0[!, :runtime_mean]
        x = classical_0[!, :problem_size]
        # v = classical_0[!,4]
        plot!(
            plts[end],
            x,
            y,
            fillalpha = 0.3,
            formatter = :plain,
            dpi = 300,
            fontfamily = "Times New Roman Bold",
            legend = :right,
            linewidth = 2,
            markersize = 3,
            markershape = :circle,
            linealpha = 0.5,
            size = (600, 300),
            label = algo,
        )
    end
    return plts[1]
end

function resultpvalue(df::DataFrames.DataFrame)
    algo = keys(groupby(df, :algorithm))
    problem_size = keys(groupby(df, :problem_size))

    for item in problem_size
        println("Problem Size: ", item["problem_size"])
        data2 = groupby(
            filter(row -> row["problem_size"] == item["problem_size"], df),
            :algorithm,
        )
        for item in data2
            for item2 in data2
                p_value = pvalue(MannWhitneyUTest(item["runtime"], item2["runtime"]))
                println(item[1][1], " vs ", item2[1][1], ", Pvalue = ", p_value)
            end
        end
    end
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export daq_leddisplay

"""
    daq_leddisplay(;kwargs...)

A LEDDisplay component.
A 7-bar LED display component.
Keyword arguments:
- `id` (String; optional): The ID used to identify the display in Dash callbacks
- `backgroundColor` (String; optional): Color of the display's background
- `className` (String; optional): Class to apply to the root component element
- `color` (String; optional): Color of the display
- `label` (optional): Description to be displayed alongside the control. To control styling,
pass an object with label and style properties. label has the following type: String | lists containing elements 'style', 'label'.
Those elements have the following types:
  - `style` (Dict; optional)
  - `label` (String; optional)
- `labelPosition` (a value equal to: 'top', 'bottom'; optional): Where the display label is positioned
- `size` (Real; optional): Size of the display
- `style` (Dict; optional): Style to apply to the root component element
- `theme` (Dict; optional): Theme configuration to be set by a ThemeProvider
- `value` (Real | String; optional): Value to be displayed. A number or a string
containing only digits (0-9), periods, and colons,
and possibly starting with a minus sign.
"""
function daq_leddisplay(; kwargs...)
        available_props = Symbol[:id, :backgroundColor, :className, :color, :label, :labelPosition, :size, :style, :theme, :value]
        wild_props = Symbol[]
        return Component("daq_leddisplay", "LEDDisplay", "dash_daq", available_props, wild_props; kwargs...)
end


text:
#=
   follows the algorithms as given in
   "Some Highly Accurate Basic Linear Algebra Subroutines"
   by Yuhang Zhao (c) 2010
   McMaster University, Hamilton, Canada

   except for Round3, where I have replaced
     a call to frexp with reinterpret logic

   comments are as they appear in the thesis
=#

rc = Ref{Int}(0) # indicates if a recursive call of iFastSumAlgorithm occurs

function iFastSum{T<:Real}(x::Array{T,1})
    global rc;
    rc[] = 0
    n = length(x)
    xs = Array{T,1}(n) # iFastSumAlgorithm is destructive
    copy!(xs, x)       # thanks to Kristoffer Carlsson
    iFastSumAlgorithm(xs, n)
end

@inline function AddTwo{T<:Real}(a::T, b::T)
    x = a+b
    z = x-a
    y = (a-(x-z))+(b-z)
    (x,y)
end

@inline function Round3{T<:Real}(s0::T, s1::T, s2::T)
    # s,e = frexp(s1)
    # if (s!=0.5 || (sign(s1)!=sign(s2)))
    if (((reinterpret(UInt64,s1) & 0x000fffffffffffff) != zero(UInt64)) || (sign(s1) != sign(s2)))
        s0
    else
        1.1*s1 + s0
    end
end

# output: the correctly rounded sum of x
# for initial call, n is length of x
function iFastSumAlgorithm{T<:Real}(x::Array{T,1},n::Int)
    global rc
    s = zero(T); loop = 1; # loop counts the number of loops

    for i in 1:n # accumlate first approximation
        s,x[i] = AddTwo(s, x[i])
    end

    # main loop
    while true # loop forever
        # (1)
        count=1; st=zero(T); loop=loop+1; sm=zero(T)
        # count points to the next position in x to store the local error
        # st is the temporary sum
        # (2)
        for i in 1:n
            st,x[count] = AddTwo(st,x[i]) # x[count] is local error
            if x[count] != zero(T)
                count += 1
                sm = max(sm,abs(st))
            end
        end
        # (3)
        em = (count-1)*eps(sm)*0.5
        # each local error <= HalfUlp(sm)
        # em is the weak upper bound on magnitude of the sum of the errors
        # (4)
        s, st = AddTwo(s,st)
        st = x[count]
        n = count
        # (5)
        if (em==zero(T) || em < eps(s)*0.5)
            if rc[] > 0
                return s # return s if it is a recursive call
            end
            w1, e1 = AddTwo(st, em)
            w2, e2 = AddTwo(st, -em)
            if ((w1+s != s) | (w2+s != s)) || (Round3(s,w1,e1) != s) || (Round3(s,w2,e2) != s)
                rc[]=1
                s1 = iFastSum(x,n) # first recursive call
                s,s1 = AddTwo(s,s1)
                s2 = iFastSum(x,n) # second recursive call
                rc[] = 0
                s = Round3(s,s1,s2)
             end
             return s
        end
    end
end


text:
#--------------------------------------------------------------------
# Spacetime Discretization methods in Julia
# Soham 02-2018
#--------------------------------------------------------------------

function projectonPatchBndbyRestriction(fn::Function, Nx::Int, M::Int, loc::Int)::Array{Float64,1}
    f2Npatch = Float64[fn(x) for x in chebgrid(2Nx, M, loc)]
    invndmx  = inv(vandermonde(2Nx))
    pvndmx   = pseudovandermonde(2Nx, chebgrid(Nx))
    fNpatch  = pvndmx*restrictmodes!(invndmx*f2Npatch, Nx)
    return fNpatch
end

function projectonPatchbyRestriction(fn::Function, Nx::Int, Ny::Int, M::Int, loc::Array{Int,1})::Array{Float64,2}
    f2Npatch = Float64[fn(x,y) for x in chebgrid(2Nx, M, loc[1]), y in chebgrid(2Ny, M, loc[2])]
    invndmx  = inv(vandermonde(2Nx))
    invndmy  = inv(vandermonde(2Ny))
    pvndmx   = pseudovandermonde(2Nx, chebgrid(Nx))
    pvndmy   = pseudovandermonde(2Ny, chebgrid(Ny))
    fNpatch  = pvndmx*restrictmodes!(invndmx*f2Npatch*invndmy', Nx, Ny)*pvndmy'
    return fNpatch
end


text:
const GaussChebyshevT{T} = QuadratureMeasure{T,<:GridArrays.ChebyshevTNodes,<:GridArrays.ChebyshevTWeights}
const GaussChebyshev = GaussChebyshevT

GaussChebyshevT(n) = GaussChebyshevT{Float64}(n)
GaussChebyshevT{T}(n) where {T} = QuadratureMeasure{T}(GridArrays.gausschebyshev(T,n)...)

const GaussChebyshevU{T} = QuadratureMeasure{T,<:GridArrays.ChebyshevUNodes,<:GridArrays.ChebyshevUWeights}

GaussChebyshevU(n) = GaussChebyshevU{Float64}(n)
GaussChebyshevU{T}(n) where {T} = QuadratureMeasure{T}(GridArrays.gausschebyshevu(T,n)...)

const GaussLegendre{T} = QuadratureMeasure{T,<:GridArrays.LegendreNodes,<:GridArrays.LegendreWeights}

GaussLegendre(n) = GaussLegendre{Float64}(n)
GaussLegendre{T}(n) where {T} = QuadratureMeasure{T}(GridArrays.gausslegendre(T,n)...)

const GaussLaguerre{T} = QuadratureMeasure{T,<:GridArrays.LaguerreNodes,<:GridArrays.LaguerreWeights}

GaussLaguerre(n, α::T) where {T <: AbstractFloat} = GaussLaguerre{T}(n, α)
GaussLaguerre(n, α) = GaussLaguerre{Float64}(n, α)
GaussLaguerre{T}(n, α) where {T} = QuadratureMeasure{T}(GridArrays.gausslaguerre(T, n, α)...)

laguerre_α(μ::GaussLaguerre) = μ.points.α

const GaussHermite{T} = QuadratureMeasure{T,<:GridArrays.HermiteNodes,<:GridArrays.HermiteWeights}

GaussHermite(n) = GaussHermite{Float64}(n)
GaussHermite{T}(n) where {T} = QuadratureMeasure{T}(GridArrays.gausshermite(T,n)...)

const GaussJacobi{T} = QuadratureMeasure{T,<:GridArrays.JacobiNodes,<:GridArrays.JacobiWeights}

GaussJacobi(n, α, β) = GaussJacobi(n, promote(α, β)...)
GaussJacobi(n, α::T, β::T) where {T} = GaussJacobi{float(T)}(n, α, β)
GaussJacobi{T}(n, α, β) where {T} = QuadratureMeasure{T}(GridArrays.gaussjacobi(T, n, α, β)...)

jacobi_α(μ::GaussJacobi) = μ.points.α
jacobi_β(μ::GaussJacobi) = μ.points.β


text:
using Test, MusicManipulations

let

rnotes = randomnotes(1000)
coun = Int[]
pit = pitches(rnotes)
pitchdict = Dict(pitch => count(x->x == pitch, pit) for pitch in pit)

@testset "allowed pitches" begin

    allowed = unique(pit)[1:2]
    anotes = filterpitches(rnotes, allowed)
    @test length(anotes) > 0
    @test length(anotes) == pitchdict[allowed[1]] + pitchdict[allowed[2]]
    for note ∈ anotes
        @test note.pitch ∈ allowed
    end

end
@testset "separate pitches" begin
    sep = separatepitches(rnotes)
    @test sort(collect(keys(sep))) == sort(unique(pitches(rnotes)))

    @test sum(length(n) for n in values(sep)) == length(rnotes)

end
end

@testset "estimate delay" begin
    for midi in (readMIDIFile("serenade_full.mid"), readMIDIFile(testmidi()))
        piano = getnotes(midi, 4)

        d = estimate_delay(piano, 0:1//3:1)
        @test abs(d) < 30

        d2 = estimate_delay_recursive(piano, 0:1//3:1, 5)

        @test d2 ≥ round(Int, d)
        @test abs(d2) < 30
    end
end

@testset "combine" begin
    rnotes = randomnotes(1000)
    # test with Array
    shuf = [rnotes[k:5:end] for k in 1:3]
    com = combine(shuf)
    @test issorted(com, by = x -> x.position)
    # test with Dict
    sepa = separatepitches(rnotes)
    com = combine(sepa)
    @test com.notes == rnotes.notes
    @test issorted(com, by = x -> x.position)
end


text:
module OMETIFF

using Unitful
using AxisArrays
using EzXML
using FileIO
using JSON
using ImageCore
using ImageMetadata
using DataStructures
using DocStringExtensions

include("utils.jl")
include("files.jl")
include("parsing.jl")
include("mmap.jl")
include("loader.jl")

end # module


text:
# ---
# title: 1483. Kth Ancestor of a Tree Node
# id: problem1483
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/kth-ancestor-of-a-tree-node/description/>
# hidden: true
# ---
#
# You are given a tree with `n` nodes numbered from `0` to `n-1` in the form of
# a parent array where `parent[i]` is the parent of node `i`. The root of the
# tree is node `0`.
#
# Implement the function `getKthAncestor``(int node, int k)` to return the
# `k`-th ancestor of the given `node`. If there is no such ancestor, return
# `-1`.
#
# The  _k-th  _ _ancestor_  of a tree node is the `k`-th node in the path from
# that node to the root.
#
#
#
# **Example:**
#
# **![](https://assets.leetcode.com/uploads/2019/08/28/1528_ex1.png)**
#
#
#
#     Input:
#     ["TreeAncestor","getKthAncestor","getKthAncestor","getKthAncestor"]
#     [[7,[-1,0,0,1,1,2,2]],[3,1],[5,2],[6,3]]
#
#     Output:
#     [null,1,0,-1]
#
#     Explanation:
#     TreeAncestor treeAncestor = new TreeAncestor(7, [-1, 0, 0, 1, 1, 2, 2]);
#
#     treeAncestor.getKthAncestor(3, 1);  // returns 1 which is the parent of 3
#     treeAncestor.getKthAncestor(5, 2);  // returns 0 which is the grandparent of 5
#     treeAncestor.getKthAncestor(6, 3);  // returns -1 because there is no such ancestor
#
#
#
#
# **Constraints:**
#
#   * `1 <= k <= n <= 5*10^4`
#   * `parent[0] == -1` indicating that `0` is the root node.
#   * `0 <= parent[i] < n` for all `0 < i < n`
#   * `0 <= node < n`
#   * There will be at most `5*10^4` queries.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# filesystem operations

export
    ctime,
    filemode,
    filesize,
    gperm,
    isblockdev,
    ischardev,
    isdir,
    isfifo,
    isfile,
    islink,
    ismount,
    ispath,
    issetgid,
    issetuid,
    issocket,
    issticky,
    lstat,
    mtime,
    operm,
    stat,
    uperm

struct StatStruct
    device  :: UInt
    inode   :: UInt
    mode    :: UInt
    nlink   :: Int
    uid     :: UInt
    gid     :: UInt
    rdev    :: UInt
    size    :: Int64
    blksize :: Int64
    blocks  :: Int64
    mtime   :: Float64
    ctime   :: Float64
end

StatStruct() = StatStruct(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)

StatStruct(buf::Union{Vector{UInt8},Ptr{UInt8}}) = StatStruct(
    ccall(:jl_stat_dev,     UInt32,  (Ptr{UInt8},), buf),
    ccall(:jl_stat_ino,     UInt32,  (Ptr{UInt8},), buf),
    ccall(:jl_stat_mode,    UInt32,  (Ptr{UInt8},), buf),
    ccall(:jl_stat_nlink,   UInt32,  (Ptr{UInt8},), buf),
    ccall(:jl_stat_uid,     UInt32,  (Ptr{UInt8},), buf),
    ccall(:jl_stat_gid,     UInt32,  (Ptr{UInt8},), buf),
    ccall(:jl_stat_rdev,    UInt32,  (Ptr{UInt8},), buf),
    ccall(:jl_stat_size,    UInt64,  (Ptr{UInt8},), buf),
    ccall(:jl_stat_blksize, UInt64,  (Ptr{UInt8},), buf),
    ccall(:jl_stat_blocks,  UInt64,  (Ptr{UInt8},), buf),
    ccall(:jl_stat_mtime,   Float64, (Ptr{UInt8},), buf),
    ccall(:jl_stat_ctime,   Float64, (Ptr{UInt8},), buf),
)

show(io::IO, st::StatStruct) = print(io, "StatStruct(mode=0o$(string(filemode(st), base = 8, pad = 6)), size=$(filesize(st)))")

# stat & lstat functions

macro stat_call(sym, arg1type, arg)
    return quote
        stat_buf = zeros(UInt8, ccall(:jl_sizeof_stat, Int32, ()))
        r = ccall($(Expr(:quote, sym)), Int32, ($(esc(arg1type)), Ptr{UInt8}), $(esc(arg)), stat_buf)
        if !(r in (0, Base.UV_ENOENT, Base.UV_ENOTDIR, Base.UV_EINVAL))
            throw(_UVError("stat", r))
        end
        st = StatStruct(stat_buf)
        if ispath(st) != (r == 0)
            error("stat returned zero type for a valid path")
        end
        return st
    end
end

stat(fd::OS_HANDLE)         = @stat_call jl_fstat OS_HANDLE fd
stat(path::AbstractString)  = @stat_call jl_stat  Cstring path
lstat(path::AbstractString) = @stat_call jl_lstat Cstring path
if RawFD !== OS_HANDLE
    global stat(fd::RawFD)  = stat(Libc._get_osfhandle(fd))
end
stat(fd::Integer)           = stat(RawFD(fd))

"""
    stat(file)

Returns a structure whose fields contain information about the file.
The fields of the structure are:

| Name    | Description                                                        |
|:--------|:-------------------------------------------------------------------|
| size    | The size (in bytes) of the file                                    |
| device  | ID of the device that contains the file                            |
| inode   | The inode number of the file                                       |
| mode    | The protection mode of the file                                    |
| nlink   | The number of hard links to the file                               |
| uid     | The user id of the owner of the file                               |
| gid     | The group id of the file owner                                     |
| rdev    | If this file refers to a device, the ID of the device it refers to |
| blksize | The file-system preferred block size for the file                  |
| blocks  | The number of such blocks allocated                                |
| mtime   | Unix timestamp of when the file was last modified                  |
| ctime   | Unix timestamp of when the file was created                        |

"""
stat(path...) = stat(joinpath(path...))

"""
    lstat(file)

Like [`stat`](@ref), but for symbolic links gets the info for the link
itself rather than the file it refers to.
This function must be called on a file path rather than a file object or a file
descriptor.
"""
lstat(path...) = lstat(joinpath(path...))

# some convenience functions

"""
    filemode(file)

Equivalent to `stat(file).mode`.
"""
filemode(st::StatStruct) = st.mode

"""
    filesize(path...)

Equivalent to `stat(file).size`.
"""
filesize(st::StatStruct) = st.size

"""
    mtime(file)

Equivalent to `stat(file).mtime`.
"""
mtime(st::StatStruct) = st.mtime

"""
    ctime(file)

Equivalent to `stat(file).ctime`.
"""
ctime(st::StatStruct) = st.ctime

# mode type predicates

"""
    ispath(path) -> Bool

Return `true` if a valid filesystem entity exists at `path`,
otherwise returns `false`.
This is the generalization of [`isfile`](@ref), [`isdir`](@ref) etc.
"""
ispath(st::StatStruct) = filemode(st) & 0xf000 != 0x0000

"""
    isfifo(path) -> Bool

Return `true` if `path` is a FIFO, `false` otherwise.
"""
isfifo(st::StatStruct) = filemode(st) & 0xf000 == 0x1000

"""
    ischardev(path) -> Bool

Return `true` if `path` is a character device, `false` otherwise.
"""
ischardev(st::StatStruct) = filemode(st) & 0xf000 == 0x2000

"""
    isdir(path) -> Bool

Return `true` if `path` is a directory, `false` otherwise.

# Examples
```jldoctest
julia> isdir(homedir())
true

julia> isdir("not/a/directory")
false
```

See also: [`isfile`](@ref) and [`ispath`](@ref).
"""
isdir(st::StatStruct) = filemode(st) & 0xf000 == 0x4000

"""
    isblockdev(path) -> Bool

Return `true` if `path` is a block device, `false` otherwise.
"""
isblockdev(st::StatStruct) = filemode(st) & 0xf000 == 0x6000

"""
    isfile(path) -> Bool

Return `true` if `path` is a regular file, `false` otherwise.

# Examples
```jldoctest
julia> isfile(homedir())
false

julia> f = open("test_file.txt", "w");

julia> isfile(f)
true

julia> close(f); rm("test_file.txt")
```

See also: [`isdir`](@ref) and [`ispath`](@ref).
"""
isfile(st::StatStruct) = filemode(st) & 0xf000 == 0x8000

"""
    islink(path) -> Bool

Return `true` if `path` is a symbolic link, `false` otherwise.
"""
islink(st::StatStruct) = filemode(st) & 0xf000 == 0xa000

"""
    issocket(path) -> Bool

Return `true` if `path` is a socket, `false` otherwise.
"""
issocket(st::StatStruct) = filemode(st) & 0xf000 == 0xc000

# mode permission predicates

"""
    issetuid(path) -> Bool

Return `true` if `path` has the setuid flag set, `false` otherwise.
"""
issetuid(st::StatStruct) = (filemode(st) & 0o4000) > 0

"""
    issetgid(path) -> Bool

Return `true` if `path` has the setgid flag set, `false` otherwise.
"""
issetgid(st::StatStruct) = (filemode(st) & 0o2000) > 0

"""
    issticky(path) -> Bool

Return `true` if `path` has the sticky bit set, `false` otherwise.
"""
issticky(st::StatStruct) = (filemode(st) & 0o1000) > 0

"""
    uperm(file)

Get the permissions of the owner of the file as a bitfield of

| Value | Description        |
|:------|:-------------------|
| 01    | Execute Permission |
| 02    | Write Permission   |
| 04    | Read Permission    |

For allowed arguments, see [`stat`](@ref).
"""
uperm(st::StatStruct) = UInt8((filemode(st) >> 6) & 0x7)

"""
    gperm(file)

Like [`uperm`](@ref) but gets the permissions of the group owning the file.
"""
gperm(st::StatStruct) = UInt8((filemode(st) >> 3) & 0x7)

"""
    operm(file)

Like [`uperm`](@ref) but gets the permissions for people who neither own the file nor are a member of
the group owning the file
"""
operm(st::StatStruct) = UInt8((filemode(st)     ) & 0x7)

# mode predicate methods for file names

for f in Symbol[
    :ispath,
    :isfifo,
    :ischardev,
    :isdir,
    :isblockdev,
    :isfile,
    :issocket,
    :issetuid,
    :issetgid,
    :issticky,
    :uperm,
    :gperm,
    :operm,
    :filemode,
    :filesize,
    :mtime,
    :ctime,
]
    @eval ($f)(path...)  = ($f)(stat(path...))
end

islink(path...) = islink(lstat(path...))

# samefile can be used for files and directories: #11145#issuecomment-99511194
samefile(a::StatStruct, b::StatStruct) = a.device==b.device && a.inode==b.inode
function samefile(a::AbstractString, b::AbstractString)
    infoa = stat(a)
    infob = stat(b)
    if ispath(infoa) && ispath(infob)
        samefile(infoa, infob)
    else
        return false
    end
end

"""
    ismount(path) -> Bool

Return `true` if `path` is a mount point, `false` otherwise.
"""
function ismount(path...)
    path = joinpath(path...)
    isdir(path) || return false
    s1 = lstat(path)
    # Symbolic links cannot be mount points
    islink(s1) && return false
    parent_path = joinpath(path, "..")
    s2 = lstat(parent_path)
    # If a directory and its parent are on different devices,  then the
    # directory must be a mount point
    (s1.device != s2.device) && return true
    (s1.inode == s2.inode) && return true
    false
end

max_stars_repo_name: ajozefiak/julia

max_issues_repo_name: AsForMe/julia

max_forks_repo_name: AsForMe/julia


text:
@testset "weight init" begin

    @testset "Warnings" begin
        @test_logs (:warn, r"Identity mapping not possible with nin != nout!") idmapping(2,3)
        @test_logs (:warn, r"Identity mapping not possible with nin != nout!") idmapping(1,1,3,4)
        @test_logs (:warn, r"Identity mapping requires odd kernel sizes!") idmapping(2,1,3,3)
        @test_logs (:warn, r"Identity mapping requires odd kernel sizes!") idmapping(1,2,3,3)
    end

    @testset "Dense ID mapping" begin
        l = Dense(3,3, initW = idmapping)
        indata = reshape(collect(Float32, 1:9), 3, 3)
        @test l(indata) == indata

        @test_logs (:warn, "Identity mapping not possible with nin != nout! Got nin=3, nout=4.") l = Dense(3,4, initW=idmapping)
        @test nin(l) == 3
        @test nout(l) == 4
    end

    @testset "Conv ID mapping" begin
        import NaiveNASflux: weights
        nch = 5
        indata = reshape(collect(Float32, 1:nch*7*13), 7, 13, nch, 1)

        l = Conv((3,5), nch=>nch, init=idmapping, pad=(1,2))
        @test l(indata) == indata
        @test size(weights(l)) == size(weights(Conv((3,5), nch=>nch)))

        l = Conv((1,3), nch=>nch, init=idmapping, pad=(0,1))
        @test l(indata) == indata
        @test size(weights(l)) == size(weights(Conv((1,3), nch=>nch)))

        @test_logs (:warn, "Identity mapping not possible with nin != nout! Got nin=$nch, nout=$(nch+2).") l = Conv((3,5), nch=> nch+2, init=idmapping, pad=(1,2))
        @test nin(l) == nch
        @test nout(l) == nch+2
    end

end


text:
using PyPlot
include("../src/kepler_solver.jl")

# Define a constant of 1/3:
const third = 1.0/3.0

function test_elliptic()
# This routine runs a test of the kep_elliptic function in kepler_solver_elliptic.jl
# Define the central force constant in terms of AU and days:
k = (2.0*pi/365.25)^2
# Initial position at 1 AU:
x0 = [0.0,1.0,0.0]
#r0 = norm(x0)
r0 = sqrt(x0[1]*x0[1]+x0[2]*x0[2]+x0[3]*x0[3])
# Circular velocity:
vcirc = sqrt(k/r0)
# Define initial velocity at apastron:
v0 = [0.98*vcirc,0.0,0.0]  # The eccentricity is about ~2(1-v0/vcirc).
h = 18.0 # 18-day timesteps

s0::Float64 = 0.0
x = zeros(Float64,3)
v = zeros(Float64,3)
s::Float64 = 0.0
t = 0.0
ssave = zeros(2)
#nsteps = 10000000
nsteps = 1000000
#nsteps = 1000
xsave=zeros(Float64,12,nsteps)  # First is time (1,:); next three are position (2-4,:); next three velocity (5-7,:); then r (8,:), drdt (9,:); then beta (10,:); finally s & ds (11:12,:)
# Create a variable to store the state at the end of a step:
state=zeros(Float64,12)

xsave[1,1]=0.0
for j=1:3
  xsave[j+1,1]=x0[j]
  xsave[j+4,1]=v0[j]
end
# Save beta:
xsave[8,1]=r0
beta0 = 2.0*k/r0-dot(v0,v0)
xsave[10,1]=beta0
#@inbounds for i=2:nsteps
for i=2:nsteps
#  x,v,r,drdt,s,beta,iter=kep_elliptic(x0,v0,r0,k,h,beta0,s0)
#  iter = kep_elliptic!(x0,v0,r0,k,h,beta0,s0,state)
  iter = kep_ell_hyp!(x0,v0,r0,k,h,beta0,s0,state)
  s = state[11]
  ds = state[12]
  if iter > 2
   println("iter: ",iter," ds: ",ds)
  end
# Increment time by h:
  state[1] += h
  for j=1:12
    xsave[j,i]=state[j]
  end
  if i >= 4
# Quadratic prediction for s of next step:
    s0 = 3.0*s - 3.0*ssave[1] + ssave[2]
    ssave[2] = ssave[1]
    ssave[1] = s
  else
    s0 = s
    ssave[2]=ssave[1]
    ssave[1]=s
  end
# Now, proceed to next step:
  for j=1:3
    x0[j]=state[1+j]
    v0[j]=state[4+j]
  end
  r0=state[8]
  beta0=state[10]
#  println(i,x,v)
end
#plot(xsave[1],xsave[2])
return xsave
end


text:
struct DataSetGroup
    name::AbstractString
    p::Ptr{Nothing}
    l::Int
    meta::Dict{AbstractString, AbstractString}
end
Base.length(dsg::DataSetGroup) = MDAL_G_datasetCount(dsg.p)

function DataSetGroup(p::Ptr{Nothing}, i::Int)
    @info "Getting group $i"
    group = MDAL_M_datasetGroup(p, i)
    name = unsafe_string(MDAL_G_name(group))
    l = MDAL_G_datasetCount(group)
    c = MDAL_G_metadataCount(group)
    @info "$c metadata in datasetgroup"
    m = Dict(unsafe_string(MDAL_G_metadataKey(group, i)) =>
             unsafe_string(MDAL_G_metadataValue(group, i)) for i in 0:c-1)
    DataSetGroup(name, group, l, m)
end

DataSetGroup(m::MDALMesh, i::Int) = DataSetGroup(m.mdalm, i)
function DataSetGroups(m::MDALMesh)
    c = MDAL_M_datasetGroupCount(m.mdalm)
    @info "$c groups in mesh"
    (DataSetGroup(m, i) for i in 1:c-1)
end

struct DataSet
    p::Ptr{Nothing}
end
Base.length(ds::DataSet) = MDAL_D_valueCount(ds.p)
Base.isvalid(ds::DataSet) = MDAL_D_isValid(ds.p)
Base.time(ds::DataSet) = MDAL_D_time(ds.p)

DataSet(m::DataSetGroup, i::Int) = DataSet(m.p, i)
function DataSet(p::Ptr{Nothing}, i::Int)
    ds = MDAL_G_dataset(p, i)
    DataSet(ds)
end

function DataSets(p::DataSetGroup)
    c = MDAL_G_datasetCount(p.p)
    @info "$c datasets in group"
    (DataSet(p, i) for i in 1:c-1)
end


text:
using Effects.List, Test, FunctionalCollections

@test list() do
  x = choose(1:10)
  isodd(x) || backtrack()
  x^2
end == (1:2:9).^2

@test list() do
  x = choose([1, 2, 3])
  y = choose([1, 2, 3])
  x^2 + y == 7 || backtrack()
  (x, y)
end == [(2, 3)]

triples = list() do
  N = 20
  i = choose(1:N)
  j = choose(i:N)
  k = choose(j:N)
  i^2 + j^2 == k^2 || backtrack()
  (i, j, k)
end

@test first(triples) == (3, 4, 5)
@test length(triples) == 6

function nqueens(n)
  ps = @Persistent []
  for i = 1:n
    next = choose(1:n)
    next in ps && backtrack() # same row
    any(abs(next-p) == (i-i′) for (i′, p) in enumerate(ps)) && backtrack() # same diagonal
    ps = push(ps, next)
  end
  return ps
end

@test length(list(() -> nqueens(8))) == 92


text:
# Reducing methods

for (mod, fname) in ((:Base, :sum), (:Base, :prod), (:Base, :maximum), (:Base, :minimum), (:Statistics, :mean))
    _fname = Symbol('_', fname)
    @eval begin
        # Returns a scalar
        @inline ($mod.$fname)(A::AbDimArray) = ($mod.$fname)(parent(A))
        # Returns a reduced array
        @inline ($mod.$_fname)(A::AbstractArray, dims::AllDimensions) =
            rebuild(A, ($mod.$_fname)(parent(A), dimnum(A, dims)), reducedims(A, dims))
        @inline ($mod.$_fname)(f, A::AbstractArray, dims::AllDimensions) =
            rebuild(A, ($mod.$_fname)(f, parent(A), dimnum(A, dims)), reducedims(A, dims))
        @inline ($mod.$_fname)(A::AbDimArray, dims::Union{Int,Base.Dims}) =
            rebuild(A, ($mod.$_fname)(parent(A), dims), reducedims(A, dims))
        @inline ($mod.$_fname)(f, A::AbDimArray, dims::Union{Int,Base.Dims}) =
            rebuild(A, ($mod.$_fname)(f, parent(A), dims), reducedims(A, dims))
    end
end

for (mod, fname) in ((:Statistics, :std), (:Statistics, :var))
    _fname = Symbol('_', fname)
    @eval begin
        # Returns a scalar
        @inline ($mod.$fname)(A::AbDimArray) = ($mod.$fname)(parent(A))
        # Returns a reduced array
        @inline ($mod.$_fname)(A::AbstractArray, corrected::Bool, mean, dims::AllDimensions) =
            rebuild(A, ($mod.$_fname)(A, corrected, mean, dimnum(A, dims)), reducedims(A, dims))
        @inline ($mod.$_fname)(A::AbDimArray, corrected::Bool, mean, dims::Union{Int,Base.Dims}) =
            rebuild(A, ($mod.$_fname)(parent(A), corrected, mean, dims), reducedims(A, dims))
    end
end

Statistics.median(A::AbDimArray) = Statistics.median(parent(A))
Statistics._median(A::AbstractArray, dims::AllDimensions) =
    rebuild(A, Statistics._median(parent(A), dimnum(A, dims)), reducedims(A, dims))
Statistics._median(A::AbDimArray, dims::Union{Int,Base.Dims}) =
    rebuild(A, Statistics._median(parent(A), dims), reducedims(A, dims))

Base._mapreduce_dim(f, op, nt::NamedTuple{(),<:Tuple}, A::AbstractArray, dims::AllDimensions) =
    rebuild(A, Base._mapreduce_dim(f, op, nt, parent(A), dimnum(A, dims)), reducedims(A, dims))
Base._mapreduce_dim(f, op, nt::NamedTuple{(),<:Tuple}, A::AbDimArray, dims::Union{Int,Base.Dims}) =
    rebuild(A, Base._mapreduce_dim(f, op, nt, parent(A), dimnum(A, dims)), reducedims(A, dims))
Base._mapreduce_dim(f, op, nt::NamedTuple{(),<:Tuple}, A::AbDimArray, dims::Colon) =
    Base._mapreduce_dim(f, op, nt, parent(A), dims)

# TODO: Unfortunately Base/accumulate.jl kwargs methods all force dims to be Integer.
# accumulate wont work unless that is relaxed, or we copy half of the file here.
Base._accumulate!(op, B, A, dims::AllDimensions, init::Union{Nothing, Some}) =
    Base._accumulate!(op, B, A, dimnum(A, dims), init)

Base._extrema_dims(f, A::AbstractArray, dims::AllDimensions) = begin
    dnums = dimnum(A, dims)
    rebuild(A, Base._extrema_dims(f, parent(A), dnums), reducedims(A, dnums))
end

# Dimension dropping

Base._dropdims(A::AbstractArray, dim::DimOrDimType) =
    rebuildsliced(A, Base._dropdims(A, dimnum(A, dim)), dims2indices(A, basetypeof(dim)(1)))
Base._dropdims(A::AbstractArray, dims::AbDimTuple) =
    rebuildsliced(A, Base._dropdims(A, dimnum(A, dims)),
                  dims2indices(A, Tuple((basetypeof(d)(1) for d in dims))))

# Function application

@inline Base.map(f, A::AbDimArray) = rebuild(A, map(f, parent(A)), dims(A))

Base.mapslices(f, A::AbDimArray; dims=1, kwargs...) = begin
    dimnums = dimnum(A, dims)
    data = mapslices(f, parent(A); dims=dimnums, kwargs...)
    rebuild(A, data, reducedims(A, DimensionalData.dims(A, dimnums)))
end

# This is copied from base as we can't efficiently wrap this function
# through the kwarg with a rebuild in the generator. Doing it this way
# wierdly makes it faster toeuse a dim than an integer.
if VERSION > v"1.1-"
    Base.eachslice(A::AbDimArray; dims=1, kwargs...) = begin
        if dims isa Tuple && length(dims) == 1
            throw(ArgumentError("only single dimensions are supported"))
        end
        dim = first(dimnum(A, dims))
        dim <= ndims(A) || throw(DimensionMismatch("A doesn't have $dim dimensions"))
        idx1, idx2 = ntuple(d->(:), dim-1), ntuple(d->(:), ndims(A)-dim)
        return (view(A, idx1..., i, idx2...) for i in axes(A, dim))
    end
end

for fname in (:cor, :cov)
    @eval Statistics.$fname(A::AbDimArray{T,2}; dims=1, kwargs...) where T = begin
        newdata = Statistics.$fname(parent(A); dims=dimnum(A, dims), kwargs...)
        I = dims2indices(A, dims, 1)
        newdims, newrefdims = slicedims(A, I)
        rebuild(A, newdata, (newdims[1], newdims[1]), newrefdims)
    end
end

# Reverse

@inline Base.reverse(A::AbDimArray{T,N}; dims=1) where {T,N} = begin
    dnum = dimnum(A, dims)
    # Reverse the dimension. TODO: make this type stable
    newdims = reversearray(DimensionalData.dims(A), dnum)
    # Reverse the data
    newdata = reverse(parent(A); dims=dnum)
    rebuild(A, newdata, newdims, refdims(A))
end

@inline reversearray(dimstorev::Tuple, dnum) = begin
    dim = dimstorev[end]
    if length(dimstorev) == dnum
        dim = rebuild(dim, val(dim), reversearray(grid(dim)))
    end
    (reversearray(Base.front(dimstorev), dnum)..., dim)
end
@inline reversearray(dims::Tuple{}, i) = ()

# Dimension reordering

for (pkg, fname) in [(:Base, :permutedims), (:Base, :adjoint),
                     (:Base, :transpose), (:LinearAlgebra, :Transpose)]
    @eval begin
        @inline $pkg.$fname(A::AbDimArray{T,2}) where T =
            rebuild(A, $fname(parent(A)), reverse(dims(A)), refdims(A))
    end
end

for fname in [:permutedims, :PermutedDimsArray]
    @eval begin
        @inline Base.$fname(A::AbDimArray{T,N}, perm) where {T,N} =
            rebuild(A, $fname(parent(A), dimnum(A, perm)), permutedims(dims(A), perm))
    end
end

# Indices

Base.firstindex(A::AbstractArray, d::DimOrDimType) = firstindex(A, dimnum(A, d))
Base.lastindex(A::AbstractArray, d::DimOrDimType) = lastindex(A, dimnum(A, d))

# Concatenation

Base._cat(catdims::Union{AbDim,AbDimTuple}, A::DimensionalArray...) = begin
    A1 = first(A)
    checkdims(A...)
    if all(hasdim(A1, catdims))
        dnum = dimnum(A1, catdims)
        newdims = map(zip(map(dims, A)...)) do ds
            basetypeof(catdims) <: basetypeof(ds[1]) ? vcat(ds...) : ds[1]
        end
        rebuild(A1, Base._cat(dnum, map(parent, A)...), )
    else
        add_dims = if (catdims isa Tuple)
            Tuple(d for d in catdims if !hasdim(A1, d))
        else
            (catdims,)
        end
        dnum = ndims(A1) + length(add_dims)
        newA = Base._cat(dnum, map(parent, A)...)
        newdims = (dims(A1)..., add_dims...)
        rebuild(A1, newA, formatdims(size(newA), newdims))
    end
end

Base.vcat(A::AbDimArray...) = begin
    checkdims(A...)
    rebuild(A, vcat(map(val, A)), vcat(map(A, dims)))
end
Base.vcat(dims::AbDim...) =
    rebuild(dims[1], vcat(map(val, dims)), vcat(map(grid, dims)))

Base.vcat(grids::Grid...) = first(grids)
Base.vcat(grids::BoundedGrid...) =
    rebuild(grids[1]; bounds=(bounds(grids[1])[1], bounds(grids[end])[end]))

checkdims(A::AbstractArray...) = checkdims(map(dims, A)...)
checkdims(dims::AbDimTuple...) = map(d -> checkdims(dims[1], d), dims)
checkdims(d1::AbDimTuple, d2::AbDimTuple) = map(checkdims, d1, d2)
checkdims(d1::AbDim, d2::AbDim) =
    basetypeof(d2) <: basetypeof(d1) || error("Dims differ: $(bastypeof(d1)), $(basetypeof(d2))")

# Index breaking

# TODO: change the index and traits of the reduced dimension
# and return a DimensionalArray.
Base.unique(A::AbDimArray{<:Any,1}) = unique(parent(A))
Base.unique(A::AbDimArray; dims::DimOrDimType) =
    unique(parent(A); dims=dimnum(A, dims))

# TODO cov, cor mapslices, eachslice, reverse, sort and sort! need _methods without kwargs in base so
# we can dispatch on dims. Instead we dispatch on array type for now, which means
# these aren't usefull unless you inherit from AbDimArray.


text:
function generate_problem_dictionary(path_to_parameters_file::String)::Dict{String,Any}

    # initialize -
    problem_dictionary = Dict{String,Any}()

    try

        # load the TOML parameters file -
        toml_dictionary = TOML.parsefile(path_to_parameters_file)["biophysical_constants"]

        # setup the initial condition array -
        initial_condition_array = [
            0.0 ;   # 1 mRNA
            5.0 ;   # nM and constant TODONE: gene concentration goes here -
            0.0 ;   # 3 I = we'll fill this in the execute script
        ]
        # TODO: calculate the mRNA_degradation_constant
       # s decided not to use this bc couldnt confirm a good reasoning from the papers for how to use the half life for the degradation constant
        mRNA_degradation_constant = log(2)/toml_dictionary["mRNA_half_life_in_min"] #converts mRNA to a constant
        #mRNA_degradation_constant = 3.75/3600
        # TODONE: VMAX for transcription -
        RNAPII_concentration= toml_dictionary["RNAPII_concentration"]                               # units: nM
        transcription_elongation_rate= toml_dictionary["transcription_elongation_rate"]                       # units: nt/s
        gene_length_in_nt= toml_dictionary["gene_length_in_nt"]                                   # units: nt

        VMAX_X = RNAPII_concentration*(transcription_elongation_rate/gene_length_in_nt)  # nM*(nt/s/nt) = nM/s

        # TODO: Stuff that I'm forgetting?
        # check to see whenever problem[""] is used, I think only this one wasn't defined
        problem_dictionary["inducer_cooperativity_parameter"] =toml_dictionary["inducer_cooperativity_parameter"]
        problem_dictionary["transcription_elongation_rate"] =toml_dictionary["transcription_elongation_rate"]
        # --- PUT STUFF INTO problem_dictionary ----
        problem_dictionary["transcription_time_constant"] = toml_dictionary["transcription_time_constant"]
        problem_dictionary["transcription_saturation_constant"] = toml_dictionary["transcription_saturation_constant"]
        problem_dictionary["E1"] = toml_dictionary["energy_promoter_state_1"]
        problem_dictionary["E2"] = toml_dictionary["energy_promoter_state_2"]
        problem_dictionary["inducer_dissociation_constant"] = toml_dictionary["inducer_dissociation_constant"]
        problem_dictionary["ideal_gas_constant_R"] = 0.008314  # kJ/mol-K
        problem_dictionary["temperature_K"] = (273.15+37)
        problem_dictionary["initial_condition_array"] = initial_condition_array
        problem_dictionary["mRNA_degradation_constant"] = mRNA_degradation_constant
        problem_dictionary["maximum_transcription_velocity"] = VMAX_X  # nM/s eq 7

        # return -
        return problem_dictionary

    catch error
        throw(error)
    end
end


text:
using Test
@isdefined(TestHelpers) || include("TestHelpers.jl")

TestHelpers.@timed_testset "batch_processing_helpers/" begin
    include("batch_processing_helpers/unit.jl")
    include("batch_processing_helpers/integration.jl")
end


text:
using DataStructures

#=
    Partition.

    Write code to partition a linked list around a value x,
    such that all nodes less than x come before all nodes
    greater than or equal to x. If x is contained within the list,
    the values of x can appear anywhere in the "right partition";
    it does not need to appear between the left and right partitions.

    EXAMPLE

    Input: 3 → 5 → 8 → 5 → 10 → 2 → 1 [partition = 5]
    Output: 3 → 1 → 2 → 10 → 5 → 5 → 8
=#

function partition(lst::LinkedList, value::Integer)::LinkedList
    let iter = iterate(lst)
        head, tail = nil(), nil()
        while iter != nothing
            v, s = iter

            if v < partition
                head = cons(v, head)
            else
                tail = cons(v, tail)
            end

            iter = iterate(lst, s)
        end
        new = cat(head, tail)
    end
    return new
end

partition(lst, value)

lst = list(3, 5, 8, 5, 10, 2, 1)
value = 5

max_stars_repo_name: bmoretz/Daily-Coding-Problem

max_issues_repo_name: bmoretz/Daily-Coding-Problem

max_forks_repo_name: bmoretz/Daily-Coding-Problem


text:
using Test
@isdefined(TestHelpers) || include("../TestHelpers.jl")

TestHelpers.@timed_testset "nets/" begin
    include("nets/sequential.jl")
end


text:
type MTsource{T1,T2, T3<:Real}
    mt :: Array{T1} # moment tensor
    loc :: Array{T2} # source's physical location in z and x axis
    nloc :: Array{Int64} # source's discrete location in z and x axis without PML boundary
    BDnloc :: Array{Int64} # source's discrete location in z and x axis with PML boundary. It is directly related to source type.
    ot :: T3 # source origin time
    not :: Int64 # source's discrete origin time
    waveform :: Array{Float64} # waveform for a single source
end

type SFsource{T1, T2<:Real}
    coeff :: Array{Int64} # single force [vx vz] or [vx vy vz]
    loc :: Array{T1} # source's physical location in z and x axis
    nloc :: Array{Int64} # source's discrete location in z and x axis without PML boundary
    BDnloc :: Array{Int64} # source's discrete location in z and x axis with PML boundary. It is directly related to source type.
    ot :: T2 # source origin time
    not :: Int64 # source's discrete origin time
    waveform :: Array{Float64} # waveform for a single source
end

function CreateWaveform(sn::Int64, not::Int64, nT::Int64, waveform::Array{Float64})
  SourceWaveform = zeros(nT)
  twf = length(waveform)
  if not < nT && not > 0
    if not+twf-1 <= nT
      SourceWaveform[not:not+twf-1] = SourceWaveform[not:not+twf-1] + waveform
    elseif not+twf-1 > nT
      SourceWaveform[not:nT] = SourceWaveform[not:nT] + waveform[1:nT-not+1]
    end
  elseif not >= nT
    errormessage = string("check No.", string(sn), " source's starting time")
    error(errormessage)
  elseif not == 0
    if twf <= nT
      SourceWaveform[1:twf] = SourceWaveform[1:twf] + waveform
    elseif twf > nT
      SourceWaveform[1:nT] = SourceWaveform[1:nT] + waveform[1:nT]
    end
  end
  return SourceWaveform
end

function BDnlocation(ext::Int64, snz::Int64, snx::Int64, iflag::Int64)
    # snx and snz : source discrete location in x and z direction

    if iflag == 1 #free surface
        BDnloc = [snz, snx+ext]
    elseif iflag == 2 #unlimited medium
        BDnloc = [snz+ext, snx+ext]
    end
    return BDnloc
end

function BDnlocation(ext::Int64, snz::Int64, snx::Int64, sny::Int64, iflag::Int64)
    # snx and snz : source discrete location in x and z direction

    if iflag == 1 #free surface
        BDnloc = [snz, snx+ext, sny+ext]
    elseif iflag == 2 #unlimited medium
        BDnloc = [snz+ext, snx+ext, sny+ext]
    end
    return BDnloc
end

function CreateMtSource{T1,T2,T3<:Real}(
    loc :: Array{T1},
    ot :: T2,
    mt :: Array{T3,1}, # moment tensor (mxx,mzz,mxz)
    waveform :: Array{Float64},
    medium :: Union{elastic2d,acoustic2d},
    nwf :: Union{nelwf2d,nacwf2d})

    if length(mt) != 3
      error("Three elements in 2d moment tensor (mxx, mzz, mxz)")
    end

    nloc = [Int64(round(loc[1]/medium.dz)),  Int64(round(loc[2]/medium.dx))]
    BDnloc = BDnlocation(medium.ext, nloc[1], nloc[2], medium.iflag)
    otinput = ot
    not = Int64(round(otinput/medium.dt))
    SourceWaveform = CreateWaveform(1, not, medium.nT, waveform)

    Sou = MTsource(mt, loc, nloc, BDnloc, otinput, not, SourceWaveform)
  end

  function CreateMtSource{T1,T2,T3<:Real}(
      loc::Array{T1},
      ot::T2,
      mt :: Array{T3,1}, # moment tensor (mxx,myy,mzz,mxy,mxz,myz)
      waveform::Array{Float64},
      medium::Union{elastic3d,acoustic3d},
      nwf::Union{nelwf3d,nacwf3d})

      if length(mt) != 6
        error("Six elements in 3d moment tensor (mxx, myy, mzz, mxy, mxz, myz)")
      end

      nloc = [Int64(round(loc[1]/medium.dz)),  Int64(round(loc[2]/medium.dx)), Int64(round(loc[3]/medium.dy))]
      BDnloc = BDnlocation(medium.ext, nloc[1], nloc[2], nloc[3], medium.iflag)
      otinput = ot
      not = Int64(round(otinput/medium.dt))
      SourceWaveform = CreateWaveform(1, not, medium.nT, waveform)

      Sou = MTsource(mt, loc, nloc, BDnloc, otinput, not, SourceWaveform)
    end

    function CreateSfSource{T1,T2<:Real}(
        loc :: Array{T1},
        ot :: T2,
        coeff :: Array{Int64}, # moment tensor (mxx,mzz,mxz)
        waveform :: Array{Float64},
        medium :: Union{elastic2d,acoustic2d},
        nwf :: Union{nelwf2d,nacwf2d})

        if length(coeff) != 2
          error("Two elements in 2d coefficient (vx vz).")
        end
        if (!(coeff[1]==0 || coeff[1]==1) || !(coeff[2]==0 || coeff[2]==1))==true
          error("coefficient can only be chosen from 0 and 1.")
        end

        nloc = [Int64(round(loc[1]/medium.dz)),  Int64(round(loc[2]/medium.dx))]
        BDnloc = BDnlocation(medium.ext, nloc[1], nloc[2], medium.iflag)
        otinput = ot
        not = Int64(round(otinput/medium.dt))
        SourceWaveform = CreateWaveform(1, not, medium.nT, waveform)

        Sou = SFsource(coeff, loc, nloc, BDnloc, otinput, not, SourceWaveform)
      end

      function CreateSfSource{T1,T2<:Real}(
          loc::Array{T1},
          ot::T2,
          coeff :: Array{Int64,1}, # moment tensor (mxx,myy,mzz,mxy,mxz,myz)
          waveform::Array{Float64},
          medium::Union{elastic3d,acoustic3d},
          nwf::Union{nelwf3d,nacwf3d})

          if length(coeff) != 3
            error("Six elements in 3d moment tensor (mxx, myy, mzz, mxy, mxz, myz)")
          end
          if (!(coeff[1]==0 || coeff[1]==1) || !(coeff[2]==0 || coeff[2]==1) || !(coeff[3]==0 || coeff[3]==1))==true
            error("coefficient can only be chosen from 0 and 1.")
          end

          nloc = [Int64(round(loc[1]/medium.dz)),  Int64(round(loc[2]/medium.dx)), Int64(round(loc[3]/medium.dy))]
          BDnloc = BDnlocation(medium.ext, nloc[1], nloc[2], nloc[3], medium.iflag)
          otinput = ot
          not = Int64(round(otinput/medium.dt))
          SourceWaveform = CreateWaveform(1, not, medium.nT, waveform)

          Sou = SFsource(coeff, loc, nloc, BDnloc, otinput, not, SourceWaveform)
        end

#=== MT source ===#
function initMTsource{T1,T2,T3<:Real}(
    sn::Int64,
    loc::Array{T1}, # location
    ot::Union{T2,Array{T2}}, # origin time
    mt:: Array{T3}, # type
    waveform::Array{Float64},
    model::Union{elmod2d,acmod2d,elmod3d,acmod3d})

    sou = Array{MTsource}(sn)
    for i in 1:sn
        sou[i] = CreateMtSource(
        loc[i,:],
        ot[i],
        mt[i,:],
        waveform[:,i],
        model.medium,
        model.nwf)
      end
    return sou
end

#=== SF source ===#
function initSFsource{T1,T2<:Real}(
    sn::Int64,
    loc::Array{T1}, # location
    ot::Union{T2,Array{T2}}, # origin time
    coeff:: Array{Int64}, # type
    waveform::Array{Float64},
    model::Union{elmod2d,acmod2d,elmod3d,acmod3d})

    sou = Array{SFsource}(sn)
    for i in 1:sn
        sou[i] = CreateSfSource(
        loc[i,:],
        ot[i],
        coeff[i,:],
        waveform[:,i],
        model.medium,
        model.nwf)
      end
    return sou
end


text:
using HyperGraphs, SimpleTraits

# making sure types are correctly traited
@test istrait(IsOriented{ChemicalHyperGraph})
@test istrait(IsOriented{ChemicalHyperEdge})
@test !istrait(IsOriented{HyperGraph})
@test !istrait(IsOriented{HyperEdge})
@test istrait(IsWeighted{ChemicalHyperGraph})
@test istrait(IsWeighted{ChemicalHyperEdge})
@test !istrait(IsWeighted{HyperGraph})
@test !istrait(IsWeighted{HyperEdge})

## trait functions ##

# chemical hypergraphs and hyperedges are oriented and weighted
@test isoriented(ChemicalHyperGraph()) && isoriented(ChemicalHyperGraph)
@test isoriented(ChemicalHyperEdge()) && isoriented(ChemicalHyperEdge)
@test isweighted(ChemicalHyperGraph()) && isoriented(ChemicalHyperGraph)
@test isweighted(ChemicalHyperEdge()) && isoriented(ChemicalHyperEdge)

# hypergraphs and hyperedges are unoriented and unweighted
@test !isoriented(HyperGraph()) && !isoriented(HyperGraph)
@test !isoriented(HyperEdge()) && !isoriented(HyperEdge)
@test !isweighted(HyperGraph()) && !isoriented(HyperGraph)
@test !isweighted(HyperEdge()) && !isoriented(HyperEdge)

## edge cases ##

# other types and object should not be traited
@test !isoriented(String) && !isoriented(3)
@test !isweighted(String) && !isweighted(3)
@test_throws TypeError istrait(IsOriented{Int})


text:
using DelimitedFiles
using Statistics
using Plots

include("exact.jl")

A, B, C = indexGrid(100, 100, 100)

Lx = 5
Ly = 4
Lz = 3

N = 10

t = 0

F = modeFrequency.(A, B, C, Lx, Ly, Lz)

idx = sortperm(F[:])

f = F[idx[2:(N + 1)]]

for n in 1:length(f)
    writedlm("frequencies/f_$n.csv", f[n])
end

a = A[idx[2:(N + 1)]]
b = B[idx[2:(N + 1)]]
c = C[idx[2:(N + 1)]]

norms = zeros(N)

for d in 0:(N - 1)

    linedata = readdlm("export/data_$d.csv", ',')

    femSolution = real.(
        (Float64.(linedata[2:end, 1]) .+ im * Float64.(linedata[2:end, 2])) *
        exp(im * 2π * f[d + 1] * t)
    )

    x = Float64.(linedata[2:end, 3])
    y = Float64.(linedata[2:end, 4])
    z = Float64.(linedata[2:end, 5])

    exactSolution = mode.(
        x .+ 0.5 * Lx, y .+ 0.5 * Ly, z .+ 0.5 * Lz,
        a[d + 1], b[d + 1], c[d + 1],
        Lx, Ly, Lz
    ) * cospi(2 * f[d + 1] * t)

    me, ie = findmax(exactSolution)

    femSolution = femSolution * me / femSolution[ie]

    errorField = femSolution - exactSolution

    norms[d + 1] = sqrt(sum(abs2.(errorField)))

    dBError = 20.0 * log10.(abs.(femSolution ./ exactSolution))

    writedlm(
        "error/error_$d.csv",
        ["error_field" "dB_error_field" "x" "y" "z"; [errorField dBError x y z]],
        ','
    )

end

bar(
    1:N,
    norms,
    legend=false,
    framestyle=:box,
    xlabel="Mode frequency [Hz]",
    ylabel="Norm of the Error Field [Pa]",
    yscale=:log10,
    xticks=(1:N, round.(f, digits=3))
)

max_stars_repo_name: CrocoDuckoDucks/Acoustics-Models

max_issues_repo_name: CrocoDuckoDucks/Acoustics-Models

max_forks_repo_name: CrocoDuckoDucks/Acoustics-Models


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

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

# ╔═╡ b46a1e96-d4e9-47a0-818b-d69146434db3
using Revise; using Utils

# ╔═╡ a367e7ea-59ae-11eb-3ae8-7935e1bab240
using PlutoUI, DataFramesMeta, Glob, AstroImages, HDF5, Dates, Printf, StatsPlots, CSV, CCDReduction, FITSIO, DataFrames, RecursiveArrayTools, BenchmarkTools, Statistics, PaddedViews, PlotlyBase, KernelDensity, PyCall

# ╔═╡ 78636301-9ad5-4562-8526-fb92a5e88c84
TableOfContents()

# ╔═╡ 3fc21e1c-355f-40a6-a846-2acea97bdd94
md"""### $(@bind run_GPTS_WLC_corner CheckBox()) GPTS WLC Corner Plot

Visualize GPT data stored in
"""

# ╔═╡ bac5f350-ce60-4d18-bdbd-41f51dea7837
if run_GPTS_WLC_corner
	# Load
	gpts_pkl_to_dict = load_pickle(
		"Projects/HATP23b/data_detrending/out_c/HATP23b/hp23b_170609_custom/white-light/BMA_posteriors.pkl"
	)
	println(keys(gpts_pkl_to_dict))
	const gpts_params = ["p", "t0", "P", "rho", "inc", "b", "aRs", "q1"]
	gpts_params_dict = sub_dict(gpts_pkl_to_dict, gpts_params)

	df_gpts = DataFrame(gpts_params_dict)

	let
		corner(df_gpts, Symbol.(gpts_params), size=(15_000, 30_000))
		png("/home/mango/Desktop/test")
	end
end

# ╔═╡ 2cab7acb-65e8-477b-88ab-16ed23c14c4e
md"""
## Plot configs
"""

# ╔═╡ b16c144e-064f-49a7-9812-52cd695e7b4b
theme(:dark); plotly()

# ╔═╡ 16df4db5-a788-4b5c-85d5-23937a1a890a
const pal = palette(["#f5d300", "#d55e00", "#08f7fe", "#0173b2", "#029e73"])

# ╔═╡ d2951668-8a8f-4cae-a011-4390708933ac
default(
	titlefont = "Lato",
	guidefont = "Lato",
	markerstrokewidth = 0,
	palette=pal,
)

# ╔═╡ 7cb0686e-00db-4016-8c55-a69a4c260ebe
md"""
## Packages
"""

# ╔═╡ Cell order:
# ╟─78636301-9ad5-4562-8526-fb92a5e88c84
# ╟─3fc21e1c-355f-40a6-a846-2acea97bdd94
# ╠═bac5f350-ce60-4d18-bdbd-41f51dea7837
# ╟─2cab7acb-65e8-477b-88ab-16ed23c14c4e
# ╠═b16c144e-064f-49a7-9812-52cd695e7b4b
# ╠═d2951668-8a8f-4cae-a011-4390708933ac
# ╠═16df4db5-a788-4b5c-85d5-23937a1a890a
# ╟─7cb0686e-00db-4016-8c55-a69a4c260ebe
# ╠═b46a1e96-d4e9-47a0-818b-d69146434db3
# ╠═a367e7ea-59ae-11eb-3ae8-7935e1bab240


text:
# Test of binary heaps

# auxiliary functions

function heap_values(h::MutableBinaryHeap{VT,O}) where {VT,O}
    n = length(h)
    nodes = h.nodes
    @assert length(nodes) == n
    vs = Vector{VT}(undef, n)
    for i = 1 : n
        vs[i] = nodes[i].value
    end
    vs
end

function list_values(h::MutableBinaryHeap{VT,O}) where {VT,O}
    n = length(h)
    nodes = h.nodes
    nodemap = h.node_map
    vs = Vector{VT}()
    for i = 1 : length(nodemap)
        id = nodemap[i]
        if id > 0
            push!(vs, nodes[id].value)
        end
    end
    vs
end

function verify_heap(h::MutableBinaryHeap{VT,O}) where {VT,O}
    ord = h.ordering
    nodes = h.nodes
    n = length(h)
    m = div(n,2)
    for i = 1 : m
        v = nodes[i].value
        lc = i * 2
        if lc <= n
            if Base.lt(ord, nodes[lc].value, v)
                return false
            end
        end
        rc = lc + 1
        if rc <= n
            if Base.lt(ord, nodes[rc].value, v)
                return false
            end
        end
    end
    return true
end

@testset "MutableBinheap" begin

    vs = [4, 1, 3, 2, 16, 9, 10, 14, 8, 7]
    vs2 = collect(enumerate(vs))
    ordering = Base.Order.By(last)

    @testset "construct heap" begin
        MutableBinaryHeap{Int, Base.ForwardOrdering}()
        MutableBinaryHeap{Int, Base.ForwardOrdering}(vs)

        MutableBinaryHeap{Int, Base.ReverseOrdering}()
        MutableBinaryHeap{Int, Base.ReverseOrdering}(vs)

        MutableBinaryMinHeap{Int}()
        MutableBinaryMinHeap{Int}(vs)
        MutableBinaryMinHeap(vs)

        MutableBinaryMaxHeap{Int}()
        MutableBinaryMaxHeap{Int}(vs)
        MutableBinaryMaxHeap(vs)

        MutableBinaryHeap{eltype(vs2)}(ordering)
        MutableBinaryHeap{eltype(vs2)}(ordering, vs2)
        MutableBinaryHeap(ordering, vs2)

        @test true
    end

    @testset "basic tests" begin
        h = MutableBinaryMinHeap{Int}()

        @test length(h) == 0
        @test isempty(h)
        @test eltype(h) == Int
        @test eltype(typeof(h)) == Int
    end

    @testset "make mutable binary minheap" begin
        h = MutableBinaryMinHeap(vs)

        @test length(h) == 10
        @test !isempty(h)
        @test first(h) == 1
        @test isequal(list_values(h), vs)
        @test isequal(heap_values(h), [1, 2, 3, 4, 7, 9, 10, 14, 8, 16])
        @test sizehint!(h, 100) === h
    end

    @testset "make mutable binary maxheap" begin
        h = MutableBinaryMaxHeap(vs)

        @test length(h) == 10
        @test !isempty(h)
        @test first(h) == 16
        @test isequal(list_values(h), vs)
        @test isequal(heap_values(h), [16, 14, 10, 8, 7, 3, 9, 1, 4, 2])
        @test sizehint!(h, 100) === h
    end

    @testset "make mutable binary custom ordering heap" begin
        h = MutableBinaryHeap(ordering, vs2)

        @test length(h) == 10
        @test !isempty(h)
        @test first(h) == (2, 1)
        @test isequal(list_values(h), vs2)
        @test isequal(heap_values(h), [(2, 1), (4, 2), (3, 3), (1, 4), (10, 7), (6, 9), (7, 10), (8, 14), (9, 8), (5, 16)])
        @test sizehint!(h, 100) === h
    end

    @testset "hmin / push! / pop!" begin
        hmin = MutableBinaryMinHeap{Int}()
        @test length(hmin) == 0
        @test isempty(hmin)

        # test push!
        ss = Any[
            [4],
            [1, 4],
            [1, 4, 3],
            [1, 2, 3, 4],
            [1, 2, 3, 4, 16],
            [1, 2, 3, 4, 16, 9],
            [1, 2, 3, 4, 16, 9, 10],
            [1, 2, 3, 4, 16, 9, 10, 14],
            [1, 2, 3, 4, 16, 9, 10, 14, 8],
            [1, 2, 3, 4, 7, 9, 10, 14, 8, 16]]
        for i = 1 : length(vs)
            ia = push!(hmin, vs[i])
            @test ia == i
            @test length(hmin) == i
            @test !isempty(hmin)
            @test isequal(list_values(hmin), vs[1:i])
            @test isequal(heap_values(hmin), ss[i])
        end

        # test pop!
        @test isequal(extract_all!(hmin), [1, 2, 3, 4, 7, 8, 9, 10, 14, 16])
        @test isempty(hmin)
    end

    @testset "hmax / push! / pop!" begin
        hmax = MutableBinaryMaxHeap{Int}()
        @test length(hmax) == 0
        @test isempty(hmax)

        # test push!
        ss = Any[
            [4],
            [4, 1],
            [4, 1, 3],
            [4, 2, 3, 1],
            [16, 4, 3, 1, 2],
            [16, 4, 9, 1, 2, 3],
            [16, 4, 10, 1, 2, 3, 9],
            [16, 14, 10, 4, 2, 3, 9, 1],
            [16, 14, 10, 8, 2, 3, 9, 1, 4],
            [16, 14, 10, 8, 7, 3, 9, 1, 4, 2]]
        for i = 1 : length(vs)
            ia = push!(hmax, vs[i])
            @test ia == i
            @test length(hmax) == i
            @test !isempty(hmax)
            @test isequal(list_values(hmax), vs[1:i])
            @test isequal(heap_values(hmax), ss[i])
        end

        # test pop!
        @test isequal(extract_all!(hmax), [16, 14, 10, 9, 8, 7, 4, 3, 2, 1])
        @test isempty(hmax)
    end

    @testset "Custom ordering push! / pop!" begin
        heap = MutableBinaryHeap{Tuple{Int,Int}}(ordering)
        @test length(heap) == 0
        @test isempty(heap)

        # test push!
        ss = Any[
            [(1, 4)],
            [(2, 1), (1, 4)],
            [(2, 1), (1, 4), (3, 3)],
            [(2, 1), (4, 2), (3, 3), (1, 4)],
            [(2, 1), (4, 2), (3, 3), (1, 4), (5, 16)],
            [(2, 1), (4, 2), (3, 3), (1, 4), (5, 16), (6, 9)],
            [(2, 1), (4, 2), (3, 3), (1, 4), (5, 16), (6, 9), (7, 10)],
            [(2, 1), (4, 2), (3, 3), (1, 4), (5, 16), (6, 9), (7, 10), (8, 14)],
            [(2, 1), (4, 2), (3, 3), (1, 4), (5, 16), (6, 9), (7, 10), (8, 14), (9, 8)],
            [(2, 1), (4, 2), (3, 3), (1, 4), (10, 7), (6, 9), (7, 10), (8, 14), (9, 8), (5, 16)]]
        for i = 1 : length(vs2)
            ia = push!(heap, vs2[i])
            @test ia == i
            @test length(heap) == i
            @test !isempty(heap)
            @test isequal(list_values(heap), vs2[1:i])
            @test isequal(heap_values(heap), ss[i])
        end

        # test pop!
        @test isequal(extract_all!(heap), sort(vs2, order=ordering))
        @test isempty(heap)
    end

    @testset "hybrid push! and pop!" begin
        h = MutableBinaryMinHeap{Int}()

        push!(h, 5)
        push!(h, 10)
        @test isequal(heap_values(h), [5, 10])
        @test isequal(list_values(h), [5, 10])

        @test pop!(h) == 5
        @test isequal(heap_values(h), [10])
        @test isequal(list_values(h), [10])

        push!(h, 7)
        push!(h, 2)
        @test isequal(heap_values(h), [2, 10, 7])
        @test isequal(list_values(h), [10, 7, 2])

        @test pop!(h) == 2
        @test isequal(heap_values(h), [7, 10])
        @test isequal(list_values(h), [10, 7])

    end

    @testset "test delete!" begin
        vs=[1,2,4,6,3]
        h = MutableBinaryMinHeap(vs)

        @test isequal(heap_values(delete!(h,1)), [2, 3, 4, 6])
        @test isequal(heap_values(delete!(h,5)), [2, 6, 4])
        @test_throws BoundsError delete!(h,10)
        @test_throws BoundsError delete!(h,0)
        @test_throws BoundsError delete!(h,-5)
        @test pop!(h) == 2
        @test pop!(h) == 4
        push!(h,2)
        @test pop!(h) == 2
        @test pop!(h) == 6
        @test isempty(h)
        #bubble_up test
        vv = [1, 9, 22, 17, 11, 33, 27, 21, 19]
        hp = MutableBinaryMinHeap(vv)
        delete!(hp, 6)
        @test isequal(extract_all!(hp), [1, 9, 11, 17, 19, 21, 22, 27])
    end

    @testset "test delete! at end" begin
        h = MutableBinaryMinHeap{Int}()
        push!(h, 1)
        handle = push!(h, 2)
        delete!(h, handle)
        @test isequal(heap_values(h), [1])
    end

    @testset "test delete! at end of 3-element heap" begin
        h = MutableBinaryMinHeap{Int}()
        push!(h, 1)
        push!(h, 2)
        handle = push!(h, 3)
        delete!(h, handle)
        @test isequal(heap_values(h), [1, 2])
    end

    @testset "test update! and top_with_handle" begin
        for (hf,m) = [(MutableBinaryMinHeap,-2.0), (MutableBinaryMaxHeap,2.0)]
            xs = rand(100)
            h = hf(xs)
            @test length(h) == 100
            @test verify_heap(h)

            for t = 1:1000
                i = rand(1:100)
                v = rand()
                xs[i] = v
                update!(h, i, v)
                @test length(h) == 100
                @test verify_heap(h)
                @test isequal(list_values(h), xs)

                v, i = top_with_handle(h)
                update!(h, i, m)
                xs[i] = m
                @test length(h) == 100
                @test verify_heap(h)
                @test isequal(list_values(h), xs)
                @test top_with_handle(h) == (m, i)

                update!(h, i, -m)
                xs[i] = -m
                @test length(h) == 100
                @test verify_heap(h)
                @test isequal(list_values(h), xs)
                @test top_with_handle(h)[2] ≠ i

                update!(h, i, v)
                xs[i] = v
                @test length(h) == 100
                @test verify_heap(h)
                @test isequal(list_values(h), xs)
                @test top_with_handle(h) == (v, i)

                i = rand(1:100)
                v = rand()
                h[i] = v
                xs[i] = v
                @test length(h) == 100
                @test verify_heap(h)
                @test isequal(list_values(h), xs)
                @test v == h[i]
            end
        end
    end

    @testset "test push! and update! conversion" begin # issue 399
        h = MutableBinaryMinHeap{Float64}()
        push!(h, 3.0)
        push!(h, 5)
        push!(h, Rational(4, 8))
        push!(h, Complex(10.1, 0.0))
        @test isequal(heap_values(h), [0.5, 5.0, 3.0, 10.1])

        update!(h, 2, 20)
        @test isequal(heap_values(h), [0.5, 10.1, 3.0, 20.0])
    end
end # @testset MutableBinheap


text:
#### Precipitation component in atmosphere model
abstract type PrecipitationModel end

export NoPrecipitation, Rain

using ..Microphysics

vars_state(    ::PrecipitationModel, FT) = @vars()
vars_gradient( ::PrecipitationModel, FT) = @vars()
vars_diffusive(::PrecipitationModel, FT) = @vars()
vars_aux(      ::PrecipitationModel, FT) = @vars()

function atmos_nodal_update_aux!(::PrecipitationModel, m::AtmosModel, state::Vars, aux::Vars, t::Real);end
function flux_precipitation!(::PrecipitationModel, atmos::AtmosModel, flux::Grad, state::Vars, aux::Vars, t::Real);end
function diffusive!(::PrecipitationModel, diffusive, ∇transform, state, aux, t, ρD_t);end
function flux_diffusive!(::PrecipitationModel, flux::Grad, state::Vars, diffusive::Vars, aux::Vars, t::Real);end
function flux_nondiffusive!(::PrecipitationModel, flux::Grad, state::Vars, aux::Vars, t::Real);end
function gradvariables!(::PrecipitationModel, transform::Vars, state::Vars, aux::Vars, t::Real);end

"""
    NoPrecipitation <: PrecipitationModel

No precipitation.
"""
struct NoPrecipitation <: PrecipitationModel end

"""
    Rain <: PrecipitationModel

Precipitation model with rain only.
"""
struct Rain <: PrecipitationModel end

vars_state(    ::Rain, FT) = @vars(ρq_rain::FT)
vars_gradient( ::Rain, FT) = @vars(q_rain::FT)
vars_diffusive(::Rain, FT) = @vars(ρd_q_rain::SVector{3,FT})
vars_aux(      ::Rain, FT) = @vars(terminal_velocity::FT, src_q_rai_tot::FT)

function atmos_nodal_update_aux!(rain::Rain, atmos::AtmosModel,
                                         state::Vars, aux::Vars, t::Real)
  FT = eltype(state)
  q_rain = state.precipitation.ρq_rain/state.ρ
  if q_rain > FT(0) # TODO - need a way to prevent negative values
    aux.precipitation.terminal_velocity = terminal_velocity(q_rain, state.ρ)
  else
    aux.precipitation.terminal_velocity = FT(0)
  end

  ρ = state.ρ
  ρinv = 1/state.ρ
  p = aux.pressure
  # TODO - ensure positive definite
  q_tot = max(FT(0), state.ρq_tot*ρinv)
  q_rai = max(FT(0), state.ρq_rain*ρinv)

  # current state
  ts    = thermo_state(rain, state, aux)
  # q     = PhasePartition(q_tot, q_liq, q_ice)
  q     = PhasePartition(ts)
  # T     = air_temperature(e_int, q)
  T     = air_temperature(ts)
  # equilibrium state at current T
  # q_eq = PhasePartition_equil(T, ρ, q_tot)
  q_eq = PhasePartition_equil(ts)

  # tendency from cloud water condensation/evaporation
  # src_q_liq = conv_q_vap_to_q_liq(q_eq, q)# TODO - temporary handling ice

  # tendencies from rain
  src_q_rai_evap = conv_q_rai_to_q_vap(q_rai, q, T , p, ρ)
  src_q_rai_acnv = conv_q_liq_to_q_rai_acnv(q.liq)
  src_q_rai_accr = conv_q_liq_to_q_rai_accr(q.liq, q_rai, ρ)

  aux.precipitation.src_q_rai_tot = src_q_rai_acnv + src_q_rai_accr + src_q_rai_evap
end

function gradvariables!(rain::Rain, transform::Vars, state::Vars, aux::Vars, t::Real)
  transform.precipitation.q_rain = state.precipitation.ρq_rain/state.ρ
end

function diffusive!(rain::Rain, diffusive::Vars, ∇transform::Grad, state::Vars, aux::Vars, t::Real, ρν::Union{Real,AbstractMatrix}, D_T)
  # diffusive flux
  diffusive.precipitation.ρd_q_rain = state.ρ .* (-D_T) .* ∇transform.precipitation.q_rain
end

function flux_precipitation!(rain::Rain, atmos::AtmosModel, flux::Grad, state::Vars, aux::Vars, t::Real)
  u = state.ρu / state.ρ
  flux.precipitation.ρq_rain += state.precipitation.ρq_rain * u
end

function flux_diffusive!(rain::Rain, flux::Grad, state::Vars, diffusive::Vars, aux::Vars, t::Real, D_T)
  flux.precipitation.ρq_rain = diffusive.precipitation.ρd_q_rain
end

function boundarycondition_precipitation!(rain::Rain, stateP::Vars, diffP::Vars, auxP::Vars,
    nM, stateM::Vars, diffM::Vars, auxM::Vars, bctype::BC, t) where {BC}
  stateP.precipitation.ρq_rain = eltype(stateP)(0)
end

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text:
module API

using PCloud
# import ..PCloud: PCloudClient, OPTS, PCloudError
using ..PCloud: HTTP, JSON3

include("pcloud_api.jl")

pushfiles!(body, file) = push!(body, :files => file)
function pushfiles!(body, file::Pair)
    push!(body, :files => HTTP.Multipart(file.first, ioify(file.second)))
end
function pushfiles!(body, file::AbstractString)
    push!(body, :files => open(file, "r"))
end
function pushfiles!(body, files::AbstractVector)
    for file in files
        pushfiles!(body, file)
    end
end

ioify(elem::IO) = elem
ioify(elem) = IOBuffer(elem)
function ioify(elem::IOBuffer)
    seekstart(elem)
    return elem
end

function query(client::PCloud.PCloudClient, method, params)
    response = if haskey(params, :files)
        body = Pair[]
        for (k, v) in params
            k == :files && continue
            # We need strings in multipart. Fortunately all parameters other then
            # `files` should transform to String straightforwardly. But it is of
            # course one of the possible source of errors.
            push!(body, k => string(v))
        end
        pushfiles!(body, params[:files])
    uri = client.apiep * method
        HTTP.post(uri, ["Connection" => "Keep-Alive"], HTTP.Form(body); cookies = client.cookies)
    else
    uri = client.apiep * method * "?" * HTTP.URIs.escapeuri(params)
        HTTP.get(uri, ["Connection" => "Keep-Alive"]; cookies = client.cookies)
    end
    res = JSON3.read(String(response.body))
    if res.result == 0
        return res
    else
        throw(PCloud.PCloudError(res))
    end
end

for (f, fdoc) in PCLOUD_API
    @eval begin
        $f(client::PCloud.PCloudClient = PCloud.OPTS.client; kwargs...) = query(client, String(Symbol($f)), Dict(kwargs))
        @doc $fdoc $f
        export $f
    end
end

end # module


text:
module models
export Conv, Dense, Chain, train!, set_optim!, zero_grad_model

using Knet, Statistics, Random

function zero_grad_model(model)
    for p in params(model)
        zero_grad(p.opt)
    end

function zero_grad(optim::Adam)
    # gradient resetting
    layer1.opt.fstm = Nothing
    layer1.opt.scndm = Nothing

function set_optim!(model, optim, learning_rate)
    for p in params(model)
        p.opt = optim(lr=learning_rate)
    end

function train!(model, L)
    for p in params(model)
        g = grad(L, p)
        update!(value(p), value(g), p.opt)
    end
    return value(L)
end

# Defining the convolutional layer:
struct Conv
    w
    b
    f
    p
end
(c::Conv)(x) = c.f.(pool(conv4(c.w, dropout(x, c.p)) .+ c.b))
Conv(w1::Int, w2::Int, cx::Int, cy::Int, f = relu; pdrop = 0) = Conv(param(w1, w2, cx, cy), param0(1, 1, cy, 1), f, pdrop)

# Defining dense layer
struct Dense
    w
    b
    f
    p
end
(d::Dense)(x) = d.f.(d.w * mat(dropout(x, d.p)) .+ d.b) # mat reshapes 4-D tensor to 2-D matrix so we can use matmul
Dense(i::Int, o::Int, f = relu; pdrop = 0) = Dense(param(o, i), param0(o), f, pdrop)

struct Chain
    layers
    Chain(layers...) = new(layers)
end
(c::Chain)(x) = (for l in c.layers
    x = l(x)
end;
x)
(c::Chain)(x, y) = nll(c(x), y)

end


text:
@testset "CoreModel simple functions" begin
    cp = test_LP()
    @test n_reactions(cp) == 3
    @test n_metabolites(cp) == 4
    @test n_coupling_constraints(cp) == 0

    cp2 = test_LP()
    @test isequal(cp, cp2)
    cp2.S[1] = 1
    @test !isequal(cp, cp2)
    @test isequal(cp, copy(cp))

    cp = test_coupledLP()
    @test n_coupling_constraints(cp) == 2000
    @test isequal(cp, copy(cp))
end

@testset "Analysis utilities" begin
    model_path = download_data_file(
        "http://bigg.ucsd.edu/static/models/e_coli_core.json",
        joinpath("data", "e_coli_core.json"),
        "7bedec10576cfe935b19218dc881f3fb14f890a1871448fc19a9b4ee15b448d8",
    )

    model = load_model(StandardModel, model_path)

    # FBA
    glc = model.reactions["EX_glc__D_e"]
    optimizer = Tulip.Optimizer # quiet by default
    sol = flux_balance_analysis_dict(
        model,
        optimizer;
        modifications = [change_objective("BIOMASS_Ecoli_core_w_GAM")],
    )

    # atom tracker
    atom_fluxes = atom_exchange(sol, model)
    @test isapprox(atom_fluxes["C"], -37.1902, atol = 1e-3)

    # exchange trackers
    consuming, producing = exchange_reactions(sol, model; verbose = false)
    @test isapprox(consuming["EX_nh4_e"], -4.76532, atol = 1e-3)

    # metabolite trackers
    consuming, producing = metabolite_fluxes(sol, model)
    @test isapprox(consuming["atp_c"]["PFK"], -7.47738, atol = 1e-3)
    @test isapprox(producing["atp_c"]["PYK"], 1.75818, atol = 1e-3)

    # set bounds
    cbm = make_optimization_model(model, optimizer)
    ubs = cbm[:ubs]
    lbs = cbm[:lbs]
    glucose_index = first(indexin(["EX_glc__D_e"], reactions(model)))
    o2_index = first(indexin(["EX_o2_e"], reactions(model)))
    atpm_index = first(indexin(["ATPM"], reactions(model)))
    set_bound(glucose_index, cbm; ub = -1.0, lb = -1.0)
    @test normalized_rhs(ubs[glucose_index]) == -1.0
    @test normalized_rhs(lbs[glucose_index]) == 1.0
    set_bound(o2_index, cbm; ub = 1.0, lb = 1.0)
    @test normalized_rhs(ubs[o2_index]) == 1.0
    @test normalized_rhs(lbs[o2_index]) == -1.0
end


text:
# this many digits[bits] (are given by) that many bits[digits]
# sum([abs(i-maximin_digits(maximin_bits(i))) for i=24:4000])

const log10of2 = log10(2)
const log2of10 = log2(10)

@inline maximin_digits(nbits) = floor(Int, log10of2*nbits + 0.125)   # match on round trip
@inline maximin_bits(ndigits) = floor(Int, log2of10*ndigits + 2.625) # at least enough, at most 3 more on round trip

@inline function digits4bits(nbits)
    bits = log10of2 * nbits
    bits += 0.125
    floor(Int, bits)
end

@inline function bits4digits(ndigits)
    digs = log2of10 * ndigits
    digs += 2.625
    floor(Int, digs)
end

#     working_precision exceeds evinced_precision

# additional precision (more of the significant bits)
#    absorbs some kinds of numerical jitter
#    before any undesired resonance occurs
const BitsOfStability  =  9

# additional accuracy (more bits of the significand)
#    may compensate for 1,2,3 ulp enclosure widenings
const BitsOfAbsorption = 15

const ExtraBits = Ref(BitsOfStability + BitsOfAbsorption)

@inline workingbits(evincedbits) = evincedbits + ExtraBits.x
@inline evincedbits(workingbits) = workingbits - ExtraBits.x

# default precision
const MINIMUM_PRECISION = 24
const DEFAULT_PRECISION = Ref(workingbits(128 - ExtraBits.x))

# these typed significands have this many signficant bits

workingprecision(::Type{Mag}) = 30 # bits of significand

workingprecision(::Type{ArbFloat}) = DEFAULT_PRECISION.x
workingprecision(::Type{ArbReal}) = DEFAULT_PRECISION.x
workingprecision(::Type{ArbComplex}) = DEFAULT_PRECISION.x

workingprecision(::Type{ArbFloat{P}}) where {P} = P
workingprecision(::Type{ArbReal{P}}) where {P} = P
workingprecision(::Type{ArbComplex{P}}) where {P} = P

workingprecision(x::ArbFloat{P}) where {P} = P
workingprecision(x::ArbReal{P}) where {P} = P
workingprecision(x::ArbComplex{P}) where {P} = P

# these typed significands have this many signficant bits shown

precision(::Type{Mag}) = 30 # bits of significand
precision(::Type{ArbFloat}) = evincedbits(DEFAULT_PRECISION.x)
precision(::Type{ArbReal}) = evincedbits(DEFAULT_PRECISION.x)
precision(::Type{ArbComplex}) = evincedbits(DEFAULT_PRECISION.x)

precision(::Type{ArbFloat{P}}) where {P} = evincedbits(P)
precision(::Type{ArbReal{P}}) where {P} = evincedbits(P)
precision(::Type{ArbComplex{P}}) where {P} = evincedbits(P)

precision(x::ArbFloat{P}) where {P} = evincedbits(P)
precision(x::ArbReal{P}) where {P} = evincedbits(P)
precision(x::ArbComplex{P}) where {P} = evincedbits(P)

function setprecision(::Type{T}, n::Int) where {T<:Union{ArbFloat,ArbReal,ArbComplex}}
    global DEFAULT_PRECISION
    n < MINIMUM_PRECISION && throw(DomainError("bit precision must be >= $MINIMUM_PRECISION"))
    DEFAULT_PRECISION.x = workingbits(n)
    return n
end

function setworkingprecision(::Type{T}, n::Int) where {T<:Union{ArbFloat,ArbReal,ArbComplex}}
    global DEFAULT_PRECISION
    n < workingbits(MINIMUM_PRECISION) && throw(DomainError("working bit precision must be >= $(workingbits(MINIMUM_PRECISION))"))
    DEFAULT_PRECISION.x = n
    return n
end

extrabits() = ExtraBits.x

function setextrabits(n::Int)
    global ExtraBits, DEFAULT_PRECISION
    ExtraBits.x = max(0,n)
    DEFAULT_PRECISION.x = workingbits(128 - ExtraBits.x)
    return workingprecision(ArbFloat)
end


text:
export evalObj, evalObjAndGrad

function append(A::Tuple,B)
	return (A...,B)
end

function append(A,B)
	return (A,B)
end

function append(A,B::Tuple)
	return(A,B...)
end

function append(A::Tuple,B::Tuple)
	return (A..., B...)
end

function myMap(f::Function,Θ::AbstractArray)
    return f(Θ)
end

function myMap(f::Function,Θ::Tuple)
    fΘ = Array{Any}(undef,length(Θ))
    for k=1:length(Θ)
        fΘ[k] = myMap(f,Θ[k])
    end
    return tuple(fΘ...)
end

function evalObjAndGrad(J,Θ::Vector,parms,ps)
    # put theta into vectors
    parms = vec2param!(Θ,parms)
	# Jc = J(parms)
    # back = gradient(()->J(parms), ps)
    Jc,back = Zygote.pullback(() -> J(parms), ps)
	gc = back(Zygote.sensitivity(Jc))

    dJ = Θ .* 0.0
    cnt = 0;
    for p in ps
		if gc[p]!=nothing
        	gp = vec(gc[p])
			dJ[cnt+1:cnt+length(gp)] = gp
			cnt +=length(gp)
		else
			println("grad was nothing")
		end
    end
    return Jc,dJ
end

function evalObj(J,Θ::Vector,parms,ps)
    # put theta into vectors
    parms = vec2param!(Θ,parms)
    Jc = J(parms)
    return Jc
end

max_stars_repo_name: EmoryMLIP/MFGnet

max_issues_repo_name: EmoryMLIP/MFGnet

max_forks_repo_name: EmoryMLIP/MFGnet


text:
using Condor
using Test

@testset "Condor.jl" begin
    # Write your tests here.
end


text:
function consensus{V}(g::AbstractGraph{V}, m::Vector{Int})
    !is_directed(g) || error("graph must be undirected.")
    @graph_requires g edge_map
    w = zeros(Int, num_edges(g))
    for e in edges(g)
        i = vertex_index(source(e, g), g)
        j = vertex_index(target(e, g), g)
        if m[i] == m[j]
            w[edge_index(e, g)] = 1
        end
    end
    w
end


text:
using SbmlInterface
using Documenter

makedocs(;
    modules=[SbmlInterface],
    authors="Paul F. Lang",
    repo="https://github.com/paulflang/SbmlInterface.jl/blob/{commit}{path}#L{line}",
    sitename="SbmlInterface.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://paulflang.github.io/SbmlInterface.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Installation" => "installation.md",
        "Tutorial" => "tutorial.md",
        "Example" => "examples.md",
        "API reference" => "api_reference.md",
        "Contributing" => "contributing.md",
        "About" => "about.md"
    ],
)

deploydocs(;
    repo="github.com/paulflang/SbmlInterface.jl",
    devbranch="main"
)


text:
using EngEconomics, Roots

# Given
aInit = -200000
aAnnual = 30000
aScrap = 45000
aLife = 20

bInit = -180000
bAnnual = 24000
bScrap = 20000
bLife = 15

MARR = 0.09

# Find Annual Worth Analysis
aAW = aInit * capitalRecoveryFactor(MARR, aLife) + aAnnual + aScrap * sinkingFundFactor(MARR, aLife)

bAW = bInit * capitalRecoveryFactor(MARR, bLife) + bAnnual + bScrap * sinkingFundFactor(MARR, bLife)

if aAW > bAW
	println("A")
else
	println("B")
end


text:
using RDatasets, DataArrays, Gadfly

set_default_plot_size(6inch, 3inch)

plot(dataset("car", "SLID"), x="Wages", color="Language",
     Geom.histogram(position=:stack))


text:
using Base.Test
using Tensors

# Single mode QR

t = rand([Mode(k,2) for k in 1:4])
q,r = qr(t, 3)
@test t ≈ q*r

# Multiple modes QR

t = rand([Mode(k,2) for k in 1:4])
q,r = qr(t, Any[3,1], 5)
@test t ≈ q*r

# Single mode SVD

t = rand([Mode(k,2) for k in 1:4])
u,s,v = svd(t, 3)
@test t ≈ u*scale(s,v)

# Multiple modes SVD

t = rand([Mode(k,2) for k in 1:4])
u,s,v = svd(t, Any[4,2], 5)
@test t ≈ scale(u,s)*v


text:
module Day11

export get_inputs, get_solution1, get_solution2

using OffsetArrays

## Input getting
function get_inputs()
    test_input1 = test_input2 = read_inputs("test_input1.txt")
    test_output1 = 37
    test_output2 = 26
    data = read_inputs("input.txt")
    return (; test_input1, test_input2, test_output1, test_output2, data)
end

function read_inputs(f_name)
    lines = collect.(readlines(joinpath(@__DIR__, f_name)))
    mat = permutedims(reduce(hcat, lines), [2,1])
    return pad_matrix(mat)
end

## Solution functions
function pad_matrix(M, default='.')
    P = OffsetArray(similar(M, (size(M).+2)...), -1, -1)
    P .= default
    for i in CartesianIndices(M)
        P[i] = M[i]
    end
    return P
end

const adj_inds = [(i, j) for i in -1:1, j in -1:1 if (i!=0 || j!=0)]

get_adjacent(data, ci_tuple, sz) = (data[(ij .+ ci_tuple)...] for ij in adj_inds)

get_seen(data, ci_tuple, sz) = Base.Generator(adj_inds) do direction
    idx = ci_tuple
    elem = '.'
    while all(idx .∈ sz) && elem=='.'
        idx = idx .+ direction
        elem = data[idx...]
    end
    elem
end

function update_seats!(data, temp=copy(data), see_func=get_adjacent, occ_tolerance=4)
    temp .= data
    sz = ntuple(n -> 1:lastindex(data, n)-1, 2)
    @inbounds for ci in CartesianIndices(sz)
        ci_tuple = Tuple(ci)
        seen = see_func(data, ci_tuple, sz)
        occupied = (x == '#' for x in seen)
        if data[ci] == 'L'
            if !any(occupied)
                temp[ci] = '#'
            end
        elseif data[ci] == '#'
            if count(occupied) >= occ_tolerance
                temp[ci] = 'L'
            end
        end
    end
    data .= temp
    return data
end

function get_solution(data, see_func=get_adjacent, occ_tolerance=4)
    data = copy(data)
    temp = similar(data)
    new_data = update_seats!(copy(data), temp, see_func, occ_tolerance)
    while new_data != data
        data .= new_data
        update_seats!(new_data, temp, see_func, occ_tolerance)
    end
    return count(x == '#' for x in new_data)
end

# Part 1
get_solution1(data) = get_solution(data)

# Part 2
get_solution2(data) = get_solution(data, get_seen, 5)

end


text:
abstract type AbstractXCCalculator
end

# Use internal XC
struct XCCalculator <: AbstractXCCalculator
end

# Use Libxc
struct LibxcXCCalculator <: AbstractXCCalculator
end

# These types should contain xcfunc field


text:
export SlicedSystem, slice

"""
    SlicedSystem(F::AbstractSystem, L::LinearSubspace)

Given a polynomial system `F` with the affine linear subspace `L`
defined by ``Ax = b`` this constructs the system `[F(x); A x - b] = 0`.
See also [`slice`](@ref).
"""
struct SlicedSystem{S<:AbstractSystem,T} <: AbstractSystem
    system::S
    affine::LinearSubspace{T}
end
SlicedSystem(
    F::System,
    A::LinearSubspace;
    compile::Union{Bool,Symbol} = COMPILE_DEFAULT[],
) = SlicedSystem(fixed(F; compile = compile), A)

"""
    slice(F::System, L::LinearSubspace; compile = $(COMPILE_DEFAULT[]))
    slice(F::AbstractSystem, L::LinearSubspace)

Slice the zero set of the polynomial system `F` with the affine linear subspace `L`
defined by ``Ax = b``. This constructs the system `[F(x); A x - b] = 0`.
See also [`SlicedSystem`](@ref).
"""
slice(F::System, A::LinearSubspace; compile::Union{Bool,Symbol} = COMPILE_DEFAULT[]) =
    SlicedSystem(fixed(F; compile = compile), A)
slice(
    F::AbstractSystem,
    A::LinearSubspace;
    compile::Union{Bool,Symbol} = COMPILE_DEFAULT[],
) = SlicedSystem(F, A)
Base.intersect(F::Union{System,AbstractSystem}, A::LinearSubspace) = SlicedSystem(F, A)
Base.intersect(A::LinearSubspace, F::Union{System,AbstractSystem}) = SlicedSystem(F, A)

ModelKit.variables(F::SlicedSystem) = variables(F.system)
ModelKit.parameters(F::SlicedSystem) = parameters(F.system)
ModelKit.variable_groups(F::SlicedSystem) = nothing

linear_subspace(A::SlicedSystem) = A.affine
system(A::SlicedSystem) = A.system

function Base.size(F::SlicedSystem)
    m, n = size(F.system)
    (m + codim(F.affine), n)
end

function (F::SlicedSystem)(x::AbstractVector, p = nothing)
    u = F.system(x, p)
    append!(u, extrinsic(F.affine)(x))
    u
end

function ModelKit.evaluate!(u, F::SlicedSystem, x, p = nothing)
    evaluate!(u, F.system, x, p)

    @unpack A, b = extrinsic(F.affine)
    m = size(F.system, 1)
    n, N = size(A)
    @inbounds for i = 1:n
        u[m+i] = -b[i]
    end
    @inbounds for j = 1:N
        xj = x[j]
        for i = 1:n
            u[m+i] += A[i, j] * xj
        end
    end

    u
end

function ModelKit.evaluate_and_jacobian!(u, U, F::SlicedSystem, x, p = nothing)
    evaluate_and_jacobian!(u, U, F.system, x, p)

    @unpack A, b = extrinsic(F.affine)
    m = size(F.system, 1)
    n, N = size(A)
    @inbounds for i = 1:n
        u[m+i] = -b[i]
    end
    @inbounds for j = 1:N
        xj = x[j]
        for i = 1:n
            u[m+i] += A[i, j] * xj
        end
    end

    @inbounds for j = 1:N, i = 1:n
        U[m+i, j] = A[i, j]
    end

    nothing
end

function ModelKit.taylor!(u, v::Val, F::SlicedSystem, tx, p = nothing)
    u .= zero(eltype(u))
    taylor!(u, v, F.system, tx, p)
    u
end


text:
using Test

@testset "multithreading_enabled" begin
    @test Threads.nthreads() == 8
end


text:
using Revise
include("../src/FxRobot.jl");import .FxRobot as fx

asset = "ETHUSDT"
from = "2021-03-01"
to = "2021-03-10"
period = "1m"
out_path = "./data/csv/daily"

fx.download_daily_binance_data(asset,from,to,period,out_path)

asset = "ETHUSDT"
from = "2021-01-01"
# to = string(today())
to = "2021-11-01"
period = "1m"
out_path = "./data/csv/monthly"

fx.download_monthly_binance_data(asset,from,to,period,out_path)

file_regexp = r"(.*)(ETHUSDT-1m)(.*)(\.csv)"
df = fx.trans_df( fx.csv_to_dataframe(out_path;file_regexp=file_regexp) )

arrow_path = abspath("./data/arrow/ETHUSDT-1m-2021.arrow")
fx.arrow_write(df,arrow_path)
df = fx.arrow_read(arrow_path)
df

max_stars_repo_name: mrchaos/fxrobot

max_issues_repo_name: mrchaos/fxrobot

max_forks_repo_name: mrchaos/fxrobot


text:
using ReachabilityAnalysis, SparseArrays, JLD2

examples_dir = normpath(@__DIR__, "..", "..", "..", "examples")
HEAT01_path = joinpath(examples_dir, "Heat3D", "HEAT01.jld2")


text:
function parity(s::Int, i::Int, j::Int)
    c = 0
    for k=min(i, j)+1:max(i, j)-1
        c += (s>>k)&1
    end
    return isodd(c) ? -1. : 1.
end