text:
using Rayuela, HDF5

function get_dataset_numbers(dataset_name::String)
  nquery, nbase, ntrain = 0, 0, 0
  if dataset_name == "MNIST"
    nquery, nbase, ntrain = Int(10e3), Int(60e3), Int(60e3)
  elseif dataset_name == "labelme"
    nquery, nbase, ntrain = Int(2e3), Int(20019), Int(20019)
  elseif dataset_name in ["SIFT1M", "Convnet1M"]
    nquery, nbase, ntrain = Int(1e4), Int(1e6), Int(1e5)
  else
    error("dataset unknown")
  end
  return nquery, nbase, ntrain
end

function test_experiment(
  dataset_name::String,
  m::Integer,
  trial::Integer)

  @show dataset_name, m, trial
  bpath = joinpath("/home/julieta/Desktop/CQ/build/results/", lowercase(dataset_name) * "_$m", "trial_$trial")
  h = 256

  # Read the codes
  B = read_cq_bvecs(joinpath(bpath, "B"))
  B = convert(Matrix{Int16}, B)
  B .+= 1

  # Read the codebooks
  K = read_cq_fvecs(joinpath(bpath, "D"))
  C = Vector{Matrix{Cfloat}}(m)
  indices = splitarray(1:(m*h), m)
  for i = 1:m
    C[i] = K[:, indices[i]]
  end

  # Load the data
  nquery, nbase, ntrain = get_dataset_numbers(dataset_name)
  verbose = true
  Xt = Rayuela.read_dataset(dataset_name, ntrain)
  Xb = Rayuela.read_dataset(dataset_name * "_base", nbase)
  Xq = Rayuela.read_dataset(dataset_name * "_query", nquery, verbose)[:,1:nquery]
  gt = Rayuela.read_dataset(dataset_name * "_groundtruth", nquery, verbose)
  if dataset_name == "SIFT1M" || dataset_name == "GIST1M"
    gt = gt .+ 1
  end
  gt = convert( Vector{UInt32}, gt[1,1:nquery] )

  qerr = qerror(Xt, B, C) * nbase
  @printf("Qerror: %e\n", qerr )

  if dataset_name in ["SIFT1M", "Convnet1M"]
    # Encode the base
    B_base = convert(Matrix{Int16}, rand(1:h, m, size(Xb,2)))
    ilsiter, icmiter, randord, npert, cpp, V = 4, 4, false, 0, true, true
    B_base = Rayuela.encoding_icm(Xb, B_base, C, ilsiter, icmiter, randord, npert, cpp, V)

    qerr = qerror(Xb, B_base, C)
    @printf("%e\n", qerr)

    # Put the base codes into B
    B = B_base
  end

  # Fast search
  knn = 1000
  if verbose; print("Querying m=$m ... "); end
  @time dists, idx = linscan_cq(B, Xq, C, knn)
  if verbose; println("done"); end
  rec = eval_recall(gt, idx, knn)

  h5write(joinpath(bpath, "recall.h5"), "recall", rec)
end

function run_experiments(
  dataset_name::String,
  dictionaries_count::Int=8)

  ntrials = 10
  for trial = 1:1

    @show "trial", trial

    a = CQ_parameters()
    a.output_file_prefix = "/home/julieta/Desktop/CQ/build/results/sift_$(dictionaries_count)/trial_$trial/"
    a.dictionaries_count = dictionaries_count
    # a.max_iter=1

    if dataset_name == "SIFT1M"
      # do nothing
    elseif dataset_name == "MNIST"
      a.points_count     = Int(60e3)
      a.space_dimension  = 28 * 28
      a.mu               = 0.00001f0
      a.queries_count    = Int(10e3)
      a.groundtruth_length = 1
      a.output_file_prefix = "/home/julieta/Desktop/CQ/build/results/mnist_$(dictionaries_count)/trial_$trial/"
      a.points_file        = "/home/julieta/Desktop/CQ/build/data/mnist/mnist_learn.fvecs"
      a.queries_file       = "/home/julieta/Desktop/CQ/build/data/mnist/mnist_query.fvecs"
      a.groundtruth_file   = "/home/julieta/Desktop/CQ/build/data/mnist/mnist_groundtruth.ivecs"
    elseif dataset_name == "labelme"
      a.points_count     = 20019
      a.space_dimension  = 512
      a.mu               = 100f0
      a.queries_count    = 2000
      a.groundtruth_length = 1
      a.output_file_prefix = "/home/julieta/Desktop/CQ/build/results/labelme_$(dictionaries_count)/trial_$trial/"
      a.points_file        = "/home/julieta/Desktop/CQ/build/data/labelme/labelme_learn.fvecs"
      a.queries_file       = "/home/julieta/Desktop/CQ/build/data/labelme/labelme_query.fvecs"
      a.groundtruth_file   = "/home/julieta/Desktop/CQ/build/data/labelme/labelme_groundtruth.ivecs"
    else
      error("dataset unknown")
    end

    a.trained_dictionary_file = joinpath(a.output_file_prefix, "D")
    a.trained_binary_codes_file = joinpath(a.output_file_prefix, "B")
    a.output_retrieved_results_file = joinpath(a.output_file_prefix, "recall")

    # Prepare output folder
    if !Base.Filesystem.ispath(a.output_file_prefix)
      mkpath(a.output_file_prefix)
    end

    # Write parameters
    pname = "/home/julieta/Desktop/CQ/build/config.txt"
    dump_CQ_parameters(a, pname)

    # Run the experiment
    cmd = joinpath("/home/julieta/Desktop/CQ/build/CompositeQuantization")
    run(`$cmd $pname`)

  end
end

function dump_dataset(dataset_name, verbose=true)
  nquery, nbase, ntrain = get_dataset_numbers(dataset_name)
  Xt = Rayuela.read_dataset(dataset_name * "_base", ntrain)
  Xb = Rayuela.read_dataset(dataset_name * "_base", nbase)
  Xq = Rayuela.read_dataset(dataset_name * "_query", nquery, verbose)[:,1:nquery]
  gt = Rayuela.read_dataset(dataset_name * "_groundtruth", nquery, verbose)

  bpath = joinpath("/home/julieta/Desktop/CQ/build/data/", lowercase(dataset_name))
  if !Base.Filesystem.ispath(bpath)
    mkpath(bpath)
  end

  fvecs_write(Xt, joinpath(bpath, "$(lowercase(dataset_name))_learn.fvecs"))
  fvecs_write(Xb, joinpath(bpath, "$(lowercase(dataset_name))_base.fvecs"))
  fvecs_write(Xq, joinpath(bpath, "$(lowercase(dataset_name))_query.fvecs"))
  ivecs_write(convert(Matrix{Int32},gt.-=1), joinpath(bpath, "$(lowercase(dataset_name))_groundtruth.ivecs"))
end

# dump_dataset("MNIST")
# dump_dataset("labelme")
# dump_dataset("Convnet1M")

# @show("run_experiments")
# run_experiments("MNIST", 16)
run_experiments("labelme", 8)
# run_experiments("SIFT1M", 16)
# run_experiments("Convnet1M")

for trial = 3
  # test_experiment("labelme", 8, trial)
  # test_experiment("labelme", 16, trial)
  # test_experiment("MNIST", 8, trial)
  # test_experiment("MNIST", 16, trial)

  # test_experiment("SIFT1M", 8, trial)
  # test_experiment("SIFT1M", 16, trial)
end


text:
#!/usr/bin/env julia

"""
# Problem 3: Largest prime factor

The prime factors of 13195 are 5, 7, 13 and 29.

What is the largest prime factor of the number `n`?
"""
module Problem003

include("Primes.jl")

"""
    solve(; n::Integer = 600851475143)

Returns the solution for [`Problem003`](@ref) with the given parameters.

## Preconditions
- `n ≥ 2`
"""
function solve(; n::Integer = 600851475143)
    return last(Primes.prime_factorization(n))[1]
end

end

ext: jl

max_stars_repo_name: curtislb/ProjectEuler

max_issues_repo_name: curtislb/ProjectEuler

max_forks_repo_name: curtislb/ProjectEuler


text:
@testset "Expression" begin
    @test_throws ErrorException("Invalid outputType") Expression("name", "expr", outputType="a")
end


text:
module Nanofactory

function produce(inventory, formulae)
    targets = filter(produceable, inventory)
    while !isempty(targets)
        element, deficit = pop!(targets)
        reagents, produced = formulae[element]
        m = ceil(Int, abs(deficit) / produced)
        for (reagent, qty) in reagents
            current = get(inventory, reagent, 0)
            inventory[reagent] = current - (m * qty)
        end
        inventory[element] += (m * produced)
        targets = filter(produceable, inventory)
    end
end

function produceable(kv)
    return kv.first != "ORE" && kv.second < 0
end

function parse_equations(reactions::Array{String, 1})
    formulae = Dict{String,Pair{Array{Pair{String,Int},1},Int}}()
    for reaction in reactions
        rs, product = split(reaction, "=>")
        qty, element = split(strip(product), " ")
        reagents = parse_reagents(convert(String, rs))
        formulae[convert(String, element)] = Pair(reagents, parse(Int, qty))
    end
    return formulae
end

function parse_reagents(reagents::String)
    acc::Array{Pair{String,Int}, 1} = []
    for reagent in split(strip(reagents), ",")
        qty, element = split(strip(reagent), " ")
        push!(acc, Pair(convert(String, element), parse(Int, qty)))
    end
    return acc
end

end


text:
@testset "utils.jl" begin
    @testset "DiscreteMeasure" begin
        @test_throws ErrorException SOT.DiscreteMeasure(randn(3), rand(2))
    end

    @testset "gradient step" begin
        c(x, y) = abs(x - y)
        τ = rand()
        x = randn()
        xs = randn(100)
        ps = rand(100)
        ps ./= sum(ps)
        ν = SOT.DiscreteMeasure(xs, ps)
        v = zeros(100)

        # unregularized and regularized approach
        for ε in (nothing, abs(randn()))
            # out-of-place method
            z0 = SOT.gradient_step(c, τ, ν, x, ε)

            # in-place method with `v = 0`
            z = zero(z0)
            tmp = similar(z)
            SOT.gradient_step!(c, z, τ, v, ν, x, tmp, ε)
            @test z0 ≈ z0
        end
    end
end


text:
#--
# Threads.nthreads() displays number of available threads
# if ==1, start julia with JULIA_NUM_THREADS=4 julia
#--

#--
# If no graphics (X-windows) required, then ENV["GKSwstype"]="nul"
#--
using FDModBiotCyl
#
using CPUTime
using ProgressMeter
using Plots
using Printf
#

#--Suppressing graphics
#ENV["GKSwstype"]="nul"

@__DIR__
@__FILE__

#================
Drawing a model
================#
function drawmodel(Field_Var,nz,dz,nr,dr,LPML_r,LPML_z)
   println("Drawing a model...")
   zvec=[0:dz:(nz-1)*dz]
   rvec=[0:dr:(nr-1)*dr]
   R=(nr-1)*dr
#   plt1=heatmap(rvec,zvec,Vp,yflip=true,xlabel="R (m)",ylabel="Z (m)",ratio=1)
   plt1=heatmap(rvec,zvec,Field_Var,yflip=true,xlabel="R (m)",ylabel="Z (m)")
   plot!([0;R],ones(2)*(LPML_z-1)*dz,label="",color=:black)
   plot!([0;R],ones(2)*((nz-LPML_z+1)*dz),label="",color=:black)
   plot!(ones(2)*((nr-LPML_r+1)*dr),[0;(nz-1)*dz],label="",color=:black)
   display(plot(plt1))
   tmp_filename=@sprintf "tmp_FDmodel.png"
   tmp_filename_full=string(@__DIR__,"/",tmp_filename)
   png(tmp_filename_full)
   println("--------")
   println("Model config figure is saved in:", tmp_filename_full)
   println("--------")
end

#==========================================
Creating a model for a borehole embedded in
 an homogeneous elastic medium
==========================================#
function makemodel_homogeneous_Elastic()
   #Necessary input matrices to the main loop function are:
   #Rho: Bulk density
   #Rhof: Fluid density
   #M,C,H: Poroelastic moduli
   #G: Formation shear modulus
   #D1,D2: Poroelastic moduli in Ou's formulation
   #Flag_AC,Flag_E: flag specifying acoustic region or elastic region
   #Other input parameters to the main loop function are:
   #ir_wall: grid number in the r direction at which a borehole wall starts (single value)
   #
   # This function creates the above input parameters assuming a borehole embedded in
   # an homogeneous elastic medium

   #Model size
   nr=351 #samples
   nz=1001 #samples
   dr=0.005 #meter
   dz=0.005 #meter

   #==============================================
   Initializing grain/frame parameters
   ==============================================#
   #First creating grain/frame parameters and then converting them
   #to necessary poroelastic parameters for our FD
   #----solid phase-----
   Km=zeros(nz,nr) #frame
   Ks=zeros(nz,nr) #grain
   G=zeros(nz,nr) #bulk
   Rho=zeros(nz,nr) #bulk (weighted average of Rhos and Rhof)
   Rhos=zeros(nz,nr) #grain
   #---fluid phase----
   Kf=zeros(nz,nr)
   Rhof=zeros(nz,nr)
   Kappa0=zeros(nz,nr) #static permeability (m^2, 1D=9.869E-13m^2)
   Eta=zeros(nz,nr) #dynamic viscosity Pa.s=kg/(m.s)
   #---other parameters---
   Phi=zeros(nz,nr)
   Tot=zeros(nz,nr) #Tortuosity factor (see Sidler 2014)
   #End initializing poroelastic parameters

   #==============================================
   Fluid parameters
   ==============================================#
   #Assuming that fluid properties are identical everywhere
   Kf[:,:]=ones(nz,nr)*2.25*10^9
   Rhof[:,:]=ones(nz,nr)*1000
   Eta[:,:]=ones(nz,nr)*1.0*10^(-3) #water 1E-3 Pa.s

   #========================================
   Inclusion of a borehole (acoustic media)
   ========================================#
   Flag_AC=zeros(nz,nr)
   #Borehole wall (homogeneous radius)
   ir_wall=Int(round(0.1/dr))

   for iz=1:nz
       Flag_AC[iz,1:ir_wall]=ones(ir_wall)
       Rho[iz,1:ir_wall]=ones(ir_wall)*Rhof[1,1] #Rho=Rhof
       G[iz,1:ir_wall]=ones(ir_wall)*0.0 #G=0
       Phi[iz,1:ir_wall]=ones(ir_wall) #Phi=1
       Km[iz,1:ir_wall]=ones(ir_wall)*0.0 #Km=0 (G=0,phi=0->M=C=H=Kf, Ou's)
       Ks[iz,1:ir_wall]=ones(ir_wall)*Kf[1,1] #
       Eta[iz,1:ir_wall]=ones(ir_wall)*0.0 #eta=0 (-> D2=0, Ou's)
   end

   #==========================
   Homogeneous Elastic model
   ==========================#
   Flag_E=zeros(nz,nr)
   Vp1_elastic=2500.0 #Or, Specify K_elastic
   Vs1_elastic=2000.0
   Rho1_elastic=2500.0

#
   for iz=1:nz #Elastic media

       G_elastic=Vs1_elastic^2*Rho1_elastic
       K_elastic=Vp1_elastic^2*Rho1_elastic-4/3*G_elastic

       Flag_E[iz,ir_wall+1:end]=ones(nr-ir_wall)
       Rho[iz,ir_wall+1:end]=ones(nr-ir_wall)*Rho1_elastic #Rho=Rho_elastic
       Rhof[iz,ir_wall+1:end]=ones(nr-ir_wall)*Rho1_elastic #Rhof=Rho_elastic
       G[iz,ir_wall+1:end]=ones(nr-ir_wall)*G_elastic #G=Gs (Ou's)
       #--
       Phi[iz,ir_wall+1:end]=ones(nr-ir_wall)*0.0 #Phi=0
       Km[iz,ir_wall+1:end]=ones(nr-ir_wall)*K_elastic #Km
       Ks[iz,ir_wall+1:end]=ones(nr-ir_wall)*K_elastic #Ks=Km (-> M=Inf, Ou's)
       #--
       Kappa0[iz,ir_wall+1:end]=ones(nr-ir_wall)*0.0 #k0=0 (-> D1=D2=Inf, Ou's)
   end

   #=====================================================
   Converting grain/frame parameters into poroelastic parameters
   =====================================================#
   #Rho, Rhof, M, C, H, G, D1, D2
   #Rho as is
   #Rhof as is
   Alpha=ones(nz,nr)-Km./Ks
   M=( (Alpha-Phi)./Ks+Phi./Kf ).^(-1)
   C=M.*Alpha
   H=Km+4/3*G+M.*Alpha.^2
   #G as is
   D1=Tot.*Rhof./Phi #(1+2/m)*alpha_inf*rhof/phi
   D2=Eta./Kappa0

   #Correction of H and C for elastic media
   for ir=1:nr
      for iz=1:nz
         if(Flag_E[iz,ir]==1)
            H[iz,ir]=Ks[iz,ir]+4/3*G[iz,ir]
            C[iz,ir]=Ks[iz,ir] #Ou's
         end

         if(Flag_AC[iz,ir]==1)
#            D1[iz,ir]=1.0 #Guan's
            D1[iz,ir]=1000. #Ou's D1=rhof
            D2[iz,ir]=0. #Ou's D2=0
         end
      end
   end

   return nr,nz,dr,dz,Rho,Rhof,M,C,H,G,D1,D2,Flag_AC,Flag_E,ir_wall
end

#================
Drawing a snapshot
================#
function drawsnap(Data1,nz,dz,nr,dr,LPML_r,LPML_z,title_str)
   println("Drawing a snapshot...")
   zvec=[0:dz:(nz-1)*dz]
   rvec=[0:dr:(nr-1)*dr]
   R=(nr-1)*dr
   plt1=heatmap(rvec,zvec,Data1,yflip=true,xlabel="R (m)",ylabel="Z (m)",legend=:none)
   heatmap!(title=title_str)
   plot!([0;R],ones(2)*(LPML_z-1)*dz,label="",color=:black)
   plot!([0;R],ones(2)*((nz-LPML_z+1)*dz),label="",color=:black)
   plot!(ones(2)*((nr-LPML_r+1)*dr),[0;(nz-1)*dz],label="",color=:black)
   display(plot(plt1))
end

#----MAIN loop function with src----
function main_loop!(nr,nz,dr,dz,Rho,Rhof,M,C,H,G,D1,D2,dt,nt,T,
   vr,vz,trr,tpp,tzz,trz,
   vfr,vfz,pf,
   ir_wall,
   src_index,src_dn,
   Flag_AC,Flag_E,
   nrec,index_allrec_tii,rec_tii,
   LPML_r,LPML_z,PML_Wr,PML_Wz,PML_IWr,PML_Wr2,PML_Wz2,PML_IWr2,
   snapshots_vr,snapshots_vz,snapshots_pf,nsnap,itvec_snap,nskip)

#println("start")

#--
Flag_vf_zero=get_Flag_vf_zero(Flag_AC,Flag_E,nr,nz)

#initializing matrices of PML
Prz_T,Pzz_T,Rz_T,Srz_T,PzzPE_T,RzPE_T,
Prz_B,Pzz_B,Rz_B,Srz_B,PzzPE_B,RzPE_B,
Prr_R,Qrp_R,Pzr_R,Qzp_R,Rr_R,Rp_R,Rrz_R,PrrPE_R,RrPE_R,RpPE_R,
Prr_TR,Qrp_TR,Prz_TR,Pzr_TR,Qzp_TR,Pzz_TR,Rr_TR,Rp_TR,Rz_TR,Rrz_TR,Srz_TR,PrrPE_TR,PzzPE_TR,RrPE_TR,RpPE_TR,RzPE_TR,
Prr_BR,Qrp_BR,Prz_BR,Pzr_BR,Qzp_BR,Pzz_BR,Rr_BR,Rp_BR,Rz_BR,Rrz_BR,Srz_BR,PrrPE_BR,PzzPE_BR,RrPE_BR,RpPE_BR,RzPE_BR=
init_PML_Por(LPML_r,LPML_z,nr,nz)

memT_vr,memT_vz,memT_trr,memT_tpp,memT_tzz,memT_trz,memT_vfr,memT_vfz,memT_pf,
memB_vr,memB_vz,memB_trr,memB_tpp,memB_tzz,memB_trz,memB_vfr,memB_vfz,memB_pf,
memR_vr,memR_vz,memR_trr,memR_tpp,memR_tzz,memR_trz,memR_vfr,memR_vfz,memR_pf=
init_memory_variables_Por(LPML_r,LPML_z,nr,nz)
#return

# Making sure RightBC for stress (zero values are trp and trz)
ApplyBCRight_stress!(vr,vz, #Use with flag_zero=1 (see ApplyBCRight_stress1D01)
    trz,
    G,nr,nz,dr,dz,dt)

@showprogress for ii=1:nt

#----Time at (ii-1)*dt---(updating velocities)-----

#--PML: save velocity at previous step
PML_save_vel!(memT_vr,memT_vz,memT_vfr,memT_vfz,
              memB_vr,memB_vz,memB_vfr,memB_vfz,
              memR_vr,memR_vz,memR_vfr,memR_vfz,
              vr,vz,vfr,vfz,nr,nz,LPML_z,LPML_r)
#println("update vel")

# Main velocity update!
update_velocity_1st_Por!(vr,vz,trr,tpp,tzz,trz,vfr,vfz,pf,Rho,Rhof,D1,D2,Flag_vf_zero,nr,nz,dr,dz,dt,LPML_z,LPML_r)

#PML: update velocity
PML_update_vel!(vr,vz,trr,tpp,tzz,trz,vfr,vfz,pf,
              Rho,Rhof,D1,D2,Flag_vf_zero,nr,nz,dr,dz,dt,
              Prz_T,Pzz_T,PzzPE_T,
              Prz_B,Pzz_B,PzzPE_B,
              Prr_R,Qrp_R,Pzr_R,Qzp_R,Rr_R,Rp_R,Rrz_R,
              PrrPE_R,RrPE_R,RpPE_R,
              Prz_TR,Pzz_TR,PzzPE_TR,
              Prr_TR,Qrp_TR,Pzr_TR,Qzp_TR,Rr_TR,Rp_TR,Rrz_TR,
              PrrPE_TR,RrPE_TR,RpPE_TR,
              Prz_BR,Pzz_BR,PzzPE_BR,
              Prr_BR,Qrp_BR,Pzr_BR,Qzp_BR,Rr_BR,Rp_BR,Rrz_BR,
              PrrPE_BR,RrPE_BR,RpPE_BR,
              LPML_z,LPML_r)
#println("update vel Left")

# Making sure RightBC for velocity (zero values are vr)
ApplyBCRight_vel!(vr,vz,
      trr,tpp,tzz,trz,
      vfr,vfz,pf,
      Rho,Rhof,D1,D2,Flag_vf_zero,nr,nz,dr,dz,dt)

#---velocity B.Cs-----
#Periodic Left Edge
ApplyBCLeft_vel!(vr,vz,trr,tpp,tzz,trz,vfr,vfz,pf,
                          Rho,Rhof,D1,D2,Flag_vf_zero,
                          nr,nz,dr,dz,dt,
                          Prz_T,Pzz_T,PzzPE_T,
                          Prz_B,Pzz_B,PzzPE_B,
                          LPML_z)

#PML: update memory variables for stress (Rx and Sxx) using velocity at two time steps
PML_update_memRS!(Rz_T,Srz_T,RzPE_T,
                  Rz_B,Srz_B,RzPE_B,
                  Rr_R,Rp_R,Rrz_R,RrPE_R,RpPE_R,
                  Rr_TR,Rp_TR,Rz_TR,Rrz_TR,Srz_TR,RrPE_TR,RpPE_TR,RzPE_TR,
                  Rr_BR,Rp_BR,Rz_BR,Rrz_BR,Srz_BR,RrPE_BR,RpPE_BR,RzPE_BR,
                  memT_vr,memT_vz,memT_vfr,memT_vfz,
                  memB_vr,memB_vz,memB_vfr,memB_vfz,
                  memR_vr,memR_vz,memR_vfr,memR_vfz,
                  vr,vz,vfr,vfz,nr,nz,dr,dz,dt,LPML_z,LPML_r,PML_Wz,PML_Wr,PML_IWr,PML_Wz2,PML_Wr2,PML_IWr2)

#--src injection (velocity)
#srcamp=src_func[ii]
#srcapply!(vz,src_index,src_dn,srcamp)
#srcapply!(vr,src_index,src_dn,srcamp)

#----Time at (ii-1)*dt+dt/2----(updating stress)---

#--PML: save stress at previous step
PML_save_stress!(memT_trr,memT_tpp,memT_tzz,memT_trz,memT_pf,
                 memB_trr,memB_tpp,memB_tzz,memB_trz,memB_pf,
                 memR_trr,memR_tpp,memR_tzz,memR_trz,memR_pf,
                 trr,tpp,tzz,trz,pf,nr,nz,LPML_z,LPML_r)

#println("update stress")

# Main stress update!
update_stress_1st_Por!(vr,vz,trr,tpp,tzz,trz,vfr,vfz,pf,
     M,C,H,G,nr,nz,dr,dz,dt,LPML_z,LPML_r)
#error()

#PML: stress update
PML_update_stress!(vr,vz,trr,tpp,tzz,trz,vfr,vfz,pf,
                  M,C,H,G,nr,nz,dr,dz,dt,
                  Rz_T,Srz_T,RzPE_T,
                  Rz_B,Srz_B,RzPE_B,
                  Rr_R,Rp_R,Rrz_R,RrPE_R,RpPE_R,
                  Rr_TR,Rp_TR,Rz_TR,Rrz_TR,Srz_TR,RrPE_TR,RpPE_TR,RzPE_TR,
                  Rr_BR,Rp_BR,Rz_BR,Rrz_BR,Srz_BR,RrPE_BR,RpPE_BR,RzPE_BR,
                  LPML_z,LPML_r)

# Making sure RightBC for stress (zero values are trp and trz)
ApplyBCRight_stress!(vr,vz, #Use with flag_zero=1 (see ApplyBCRight_stress1D01)
     trz,
     G,nr,nz,dr,dz,dt)

#---stress B.Cs
#println("update stress Left")
ApplyBCLeft_stress!(vr,vz,trr,tpp,tzz,trz,vfr,vfz,pf,
                    M,C,H,G,nr,nz,dr,dz,dt,
                    Rz_T,Srz_T,RzPE_T,
                    Rz_B,Srz_B,RzPE_B,
                    LPML_z)

#--src injection (stress:monopole src)
srcamp=-src_func[ii]
srcamp=srcamp*dt/Rhof[1,1]
srcapply!(tzz,src_index,src_dn,srcamp)
srcapply!(trr,src_index,src_dn,srcamp)
srcapply!(tpp,src_index,src_dn,srcamp)

#--PML: update memory variables for velocity (Pxx and Qxx) using stress at two time steps
PML_update_memPQ!(Prz_T,Pzz_T,PzzPE_T,
                  Prz_B,Pzz_B,PzzPE_B,
                  Prr_R,Qrp_R,Pzr_R,Qzp_R,PrrPE_R,
                  Prr_TR,Qrp_TR,Prz_TR,Pzr_TR,Qzp_TR,Pzz_TR,PrrPE_TR,PzzPE_TR,
                  Prr_BR,Qrp_BR,Prz_BR,Pzr_BR,Qzp_BR,Pzz_BR,PrrPE_BR,PzzPE_BR,
                  memT_trr,memT_tpp,memT_tzz,memT_trz,memT_pf,
                  memB_trr,memB_tpp,memB_tzz,memB_trz,memB_pf,
                  memR_trr,memR_tpp,memR_tzz,memR_trz,memR_pf,
                  trr,tpp,tzz,trz,pf,nr,nz,dr,dz,dt,LPML_z,LPML_r,PML_Wz,PML_Wr,PML_IWr,PML_Wz2,PML_Wr2,PML_IWr2)

#-----End of updating field----

#println("receiver")
#Receiver field (extraction)
#getRecData_from_index!(vr,rec_vr,index_allrec_vr,nrec,ii)
#getRecData_from_index!(vz,rec_vz,index_allrec_vz,nrec,ii)
getRecData_from_index!(trr,rec_tii,index_allrec_tii,nrec,ii)

#--Snapshots
#println("check snap")
check_snap=findall(x -> x==ii,itvec_snap)
if (length(check_snap)!=0)
   println("ii=",ii)
   cnt_snap=Int(check_snap[1])
   get_snapshots!(snapshots_vz,cnt_snap,vz,nr,nz)
   get_snapshots!(snapshots_vr,cnt_snap,vr,nr,nz)
   get_snapshots!(snapshots_pf,cnt_snap,pf,nr,nz)

   drawsnap(snapshots_vz[:,:,cnt_snap],nz,dz,nr,dr,
            LPML_r,LPML_z,"Snapshots: Vz")

end

end #i Time

end #function
#------------------------------

#-----ENTRY POINT HERE----
CPUtic()
start=time()

#======================
General variables
======================#

#time samples
dt=1E-6
nt=2501
T=(nt-1)*dt
tvec=range(0.0,T,length=nt) # range object (no memory allocation)
tvec=collect(tvec) # a vector

#PML thickness in samples
LPML_r=60
LPML_z=60

#======================
Creating a model
======================#
nr,nz,dr,dz,Rho,Rhof,M,C,H,G,D1,D2,Flag_AC,Flag_E,ir_wall=makemodel_homogeneous_Elastic()
drawmodel(H,nz,dz,nr,dr,LPML_r,LPML_z)
#error()

#======================
Stability check
======================#
Vmax=maximum((H./Rho).^(0.5))
check_stability01(dt,dr,Vmax,0)
check_stability01(dt,dz,Vmax,0)
#error()
#======================
Creating a src wavelet
======================#
f0=5000 #src Freq
delay=1/f0*1.0
src_func=myricker2(tvec,f0,delay,1) #when using 1st derivative Gaussian
tmp_maxamp=maximum(map(abs,src_func))
src_func=src_func/tmp_maxamp*1E+3

display(plot(tvec,src_func[:],title="src function"))
tmpfilename_src=string(@__DIR__,"/src_function.png")
png(tmpfilename_src)
println("Source signature figure saved in ",tmpfilename_src)

#================================
Point src geometry
================================#
srcgeom=zeros(1,2) #(z,r)
srcgeom[1,1]=(nz-1)*dz/3 #z meter
srcgeom[1,2]=0*dr #r meter
#src_index,src_dn=get_srcindex_monopole(srcgeom,dr,dz)
#Gaussian point src
wsize=5 #window size (odd number)
wsigma=1 #std
src_index,src_dn=get_srcindex_pGauss(srcgeom,dr,dz,wsize,wsigma)

#error()
#==============================
Initializing field variables
==============================#
vr,vz,trr,tpp,tzz,trz,vfr,vfz,pf=init_fields_Por(nz,nr) #

# Initializing PML profile
LPML_r,LPML_z,PML_Wr,PML_Wz,
PML_IWr,PML_Wr2,PML_Wz2,PML_IWr2=init_PML_profile(LPML_r,LPML_z,Vmax,dr,dz,nr,f0)
#PML_check(LPML_r,LPML_z,Vmax,dr,dz,f0)

#===================
Receiver geometry
====================#
#--Fluid pressure at borehole center
nrec=nz
recgeom=zeros(nrec,2) #(z,r)
for irec=1:nrec
   recgeom[irec,1]=dz*(irec-1)
   recgeom[irec,2]=0.0
end
rec_tii,index_allrec_tii=init_receiver_hydrophone(recgeom,nrec,dr,dz,nr,nz,nt)

#==============================
Snapshot settings
==============================#
nskip=100 #change here to adjust time sampling (1~nt)
snapshots_vz,nsnap,itvec_snap=init_snapshots(nt,nz,nr,nskip)
snapshots_vr,nsnap,itvec_snap=init_snapshots(nt,nz,nr,nskip)
snapshots_pf,nsnap,itvec_snap=init_snapshots(nt,nz,nr,nskip)

#error()

#==============================
Start main FD Loop
==============================#

main_loop!(nr,nz,dr,dz,Rho,Rhof,M,C,H,G,D1,D2,dt,nt,T,
   vr,vz,trr,tpp,tzz,trz,
   vfr,vfz,pf,
   ir_wall,
   src_index,src_dn,
   Flag_AC,Flag_E,
   nrec,index_allrec_tii,rec_tii,
   LPML_r,LPML_z,PML_Wr,PML_Wz,PML_IWr,PML_Wr2,PML_Wz2,PML_IWr2,
   snapshots_vr,snapshots_vz,snapshots_pf,nsnap,itvec_snap,nskip)

CPUtoc()
println("elapsed real time: ", round(time() - start;digits=3)," seconds")

println("FD finished without a problem.")

println("Plotting receiver response and comparing with a result using DWI...")
offset=recgeom[:,1]-srcgeom[1]*ones(size(recgeom[:,1]))
plt1=heatmap(offset[:],tvec[:],-rec_tii,xlabel="Offset (m)",
            ylabel="time (s)",title="fluid pressure",clim=(-0.01,0.01))
heatmap!(flip=true)
#Loading data pre-caluculated by DWI
ns_dwi=1024
tvec_dwi=collect(range(0,0.003,length=ns_dwi))
data_dwi = Array{Float32}(undef, ns_dwi, 1); #
#read!("/Users/mminato/work/Git_Projects/FDModBiotCyl.jl/trunk/DWI/DWI/work/dwi_2m_5kHz.bin", data_dwi)
dwifilepath=joinpath(@__DIR__,"dwi_2m_5kHz.bin")
read!(dwifilepath, data_dwi)
plt2=plot(tvec,-rec_tii[:,735]/maximum(rec_tii[:,735]),xlabel="time (s)",label="FD")
plot!(tvec_dwi+ones(size(tvec_dwi))*0.000065,-data_dwi/maximum(data_dwi),xlabel="time (s)",
      label="DWI",title="offset@2m")
display(plot(plt1,plt2,layout=(2,1)))


text:
using SortingNetworks
using Test

function bittuple_impl(::Val{N}) where {N}
    args = [:(ithbit(x, $i)) for i in 0:N-1]
   :(tuple($(args...)))
end
ithbit(x, i) = Bool((x >> i) & true)
@generated bittuple(x, N) = bittuple_impl(N())

function bittuples(vN::Val{N}) where {N}
    @assert N <= 25
    xs = UInt64(1):UInt64(2)^N
    [bittuple(x, vN) for x in xs]
end

@noinline function test_swapsort(vN::Val)
    for cba ∈ bittuples(vN)
        abc = swapsort(cba)
        @test issorted(abc)
    end
end

test_swapsort(n) = test_swapsort(Val{n}())

@testset "test 0-1 principle implementation" begin
    N = 8
    tups = @inferred bittuples(Val{N}())
    @test length(tups) == 2^N
    @test length(unique(tups)) == length(tups)
end

@testset "0-1 principle" for n ∈ 1:25
    test_swapsort(n)
end


text:
using .CuArrays.CUDNN: @check, libcudnn, cudnnStatus_t, libcudnn_handle,
  cudnnDataType, TensorDesc, FilterDesc

using LinearAlgebra

mutable struct DropoutDesc
  ptr::Ptr{Nothing}
  states::CuVector{UInt8}
end

Base.unsafe_convert(::Type{Ptr{Nothing}}, dd::DropoutDesc) = dd.ptr

function DropoutDesc(ρ::Real; seed::Integer=0)
  d = [C_NULL]
  s = Csize_t[0]
  @check ccall((:cudnnCreateDropoutDescriptor,libcudnn), cudnnStatus_t, (Ptr{Ptr{Nothing}},), d)
  @check ccall((:cudnnDropoutGetStatesSize,libcudnn),cudnnStatus_t,(Ptr{Nothing},Ptr{Csize_t}),libcudnn_handle[],s)
  states = CuArray{UInt8}(s[]) # TODO: can we drop this when ρ=0?
  desc = DropoutDesc(d[], states)
  @check ccall((:cudnnSetDropoutDescriptor,libcudnn),cudnnStatus_t,(Ptr{Nothing},Ptr{Nothing},Cfloat,Ptr{Nothing},Csize_t,Culonglong),
    desc,libcudnn_handle[],ρ,states,length(states),seed)
  finalizer(desc) do x
    @check ccall((:cudnnDestroyDropoutDescriptor,libcudnn),cudnnStatus_t,(Ptr{Nothing},),x)
  end
  return desc
end

const RNN_RELU = 0 # Stock RNN with ReLu activation
const RNN_TANH = 1 # Stock RNN with tanh activation
const LSTM = 2     # LSTM with no peephole connections
const GRU = 3      # Using h' = tanh(r * Uh(t-1) + Wx) and h = (1 - z) * h' + z * h(t-1)

const LINEAR_INPUT = 0
const SKIP_INPUT = 1

const UNIDIRECTIONAL = 0
const BIDIRECTIONAL = 1

const RNN_ALGO_STANDARD = 0
const RNN_ALGO_PERSIST_STATIC = 1
const RNN_ALGO_PERSIST_DYNAMIC = 2

# param layout:
# RNN: [weight, bias] × [input, hidden]
# GRU: [weight, bias] × [input, hidden] × [reset, update, newmem]
# LSTM: [weight, bias] × [input, hidden] × [input, forget, newmem, output]

function params(w::CuVector, input, hidden, n = 1)
  slice(offset, shape) = reshape(view(w, offset.+(1:prod(shape))), shape)
  wx = slice(0, (input, hidden*n))
  wh = slice(length(wx), (hidden, hidden*n))
  bias = view(w, length(wx)+length(wh) .+ (1:hidden*n))
  (wx, wh), bias
end

mutable struct RNNDesc{T}
  mode::Int
  input::Int
  hidden::Int
  params::CuVector{T}
  weights::NTuple{2,CuMatrix{T}}
  bias::CuVector{T}
  ptr::Ptr{Nothing}
end

Base.unsafe_convert(::Type{Ptr{Nothing}}, d::RNNDesc) = d.ptr

function rnnParamSize(T, r, input)
  size = Csize_t[0]
  @check ccall((:cudnnGetRNNParamsSize, libcudnn), cudnnStatus_t, (Ptr{Nothing},Ptr{Nothing},Ptr{Nothing},Ptr{Csize_t},Cint),
    libcudnn_handle[], r, TensorDesc(T, (1,input,1)), size, cudnnDataType(T))
  return Int(size[])÷sizeof(T)
end

ngates(mode) = [1, 1, 4, 3][mode+1]
ngates(r::RNNDesc) = ngates(r.mode)

function RNNDesc{T}(mode::Int, input::Int, hidden::Int; layers = 1) where T
  d = [C_NULL]
  @check ccall((:cudnnCreateRNNDescriptor,libcudnn),cudnnStatus_t,(Ptr{Ptr{Nothing}},),d)

  dropoutDesc = DropoutDesc(0)
  inputMode = LINEAR_INPUT
  direction = UNIDIRECTIONAL
  algo = RNN_ALGO_STANDARD
  @check ccall((:cudnnSetRNNDescriptor_v6,libcudnn), cudnnStatus_t, (Ptr{Nothing},Ptr{Nothing},Cint,Cint,Ptr{Nothing},Cint,Cint,Cint,Cint,Cint),
    libcudnn_handle[],d[],hidden,layers,dropoutDesc,inputMode,direction,mode,algo,cudnnDataType(T))

  w = cuzeros(T, rnnParamSize(T, d[], input))
  # TODO: avoid reserve allocation here
  rd = RNNDesc{T}(mode, input, hidden, w, params(w, input, hidden, ngates(mode))..., d[])
  finalizer(rd) do x
    @check ccall((:cudnnDestroyRNNDescriptor,libcudnn),cudnnStatus_t,(Ptr{Nothing},),x)
  end
  return rd
end

function rnnWorkspaceSize(r::RNNDesc, seqlen, xdesc)
  size = Csize_t[0]
  @check ccall((:cudnnGetRNNWorkspaceSize, libcudnn), cudnnStatus_t, (Ptr{Nothing},Ptr{Nothing},Cint,Ptr{Ptr{Nothing}},Ptr{Csize_t}),
    libcudnn_handle[], r, seqlen, xdesc, size)
  return Int(size[])
end

const workspace = [CuVector{UInt8}(1)]

getworkspace(bytes) =
  length(workspace[]) ≥ bytes ?
    workspace[] :
    (workspace[] = CuVector{UInt8}(bytes))

getworkspace(r::RNNDesc, seqlen, xdesc) =
  getworkspace(rnnWorkspaceSize(r, seqlen, xdesc))

function rnnTrainingReserveSize(r::RNNDesc, seqlen, xdesc)
  size = Csize_t[0]
  @check ccall((:cudnnGetRNNTrainingReserveSize,libcudnn), cudnnStatus_t, (Ptr{Nothing}, Ptr{Nothing}, Cint, Ptr{Ptr{Nothing}}, Ptr{Csize_t}),
    libcudnn_handle[], r, seqlen, xdesc, size)
  return Int(size[])
end

function cudnnRNNForward(rnn::RNNDesc{T}, seqlen, xd, x, hd, h, cd, c, wd, w, yd, y, hod, ho, cod, co,
                         workspace, reserve=nothing) where T
  if reserve == nothing
    @check ccall((:cudnnRNNForwardInference, libcudnn), cudnnStatus_t,
                 (Ptr{Nothing}, Ptr{Nothing}, Cint,
                  Ptr{Ptr{Nothing}}, Ptr{T}, Ptr{Nothing}, Ptr{T}, Ptr{Nothing}, Ptr{T},
                  Ptr{Nothing}, Ptr{T}, Ptr{Ptr{Nothing}}, Ptr{T}, Ptr{Nothing}, Ptr{T},
                  Ptr{Nothing}, Ptr{T},
                  Ptr{Nothing}, Csize_t),
                 libcudnn_handle[], rnn, seqlen,
                 xd, x, hd, h, cd, c, wd, w, yd, y, hod, ho, cod, co,
                 workspace, length(workspace))
  else
    @check ccall((:cudnnRNNForwardTraining, libcudnn), cudnnStatus_t,
                 (Ptr{Nothing}, Ptr{Nothing}, Cint,
                  Ptr{Ptr{Nothing}}, Ptr{T}, Ptr{Nothing}, Ptr{T}, Ptr{Nothing}, Ptr{T}, Ptr{Nothing}, Ptr{T}, Ptr{Ptr{Nothing}}, Ptr{T}, Ptr{Nothing}, Ptr{T}, Ptr{Nothing}, Ptr{T},
                  Ptr{Nothing}, Csize_t, Ptr{Nothing}, Csize_t),
                 libcudnn_handle[], rnn, seqlen,
                 xd, x, hd, h, cd, c, wd, w, yd, y, hod, ho, cod, co,
                 workspace, length(workspace), reserve, length(reserve))
  end
end

xDesc(x) = [TensorDesc(eltype(x), (1, size(x, 1), size(x, 2)))]

hDesc(h::Nothing) = C_NULL, C_NULL
hDesc(x::Integer) = (@assert x == 0; hDesc(nothing))
function hDesc(h::CuArray)
  TensorDesc(eltype(h), (size(h, 1), size(h, 2), 1)), h
end

# TODO: can we just manipulate strides here?
# TODO: should use repmat, but this isn't implemented.
hBatch(x::AbstractVector, h::CuVector) = h
hBatch(x::AbstractMatrix, h::CuVector) = h .* cuones(1, size(x, 2))
hBatch(x::AbstractMatrix, h::CuMatrix) = h .* cuones(1, size(h,2) == 1 ? size(x,2) : 1)

function forward(rnn::RNNDesc{T}, x::CuArray{T}, h_::CuArray{T}, c_ = nothing, train = Val{false}) where T
  h = hBatch(x, h_)
  c = c_ == nothing ? nothing : hBatch(x, c_)
  @assert size(x, 1) == rnn.input
  @assert size(h, 1) == rnn.hidden
  @assert size(x, 2) == size(h, 2)
  seqLength = 1
  xdesc = xDesc(x)
  y = x isa AbstractVector ? similar(x, rnn.hidden) : similar(x, rnn.hidden, size(x, 2))
  ho = similar(h)
  ydesc = xDesc(y)
  workspace = getworkspace(rnn, seqLength, xdesc)
  reserve = train == Val{true} ?
    CuVector{UInt8}(rnnTrainingReserveSize(rnn, seqLength, xdesc)) :
    nothing
  co = c == nothing ? c : similar(c)
  cudnnRNNForward(rnn, seqLength,
                  xdesc, x,
                  hDesc(h)...,
                  hDesc(c)...,
                  FilterDesc(T, (1, 1, length(rnn.params))), rnn.params,
                  ydesc, y,
                  hDesc(ho)...,
                  hDesc(co)...,
                  workspace, reserve)
  result = c == nothing ? (y, ho) : (y, ho, co)
  return train == Val{true} ? (reserve, result) : result
end

forwardTrain(rnn::RNNDesc{T}, x::CuArray{T}, h::CuArray{T}, c = nothing) where T =
  forward(rnn, x, h, c, Val{true})

function cudnnRNNBackwardData(rnn::RNNDesc{T}, seqlen, yd, y, dyd, dy, dhod, dho, dcod, dco,
                              wd, w, hd, h, cd, c, dxd, dx, dhd, dh, dcd, dc, ws, rs) where T
  @check ccall((:cudnnRNNBackwardData,libcudnn),cudnnStatus_t,
               (Ptr{Nothing}, Ptr{Nothing}, Cint,
                Ptr{Ptr{Nothing}}, Ptr{T}, Ptr{Ptr{Nothing}}, Ptr{T}, Ptr{Nothing}, Ptr{T},
                Ptr{Nothing}, Ptr{T}, Ptr{Nothing}, Ptr{T}, Ptr{Nothing}, Ptr{T}, Ptr{Nothing},
                Ptr{T}, Ptr{Ptr{Nothing}}, Ptr{T}, Ptr{Nothing}, Ptr{T}, Ptr{Nothing}, Ptr{T},
                Ptr{Nothing}, Csize_t, Ptr{Nothing}, Csize_t),
               libcudnn_handle[], rnn, seqlen, yd, y, dyd, dy, dhod, dho, dcod, dco,
               wd, w, hd, h, cd, c, dxd, dx, dhd, dh, dcd, dc, ws, length(ws), rs, length(rs))
end

function backwardData(rnn::RNNDesc{T}, y, dy_, dho, dco, h, c, reserve) where T
  # Same as above, any more efficient way?
  dy = dy_ isa Integer ? zero(y) : dy_
  yd = xDesc(y)
  dx = y isa AbstractVector ? similar(dy, rnn.input) : similar(dy, rnn.input, size(dy, 2))
  dh = similar(h)
  dc = c == nothing ? nothing : similar(c)
  cudnnRNNBackwardData(rnn, 1,
    yd, y, yd, dy, hDesc(dho)..., hDesc(dco)...,
    FilterDesc(T, (1, 1, length(rnn.params))), rnn.params,
    hDesc(h)..., hDesc(c)..., xDesc(dx), dx, hDesc(dh)..., hDesc(dc)...,
    workspace[], reserve)
  return c == nothing ? (dx, dh) : (dx, dh, dc)
end

backwardData(rnn, y, dy, dho, hx, reserve) =
  backwardData(rnn, y, dy, dho, nothing, hx, nothing, reserve)

function cudnnRNNBackwardWeights(rnn::RNNDesc{T}, seqlen, xd, x, hd, h, yd, y, dwd, dw,
                                 workspace, reserve) where T
  @check ccall((:cudnnRNNBackwardWeights,libcudnn), cudnnStatus_t,
               (Ptr{Nothing}, Ptr{Nothing}, Cint,  # handle, rnnDesc, seqLength
                Ptr{Ptr{Nothing}}, Ptr{T}, #x
                Ptr{Nothing}, Ptr{T}, #hx
                Ptr{Ptr{Nothing}}, Ptr{T}, #y
                Ptr{Nothing}, Csize_t, #ws
                Ptr{Nothing}, Ptr{T}, #dw
                Ptr{Nothing}, Csize_t), #rs
               libcudnn_handle[], rnn, seqlen, xd, x, hd, h, yd, y,
               workspace, length(workspace), dwd, dw, reserve, length(reserve))
end

function backwardWeights(rnn::RNNDesc{T}, x, h, y, reserve) where T
  dw = zero(rnn.params)
  cudnnRNNBackwardWeights(rnn, 1,
    xDesc(x), x, hDesc(h)..., xDesc(y), y,
    FilterDesc(T, (1, 1, length(dw))), dw,
    workspace[], reserve)
  return params(dw, rnn.input, rnn.hidden, ngates(rnn))
end

# Interface

import ..Flux: Flux, relu
import ..Tracker: TrackedArray
using .CuArrays.CUDAnative
using .CuArrays: @cuindex, cudims

function LinearAlgebra.copy_transpose!(dst::CuArray, src::CuArray)
  function kernel(dst, src)
    I = @cuindex dst
    dst[I...] = src[reverse(I)...]
    return
  end
  blk, thr = cudims(dst)
  @cuda blocks=blk threads=thr kernel(dst, src)
  return dst
end

CuParam{T,N} = Union{CuArray{T,N},TrackedArray{T,N,CuArray{T,N}}}
CuRNN{T} = Flux.RNNCell{<:Union{typeof(tanh),typeof(relu)},<:CuParam{T,2},<:CuParam{T,1}}
CuGRU{T} = Flux.GRUCell{<:CuParam{T,2},<:CuParam{T,1}}
CuLSTM{T} = Flux.LSTMCell{<:CuParam{T,2},<:CuParam{T,1}}
CuRNNs{T} = Union{CuRNN{T},CuGRU{T},CuLSTM{T}}

function copyparams!(m::CuRNNs, d::RNNDesc)
  Wi, Wh = d.weights
  copy_transpose!(Wi, Flux.data(m.Wi))
  copy_transpose!(Wh, Flux.data(m.Wh))
  copy_transpose!(d.bias, Flux.data(m.b))
  return
end

function RNNDesc(m::CuRNNs{T}) where T
  h, i = length(m.h), size(m.Wi, 2)
  mode = m isa CuRNN ?
    (m.σ == tanh ? RNN_TANH : RNN_RELU) :
    m isa CuGRU ? GRU : LSTM
  r = RNNDesc{T}(mode, i, h)
  return r
end

const descs = WeakKeyDict()

function desc(rnn)
  d = haskey(descs, rnn) ? descs[rnn] : (descs[rnn] = RNNDesc(rnn))
  copyparams!(rnn, d)
  return d
end

import Flux.Tracker
import Flux.Tracker: data, istracked, track, unbroadcast, @grad, nobacksies

istrain(m::CuRNNs, args...) = any(x -> x isa TrackedArray, (m.Wi, m.Wh, m.b, args...))

function (m::CuRNN{T})(h::CuParam{T}, x::CuParam{T}) where T <: Union{Float32,Float64}
  result = istrain(m, h, x) ?
    track(m, x, h, m.Wi, m.Wh, m.b) :
    forward(desc(m), x, h)
  return result[2], result[1]
end

function (m::CuGRU{T})(h::CuParam{T}, x::CuParam{T}) where T <: Union{Float32,Float64}
  result = istrain(m, h, x) ?
    track(m, x, h, m.Wi, m.Wh, m.b) :
    forward(desc(m), x, h)
  return result[2], result[1]
end

function (m::CuLSTM{T})(h::NTuple{2,CuParam{T}}, x::CuParam{T}) where T <: Union{Float32,Float64}
  result = istrain(m, h, x) ?
    track(m, x, h[1], h[2], m.Wi, m.Wh, m.b) :
    forward(desc(m), x, h[1], h[2])
  return (result[2], result[3]), result[1]
end

(m::CuRNN{T})(h::CuParam{T}, x) where T <: Union{Float32,Float64} = m(h, CuArray{T}(x))
(m::CuGRU{T})(h::CuParam{T}, x) where T <: Union{Float32,Float64} = m(h, CuArray{T}(x))
(m::CuLSTM{T})(h::NTuple{2,CuParam{T}}, x) where T <: Union{Float32,Float64} = m(h, CuArray{T}(x))

@grad function (m::Union{CuRNN,CuGRU})(x, h, Wi, Wh, b)
  reserve, result = forwardTrain(desc(m), data(x), data(h))
  result, function (Δ)
    y, ho = result
    dy, dho = Δ
    h_ = hBatch(x, data(h))
    dx, dh = backwardData(descs[m], y, dy, dho, h_, reserve)
    (dWi, dWh), db = backwardWeights(descs[m], data(x), h_, y, reserve)
    nobacksies(:RNN, (dx, unbroadcast(h, dh), transpose(dWi), transpose(dWh), db))
  end
end

@grad function (m::CuLSTM)(x, h, c, Wi, Wh, b)
  reserve, result = forwardTrain(desc(m), data.((x, h, c))...)
  result, function (Δ)
    y, ho = result
    dy, dho, dco = Δ
    h_ = hBatch(x, data(h))
    c_ = hBatch(x, data(c))
    dx, dh, dc = backwardData(descs[m], y, dy, dho, dco, h_, c_, reserve)
    (dWi, dWh), db = backwardWeights(descs[m], data(x), h_, y, reserve)
    nobacksies(:RNN,
      (dx, unbroadcast(h, dh), unbroadcast(c, dc),
       transpose(dWi), transpose(dWh), db))
  end
end


text:
# CUFFT API function return values
@enum(cufftResult,
    CUFFT_SUCCESS        = 0,  #  The cuFFT operation was successful
    CUFFT_INVALID_PLAN   = 1,  #  cuFFT was passed an invalid plan handle
    CUFFT_ALLOC_FAILED   = 2,  #  cuFFT failed to allocate GPU or CPU memory
    CUFFT_INVALID_TYPE   = 3,  #  No longer used
    CUFFT_INVALID_VALUE  = 4,  #  User specified an invalid pointer or parameter
    CUFFT_INTERNAL_ERROR = 5,  #  Driver or internal cuFFT library error
    CUFFT_EXEC_FAILED    = 6,  #  Failed to execute an FFT on the GPU
    CUFFT_SETUP_FAILED   = 7,  #  The cuFFT library failed to initialize
    CUFFT_INVALID_SIZE   = 8,  #  User specified an invalid transform size
    CUFFT_UNALIGNED_DATA = 9,  #  No longer used
    CUFFT_INCOMPLETE_PARAMETER_LIST = 10, #  Missing parameters in call
    CUFFT_INVALID_DEVICE = 11, #  Execution of a plan was on different GPU than plan creation
    CUFFT_PARSE_ERROR    = 12, #  Internal plan database error
    CUFFT_NO_WORKSPACE   = 13,  #  No workspace has been provided prior to plan execution
    CUFFT_NOT_IMPLEMENTED = 14, # Function does not implement functionality for parameters given.
    CUFFT_LICENSE_ERROR  = 15, # Used in previous versions.
    CUFFT_NOT_SUPPORTED  = 16  # Operation is not supported for parameters given.
)

const cufftReal = Float32
const cufftDoubleReal = Float64

const cufftComplex = Complex64
const cufftDoubleComplex = Complex128

# CUFFT transform directions
const CUFFT_FORWARD = -1 # Forward FFT
const CUFFT_INVERSE =  1 # Inverse FFT

# CUFFT supports the following transform types
const cufftType = Cint
const CUFFT_R2C = 0x2a     # Real to Complex
const CUFFT_C2R = 0x2c     # Complex to Real
const CUFFT_C2C = 0x29     # Complex to Complex
const CUFFT_D2Z = 0x6a     # Double to Double-Complex
const CUFFT_Z2D = 0x6c     # Double-Complex to Double
const CUFFT_Z2Z = 0x69     # Double-Complex to Double-Complex

const cufftCompatibility = Cint
const   CUFFT_COMPATIBILITY_NATIVE          = 0x00
const   CUFFT_COMPATIBILITY_FFTW_PADDING    = 0x01
const   CUFFT_COMPATIBILITY_FFTW_ASYMMETRIC = 0x02
const   CUFFT_COMPATIBILITY_FFTW_ALL        = 0x03

const cufftHandle = Cint


text:
const FLAG_SIMULATION_POISSON_VARY_RANK_SMALL                 = false
const FLAG_SIMULATION_POISSON_VARY_RANK_MEDIUM                = true
const FLAG_SIMULATION_POISSON_VARY_RANK_LARGE                 = false

const FLAG_SIMULATION_POISSON_VARY_MISSING_PERCENTAGE_SMALL   = false
const FLAG_SIMULATION_POISSON_VARY_MISSING_PERCENTAGE_MEDIUM  = false
const FLAG_SIMULATION_POISSON_VARY_MISSING_PERCENTAGE_LARGE   = false

FLAG_SIMULATION_POISSON_VARY_RANK_SMALL ?
  include("simulation_poisson_vary_rank_small.jl")  : nothing

FLAG_SIMULATION_POISSON_VARY_RANK_MEDIUM ?
  include("simulation_poisson_vary_rank_medium.jl")  : nothing

# @testset "$(format("ADMM Algorithm: Small Input Simulation [Poisson 400 x 400]"))" begin
#   if SIMULATION_STATUS_POISSON_VARY_RANK_SMALL == false
#     let
#       Random.seed!(65536)
#       ROW = 400
#       COL = 400
#       for input_rank in 1:2:400
#         @printf("small case: rank = %d\n", input_rank)
#         dist = Poisson()
#         timer = TimerOutput()
#         Base.Filesystem.mkpath("./test_result/poisson/small(400x400)")
#         io = open("./test_result/poisson/small(400x400)/rank"*string(input_rank)*".log", "w")
#         truth_matrix        = rand([(FixedRankMatrix(Distributions.Poisson(5), rank = input_rank), ROW, COL)])
#         sample_model        = provide(Sampler{BernoulliModel}(), rate = 0.8)
#         input_matrix        = sample_model.draw(truth_matrix)
#         manual_type_matrix  = Array{Symbol}(undef, ROW, COL)
#         manual_type_matrix .= :Poisson
#         @timeit timer  "Poisson(400x400)" * "| rank="*string(input_rank) begin
#           completed_matrix, type_tracker, tracker = complete(A               = input_matrix,
#                                                              maxiter         = 200,
#                                                              ρ               = 0.3,
#                                                              use_autodiff    = false,
#                                                              gd_iter         = 3,
#                                                              debug_mode      = false,
#                                                              project_rank    = nothing,
#                                                              io              = io,
#                                                              type_assignment = manual_type_matrix)
#         end

#         log_simulation_result(Poisson(), completed_matrix, truth_matrix, type_tracker,tracker, io = io)
#         show(io, MIME("text/plain"), timer)
#       end
#     end
#   end
# end

# @testset "$(format("ADMM Algorithm: Medium Input Simulation [Poisson 2000 x 2000]"))" begin
#   if SIMULATION_STATUS_POISSON_VARY_RANK_MEDIUM == false
#     let
#       Random.seed!(65536)
#       ROW = 2000
#       COL = 2000
#       for input_rank in union(1, collect(10:10:500))
#         @printf("medium case: rank = %d\n", input_rank)
#         dist = Poisson()
#         timer = TimerOutput()
#         Base.Filesystem.mkpath("./test_result/poisson/medium(2000x2000)")
#         io = open("./test_result/poisson/medium(2000x2000)/rank"*string(input_rank)*".log", "w")
#         truth_matrix        = rand([(FixedRankMatrix(Distributions.Poisson(5), rank = input_rank), ROW, COL)])
#         sample_model        = provide(Sampler{BernoulliModel}(), rate = 0.8)
#         input_matrix        = sample_model.draw(truth_matrix)
#         manual_type_matrix  = Array{Symbol}(undef, ROW, COL)
#         manual_type_matrix .= :Poisson
#         @timeit timer  "Poisson(2000x2000)" * "| rank="*string(input_rank) begin
#           completed_matrix, type_tracker, tracker = complete(A               = input_matrix,
#                                                              maxiter         = 200,
#                                                              ρ               = 0.3,
#                                                              use_autodiff    = false,
#                                                              gd_iter         = 3,
#                                                              debug_mode      = false,
#                                                              project_rank    = input_rank,
#                                                              io              = io,
#                                                              type_assignment = manual_type_matrix)
#         end
#         log_simulation_result(Poisson(), completed_matrix, truth_matrix, type_tracker,tracker, io = io)
#         show(io, MIME("text/plain"), timer)
#       end
#     end
#   end
# end


text:
#=

This file specifies how to compute aerosol optical properties using the Siewert-NAI2 method

=#

"""
    $(FUNCTIONNAME)(model::MieModel{FDT}) where FDT<:NAI2

Reference: Suniti Sanghavi 2014, https://doi.org/10.1016/j.jqsrt.2013.12.015

Compute the aerosol optical properties using the Siewert-NAI2 method
Input: MieModel, holding all computation and aerosol properties
Output: AerosolOptics, holding all Greek coefficients and Cross-Sectional information
"""
function compute_aerosol_optical_properties(model::MieModel{FDT}, FT2::Type=Float64) where FDT <: NAI2

    # Unpack the model
    @unpack computation_type, aerosol, λ, polarization_type, truncation_type, r_max, nquad_radius, wigner_A, wigner_B = model

    # Extract variables from aerosol struct:
    @unpack size_distribution, nᵣ, nᵢ = aerosol

    # Imaginary part of the refractive index must be ≥ 0
    @assert nᵢ ≥ 0

    # Get the refractive index's real part type
    FT = eltype(nᵣ);
    # @assert FT == Float64 "Aerosol computations require 64bit"
    # Get radius quadrature points and weights (for mean, thus normalized):
    #

    # Just sample from 0.25%ile to 99.75%ile:
    start,stop = quantile(size_distribution,[0.0025,0.9975])
    #r, wᵣ = gauleg(nquad_radius, 0.0, r_max ; norm=true)
    r, wᵣ = gauleg(nquad_radius, start, min(stop,r_max) ; norm=true)

    # Wavenumber
    k = 2π / λ

    # Size parameter
    x_size_param = k * r # (2πr/λ)

    # Compute n_max for largest size:
    n_max = get_n_max(maximum(x_size_param))

    # Determine max amount of Gaussian quadrature points for angle dependence of
    # phase functions:
    n_mu = 2n_max - 1;

    # Obtain Gauss-Legendre quadrature points and weights for phase function
    μ, w_μ = gausslegendre(n_mu)

    # Compute π and τ functions
    leg_π, leg_τ = compute_mie_π_τ(μ, n_max)

    # Pre-allocate arrays:
    S₁    = zeros(Complex{FT}, n_mu, nquad_radius)
    S₂    = zeros(Complex{FT}, n_mu, nquad_radius)
    f₁₁   = zeros(FT, n_mu, nquad_radius)
    f₃₃   = zeros(FT, n_mu, nquad_radius)
    f₁₂   = zeros(FT, n_mu, nquad_radius)
    f₃₄   = zeros(FT, n_mu, nquad_radius)
    C_ext = zeros(FT, nquad_radius)
    C_sca = zeros(FT, nquad_radius)

    # Standardized weights for the size distribution:
    wₓ = compute_wₓ(size_distribution, wᵣ, r, r_max)

    # Loop over size parameters
    @showprogress 1 "Computing PhaseFunctions Siewert NAI-2 style ..." for i = 1:length(x_size_param)

        # Maximum expansion (see eq. A17 from de Rooij and Stap, 1984)
        n_max = get_n_max(x_size_param[i])

        # In Domke methods, we want to pre-allocate these as 2D outside of this loop.
        an = (zeros(Complex{FT}, n_max))
        bn = (zeros(Complex{FT}, n_max))

        # Weighting for sums of 2n+1
        n_ = collect(FT, 1:n_max);
        n_ = 2n_ .+ 1

        # Pre-allocate Dₙ  :
        y = x_size_param[i] * (aerosol.nᵣ - aerosol.nᵢ);
        nmx = round(Int, max(n_max, abs(y)) + 51)
        Dₙ  = zeros(Complex{FT}, nmx)

        # Compute aₙ,bₙ and S₁,S₂
        compute_mie_ab!(x_size_param[i], aerosol.nᵣ + aerosol.nᵢ * im, an, bn, Dₙ )
        compute_mie_S₁S₂!(an, bn, leg_π, leg_τ, view(S₁, :, i), view(S₂, :, i))

        # Compute Extinction and scattering cross sections:
        C_sca[i] = 2π / k^2 * (n_' * (abs2.(an) + abs2.(bn)))
        C_ext[i] = 2π / k^2 * (n_' * real(an + bn))

        # Compute scattering matrix components per size parameter:
        f₁₁[:,i] =  0.5 / x_size_param[i]^2  * real(abs2.(S₁[:,i]) + abs2.(S₂[:,i]));
        f₃₃[:,i] =  0.5 / x_size_param[i]^2  * real(S₁[:,i] .* conj(S₂[:,i]) + S₂[:,i] .* conj(S₁[:,i]));
        f₁₂[:,i] = -0.5 / x_size_param[i]^2  * real(abs2.(S₁[:,i]) - abs2.(S₂[:,i]));
        f₃₄[:,i] = -0.5 / x_size_param[i]^2  * imag(S₁[:,i] .* conj(S₂[:,i]) - S₂[:,i] .* conj(S₁[:,i]));

    end

    # Calculate bulk scattering and extinction coeffitientcs
    bulk_C_sca =  sum(wₓ .* C_sca)
    bulk_C_ext =  sum(wₓ .* C_ext)

    # Compute bulk scattering
    wr = (4π * r.^2 .*  wₓ)
    bulk_f₁₁   =  f₁₁ * wr
    bulk_f₃₃   =  f₃₃ * wr
    bulk_f₁₂   =  f₁₂ * wr
    bulk_f₃₄   =  f₃₄ * wr

    # Normalize Phase function with bulk scattering cross section.
    bulk_f₁₁ /= bulk_C_sca
    bulk_f₃₃ /= bulk_C_sca
    bulk_f₁₂ /= bulk_C_sca
    bulk_f₃₄ /= bulk_C_sca

    # Range of l-values
    l_max = length(μ);

    # Get legendre polynomials for l-max
    P, P², R², T² = compute_legendre_poly(μ, l_max)

    # Compute Greek coefficients:
    α = zeros(FT, l_max)
    β = zeros(FT, l_max)
    δ = zeros(FT, l_max)
    γ = zeros(FT, l_max)
    ϵ = zeros(FT, l_max)
    ζ = zeros(FT, l_max)

    # Compute Greek coefficients from bulk scattering matrix elements (spherical only here!)
    for l = 0:length(β) - 1

        # pre-factor:
        fac = l ≥ 2 ? (2l + 1) / 2 * sqrt(1 / ((l - 1) * (l) * (l + 1) * (l + 2))) : 0

        # Compute Greek coefficients
        # Eq 17 in Sanghavi 2014

        δ[l + 1] = (2l + 1) / 2 * w_μ' * (bulk_f₃₃ .* P[:,l + 1])
        β[l + 1] = (2l + 1) / 2 * w_μ' * (bulk_f₁₁ .* P[:,l + 1])
        γ[l + 1] = fac      * w_μ' * (bulk_f₁₂ .* P²[:,l + 1])
        ϵ[l + 1] = fac      * w_μ' * (bulk_f₃₄ .* P²[:,l + 1])
        ζ[l + 1] = fac      * w_μ' * (bulk_f₃₃ .* R²[:,l + 1] + bulk_f₁₁ .* T²[:,l + 1])
        α[l + 1] = fac      * w_μ' * (bulk_f₁₁ .* R²[:,l + 1] + bulk_f₃₃ .* T²[:,l + 1])
    end

    # Check whether this is a Dual number (if so, don't do any conversions)
    if FT2 <: AbstractFloat
        # Create GreekCoefs object with α, β, γ, δ, ϵ, ζ
        greek_coefs = GreekCoefs(convert.(FT2, α),
                                 convert.(FT2, β),
                                 convert.(FT2, γ),
                                 convert.(FT2, δ),
                                 convert.(FT2, ϵ),
                                 convert.(FT2, ζ))
        # Return the packaged AerosolOptics object
        return AerosolOptics(greek_coefs=greek_coefs, ω̃=FT2(bulk_C_sca / bulk_C_ext), k=FT2(bulk_C_ext), fᵗ=FT2(1))

    else
        greek_coefs = GreekCoefs(α, β, γ,δ,ϵ,ζ)
        return AerosolOptics(greek_coefs=greek_coefs, ω̃=(bulk_C_sca / bulk_C_ext), k=(bulk_C_ext), fᵗ=FT(1))
    end
end

function compute_ref_aerosol_extinction(model::MieModel{FDT}, FT2::Type=Float64) where FDT <: NAI2

    # Unpack the model
    @unpack computation_type, aerosol, λ, polarization_type, r_max, nquad_radius = model

    # Extract variables from aerosol struct:
    @unpack size_distribution, nᵣ, nᵢ = aerosol

    # Imaginary part of the refractive index must be ≥ 0
    @assert nᵢ ≥ 0 "Imaginary part of the refractive index must be ≥ 0 (definition)"

    # Get the refractive index's real part type
    FT = eltype(nᵣ);
    #@show FT
    #@assert FT == Float64 "Aerosol computations require 64bit"
    # Get radius quadrature points and weights (for mean, thus normalized):
    r, wᵣ = gauleg(nquad_radius, 0.0, r_max ; norm=true)

    # Wavenumber
    k = 2π / λ

    # Size parameter
    x_size_param = k * r # (2πr/λ)

    # Compute n_max for largest size:
    n_max = get_n_max(maximum(x_size_param))

    # Determine max amount of Gaussian quadrature points for angle dependence of
    # phase functions:
    n_mu = 2n_max - 1;

    # Obtain Gauss-Legendre quadrature points and weights for phase function
    μ, w_μ = gausslegendre(n_mu)

    # Compute π and τ functions
    leg_π, leg_τ = compute_mie_π_τ(μ, n_max)

    # Pre-allocate arrays:
    C_ext = zeros(FT, nquad_radius)
    C_sca = zeros(FT, nquad_radius)

    # Standardized weights for the size distribution:
    wₓ = compute_wₓ(size_distribution, wᵣ, r, r_max)

    # Loop over size parameters
    @showprogress 1 "Computing extinction XS at reference wavelength ..." for i = 1:length(x_size_param)

        # Maximum expansion (see eq. A17 from de Rooij and Stap, 1984)
        n_max = get_n_max(x_size_param[i])

        # In Domke methods, we want to pre-allocate these as 2D outside of this loop.
        an = (zeros(Complex{FT}, n_max))
        bn = (zeros(Complex{FT}, n_max))

        # Weighting for sums of 2n+1
        n_ = collect(FT, 1:n_max);
        n_ = 2n_ .+ 1

        # Pre-allocate Dn:
        y = x_size_param[i] * (aerosol.nᵣ - aerosol.nᵢ);
        nmx = round(Int, max(n_max, abs(y)) + 51)
        Dn = zeros(Complex{FT}, nmx)

        # Compute an,bn and S₁,S₂
        compute_mie_ab!(x_size_param[i], aerosol.nᵣ + aerosol.nᵢ * im, an, bn, Dn)

        # Compute Extinction and scattering cross sections:
        C_ext[i] = 2π / k^2 * (n_' * real(an + bn))
    end

    # Calculate bulk extinction coeffitient
    bulk_C_ext =  sum(wₓ .* C_ext)

    # Return the bulk extinction coeffitient
    return bulk_C_ext
end

"""
    $(FUNCTIONNAME)(aerosol::Aerosol, λ)

Compute phase function from aerosol distribution with log-normal mean μ [µm] and σ
Output: μ, w_μ, P, C_ext, C_sca, g
"""
function phase_function(aerosol::Aerosol, λ, r_max, nquad_radius)

    # Extract variables from aerosol struct:
    @unpack size_distribution, nᵣ, nᵢ = aerosol

    # Imaginary part of the refractive index must be ≥ 0
    @assert nᵢ ≥ 0 "Imaginary part of the refractive index must be ≥ 0 (definition)"

    # Get the refractive index's real part type
    FT = eltype(nᵣ);
    # @assert FT == Float64 "Aerosol computations require 64bit"
    # Get radius quadrature points and weights (for mean, thus normalized):
    r, wᵣ = gauleg(nquad_radius, 0.0, r_max ; norm=true)

    # Wavenumber
    k = 2π / λ

    # Size parameter
    x_size_param = k * r # (2πr/λ)

    # Compute n_max for largest size:
    n_max = get_n_max(maximum(x_size_param))

    # Determine max amount of Gaussian quadrature points for angle dependence of
    # phase functions:
    n_mu = 2n_max - 1;

    # Obtain Gauss-Legendre quadrature points and weights for phase function
    μ, w_μ = gausslegendre(n_mu)

    # Compute π and τ functions
    leg_π, leg_τ = compute_mie_π_τ(μ, n_max)

    # Pre-allocate arrays:
    S₁    = zeros(Complex{FT}, n_mu, nquad_radius)
    S₂    = zeros(Complex{FT}, n_mu, nquad_radius)
    f₁₁   = zeros(FT, n_mu, nquad_radius)
    C_ext = zeros(FT, nquad_radius)
    C_sca = zeros(FT, nquad_radius)

    # Standardized weights for the size distribution:
    wₓ = compute_wₓ(size_distribution, wᵣ, r, r_max)

    # Loop over size parameters
    @showprogress 1 "Computing PhaseFunctions Siewert NAI-2 style ..." for i = 1:length(x_size_param)

        # Maximum expansion (see eq. A17 from de Rooij and Stap, 1984)
        n_max = get_n_max(x_size_param[i])

        # In Domke methods, we want to pre-allocate these as 2D outside of this loop.
        an = (zeros(Complex{FT}, n_max))
        bn = (zeros(Complex{FT}, n_max))

        # Weighting for sums of 2n+1
        n_ = collect(FT, 1:n_max);
        n_ = 2n_ .+ 1

        # Pre-allocate Dn:
        y = x_size_param[i] * (aerosol.nᵣ - aerosol.nᵢ);
        nmx = round(Int, max(n_max, abs(y)) + 51)
        Dn = zeros(Complex{FT}, nmx)

        # Compute an,bn and S₁,S₂
        compute_mie_ab!(x_size_param[i], aerosol.nᵣ + aerosol.nᵢ * im, an, bn, Dn)
        compute_mie_S₁S₂!(an, bn, leg_π, leg_τ, view(S₁, :, i), view(S₂, :, i))

        # Compute Extinction and scattering cross sections:
        C_sca[i] = 2π / k^2 * (n_' * (abs2.(an) + abs2.(bn)))
        C_ext[i] = 2π / k^2 * (n_' * real(an + bn))
        # @show r[i], x_size_param[i], C_ext[i], C_ext[i]/(4π*r[i]^2), C_ext[i]*1e-8
        # Compute scattering matrix components per size parameter (might change column/row ordering):
        f₁₁[:,i] =  0.5 / x_size_param[i]^2  * real(abs2.(S₁[:,i]) + abs2.(S₂[:,i]));
    end

    # Calculate bulk scattering and extinction coeffitientcs
    bulk_C_sca =  sum(wₓ .* C_sca)
    bulk_C_ext =  sum(wₓ .* C_ext)

    # Compute bulk scattering
    wr = (4π * r.^2 .*  wₓ)
    bulk_f₁₁   =  f₁₁ * wr

    # Normalize Phase function with bulk scattering cross section.
    bulk_f₁₁ /= bulk_C_sca

    # Assymetry factor g = 0.5\int_{-1}^1 μ P(μ) dμ
    g = 1/2 * w_μ'*(μ .*  bulk_f₁₁ )
    return μ, w_μ, bulk_f₁₁, bulk_C_ext, bulk_C_sca, g
end

"""
    $(FUNCTIONNAME)(r::FT, λ::FT, nᵣ::FT, nᵢ::FT)

Compute phase function from mono-modal aerosol with radius `r` at wavelength `λ`, both in `[μm]`
Output: μ, w_μ, P, C_ext, C_sca, g
"""
function phase_function(r::FT, λ::FT, nᵣ::FT, nᵢ::FT) where {FT<:AbstractFloat}
    # Imaginary part of the refractive index must be ≥ 0 (definition)
    @assert nᵢ ≥ 0 "Imaginary part of the refractive index must be ≥ 0 (definition)"
    # Wavenumber
    k = 2π / λ

    # Size parameter
    size_param = k * r # (2πr/λ)

    # Compute n_max (doubled here to make plots smoother):
    n_max = 2get_n_max(size_param)

    # Determine max amount of Gaussian quadrature points for angle dependence of
    # phase functions:
    n_mu = 2n_max - 1;

    # Obtain Gauss-Legendre quadrature points and weights for phase function
    μ, w_μ = gausslegendre(n_mu)

    # Compute π and τ functions
    leg_π, leg_τ = compute_mie_π_τ(μ, n_max)

    # Pre-allocate arrays:
    S₁    = zeros(Complex{FT}, n_mu)
    S₂    = zeros(Complex{FT}, n_mu)

    # In Domke methods, we want to pre-allocate these as 2D outside of this loop.
    an = (zeros(Complex{FT}, n_max))
    bn = (zeros(Complex{FT}, n_max))

    # Weighting for sums of 2n+1
    n_ = collect(FT, 1:n_max);
    n_ = 2n_ .+ 1

    # Pre-allocate Dn:
    y = size_param * (nᵣ - nᵢ);
    nmx = round(Int, max(n_max, abs(y)) + 51)
    Dn = zeros(Complex{FT}, nmx)

    # Compute an,bn and S₁,S₂
    compute_mie_ab!(size_param, nᵣ + nᵢ * im, an, bn, Dn)
    compute_mie_S₁S₂!(an, bn, leg_π, leg_τ, S₁, S₂)

    # Compute Extinction and scattering cross sections:
    C_sca = 2π / k^2 * (n_' * (abs2.(an) + abs2.(bn)))
    C_ext = 2π / k^2 * (n_' * real(an + bn))

    # Compute scattering matrix components per size parameter (might change column/row ordering):
    f₁₁ =  0.5 / size_param^2  * real(abs2.(S₁) + abs2.(S₂));
    f₁₁ *= 4π * r.^2
    f₁₁ /= C_sca

    # Asymmetry factor g
    g = 1/2 * w_μ'*(μ .*  f₁₁ )
    return μ, w_μ, f₁₁, C_ext, C_sca, g
end


text:
mutable struct IdealCamera <: AbstractCamera
    landmarks_::Vector{Landmark}
    last_observation_::Vector{Vector{Float64}}
    distance_range::Vector{Float64}
    direction_range::Vector{Float64}
    function IdealCamera(
        landmark::Vector{Landmark},
        distance_range = [0.5, 6.0],
        direction_range = [-pi / 3, pi / 3],
    )
        new(landmark, Vector{Vector{Float64}}[], distance_range, direction_range)
    end
end

function visible(camera::IdealCamera, polarpos::Vector{Float64})
    return camera.distance_range[1] <= polarpos[1] <= camera.distance_range[2] &&
           camera.direction_range[1] <= polarpos[2] <= camera.direction_range[2]
end

function observation_function(camera_pose::Vector{Float64}, landmark_pos::Vector{Float64})
    diff = landmark_pos .- camera_pose[1:2]
    ϕ = atan(diff[2], diff[1]) - camera_pose[3]
    while ϕ >= pi
        ϕ -= 2 * pi
    end
    while ϕ < -pi
        ϕ += 2 * pi
    end
    return [sqrt(sum(diff .* diff)), ϕ]
end

function observations(camera::Union{IdealCamera,Nothing}, camera_pose::Vector{Float64})
    if typeof(camera) == Nothing
        return nothing
    end
    n = size(camera.landmarks_)[1]
    observed = [[1.0]]
    pop!(observed)
    for i = 1:n
        z = observation_function(camera_pose, camera.landmarks_[i].pos_)
        if visible(camera, z)
            push!(observed, [z..., camera.landmarks_[i].id])
        end
    end
    camera.last_observation_ = copy(observed)
    return observed
end

function draw(
    camera::Union{IdealCamera,Nothing},
    camera_pose::Vector{Float64},
    p::Plot{T},
) where {T}
    if camera == nothing
        return
    end
    for obsv in camera.last_observation_
        x, y, θ = camera_pose
        distance, direction = obsv[1], obsv[2]
        lx = x + distance * cos(direction + θ)
        ly = y + distance * sin(direction + θ)
        p = plot!([x, lx], [y, ly], color = "pink", legend = nothing)
    end
end

struct Uniform2D
    low::Vector{Float64}
    upp::Vector{Float64}
    uni::Uniform{Float64}
    function Uniform2D(x::Vector{Float64}, y::Vector{Float64})
        new([x[1], y[1]], [x[2], y[2]], Uniform())
    end
end

mutable struct RealCamera <: AbstractCamera
    landmarks_::Vector{Landmark}
    last_observation_::Vector{Vector{Float64}}
    distance_range::Vector{Float64}
    direction_range::Vector{Float64}
    distance_noise_rate::Float64
    direction_noise::Float64
    distance_bias_rate_std::Float64
    direction_bias::Float64
    phantom_prob::Float64
    phantom_distrib::Uniform2D
    overlook_prob::Float64
    occlusion_prob::Float64
    function RealCamera(
        landmarks::Vector{Landmark},
        distance_range = [0.5, 6.0],
        direction_range = [-pi / 3, pi / 3];
        distance_noise_rate = 0.1,
        direction_noise = pi / 90,
        distance_bias_rate_stddev = 0.1,
        direction_bias_stddev = pi / 90,
        phantom_prob = 0.0,
        phantom_range_x = [-5.0, 5.0],
        phantom_range_y = [-5.0, 5.0],
        overlook_prob = 0.1,
        occlusion_prob = 0.0,
    )

        new(
            landmarks,
            Vector{Vector{Float64}}(undef, 0),
            distance_range,
            direction_range,
            distance_noise_rate,
            direction_noise,
            rand(Normal(0.0, distance_bias_rate_stddev)),
            rand(Normal(0.0, direction_bias_stddev)),
            phantom_prob,
            Uniform2D(phantom_range_x, phantom_range_y),
            overlook_prob,
            occlusion_prob,
        )
    end
end

function visible(camera::AbstractCamera, polarpos::Vector{Float64})
    return camera.distance_range[1] <= polarpos[1] <= camera.distance_range[2] &&
           camera.direction_range[1] <= polarpos[2] <= camera.direction_range[2]
end

function apply_noise(camera::AbstractCamera, z::Vector{Float64})
    d = z[1]
    ϕ = z[2]
    errored_d = rand(Normal(d, d * camera.distance_noise_rate))
    errored_ϕ = rand(Normal(ϕ, camera.direction_noise))
    return [errored_d, errored_ϕ]
end

function apply_bias(camera::AbstractCamera, z::Vector{Float64})
    d = z[1]
    ϕ = z[2]
    d += d * camera.distance_bias_rate_std
    ϕ += camera.direction_bias
    return [d, ϕ]
end

function gen_phantom(dist::Uniform2D)
    c = dist.upp - dist.low
    x = [rand(dist.uni) for i = 1:2]
    return dist.low .+ (x .* c)
end

function if_overlook(camera::AbstractCamera)
    if rand(Uniform()) < camera.overlook_prob
        return true
    else
        return false
    end
end

function apply_occlusion(camera::AbstractCamera, z::Vector{Float64})
    if rand(Uniform()) < camera.occlusion_prob
        return [z[1] + rand(Uniform()) * (camera.distance_range[2] - z[1]), z[2]]
    else
        return z
    end
end

function observations(
    camera::RealCamera,
    camera_pose::Vector{Float64};
    noise = false,
    bias = false,
    phantom = false,
    overlook = false,
    occlusion = false,
)
    n = size(camera.landmarks_)[1]
    observed = [[1.0]]
    pop!(observed)
    for i = 1:n
        z = observation_function(camera_pose, camera.landmarks_[i].pos_)
        if phantom
            if rand(Uniform()) < camera.phantom_prob
                phantom_pos = gen_phantom(camera.phantom_distrib)
                z = observation_function(camera_pose, phantom_pos)
            end
        end
        if overlook && if_overlook(camera)
            # not added
            continue
        end
        if occlusion
            z = apply_occlusion(camera, z)
        end
        if visible(camera, z)
            if bias
                z = apply_bias(camera, z)
            end
            if noise
                z = apply_noise(camera, z)
            end
            push!(observed, [z[1], z[2], camera.landmarks_[i].id])
        end
    end
    camera.last_observation_ = copy(observed)
    return observed
end

mutable struct PsiCamera <: AbstractCamera
    landmarks_::Vector{Landmark}
    last_observation_::Vector{Vector{Float64}}
    distance_range::Vector{Float64}
    direction_range::Vector{Float64}
    distance_noise_rate::Float64
    direction_noise::Float64
    distance_bias_rate_std::Float64
    direction_bias::Float64
    phantom_prob::Float64
    phantom_distrib::Uniform2D
    overlook_prob::Float64
    occlusion_prob::Float64
    function PsiCamera(
        landmarks::Vector{Landmark},
        distance_range = [0.5, 6.0],
        direction_range = [-pi / 3, pi / 3];
        distance_noise_rate = 0.1,
        direction_noise = pi / 90,
        distance_bias_rate_stddev = 0.1,
        direction_bias_stddev = pi / 90,
        phantom_prob = 0.0,
        phantom_range_x = [-5.0, 5.0],
        phantom_range_y = [-5.0, 5.0],
        overlook_prob = 0.1,
        occlusion_prob = 0.0,
    )

        new(
            landmarks,
            Vector{Vector{Float64}}(undef, 0),
            distance_range,
            direction_range,
            distance_noise_rate,
            direction_noise,
            rand(Normal(0.0, distance_bias_rate_stddev)),
            rand(Normal(0.0, direction_bias_stddev)),
            phantom_prob,
            Uniform2D(phantom_range_x, phantom_range_y),
            overlook_prob,
            occlusion_prob,
        )
    end
end

function observations(camera::PsiCamera, camera_pose::Vector{Float64})
    orientation_noise = pi / 90.0
    n = size(camera.landmarks_)[1]
    observed = Vector{Vector{Float64}}(undef, 0)
    for i = 1:n
        landmark = camera.landmarks_[i]
        ψ = rand(
            Normal(
                atan(camera_pose[2] - landmark.pos_[2], camera_pose[1] - landmark.pos_[1]),
                orientation_noise,
            ),
        )
        z = observation_function(camera_pose, landmark.pos_)

        if rand(Uniform()) < camera.phantom_prob
            phantom_pos = gen_phantom(camera.phantom_distrib)
            z = observation_function(camera_pose, phantom_pos)
        end

        if if_overlook(camera)
            # not added
            continue
        end

        if visible(camera, z)
            z = apply_bias(camera, z)
            z = apply_noise(camera, z)
            push!(observed, [z[1], z[2], ψ, landmark.id])
        end
    end
    camera.last_observation_ = copy(observed)
    return observed
end

function draw(camera::AbstractCamera, camera_pose::Vector{Float64}, p::Plot{T}) where {T}
    for obsv in camera.last_observation_
        x, y, θ = camera_pose
        distance, direction = obsv[1], obsv[2]
        lx = x + distance * cos(direction + θ)
        ly = y + distance * sin(direction + θ)
        p = plot!([x, lx], [y, ly], color = "pink", legend = nothing)
    end
end

max_stars_repo_name: soblin/JuliaProbo

max_issues_repo_name: soblin/JuliaProbo

max_forks_repo_name: soblin/JuliaProbo


text:
x = 1
y = 1
z = 1

# MIDI Mapping
# https://github.com/Ableton/push-interface/blob/master/doc/AbletonPush2MIDIDisplayInterface.asc#23-midi-mapping
function midi_event(msg)
    global x, y, z
    if 0xb0 == first(msg) # knobs
        length(msg) < 3 && return
        knob = msg[2]
        val = 0x7f == msg[3] ? -1 : 1
        if 0x47 == knob     # knob 1
            x += val
        elseif 0x48 == knob # knob 2
            y += val
        elseif 0x49 == knob # knob 3
            z += val
        end
        Atom.msg("updateWorkspace")
    end
end

function midi_callback(msg)
    # @info :midi_callback msg
    midi_event(msg)
end

import PushInterface: RTMIDI
import .RTMIDI: rtmidi_in_create_default, rtmidi_open_port, rtmidi_in_set_callback
import .RTMIDI: destroy, EventCB, rtmidi_callback_func

device = rtmidi_in_create_default()
rtmidi_open_port(device, 1, "Ableton Push 2 User Port")

cb_ptr = @cfunction rtmidi_callback_func Cvoid (Cdouble, Ptr{Cuchar}, Ptr{EventCB})
cond = Base.AsyncCondition()
handle = Base.unsafe_convert(Ptr{Cvoid}, cond)
ecb = EventCB(handle, 0, C_NULL)
r_ecb = Ref(ecb)
rtmidi_in_set_callback(device, cb_ptr, r_ecb)
finalizer(destroy, device)

function _callback_async_loop(cond, r_ecb)
    while isopen(cond)
        wait(cond)
        msg = codeunits(unsafe_string(r_ecb[].message))
        midi_callback(msg)
    end
end

cbloop = Task() do
    _callback_async_loop(cond, r_ecb)
end
schedule(cbloop)


text:
function get_bgd(obj::GalileoOrbitalElements)
    return jcall(obj, "getBGD", jdouble, ())
end

function get_bgde1e5a(obj::GalileoOrbitalElements)
    return jcall(obj, "getBGDE1E5a", jdouble, ())
end

function get_bgde5b_e1(obj::GalileoOrbitalElements)
    return jcall(obj, "getBGDE5bE1", jdouble, ())
end

function get_iod_nav(obj::GalileoOrbitalElements)
    return jcall(obj, "getIODNav", jint, ())
end


text:
@testset "Cost.jl" begin
    filename = "data\\test.cost"
    df = DDR2import.Cost.read(filename)
    @test df[1, :FLIGHTID] == 139729486
    @test df[2, :COUNTRYCODE] == "LB"
    @test df[3, :COST] == 86.7351
    #LAST TEST
    @test df[7, :COUNTRYCODE] == "LY"
end


text:
# combinations(a, k) = Combinatorics.combinations(a, k)

Combinatorics.combinations(a, k::Integer, ::Type{Val{:unique}}) = combinations(a, k)
Combinatorics.combinations(a, k::Integer, ::Type{Val{:repeated}}) = multicombinations(a, k)