text:
using Test, DataFrames, Suppressor, HMCResearchRandomWalks.Enz

@testset "Setter for simulation control" begin
    set_THICK_ENZ()
    @test Enz.THICK_ENZ == true
    set_THICK_ENZ(false)
    @test Enz.THICK_ENZ == false
    set_THICK_ENZ(true)
    @test Enz.THICK_ENZ == true
end

@testset "fast custom enz_sim" begin
    set_THICK_ENZ(false)
    set_NUMBER_OF_WALKS(10)
    set_MAX_STEPS_PER_WALK(1000000)
    set_ENZ_STEP_SIZE(0.07)
    output = @capture_out begin enz_sim(5534) end
    correct_output = string(
    "############################\n",
    "       Enz Experiment       \n",
    "############################\n",
    "_________Parameters_________\n",
    "Enzyme Thickness:\t0.2\n",
    "# of trials:\t\t10\n",
    "# of steps (max):\t1000000\n",
    "Step size, water:\t0.1\n",
    "Step size, enz:\t\t0.07\n",
    "Random seed:\t\t5534\n",
    "_________Results____________\n",
    "# in sensor:\t\t10\n",
    "# of escaped:\t\t0\n",
    "# unresolved:\t\t0\n"
    )
    @test output === correct_output
end

@testset "enz_print" begin
    n = 12312020
    m = 1225
    seed = 4321
    set_NUMBER_OF_WALKS(n)
    set_MAX_STEPS_PER_WALK(m)
    output = @capture_out begin Enz.enz_print(seed) end
    correct_output = string(
    "############################\n",
    "   Compare Enz Step Sizes   \n",
    "############################\n",
    "_________Parameters_________\n",
    "# of trials:\t\t$n\n",
    "# of steps (max):\t$m\n",
    "Step size, water:\t0.1\n",
    "Random seed:\t\t4321\n"
    )
    @test output === correct_output
end

@testset "getdata_enz.." begin
    set_NUMBER_OF_WALKS(5)
    set_MAX_STEPS_PER_WALK()
    @suppress begin
        global getdata_enz_test_df = getdata_enzstepsize(
            [0.003, 0.005, 0.007],
            1234
        )
    end
    thk = "thick, 2 um"
    thn = "thin, 0.2 um"
    correct_df = DataFrame(
        thick_or_thin = [thk, thn, thk, thn, thk, thn],
        enz_step_size = [0.003, 0.003, 0.005, 0.005, 0.007, 0.007],
        sensor_yield = [1,5,3,5,2,5],
        escaped = [4,0,2,0,3,0],
        unresolved = [0,0,0,0,0,0]
    )
    @test getdata_enz_test_df == correct_df
end

max_stars_repo_name: s-y-wu/HMCResearchRandomWalks

max_issues_repo_name: s-y-wu/HMCResearchRandomWalks

max_forks_repo_name: s-y-wu/HMCResearchRandomWalks


text:
# Autogenerated wrapper script for FLINT_jll for i686-linux-musl
export libflint

using GMP_jll
using MPFR_jll
JLLWrappers.@generate_wrapper_header("FLINT")
JLLWrappers.@declare_library_product(libflint, "libflint.so.16")
function __init__()
    JLLWrappers.@generate_init_header(GMP_jll, MPFR_jll)
    JLLWrappers.@init_library_product(
        libflint,
        "lib/libflint.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
      if !Sys.iswindows() && !(get(ENV, "NEMO_THREADED", "") == "1")
    #to match the global gmp ones
    fm = dlsym(libflint_handle, :__flint_set_memory_functions)
    ccall(fm, Nothing,
      (Ptr{Nothing},Ptr{Nothing},Ptr{Nothing},Ptr{Nothing}),
        cglobal(:jl_malloc),
        cglobal(:jl_calloc),
        cglobal(:jl_realloc),
        cglobal(:jl_free))
  end

end  # __init__()

ext: jl


text:
#=
   From Eugene Charniak
   "Bayesian Networks without Tears"
   page 51

   """
   Suppose when I go home at night, I want to know if my family is home before I try the doors.
   (Perhaps the most convenient door to enter is double locked when nobody is home.)
   Now, often when my wife leaves the house, she turns on an outdoor light. However,
   she sometimes turns on this light if she is expecting a guest. Also, we have a dog.
   When nobody is home, the dog is put in the back yard. The same is true if the dog
   has bowel troubles. Finally, if the dog is in the backyard, I will probably hear her barking (or
   what I think is her barking), but sometimes I can be confused by other dogs barking. This
   example, partially inspired by Pearl’s (1988) earthquake example, is illustrated in figure 1.
   There we find a graph not unlike many we see in AI. We might want to use such diagrams to
   predict what will happen (if my family goes out, the dog goes out) or to infer causes from
   observed effects (if the light is on and the dog is out, then my family is probably out).

   The important thing to note about this example is that the causal connections are
   not absolute. Often, my family will have left without putting out the dog or turning on a
   light. Sometimes we can use these diagrams anyway, but in such cases, it is hard to know
   what to infer when not all the evidence points the same way. Should I assume the family is
   out if the light is on, but I do not hear the dog? What if I hear the dog, but the light is
   out? Naturally, if we knew the relevant probabilities, such as P(family-out | light-on, ¬ hear-
   bark), then we would be all set. However, typically, such numbers are not available for
   all possible combinations of circumstances. Bayesian networks allow us to calculate them
   from a small set of probabilities, relating only neighboring nodes.
   """

  Summary Statistics
     parameters      mean       std   naive_se      mcse          ess      rhat   ess_per_sec
         Symbol   Float64   Float64    Float64   Float64      Float64   Float64       Float64

     family_out    0.5434    0.4981     0.0050    0.0238     412.4049    1.0003       43.8169
  bowel_problem    0.0045    0.0669     0.0007    0.0027     617.6361    1.0026       65.6222
       light_on    0.9996    0.0200     0.0002    0.0004    2332.5805    1.0003      247.8305
        dog_out    0.4205    0.4937     0.0049    0.0206     564.5462    1.0002       59.9815
      hear_bark    0.0001    0.0100     0.0001    0.0001   10004.0032    1.0000     1062.8988

   Cf ~/webppl/family_out_problem.wppl
=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function family_out_problem()
    family_out    ~ flip(0.15)
    bowel_problem ~ flip(0.01)
    light_on      ~ family_out ? flip(0.6) : flip(0.05)
    dog_out       ~
        (family_out==true && bowel_problem==true)   ? flip(0.99) :
        (family_out==true && bowel_problem==false)  ? flip(0.90) :
        (family_out==false && bowel_problem==true)  ? flip(0.97) :
        (family_out==false && bowel_problem==false) ? flip(0.3)  : Dirac(false)

    hear_bark ~ dog_out ? flip(0.7) : flip(0.01);

    #=
      To take the earlier example, if I observe that the light is on (light-on = true) but
      do not hear my dog (hear-bark = false), I can calculate the conditional probability
      of family-out given these pieces of evidence. (For this case, it is .5.)
    =#
    true ~ Dirac(light_on == true)
    true ~ Dirac(hear_bark == false)

end

model = family_out_problem()

num_chains = 4

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(5), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns)
# display(plot(chns))

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
# replace underscores with spaces to LaTeX is happy
replace_(x) = replace(String(x), "_" => " ")
sumby(v, n) = [sum(@view(v[i:min(i+n-1, length(v))])) for i in 1:n:length(v)]
ssum(x, s) = sum(sum.(x[s]))

# Main workhorse function for line plots
function make_plot(
        data::SocketCounterRecord,
        names;
        sumover = 1,
        plotsum = false,
        xlabel = "Time (s)",
        ylabel = "",
        ymax = nothing,
        modifier = identity,
        title = "",
        reducer = sum,
        linewidth = "4pt",
    )

    selected_data = Dict(k => reducer.(retrieve(data, k)) for k in names)

    if plotsum
        selected_data[:total_sum] = map(sum, zip(values(selected_data)...))
        push!(names, :total_sum)
    end
    _sums = Dict(k => sum(v) for (k,v) in selected_data)
    display(_sums)

    # Plot the results
    plots = []
    for name in names
        v = selected_data[name]
        coords = Coordinates(1:sumover:length(v), sumby(modifier(v), sumover))
        push!(plots,
            @pgf(PlotInc(
                {
                    line_width = linewidth,
                },
                coords
            )),
            @pgf(LegendEntry(replace_(name)))
        )
    end

    empty!(PGFPlotsX.CUSTOM_PREAMBLE)
    push!(PGFPlotsX.CUSTOM_PREAMBLE, "\\usepgfplotslibrary{colorbrewer}")
    plt = TikzDocument()
    push!(plt, "\\pgfplotsset{cycle list/Paired}")

    tikz = TikzPicture()

    opts = []
    isnothing(ymax) || push!(opts, "ymax = $ymax")

    axs = @pgf Axis(
        {
            width = "8cm",
            height = "8cm",
            ultra_thick,
            legend_style = {
                at = Coordinate(0.02, 1.02),
                anchor = "south west",
            },
            legend_columns = 2,
            grid = "both",
            xlabel = xlabel,
            ylabel = ylabel,
            title = title,
            opts...
        },
        plots...
    )
    push!(tikz, axs)
    push!(plt, tikz)

    return plt
end

labelize(x) = x
labelize(x::String) = replace_(x)
labelize(x::Tuple) = join(labelize.(x), " ")

# Now bar plots!
function bar_plot(
        # The actual data for each item of interest.
        data::Vector{<:SocketCounterRecord},
        # The names of the counters to retrieve,
        counters,
        # The series lables
        labels;
        xlabel = "",
        ylabel = "",
        title = "",
        aggregator = sum,
        reducer = sum,
        ymax = nothing,
        width = "8cm",
        height = "8cm",
    )

    @assert length(data) == length(labels)

    # First step - pull out the data we want into dictionary form.
    selected_data = map(data) do d
        # For each counters, aggregate the counters according to `aggregate`,
        # then reduce the resulting vector to a single number using `reducer`.
        return map(counters) do counter
            aggregated = aggregator.(retrieve(d, counter))
            @show sum(aggregated)
            reduced_data = reducer(aggregated)
            return counter => reduced_data
        end |> Dict
    end

    # Plot the results
    plots = []
    for counter in counters
        y_coords = getindex.(selected_data, counter)
        x_coords = labelize.(labels)

        coords = Coordinates(x_coords, y_coords)
        push!(plots,
            @pgf(PlotInc(coords)),
            @pgf(LegendEntry(replace_(counter)))
        )
    end

    empty!(PGFPlotsX.CUSTOM_PREAMBLE)
    push!(PGFPlotsX.CUSTOM_PREAMBLE, "\\usepgfplotslibrary{colorbrewer}")
    plt = TikzDocument()
    push!(plt, "\\pgfplotsset{cycle list/Paired}")

    tikz = TikzPicture()

    opts = []
    isnothing(ymax) || push!(opts, "ymax = $ymax")

    axs = @pgf Axis(
        {
            width = width,
            height = height,
            ultra_thick,
            ybar,
            ymin = 0,
            legend_style = {
                at = Coordinate(0.98, 1.02),
                anchor = "south east",
            },
            legend_columns = 2,
            grid = "both",

            # Set up symbolic coordinates.
            symbolic_x_coords = labelize.(labels),
            nodes_near_coords_align={vertical},
            xtick = "data",
            xticklabel_style = {
                rotate = -15,
            },
            #enlarge_x_limits = 0.3,

            # Setup labels
            xlabel = xlabel,
            ylabel = ylabel,
            title = title,
            opts...
        },
        plots...
    )
    push!(tikz, axs)
    push!(plt, tikz)

    return plt
end

# Look in the `data` directory - deserialize all files that match `prefix`.
# Furthermore - do a pairwise merging of the named tuple entries as long as the difference
# in lengths between the everything is with `max_truncate`.
function load(
        file,
        params = NamedTuple();
        # keyword arguments
        max_truncate = 5,
        socket = 2,
    )

    database = deserialize(joinpath(DATADIR, file))

    # Filter out entries that match the request.
    filtered_database = database[params]
    if length(filtered_database) == 0
        println(database)
        throw(error("No results matching your query!"))
    elseif length(filtered_database) > 1
        println(filtered_database)
        throw(error("Found $(length(filtered_database)) entries"))
    end

    # Get the actual payload from the data.
    data = first(filtered_database).data[Symbol("socket_$(socket-1)")]::SocketCounterRecord

    # Check if number of samples is about the same.
    min, max = extrema(length, walkleaves(data))
    if max - min > max_truncate
        error("Sizes of data not within $max_truncate. Sizes are between $min and $max.")
    end

    # Funtion to resize all arrays
    resize!.(walkleaves(data), min)
    return data
end

function load_multiple(
        file,
        params::Union{<:NamedTuple, <:Vector},
        # The field that should be different across the selected matches
        should_vary;
        socket = 2,
    )

    db = deserialize(joinpath(DATADIR, file))

    # Filter out entries that match the request.
    #
    # If `params` is a NamedTuple, find all matches.
    # If `params` is a vector of named tuples, find a match for each.
    if isa(params, NamedTuple)
        filtered_database = db[params]
    elseif isa(params, Vector)
        filtered_database = Slab()
        for p in params
            @show p
            x = db[p]
            if length(x) > 1
                println(x)
                error()
            end
            filtered_database[p] = x[1].data
        end
    end

    if isnothing(filtered_database)
        println(db)
        throw(error("No results matching your query!"))
    end

    # Ensure that the `should_vary` field of the filtered_database does infact have
    # all unique entries.
    @assert allunique(filtered_database[should_vary])

    # Now, extract the socket records and return!
    key = Symbol("socket_$(socket - 1)")
    data = [entry.data[key]::SocketCounterRecord for entry in filtered_database]
    return data, filtered_database[should_vary]
end

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
"""
    Comrade
Radio Observation Sampling Exploration
"""
module Comrade

using AbstractFFTs
using AbstractMCMC
using Accessors: @set
using ArgCheck: @argcheck
using BasicInterpolators
using DocStringExtensions
using ChainRulesCore
using ComradeBase
using FITSIO
using FileIO
using ForwardDiff
using FFTW: fft, fftfreq, fftshift, ifft, ifft!, ifftshift, plan_fft
using FLoops
#using MappedArrays: mappedarray
using MeasureBase, MeasureTheory
using NFFT: nfft, plan_nfft
using PaddedViews
using PyCall: pyimport, PyNULL, PyObject
using SpecialFunctions
using Reexport
using Requires: @require
using StructArrays: StructVector, StructArray, append!!
using Tables
using UUIDs
# Write your package code here.

@reexport using ComradeBase
export AD
using ComradeBase: visanalytic, imanalytic

export SequentialEx, ThreadedEx

const ehtim = PyNULL()

"""
    $(SIGNATURES)
Loads the [eht-imaging](https://github.com/achael/eht-imaging) library and stores it in the
`ehtim` variable.

# Notes
This will fail if ehtim isn't installed in the python installation that PyCall references.
"""
function load_ehtim()
    #try
    copy!(ehtim, pyimport("ehtim"))
    #catch
    #    @warn "No ehtim installation found in python path. Some data functionality will not work"
    #end
end

export rad2μas, μas2rad, ehtim, load_ehtim, logdensity_def, logdensityof

"""
    rad2μas(x)
Converts a number from radians to μas
"""
@inline rad2μas(x) = 180/π*3600*1e6*x

"""
    μas2rad(x)
Converts a number from μas to rad
"""
@inline μas2rad(x) = x/(180/π*3600*1e6)

#include("interface.jl")
#include("images/images.jl")
import ComradeBase: flux, radialextent, intensitymap, intensitymap!
export create_cache
include("observations/observations.jl")
include("models/models.jl")
include("distributions/distributions.jl")
include("distributions/radiolikelihood.jl")
include("visualizations/visualizations.jl")
include("bayes/bayes.jl")
include("inference/inference.jl")
include("calibration/gains.jl")

function __init__()
    # @require AdvancedHMC="0bf59076-c3b1-5ca4-86bd-e02cd72cde3d" include("inference/advancedhmc.jl")
    # @require AdaptiveMCMC="717c3277-546c-407d-8270-09a39a0919a0" include("inference/adaptivemcmc.jl")
    # @require NestedSamplers="41ceaf6f-1696-4a54-9b49-2e7a9ec3782e" include("inference/nested.jl")
    # @require Dynesty="eb527566-0f3e-4aab-bb5f-9d2e403dba70" include("inference/dynesty.jl")
    # @require Pathfinder="b1d3bc72-d0e7-4279-b92f-7fa5d6d2d454" include("inference/pathfinder.jl")
    # @require GalacticOptim="a75be94c-b780-496d-a8a9-0878b188d577" include("inference/galacticoptim.jl")
    # FIX THIS
    del_format(format"FITS")
    add_format(format"FITS",
        # See https://www.loc.gov/preservation/digital/formats/fdd/fdd000317.shtml#sign
        [0x53,0x49,0x4d,0x50,0x4c,0x45,0x20,0x20,0x3d,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x20,0x54],
        [".fit", ".fits", ".fts", ".FIT", ".FITS", ".FTS", ".fit",],
        [:FITSIO => UUID("525bcba6-941b-5504-bd06-fd0dc1a4d2eb")],
        [:AstroImages => UUID("fe3fc30c-9b16-11e9-1c73-17dabf39f4ad")],
        [:Comrade => UUID("99d987ce-9a1e-4df8-bc0b-1ea019aa547b")]
    )

end

end


text:
module EmpiricalCovarianceOperatorsTests
    using EmpiricalCovarianceOperators
    using Test
    using Random
    using Statistics

    @testset "CovarianceOperator Tests" begin
        x = randn(500)
        X = randn(20,500)
        C = cov(X,X)
        A = CovarianceOperator(X)
        y = C*x
        Y = randn(20,300)
        A2 = CrossCovarianceOperator(X,Y)
        C2 = cov(X,Y)
        x2 = randn(300)
        y2 = C2*x2
        y1p5 = C2'*x
        @test isapprox(C, convert(Array{Float64}, A))
        @test isapprox(y, A*x)
        @test isapprox(convert(Array{Float64}, A), convert(Array{Float64}, A'))
        @test_broken isapprox(C2, convert(Array{Float64}, A2))
        @test isapprox(y2, A2*x2)
        @test isapprox(y1p5, A2'*x)
        @test_broken isapprox(x'*B*x2, x2'*B'*x)
    end
end


text:
# Required for BSON loading output.bson
using Cytof5, Distributions

using Random

# TODO: Remove this dependency
using RCall

# TODO: remove these in favor of BSON
# using JLD2, FileIO
using BSON

include("util.jl")
#include("CytofImg.jl")

"""
# Example

```julia
post_process("path/to/output.bson")
```
"""
function post_process(PATH_TO_OUTPUT, thresh=0.99, min_presences=[0, .01, .03, .05]) # bla/output.bson
  OUTDIR = join(split(PATH_TO_OUTPUT, "/")[1:end-1], "/")
  IMGDIR = "$OUTDIR/img/"
  mkpath(IMGDIR)

  println("Loading Data ...")
  BSON.@load "$PATH_TO_OUTPUT" out ll lastState c metrics init dden simdat
  y_dat = simdat[:y]

  I, K = size(simdat[:W])
  K_MCMC = size(lastState.W, 2)
  J = size(lastState.Z, 1)
  MCMC_ITER = length(out[1])
  BURN = length(ll) - MCMC_ITER

  # Plot loglikelihood
  util.plotPdf("$(IMGDIR)/ll.pdf")
  util.plot(ll[(end-MCMC_ITER):end], ylab="log-likelihood", xlab="MCMC iteration (post-burn)",
            typ="l");
  util.devOff()

  util.plotPdf("$(IMGDIR)/ll_entire_history.pdf")
  util.plot(ll, ylab="log-likelihood", xlab="MCMC iteration", typ="l");
  util.abline(v=BURN);
  util.devOff()

  function addGridLines(J::Int, K::Int, col="grey")
    util.abline(v=(1:K) .+ .5, h=(1:J) .+ .5, col=col);
  end

  # Plot Z
  Zpost = util.getPosterior(:Z, out[1])
  Zmean = mean(Zpost)

  util.plotPdf("$IMGDIR/Z_mean.pdf")
  util.myImage(Zmean, xlab="Features", ylab="Markers", addL=true, col=util.greys(11),
               f=Z->addGridLines(J,K_MCMC));
  util.devOff()

  util.plotPdf("$IMGDIR/Z_mean_est_leftordered.pdf")
  util.myImage(Cytof5.Model.leftOrder((Zmean .> .5)*1),
               xlab="Features", ylab="Markers", addL=true, col=util.greys(11),
               f=Z->addGridLines(J,K_MCMC));
  util.devOff()

  util.plotPdf("$IMGDIR/Z_true.pdf")
  util.myImage(simdat[:Z], xlab="Features", ylab="Markers");
  addGridLines(J, K)
  util.devOff()

  # Plot W
  Wpost = util.getPosterior(:W, out[1])
  Wpost = cat(Wpost..., dims=3) # I x K_MCMC x NMCMC

  util.plotPdf("$IMGDIR/W.pdf")
  for i in 1:I
    util.boxplot(Wpost[i, :, :]', xlab="Features", ylab="W$(i)", );
    # Plot data
    util.abline(h=simdat[:W][i, :], lty=2, col="grey")
  end
  util.devOff()

  # Plot alpha
  alphaPost = util.getPosterior(:alpha, out[1])
  println("Making alpha ...")
  util.plotPdf("$(IMGDIR)/alpha.pdf")
  util.plotPost(alphaPost, ylab="density", xlab="alpha", main="");
  util.devOff()

  # Posterior of eps
  println("Making eps...")
  epsPost = hcat(util.getPosterior(:eps, out[1])...)'
  util.plotPdf("$IMGDIR/eps.pdf")
  util.plotPosts(epsPost);
  util.devOff()

  # Get lam
  lamPost = util.getPosterior(:lam, out[1])
  # unique(lamPost)

  # Get mus
  delta0Post = cat([m[:delta][0] for m in out[1]]..., dims=2) # L0 x NMCMC
  delta1Post = cat([m[:delta][1] for m in out[1]]..., dims=2) # L1 x NMCMC
  mus0Post = -cumsum(delta0Post, dims=1)
  mus1Post = cumsum(delta1Post, dims=1)
  musPost = [ mus0Post; mus1Post ]

  util.plotPdf("$IMGDIR/mus.pdf")
  R"boxplot"(musPost', ylab="mu*", xlab="", xaxt="n", col="steelblue", pch=20, cex=0);
  #util.plot(1:size(musPost, 2), mean(musPost, dims=1), typ="n", ylab="μ*", xlab="", xaxt="n")
  #util.addErrbar(R"t(apply($musPost, 2, quantile, c(.025, .975)))",
  #               x=1:size(musPost, 2), ylab="μ*", xlab="", xaxt="n", col="blue", lend=1, lwd=10);
  util.abline(h=0, v=size(mus0Post, 2) + .5, col="grey30", lty=1);
  util.abline(h=simdat[:mus][0], lty=2, col="steelblue");
  util.abline(h=simdat[:mus][1], lty=2, col="steelblue");
  util.devOff()

  # Get sig2
  sig2Post = hcat(util.getPosterior(:sig2, out[1])...) # I x NMCMC

  util.plotPdf("$IMGDIR/sig2.pdf")
  util.boxplot(sig2Post', xlab="samples", ylab="sig2")
  util.abline(h=simdat[:sig2], lty=2, col="grey")
  util.devOff()

  # Plot posterior predictive of expressed cells (blue)
  # Compared against observed positive expressions, with kde bw=0.1
  function plot_dden_expressed(; y_dat, out1, bw=.1, ygridsize=200, ygrid_range=[-1, 6])
    mkpath("$(IMGDIR)/dden-expressed")
    # Get eta post
    eta1Post = cat([m[:eta][1] for m in out1]..., dims=4) # I x J x L1 x NMCMC
    ygrid = collect(range(ygrid_range[1], stop=ygrid_range[2], length=ygridsize))
    ygrid = reshape(ygrid, 1, ygridsize)
    for i in 1:I
      for j in 1:J
        print("\ri: $i, j: $j   ")
        ddij = [begin
                  eta1_ij = eta1Post[i, j, :, b:b]
                  comp2 = pdf.(Normal.(mus1Post[:, b:b], sqrt(sig2Post[i, b])), ygrid)
                  vec(sum(eta1_ij .* comp2, dims=1))
                end for b in 1:MCMC_ITER]
        ddij = hcat(ddij...) # ygridsize x NMCMC
        ddij_mean = mean(ddij, dims=2)
        ddij_lower = [quantile(ddij[g, :], .025) for g in 1:ygridsize]
        ddij_upper = [quantile(ddij[g, :], .975) for g in 1:ygridsize]

        util.plotPdf("$(IMGDIR)/dden-expressed/dden_i$(i)_j$(j).pdf")
        util.plot(ygrid, ddij_mean, lwd=2, col="blue",
                  typ="l", xlab="y", ylab="density", main="",
                  xlim=[-1, 6], ylim=[0, maximum(ddij_upper)])
        util.colorBtwn(ygrid, ddij_lower, ddij_upper, from=-10, to=10,
                       col=util.rgba("blue", .3))
        yij = y_dat[i][:, j]
        yij_pos = yij[yij .> 0]
        util.lines(util.density(yij_pos, bw=bw))
        util.devOff()
      end
    end
    println()
  end
  plot_dden_expressed(y_dat=simdat[:y], out1=out[1])

  # Posterior of y_imputed
  y_imputed = [ o[:y_imputed] for o in out[2] ]
  util.plotPdf("$(IMGDIR)/ydatPost.pdf")
  R"par(mfrow=c(4,2))"
  for i in 1:I
    for j in 1:J
      util.plot(util.density(simdat[:y][i][:, j], na=true), col="red", xlim=[-8,8],
                main="Y sample: $(i), marker: $(j)", bty="n", fg="grey")
      for iter in 1:length(y_imputed)
        yimp = y_imputed[iter]
        util.lines(util.density(yimp[i][:, j]), col=util.rgba("blue", .5))
      end
      util.lines(util.density(simdat[:y_complete][i][:, j]), col="grey")
    end
  end
  R"par(mfrow=c(1,1))"
  util.devOff()

  idx_missing = [ findall(isnan.(y_dat[i])) for i in 1:I ]
  idx = idx_missing[2][1]
  util.plotPdf("$(IMGDIR)/y_trace.pdf")
  util.hist([ y_imputed[b][2][idx] for b in 1:length(y_imputed) ], col="blue", border="transparent")
  util.plot([ y_imputed[b][2][idx] for b in 1:length(y_imputed) ], typ="l")
  util.devOff()

  # ARI - adjusted Rand Index ∈  (0, 1). Metric for clustering.
  # Higher is better.
  open("$IMGDIR/ari.txt", "w") do file
    ariCytof = [ x[1] for x in util.ari.(simdat[:lam], lastState.lam) ]
    write(file, "ARI lam: $ariCytof\n")
  end

  println("Make png...")

  # Set png resolution settings
  s_png = 10
  fy(clus) = util.addCut(clus, s_png)
  fZ(Z) = util.addGridLines(Z, s_png)

  for i in 1:I
    util.plotPng("$IMGDIR/y_imputed$(i).png", s=s_png)
    util.yZ_inspect(out[1], i=i, lastState.y_imputed, zlim=[-4,4], using_zero_index=false,
                    thresh=thresh, fy=fy, fZ=fZ, col=util.blueToRed(9))
    util.devOff()

    util.plotPng("$IMGDIR/y_dat$(i).png", s=s_png)
    util.yZ_inspect(out[1], i=i, simdat[:y], zlim=[-4,4], na="black", using_zero_index=false,
                    thresh=thresh, fy=fy, fZ=fZ, col=util.blueToRed(9))
    util.devOff()
  end
  println("Done with png...")

  #= Plots.jl
  for i in 1:I
    CytofImg.yZ_inspect(out[1], lastState.y_imputed, i, thresh=.9)
    CytofImg.Plots.savefig("$IMGDIR/y_imputed$(i).png")

    CytofImg.yZ_inspect(out[1], simdat[:y], i, thresh=.9)
    CytofImg.Plots.savefig("$IMGDIR/y_dat$(i).png")
  end
  =#

  # Missing Mechanism
  open("$IMGDIR/beta.txt", "w") do file
    for i in 1:I
      bi = join(c.beta[:, i], ", ")
      write(file, "for i=$i, beta = $bi \n")
    end
  end

  # Plot missing mechanism
  util.plotPdf("$IMGDIR/prob_miss.pdf")
  R"par(mfrow=c($(I), 1))"
  for i in 1:I
    util.plotProbMiss(c.beta, i)
  end
  R"par(mfrow=c(1,1))"
  util.devOff()

  # TODO: DO THIS
  # Plot QQ
  # util.plotPdf("$IMGDIR/qq.pdf")
  # R"par(mfrow=c(3, 3), mar=c(5.1, 4, 2, 1))"
  # y_obs_range = util.y_obs_range(y_dat)
  # for i in 1:I
  #   for j in 1:J
  #     print("\ri: $i, j: $j")
  #     # QQ of observed expression levels
  #     y_obs, y_pp = util.qq_yobs_postpred(y_dat, i, j, out)
  #     util.myQQ(y_obs, y_pp, pch=20, ylab="post pred quantiles",
  #               xlab="y (observed) quantiles", main="i: $i, j: $j", xlim=y_obs_range,
  #               ylim=y_obs_range)
  #   end
  # end
  # R"par(mfrow=c(1, 1), mar=mar.default())"
  # println()
  # util.devOff()

  # TODO: consider removing this
  # # QQ with observed values
  # util.plotPdf("$IMGDIR/qq_truedat.pdf")
  # R"par(mfrow=c(3, 3), mar=c(5.1, 4, 2, 1))"
  # y_obs_range = util.y_obs_range(simdat[:y_complete])
  # for i in 1:I
  #   for j in 1:J
  #     println("i: $i, j: $j")
  #     # QQ of observed expression levels
  #     y_obs, y_pp = util.qq_yobs_postpred(simdat[:y_complete], i, j, out)
  #     util.myQQ(y_obs, y_pp, pch=20, ylab="post pred quantiles", xlab="y
  #               (observed) quantiles", main="i: $i, j: $j", xlim=y_obs_range,
  #               ylim=y_obs_range)
  #   end
  # end
  # R"par(mfrow=c(1, 1), mar=mar.default())"
  # util.devOff()

  # Plot Posterior Density for i, j observed
  mkpath("$(IMGDIR)/dden")
  for i in 1:I
    for j in 1:J
      println("Plot posterior density for (i: $i, j: $j) observed ...")
      util.plotPdf("$(IMGDIR)/dden/dden_i$(i)_j$(j).pdf")
      dyij = util.density(filter(yij -> !isnan(yij), y_dat[i][:, j]), bw=.1)
      dd_ij = hcat([dd[i, j] for dd in dden]...)
      pdyij_mean = R"rowMeans($dd_ij)" .+ 0
      pdyij_lower = R"apply($dd_ij, 1, quantile, .025)" .+ 0
      pdyij_upper = R"apply($dd_ij, 1, quantile, .975)" .+ 0
      h = maximum([dyij[:y] .+ 0; pdyij_upper])
      util.plot(c.y_grid, pdyij_mean, xlab="y", ylab="density", ylim=[0, h],
                main="i: $i, j: $j", col="blue", lwd=2, typ="l")
      util.colorBtwn(c.y_grid, pdyij_lower, pdyij_upper, from=-10, to=10,
                     col=util.rgba("blue", .3))
      util.lines(dyij)
      util.devOff()
    end
  end

  # Separate graphs
  mkpath("$IMGDIR/sep")
  for i in 1:I
    println("Separate graphs for yZ $(i)...")
    idx_best = R"estimate_ZWi_index($(out[1]), $i)"[1]
    Zi = out[1][idx_best][:Z]
    Wi = out[1][idx_best][:W][i,:]

    # Point Est for Wi and Zi. txt and pdf.
    open("$IMGDIR/sep/W_$(i)_hat.txt", "w") do file
      write(file, "$(join(Wi, "\n"))\n")
    end

    open("$IMGDIR/sep/W$(i)_hat_ordered_cumsum.txt", "w") do file
      ord = sortperm(Wi, rev=true)
      cs_wi_sorted = cumsum(Wi[ord])
      write(file, "num_features,k,wi,cumprop\n")
      for k in 1:K_MCMC
        write(file, "$(k),$(ord[k]),$(Wi[ord][k]),$(cs_wi_sorted[k])\n")
      end
    end

    open("$IMGDIR/sep/Z$(i)_hat.txt", "w") do file
      for j in 1:J
        zj = join(Int.(Zi[j, :]), ",")
        write(file, "$(zj)\n")
      end
    end

    lami = out[1][idx_best][:lam][i]
    ord = sortperm(Wi, rev=true)
    lami = util.reorder_lami(ord, lami)

    for min_presence in min_presences
      common_celltypes = util.get_common_celltypes(Wi, thresh=min_presence,
                                                   filter_by_min_presence=true)

      if all(isnan.(lami))
        println("Sample too noisy.")
      else
        println("common celltypes (min_presence > $min_presence): $common_celltypes")
        K_trunc = length(common_celltypes)

        util.plotPng("$IMGDIR/sep/y_dat$(i)_only_minpresence$(min_presence).png", s=s_png)
        ord_yi = sortperm(lami)
        function f(y)
          util.addCut(lami, s_png)
          R"""
          axis(1, at=1:$J, label=1:$J, las=2, fg="grey", cex.axis=1.5)
          """
        end
        util.myImage(y_dat[i][ord_yi[1 .<= lami[ord_yi] .<= K_trunc], :],
                     addL=true, f=f,
                     zlim=[-4,4], col=util.blueToRed(9), na="black", xlab="markers",
                     ylab="cells", xaxt="n", var"cex.lab"=1.5, var"cex.y.leg"=1.5);
        util.devOff()

        util.plotPdf("$IMGDIR/sep/Z_hat$(i)_minpresence$(min_presence).pdf", w=5, h=10)
        util.myImage(Zi[:, common_celltypes], addL=false, ylab="markers", yaxt="n",
                     f=Z->addGridLines(J, K_trunc), xaxt="n", xlab="celltypes",
                     var"cex.lab"=1.5);

        perc = string.(round.(Wi[common_celltypes] * 100, digits=2), "%")
        R"""
        axis(3, at=1:$K_trunc, label=$(perc), las=2, fg="grey", cex.axis=1)
        axis(1, at=1:$K_trunc, label=$(common_celltypes), las=1,
             fg="grey", cex.axis=1)
        axis(2, at=1:$J, label=1:$J, las=2, fg="grey", cex.axis=1.3)
        """
        util.devOff()

        util.plotPdf("$IMGDIR/sep/ZT_hat$(i)_minpresence$(min_presence).pdf")
        perc = string.(round.(Wi[common_celltypes] * 100, digits=2), "%")
        function fz(z)
          addGridLines(K_trunc, J)
          R"""
          axis(4, at=1:$K_trunc, label=$(perc), las=1, fg="grey", cex.axis=1.3)
          axis(2, at=1:$K_trunc, label=$(common_celltypes), las=1,
               fg="grey", cex.axis=1.3)
          axis(1, at=1:$J, label=1:$J, las=2, fg="grey", cex.axis=1.3)
          """
        end
        R"par(mar=c(5.1, 4, 2.1, 5.1))"
        util.myImage(Zi[:, common_celltypes]', addL=false, ylab="celltypes", yaxt="n",
                     f=fz, xaxt="n", xlab="markers", var"cex.lab"=1.5);
        R"par(mar=c(5, 4, 4, 2) + .1)"
        util.devOff()
      end
    end
  end # separate graphs

end


text:
include("setup_1d.jl")

println("dwt! by lifting (N=",N,"), ", L, " levels")
f(x0, L, wl) = for i = 1:tn; dwt!(x0, wl, L); end
f(x0, L, wl);
@time f(x0, L, wl);


text:
using Revise # reduce recompilation time
using LinearAlgebra
using SparseArrays
using BenchmarkTools
using UnPack
using ToeplitzMatrices
using Test
using DelimitedFiles

using KernelAbstractions
using CUDA
using CUDAapi
CUDA.allowscalar(false)

push!(LOAD_PATH, "./src")
using CommonUtils
using Basis1D
using Basis2DTri
using UniformTriMesh

push!(LOAD_PATH, "./examples/EntropyStableEuler")
using EntropyStableEuler

using SetupDG

const Nd = 5
S_N(x) = @. sin(pi*x/h)/(2*pi/h)/tan(x/2)
"""
Vandermonde matrix of sinc basis functions determined by h,
evaluated at r
"""
function vandermonde_Sinc(h,r)
    N = convert(Int, 2*pi/h)
    V = zeros(length(r),N)
    for n = 1:N
        V[:,n] = S_N(r.-n*h)
    end
    V[1,1] = 1
    V[end,end] = 1
    return V
end

vector_norm(U) = CuArrays.sum((x->x.^2).(U))

# ================================= #
# ============ Routines =========== #
# ================================= #
function u_to_v!(VU,Q,Nd,Nq)
    index = (blockIdx().x-1)*blockDim().x + threadIdx().x
    stride = blockDim().x*gridDim().x

    @inbounds begin
    for i = index:stride:ceil(Int,length(VU)/Nd)
        k = div(i-1,Nq) # Current element
        n = mod1(i,Nq) # Current quad node

        idx = k*Nd*Nq+n
        rho = Q[idx]
        rhou = Q[idx+Nq]
        rhov = Q[idx+2*Nq]
        rhow = Q[idx+3*Nq]
        E = Q[idx+4*Nq]
        rhoUnorm = rhou^2+rhov^2+rhow^2
        rhoe = E-.5*rhoUnorm/rho
        sU = CUDA.log((gamma-1.0f0)*rhoe/CUDA.exp(gamma*CUDA.log(rho)))

        VU[idx] = (-E+rhoe*(gamma+1.0f0-sU))/rhoe
        VU[idx+Nq] = rhou/rhoe
        VU[idx+2*Nq] = rhov/rhoe
        VU[idx+3*Nq] = rhow/rhoe
        VU[idx+4*Nq] = -rho/rhoe
    end
    end
end

function v_to_u!(Qh,Nd,Nh)
    index = (blockIdx().x-1)*blockDim().x + threadIdx().x
    stride = blockDim().x*gridDim().x

    @inbounds begin
    for i = index:stride:ceil(Int,length(Qh)/Nd)
        k = div(i-1,Nh) # Current element
        n = mod1(i,Nh) # Current hybridized node

        idx = k*Nd*Nh+n
        v1 = Qh[idx     ]
        vu = Qh[idx+  Nh]
        vv = Qh[idx+2*Nh]
        vw = Qh[idx+3*Nh]
        vE = Qh[idx+4*Nh]
        vUnorm = vu^2+vv^2+vw^2
        rhoeV = CUDA.exp(1.0f0/(gamma-1.0f0)*CUDA.log((gamma-1.0f0)/CUDA.exp(gamma*CUDA.log(-vE))))*CUDA.exp(-(gamma-v1+vUnorm/(2.0f0*vE))/(gamma-1.0f0))

        Qh[idx     ] = -rhoeV*vE
        Qh[idx+  Nh] = rhoeV*vu
        Qh[idx+2*Nh] = rhoeV*vv
        Qh[idx+3*Nh] = rhoeV*vw
        Qh[idx+4*Nh] = rhoeV*(1.0f0-vUnorm/(2.0f0*vE))
    end
    end
end

#TODO: combine with v_to_u
function u_to_primitive!(Qh,Nd,Nh)
    index = (blockIdx().x-1)*blockDim().x + threadIdx().x
    stride = blockDim().x*gridDim().x

    @inbounds begin
    for i = index:stride:ceil(Int,length(Qh)/Nd)
        k = div(i-1,Nh) # Current element
        n = mod1(i,Nh) # Current hybridized node
        idx = k*Nd*Nh+n

        rho = Qh[idx]
        rhou = Qh[idx+Nh]
        rhov = Qh[idx+2*Nh]
        rhow = Qh[idx+3*Nh]
        E = Qh[idx+4*Nh]

        beta = rho/(2.0f0*(gamma-1.0f0)*(E-.5*(rhou^2+rhov^2+rhow^2)/rho))

        Qh[idx+Nh] = rhou/rho
        Qh[idx+2*Nh] = rhov/rho
        Qh[idx+3*Nh] = rhow/rho
        Qh[idx+4*Nh] = beta
    end
    end
end

function extract_face_val_conservative!(Uf,Qh,Nfp_P,Nq_P,Nh_P)
    index = (blockIdx().x-1)*blockDim().x + threadIdx().x
    stride = blockDim().x*gridDim().x

    @inbounds begin
    for i = index:stride:length(Uf)
        k = div(i-1,Nfp_P)
        n = mod1(i,Nfp_P)
        Uf[i] = Qh[k*Nh_P+Nq_P+n]
    end
    end
end

function extract_face_val!(QM,Qh,Nfp_P,Nq_P,Nh_P)
    index = (blockIdx().x-1)*blockDim().x + threadIdx().x
    stride = blockDim().x*gridDim().x

    @inbounds begin
    for i = index:stride:length(QM)
        k = div(i-1,Nfp_P)
        n = mod1(i,Nfp_P)
        QM[i] = Qh[k*Nh_P+Nq_P+n]
    end
    end
end

function CU_logmean(uL,uR,logL,logR)
    da = uR-uL
    aavg = .5f0*(uL+uR)
    f = da/aavg
    if CUDA.abs(f)<1e-4
        v = f^2
        return aavg*(1.0f0 + v*(-.2f0-v*(.0512f0 - v*0.026038857142857f0)))
    else
        return -da/(logL-logR)
    end
end

function CU_euler_flux(rhoM,uM,vM,wM,betaM,rhoP,uP,vP,wP,betaP,rhologM,betalogM,rhologP,betalogP)
    rholog  = CU_logmean(rhoM,rhoP,rhologM,rhologP)
    betalog = CU_logmean(betaM,betaP,betalogM,betalogP)

    rhoavg = .5f0*(rhoM+rhoP)
    uavg   = .5f0*(uM+uP)
    vavg   = .5f0*(vM+vP)
    wavg   = .5f0*(wM+wP)

    unorm = uM*uP+vM*vP+wM*wP
    pa = rhoavg/(betaM+betaP)
    E_plus_p = rholog/(2.0f0*(gamma-1.0f0)*betalog) + pa + .5f0*rholog*unorm

    FxS1 = (@. rholog*uavg)
    FxS2 = (@. FxS1*uavg + pa)
    FxS3 = (@. FxS1*vavg) # rho * u * v
    FxS4 = (@. FxS1*wavg) # rho * u * w
    FxS5 = (@. E_plus_p*uavg)
    FyS1 = (@. rholog*vavg)
    FyS2 = (@. FxS3) # rho * u * v
    FyS3 = (@. FyS1*vavg + pa)
    FyS4 = (@. FyS1*wavg) # rho * v * w
    FyS5 = (@. E_plus_p*vavg)
    FzS1 = (@. rholog*wavg)
    FzS2 = (@. FxS4) # rho * u * w
    FzS3 = (@. FyS4) # rho * v * w
    FzS4 = (@. FzS1*wavg + pa) # rho * w^2 + p
    FzS5 = (@. E_plus_p*wavg)

    return FxS1,FxS2,FxS3,FxS4,FxS5,FyS1,FyS2,FyS3,FyS4,FyS5,FzS1,FzS2,FzS3,FzS4,FzS5
end

function construct_lam!(lam,Uf,nxJ,nyJ,sJ,Nfp,K,Nd)
    index = (blockIdx().x-1)*blockDim().x + threadIdx().x
    stride = blockDim().x*gridDim().x

    @inbounds begin
    for i = index:stride:Nfp*K
        k = div(i-1,Nfp)
        n = mod1(i,Nfp)

        rhoM  = Uf[k*Nd*Nfp+n      ]
        rhouM = Uf[k*Nd*Nfp+n+  Nfp]
        rhovM = Uf[k*Nd*Nfp+n+2*Nfp]
        rhowM = Uf[k*Nd*Nfp+n+3*Nfp]
        EM    = Uf[k*Nd*Nfp+n+4*Nfp]
        nxJ_val = nxJ[k*Nfp+n]
        nyJ_val = nyJ[k*Nfp+n]
        sJ_val  = sJ[k*Nfp+n]
        rhouM_n = (rhouM*nxJ_val+rhovM*nyJ_val)/sJ_val
        lam[k*Nfp+n] = CUDA.sqrt(CUDA.abs(rhouM_n/rhoM))+CUDA.sqrt(gamma*(gamma-1.0f0)*(EM-0.5f0*rhouM_n^2/rhoM)/rhoM)
    end
    end
end

function surface_kernel!(flux,QM,QP,Uf,UfP,LFc,nxJ,nyJ,Nfp,K,Nd)
    index = (blockIdx().x-1)*blockDim().x + threadIdx().x
    stride = blockDim().x*gridDim().x
    @inbounds begin
    for i = index:stride:Nfp*K
        k = div(i-1,Nfp)
        n = mod1(i,Nfp)

        rhoM  = QM[k*Nd*Nfp+n      ]
        uM    = QM[k*Nd*Nfp+n+  Nfp]
        vM    = QM[k*Nd*Nfp+n+2*Nfp]
        wM    = QM[k*Nd*Nfp+n+3*Nfp]
        betaM = QM[k*Nd*Nfp+n+4*Nfp]
        rhoP  = QP[k*Nd*Nfp+n      ]
        uP    = QP[k*Nd*Nfp+n+  Nfp]
        vP    = QP[k*Nd*Nfp+n+2*Nfp]
        wP    = QP[k*Nd*Nfp+n+3*Nfp]
        betaP = QP[k*Nd*Nfp+n+4*Nfp]
        nxJ_val = nxJ[k*Nfp+n]
        nyJ_val = nyJ[k*Nfp+n]

        rhouM = Uf[k*Nd*Nfp+n+  Nfp]
        rhovM = Uf[k*Nd*Nfp+n+2*Nfp]
        rhowM = Uf[k*Nd*Nfp+n+3*Nfp]
        EM    = Uf[k*Nd*Nfp+n+4*Nfp]
        rhouP = UfP[k*Nd*Nfp+n+  Nfp]
        rhovP = UfP[k*Nd*Nfp+n+2*Nfp]
        rhowP = UfP[k*Nd*Nfp+n+3*Nfp]
        EP    = UfP[k*Nd*Nfp+n+4*Nfp]
        LFc_val = LFc[k*Nfp+n]

        rhologM = CUDA.log(rhoM)
        rhologP = CUDA.log(rhoP)
        betalogM = CUDA.log(betaM)
        betalogP = CUDA.log(betaP)

        FxS1,FxS2,FxS3,FxS4,FxS5,FyS1,FyS2,FyS3,FyS4,FyS5,_,_,_,_,_ = CU_euler_flux(rhoM,uM,vM,wM,betaM,rhoP,uP,vP,wP,betaP,rhologM,betalogM,rhologP,betalogP)

        flux[k*Nd*Nfp+n      ] = nxJ_val*FxS1+nyJ_val*FyS1-LFc_val*(rhoP-rhoM)
        flux[k*Nd*Nfp+n+  Nfp] = nxJ_val*FxS2+nyJ_val*FyS2-LFc_val*(rhouP-rhouM)
        flux[k*Nd*Nfp+n+2*Nfp] = nxJ_val*FxS3+nyJ_val*FyS3-LFc_val*(rhovP-rhovM)
        flux[k*Nd*Nfp+n+3*Nfp] = nxJ_val*FxS4+nyJ_val*FyS4-LFc_val*(rhowP-rhowM)
        flux[k*Nd*Nfp+n+4*Nfp] = nxJ_val*FxS5+nyJ_val*FyS5-LFc_val*(EP-EM)
    end
    end
end

function flux_differencing_xy_kernel!(gradfh,Qh,rxJ,sxJ,ryJ,syJ,Qrh_skew,Qsh_skew,Nh,Nh_P,Np_F,K,Nd,Nq_P,num_threads,Nk_max,Nk_tri_max)
    index = (blockIdx().x-1)*blockDim().x + threadIdx().x
    stride = blockDim().x*gridDim().x
    #TODO: using stride to avoid possible insufficient blocks?
    tid_start = (blockIdx().x-1)*blockDim().x + 1 # starting threadid of current block
    tid_end = (blockIdx().x-1)*blockDim().x + num_threads # ending threadid of current block
    k_start = div(tid_start-1,Nh) # starting elementid of current block
    k_end = tid_end < Nh*K ? div(tid_end-1,Nh) : K-1 # ending elementid of current block
    Nk = k_end-k_start+1 # Number of elements in current block
    k_tri_start = div(tid_start-1,Nh_P) # starting triangle id of current block
    k_tri_end = tid_end < Nh*K ? div(tid_end-1,Nh_P) : Np_F*K-1 # ending triangle id of current block
    Nk_tri = k_tri_end-k_tri_start+1 # Number of triangles in current block

    # Parallel read
    Qh_shared = @cuDynamicSharedMem(NUM_TYPE,Nh_P*Nd*Nk_tri_max)
    load_size = div(Nh_P*Nd*Nk_tri-1,num_threads)+1 # size of read of each thread
    @inbounds begin
    for i in (threadIdx().x-1)*load_size+1:min(threadIdx().x*load_size,Nh_P*Nd*Nk_tri)
        offset_tri = div(i-1,Nh_P*Nd) # distance with k_tri_start
        n_local = mod1(i,Nh_P*Nd)
        n_d = div(n_local-1,Nh_P) # Current component
        m_xy = mod1(n_local,Nh_P) # Current hybridized node index at xy slice
        i_tri = k_tri_start+offset_tri # Current triangle
        k = div(i_tri,Np_F) # Current element
        i_xy = mod(i_tri,Np_F) # Current local xy slice id on the elment
        Qh_shared[i] = Qh[k*Nd*Nh+i_xy*Nh_P+m_xy+n_d*Nh]
    end
    end

    sync_threads()

    @inbounds begin
    if index <= Nh*K
        i = index
        k = div(i-1,Nh)   # Current element
        m = mod1(i,Nh)    # Current node on element
        i_xy = div(m-1,Nh_P) # Current local xy slice id on the element
        m_xy = mod1(i,Nh_P)  # Current hybridized node index at x-y slice

        k_tri_offset = div(i-1,Nh_P)-k_tri_start
        rhoL  = Qh_shared[k_tri_offset*Nd*Nh_P+m_xy       ]
        uL    = Qh_shared[k_tri_offset*Nd*Nh_P+m_xy+  Nh_P]
        vL    = Qh_shared[k_tri_offset*Nd*Nh_P+m_xy+2*Nh_P]
        wL    = Qh_shared[k_tri_offset*Nd*Nh_P+m_xy+3*Nh_P]
        betaL = Qh_shared[k_tri_offset*Nd*Nh_P+m_xy+4*Nh_P]
        rhologL  = CUDA.log(rhoL)
        betalogL = CUDA.log(betaL)

        rho_sum = 0.0f0
        u_sum = 0.0f0
        v_sum = 0.0f0
        w_sum = 0.0f0
        beta_sum = 0.0f0

        # Assume Affine meshes

        rxJ_val = rxJ[k+1]
        sxJ_val = sxJ[k+1]
        ryJ_val = ryJ[k+1]
        syJ_val = syJ[k+1]

        for n_xy = 1:Nh_P
            if m_xy <= Nq_P || n_xy <= Nq_P
                rhoR  = Qh_shared[k_tri_offset*Nd*Nh_P+n_xy       ]
                uR    = Qh_shared[k_tri_offset*Nd*Nh_P+n_xy+  Nh_P]
                vR    = Qh_shared[k_tri_offset*Nd*Nh_P+n_xy+2*Nh_P]
                wR    = Qh_shared[k_tri_offset*Nd*Nh_P+n_xy+3*Nh_P]
                betaR = Qh_shared[k_tri_offset*Nd*Nh_P+n_xy+4*Nh_P]
                rhologR = CUDA.log(rhoR)
                betalogR = CUDA.log(betaR)

                FxS1,FxS2,FxS3,FxS4,FxS5,FyS1,FyS2,FyS3,FyS4,FyS5,FzS1,FzS2,FzS3,FzS4,FzS5 = CU_euler_flux(rhoL,uL,vL,wL,betaL,rhoR,uR,vR,wR,betaR,rhologL,betalogL,rhologR,betalogR)

                Qx_val = 2.0f0*(rxJ_val*Qrh_skew[m_xy,n_xy]+sxJ_val*Qsh_skew[m_xy,n_xy])
                Qy_val = 2.0f0*(ryJ_val*Qrh_skew[m_xy,n_xy]+syJ_val*Qsh_skew[m_xy,n_xy])
                rho_sum  += Qx_val*FxS1+Qy_val*FyS1
                u_sum    += Qx_val*FxS2+Qy_val*FyS2
                v_sum    += Qx_val*FxS3+Qy_val*FyS3
                w_sum    += Qx_val*FxS4+Qy_val*FyS4
                beta_sum += Qx_val*FxS5+Qy_val*FyS5
            end
        end

        gradfh[k*Nd*Nh+m     ] = rho_sum
        gradfh[k*Nd*Nh+m+  Nh] = u_sum
        gradfh[k*Nd*Nh+m+2*Nh] = v_sum
        gradfh[k*Nd*Nh+m+3*Nh] = w_sum
        gradfh[k*Nd*Nh+m+4*Nh] = beta_sum
    end
    end
    return
end

function flux_differencing_z_kernel!(gradfh,Qh,wq,h,J,Qth,Nh,Nq,Nh_P,Np_F,K,Nd,Nq_P,num_threads,Nk_max,Nk_fourier_max)
    index = (blockIdx().x-1)*blockDim().x + threadIdx().x
    stride = blockDim().x*gridDim().x
    #TODO: using stride to avoid possible insufficient blocks?
    tid_start = (blockIdx().x-1)*blockDim().x + 1 # starting threadid of current block
    tid_end = (blockIdx().x-1)*blockDim().x + num_threads # ending threadid of current block
    k_start = div(tid_start-1,Nq) # starting elementid of current block
    k_end = tid_end < Nq*K ? div(tid_end-1,Nq) : K-1 # ending elementid of current block
    Nk = k_end-k_start+1 # Number of elements in current block
    k_fourier_start = div(tid_start-1,Np_F)
    k_fourier_end = tid_end < Nq*K ? div(tid_end-1,Np_F) : Nq_P*K-1
    Nk_fourier = k_fourier_end-k_fourier_start+1

    Qh_shared = @cuDynamicSharedMem(NUM_TYPE,Np_F*Nd*Nk_fourier_max)
    load_size = div(Np_F*Nd*Nk_fourier-1,num_threads)+1 # size of read of each thread
    @inbounds begin
    for i in (threadIdx().x-1)*load_size+1:min(threadIdx().x*load_size,Np_F*Nd*Nk_fourier)
        offset_fourier = div(i-1,Np_F*Nd) # distance with k_fourier_start
        n_local = mod1(i,Np_F*Nd)
        n_d = div(n_local-1,Np_F) # Current component
        m_z = mod1(n_local,Np_F) # Current local node index at fourier slice
        i_fourier = k_fourier_start+offset_fourier
        k = div(i_fourier,Nq_P) # Current element
        i_z = mod(i_fourier,Nq_P) # Current local z slice id on the element
        Qh_shared[i] = Qh[k*Nd*Nh+(m_z-1)*Nh_P+i_z+1+n_d*Nh]
    end
    end
    sync_threads()

    @inbounds begin
    if index <= Nq*K
        i = index
        k = div(i-1,Nq) # Current element
        m = mod1(i,Nq) # local node number in current element
        i_z = div(m-1,Np_F) # Fourier element id on current element
        m_z = mod1(m,Np_F) # local node number on Fourier element

        k_fourier_offset = div(i-1,Np_F)-k_fourier_start
        rhoL  = Qh_shared[k_fourier_offset*Nd*Np_F+m_z       ]
        uL    = Qh_shared[k_fourier_offset*Nd*Np_F+m_z+  Np_F]
        vL    = Qh_shared[k_fourier_offset*Nd*Np_F+m_z+2*Np_F]
        wL    = Qh_shared[k_fourier_offset*Nd*Np_F+m_z+3*Np_F]
        betaL = Qh_shared[k_fourier_offset*Nd*Np_F+m_z+4*Np_F]
        rhologL = CUDA.log(rhoL)
        betalogL = CUDA.log(betaL)

        rho_sum = 0.0f0
        u_sum = 0.0f0
        v_sum = 0.0f0
        w_sum = 0.0f0
        beta_sum = 0.0f0

        for n_z = 1:Np_F
            rhoR  = Qh_shared[k_fourier_offset*Nd*Np_F+n_z       ]
            uR    = Qh_shared[k_fourier_offset*Nd*Np_F+n_z+  Np_F]
            vR    = Qh_shared[k_fourier_offset*Nd*Np_F+n_z+2*Np_F]
            wR    = Qh_shared[k_fourier_offset*Nd*Np_F+n_z+3*Np_F]
            betaR = Qh_shared[k_fourier_offset*Nd*Np_F+n_z+4*Np_F]

            rhologR = CUDA.log(rhoR)
            betalogR = CUDA.log(betaR)

            _,_,_,_,_,_,_,_,_,_,FzS1,FzS2,FzS3,FzS4,FzS5 = CU_euler_flux(rhoL,uL,vL,wL,betaL,rhoR,uR,vR,wR,betaR,rhologL,betalogL,rhologR,betalogR)

            Qz_val = 2.0f0/J*wq[i_z+1]*Qth[m_z,n_z]
            rho_sum += Qz_val*FzS1
            u_sum += Qz_val*FzS2
            v_sum += Qz_val*FzS3
            w_sum += Qz_val*FzS4
            beta_sum += Qz_val*FzS5
        end
        gradfh[k*Nd*Nh+(m_z-1)*Nh_P+i_z+1     ] += rho_sum
        gradfh[k*Nd*Nh+(m_z-1)*Nh_P+i_z+1+  Nh] += u_sum
        gradfh[k*Nd*Nh+(m_z-1)*Nh_P+i_z+1+2*Nh] += v_sum
        gradfh[k*Nd*Nh+(m_z-1)*Nh_P+i_z+1+3*Nh] += w_sum
        gradfh[k*Nd*Nh+(m_z-1)*Nh_P+i_z+1+4*Nh] += beta_sum
    end
    end
    return
end

function volume_kernel!(gradfh,Qh,rxJ,sxJ,ryJ,syJ,wq,h,J,Qrh_skew,Qsh_skew,Qth,Nh,Nh_P,Np_F,K,Nd,Nq_P,num_threads,Nk_max)
    index = (blockIdx().x-1)*blockDim().x + threadIdx().x
    stride = blockDim().x*gridDim().x
    #TODO: using stride to avoid possible insufficient blocks?
    tid_start = (blockIdx().x-1)*blockDim().x + 1 # starting threadid of current block
    tid_end = (blockIdx().x-1)*blockDim().x + num_threads # ending threadid of current block
    k_start = div(tid_start-1,Nh) # starting elementid of current block
    k_end = tid_end < Nh*K ? div(tid_end-1,Nh) : K-1 # ending elementid of current block
    Nk = k_end-k_start+1 # Number of elements in current block

    # Parallel read
    Qh_shared = @cuDynamicSharedMem(NUM_TYPE,Nh*Nd*Nk_max)
    load_size = div(Nh*Nd*Nk-1,num_threads)+1 # size of read of each thread

    @inbounds begin
    for r_id in (threadIdx().x-1)*load_size+1:min(threadIdx().x*load_size,Nh*Nd*Nk)
        Qh_shared[r_id] = Qh[k_start*Nd*Nh+r_id]
    end
    end
    sync_threads()

    @inbounds begin
    if index <= Nh*K
        i = index
        k = div(i-1,Nh)
        m = mod1(i,Nh)

        k_offset = k-k_start

        rhoL  = Qh_shared[k_offset*Nd*Nh+m     ]
        uL    = Qh_shared[k_offset*Nd*Nh+m+  Nh]
        vL    = Qh_shared[k_offset*Nd*Nh+m+2*Nh]
        wL    = Qh_shared[k_offset*Nd*Nh+m+3*Nh]
        betaL = Qh_shared[k_offset*Nd*Nh+m+4*Nh]
        rhologL = CUDA.log(rhoL)
        betalogL = CUDA.log(betaL)

        t = div(m-1,Nh_P) # Current x-y slice
        s = mod1(m,Nh_P) # Current hybridized node index at x-y slice
        xy_idx = t*Nh_P+1:(t+1)*Nh_P # Nonzero index for Qrh, Qsh
        z_idx = s:Nh_P:s+(Np_F-1)*Nh_P

        rho_sum = 0.0f0
        u_sum = 0.0f0
        v_sum = 0.0f0
        w_sum = 0.0f0
        beta_sum = 0.0f0

        # Assume Affine meshes

        rxJ_val = rxJ[k+1]
        sxJ_val = sxJ[k+1]
        ryJ_val = ryJ[k+1]
        syJ_val = syJ[k+1]

        # TODO: better way to indexing
        for n = 1:Nh
            if n in xy_idx || n in z_idx

                rhoR  = Qh_shared[k_offset*Nd*Nh+n     ]
                uR    = Qh_shared[k_offset*Nd*Nh+n+  Nh]
                vR    = Qh_shared[k_offset*Nd*Nh+n+2*Nh]
                wR    = Qh_shared[k_offset*Nd*Nh+n+3*Nh]
                betaR = Qh_shared[k_offset*Nd*Nh+n+4*Nh]
                rhologR = CUDA.log(rhoR)
                betalogR = CUDA.log(betaR)

                FxS1,FxS2,FxS3,FxS4,FxS5,FyS1,FyS2,FyS3,FyS4,FyS5,FzS1,FzS2,FzS3,FzS4,FzS5 = CU_euler_flux(rhoL,uL,vL,wL,betaL,rhoR,uR,vR,wR,betaR,rhologL,betalogL,rhologR,betalogR)

                if n in xy_idx
                    col_idx = mod1(n,Nh_P)
                    Qx_val = 2.0f0*(rxJ_val*Qrh_skew[s,col_idx]+sxJ_val*Qsh_skew[s,col_idx])
                    Qy_val = 2.0f0*(ryJ_val*Qrh_skew[s,col_idx]+syJ_val*Qsh_skew[s,col_idx])
                    rho_sum += Qx_val*FxS1+Qy_val*FyS1
                    u_sum += Qx_val*FxS2+Qy_val*FyS2
                    v_sum += Qx_val*FxS3+Qy_val*FyS3
                    w_sum += Qx_val*FxS4+Qy_val*FyS4
                    beta_sum += Qx_val*FxS5+Qy_val*FyS5
                end
                #=
                if n in z_idx && s <= Nq_P
                    col_idx = div(n-1,Nh_P)+1
                    wqn = 2.0f0/J*wq[s]
                    Qz_val = wqn*Qth[t+1,col_idx]
                    rho_sum += Qz_val*FzS1
                    u_sum += Qz_val*FzS2
                    v_sum += Qz_val*FzS3
                    w_sum += Qz_val*FzS4
                    beta_sum += Qz_val*FzS5
                end
                =#
            end
        end

        gradfh[k*Nd*Nh+m     ] = rho_sum
        gradfh[k*Nd*Nh+m+  Nh] = u_sum
        gradfh[k*Nd*Nh+m+2*Nh] = v_sum
        gradfh[k*Nd*Nh+m+3*Nh] = w_sum
        gradfh[k*Nd*Nh+m+4*Nh] = beta_sum
    end
    end
    return
end

# ================================= #
# ============ Routines =========== #
# ================================= #

function rhs(Q,ops,mesh,param,num_threads,compute_rhstest,enable_test)
    Vq,wq,Qrh_skew,Qsh_skew,Qth,Ph,LIFTq,VPh,Wq = ops
    rxJ,sxJ,ryJ,syJ,nxJ,nyJ,sJ,J,h,mapP,mapP_d = mesh
    K,Np_P,Nq_P,Nfp_P,Nh_P,Np_F,Nq,Nfp,Nh,Nk_max,Nk_tri_max,Nk_fourier_max = param

    VU = CUDA.fill(CUDA.zero(NUM_TYPE),Nq*Nd*K)
    Uf = CUDA.fill(CUDA.zero(NUM_TYPE),Nfp*Nd*K)
    QM = CUDA.fill(CUDA.zero(NUM_TYPE),Nfp*Nd*K)
    flux = CUDA.fill(CUDA.zero(NUM_TYPE),Nfp*K*Nd)
    gradfh = CUDA.fill(CUDA.zero(NUM_TYPE),Nh*Nd*K)
    lam = CUDA.fill(CUDA.zero(NUM_TYPE),Nfp*K)
    LFc = CUDA.fill(CUDA.zero(NUM_TYPE),Nfp*K)

    # Entropy Projection
    @cuda threads=num_threads blocks=ceil(Int,Nq*K/num_threads) u_to_v!(VU,Q,Nd,Nq)
    synchronize()
    Qh = reshape(Ph*reshape(VU,Nq_P,Np_F*Nd*K),Nh*Nd*K)
    @cuda threads=num_threads blocks=ceil(Int,Nh*K/num_threads) v_to_u!(Qh,Nd,Nh)
    synchronize()
    @cuda threads=num_threads blocks=ceil(Int,Nfp*Nd*K/num_threads) extract_face_val_conservative!(Uf,Qh,Nfp_P,Nq_P,Nh_P)
    UfP = Uf[mapP_d]
    synchronize()
    @cuda threads=num_threads blocks=ceil(Int,Nh*K/num_threads) u_to_primitive!(Qh,Nd,Nh)
    synchronize()

    # Compute Surface values
    @cuda threads=num_threads blocks=ceil(Int,Nfp*Nd*K/num_threads) extract_face_val!(QM,Qh,Nfp_P,Nq_P,Nh_P)
    synchronize()
    QP = QM[mapP_d]

    # LF dissipation
    @cuda threads=num_threads blocks=ceil(Int,Nfp*K/num_threads) construct_lam!(lam,Uf,nxJ,nyJ,sJ,Nfp,K,Nd)
    synchronize()
    LFc .= .5f0*CUDA.max.(lam,lam[mapP]).*sJ

    # Surface kernel
    @cuda threads=num_threads blocks=ceil(Int,Nfp*K/num_threads) surface_kernel!(flux,QM,QP,Uf,UfP,LFc,nxJ,nyJ,Nfp,K,Nd)
    synchronize()
    flux = reshape(LIFTq*reshape(flux,Nfp_P,Np_F*Nd*K),Nq*Nd*K)

    # Volume kernel
    @cuda threads=num_threads blocks=ceil(Int,Nh*K/num_threads) shmem=sizeof(NUM_TYPE)*Nh_P*Nd*Nk_tri_max flux_differencing_xy_kernel!(gradfh,Qh,rxJ,sxJ,ryJ,syJ,Qrh_skew,Qsh_skew,Nh,Nh_P,Np_F,K,Nd,Nq_P,num_threads,Nk_max,Nk_tri_max)
    @cuda threads=num_threads blocks=ceil(Int,Nq*K/num_threads) shmem=sizeof(NUM_TYPE)*Np_F*Nd*Nk_fourier_max flux_differencing_z_kernel!(gradfh,Qh,wq,h,J,Qth,Nh,Nq,Nh_P,Np_F,K,Nd,Nq_P,num_threads,Nk_max,Nk_fourier_max)
    #@cuda threads=num_threads blocks=ceil(Int,Nh*K/num_threads) shmem=sizeof(NUM_TYPE)*Nh*Nd*Nk_max volume_kernel!(gradfh,Qh,rxJ,sxJ,ryJ,syJ,wq,h,J,Qrh_skew,Qsh_skew,Qth,Nh,Nh_P,Np_F,K,Nd,Nq_P,num_threads,Nk_max)
    synchronize()
    gradf = reshape(VPh*reshape(gradfh,Nh_P,Np_F*Nd*K),Nq*Nd*K)

    # Combine
    rhsQ = -(gradf+flux)

    # Compute rhstest
    rhstest = 0
    if compute_rhstest
        rhstest = CUDA.sum(Wq*reshape(VU.*rhsQ,Nq_P,Np_F*Nd*K))
    end

    if enable_test
        @show CUDA.maximum(VU)
        @show CUDA.minimum(VU)
        @show CUDA.sum(VU)
        @show CUDA.maximum(Qh)
        @show CUDA.minimum(Qh)
        @show CUDA.sum(Qh)
        @show CUDA.maximum(QM)
        @show CUDA.minimum(QM)
        @show CUDA.sum(QM)
        @show CUDA.maximum(QP)
        @show CUDA.minimum(QP)
        @show CUDA.sum(QP)
        @show CUDA.maximum(lam)
        @show CUDA.minimum(lam)
        @show CUDA.sum(lam)
        @show CUDA.maximum(LFc)
        @show CUDA.minimum(LFc)
        @show CUDA.sum(LFc)
        @show CUDA.maximum(flux)
        @show CUDA.minimum(flux)
        @show CUDA.sum(flux)
        @show CUDA.maximum(gradf)
        @show CUDA.minimum(gradf)
        @show CUDA.sum(gradf)
        @show CUDA.maximum(rhsQ)
        @show CUDA.minimum(rhsQ)
        @show CUDA.sum(rhsQ)
        @show rhstest
        @show typeof(Q)
        @show typeof(VU)
        @show typeof(Qh)
        @show typeof(QM)
        @show typeof(QP)
        @show typeof(flux)
        @show typeof(gradf)
        @show typeof(rhsQ)
    end
    return rhsQ,rhstest
end

"Constants"
const sp_tol      = 1e-12
const shmem_limit = 1024  # Max number of elements in shared mem
const enable_test = false
const errTol      = 1e-5  # Error tolerance in time stepping
const gamma       = 1.4f0
const NUM_TYPE    = Float64
"Program parameters"
const compute_L2_err          = false
const add_LF_dissipation      = true

N_P_arr = [3]
Np_F_arr = [2;4;6;8;10;12;14;16;18]
K1D_arr = [20]

errArr = zeros(length(K1D_arr),length(Np_F_arr),length(N_P_arr))

for N_P_idx = 1:length(N_P_arr)
    for Np_F_idx = 1:length(Np_F_arr)
        for K1D_arr_idx = 1:length(K1D_arr)
            "Approximation Parameters"
            N_P   = N_P_arr[N_P_idx];    # The order of approximation in polynomial dimension
            Np_P  = Int((N_P+1)*(N_P+2)/2)
            Np_F  = Np_F_arr[Np_F_idx];    # The order of approximation in Fourier dimension
            K1D   = K1D_arr[K1D_arr_idx];   # Number of elements in polynomial (x,y) dimension
            CFL   = 0.5;
            T     = 1.0;  # End time

            L_fourier = 10
            L_X = 10
            L_Y = 10

            println("N_P = $N_P, Np_F = $Np_F, K1D = $K1D")
            "Time integration Parameters"
            rk4a,rk4b,rk4c = rk45_coeffs()
            CN             = (N_P+1)*(N_P+2)*3/2  # estimated trace constant for CFL
            dt             = CFL * 2 / CN / K1D
            Nsteps         = convert(Int,ceil(T/dt))
            dt             = T/Nsteps

            "Initialize Reference Element in Fourier dimension"
            h       = 2*pi/Np_F
            column  = [0; .5*(-1).^(1:Np_F-1).*cot.((1:Np_F-1)*h/2)]
            column2 = [-pi^2/3/h^2-1/6; -((-1).^(1:Np_F-1)./(2*(sin.((1:Np_F-1)*h/2)).^2))]
            Dt      = Array{Float64,2}(Toeplitz(column,column[[1;Np_F:-1:2]]))
            D2t     = Array{Float64,2}(Toeplitz(column2,column2[[1;Np_F:-1:2]]))
            t       = LinRange(h,2*pi,Np_F)

            "Initialize Reference Element in polynomial dimension"
            rd = init_reference_tri(N_P);
            @unpack fv,Nfaces,r,s,VDM,V1,Dr,Ds,rf,sf,wf,nrJ,nsJ,rq,sq,wq,Vq,M,Pq,Vf,LIFT = rd
            Lq    = LIFT

            "Problem parameters"
            Nq_P   = length(rq)
            Nfp_P  = length(rf)
            Nh_P   = Nq_P+Nfp_P # Number of hybridized points
            Nq     = Nq_P*Np_F
            Nh     = Nh_P*Np_F
            Nfp    = Nfp_P*Np_F

            # Adpatively determine number of threads
            thread_count = 2^8
            for n = 8:-1:5
                thread_count = 2^n
                Nk_tri_max = div(thread_count-2,Nh_P)+2 # Max amount of triangles in a block
                Nk_fourier_max = div(thread_count-2,Np_F)+2 # Max amount of fourier slides in a block
                if max(Nd*Nh_P*Nk_tri_max,Nd*Np_F*Nk_fourier_max) < shmem_limit
                    break
                end
            end
            num_threads = thread_count
            Nk_max = num_threads == 1 ? 1 : div(num_threads-2,Nh)+2 # Max amount of element in a block
            Nk_tri_max = num_threads == 1 ? 1 : div(num_threads-2,Nh_P)+2 # Max amount of triangles in a block
            Nk_fourier_max = num_threads == 1 ? 1 : div(num_threads-2,Np_F)+2 # Max amount of fourier slides in a block

            "Mesh related variables"
            # Initialize 2D triangular mesh
            VX,VY,EToV = uniform_tri_mesh(K1D,K1D)
            @. VX = (1+VX)*L_X/2
            @. VY = (1+VY)*L_Y/2
            md   = init_mesh((VX,VY),EToV,rd)
            # Intialize triangular prism
            VX   = repeat(VX,2)
            VY   = repeat(VY,2)
            VZ   = [L_fourier/Np_F*ones((K1D+1)*(K1D+1),1); L_fourier*ones((K1D+1)*(K1D+1),1)]
            EToV = [EToV EToV.+(K1D+1)*(K1D+1)]

            # Make domain periodic
            @unpack Nfaces,Vf = rd
            @unpack xf,yf,K,mapM,mapP,mapB = md
            LX,LY = (x->maximum(x)-minimum(x)).((VX,VY)) # find lengths of domain
            mapPB = build_periodic_boundary_maps(xf,yf,LX,LY,Nfaces*K,mapM,mapP,mapB)
            mapP[mapB] = mapPB

            # Initialize 3D mesh
            @unpack x,y,xf,yf,xq,yq,rxJ,sxJ,ryJ,syJ,J,nxJ,nyJ,sJ,mapM,mapP,mapB = md
            x,y,xf,yf,xq,yq,rxJ,sxJ,ryJ,syJ,J,nxJ,nyJ,sJ = (x->reshape(repeat(x,inner=(1,Np_F)),size(x,1),Np_F,K)).((x,y,xf,yf,xq,yq,rxJ,sxJ,ryJ,syJ,J,nxJ,nyJ,sJ))
            z,zq,zf = (x->reshape(repeat(collect(L_fourier/Np_F:(L_fourier/Np_F):L_fourier),inner=(1,x),outer=(K,1))',x,Np_F,K)).((Np_P,Nq_P,Nfp_P))
            mapM = reshape(1:Nfp_P*Np_P*K,Nfp_P,Np_P,K)
            mapP_2D = (x->mod1(x,Nfp_P)+div(x-1,Nfp_P)*Nfp_P*Np_F).(mapP)
            mapP = reshape(repeat(mapP_2D,inner=(1,Np_F)),Nfp_P,Np_F,K)
            for j = 1:Np_F
                mapP[:,j,:] = mapP[:,j,:].+(j-1)*Nfp_P
            end
            mapP_d = zeros(Int32,Nfp_P,Np_F,Nd,K)
            for k = 1:K
                for j = 1:Np_F
                    for i = 1:Nfp_P
                        val = mapP[i,j,k]
                        elem = div(val-1,Nfp)
                        n = mod1(val,Nfp)
                        mapP_d[i,j,1,k] = elem*Nd*Nfp+n
                        mapP_d[i,j,2,k] = elem*Nd*Nfp+n+Nfp
                        mapP_d[i,j,3,k] = elem*Nd*Nfp+n+2*Nfp
                        mapP_d[i,j,4,k] = elem*Nd*Nfp+n+3*Nfp
                        mapP_d[i,j,5,k] = elem*Nd*Nfp+n+4*Nfp
                    end
                end
            end
            mapP = mapP[:]
            mapP_d = mapP_d[:]

            #############################################
            ######  Construct hybridized operators ######
            #############################################

            # scale by Fourier dimension
            M = h*M
            wq = h*wq
            wf = h*wf

            # Hybridized operators
            Vh = [Vq;Vf]
            rxJ,sxJ,ryJ,syJ = (x->mapslices((y->Vh*y),x,dims=(1,2))).((rxJ,sxJ,ryJ,syJ))
            Ef = Vf*Pq
            Br,Bs = (x->diagm(wf.*x)).((nrJ,nsJ))
            Qr,Qs = (x->Pq'*M*x*Pq).((Dr,Ds))
            Qrh,Qsh = (x->1/2*[x[1]-x[1]' Ef'*x[2];
                               -x[2]*Ef   x[2]]).(((Qr,Br),(Qs,Bs)))
            Qrh_skew,Qsh_skew = (x->1/2*(x-x')).((Qrh,Qsh))
            Qt = Dt
            Qth = Qt # Not the SBP operator, weighted when flux differencing
            Ph = [Vq;Vf]*Pq

            # TODO: assume mesh uniform affine, so Jacobian are constants
            # TODO: fix other Jacobian parts
            JP = L_X*L_Y/K1D/K1D/4
            JF = L_fourier/pi/2
            J = JF*JP
            wq = J*wq
            wf = JF*wf
            Lq = 1/JP*Lq
            Qrh = JF*Qrh
            Qsh = JF*Qsh
            Qth = JP*Qth
            Qrh_skew = 1/2*(Qrh-Qrh')
            Qsh_skew = 1/2*(Qsh-Qsh')
            LIFTq = Vq*Lq
            VPh = Vq*[Pq Lq]*diagm(1 ./ [wq;wf])
            Winv = diagm(1 ./ [wq;wf])
            Wq = diagm(wq)

            # Assume affine meshes
            rxJ = [rxJ[1,1,k] for k in 1:K]
            sxJ = [sxJ[1,1,k] for k in 1:K]
            ryJ = [ryJ[1,1,k] for k in 1:K]
            syJ = [syJ[1,1,k] for k in 1:K]

            # convert precision
            Vq,wq,Qrh_skew,Qsh_skew,Qth,Ph,LIFTq,VPh,Wq,rxJ,sxJ,ryJ,syJ,nxJ,nyJ,sJ,J,h,rk4a,rk4b = (A->convert.(NUM_TYPE,A)).((Vq,wq,Qrh_skew,Qsh_skew,Qth,Ph,LIFTq,VPh,Wq,rxJ,sxJ,ryJ,syJ,nxJ,nyJ,sJ,J,h,rk4a,rk4b))

            # TODO: refactor
            Vq,Vf,wq,wf,Pq,Lq,Qrh_skew,Qsh_skew,Qth,Ph,LIFTq,VPh,Wq,Winv,rxJ,sxJ,ryJ,syJ = (x->CuArray(x)).((Vq,Vf,wq,wf,Pq,Lq,Qrh_skew,Qsh_skew,Qth,Ph,LIFTq,VPh,Wq,Winv,rxJ,sxJ,ryJ,syJ))
            nxJ = CuArray(nxJ[:])
            nyJ = CuArray(nyJ[:])
            sJ  = CuArray(sJ[:])

            ops = (Vq,wq,Qrh_skew,Qsh_skew,Qth,Ph,LIFTq,VPh,Wq)
            mesh = (rxJ,sxJ,ryJ,syJ,nxJ,nyJ,sJ,J,h,mapP,mapP_d)
            param = (K,Np_P,Nq_P,Nfp_P,Nh_P,Np_F,Nq,Nfp,Nh,Nk_max,Nk_tri_max,Nk_fourier_max)

            xq,yq,zq = (x->reshape(x,Nq_P,Np_F,K)).((xq,yq,zq))
            c1 = 0
            c2 = 5.0
            c3 = 4.5
            p0 = 1/gamma
            P_max = 0.4
            r1_exact(x,y,z,t) = zeros(size(x))
            r2_exact(x,y,z,t) = @. -(z-c3-t)
            r3_exact(x,y,z,t) = @. y-c2
            P_exact(x,y,z,t) = P_max*exp.((1 .-r1_exact(x,y,z,t).^2-r2_exact(x,y,z,t).^2-r3_exact(x,y,z,t).^2)/2)
            ρ_exact(x,y,z,t) = (1 .-(gamma-1)/2*P_exact(x,y,z,t).^2).^(1/(gamma-1))
            u_exact(x,y,z,t) = P_exact(x,y,z,t).*r1_exact(x,y,z,t)
            v_exact(x,y,z,t) = P_exact(x,y,z,t).*r2_exact(x,y,z,t)
            w_exact(x,y,z,t) = P_exact(x,y,z,t).*r3_exact(x,y,z,t).+1
            E_exact(x,y,z,t) = p0/(gamma-1)*(1 .-(gamma-1)/2*P_exact(x,y,z,t).^2).^(gamma/(gamma-1)).+ρ_exact(x,y,z,t)/2 .*(u_exact(x,y,z,t).^2 .+v_exact(x,y,z,t).^2 .+w_exact(x,y,z,t).^2)
            p_exact(x,y,z,t) = pfun(ρ_exact(x,y,z,t),(ρ_exact(x,y,z,t).*u_exact(x,y,z,t),ρ_exact(x,y,z,t).*v_exact(x,y,z,t),ρ_exact(x,y,z,t).*w_exact(x,y,z,t)),E_exact(x,y,z,t))
            ρ = ρ_exact(xq,yq,zq,0)
            u = u_exact(xq,yq,zq,0)
            v = v_exact(xq,yq,zq,0)
            w = w_exact(xq,yq,zq,0)
            E = E_exact(xq,yq,zq,0)
            p = p_exact(xq,yq,zq,0)

            Q_exact(x,y,z,t) = (ρ_exact(x,y,z,t),u_exact(x,y,z,t),v_exact(x,y,z,t),w_exact(x,y,z,t),p_exact(x,y,z,t))

            Q = primitive_to_conservative(ρ,u,v,w,p)
            Q = collect(Q)
            Q_vec = zeros(Nq_P,Np_F,Nd,K)
            Q_ex_vec = zeros(Nq_P,Np_F,Nd,K)
            # TODO: clean up
            rq2,sq2,wq2 = quad_nodes_2D(N_P+2)
            Vq2 = vandermonde_2D(N_P,rq2,sq2)/VDM
            xq2,yq2,zq2 = (x->Vq2*reshape(x,Np_P,Np_F*K)).((x,y,z))
            Q_ex = Q_exact(xq2,yq2,zq2,T)

            for k = 1:K
                for d = 1:5
                    @. Q_vec[:,:,d,k] = Q[d][:,:,k]
                end
            end
            Q = Q_vec[:]
            Q = convert.(NUM_TYPE,Q)

            Q = CuArray(Q)
            resQ = CUDA.fill(0.0f0,Nq*Nd*K)

            println("Start timestepping")
            ################################
            ######   Time stepping   #######
            ################################
            # ===================== #
            # ======= RK45 ======== #
            # ===================== #
            @time begin

            @inbounds begin
            for i = 1:Nsteps
                rhstest = 0.0

                for INTRK = 1:5
                    if enable_test
                        @show "==============="
                        @show INTRK
                        @show "==============="
                    end
                    compute_rhstest = INTRK==5
                    rhsQ,rhstest = rhs(Q,ops,mesh,param,num_threads,compute_rhstest,enable_test)
                    resQ .= rk4a[INTRK]*resQ+dt*rhsQ
                    Q .= Q + rk4b[INTRK]*resQ
                end

                if i%10==0 || i==Nsteps
                    println("Time step: $i out of $Nsteps with rhstest = $rhstest")
                end
            end
            end # end @inbounds

            end # end @time

            # Calculate L2 error
            Q = Array(Q)
            Q = reshape(Q,Nq_P,Np_F,Nd,K)
            Q_mat = [zeros(Nq_P,Np_F,K),zeros(Nq_P,Np_F,K),zeros(Nq_P,Np_F,K),zeros(Nq_P,Np_F,K),zeros(Nq_P,Np_F,K)]
            for k = 1:K
                for d = 1:5
                    Q_mat[d][:,:,k] = Q[:,:,d,k]
                end
            end
            Pq = Array(Pq)
            Q = (x->Vq2*Pq*reshape(x,Nq_P,Np_F*K)).(Q_mat)
            p = pfun(Q[1],(Q[2],Q[3],Q[4]),Q[5])
            Q = (Q[1],Q[2]./Q[1],Q[3]./Q[1],Q[4]./Q[1],p)
            Q_ex = (x->reshape(x,length(rq2),Np_F*K)).(Q_ex)
            L2_err = 0.0
            for fld in 1:Nd
                L2_err += sum(h*J*wq2.*(Q[fld]-Q_ex[fld]).^2)
            end
            println("L2err at final time T = $T is $L2_err\n")
            errArr[K1D_arr_idx,Np_F_idx,N_P_idx] = L2_err
        end
    end
end

open("array_conv.txt","w") do io
    writedlm(io,Array(errArr))
end

max_stars_repo_name: yiminllin/ESDG-Fourier

max_issues_repo_name: yiminllin/ESDG-Fourier

max_forks_repo_name: yiminllin/ESDG-Fourier


text:
"""
Gas solubility.

* [`solub_co2`](@ref)
* [`solub_ocs`](@ref)
"""
module Solub

include("../docstring_style.jl")

using Canopy.Constants: atm
using Canopy.Air: air_molar

export solub_co2, solub_ocs

"""
Calculate CO2 solubility in natural water from temperature [K] and salinity [g
kg^-1 seawater].

If the keyword argument `bunsen` is true (default), return the dimensionless
Bunsen solubility coefficient, otherwise return the Henry solubility
coefficient defined via concentration and partial pressure [mol L^-1 atm^-1].

# References

* Murray, C. N., & Riley, J. P. (1971). The solubility of gases in distilled
  water and sea water--IV. Carbon dioxide. _Deep Sea Research_, _18_(5),
  533--541. <https://doi.org/10.1016/0011-7471(71)90077-5>
* Weiss, R. F. (1974). Carbon dioxide in water and seawater: The solubility of
  a non-ideal gas. _Marine Chemistry_, _2_(3), 203--215.
  <https://doi.org/10.1016/0304-4203(74)90015-2>

# Examples

```jldoctest
julia> solub_co2(298.15)
0.8310046248555008
```
"""
function solub_co2(temp, salinity = 0.0; bunsen::Bool = true)
    t = temp * 1e-2  # a transferred scale
    # Henry solubility [mol L^-1 atm^-1]
    kcp_co2 = exp(
        -58.0931 +
        90.5069 / t +
        22.2940 * log(t) +
        salinity * (0.027766 - 0.025888 * t + 0.0050578 * t * t),
    )
    return bunsen ? kcp_co2 * 1e3 / air_molar(temp, atm) : kcp_co2
end

"""
Calculate carbonyl sulfide (OCS) solubility in pure water at a given
temperature [K].

If the keyword argument `bunsen` is true (default), return the dimensionless
Bunsen solubility coefficient, otherwise return the Henry solubility
coefficient defined via concentration and partial pressure [mol L^-1 atm^-1].

# References

* Elliott, S., Lu, E., & Rowland, F. S. (1989). Rates and mechanisms for the
  hydrolysis of carbonyl sulfide in natural waters. _Environmental Science &
  Technology_, _23_(4), 458--461. <https://doi.org/10.1021/es00181a011>
* Sun, W., Maseyk, K., Lett, C., & Seibt, U. (2015). A soil diffusion--reaction
  model for surface COS flux: COSSM v1. _Geoscientific Model Development_,
  _8_(10), 3055--3070. <https://doi.org/10.5194/gmd-8-3055-2015>

# Examples

```jldoctest
julia> solub_ocs(298.15)
0.48759826544380036

julia> solub_ocs(298.15; bunsen=false)
0.019930113211735698
```
"""
function solub_ocs(temp; bunsen::Bool = true)
    # Bunsen solubility [dimensionless]
    k_cos = temp * exp(4050.32 / temp - 20.0007)
    return bunsen ? k_cos : k_cos * air_molar(temp, atm) * 1e-3
end

solub_gas(temp) = error("Not implemented!")

end  # module


text:
using OpenTelemetry
using BenchmarkTools
using Logging

Logging.disable_logging(Logging.Error)

suite = BenchmarkGroup()

#####

trace_suite = BenchmarkGroup()
suite["Trace"] = trace_suite

const TRACER = Tracer(
    provider = TracerProvider(;
        sampler = DEFAULT_ON,
        resource = Resource(
            attributes = StaticAttrs((;
                Symbol("service.name") => "A123456789",
                Symbol("service.version") => "1.34567890",
                Symbol("service.instance.id") => "123ab456-a123-12ab-12ab-12340a1abc12",
            )),
        ),
    ),
)

function create_span(t)
    with_span(
        Span(
            "benchmarkedSpan",
            t;
            attributes = Dict{String,TAttrVal}("long.attribute" => -10000000001000000000),
        ),
    ) do
        push!(current_span(), OpenTelemetry.Event(name = "benchmarkEvent"))
    end
end

trace_suite["Create Span"] = @benchmarkable create_span($TRACER)

#####

meter_suite = BenchmarkGroup()
suite["Meter"] = meter_suite

const PROVIDER = MeterProvider(;
    views = [
        View("benchmark_counter"),
        View(
            "benchmark_histogram";
            aggregation = HistogramAgg{Int}(; boundaries = Tuple(100.0:100.0:900.0)),
        ),
    ],
)
const METER = Meter("benchmark_meter"; provider = PROVIDER)
const COUNTER = Counter{Int}("benchmark_counter", METER)

function update_counter(c, n)
    c(n; name = "a", code = 2, msg = "hi")
end

meter_suite["Update Counter"] =
    @benchmarkable update_counter(c, n) setup = (c = $COUNTER; n = rand(1:100))

const HISTOGRAM = Histogram{Int}("benchmark_histogram", METER)

function update_histogram(h, n)
    h(n; name = "a", code = 2, msg = "hi")
end

meter_suite["Update Histogram"] =
    @benchmarkable update_histogram(h, n) setup = (h = $HISTOGRAM; n = rand(1:1_000))

#####

tune!(suite)
results = run(suite)
println(results)
BenchmarkTools.save(joinpath(@__DIR__, "output.json"), median(results))


text:
module HWunconstrained

	# imports: which packages are we going to use in this module?
	using Distributions, Optim, Plots, DataFrames
	using Random
	using Statistics
	using LinearAlgebra

    export maximize_like_grad, makeData, loglik, grad!, maximize_like, plotLike, plotGrad

	# methods/functions
	# -----------------

	# data creator
	# should/could return a dict with beta,numobs,X,y,norm)
	# true coeff vector, number of obs, data matrix X, response vector y, and a type of parametric distribution for G.
	function makeData(n=10_000)
	           beta = [ 1; 1.5; -0.5 ]
	           Random.seed!(54321)
	           numobs = n
	           X = randn(numobs,3)    # n,k
	           epsilon = randn(numobs)
	           Y = X * beta + epsilon
	           y = 1.0 * (Y .> 0)
	           norm = Normal(0,1)    # create a normal distribution object with mean 0 and variance 1
	           return Dict("beta"=>beta,"n"=>numobs,"X"=>X,"y"=>y,"dist"=>norm)
	       end

	function loglik(betas::Vector, d::Dict)
		lh = log.(cdf.(d["dist"],d["X"]*betas) .^ d["y"] .* (1 .- cdf.(d["dist"],d["X"]*betas)) .^ (1 .- d["y"]))
		return (-1) * sum(lh) / d["n"]
	end

	# gradient of the likelihood at x
	function grad!(storage::Vector,betas::Vector,d)

		for j=1:3
			temp = zeros(length(d["y"]))
 			for i=1:d["n"]
				temp[i] = (d["y"][i] - cdf.(d["dist"], d["X"][i,:]' * betas)) * d["X"][i,j]
			end
			storage[j] = sum(temp)
		end
	end

	# hessian of the likelihood at x
	function hessian!(storage::Matrix,betas::Vector,d)

		#call the gradient we just made

		# Compute the hessian
		.*(storage, storage')
		#change it to matrix

	end

	function info_mat(betas::Vector,d)

	end

	function inv_Info(betas::Vector,d)

	end

	"""
	inverse of observed information matrix
	"""
	function inv_observedInfo(betas::Vector,d)

	end

	"""
	standard errors
	"""
	function se(betas::Vector,d::Dict)

	end

	# function that maximizes the log likelihood without the gradient
	# with a call to `optimize` and returns the result

# Check i( is intentional that the maximization is not so good)
	function maximize_like(x0=[0.8, 1.0, -1.0],meth=NelderMead())
		maximum_like(betas) = loglik(betas, makeData())
		#maximum_like(betas) = (-1) * sum(log.(cdf.(makeData()["dist"], (makeData()["X"] * betas)) ^makeData()["y"]) .+
		#		(log.(1 .- cdf.(makeData()["dist"], (makeData()["X"] * betas)) ^(1 .- makeData()["y"]))))
		nm = optimize(maximum_like, x0, meth)
		return(nm)
	end

	function maximize_like_helpNM(x0=[ 1; 1.5; -0.5 ],meth=NelderMead())

	end

	# function that maximizes the log likelihood with the gradient
	# with a call to `optimize` and returns the result
	function maximize_like_grad(x0=[0.8,1.0,-0.1],meth=BFGS())

	end

	function maximize_like_grad_hess(x0=[0.8,1.0,-0.1],meth=Newton())

	end

	function maximize_like_grad_se(x0=[0.8,1.0,-0.1],meth=BFGS())

	end

	# visual diagnostics
	# ------------------

	# function that plots the likelihood
	# we are looking for a figure with 3 subplots, where each subplot
	# varies one of the parameters, holding the others fixed at the true value
	# we want to see whether there is a global minimum of the likelihood at the the true value
	function plotLike()

		# You want to call the loglik function for different vectors of beta and keep loglik in mind
		# increments = 0.01
		totest_values = collect(-3.0:0.01:3.0)
		#for the first element
		j = zeros(length(totest_values))

		d = makeData()
		b1 = d["beta"][1]
		b2 = d["beta"][2]
		b3 = d["beta"][3]
		for i=1:length(totest_values)
			betas_calc = [totest_values[i], b2, b3]
			j[i] = loglik(betas_calc, d)
		end
		p1 = plot(totest_values, j,
			lab="B2 = 1.5 ; B3 = -0.5",
			xlabel = "B1")

		#for the second element
		j = zeros(length(totest_values))
		for i=1:length(totest_values)
			j[i] = loglik([makeData()["beta"][1], totest_values[i], makeData()["beta"][3]],
					makeData())
		end
		p2 = plot(totest_values, j,
			lab="B1 = 1.0 ; B3 = -0.5",
			xlabel = "B2", ylabel = "Log Likelihood function")

		#for the second element
		j = zeros(length(totest_values))
		for i=1:length(totest_values)
			j[i] = loglik([makeData()["beta"][1], makeData()["beta"][2], totest_values[i]],
					makeData())
		end
		p3 = plot(totest_values, j,
			lab="B1 = 1.0 ; B2 = 1.5",
			xlabel = "B3")

		# Plot final graph
		l = @layout [a; b; c]
		return(plot(p1, p2, p3, layout = l))
	end

	function plotGrad()
			# You want to call the loglik function for different vectors of beta and keep loglik in mind
			# increments = 0.01
			totest_values = collect(-2.0:0.01:2.0)

# Element 1
		# reinitialize storage to true values
			betas = makeData()["beta"]
			storage = zeros(3)
			j = zeros(length(totest_values))
		# run the function grad!() to get the gradient & extract the coef of interest
			for i=1:length(totest_values)
				grad!(storage, [totest_values[i], betas[2], betas[3]], makeData())
				j[i] = storage[1]
			end
			p1 = plot(totest_values, j,
				lab="B2 = 1.5 ; B3 = -0.5",
				xlabel = "B1")

# Element 2
			betas = makeData()["beta"]
			storage = zeros(3)
			j = zeros(length(totest_values))
			for i=1:length(totest_values)
				grad!(storage, [betas[1], totest_values[i], betas[3]], makeData())
				j[i] = storage[2]
			end
			p2 = plot(totest_values, j,
				lab="B1 = 1.0 ; B3 = -0.5",
				xlabel = "B2", ylabel = "Log likelihood gradient")

# Element 3
			betas = makeData()["beta"]
			storage = zeros(3)
			j = zeros(length(totest_values))
			for i=1:length(totest_values)
				grad!(storage, [betas[1], betas[2], totest_values[i]], makeData())
				j[i] = storage[3]
			end
			p3 = plot(totest_values, j,
				lab="B1 = 1.0 ; B2 = 1.5",
				xlabel = "B3")
				l = @layout [a; b; c]
			return(plot(p1, p2, p3, layout = l))
		end

	function saveplots()
		p1 = plotLike()
		savefig(p1,joinpath(dirname(@__FILE__),"..","likelihood.png"))
		p2 = plotGrad()
		savefig(p2,joinpath(dirname(@__FILE__),"..","gradient.png"))
		@info("saved both plots.")
	end

end


text:
const _fdm = central_fdm(5, 1; max_range=1e-2)
const TEST_INFERRED = Ref(true)
const TRANSFORMS_TO_ALT_TANGENTS = Function[x->@thunk(x)] # e.g. [_ -> ZeroTangent(), x -> rebasis(x)]

"sets up TEST_INFERRED based ion enviroment variables"
function init_test_inferred_setting!()
    TEST_INFERRED[] = if haskey(ENV, "CHAINRULES_TEST_INFERRED")
        parse(Bool, "CHAINRULES_TEST_INFERRED")
    else
        !parse(Bool, get(ENV, "JULIA_PKGEVAL", "false"))
    end

    !TEST_INFERRED[] && @warn "inference tests have been disabled"
end


text:
#=
GasLaws.jl:
- Julia version: 1.5.1
- Author: tbald
- Date: 2020-12-20
=#

module GasLaws
using Unitful
using Unitful: 𝐋,𝐍

const R=.0831446261815324u"bar*L*K^-1*mol^-1"
@derived_dimension MolarVolume 𝐋^3/𝐍

#Returns amount in desired unit
function ideal(P::Unitful.Pressure,V::Unitful.Volume,
               n::Unitful.AmountUnits,T::Unitful.Temperature)::Unitful.Amount
    newT= T |> u"K"
    units=Unitful.unit(P)*Unitful.unit(V)/(u"K"*n)
    newR=R |> units
    newN=P*V/(newR*newT)
end

#Returns Temperature in desired unit
function ideal(P::Unitful.Pressure,V::Unitful.Volume,
               n::Unitful.Amount,T::Unitful.TemperatureUnits)::Unitful.Temperature
    units=Unitful.unit(P)*Unitful.unit(V)/(u"K"*Unitful.unit(n))
    newR=R |> units
    newT=P*V/(newR*n)
    newT |> T
end

#Returns Volume in desired unit
function ideal(P::Unitful.Pressure,V::Unitful.VolumeUnits,
               n::Unitful.Amount,T::Unitful.Temperature)::Unitful.Volume
    newT= T |> u"K"
    units=Unitful.unit(P)*V/(u"K"*Unitful.unit(n))
    newR=R |> units
    newV=(n*newR*newT)/P
end

#Returns Pressure in desired unit
function ideal(P::Unitful.PressureUnits,V::Unitful.Volume,
               n::Unitful.Amount,T::Unitful.Temperature)::Unitful.Pressure
    newT= T |> u"K"
    units=P*Unitful.unit(V)/(u"K"*Unitful.unit(n))
    newR=R |> units
    newP=(n*newR*newT)/V
end

#Solves for pressure using Van der Waals
#Default a and b for CO2
function Van_der_Waals(Vm::MolarVolume,T::Unitful.Temperature;a=3.658,b=0.04286)::typeof(1.0u"bar")
    full_a=a*1.0u"bar*L^2*mol^-2"
    full_b=b*1.0u"L*mol^-1"
    newT= T |> u"K"

    p=R*newT/(Vm-full_b)
    p+=-full_a/Vm^2
end

#Solves for pressure using Redlich-Kwong equation
##Default a and b for CO2
function Redlich_Kwong(Vm::MolarVolume,T::Unitful.Temperature;a=64.43,b=0.02963)::typeof(1.0u"bar")
    full_a=a*1.0u"bar*L^2*mol^-2*K^(1/2)"
    full_b=b*1.0u"L*mol^-1"
    newT= T |> u"K"

    p=R*newT/(Vm-full_b)
    p+=-full_a/(newT^(1/2)*Vm*(Vm+full_b))
end

##Could add reduced versions of the real gas laws

end

max_stars_repo_name: TyBalduf/IdealGas

max_issues_repo_name: TyBalduf/IdealGas

max_forks_repo_name: TyBalduf/IdealGas


text:
"""Lateral flow in a HF135 nitrocellulose membrane."""
module ExampleHF135Richards

using Fronts
using Fronts.PorousModels
using Plots

# Wetting of an HF135 membrane, Van Genuchten model
# Data from Buser (PhD thesis, 2016)
# http://hdl.handle.net/1773/38064
k = 5.50e-13  # m^2
α = 0.2555  # 1/m
n = 2.3521

θr = 0.0473
θs = 0.945

hi = -30.591486389337
hb = 0

model = VanGenuchten(n=n, α=α, k=k, θr=θr, θs=θs)

prob = DirichletProblem(RichardsEquation(model), i=hi, b=hb)

h = solve(prob)

r = range(0, 0.05, length=500)

plt = plot(r, r -> h(r,60), label="t = 60 s", xguide="r [m]", yguide="h [m]")
display(plt)

plt = plot(r, r -> θh(model, h(r,60)), label="t = 60 s", xguide="r [m]", yguide="θ")
display(plt)

plt = plot(r[2:end], r -> flux(h, r, 60), label="t = 60 s", xguide="r [m]", yguide="U [m/s]")
display(plt)

end


text:
using Devectorize

function async_gs(m,b,x,r,eps)
###
# Solves the system m * x = b using distributed async Gauss-Seidel method
# m is a distributed matrix
# b,x,r are distributed vectors
# eps is the stopping criterium
###
    res = Inf
    I = localindexes(x)[1]
    localx = convert(Array,x)
    n = size(m,2)

    updatingx = false
    calculatingr = false

    iter = 0

    while res > eps

        # Updating x

        for i in I
            s = 0
            for j = 1:n
                if i != j
                    s += m[i,j] * localx[j]
                end
            end

            localx[i] = (b[i] - s) / m[i,i]
            localpart(x)[i - I.start + 1] = localx[i]
        end

        # Stopping criteria

        res = 0

        for i in I

            s = b[i]

            for j = 1:n

                s -= m[i,j] * localx[j]

            end

            localpart(r)[i - I.start + 1] = abs(s)

            res = max(res,s)

        end

        if res <= eps && !calculatingr

            calculatingr = true
            @async begin
                res = norm(r,Inf)
                calculatingr = false
            end

        end

        if !updatingx

            updatingx = true
            @async begin
                copy!(localx,x)
                updatingx = false
            end

        end

        iter += 1
    end
    return iter
end

function shared_gs(m,b,x,r,eps)
    ###
    # Solves the system m * x = b using a shared memory parallel
    # Gauss-Seidel method
    ###
    res = norm(r,Inf)
    I = localindexes(x)
    n = size(m,2)

    iter = 0
    while res > eps

        # Updating x

        for i in I
            s = 0
            for j = 1:n
                if i != j
                    s += m[i,j] * x[j]
                end
            end

            x[i] = (b[i] - s) / m[i,i]

        end

        # Stopping criterium

        res = 0

        for i in I

            s = b[i]

            for j = 1:n

                s -= m[i,j] * x[j]

            end

            r[i] = abs(s)

            res = max(res,s)

        end

        if res <= eps

            res = norm(r,Inf)

        end

        sleep(rand() / 1000)

        iter += 1
    end
    return iter
end

function serial_gs(m,b,x,r,eps)

    res = Inf
    n = size(m,1)
    iter = 0

    while res > eps
        for i in 1:n
            s = b[i]
            for j = 1:n
                if i != j
                    s -= m[i,j] * x[j]
                end
            end
            x[i] = s / m[i,i]
        end

        res = 0

        for i = 1:n
            s = b[i] - BLAS.dot(n,m[i,:],1,x,1)
#            for j = 1:n
#                s -= m[i,j] * x[j]
#            end
            res = max(res,abs(s))
        end

        iter += 1
    end

    return iter
end

function gs(m,b,x,N)

    for k = 1:N
        I = localindexes(x)[1]
        for i in I
            s = 0
            for j = 1:size(m,2)
                if i != j
                    s += m[i,j] * x[j]
                end
            end
            localpart(x)[i - I.start + 1] = (b[i] - s) / m[i,i]
        end
    end
end

function initializeRandomDistributedDDMatrix(I,l,u)
    m = l + (u - l) * rand(size(I[1],1),size(I[2],1))
    for i in 1:size(I[1],1)
        @devec m[i,I[1][i]] = sum(abs(m[i,:]))
    end
    return m
end

function generateRandomDistributedSystem(n,sol=ones((n,1)),w=workers(),
                                         l=-100,u=100)

    dist = [size(w,1),1]

    M = DArray((I)->initializeRandomDistributedDDMatrix(I,l,u),(n,n),w,dist);

    b = dzeros((size(M,1),1),w,dist)
    r = similar(b)

    for i in w
        @spawnat i localpart(b)[:,1] = localpart(M) * sol
        @spawnat i localpart(r)[:,1] = b[localindexes(b)[1],1]
    end

    x0 = dzeros((size(M,1),1),w,dist)

    return (M,b,x0,r)
end

# Para rodar
# require("gs_distribuido.jl")
#
# (M,b,x,r) = generateRandomDistributedSystem(20)
#
# for i in workers()
#     @spawnat i async_gs(M,b,x,r,0.1)
# end
#
# Ou, a versao de 1 linha (imprime o tempo e o total de iteracoes):
#
# @time reduce(+,pmap(fetch,{(@spawnat i async_gs(M,b,x,r,0.1)) for i in workers()}))

max_stars_repo_name: fsobral/ccs

max_issues_repo_name: fsobral/ccs

max_forks_repo_name: fsobral/ccs


text:
@enum CLBlastLayout CLBlastLayoutRowMajor=101 CLBlastLayoutColMajor=102
@enum CLBlastTranspose CLBlastTransposeNo=111 CLBlastTransposeYes=112 CLBlastTransposeConjugate=113
@enum CLBlastTriangle CLBlastTriangleUpper=121 CLBlastTriangleLower=122
@enum CLBlastDiagonal CLBlastDiagonalNonUnit=131 CLBlastDiagonalUnit=132
@enum CLBlastSide CLBlastSideLeft=141 CLBlastSideRight=142
@enum CLBlastPrecision CLBlastPrecisionHalf=16 CLBlastPrecisionSingle=32 CLBlastPrecisionDouble=64 CLBlastPrecisionComplexSingle=3232 CLBlastPrecisionComplexDouble=6464