text:
"""
  IsingJ

  ver 0.2.2 Released Sep 2020
  This code is distributed under the constitution
  of MIT LICENSE.

  module Initial

  2019/04/24  Released by PearCandy                 """

module MCsampling
    export mcsampling

""" Algorithms for solving classical Ising model
    Args:
        spin:
            spin
        sweepnum:
            toal  steps for sampling
        warm_up:
            warm up steps before sampling
        N:
            simulation (Matrix) size
        j0:
            exchange coupling
        ext:
            magnetic field
        temp:
            temperature (Tc=2.26918)                """

    function mcsampling(spin,sweepnum,warmup,N,j0,ext,temp)
        mag=0
        eng=0
        mag2=0
        eng2=0

        for io in 1:sweepnum+warmup    # montecarlo sweep
            for i in 1:N
                for j in 1:N
                    E_old=get_energy(spin,i,j,j0,N,ext)
                    spin[i,j]=(-1.0)*spin[i,j]          # spin flip at j-site
                    E_new=get_energy(spin,i,j,j0,N,ext)
                    delta=E_new-E_old                   # energy diff after spin flip
                    if(delta <= 0)                      # judgement of spin flip or not
                        spin[i,j]=spin[i,j]
                    else
                        if(rand() <= exp((-1.0)*delta/temp))
                            spin[i,j]=spin[i,j]
                        else
                            spin[i,j]=(-1.0)*spin[i,j]
                        end
                    end
                end
            end

            if(io > warmup)
                # 1. total spin
                total_spin=sampling_magnetization(spin,N)
                mag=mag+total_spin/sweepnum
                mag2=mag2+(total_spin*total_spin)/sweepnum

                # 2. total energy
                total_energy=sampling_energy(spin,N,j0,ext)
                eng=eng+total_energy/sweepnum
                eng2=eng2+(total_energy*total_energy)/sweepnum
            else
            end
        end
        # cv, chi
        cv=(eng2-eng*eng)/(temp*temp)
        chi=(mag2-mag*mag)/(temp)
        return spin,mag,eng,cv,chi
    end

    # sampling energy
    function sampling_energy(spin,N,j0,ext)
        total_energy=0
        for i in 1:N
            for j in 1:N
                if(i==N)
                    right=spin[1,j]
                else
                    right=spin[i+1,j]
                end
                if(j==N)
                    down=spin[i,1]
                else
                    down=spin[i,j+1]
                end
                total_energy=total_energy+(-j0)*spin[i,j]*(right+down)-ext*spin[i,j]
            end
        end
        return total_energy/(N*N)
    end

   # sampling magnetization
    function sampling_magnetization(spin,N)
        total_spin=0
        for i in 1:N
            for j in 1:N
                total_spin=total_spin+spin[i,j]
            end
        end
        return total_spin/(N*N)
    end

  # get energy
    function get_energy(spin,i,j,j0,N,ext)
        if(i==1)
            left  = spin[N,j]
            right = spin[2,j]
        elseif(i==N)
            left  = spin[N-1,j]
            right = spin[1,j]
        else
            left  = spin[i-1,j]
            right = spin[i+1,j]
        end
        if(j==1)
            up    = spin[i,N]
            down  = spin[i,2]
        elseif(j==N)
            up    = spin[i,N-1]
            down  = spin[i,1]
        else
            up    = spin[i,j-1]
            down  = spin[i,j+1]
        end

        deltaE = (-1.0)*j0*(spin[i,j]*(left+right+up+down))-spin[i,j]*ext
        return deltaE
    end
end #module


text:
# Sample a gas fitness DGP and estimate the parameters for different lenghts of
# time series
using JLD2
using Hwloc

counts = Hwloc.histmap(Hwloc.topology_load())
ncores = counts[:Core]
if nworkers()<ncores
    addprocs(ncores - 1  - nworkers())
end
using DynNets
@everywhere using DynNets
#
N_est = 150
N = 20 ;  NTV =1
indTvPars = falses(N);indTvPars[1:NTV] = true
Tvals = [100,250,500,750]
N_T = length(Tvals)
est_times = SharedArray{Float64}(zeros(N_T,N_est))
est_w = SharedArray{Float64}(zeros(N_T,N_est,N))
est_A = SharedArray{Float64}(zeros(N_T,N_est,NTV))
est_B = SharedArray{Float64}(zeros(N_T,N_est,NTV))
est_conv_flag = SharedArray{Bool,2}(N_T,N_est)
est_flag = SharedArray{Bool,2}(N_T,N_est)

real_A = SharedArray{Float64}(zeros(N_T,N_est,NTV))
real_B = SharedArray{Float64}(zeros(N_T,N_est,NTV))
 @progress "T"   for ind = 1:N_T
    T = Tvals[ind]
    @sync @parallel    for n = 1:N_est
        sim_A = rand((0.5:0.00001:0.9),NTV)
        sim_B = rand((0.4:0.00001:0.999),NTV)
        real_A[ind,n,1:NTV] = sim_A
        real_B[ind,n,1:NTV] = sim_B
        Y_T,fit_T =  DynNets.sampleParDynNetsBin1(T,"Fit-Multi-gas",
                                                    indTvPar = [trues(NTV);falses(N-NTV)],
                                                    α = 5,β=1,
                                                    Agas =sim_A , Bgas = sim_B  )
        degsT = squeeze(sum(Y_T,2),2)
        model_data = DynNets.SdErgm1(degsT,[zeros(i)for i in [N,NTV,NTV]], indTvPars)
        #estimate
        est_times[ind,n] = @elapsed hat_pars ,hat_inds,conv_flag = DynNets.estimate(model_data)

        est_w[ind,n,:] = hat_pars[1]
        est_A[ind,n,:] = hat_pars[2]
        est_B[ind,n,:] = hat_pars[3]
        est_conv_flag[ind,n] = conv_flag
        est_flag[ind,n] = true
        #println([T,n])
    end

    # save
    save_path = "./data/estimatesTest/EstVarT.jld"
    @save(save_path)
    println([T,sum(est_times[ind,:])])
end

ext: jl


text:
using Revise;
using SortingLab;

N = 100_000_000;
K = 100;

svec = rand("id".*dec.(1:N÷K, 10), N);

@time svec_sorted = radixsort(svec, false, (22, 0x3fffff));
issorted(svec_sorted) # true

@time svec_sorted = radixsort(svec, false, (11, 0x7ff));
issorted(svec_sorted) # true

@time svec_sorted = radixsort(svec); # 16； fastest
issorted(svec_sorted) # true

function quirkysort(svec)
    fsortperm()


text:
using GalaxyWarp
using Test
using Parameters
using PCHIPInterpolation
using Optim
using BlackBoxOptim
using CairoMakie
using Distributions: Normal

@testset "1D Signal Fitting" begin
    vv=collect(-400.0:20.0:400.0)

    test_parss=[[1.2,-100.0,22.0],
        [1.2,-100.0,22.0,  2.0,  11.0,30.0],
        [0.5,-100.0,22.0,  0.6,-170.0,26.0,  2.0,11.0,30.0],
        [0.5,-100.0,22.0,  0.5,-130.0,16.0,  2.0,11.0,30.0]
    ]
    for test_pars in test_parss
        test_signal=GalaxyWarp.gmodel(vv,test_pars)
        result_pars = GalaxyWarp.fit_pixel(test_signal,0.1,vv,0.1;smin=0.8,smax=100.0,svmin=15,svmax=34.0)[1]

        AA=test_pars[1:3:end]
        VVc=test_pars[2:3:end]
        SSv=test_pars[3:3:end]

        fAA=result_pars[1:3:end]
        fVVc=result_pars[2:3:end]
        fSSv=result_pars[3:3:end]
        @testset "$(size(AA)[1]) gaussians" begin
            @test isapprox(size(AA)[1] ,size(fAA)[1], atol=1)
            if size(AA)[1] == size(fAA)[1]
                    @test isapprox(minimum(abs.(sort(AA) .- sort(fAA))), 0.0, atol=0.5)
                    @test isapprox(minimum(abs.(sort(VVc) .- sort(fVVc))), 0.0, atol=15.0)
                    @test isapprox(minimum(abs.(sort(SSv) .- sort(fSSv))), 0.0, atol=10.0)
            end
        end
    end

    for test_pars in test_parss
        for rms in [0.1,0.15,0.2]
            test_signal=GalaxyWarp.gmodel(vv,test_pars)
            test_signal=rand.(GalaxyWarp.Normal.(test_signal,rms))
            result_pars = GalaxyWarp.fit_pixel(test_signal,rms,vv,rms;smin=0.8,smax=100.0,svmin=15,svmax=34.0)[1]

            AA=test_pars[1:3:end]
            VVc=test_pars[2:3:end]
            SSv=test_pars[3:3:end]

            fAA=result_pars[1:3:end]
            fVVc=result_pars[2:3:end]
            fSSv=result_pars[3:3:end]
            @testset "$(size(AA)[1]) gaussians" begin
                @test isapprox(size(AA)[1] ,size(fAA)[1], atol=1)
                if size(AA)[1] == size(fAA)[1]
                        @test isapprox(minimum(abs.(sort(AA) .- sort(fAA))), 0.0, atol=0.5)
                        @test isapprox(minimum(abs.(sort(VVc) .- sort(fVVc))), 0.0, atol=15.0)
                        @test isapprox(minimum(abs.(sort(SSv) .- sort(fSSv))), 0.0, atol=10.0)
                end
            end
        end
    end
end

@with_kw mutable struct p_test
    Fit::Array{Symbol}=[
        :pa0,:pa15,:pa5,
        :i0,:i15,:i5,
        :v0,:a,:n]

    pa0::Float64=0.0
    pa15::Float64=20.0
    pa5::Float64=30.0

    PA::Array{Symbol}=[:pa0,:pa15,:pa5]
    R_pa::Vector{Float64}=[0.0,1.5,5.0]
    ppa = Interpolator

    i0::Float64=70.0
    i15::Float64=74.0
    i5::Float64=70.0

    I::Array{Symbol}=[:i0,:i15,:i5]
    R_i::Vector{Float64}=[0.0,1.5,5.0]
    ppi = Interpolator

    v0::Float64=300.0
    a::Float64= 2.0
    n::Float64= 1.0

    drc=0.05
    Rd=collect(0.01:drc:5.0)
    dR =0.07
    dV0=300.0
    dv=32.0
    V::Array{Symbol}=[:v0,:a,:n]

    Vcirc::Function = (r::Float64,v0::Float64,a::Float64,n::Float64) -> v0*r/(r^2.0+a^2.0)^n

end

@testset "Simulated data Fitting" begin
    filename="$(@__DIR__)/simulated1/"
    run(`mkdir -p $(filename)`)
    pixel_fit_steps=800
    kinematics_fittting_steps=5000

    @time simdata=GalaxyWarp.sim_data("simulated1";theta_dist="uniform", Vcirc = x-> 320.0*x/(x^2.0+1.5^2.0)^1.2,I=x->70.0,PA=x->6.0*x, Sb = x->100.0, sigma_inst=x->10.0, cloud_mass=1e5 ,xl=3.0,yl=5.0,vl=350.0,dx=0.05,dy=0.05,dv=20.0,beam=[0.08,0.07],Rmin=0.1,dR=0.1,Rmax = 4.5,a=0.01,rms=0.05,Rcmax=2.1,Rcmap=1.7,k=3)

    @test length(simdata["Xc"])>100000
    @test length(simdata["Yc"])==length(simdata["Xc"])
    @test length(simdata["Vc"])==length(simdata["Xc"])
    @test length(simdata["Zc"])==length(simdata["Xc"])

    clouds_sim=GalaxyWarp.fit_pixelst(simdata;x0=-5.0,x1=5.0,y0=-5.5,y1=5.5,sigma=3.0, ncrit=1,N0=1,Nmax=7,vmin=-400.0, vmax=400.0,svmin=15.0, svmax=32.0,maxsteps=pixel_fit_steps, popsize=16)

    @time GalaxyWarp.cloud_fitting_diagnostics(simdata,clouds_sim,"$(filename)"*"test")

    lenX0=length(clouds_sim.Xp)
    @test lenX0>2000
    @test length(clouds_sim.Yp)==lenX0
    @test length(clouds_sim.Vp)==lenX0
    @test length(clouds_sim.belongs)==0
    @test length(clouds_sim.belongs)==0

    clouds_sim_f1=GalaxyWarp.filter_clouds(clouds_sim,[Dict("x0"=>-1.0)])

    lenX1=length(clouds_sim_f1.Xp)
    @test lenX1<lenX0
    @test length(clouds_sim_f1.Vp)==lenX1
    @test length(clouds_sim_f1.Yp)==lenX1

    @time GalaxyWarp.merge_pixel_clouds!(clouds_sim;dr=0.1,dv=50.0,mp=7,df=3.0*simdata["rms"],roll_d=8,max_d=1.0)

    lenX1m=length(clouds_sim.Xc)
    @test lenX1m>10
    @test lenX1m<lenX1
    @test length(clouds_sim.Yc)==lenX1m
    @test length(clouds_sim.Vc)==lenX1m

    clouds_sim_f2=GalaxyWarp.filter_clouds(clouds_sim,[Dict("x0"=>-1.0)])
    lenX1_2=length(clouds_sim_f2.Xp)
    @test lenX1_2 == lenX1
    @test length(clouds_sim_f2.Vp)==lenX1_2
    @test length(clouds_sim_f2.Yp)==lenX1_2

    lenX1_2c=length(clouds_sim_f2.Xc)
    @test lenX1_2c < length(clouds_sim.Xc)
    @test length(clouds_sim_f2.Vc)==lenX1_2c
    @test length(clouds_sim_f2.Yc)==lenX1_2c

    GalaxyWarp.cloud_fitting_diagnostics(simdata,clouds_sim_f1,"$(filename)"*"test")

    Pars=p_test()
    par0=copy([getfield(Pars,s) for s in Pars.Fit])
    GalaxyWarp.update_parameters!(Pars,par0)
    disks0 = GalaxyWarp.make_disks(Pars.Rd,Pars)

    function logposterior(pars::Vector{Float64})::Float64
		GalaxyWarp.update_parameters!(Pars,pars)
		likedV=GalaxyWarp.likelihood(clouds_sim,Pars;dv=30.0)
		#prior=GalaxyWarp.logprior(Pars,Priors)
		return likedV
	end

    #res =optimize(logposterior, par0, NelderMead(),Optim.Options(iterations = 3600, store_trace = true, show_trace = true,g_tol=1.0e-3))

	#MAP=Optim.minimizer(res)

    plower=[0.0,0.0,0.0,
    0.0,50.0,50.0,
    50.0,0.5,0.3]

    pupper=[180.0,180.0,50.0,
    90.0,90.0,90.0,
    500.0,4.0,3.0]

    opt_L = bbsetup(logposterior;
		SearchRange = collect(zip(plower, pupper)),
		PopulationSize=120,MaxSteps = kinematics_fittting_steps
		)
	res_L = bboptimize(opt_L)
	MAP=best_candidate(res_L)

	@show MAP
    @show logposterior(MAP)

    GalaxyWarp.update_parameters!(Pars,MAP)
    disks = GalaxyWarp.make_disks(Pars.Rd,Pars)
	GalaxyWarp.cloud_geometry!(clouds_sim,disks,Pars)

    fsky=GalaxyWarp.plot_sky(clouds_sim,disks,msize=7)
	save("$(filename)"*"sky.png",fsky)

    ## Filtering by dV
    clouds_sim_res=GalaxyWarp.filter_clouds(clouds_sim,[Dict("dv"=>10.0)])
    GalaxyWarp.cloud_fitting_diagnostics(simdata,clouds_sim_res,"$(filename)"*"test_res")

    lenX1_3=length(clouds_sim_res.Xp)
    @test lenX1_3 < length(clouds_sim.Xp)
    @test length(clouds_sim_res.Vp)==lenX1_3
    @test length(clouds_sim_res.Yp)==lenX1_3

    ## Checking how good is the fit
    pa_real0=0.0
    pa_real15=6.0*1.5
    pa_real5=6.0*5.0

    i_real0=70.0
    i_real15=70.0
    i_real5=70.0

    v0_real=320.0
    a_real=1.5
    n_real=1.2

    dpa=7.0
    di=7.0
    @test isapprox(pa_real0,Pars.pa0,atol=dpa)
    @test isapprox(pa_real15,Pars.pa15,atol=dpa)
    @test isapprox(pa_real5,Pars.pa5,atol=dpa)
    @test isapprox(i_real0,Pars.i0,atol=di)
    @test isapprox(i_real15,Pars.i15,atol=di)
    @test isapprox(i_real5,Pars.i5,atol=di)
    @test isapprox(v0_real,Pars.v0,atol=20.0 )
    @test isapprox(a_real,Pars.a,atol=0.2 )
    @test isapprox(n_real,Pars.n,atol=0.2 )

    REAL=[pa_real0,pa_real15,pa_real5,i_real0,i_real15,i_real5,v0_real,a_real,n_real]

    @show REAL
    @show logposterior(REAL)

    GalaxyWarp.update_parameters!(Pars,REAL)
    disks_real = GalaxyWarp.make_disks(Pars.Rd,Pars)
	GalaxyWarp.cloud_geometry!(clouds_sim,disks_real,Pars)
    fsky=GalaxyWarp.plot_sky(clouds_sim,disks,msize=7)
	save("$(filename)"*"sky_real.png",fsky)

	#GalaxyWarp.cloud_geometry!(clouds_sim,disks,Pars)

	fdisks=GalaxyWarp.plot_disks(Dict("MAP"=>disks,"d0"=>disks0,"REAL"=>disks_real))
	save("$(filename)"*"disks.png",fdisks)

    for ang in [0.0,15.0,30.0]
		pvdf=GalaxyWarp.plot_pvd(simdata,ang;disks=Dict("MAP"=>disks,"d0"=>disks0,"REAL"=>disks_real),clouds=clouds_sim,merged=true)
		save("$(filename)"*"pvd_$(ang).png",pvdf)
	end
end

@testset "Noisy Simulated data Fitting" begin
    simdata=GalaxyWarp.sim_data("noisy_cube";theta_dist="uniform", Vcirc = x-> 320.0*x/(x^2.0+1.5^2.0)^1.2,I=x->70.0,PA=x->6.0*x, Sb = x->100.0, sigma_inst=x->10.0, cloud_mass=1e5 ,xl=3.0,yl=5.0,vl=350.0,dx=0.05,dy=0.05,dv=20.0,beam=[0.08,0.07],Rmin=0.1,dR=0.1,Rmax = 4.5,a=0.01,rms=0.05,Rcmax=2.1,Rcmap=1.7,k=3)

    simdata["data"] .=rand.(Normal.(simdata["data"], abs.(0.1 .* simdata["data"])))

    @show size(simdata["data"])
    @test length(simdata["Xc"])>100000
    @test length(simdata["Yc"])==length(simdata["Xc"])
    @test length(simdata["Vc"])==length(simdata["Xc"])
    @test length(simdata["Zc"])==length(simdata["Xc"])

    clouds_sim=GalaxyWarp.fit_pixelst(simdata;x0=-5.0,x1=5.0,y0=-5.5,y1=5.5,sigma=3.0, ncrit=1,N0=1,Nmax=7,vmin=-400.0, vmax=400.0,svmin=15.0, svmax=32.0,maxsteps=800, popsize=16)

    GalaxyWarp.cloud_fitting_diagnostics(simdata,clouds_sim,"$(@__DIR__)/test")

    lenX0=length(clouds_sim.Xp)
    @test lenX0>2000
    @test length(clouds_sim.Yp)==lenX0
    @test length(clouds_sim.Vp)==lenX0
    @test length(clouds_sim.belongs)==0
    @test length(clouds_sim.belongs)==0

    clouds_sim_f1=GalaxyWarp.filter_clouds(clouds_sim,[Dict("x0"=>-1.0)])

    lenX1=length(clouds_sim_f1.Xp)
    @test lenX1<lenX0
    @test length(clouds_sim_f1.Vp)==lenX1
    @test length(clouds_sim_f1.Yp)==lenX1

    @time GalaxyWarp.merge_pixel_clouds!(clouds_sim;dr=0.1,dv=50.0,mp=7,df=3.0*simdata["rms"],roll_d=8,max_d=1.0)

    GalaxyWarp.cloud_fitting_diagnostics(simdata,clouds_sim_f1,"$(@__DIR__)/test")

    lenX1m=length(clouds_sim.Xc)
    @test lenX1m>10
    @test lenX1m<lenX1
    @test length(clouds_sim.Vc)==lenX1m
    @test length(clouds_sim.Yc)==lenX1m

    Pars=p_test()
    par0=copy([getfield(Pars,s) for s in Pars.Fit])
    GalaxyWarp.update_parameters!(Pars,par0)
    disks0 = GalaxyWarp.make_disks(Pars.Rd,Pars)

    function logposterior(pars::Vector{Float64})::Float64
		GalaxyWarp.update_parameters!(Pars,pars)
		likedV=GalaxyWarp.likelihood(clouds_sim,Pars;dv=30.0)
		#prior=GalaxyWarp.logprior(Pars,Priors)
		return likedV
	end

    #res =optimize(logposterior, par0, NelderMead(),Optim.Options(iterations = 3600, store_trace = true, show_trace = true,g_tol=1.0e-3))

	#MAP=Optim.minimizer(res)

    plower=[0.0,0.0,0.0,
    0.0,50.0,50.0,
    50.0,0.5,0.3]

    pupper=[180.0,180.0,50.0,
    90.0,90.0,90.0,
    500.0,4.0,3.0]

    steps=10000
    opt_L = bbsetup(logposterior;
		SearchRange = collect(zip(plower, pupper)),
		PopulationSize=120,MaxSteps = steps
		)
	res_L = bboptimize(opt_L)
	MAP=best_candidate(res_L)

	@show MAP
    @show logposterior(MAP)

    GalaxyWarp.update_parameters!(Pars,MAP)
    disks = GalaxyWarp.make_disks(Pars.Rd,Pars)
	GalaxyWarp.cloud_geometry!(clouds_sim,disks,Pars)

    fsky=GalaxyWarp.plot_sky(clouds_sim,disks,msize=7)
	save("$(@__DIR__)/sky.png",fsky)

    pa_real0=0.0
    pa_real15=6.0*1.5
    pa_real5=6.0*5.0

    i_real0=70.0
    i_real15=70.0
    i_real5=70.0

    v0_real=320.0
    a_real=1.5
    n_real=1.2

    dpa=7.0
    di=7.0
    @test isapprox(pa_real0,Pars.pa0,atol=dpa)
    @test isapprox(pa_real15,Pars.pa15,atol=dpa)
    @test isapprox(pa_real5,Pars.pa5,atol=dpa)
    @test isapprox(i_real0,Pars.i0,atol=di)
    @test isapprox(i_real15,Pars.i15,atol=di)
    @test isapprox(i_real5,Pars.i5,atol=di)
    @test isapprox(v0_real,Pars.v0,atol=20.0 )
    @test isapprox(a_real,Pars.a,atol=0.2 )
    @test isapprox(n_real,Pars.n,atol=0.2 )

    REAL=[pa_real0,pa_real15,pa_real5,i_real0,i_real15,i_real5,v0_real,a_real,n_real]

    @show REAL
    @show logposterior(REAL)

    GalaxyWarp.update_parameters!(Pars,REAL)
    disks_real = GalaxyWarp.make_disks(Pars.Rd,Pars)
	GalaxyWarp.cloud_geometry!(clouds_sim,disks_real,Pars)
    fsky=GalaxyWarp.plot_sky(clouds_sim,disks,msize=7)
	save("$(@__DIR__)/sky_real.png",fsky)

	#GalaxyWarp.cloud_geometry!(clouds_sim,disks,Pars)

	fdisks=GalaxyWarp.plot_disks(Dict("MAP"=>disks,"d0"=>disks0,"REAL"=>disks_real))
	save("$(@__DIR__)/disks.png",fdisks)

    for ang in [0.0,15.0,30.0]
		pvdf=GalaxyWarp.plot_pvd(simdata,ang;disks=Dict("MAP"=>disks,"d0"=>disks0,"REAL"=>disks_real),clouds=clouds_sim,merged=true)
		save("$(@__DIR__)/pvd_$(ang).png",pvdf)
	end
end

const G = 4.304574992e-06

@with_kw mutable struct p_mrk79_
    Fit::Array{Symbol}=[
        :pa0,:pa15,:pa5,
        :i0,:i15,:i5,
        :v0,:a,:n]

    pa0::Float64=130.0
    pa15::Float64=130.0
    pa5::Float64=130.0

    PA::Array{Symbol}=[:pa0,:pa15,:pa5]
    R_pa::Vector{Float64}=[0.0,5.0,15.0]
    ppa = Interpolator

    i0::Float64=50.0
    i15::Float64=50.0
    i5::Float64=50.0

    I::Array{Symbol}=[:i0,:i15,:i5]
    R_i::Vector{Float64}=[0.0,5.0,15.0]
    ppi = Interpolator

    Mbh=1.0e8
    v0::Float64=300.0
    a::Float64= 2.0
    n::Float64= 1.0

    drc=0.1
    Rd=collect(0.01:drc:15.0)
    dR =0.1
    dV0=300.0
    dv=32.0
    V::Array{Symbol}=[:v0,:a,:n]

    Vcirc::Function = (r::Float64,v0::Float64,a::Float64,n::Float64) -> sqrt(
		G * Mbh / r + (v0*r/(r^2.0+a^2.0)^n)^2.0)

end

pwd()
#@testset "real data Fitting" begin
    data=GalaxyWarp.load_cube("$(@__DIR__)/data/COcube.fits", "CO21"; z=0.022296, rms=0.0012, centerx="median", centery="median",zaxis="vel", line_rest=230.538e9, x0=-15.5, x1=15.5, y0=-15.5, y1=15.5, v0=-550.0, v1=550.0)

    figmom21=GalaxyWarp.moments(data,1400,800;x0=-10.0,x1=10.0,y0=-10.0,y1=10.0,v0=-300.0,v1=300.0,dx=1.0,dv=5.0,vmax=80.0,denoise=4.5)

    #@show size(data["data"])
    @test length(data["Xc"])>100000
    @test length(data["Yc"])==length(data["Xc"])
    @test length(data["Vc"])==length(data["Xc"])
    @test length(data["Zc"])==length(data["Xc"])

    clouds_mrk79=GalaxyWarp.fit_pixelst(data;x0=-10.0,x1=10.0,y0=-10.5,y1=10.5,sigma=3.0, ncrit=1,N0=1,Nmax=7,vmin=-400.0, vmax=400.0,svmin=15.0, svmax=32.0,maxsteps=800, popsize=16)

    GalaxyWarp.cloud_fitting_diagnostics(data,clouds_mrk79,"$(@__DIR__)/mrk79")

    lenX0=length(clouds_mrk79.Xp)
    @test lenX0>2000
    @test length(clouds_mrk79.Yp)==lenX0
    @test length(clouds_mrk79.Vp)==lenX0
    @test length(clouds_mrk79.belongs)==0
    @test length(clouds_mrk79.belongs)==0

    data["beam"]
	pixels_per_beam=Int(round(data["beam"][1]*data["beam"][2]/data["dx"]^2.0))

    @time GalaxyWarp.merge_pixel_clouds!(clouds_mrk79;dr=1.5,dv=50.0,mp=10,df=3.0*data["rms"],roll_d=8,max_d=1.0)

    GalaxyWarp.cloud_fitting_diagnostics(data,clouds_mrk79,"$(@__DIR__)/mrk79_m";slit=1.0)

    lenX1m=length(clouds_mrk79.Xc)
    @test length(clouds_mrk79.Vc)==lenX1m
    @test length(clouds_mrk79.Yc)==lenX1m

    Pars=p_mrk79_()
    par0=copy([getfield(Pars,s) for s in Pars.Fit])
    GalaxyWarp.update_parameters!(Pars,par0)
    disks0 = GalaxyWarp.make_disks(Pars.Rd,Pars)

    function logposterior(pars::Vector{Float64})::Float64
		GalaxyWarp.update_parameters!(Pars,pars)
		likedV=GalaxyWarp.likelihood(clouds_mrk79,Pars;dv=30.0)
		#prior=GalaxyWarp.logprior(Pars,Priors)
		return likedV
	end

    #res =optimize(logposterior, par0, NelderMead(),Optim.Options(iterations = 3600, store_trace = true, show_trace = true,g_tol=1.0e-3))

	#MAP=Optim.minimizer(res)

    plower=[0.0,0.0,0.0,
    0.0,0.0,0.0,
    50.0,0.5,0.3]

    pupper=[180.0,180.0,50.0,
    90.0,90.0,90.0,
    500.0,4.0,3.0]

    steps=10000
    opt_L = bbsetup(logposterior;
		SearchRange = collect(zip(plower, pupper)),
		PopulationSize=120,MaxSteps = steps
		)
	res_L = bboptimize(opt_L)
	MAP=best_candidate(res_L)

	@show MAP
    @show logposterior(MAP)

    GalaxyWarp.update_parameters!(Pars,MAP)
    disks = GalaxyWarp.make_disks(Pars.Rd,Pars)
	GalaxyWarp.cloud_geometry!(clouds_mrk79,disks,Pars)

    fsky=GalaxyWarp.plot_sky(clouds_mrk79,disks;msize=10,zoom_image=[-10.0,10.0,-10.0,10.0,2.0,2.0])
	save("$(@__DIR__)/sky.png",fsky)

    # pa_real0=0.0
    # pa_real15=6.0*1.5
    # pa_real5=6.0*5.0

    # i_real0=70.0
    # i_real15=70.0
    # i_real5=70.0

    # v0_real=320.0
    # a_real=1.5
    # n_real=1.2

    # dpa=7.0
    # di=7.0
    # @test isapprox(pa_real0,Pars.pa0,atol=dpa)
    # @test isapprox(pa_real15,Pars.pa15,atol=dpa)
    # @test isapprox(pa_real5,Pars.pa5,atol=dpa)
    # @test isapprox(i_real0,Pars.i0,atol=di)
    # @test isapprox(i_real15,Pars.i15,atol=di)
    # @test isapprox(i_real5,Pars.i5,atol=di)
    # @test isapprox(v0_real,Pars.v0,atol=20.0 )
    # @test isapprox(a_real,Pars.a,atol=0.2 )
    # @test isapprox(n_real,Pars.n,atol=0.2 )

    # REAL=[pa_real0,pa_real15,pa_real5,i_real0,i_real15,i_real5,v0_real,a_real,n_real]

    # @show REAL
    # @show logposterior(REAL)

    # GalaxyWarp.update_parameters!(Pars,REAL)
    # disks_real = GalaxyWarp.make_disks(Pars.Rd,Pars)
	# GalaxyWarp.cloud_geometry!(clouds_sim,disks_real,Pars)
    # fsky=GalaxyWarp.plot_sky(clouds_sim,disks,msize=7)
	# save("$(@__DIR__)/sky_real.png",fsky)

	#GalaxyWarp.cloud_geometry!(clouds_sim,disks,Pars)

	fdisks=GalaxyWarp.plot_disks(Dict("MAP"=>disks,"d0"=>disks0))
	save("$(@__DIR__)/mrk_79_disks.png",fdisks)

    for ang in [130.0]
		pvdf=GalaxyWarp.plot_pvd(data,ang;disks=Dict("MAP"=>disks,"d0"=>disks0),clouds=clouds_mrk79,merged=true,slit=1.0)
		save("$(@__DIR__)/mrk79_pvd_$(ang).png",pvdf)
	end
#end

const G = 4.304574992e-06
const K1 =0.0013329419990568908
function VNFW(R::Float64, M_vir::Float64, c::Float64)
	K2=119.73662477388707 * R * c
	return K1 *sqrt(M_vir) *
		sqrt((log(1.0 + K2 / M_vir^0.333333333) - K2 / (M_vir^0.3333333 *(1.0 + K2 / M_vir^0.333333333))
		) / (-c / (c + 1.0) + log(c + 1.0))) / sqrt(R)
end

@with_kw mutable struct p_5055_
    Rmax=7.5
    Fit::Array{Symbol}=[
        :pa0,:pa15,:pa5,
        :i0,:i15,:i5,
        :v0,:a,:n]

    pa0::Float64=130.0
    pa15::Float64=130.0
    pa5::Float64=130.0

    PA::Array{Symbol}=[:pa0,:pa15,:pa5]
    R_pa::Vector{Float64}=[0.0,Rmax/2,Rmax]
    ppa = Interpolator

    i0::Float64=50.0
    i15::Float64=50.0
    i5::Float64=50.0

    I::Array{Symbol}=[:i0,:i15,:i5]
    R_i::Vector{Float64}=[0.0,Rmax/2,Rmax]
    ppi = Interpolator

    Mbh=1.0e8
    v0::Float64=300.0
    a::Float64= 2.0
    n::Float64= 1.0

    drc=0.2
    Rd=collect(0.01:drc:Rmax)
    dR =0.1
    dV0=300.0
    dv=32.0
    V::Array{Symbol}=[:v0,:a,:n]

    Vcirc::Function = (r::Float64,v0::Float64,a::Float64,n::Float64) -> sqrt(
		G * Mbh / r + (v0*r/(r^2.0+a^2.0)^n)^2.0)

end

using GalaxyWarp
pwd()
#@testset "real data Fitting" begin
    run(`mkdir -p test/NGC5055/`)
    data_5055_HI=GalaxyWarp.load_cube("$(@__DIR__)/data/NGC5055_HI_lab.fits", "HI"; z=0.001678, rms=0.00035, centerx="median", centery="median",beam=[0.1,0.1,0.0],zaxis="vel", x0=-15.5, x1=15.5, y0=-15.5, y1=15.5, v0=-550.0, v1=550.0)

    figmom=GalaxyWarp.moments(data_5055_HI,1400,800;x0=-10.0,x1=10.0,y0=-10.0,y1=10.0,v0=-300.0,v1=300.0,dx=1.0,dv=5.0,vmax=80.0,denoise=4.5)

    clouds_5055_HI=GalaxyWarp.fit_pixelst(data_5055_HI;x0=-10.0,x1=10.0,y0=-10.5,y1=10.5,sigma=3.0, ncrit=1,N0=1,Nmax=7,vmin=-400.0, vmax=400.0,svmin=15.0, svmax=32.0,maxsteps=800, popsize=16)

    ## CO21

    data_5055_CO=GalaxyWarp.load_cube("$(@__DIR__)/data/NGC5055_CO21_lab.fits", "CO21"; z=0.001678, rms=0.1, centerx="median", centery="median",zaxis="vel", x0=-7.5, x1=7.5, y0=-5.5, y1=5.5, v0=-550.0, v1=550.0) #beam=[0.1,0.1,0.0],)

    figmom=GalaxyWarp.moments(data_5055_CO,1400,800;x0=-7.0,x1=7.0,y0=-5.0,y1=5.0,v0=-300.0,v1=300.0,dx=1.0,dv=5.0,vmax=80.0,denoise=4.0)
    save("$(@__DIR__)/NGC5055/CO_moments.png",figmom)

    clouds_5055_CO=GalaxyWarp.fit_pixelst(data_5055_CO;x0=-7.0,x1=7.0,y0=-5.0,y1=5.0,sigma=3.0, ncrit=1,N0=1,Nmax=5,vmin=-340.0, vmax=340.0,svmin=3.0, svmax=55.0,maxsteps=1000, popsize=16)

    GalaxyWarp.cloud_fitting_diagnostics(data_5055_CO,clouds_5055_CO,"$(@__DIR__)/NGC5055/NGC5055")

    clouds_5055_CO_f=GalaxyWarp.filter_clouds(clouds_5055_CO,[Dict("f0"=>3.5)])

    GalaxyWarp.cloud_fitting_diagnostics(data_5055_CO,clouds_5055_CO_f,"$(@__DIR__)/NGC5055/NGC5055_filt")

    data_5055_CO["beam"]
	pixels_per_beam=Int(round(data_5055_CO["beam"][1]*data_5055_CO["beam"][2]/data_5055_CO["dx"]^2.0))

    @time GalaxyWarp.merge_pixel_clouds!(clouds_5055_CO;dr=0.5,dv=50.0,mp=10,df=3.0*data_5055_CO["rms"],roll_d=8,max_d=1.0)

    GalaxyWarp.cloud_fitting_diagnostics(data_5055_CO,clouds_5055_CO,"$(@__DIR__)/NGC5055/NGC5055_m";slit=0.1)

    lenX1m=length(clouds_5055_CO.Xc)
    @test length(clouds_5055_CO.Vc)==lenX1m
    @test length(clouds_5055_CO.Yc)==lenX1m

    Pars=p_5055_()
    par0=copy([getfield(Pars,s) for s in Pars.Fit])
    GalaxyWarp.update_parameters!(Pars,par0)
    disks0 = GalaxyWarp.make_disks(Pars.Rd,Pars)

    function logposterior(pars::Vector{Float64})::Float64
		GalaxyWarp.update_parameters!(Pars,pars)
		likedV=GalaxyWarp.likelihood(clouds_5055_CO,Pars;dv=30.0)
		#prior=GalaxyWarp.logprior(Pars,Priors)
		return likedV
	end

    #res =optimize(logposterior, par0, NelderMead(),Optim.Options(iterations = 3600, store_trace = true, show_trace = true,g_tol=1.0e-3))

	#MAP=Optim.minimizer(res)

    plower=[0.0,0.0,0.0,
    0.0,0.0,0.0,
    50.0,0.1,0.3]

    pupper=[360.0,360.0,360.0,
    90.0,90.0,90.0,
    500.0,4.0,3.0]

    steps=10000
    opt_L = bbsetup(logposterior;
		SearchRange = collect(zip(plower, pupper)),
		PopulationSize=120,MaxSteps = steps
		)
	res_L = bboptimize(opt_L)
	MAP=best_candidate(res_L)

	@show MAP
    @show logposterior(MAP)

    GalaxyWarp.update_parameters!(Pars,MAP)
    disks = GalaxyWarp.make_disks(Pars.Rd,Pars)
	GalaxyWarp.cloud_geometry!(clouds_5055_CO,disks,Pars)

    fsky=GalaxyWarp.plot_sky(clouds_5055_CO,disks;msize=10,zoom_image=[-10.0,10.0,-10.0,10.0,2.0,2.0])
	save("$(@__DIR__)/sky.png",fsky)

    # pa_real0=0.0
    # pa_real15=6.0*1.5
    # pa_real5=6.0*5.0

    # i_real0=70.0
    # i_real15=70.0
    # i_real5=70.0

    # v0_real=320.0
    # a_real=1.5
    # n_real=1.2

    # dpa=7.0
    # di=7.0
    # @test isapprox(pa_real0,Pars.pa0,atol=dpa)
    # @test isapprox(pa_real15,Pars.pa15,atol=dpa)
    # @test isapprox(pa_real5,Pars.pa5,atol=dpa)
    # @test isapprox(i_real0,Pars.i0,atol=di)
    # @test isapprox(i_real15,Pars.i15,atol=di)
    # @test isapprox(i_real5,Pars.i5,atol=di)
    # @test isapprox(v0_real,Pars.v0,atol=20.0 )
    # @test isapprox(a_real,Pars.a,atol=0.2 )
    # @test isapprox(n_real,Pars.n,atol=0.2 )

    # REAL=[pa_real0,pa_real15,pa_real5,i_real0,i_real15,i_real5,v0_real,a_real,n_real]

    # @show REAL
    # @show logposterior(REAL)

    # GalaxyWarp.update_parameters!(Pars,REAL)
    # disks_real = GalaxyWarp.make_disks(Pars.Rd,Pars)
	# GalaxyWarp.cloud_geometry!(clouds_sim,disks_real,Pars)
    # fsky=GalaxyWarp.plot_sky(clouds_sim,disks,msize=7)
	# save("$(@__DIR__)/sky_real.png",fsky)

	#GalaxyWarp.cloud_geometry!(clouds_sim,disks,Pars)

	fdisks=GalaxyWarp.plot_disks(Dict("MAP"=>disks,"d0"=>disks0))
	save("$(@__DIR__)/mrk_79_disks.png",fdisks)

    for ang in [0.0,70.0,130.0]
		pvdf=GalaxyWarp.plot_pvd(data_5055_CO,ang;disks=Dict("MAP"=>disks,"d0"=>disks0),clouds=clouds_5055_CO,merged=true,slit=0.2)
		save("$(@__DIR__)/NGC5055/NGC5055_pvd_$(ang).png",pvdf)
	end
#end


text:
export weightednorm
export entropyintegral

using LinearAlgebra: norm, dot

weightednorm(dg, q, p = 2; componentwise=false) =
  weightednorm(q, Val(p), dg.MJ, componentwise)

function weightednorm(q, ::Val{p}, MJ, componentwise) where {p}
  n = map(components(q)) do f
    sum(MJ .* abs.(f) .^ p) ^ (1 // p)
  end
  componentwise ? n : norm(n, p)
end

function weightednorm(q, ::Val{Inf}, _, componentwise)
  n = map(components(q)) do f
    maximum(abs.(f))
  end
  componentwise ? n : norm(n, Inf)
end

entropyintegral(dg, q) = sum(dg.MJ .* entropy.(Ref(dg.law), q, dg.auxstate))
entropyproduct(dg, p, q) = sum(dg.MJ .* dot.(entropyvariables.(Ref(dg.law), p, dg.auxstate), q))


text:
function fpBGrd(x,Param)
if lower(Param.ProfpBGrd)== "baldaufcart"
    delta=Param.Grav/(Param.Rd*Param.T0);

    p=Param.p0*exp(-delta*x[3]);
else
    p=0;
end
return p
end


text:
"""
      SeisPatchProcess(in,out;<keyword arguments>)
Read data from disk, split into multidimensional overlapping patches,
apply processes, merge patches with tapers, and write to disk.

Processing of patches can be done in parallel by running the code using
(for example): julia -p 2 script_name.jl

Important Notice: You must make declare global variables f and f_param on
every processor. You can do this in your main function by typing (for
example):
@everywhere global f_param = ["style"=>"mxmyhxhy", "Niter"=> 100, "alpha"=>1,"fmax"=>80]
@everywhere global f = [SeisPOCS]

f is an array of functions that have the following syntax:
d2, h2 = f(d1,h1,f_param), where
param is a dictionary (Dict) of parameters for the function.

note that param should contain parameters for the patching and unpatching operations.

to execute the code type:  julia -p 4 my_code.jl where 4 can be replaced
    with the number of processors you wish to use.
*Credits: A. Stanton*
"""
function SeisPatchProcess(in,out,functions,param_patch=Dict(),param_f=Dict();patch_file="patch",patch_out="patch_function",rm=1)

extent = SeisMain.ReadTextHeader(in)
println(extent.n1)

    patches,npatch = SeisPatch(in,patch_file;param_patch...)

	list = Vector{patch_list}(undef,npatch)
	patches_out = Vector{String}(undef,npatch)

	for i=1:npatch

 		p_out=join([patch_out i])
		list[i] = patch_list(patches[i],p_out,param_f,functions)
		patches_out[i] = p_out
	end

	a = pmap(MyProcess,list)

    SeisUnPatch(patches_out,out;nt=extent.n1,param_patch...)

	if rm==1
		for ipatch = 1 : npatch
		SeisRemove(patches[ipatch])
		SeisRemove(patches_out[ipatch])
		end
	end

end

function MyProcess(params)

	f = params.functions
	p = params.param_f

	d1,h1,e1 = SeisRead(params.p_file)

	for ifunc = 1 : length(f)
		func = f[ifunc]
		d2 = func(d1;p[ifunc]...)
		d1 = copy(d2)
	#h1 = copy(h2)
	end

	SeisWrite(params.p_out,d1,h1,e1)
	return(1)

end

mutable struct patch_list
       p_file
       p_out
       param_f
       functions
end


text:
# This is a global github authentication token that is set the first time
# we authenticate and then reused
const _github_auth = Ref{GitHub.Authorization}()

function github_auth(;allow_anonymous::Bool=true)
    if !isassigned(_github_auth) || !allow_anonymous && isa(_github_auth[], GitHub.AnonymousAuth)
        # If the user is feeding us a GITHUB_TOKEN token, use it!
        if length(get(ENV, "GITHUB_TOKEN", "")) == 40
            _github_auth[] = GitHub.authenticate(ENV["GITHUB_TOKEN"])
        else
            if allow_anonymous
                _github_auth[] = GitHub.AnonymousAuth()
            else
                # If we're not allowed to return an anonymous authorization,
                # then let's create a token.
                _github_auth[] = GitHub.authenticate(obtain_token())
            end
        end
    end

    if !isa(_github_auth[], GitHub.AnonymousAuth)
        # Also shove this into Registrator, so it uses the same token.
        Registrator.CommentBot.CONFIG["github"] = Dict("token" => _github_auth[].token)
    end
    return _github_auth[]
end

function obtain_token(; ins=stdin, outs=stdout, github_api=GitHub.DEFAULT_API)
    println(outs)
    printstyled(outs, "Creating a github access token\n", bold=true)
    println(outs, """
    To continue, we need to create a GitHub access token, so that we can do
    things like fork Yggdrasil and create pull requests against it.  To create
    an access token, you will be asked to enter your GitHub credentials:
    """)

    params = Dict(
        # The only thing we really need to do is write to public repos
        "scopes" => [
            "public_repo",
        ],
        "note" => "BinaryBuilder.jl generated token",
        "note_url" => "https://github.com/JuliaPackaging/BinaryBuilder.jl",
        "fingerpint" => randstring(40)
    )

    while true
        user = nonempty_line_prompt("Username", "GitHub username:", ins=ins, outs=outs)
        password = nonempty_line_prompt("Password", "GitHub password:"; ins=ins, outs=outs, echo=false)

        # Shuffle this junk off to the GH API
        headers = Dict{String, String}("User-Agent"=>"BinaryBuilder-jl")
        auth = GitHub.UsernamePassAuth(user, password)
        resp = GitHub.gh_post(
            github_api,
            "/authorizations";
            headers=headers,
            params=params,
            handle_error=false,
            auth=auth,
        )
        if resp.status == 401 && startswith(strip(HTTP.getkv(resp.headers, "X-GitHub-OTP", "")), "required")
            otp_code = nonempty_line_prompt("Authorization code", "Two-factor authentication in use, enter auth code:")
            resp = GitHub.gh_post(
                github_api,
                "/authorizations";
                headers=merge(headers, Dict("X-GitHub-OTP" => otp_code)),
                params=params,
                handle_error=false,
                auth = auth,
            )
        end
        if resp.status == 401
            printstyled(outs, "Invalid credentials!", color=:red)
            println(outs)
            continue
        end
        if resp.status != 200
            GitHub.handle_response_error(resp)
        end
        token = JSON.parse(HTTP.payload(resp, String))["token"]

        print(outs, strip("""
        Successfully obtained GitHub authorization token!
        This token will be used for the rest of this BB session, however if you
        wish to use this token permanently, you may set it in your environment
        as via the following line in your """))
        printstyled(outs, "~/.julia/config/startup.jl"; bold=true)
        println(outs, "file:")
        println(outs)

        printstyled(outs, "    ENV[\"GITHUB_TOKEN\"] = \"$(token)\""; bold=true)
        println(outs)

        println(outs, "This token is sensitive, so only do this in a computing environment you trust.")
        println(outs)

        return token
    end
end

function create_or_update_pull_request(repo, params; auth=github_auth())
    try
        return create_pull_request(repo; params=params, auth=auth)
    catch ex
        # If it was already created, search for it so we can update it:
        if Registrator.CommentBot.is_pr_exists_exception(ex)
            prs, _ = pull_requests(repo; auth=auth, params=Dict(
                "state" => "open",
                "base" => params["base"],
                "head" => string(split(repo, "/")[1], ":", params["head"]),
            ))
            return update_pull_request(repo, first(prs).number; auth=auth, params=params)
        else
            rethrow(ex)
        end
    end
end


max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using JLBoost
using Test
using DataFrames
using Random

# @testset "smoke test" begin
#     df = DataFrame(x=[1,1,1,0,0], y = [1,1,1,0,0])
#     jlboost(df, :y; nrounds=4)
# end

@testset "smoke test" begin
    X = rand(100, 50)

    #y = 1.0*(rand(100) .> 0.5)
    y = (sum(X', dims=[1]).>25)[1,:];

    df = DataFrame(X, :auto);
    df[!,"y"] = y;

    xgb = jlboost(df, "y"; verbose=false, max_depth = 5);

    y_hat = predict(xgb, DataFrame(X, :auto))
    println(y_hat)

    y_hat = (y_hat .>=0)
    accuracy = sum(y .== y_hat)/length(y)

    ac = AUC(y_hat, df[!,"y"])
    println("Accuracy $(accuracy)")
    println("AUC: $(ac)")
end


text:
const TypedLength{N,T} = Union{NTuple{N,T}, StaticArray{Tuple{N},T}}

to_complex(::Type{A}) where {A <: StaticArray} = SVector
to_complex(::Type{A}) where {A <: NTuple} = NTuple

@generated function fft(x::A) where {N, T, A <: TypedLength{N,T}}
    fft_expr = fft_expression(N,T, to_complex(A))
    quote
    #    $(Expr(:meta, :inline))
        @inbounds begin
            Base.Cartesian.@nextract $N x x
        end
        $fft_expr
    end
end

max_issues_repo_name: chriselrod/StaticFFTs

max_forks_repo_name: chriselrod/StaticFFTs


text:
function dns_service_browse(service_type)

    result = Channel{DNSService}(0)

    p = open(`avahi-browse --no-fail --parsable --resolve $service_type`)

    atexit(()->kill(p, Base.SIGKILL))

    @async try
        while isopen(result) &&
              (process_running(p) || bytesavailable(p) > 0)

            l = readline(p)
            if !startswith(l, "=")
                continue
            end
            l = split(l, ";")
            if l[3] == "IPv6" && l[5] == service_type
                tag, net, ip, name, type, domain, target, address, port = l
                push!(result, name => (target, parse(Int16, port)))
            end
        end
    catch err
        if !isopen(result) &&
           err isa InvalideStateException &&
           err.state == closed
            # Ignore push to closed Channel error.
        else
            exception=(err, catch_backtrace())
            @error "Error reading `avahi-browse` output." exception
        end
    finally
        close(result)
        kill(p)
    end

    return result
end

function register_dns_service(name, service_type, port)
    global dns_process
    dns_process[(name, service_type)] =
    open(`avahi-publish --no-fail --service $name _tinyrpc._tcp $port`)
end


text:
module Sphere

function get_fn(argd::Dict)
    center = [float(x) for x in argd["center"]]
    radius = float(argd["radius"])
    k = float(argd["k"])

    function fn_compute(coordinates)
        d = norm(coordinates - center)
        if d > radius
            return k * d
        else
            return 0.0
        end
    end

    return fn_compute
end

export get_fn
end

max_stars_repo_name: yuhangwang/SimShape

max_issues_repo_name: yuhangwang/SimShape

max_forks_repo_name: yuhangwang/SimShape


text:
# The indentity function
function identity_sampler(rawdata)
    return rawdata
end


text:
using BlackBoxOptim

# From Hans W. Borchers in
# https://groups.google.com/forum/#!topic/julia-opt/HltM-OGjueo
function trefethen(p)
  x = p[1]; y = p[2]
  return exp(sin(50 * x)) + sin(60 * exp(y)) + sin(70 * sin(x)) +
             sin(sin(80 * y)) - sin(10 * (x + y)) + (x^2 + y^2)/4
end

p = BlackBoxOptim.fixeddim_problem(trefethen; range = (-1.0, 1.0), dims = 2, name = "trefethen")

# To just run one optimization:
bboptimize(p; max_time = 2.0, method = :generating_set_search,
  parameters = {:SaveFitnessTraceToCsv => true})

# To compare several methods in several repetitions:
#BlackBoxOptim.repeated_bboptimize(5, p, 2, [
#  :generating_set_search,
#  :adaptive_de_rand_1_bin_radiuslimited,
#  :random_search,
#  :separable_nes],
#  5.0)


text:
include("defects.jl")

function Matrix_D(position::Location, z, potential::Vector{LocalPotential})
    D = inv(
        (diag_V(potential) |> inv) -
        propagator_matrix(z, map(x -> x.loc, potential)),
    )
end

function δρ_calc(position::Location, z, potential::Vector{LocalPotential})
    D = Matrix_D(position, z, potential)
    PropVectorR = map(x -> Ξ(x.loc - position, z), potential)
    output = transpose(PropVectorR) * D * PropVectorR
    return output
end

function δρR(R::Location, s::AtomsSystem)
    μ = s.μ
    potential = s.potential
    f_int(x) = (1 / 2 * π) * δρ_calc(R, μ + 1im * x, potential)
    res = quadgk(f_int, 0, Inf, maxevals = NumEvals, rtol = ν)
    return 2 * real(res[1][1])
end

function spectral_bulk(ω, R::Location, s::AtomsSystem)
    μ = s.μ
    potential = s.potential
    pristine_spectral = -Ξ(Location(0, 0), ω + 1im * η) / π |> imag
    correction_spectral = -δρ_calc(R, ω + 1im * η, potential) / π |> imag
    return (pristine_spectral + correction_spectral)
end


text:
using JET.JETInterface
using Base.Core
using StructTypes
import JET:
    JET,
    @invoke,
    isexpr

using HTTP: Stream
using FilePathsBase: AbstractPath
import Parsers

const CC = Core.Compiler

# avoid kwargs in write due as it makes the analysis more complicated
# https://github.com/JuliaLang/julia/issues/9551
# https://discourse.julialang.org/t/untyped-keyword-arguments/24228
# https://discourse.julialang.org/t/closure-over-a-function-with-keyword-arguments-while-keeping-access-to-the-keyword-arguments/15574

function construct_error(T::DataType, d)
    struct_keys = collect(fieldnames(T))
    data_keys = collect(keys(d))
    ks = Symbol[]

    for k in struct_keys
        if !(k in data_keys)
            push!(ks, k)
        end
    end

    if isempty(ks)
        return nothing
    else
        return DataMissingKey(T,
            sort!(struct_keys),
            sort!(data_keys),
        )
    end
end

write(stream::Stream{<:Request}, data, status_code=ResponseCodes.Default()) = write(stream.message.response, data, status_code)

function write(res::Response, headers::Headers{T}, status_code=ResponseCodes.Default()) where {T}
    val = headers.val
    if !isnothing(val)
        for (header, value) in zip(fieldnames(val), fieldvalues(val))
            HTTP.setheader(res, headerize(header) => value)
        end
    end
end

function write(res::Response, data::Body{T}, status_code=ResponseCodes.Default()) where {T}
    if StructTypes.StructType(T) == StructTypes.NoStructType()
        error("Unsure how to write type $T to stream")
    else
        if !isnothing(data.val)
            b = IOBuffer()
            JSON3.write(b, data.val, allow_inf=true)
            res.body = take!(b)
        end

        m = mime_type(T)
        if !isnothing(m)
            write(res, Headers(content_type=m), status_code)
        end
    end
    res.status = Int(status_code)
end

function write(res::Response, path::AbstractPath, status_code=ResponseCodes.Default())
    body = Base.read(path)
    res.body = body
    m = mime_type(path)
    if !isnothing(m)
        write(res, Headers(content_type=m))
    end
end

function write(stream::Response, data::T, status_code=ResponseCodes.Default()) where {T}
    if StructTypes.StructType(T) != StructTypes.NoStructType()
        write(stream, Body(T), status_code)
    elseif T isa Exception
        write(stream, Body(string(data)), status_code)
    else
        error("Unable in infer correct write location please wrap in Body or Headers")
    end
end

function write(res::Response, status_code::Union{ResponseCodes.ResponseCode,Int})
    res.status = Int(status_code)
end

read(stream::Stream{<:Request}, b::Body{T}) where {T} = read(stream.message, b)
read(stream::Stream{A,B}, b) where {A<:Request,B} = read(stream.message, b)

function read(req::Request, ::Body{T}) where {T}
    d = JSON3.read(req.body)
    try
        return StructTypes.constructfrom(T, d)
    catch e
        maybe_e = construct_error(T, d)
        if isnothing(e)
            rethrow(e)
        else
            throw(maybe_e)
        end
    end
end

function construct_data(data, T::DataType)
    convert_numbers!(data, T)
    StructTypes.constructfrom(T, data)
end

function read(req::Request, ::Params{T}) where {T}
    try
        if hasfield(Request, :context)
            d = req.context[:params]
            return construct_data(d, T)
        else
            error("Params not supported on this version of HTTP")
        end
    catch e
        @debug "Failed to convert path params into $T"
        rethrow(e)
    end
end

function convert_numbers!(data::AbstractDict, T)
    for (k, t) in zip(fieldnames(T), fieldtypes(T))
        if t <: Union{Number, Missing, Nothing}
            data[k] = Parsers.parse(Float64, data[k])
        end
    end
    data
end

function read(req::Request, ::Query{T}) where {T}
    try
        q::Dict{Symbol,Any} = Dict(Symbol(k) => v for (k, v) in queryparams(req.url))
        return construct_data(q, T)
    catch e
        @debug "Failed to convert query into $T"
        rethrow(e)
    end
end

function read(req::Request, ::Headers{T}) where {T}
    fields = fieldnames(T)
    d = Dict{Symbol,Any}()

    for i in fields
        h = headerize(i)
        if HTTP.hasheader(req, h)
            d[i] = HTTP.header(req, h)
        else
            d[i] = missing
        end
    end

    try
        return construct_data(d, T)
    catch e
        rethrow(e)
    end
end

struct DispatchAnalyzer{T} <: AbstractAnalyzer
    state::AnalyzerState
    opts::BitVector
    frame_filter::T
    __cache_key::UInt
end

function DispatchAnalyzer(;
    ## a predicate, which takes `CC.InfernceState` and returns whether we want to analyze the call or not
    frame_filter=x::Core.MethodInstance -> true,
    jetconfigs...)
    state = AnalyzerState(; jetconfigs...)
    ## we want to run different analysis with a different filter, so include its hash into the cache key
    cache_key = state.param_key
    cache_key = hash(frame_filter, cache_key)
    return DispatchAnalyzer(state, BitVector(), frame_filter, cache_key)
end

## AbstractAnalyzer API requirements
JETInterface.AnalyzerState(analyzer::DispatchAnalyzer) = analyzer.state
JETInterface.AbstractAnalyzer(analyzer::DispatchAnalyzer, state::AnalyzerState) = DispatchAnalyzer(state, analyzer.opts, analyzer.frame_filter, analyzer.__cache_key)
JETInterface.ReportPass(analyzer::DispatchAnalyzer) = DispatchAnalysisPass()
JETInterface.get_cache_key(analyzer::DispatchAnalyzer) = analyzer.__cache_key

struct DispatchAnalysisPass <: ReportPass end
## ignore all reports defined by JET, since we'll just define our own reports
(::DispatchAnalysisPass)(T::Type{<:InferenceErrorReport}, @nospecialize(_...)) = return

function CC.finish!(analyzer::DispatchAnalyzer, frame::Core.Compiler.InferenceState)

    caller = frame.result

    ## get the source before running `finish!` to keep the reference to `OptimizationState`
    src = caller.src
    ## run `finish!(::AbstractAnalyzer, ::CC.InferenceState)` first to convert the optimized `IRCode` into optimized `CodeInfo`
    ret = @invoke CC.finish!(analyzer::AbstractAnalyzer, frame::CC.InferenceState)

    if analyzer.frame_filter(frame.linfo)
        ReportPass(analyzer)(IoReport, analyzer, caller, src)
    end

    return ret
end

@reportdef struct IoReport <: InferenceErrorReport
    slottypes
end

JETInterface.print_report(::IO, ::IoReport) =  "detected io"

function (::DispatchAnalysisPass)(::Type{IoReport}, analyzer::DispatchAnalyzer, caller::CC.InferenceResult, opt::CC.OptimizationState)
    (; src, linfo, slottypes, sptypes) = opt

    fn = get(slottypes, 1, nothing)
    if fn == Core.Const((@__MODULE__).read)

        add_new_report!(analyzer, caller, IoReport(caller, slottypes))

    elseif fn == Core.Const((@__MODULE__).write)

        status_code = get(slottypes, 4, nothing)
        if !(status_code isa Core.Const)
            return
        end

        data = get(slottypes, 3, nothing)
        if !(data isa Type)
            return
        end

        add_new_report!(analyzer, caller, IoReport(caller, slottypes))
    end
end

extract_type(::Type{T}) where {T} = T

function handler_writes(@nospecialize(handler))
    calls = JET.report_call(handler, Tuple{Stream}, analyzer=DispatchAnalyzer)
    reports = JET.get_reports(calls)
    fn = Core.Const((@__MODULE__).write)
    filter!(x -> x.slottypes[1] == fn, reports)
    l = map(reports) do r
        res_type = r.slottypes[3]
        res_code = r.slottypes[4].val
        (extract_type(res_type), res_code)
    end
    filter!(x -> x[1] != Any, l)
    unique!(l)
end

function handler_reads(@nospecialize(handler))
    calls = JET.report_call(handler, Tuple{Stream}, analyzer=DispatchAnalyzer)
    reports = JET.get_reports(calls)
    fn = Core.Const((@__MODULE__).read)
    filter!(x -> x.slottypes[1] == fn, reports)

    l = map(reports) do r
        res_type = r.slottypes[3]
        # extracts the type from Core.Const
        res_type = CC.widenconst(res_type)
        return res_type
    end
    unique!(l)
end

handler_reads(handler::AbstractHandler) = handler_reads(handler.fn)


text:
"""
# module Integrators

# Examples

```jldoctest
julia>
```
"""
module Integrators

export velocity_verlet

function velocity_verlet(x::Float64, v::Float64, f::Function, dt::Float64, atomic_mass::Float64)
    nextx = x + v * dt + f(x) / (2 * atomic_mass) * dt^2
    nextv = v + (f(nextx) + f(x)) / (2 * atomic_mass) * dt
    nextx, nextv
end

function velocity_verlet(x::Float64, v::Float64, f::Function, dt::Float64, atomic_mass::Float64, steps_amount::Int)
    for i in 1:steps_amount
        x, v = velocity_verlet(x, v, f, dt, atomic_mass)
    end
    x, v
end

end


text:
if !@isdefined TEST
    using Revise
end
using Mitosis
using Random, Test, LinearAlgebra, Statistics

# Define some vectors and matrices
ξ0 = [1., 0.]
x = [1.2, 0.1]
P0 = Matrix(1.0I, 2, 2)

Φ = [0.8 0.5; -0.1 0.8]
β = [0.1, 0.2]
Q = [0.2 0.0; 0.0 1.0]

H = [1.0 0.0]
R = Matrix(1.0I, 1, 1)

# and a nonlinear vector function
f(x) = [0.8 0.5; -0.1 0.8]*[atan(x[1]), atan(x[2])] + [0.1, 0.2]

# Define some transition kernels.

# We define the equivalent of the Soss model
:(m = @model ξ0 begin
          x0 ~ MvNormal(ξ0, P0) # priortransition
          y0 ~ MvNormal(H*x0, R) # partialobservation
          x1 ~ MvNormal(f(x0), Q) # nonlineartransition
          y1 ~ MvNormal(H*x1, R) # partial observations
          x2 ~ MvNormal(x1, P0) # full observation
          # return y0, y1, x2
      end;)

# We use AffineMap, LinearMap and ConstantMap
# For example
@test AffineMap(Φ, β)(x) == Φ*x + β

# Prior
# X′ ~ N(X, P0)
fullobservation = priortransition = Mitosis.kernel(Gaussian; μ=LinearMap(I(2)), Σ=ConstantMap(P0))

# Nonlinear transition and linear approximation
# X′ ~ N(f(X), Q)
nonlineartransition = Mitosis.kernel(Gaussian; μ=f, Σ=ConstantMap(Q))
# X′ ~ N(ΦX + β, Q),  with Φ*x + β ≈ f(x)
linearizedtransition = Mitosis.kernel(Gaussian; μ=AffineMap(Φ, β), Σ=ConstantMap(Q))

# Partial observation
# X′ ~ N(H*X, R)
partialobservation = Mitosis.kernel(Gaussian; μ=LinearMap(H), Σ=ConstantMap(R))

# And a deterministic copy kernel
cp2 = Mitosis.Copy{2}() # copy kernel
@test cp2(1.1) == (1.1, 1.1)

# Forward sample a Bayes net
x0 = rand(priortransition(ξ0))
y0 = rand(partialobservation(x0))
            # forward model with nonlinear transition
x1 = rand(nonlineartransition(x0))
y1 = rand(partialobservation(x1))
x2 = rand(fullobservation(x1))

# We actually want to write down the
# forward sample with explicit copies
# so every state is only used as input of a single kernel
x0 = rand(priortransition(ξ0))
x0a, x0b = cp2(x0)
y0 = rand(partialobservation(x0b))
x1 = rand(nonlineartransition(x0a))
x1a, x1b = cp2(x1)
y1 = rand(partialobservation(x1b))
x2 = rand(fullobservation(x1a))

# thus we have obtained some observations
(y0, y1, x2)

# backward steps trace forward step in reverse order
# starting from the observations

m1b, p1b = backward(BFFG(), partialobservation, y1; unfused=true)
m1a, p1a = backward(BFFG(), fullobservation, x2; unfused=true)
m1, p1 = backward(BFFG(), cp2, p1a, p1b) # reverse of copy, call each child transition with unfused=true
m0b, p0b = backward(BFFG(), partialobservation, y0; unfused=true)
            # backward filter just uses linearization
m0a, p0a = backward(BFFG(), linearizedtransition, p1; unfused=true)
m0, p0 = backward(BFFG(), cp2, p0a, p0b)
m, evidence = backward(BFFG(), priortransition, p0) # not a child of copy

# this creates messages m, m0, m1, ... and an evidence approximation

# forward sampler, requires the messages from the previous pass and
# a weighted starting term ξ0

x0 = rand(forward(BFFG(), priortransition, m, weighted(ξ0)))
x0a, x0b = rand(forward(BFFG(), cp2, m0, x0))
y0 = rand(forward(BFFG(), partialobservation, m0b, x0b))
            # guided sampler targets nonlineartransition, needs to know linearizedtransition
x1 = rand(forward(BFFG(), (nonlineartransition, linearizedtransition), m0a, x0a))
x1a, x1b = rand(forward(BFFG(), cp2, m1, x1))
y1 = rand(forward(BFFG(), partialobservation, m1b, x1b))
x2 = rand(forward(BFFG(), fullobservation, m1a, x1a))

# weighted joint posterior sample (x0, x1 random)
(x0, y0, x1, y1, x2)

# final importance weight as sum of weights of _leaves_ of guided pass
# the importance weight accounts for mismatch of linear backward filter
# and nonlinear forward model
w = x2.ll + y1.ll + y0.ll


text:
module BibUtils

import BibTeX

import OrderedCollections
import OrderedCollections: OrderedDict

using StringRepFilter

include("Ref.jl")
include("bibref.jl")
include("filtering.jl")
include("globals.jl")
include("utils.jl")

export bibpaths, add_bibpath
export Ref
export get_refdict, set_refdict!
export get_bibkey, set_bibkey!
export get_year, set_year!
export get_title, set_title!
export get_keywords, set_keywords!
export get_authors, set_authors!
export get_abstract, set_abstract!
export get_files
export bibrefs, filter_bibref

end # module


text:
import Base.length
import Base.reshape
import AxisArrays.AxisArray
import JuMP.value
import LinearAlgebra: transpose
import MutableArithmetics
using AxisArrays
using JuMP
using Printf

function string_dictkeys_tosymbols(d::Dict)
    d2 = Dict()
    for (k, v) in d
        d2[Symbol(k)] = v
    end
    return d2
end

function filter_dict_to_match_struct_field_names(d::Dict, s::DataType)
    f = fieldnames(s)
    d2 = Dict()
    for k in f
        if haskey(d, k)
            d2[k] = d[k]
        else
            @sprintf "\nWARNING: dict is missing struct field %s !\n" k
        end
    end
    return d2
end

Base.@kwdef struct Parameter
	 ###  Sets  ###
	 Storage::Array{String,1}      # Set B in math; new; B = "Elec","HotThermal","ColdThermal"
     TechClass::Array{String,1}    # Set C
	 DemandBin::UnitRange{Int64}   # Set E
	 FuelType::Array{String,1}     # Set F; new; F = {"NaturalGas"} for CHP and whatever fuel source is used for Boiler
     TimeStep::UnitRange{Int64}    # Set H
     TimeStepBat::UnitRange{Int64} # Set H union {0}
	 Subdivision::Array{String,1}	# Set K; new; elements = { "CapCost", "FuelBurn" }
	 Month::UnitRange{Int64} 	    # Set M
	 DemandMonthsBin::UnitRange{Int64}	# Set N
	 Ratchets::UnitRange{Int64}	   # Set R
     Seg::UnitRange{Int64}	       # Set S
	 Tech::Array{String,1}         # Set T in math
	 FuelBin::UnitRange{Int64}	   # To be removed
	 PricingTier::UnitRange{Int64}  # Set U: Pricing Tiers (proposed revision) #new
	 NMILRegime::Array{String,1}	# Set V: Net-metering Regimes

	 ###  Subsets and Indexed Sets  ####
	 ElecStorage::Array{String,1}  # B^{e} \subset B: Electrical energy storage systems
	 #HotTES::Array{String,1}  # B^{h} \subset B: Hot thermal energy storage systems (IGNORE)
	 #ColdTES::Array{String,1}  # B^{c} \subset B: Cold thermal energy storage systems (IGNORE)
	 #ThermalStorage::Array{String,1}  # B^{th} \subset B: Thermal energy storage systems (IGNORE)
	 TimeStepRatchetsMonth::Array{Array{Int64,1},1}   #  H_m: Time steps in month m
	 TimeStepRatchets::Array{Array{Int64,1},1}   #  H_r: Time steps in ratchet r
     TimeStepsWithGrid::Array{Int64,1}  # H_g: Time steps with grid connection
     TimeStepsWithoutGrid::Array{Int64,1}	 # H \setminus H_g: Time steps without grid connection
	 DemandLookbackMonths::Array{Any,1}   # M^{lb}: fixed Look back months considered for peak pricing
	 DemandLookbackRange::Int   # number of Look back months considered for peak pricing
	 SegByTechSubdivision::AxisArray # S_{kt}: System size segments from segmentation k applied to technology t
	 TechsInClass::AxisArray # T_c \subset T: Technologies that are in class c
	 TechsByFuelType::AxisArray # T_f \subset T: Technologies that burn fuel type f
	 TechsByNMILRegime::AxisArray # T_v \subset T: Technologies that may acess net-meterng regime v
	 #AbsorptionChillers::Array{String,1}  # T^{ac} \subset T: Absorption Chillers (IGNORE)
	 #CHPTechs::Array{String,1}  # T^{CHP} \subset T: CHP technologies (IGNORE)
	 #CoolingTechs::Array{String,1}  # T^{cl} \subset T: Cooling technologies (IGNORE)
	 ElectricTechs::Array{String,1}  # T^{e} \subset T: Electricity-producing technologies
	 #ElectricChillers::Array{String,1}  # T^{ec} \subset T: Electric chillers  (IGNORE)
	 FuelBurningTechs::Array{String,1}  # T^{f} \subset T: Fuel-burning technologies
	 #HeatingTechs::Array{String,1}  # T^{ht} \subset T: Heating technologies (IGNORE)
	 TechsNoTurndown::Array{String,1}  # T^{ac} \subset T: Technologies that cannot turn down, i.e., PV and wind

	 ###  Parameters and Tables supporting Indexed Sets ###
	 TechToNMILMapping::AxisArray  # Defines set T_v: Technologies that may be access net-metering regime v

	 ###  Scaling Parameters ###
	 TimeStepScaling::Float64  # \Delta: Time step scaling [h]

	 ###  Parameters for Costs and their Functional Forms ###
     AnnualMinCharge::Float64    # c^{amc}: Utility annual minimum charge
     MonthlyMinCharge::Float64    # c^{mmc}: Utility monthly minimum charge  (not in math; will use this in min charge calculation)
	 FixedMonthlyCharge::Float64  # c^{fmc}: Utility monthly fixed charge
	 FuelCost::AxisArray # c^{u}_{f}: Unit cost of fuel type f [$/MMBTU]  in math  (NEW)
	 ElecRate::Array{Float64, 2}  #   c^{g}_{uh}: Grid energy cost in energy demand tier u during time step h  (NEW)
	 OMperUnitSize::AxisArray # c^{om}_{t}: Operation and maintenance cost of technology t per unit of system size [$/kW]
     OMcostPerUnitProd::AxisArray

	 GridExportRates::Array{Float64, 2}  # c^{e}_{uh}: Export rate for energy in energy pricing tier u in time step h   (NEW)
	 CapCostSlope::AxisArray   # c^{cm}_{ts}: Slope of capital cost curve for technology t in segment s
     CapCostYInt::AxisArray  # c^{cb}_{ts}: Y-Intercept of capital cost curve for technology t in segment s
     CapCostX::AxisArray    # X-value of inflection point (will be changed)
	 #For the replacement of CapCostX, see new parameters SegmentLB and SegmentUB in section "System size and fuel limit parameters"
	 DemandRates::Array{Float64, 2}  # c^{r}_{re}: Cost per unit peak demand in tier e during ratchet r
	 DemandRatesMonth::Array{Float64, 2}   # c^{rm}_{mn}: Cost per unit peak demand in tier n during month m

	 ###  Demand Parameters ###
	 ElecLoad::Array{Float64,1}  # \delta^{d}_{h}: Electrical load in time step h   [kW]
	 # HeatingLoad::Array{Float64,1}  # \delta^{bo}_{h}: Heating load in time step h   [MMBTU/hr]
	 # CoolingLoad::Array{Float64,1}  # \delta^{c}_{h}: Cooling load in time step h   [kW]
     DemandLookbackPercent::Float64    # \delta^{lp}: Demand Lookback proportion [fraction]
     MaxDemandInTier::Array{Float64,1}  # \delta^{t}_{e}: Maximum power demand in ratchet e
     MaxDemandMonthsInTier::Array{Float64,1}   # \delta^{mt}_{n}: Maximum monthly power demand in tier n
	 MaxGridSales::Array{<:Real, 1}   # \delta^{gs}_{u}: Maximum allowable energy sales in tier u in math; equal to sum of LoadProfile["1R",ts] on set TimeStep for tier 1 (analogous "1W") and unlimited for "1X"
     MaxUsageInTier::Array{Float64,1}   # \delta^{tu}_{u}: Maximum monthly energy demand in tier u

	 ###  Incentive Parameters ###
	 NMILLimits::AxisArray   # i^{n}_{v}: Net metering and interconnect limits in net metering regime v [kW]
     MaxProdIncent::AxisArray      # \bar{i}_t: Upper incentive limit for technology t [$]
	 ProductionIncentiveRate::AxisArray  # i^{r}_{t}: Incentive rate for technology t [$/kWh] (NEW)
	 MaxSizeForProdIncent::AxisArray  # \bar{i}^{\sigma}_t: Maximum system size to obtain production incentive for technology t [kW]

	 ###  Technology-specific Time-series Factor Parameters ###
	 ProductionFactor::AxisArray    #f^{p}_{th}  Production factor of technology t and time step h  [unitless]  (NEW)
     # f^{fa}_{th}: Fuel burn ambient correction factor of technology t at time step h [unitless]
	 # f^{ha}_{th}: Hot water ambient correction factor of technology t at time step h [unitless]
	 # f^{ht}_{th}: Hot water thermal grade correction factor t correction factor of technology t at time step h [unitless]
	 # f^{ed}_{th}: Fuel burn ambient correction factor of technology t at time step h [unitless]

	 ###  Technology-specific Factor Parameters ###
	 TurbineDerate::AxisArray  # f^{d}_{t}: Derate factor for turbine technologyt [unitless]
     MinTurndown::AxisArray    # f^{td}_{t}:  Minimum turn down for technology t [unitless]
     pwf_prod_incent::AxisArray   # f^{pi}_t: Present worth factor for incentives for technology t [unitless]
	 LevelizationFactor::AxisArray    # f^{l}_{t}: Levelization factor of technology t [unitless]

	 ###  Generic Factor Parameters ###
	 pwf_om::Float64  # f^{om}: Operations and maintenance present worth factor [unitless]
     pwf_e::Float64   # f^{e}: Energy present worth factor [unitless]
	 r_tax_owner::Float64      # f^{tow}: Tax rate factor for owner [fraction]
     r_tax_offtaker::Float64   # f^{tot}: Tax rate factor for offtaker [fraction]

	 ###  System Size and Fuel Limit Parameters ###
	 TechClassMinSize::AxisArray   #  \ubar{b}^{\sigma}_{c}: Minimum system size for technology class c [kW]
	 MaxSize::AxisArray    #  \bar{b}^{\sigma}_{t}: Maximum system size for technology t [kW]
	 SegmentMinSize::AxisArray # \ubar{b}^{\sigma s}_{tks}: Minimum system size for technology t, subdivision k, segments
	 SegmentMaxSize::AxisArray  # \bar{b}^{\sigma s}_{tks}: Maximum system size for technology t, subdivision k, segments
	 FuelLimit::AxisArray # b^{fa}_{f}: Amount of available fuel for type f [MMBTU]   (NEW)

	 ###  Efficiency Parameters ###
	 ChargeEfficiency::AxisArray  # \eta^{esi}_{bt}: Efficiency of charging storage system b using technology t  [fraction] (NEW)
	 GridChargeEfficiency::Float64   # \eta^{esig}: Efficiency of charging electrical storage using grid power [fraction] (NEW)
     DischargeEfficiency::AxisArray  # \eta^{eso}_{b}: Efficiency of discharging storage system b [fraction] (NEW)
	 # \eta^{bo}: Boiler efficiency [fraction]
	 # \eta^{ecop}: Electric chiller efficiency [fraction]
	 # \eta^{acop}: Absorption chiller efficiency [fraction]

	 ###  Storage Parameters ###   # \ubar{w}^{bkW}_{b}: Minimum power capacity of storage system b (needs to be indexed on b )
     StorageMinChargePcent::Float64     #  \ubar{w}^{mcp}_{b}: Minimum state of charge of storage system b
     InitSOC::Float64    # w^{i}_{b} Initial percent state of charge for storage system b
     StorageMinSizeEnergy::AxisArray     # \bar{w}^{bkWh}_{b}: Maximum energy capacity of storage system b [kWh]
     StorageMaxSizeEnergy::AxisArray     # \ubar{w}^{bkWh}_{b}: Minimum energy capacity of storage system b [kWh]
     StorageMinSizePower::AxisArray     # \bar{w}^{bkW}_{b}: Maximum power capacity of storage system b [kW]
     StorageMaxSizePower::AxisArray     # \ubar{w}^{bkW}_{b}: Minimum power capacity of storage system b [kW]
     StorageMinSOC::AxisArray     #  \ubar{w}^{mcp}_{b}: Minimum state of charge of storage system b [fraction]
     StorageInitSOC::AxisArray  #Initial state of charge of storage system b [fraction]
     StorageCanGridCharge::Bool  # Boolean for storage system [fraction]

	 ###  Fuel Burn Parameters ###
	 FuelBurnSlope::AxisArray # m^\text{fm}_{t}: Fuel burn rate slope parameter for technology t
	 FuelBurnYInt::AxisArray # m^\text{fb}_{t}: Fuel burn rate slope parameter for technology t

	 ### Not used or used for calculation of other parameters ###
	 two_party_factor::Float64
     analysis_years::Int64     # Used to calculate present worth factors maybe?
     AnnualElecLoad::Float64   # Not used anymore (can just sum LoadProfile["1R",h] for all h in TimeStep
     CapCostSegCount::Int64    # Size of set S
     FuelBinCount::Int64       # Size of set F
     DemandBinCount ::Int64    # Size of set E
     DemandMonthsBinCount::Int64   # Size of set N
     TimeStepCount::Int64          # Size of set H
     Points::UnitRange{Int64}      # CapCostSegCount+1; this is going to be the size of set S^{c} now
     PricingTierCount::Int64    # Size of set U

     # new parameters for reformulation
     StorageCostPerKW::AxisArray
     StorageCostPerKWH::AxisArray
	 ElectricDerate::AxisArray

     # New Sets
     SalesTiers::UnitRange
     StorageSalesTiers::Array{Int, 1}
     NonStorageSalesTiers::Array{Int, 1}
	 SalesTiersByTech::AxisArray
	 TechsBySalesTier::AxisArray
	 CurtailmentTiers::Array{Int, 1}

    # Feature Additions
     TechToLocation::AxisArray
     MaxSizesLocation::Array{Float64, 1}
     Location::UnitRange
end

function Parameter(d::Dict)
    can_be_empty = (
        "MaxSize",
        "OMperUnitSize",
        "OMcostPerUnitProd",
        "MaxProdIncent",
        "MaxSizeForProdIncent",
        "TurbineDerate",
        "MinTurndown",
        "pwf_prod_incent",
        "LevelizationFactor",
        "NMILLimits",
        "TechClassMinSize",
		"TechsBySalesTier",
		"SalesTiersByTech",
		"NMILRegime",
		"TechsByNMILRegime"
     )
    if typeof(d["Tech"]) === Array{Any, 1}  # came from Python as empty array
        d["Tech"] = convert(Array{String, 1}, d["Tech"])
    end
    for x in can_be_empty
        if typeof(x) === Array{Any, 1}  # came from Python as empty array
            d[x] = convert(Array{Float64, 1}, d[x])
        end
    end
    d = convert(Dict{Union{String, Symbol}, Any}, d)

    d[:Seg] = 1:d["CapCostSegCount"]
    d[:Points] = 0:d["CapCostSegCount"]
    d[:Month] = 1:12
    d[:Ratchets] = 1:d["NumRatchets"]
    d[:FuelBin] = 1:d["FuelBinCount"]
	d[:PricingTier] = 1:d["PricingTierCount"]
    d[:DemandBin] = 1:d["DemandBinCount"]
    d[:DemandMonthsBin] = 1:d["DemandMonthsBinCount"]
    d[:TimeStep] = 1:d["TimeStepCount"]
    d[:TimeStepBat] = 0:d["TimeStepCount"]
	d[:SalesTiers] = 1:d["SalesTierCount"]
    n_location = length(d["MaxSizesLocation"])
    d[:Location] = 1:n_location

    # the following array manipulation may have to adapt once length(d["Subdivision"]) > 1
    seg_min_size_array = reshape(transpose(reshape(d["SegmentMinSize"], length(d[:Seg]), length(d["Tech"]))),
                                 length(d["Tech"]), length(d["Subdivision"]), length(d[:Seg]))
    seg_max_size_array = reshape(transpose(reshape(d["SegmentMaxSize"], length(d[:Seg]), length(d["Tech"]))),
                                length(d["Tech"]), length(d["Subdivision"]), length(d[:Seg]))

    # convert vectors to AxisArray's with axes for REopt JuMP model
    d["TurbineDerate"] = AxisArray(d["TurbineDerate"], d["Tech"])
    d["TechToLocation"] = vector_to_axisarray(d["TechToLocation"], d["Tech"], d[:Location])
    d["pwf_prod_incent"] = AxisArray(d["pwf_prod_incent"], d["Tech"])
    d["LevelizationFactor"] = AxisArray(d["LevelizationFactor"], d["Tech"])
    d["OMperUnitSize"] = AxisArray(d["OMperUnitSize"], d["Tech"])
    d["CapCostSlope"] = vector_to_axisarray(d["CapCostSlope"], d["Tech"], d[:Seg])
    d["CapCostYInt"] = vector_to_axisarray(d["CapCostYInt"], d["Tech"], d[:Seg])
    d["CapCostX"] = vector_to_axisarray(d["CapCostX"], d["Tech"], d[:Points])
    d["MaxProdIncent"] = AxisArray(d["MaxProdIncent"], d["Tech"])
    d["MaxSizeForProdIncent"] = AxisArray(d["MaxSizeForProdIncent"], d["Tech"])
    d["MaxSize"] = AxisArray(d["MaxSize"], d["Tech"])
    d["TechClassMinSize"] = AxisArray(d["TechClassMinSize"], d["TechClass"])
    d["MinTurndown"] = AxisArray(d["MinTurndown"], d["Tech"])
    d["DemandRates"] = transpose(reshape(d["DemandRates"], d["DemandBinCount"], d["NumRatchets"]))
    d["DemandRatesMonth"] = transpose(reshape(d["DemandRatesMonth"], d["DemandMonthsBinCount"], 12))
    d["NMILLimits"] = AxisArray(d["NMILLimits"], d["NMILRegime"])
    d["TechToNMILMapping"] = vector_to_axisarray(d["TechToNMILMapping"], d["Tech"], d["NMILRegime"])
    d["OMcostPerUnitProd"] = AxisArray(d["OMcostPerUnitProd"], d["Tech"])

    # Reformulation additions
    d["StorageCostPerKW"] = AxisArray([d["StorageCostPerKW"]], d["Storage"])
    d["StorageCostPerKWH"] = AxisArray([d["StorageCostPerKWH"]], d["Storage"])
    d["FuelCost"] = AxisArray(d["FuelCost"], d["FuelType"])
    d["ElecRate"] = transpose(reshape(d["ElecRate"], d["TimeStepCount"], d["PricingTierCount"]))
    d["GridExportRates"] = transpose(reshape(d["GridExportRates"], d["TimeStepCount"], d["SalesTierCount"]))
    d["FuelBurnSlope"] = AxisArray(d["FuelBurnSlope"], d["FuelBurningTechs"])
    d["FuelBurnYInt"] = AxisArray(d["FuelBurnYInt"], d["FuelBurningTechs"])
    d["ProductionFactor"] = vector_to_axisarray(d["ProductionFactor"], d["Tech"], d[:TimeStep])
    d["ProductionIncentiveRate"] = AxisArray(d["ProductionIncentiveRate"], d["Tech"])
    d["FuelLimit"] = AxisArray(d["FuelLimit"], d["FuelType"])
    d["ChargeEfficiency"] = vector_to_axisarray(d["ChargeEfficiency"], d["Tech"], d["Storage"])
    d["DischargeEfficiency"] = AxisArray(d["DischargeEfficiency"], d["Storage"])
    d["StorageMinSizeEnergy"] = AxisArray([d["StorageMinSizeEnergy"]], d["Storage"])
    d["StorageMaxSizeEnergy"] = AxisArray([d["StorageMaxSizeEnergy"]], d["Storage"])
    d["StorageMinSizePower"] = AxisArray([d["StorageMinSizePower"]], d["Storage"])
    d["StorageMaxSizePower"] = AxisArray([d["StorageMaxSizePower"]], d["Storage"])
    d["StorageMinSOC"] = AxisArray([d["StorageMinSOC"]], d["Storage"])
    d["StorageInitSOC"] = AxisArray([d["StorageInitSOC"]], d["Storage"])
    d["SegmentMinSize"] = AxisArray(seg_min_size_array, d["Tech"], d["Subdivision"], d[:Seg])
    d["SegmentMaxSize"] = AxisArray(seg_max_size_array, d["Tech"], d["Subdivision"], d[:Seg])
	d["ElectricDerate"] = vector_to_axisarray(d["ElectricDerate"], d["Tech"], d[:TimeStep])
    d["MaxGridSales"] = [d["MaxGridSales"]]

    # Indexed Sets
    if isempty(d["FuelType"])
        d["TechsByFuelType"] = []  # array of arrays is not empty, but needs to be for AxisArray conversion
    end
    d["SegByTechSubdivision"] = vector_to_axisarray(d["SegByTechSubdivision"], d["Subdivision"], d["Tech"])
    d["TechsByFuelType"] = AxisArray(d["TechsByFuelType"], d["FuelType"])
    d["TechsInClass"] = AxisArray(d["TechsInClass"], d["TechClass"])
	d["SalesTiersByTech"] = AxisArray(d["SalesTiersByTech"], d["Tech"])
	d["TechsBySalesTier"] = AxisArray(d["TechsBySalesTier"], d[:SalesTiers])
	d["TechsByNMILRegime"] = AxisArray(d["TechsByNMILRegime"], d["NMILRegime"])

    d = string_dictkeys_tosymbols(d)
    d = filter_dict_to_match_struct_field_names(d, Parameter)
    param = Parameter(;d...)
end

# Additional dispatches to make things easier
function length(::Symbol)
    return 1
end

JuMP.value(::Val{false}) = 0

JuMP.value(x::Number) = x  # for @expression(REopt, ExportBenefitYr1, 0.0) and similar

JuMP.value(::MutableArithmetics.Zero) = 0

function vector_to_axisarray(v::Array{<:Any, 1}, ax1::Array{<:Any, 1}, ax2::Union{UnitRange, Array{<:Any, 1}})
    l1 = length(ax1)
    l2 = length(ax2)
    a = transpose(reshape(v, l2, l1))
    return AxisArray(a, ax1, ax2)
end


text:
using OrdinaryDiffEq, Test
f(u,p,t) = 2*u
u_positive(u,t,integrator) = u > 0
terminate_if_u_pos = DiscreteCallback(u_positive, terminate!)

prob = ODEProblem(f, 1.0, (0.0, 1.0)) # positive initial condition > positive u > :Terminated
integrator = init(prob, Tsit5(), callback=terminate_if_u_pos)
sol1 = solve!(integrator)
@test sol1.retcode === :Terminated

reinit!(integrator, -1.0) # negative initial condition > negative u > :Success!
sol2 = solve!(integrator)
@test sol2.retcode == :Success

prob = ODEProblem(f, -1.0,(0.0, 1.0)) # negative initial condition > negative u > :Success!
integrator = init(prob, Tsit5(), callback=terminate_if_u_pos)
sol3 = solve!(integrator)
@test sol3.retcode == :Success


text:
#=
http://jsxgraph.org/docs/symbols/Chart.html
=#


text:
# Type and constructors
# Default constructor is for the line through two points
"""
	(type) Line{T<:AnyComplex} in the complex plane

Each `Line` type is parameterized according to the common type of its complex input arguments.
"""
struct Line{T<:AnyComplex} <: AbstractClosedCurve
	base::T
	direction::T
	Line{T}(a,b) where T = new(a,sign(b-a))
end
"""
	Line(a,b)

Construct the line passing through the two given points.

	Line(a,direction=z)

Construct the line through the given point and parallel to the complex sign of the given `direction` value. """
function Line(a::Number,b::Number)
	a,b = promote(complex(float(a)),complex(float(b)))
	Line{typeof(a)}(a,b)
end

# Construct by giving a keyword
Line(z::Number;direction) = Line(z,z+direction)
#Line(z::Number;angle) = Line(z,z+exp(complex(0,angle)))

# Complex type converters
for ctype in [:Spherical,:Polar,:Complex]
	@eval begin
		function $ctype(L::Line{T}) where T<:AnyComplex
			Line($ctype(L.base),$ctype(L.base+L.direction))
		end
	end
end

# Required methods
arclength(::Line) = Inf
point(L::Line,t::Real) = L.base + (2t-1)/(t-t^2)*L.direction
(C::Line)(t::Real) = point(C,t)

ispositive(C::Line) = true  # this seems arbitrary...

"""
	arg(L::Line,z)

Find a parameter argument `t` such that `L(t)==z` is true. For an infinite `z`, return zero (but note that `L(1)` is also infinity).

This gives undefined results if `z` is not actually on the line.
"""
function arg(L::Line,z::Number)
	isinf(z) && return float(0)
	del = z-L.base
	if abs(real(del)) > abs(imag(del))
		del = real(del)
		α = real(L.direction)
	else
		del = imag(del)
		α = imag(L.direction)
	end
	for t in realroots(del,2α-del,-α)
		0-eps(del) ≤ t ≤ 1+eps(del) && return t
	end
	return []
end
unittangent(L::Line,t::Real=0) = L.direction
function tangent(L::Line{T},t::Real) where T <: AnyComplex
	T( (1/t^2+1/(1-t)^2)*L.direction )
end

+(L::Line,z::Number) = Line(L.base+z,direction=L.direction)
-(L::Line) = Line(-L.base,direction=-L.direction)
*(L::Line,z::Number) = Line(L.base*z,direction=L.direction*sign(z))

"""
	inv(L)
Invert a line `L` through the origin. In general the inverse is a `Circle` through the inverse of any three points on the line.
"""
function inv(L::Line)
	w = 1 ./ point(L,[0.25,0.5,0.75])
	Circle(w...)
end

# Other methods
isfinite(::Line) = false
slope(L::Line) = imag(L.direction)/real(L.direction)
angle(L::Line) = angle(L.direction)
conj(L::Line) = Line(conj(L.base),direction=conj(L.direction))
reverse(L::Line) = Line(L.base,direction=-L.direction)

"""
	isapprox(L1::Line,L2::Line; tol=<default>)
	L1 ≈ L2

Determine if `L1` and `L2` represent the same line, irrespective of the type or values of its parameters. Identity is determined by agreement within `tol`, which is interpreted as the weaker of absolute and relative differences.
"""
function isapprox(L1::Line,L2::Line;tol=DEFAULT[:tol])
	dz = L1.base - L2.base
	w1,w2 = L1.direction,L2.direction
	return isapprox(real(w1)*imag(w2),imag(w1)*real(w2),rtol=tol,atol=tol) &&
		isapprox(real(w1)*imag(dz),imag(w1)*real(dz),rtol=tol,atol=tol)
end

"""
	isleft(z,L::Line)

Determine whether the number `z` lies "to the left" of line `L`. This means that the angle it makes with `tangent(L)` is in the interval (0,π).

Note that `isleft` and `isright` are *not* logical opposites; a point on the curve should give `false` in both cases.
"""
isleft(z::Number,L::Line) = π > angle((z-L.base)/L.direction) > 0

"""
	isright(z,L::Line)

Determine whether the number `z` lies "to the right" of line `L`. This means that the angle it makes with `tangent(L)` is in the interval (-π,0).

Note that `isleft` and `isright` are *not* logical opposites; a point on the curve should give `false` in both cases.
"""
isright(z::Number,L::Line) = -π < angle((z-L.base)/L.direction) < 0

winding(L::Line,z::Number) = isleft(z,L) ? 1 : 0

"""
	dist(z,L::Line)

Compute the distance from number `z` to the line `L`.
"""
dist(z::Number,L::Line) = imag( (z-L.base)/L.direction )
"""
	closest(z,L::Line)

Find the point on line `L` that lies closest to `z`.
"""
function closest(z::Number,L::Line)
	s = L.direction
	L.base + real((z-L.base)/s)*s
end
"""
	reflect(z,L::Line)

Reflect the value `z` across the line `L`. (For reflection of a line through a point, use translation and negation.)
"""
function reflect(z::Number,L::Line)
	ζ = z - L.base
	L.base + L.direction*conj(ζ/L.direction)
end

function show(io::IO,L::Line)
	print(IOContext(io,:compact=>true),"Line(...",L(0.5),"...",L((sqrt(5)-1)/2),"...)")
end
function show(io::IO,::MIME"text/plain",L::Line{T}) where {T}
	print(io,"Line{$T} in the complex plane:\n   through (",L.base,") parallel to (",L.direction,")")
end

# In the plane (but not on the sphere), two points are enough to draw a line, and we want to avoid infinity.
plotdata(L::Line{T}) where T<:Union{Complex,Polar} = point(L,[0.1,0.9])


text:
using Example
using Test

my_f(1, 2)

derivative(1, 2)

@testset "Example.jl" begin
    @test my_f(2, 3) == 13
    @test my_f(1, 2) == 8
end

@testset "derivative" begin
    @test derivative(1,2) == 2
end


text:
@testset "Multithread for solving steady states: single value ODE function" begin

atol = 1e-10
# Linear ODE which starts at 0.5 and solves from t=0.0 to t=1.0
de_func(u,p,t) = 0.9u*p

us = [1,2,3] # vec
p_ = 1
u_ss_real = zeros(length(us))

sim = solve_SSODE_threads(de_func, us, p_)
res = FindSteadyStates.flatten(get_sol2array(sim))

@test isapprox(res, u_ss_real; atol=atol)

end

@testset "Multithread for solving steady states: 2D ODE functions" begin

function f(du,u,p,t)
  du[1] = p[1] * u[1]
  du[2] = p[2] * u[2]
end

us = [[1.,2.], [3.,4.], [5.,6.] ]
p = [0.9,0.8]

ode = DEsteady(func=f, p = p , u0=us, method=SSRootfind())

@time sim = solve_SSODE_threads(f, us, p; method= Tsit5())
@time sim_ = solve_SSODE_threads(ode) #Rootfind method

@testset "compare solutions" for us_ in sim.u
      s = us_.u
      @test isapprox(s, [0.,0.]; atol=1e-10)
  end

@testset "compare solutions" for us_ in sim_.u
    s = us_.u
    @test isapprox(s, [0.,0.]; atol=1e-10)
end

end


text:
using Random
function jump_sensorloc(solver, seed, n; verbose = false, test = false)

    m, x_true, a, d, d_bar = sensorloc_data(seed, n)

    model = Model(with_optimizer(solver))

    # Build SDP problem
    @variable(model, X[1:n+2, 1:n+2], PSD)

    # Constraint with distances from anchors to sensors
    for j in 1:n, k in 1:m
        # e = zeros(n, 1)
        # e[j] = -1.0
        # v = vcat(a[k], e)
        # V = v * v'
        # @constraint(model, sum(V .* X) == d_bar[k, j]^2)
        @constraint(model, X[1,1]*a[k][1]*a[k][1] + X[2,2]*a[k][2]*a[k][2]
                                - 2 * X[1, j+2] * a[k][1]
                                - 2 * X[2, j+2] * a[k][2]
                                + X[j+2, j+2]
                            == d_bar[k, j]^2)
    end

    # Constraint with distances from sensors to sensors
    count, count_all = 0, 0
    rng = Random.MersenneTwister(seed)
    has_ctr = zeros(Bool,n,n)
    for i in 1:n, j in 1:i - 1
        count_all += 1
        if Random.rand(rng) > 0.9
            count += 1
            has_ctr[i,j] = true
            # e = zeros(n, 1)
            # e[i] = 1.0
            # e[j] = -1.0
            # v = vcat(zeros(2, 1), e)
            # V = v * v'
            # @constraint(model, sum(V .* X) == d[i, j]^2)
            @constraint(model, X[i+2,i+2] + X[j+2,j+2] - 2*X[i+2,j+2] == d[i, j]^2)
        end
    end
    if verbose
        @show count_all, count
    end

    @constraint(model, X[1, 1] == 1.0)
    @constraint(model, X[1, 2] == 0.0)
    @constraint(model, X[2, 1] == 0.0)
    @constraint(model, X[2, 2] == 1.0)

    # Feasibility objective function
    @objective(model, Min, 0.0 * X[1, 1] + 0.0 * X[2, 2])

    teste = @time optimize!(model)

    XX = value.(X)

    verbose && sensorloc_eval(n, m, x_true, XX)

    objval = objective_value(model)
    stime = MOI.get(model, MOI.SolveTime())

    # @show tp = typeof(model.moi_backend.optimizer.model.optimizer)
    # @show fieldnames(tp)
    rank = -1
    try
        @show rank = model.moi_backend.optimizer.model.optimizer.sol.final_rank
    catch
    end
    status = 0
    if JuMP.termination_status(model) == MOI.OPTIMAL
        status = 1
    end

    # SDP constraints
    max_spd_violation = minimum(eigen(XX).values)

    # test violations of linear constraints
    max_lin_viol = 0.0

    # Constraint with distances from anchors to sensors
    for j in 1:n
        for k in 1:m
            val =  abs(XX[1,1]*a[k][1]*a[k][1] + XX[2,2]*a[k][2]*a[k][2] -
                2 * XX[1, j+2] * a[k][1] -
                2 * XX[2, j+2] * a[k][2] +
                XX[j+2, j+2] -
                d_bar[k, j]^2)
            if val > 0.0
                if val > max_lin_viol
                    max_lin_viol = val
                end
            end
        end
    end

    # Constraint with distances from sensors to sensors
    count, count_all = 0, 0
    rng = Random.MersenneTwister(seed)
    has_ctr = zeros(Bool,n,n)
    for i in 1:n
        for j in 1:i - 1
            if has_ctr[i,j]
                val =  abs(XX[i+2,i+2] + XX[j+2,j+2] - 2*XX[i+2,j+2] - d[i, j]^2)
                if val > max_lin_viol
                    max_lin_viol = val
                end
            end
        end
    end
    max_lin_viol = max(max_lin_viol, abs(XX[1, 1] - 1.0))
    max_lin_viol = max(max_lin_viol, abs(XX[1, 2] - 0.0))
    max_lin_viol = max(max_lin_viol, abs(XX[2, 1] - 0.0))
    max_lin_viol = max(max_lin_viol, abs(XX[2, 2] - 1.0))

    return (objval, stime, rank, status, max_lin_viol, max_spd_violation)
    # return (objval, stime)
end


text:
module Autosampler

using CSV, DataFrames, ImagineFormat, ImagineInterface
using Random, Mmap

export stim_randomize, update_imagine

"""
    stim_randomize(filein, trials::Int, fileout=insertfilename(filein, "_sequence"))

Generate a randomize sequence of stimuli. `filein` is a string containing the name
of a CSV file that describes the unique stimuli (see format described in the README).
`trials` specifies the number of pseudo-random sequences used to generate the
final stimulus sequence.
`fileout` is an optional argument used to control the name of the output file;
the default is to insert "_sequence" right before the extension of `filein`.
"""
function stim_randomize(filein::AbstractString, trials::Int; fileout=insertfilename(filein, "_sequence"), kwargs...)
    open(fileout, "w") do io
        stim_randomize(io, filein, trials; kwargs...)
    end
    return fileout
end

function stim_randomize(io::IO, filein::AbstractString, trials::Int; header=false, kwargs...)
    df = CSV.File(filein) |> DataFrame
    println("Headers found: ", String.(names(df)))
    n = size(df, 1)
    println(n, " stimuli found")
    # $ is a disallowed character since we will use it as a separator
    for i = 1:n
        occursin('\$', df[i,1]) && error("\$ is disallowed in stimulus names")
    end
    p = Int[]
    for i = 1:trials
        append!(p, randperm(n))
    end
    CSV.write(io, df[p, :]; header=header, kwargs...)
    flush(io)
end

function update_imagine(imaginefile, sequencefile; um_per_pixel=nothing, aifile=replaceext(imaginefile, ".ai"), difile=replaceext(imaginefile, ".di"), updatedfile=imaginefile, csvheader=0, kwargs...)
    header = ImagineFormat.parse_header(imaginefile)
    haskey(header, "stimulus sequence") && error(imaginefile, " already has a `stimulus sequence` entry")
    if um_per_pixel === nothing
        um_per_pixel = header["um per pixel"]
    end
    um_per_pixel < 0 && error("must specify um_per_pixel")
    header["um per pixel"] = um_per_pixel

    df = CSV.File(sequencefile; header=csvheader, kwargs...) |> DataFrame
    n = size(df, 1)

    # Scan the AI file for stimulus triggers
    ai = parse_ai(aifile, header)
    aistim = getname(ai, "stimuli")
    stimhi  = find_pulse_starts(aistim; sampmap=:volts)
    stimlo  = find_pulse_stops(aistim; sampmap=:volts)
    if occursin("camera frame TTL", header["label list"])
        aiframe = getname(ai, "camera frame TTL")
        framestarts = find_pulse_starts(aiframe; sampmap=:volts)
    elseif occursin("camera1 frame monitor", header["di label list"])
        di = parse_di(difile, header)
        diframe = getname(di, "camera1 frame monitor")
        framestarts = find_pulse_starts(diframe) #; sampmap=:volts)
    end

    # Record as frameidx after the stimulus trigger
    fps = header["frames per stack"]
    framehi, framelo = Int[], Int[]
    for (frameidxs, scanidxs) in ((framehi, stimhi), (framelo, stimlo))
        for i in scanidxs
            idx = last(searchsorted(framestarts, i))
            push!(frameidxs, idx)
        end
    end

    # Add to header
    header["stimulus sequence"] = join(df[:,1], '\$')
    header["stimulus scan hi"] = stimhi
    header["stimulus scan lo"] = stimlo
    header["stimulus frame hi"] = framehi
    header["stimulus frame lo"] = framelo

    if updatedfile == imaginefile
        mv(imaginefile, imaginefile*".orig")
        sleep(0.25)
        isfile(imaginefile) && error("failed to move the old file")
    end
    ImagineFormat.save_header(updatedfile, header; misc=(
        "stimulus sequence", "stimulus scan hi", "stimulus scan lo",
        "stimulus frame hi", "stimulus frame lo"))
    return header
end

function insertfilename(filename, tail)
    basename, ext = splitext(filename)
    return basename*tail*ext
end

function replaceext(filename, ext)
    basename, _ = splitext(filename)
    return basename*ext
end

end # module


text:
function Dfp(arg0::Dfp)
    return Dfp((Dfp,), arg0)
end

function abs(obj::Dfp)
    return jcall(obj, "abs", Dfp, ())
end

function acos(obj::Dfp)
    return jcall(obj, "acos", Dfp, ())
end

function acosh(obj::Dfp)
    return jcall(obj, "acosh", Dfp, ())
end

function add(obj::Dfp, arg0::Dfp)
    return jcall(obj, "add", Dfp, (Dfp,), arg0)
end

function add(obj::Dfp, arg0::jdouble)
    return jcall(obj, "add", Dfp, (jdouble,), arg0)
end

function asin(obj::Dfp)
    return jcall(obj, "asin", Dfp, ())
end

function asinh(obj::Dfp)
    return jcall(obj, "asinh", Dfp, ())
end

function atan(obj::Dfp)
    return jcall(obj, "atan", Dfp, ())
end

function atan2(obj::Dfp, arg0::Dfp)
    return jcall(obj, "atan2", Dfp, (Dfp,), arg0)
end

function atanh(obj::Dfp)
    return jcall(obj, "atanh", Dfp, ())
end

function cbrt(obj::Dfp)
    return jcall(obj, "cbrt", Dfp, ())
end

function ceil(obj::Dfp)
    return jcall(obj, "ceil", Dfp, ())
end

function classify(obj::Dfp)
    return jcall(obj, "classify", jint, ())
end

function copy_sign(obj::Dfp, arg0::Dfp)
    return jcall(obj, "copySign", Dfp, (Dfp,), arg0)
end

function copy_sign(obj::Dfp, arg0::jdouble)
    return jcall(obj, "copySign", Dfp, (jdouble,), arg0)
end

function copysign(::Type{Dfp}, arg0::Dfp, arg1::Dfp)
    return jcall(Dfp, "copysign", Dfp, (Dfp, Dfp), arg0, arg1)
end

function cos(obj::Dfp)
    return jcall(obj, "cos", Dfp, ())
end

function cosh(obj::Dfp)
    return jcall(obj, "cosh", Dfp, ())
end

function divide(obj::Dfp, arg0::Dfp)
    return jcall(obj, "divide", Dfp, (Dfp,), arg0)
end

function divide(obj::Dfp, arg0::jdouble)
    return jcall(obj, "divide", Dfp, (jdouble,), arg0)
end

function divide(obj::Dfp, arg0::jint)
    return jcall(obj, "divide", Dfp, (jint,), arg0)
end

function dotrap(obj::Dfp, arg0::jint, arg1::JString, arg2::Dfp, arg3::Dfp)
    return jcall(obj, "dotrap", Dfp, (jint, JString, Dfp, Dfp), arg0, arg1, arg2, arg3)
end

function equals(obj::Dfp, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function exp(obj::Dfp)
    return jcall(obj, "exp", Dfp, ())
end

function expm1(obj::Dfp)
    return jcall(obj, "expm1", Dfp, ())
end

function floor(obj::Dfp)
    return jcall(obj, "floor", Dfp, ())
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_exponent(obj::Dfp)
    return jcall(obj, "getExponent", jint, ())
end

function get_field(obj::Dfp)
    return jcall(obj, "getField", DfpField, ())
end

function get_one(obj::Dfp)
    return jcall(obj, "getOne", Dfp, ())
end

function get_radix_digits(obj::Dfp)
    return jcall(obj, "getRadixDigits", jint, ())
end

function get_real(obj::Dfp)
    return jcall(obj, "getReal", jdouble, ())
end

function get_two(obj::Dfp)
    return jcall(obj, "getTwo", Dfp, ())
end

function get_zero(obj::Dfp)
    return jcall(obj, "getZero", Dfp, ())
end

function greater_than(obj::Dfp, arg0::Dfp)
    return jcall(obj, "greaterThan", jboolean, (Dfp,), arg0)
end

function hash_code(obj::Dfp)
    return jcall(obj, "hashCode", jint, ())
end

function hypot(obj::Dfp, arg0::Dfp)
    return jcall(obj, "hypot", Dfp, (Dfp,), arg0)
end

function int_log10(obj::Dfp)
    return jcall(obj, "intLog10", jint, ())
end

function int_value(obj::Dfp)
    return jcall(obj, "intValue", jint, ())
end

function is_infinite(obj::Dfp)
    return jcall(obj, "isInfinite", jboolean, ())
end

function is_na_n(obj::Dfp)
    return jcall(obj, "isNaN", jboolean, ())
end

function is_zero(obj::Dfp)
    return jcall(obj, "isZero", jboolean, ())
end

function less_than(obj::Dfp, arg0::Dfp)
    return jcall(obj, "lessThan", jboolean, (Dfp,), arg0)
end

function linear_combination(obj::Dfp, arg0::Dfp, arg1::Dfp, arg2::Dfp, arg3::Dfp)
    return jcall(obj, "linearCombination", Dfp, (Dfp, Dfp, Dfp, Dfp), arg0, arg1, arg2, arg3)
end

function linear_combination(obj::Dfp, arg0::Dfp, arg1::Dfp, arg2::Dfp, arg3::Dfp, arg4::Dfp, arg5::Dfp)
    return jcall(obj, "linearCombination", Dfp, (Dfp, Dfp, Dfp, Dfp, Dfp, Dfp), arg0, arg1, arg2, arg3, arg4, arg5)
end

function linear_combination(obj::Dfp, arg0::Dfp, arg1::Dfp, arg2::Dfp, arg3::Dfp, arg4::Dfp, arg5::Dfp, arg6::Dfp, arg7::Dfp)
    return jcall(obj, "linearCombination", Dfp, (Dfp, Dfp, Dfp, Dfp, Dfp, Dfp, Dfp, Dfp), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function linear_combination(obj::Dfp, arg0::Vector{Dfp}, arg1::Vector{Dfp})
    return jcall(obj, "linearCombination", Dfp, (Vector{Dfp}, Vector{Dfp}), arg0, arg1)
end

function linear_combination(obj::Dfp, arg0::Vector{jdouble}, arg1::Vector{Dfp})
    return jcall(obj, "linearCombination", Dfp, (Vector{jdouble}, Vector{Dfp}), arg0, arg1)
end

function linear_combination(obj::Dfp, arg0::jdouble, arg1::Dfp, arg2::jdouble, arg3::Dfp)
    return jcall(obj, "linearCombination", Dfp, (jdouble, Dfp, jdouble, Dfp), arg0, arg1, arg2, arg3)
end

function linear_combination(obj::Dfp, arg0::jdouble, arg1::Dfp, arg2::jdouble, arg3::Dfp, arg4::jdouble, arg5::Dfp)
    return jcall(obj, "linearCombination", Dfp, (jdouble, Dfp, jdouble, Dfp, jdouble, Dfp), arg0, arg1, arg2, arg3, arg4, arg5)
end

function linear_combination(obj::Dfp, arg0::jdouble, arg1::Dfp, arg2::jdouble, arg3::Dfp, arg4::jdouble, arg5::Dfp, arg6::jdouble, arg7::Dfp)
    return jcall(obj, "linearCombination", Dfp, (jdouble, Dfp, jdouble, Dfp, jdouble, Dfp, jdouble, Dfp), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function log(obj::Dfp)
    return jcall(obj, "log", Dfp, ())
end

function log10(obj::Dfp)
    return jcall(obj, "log10", Dfp, ())
end

function log10k(obj::Dfp)
    return jcall(obj, "log10K", jint, ())
end

function log1p(obj::Dfp)
    return jcall(obj, "log1p", Dfp, ())
end

function multiply(obj::Dfp, arg0::Dfp)
    return jcall(obj, "multiply", Dfp, (Dfp,), arg0)
end

function multiply(obj::Dfp, arg0::jdouble)
    return jcall(obj, "multiply", Dfp, (jdouble,), arg0)
end

function multiply(obj::Dfp, arg0::jint)
    return jcall(obj, "multiply", Dfp, (jint,), arg0)
end

function negate(obj::Dfp)
    return jcall(obj, "negate", Dfp, ())
end

function negative_or_null(obj::Dfp)
    return jcall(obj, "negativeOrNull", jboolean, ())
end

function new_instance(obj::Dfp)
    return jcall(obj, "newInstance", Dfp, ())
end

function new_instance(obj::Dfp, arg0::Dfp)
    return jcall(obj, "newInstance", Dfp, (Dfp,), arg0)
end

function new_instance(obj::Dfp, arg0::DfpField, arg1::DfpField_RoundingMode)
    return jcall(obj, "newInstance", Dfp, (DfpField, DfpField_RoundingMode), arg0, arg1)
end

function new_instance(obj::Dfp, arg0::JString)
    return jcall(obj, "newInstance", Dfp, (JString,), arg0)
end

function new_instance(obj::Dfp, arg0::jbyte)
    return jcall(obj, "newInstance", Dfp, (jbyte,), arg0)
end

function new_instance(obj::Dfp, arg0::jbyte, arg1::jbyte)
    return jcall(obj, "newInstance", Dfp, (jbyte, jbyte), arg0, arg1)
end

function new_instance(obj::Dfp, arg0::jdouble)
    return jcall(obj, "newInstance", Dfp, (jdouble,), arg0)
end

function new_instance(obj::Dfp, arg0::jint)
    return jcall(obj, "newInstance", Dfp, (jint,), arg0)
end

function new_instance(obj::Dfp, arg0::jlong)
    return jcall(obj, "newInstance", Dfp, (jlong,), arg0)
end

function next_after(obj::Dfp, arg0::Dfp)
    return jcall(obj, "nextAfter", Dfp, (Dfp,), arg0)
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function positive_or_null(obj::Dfp)
    return jcall(obj, "positiveOrNull", jboolean, ())
end

function pow(obj::Dfp, arg0::Dfp)
    return jcall(obj, "pow", Dfp, (Dfp,), arg0)
end

function pow(obj::Dfp, arg0::jdouble)
    return jcall(obj, "pow", Dfp, (jdouble,), arg0)
end

function pow(obj::Dfp, arg0::jint)
    return jcall(obj, "pow", Dfp, (jint,), arg0)
end

function power10(obj::Dfp, arg0::jint)
    return jcall(obj, "power10", Dfp, (jint,), arg0)
end

function power10k(obj::Dfp, arg0::jint)
    return jcall(obj, "power10K", Dfp, (jint,), arg0)
end

function reciprocal(obj::Dfp)
    return jcall(obj, "reciprocal", Dfp, ())
end

function remainder(obj::Dfp, arg0::Dfp)
    return jcall(obj, "remainder", Dfp, (Dfp,), arg0)
end

function remainder(obj::Dfp, arg0::jdouble)
    return jcall(obj, "remainder", Dfp, (jdouble,), arg0)
end

function rint(obj::Dfp)
    return jcall(obj, "rint", Dfp, ())
end

function root_n(obj::Dfp, arg0::jint)
    return jcall(obj, "rootN", Dfp, (jint,), arg0)
end

function round(obj::RealFieldElement)
    return jcall(obj, "round", jlong, ())
end

function scalb(obj::Dfp, arg0::jint)
    return jcall(obj, "scalb", Dfp, (jint,), arg0)
end

function signum(obj::Dfp)
    return jcall(obj, "signum", Dfp, ())
end

function sin(obj::Dfp)
    return jcall(obj, "sin", Dfp, ())
end

function sin_cos(obj::CalculusFieldElement)
    return jcall(obj, "sinCos", FieldSinCos, ())
end

function sinh(obj::Dfp)
    return jcall(obj, "sinh", Dfp, ())
end

function sqrt(obj::Dfp)
    return jcall(obj, "sqrt", Dfp, ())
end

function strictly_negative(obj::Dfp)
    return jcall(obj, "strictlyNegative", jboolean, ())
end

function strictly_positive(obj::Dfp)
    return jcall(obj, "strictlyPositive", jboolean, ())
end

function subtract(obj::Dfp, arg0::Dfp)
    return jcall(obj, "subtract", Dfp, (Dfp,), arg0)
end

function subtract(obj::Dfp, arg0::jdouble)
    return jcall(obj, "subtract", Dfp, (jdouble,), arg0)
end

function tan(obj::Dfp)
    return jcall(obj, "tan", Dfp, ())
end

function tanh(obj::Dfp)
    return jcall(obj, "tanh", Dfp, ())
end

function to_degrees(obj::Dfp)
    return jcall(obj, "toDegrees", Dfp, ())
end

function to_double(obj::Dfp)
    return jcall(obj, "toDouble", jdouble, ())
end

function to_radians(obj::Dfp)
    return jcall(obj, "toRadians", Dfp, ())
end

function to_split_double(obj::Dfp)
    return jcall(obj, "toSplitDouble", Vector{jdouble}, ())
end

function to_string(obj::Dfp)
    return jcall(obj, "toString", JString, ())
end

function unequal(obj::Dfp, arg0::Dfp)
    return jcall(obj, "unequal", jboolean, (Dfp,), arg0)
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end