run.jl

using Pkg
Pkg.activate(".")

using Revise, Statistics, Dierckx, DataFrames
using QSFit, GFit, Gnuplot, GFitViewer, MyAstroUtils
using CL_1ES_1927p654
using Dates
using TextParse
import StatsBase, JLD2

if readlines(`hostname`)[1] != "giorgio-pc"
    Gnuplot.options.term="sixelgd size 800,600"
end

function read_epochs()
    out = DataFrame(id=Int[], filename=String[], date=String[], scale=Float64[], OIII_fwhm=Float64[])
    for (root, dirs, files) in walkdir("AT2018zf")
        for file in files
            if file[end-3:end] == ".dat"
                filename = root * "/" * file
                date = filename[27:end-4]
                push!(out, (nrow(out)+1, filename, date, NaN, NaN))
            end
        end
    end
    @assert issorted(out.date)
    return out
end

function read_spec(id; kw...)
    if id == "B03"
        file = "boller2003.txt"
    else
        irow = findfirst(epochs.id .== id)
        file = epochs[irow, :filename]
    end
    println(file)
	spec = Spectrum(Val(:ASCII), file, columns=[1,2]; label="$id: $file", kw...);
    if isa(id, Number)
        spec.flux ./= epochs[irow, :scale]
        spec.err  ./= epochs[irow, :scale]
    end
    return spec
end


function analyze_boller03()
    file = opt.path * "/boller.dat"
    if !isfile(file)
        source = QSO{q1927p654}("1ES 1927+654 (Boller03)", opt.z, ebv=opt.ebv);
        source.options[:min_spectral_coverage][:OIII_5007] = 0.35
        # source.options[:host_template] = "/home/gcalderone/Mbi1.30Zm1.49T00.0700_iTp0.00_baseFe_linear_FWHM_2.51"
        #=
        Resolution: from  Sect 3 it is  6 / 5007 * 3e5 ~ 360 km/s
        however with this value the [OIII] FWHM is 100 km/s (limit)
        I tried 180, 200 and 250, and the [OIII] luminosity do not change
        =#
        spec = read_spec("B03", resolution=200.)
        add_spec!(source, spec);
        res = qsfit(source);
        viewer(res, showcomps=[:qso_cont, :galaxy, :balmer],
               filename="$(opt.path)/results_boller.html")
        JLD2.save_object(file, res)
    else
        res = JLD2.load_object(file)
    end
    @assert res.model[:OIII_5007].norm.val == 0.00029042977121756283
    @assert res.model[:OIII_5007].norm.val / res.model[:galaxy].norm.val == 60.07373299362222
    return res
end


function estimate_SE_scale!()
    file = opt.path * "/" * "SE_scale.dat"
    if !isfile(file)
        for irow in 1:nrow(epochs)
            @info irow
            id = epochs[irow, :id]
            epochs[irow, :scale] = 1.
            spec = read_spec(id)
            epochs[irow, :scale] = median(spec.flux) / 25000

            source = QSO{q1927p654}("1ES 1927+654 ($id)", opt.z, ebv=opt.ebv);
            source.options[:min_spectral_coverage][:OIII_5007] = 0.5
            # source.options[:host_template] = "/home/gcalderone/Mbi1.30Zm1.49T00.0700_iTp0.00_baseFe_linear_FWHM_2.51"
            spec = read_spec(id)
            add_spec!(source, spec);
            res = qsfit(source);
            viewer(res, showcomps=[:qso_cont, :galaxy, :balmer],
                   filename="$(opt.path)/results_$(id).html")

            # Update all_scale and store FWHM
            @info res.model[:OIII_5007].norm.val
            epochs[irow, :scale] *= res.model[:OIII_5007].norm.val
            epochs[irow, :OIII_fwhm] = res.model[:OIII_5007].fwhm.val
        end
        @assert all(epochs[.!isfinite.(epochs.OIII_fwhm), :scale] .== NaN)
        JLD2.save_object(file, epochs)
    else
        empty!(epochs)
        append!(epochs, JLD2.load_object(file))
    end
end


function analyze_job(job; Nloop = 6)
    out = fill(NaN, length(epoch_ids[job]))

    estimate_SE_scale!()
    for loop in 1:Nloop
        file = "$(opt.path)/results_$(job)_$(loop).dat"
        if !isfile(file)
            source = QSO{q1927p654}("1ES 1927+654", opt.z, ebv=opt.ebv);
            source.options[:min_spectral_coverage][:OIII_5007] = 0.5
            source.options[:wavelength_range] = [3500, 6900]
            @gp :zoom "set grid" :-
            for id in epoch_ids[job]
                spec = read_spec(id, resolution=get(resolution, job, NaN))
                add_spec!(source, spec);
                @gp :- :zoom xr=[4750,5150] spec.λ ./ (1 + source.z) spec.flux "w l t '$(id)'"
            end
            res = qsfit_multi(source);
            JLD2.save_object(file, res)
            viewer(res, showcomps=[:qso_cont, :galaxy, :balmer], filename="$(opt.path)/results_$(job)_$(loop).html")
            viewer(res, showcomps=[:qso_cont, :galaxy, :balmer], filename="$(opt.path)/results_$(job)_$(loop)_rebin4.html", rebin=4)
        else
            @info file
            res = JLD2.load_object(file);
            println(loop, "  ", res.fitres.gofstat / res.fitres.dof, "  ", res.fitres.elapsed / 3600.)
        end

        OIII_norm = Float64[]
        for i in 1:length(epoch_ids[job])
            id = epoch_ids[job][i]
            push!(OIII_norm, res.multi[i][:OIII_5007].norm.val)
            epochs[epochs.id .== id, :scale] *= ((1 + res.multi[i][:OIII_5007].norm.val) / 2)
        end
        out = hcat(out, OIII_norm)
    end
    out = out[:,2:end]

    @gp "set grid" :-
    for i in 1:length(epoch_ids[job])
        id = epoch_ids[job][i]
        @gp :- 1:Nloop out[i, :] "w lp t '$(id)'" :-
    end
    @gp

    return out
end




const opt = (
z = 0.01942,
ebv = 0.077,
path = "output")

try;  mkdir(opt.path);  catch; end
const epochs = read_epochs()
analyze_boller03()
estimate_SE_scale!()

@gp :prenorm "set key bottom right" "set grid" xlab="Epoch" ylab="Value" :-
@gp :- :prenorm epochs.scale./1e-18 "w lp t '[OIII] norm.'" :-
@gp :- :prenorm epochs.OIII_fwhm./1000 "w lp t '[OIII] FWHM (x1000 km/s)"
save(:prenorm, term="png size 800,600", output="$(opt.path)/oiii_norm_fwhm.png")

w = 10 .^(2:0.01:log10(400));
r = sqrt.(400^2 .- w.^2) ./ 2.355;
@gp w r "w l t '400 km/s'"

w = 10 .^(2:0.01:log10(900));
r = sqrt.(900^2 .- w.^2) ./ 2.355;
@gp :- w r "w l t '900 km/s'"
@gp :- "set grid" xlab="Actual FWHM" ylab="Instr. resolution"


# Ensure no sample is negative
@gp :allepochs "set grid" xr=[3e3,1e4] cbr=[1,29] :-
@gp :- :allepochs cblabel="Epoch" xlab="Rest frame wavelength[A]" ylab="Flux density [arb.units]"
for id in epochs.id
	s = read_spec(id)
	x = s.λ;
	y = s.flux;
    y .-= median(y)
    y ./= maximum(y)
    @gp :- :allepochs x ./ (1 + opt.z) y "u 1:(\$2+$id):($id) w l notit lc pal" :-
end
@gp :- :allepochs yr=[1, 29]
save(:allepochs, term="png size 800,600", output="$(opt.path)/all_epochs.png")
@gp :- :allepochs xr=[4900, 5100]
save(:allepochs, term="png size 800,2000", output="$(opt.path)/all_epochs_zoom.png")



epoch_ids = Dict(
    :vecchi   => [17, 10, 21, 23],
    :lowres   => [8,9,10,11,12,14,15,16,18,21,22,23,24,25],
    :highres  => [5, 7, 13, 17, 20],
    :test     => [7, 13],
    :limited  => [1,2,3,13],
    :all      => collect(5:26),
    :floyds   => [6, 8, 9, 10, 11, 12, 14, 15, 16, 18, 19, 21, 22, 23])
# deleteat!(epoch_ids[:all], 4)  #insufficient coverage in na_Hg

# Can't use resolution smaller than 150 km / s otherwise some line
# will be neglected because of spectral coverage
resolution = Dict()
#    :lowres  => 350.,
#    :highres => 150.)

analyze_job(:test)
analyze_job(:highres)
analyze_job(:all)
analyze_job(:floyds)
ddd



function cont_at(bestfit, λ)
    A1 = bestfit[:qso_cont].norm.val - bestfit[:qso_cont].norm.unc
    A  = bestfit[:qso_cont].norm.val
    A2 = bestfit[:qso_cont].norm.val + bestfit[:qso_cont].norm.unc
    if A1 > A2
        A1, A2 = A2, A1
    end
    @assert A1 <= A2

    B1 = (λ / bestfit[:qso_cont].x0.val)^(bestfit[:qso_cont].alpha.val - bestfit[:qso_cont].alpha.unc)
    B  = (λ / bestfit[:qso_cont].x0.val)^(bestfit[:qso_cont].alpha.val)
    B2 = (λ / bestfit[:qso_cont].x0.val)^(bestfit[:qso_cont].alpha.val + bestfit[:qso_cont].alpha.unc)
    if B1 > B2
        B1, B2 = B2, B1
    end
    @assert B1 <= B2

    dd = extrema([A1*B1, A1*B2, A2*B1, A2*B2])

    return λ * A * B, λ * (dd[2]-dd[1])
end

# Wrap-up
chosen_job = :all
chosen_iloop = 5

res = JLD2.load_object("output/results_$(chosen_job)_$(chosen_iloop).dat")

tab = DataFrame(epoch=Int[], date=String[], instr=String[], galaxy=Float64[],
                contslope=Float64[],
                c3500=Float64[], c5100=Float64[],
                Ha_bb_norm=Float64[], Ha_bb_fwhm=Float64[], Ha_bb_voff=Float64[],
                Ha_br_norm=Float64[], Ha_br_fwhm=Float64[], Ha_br_voff=Float64[],
                Ha_na_norm=Float64[], Ha_na_fwhm=Float64[], Ha_na_voff=Float64[],
                Hb_bb_norm=Float64[], Hb_bb_fwhm=Float64[], Hb_bb_voff=Float64[],
                Hb_br_norm=Float64[], Hb_br_fwhm=Float64[], Hb_br_voff=Float64[],
                Hb_na_norm=Float64[], Hb_na_fwhm=Float64[], Hb_na_voff=Float64[],
                oiii_fwhm=Float64[],  oiii_voff=Float64[])

for id in 1:length(res.multi)
    push!(tab,
          (epoch_ids[chosen_job][id], string(split(res.source.specs[id].label, "/")[2])[18:end-4], "",
           res.multi[id][:galaxy].norm.val,
           res.multi[id][:qso_cont].alpha.val,
           cont_at(res.multi[id], 3500)[1],
           cont_at(res.multi[id], 5100)[1],
           res.multi[id][:Ha_bb].norm.patched, res.multi[id][:Ha_bb].fwhm.patched, res.multi[id][:Ha_bb].voff.patched,
           res.multi[id][:Ha_br].norm.patched, res.multi[id][:Ha_br].fwhm.patched, res.multi[id][:Ha_br].voff.patched,
           res.multi[id][:Ha_na].norm.patched, res.multi[id][:Ha_na].fwhm.patched, res.multi[id][:Ha_na].voff.patched,
           res.multi[id][:Hb_bb].norm.patched, res.multi[id][:Hb_bb].fwhm.patched, res.multi[id][:Hb_bb].voff.patched,
           res.multi[id][:Hb_br].norm.patched, res.multi[id][:Hb_br].fwhm.patched, res.multi[id][:Hb_br].voff.patched,
           res.multi[id][:Hb_na].norm.patched, res.multi[id][:Hb_na].fwhm.patched, res.multi[id][:Hb_na].voff.patched,
           res.multi[id][:OIII_5007].fwhm.patched, res.multi[id][:OIII_5007].voff.patched))
end

tab[!, :l3500] .= [3500 .* Spline1D(domain(res.multi[id])[:], res.multi[id]())(3500.) for id in 1:length(res.multi)]
tab[!, :l5100] .= [5100 .* Spline1D(domain(res.multi[id])[:], res.multi[id]())(5100.) for id in 1:length(res.multi)]
tab[!, :Ha] .= tab.Ha_br_norm .+ tab.Ha_bb_norm .+ tab.Ha_na_norm
tab[!, :Hb] .= tab.Hb_br_norm .+ tab.Hb_bb_norm .+ tab.Hb_na_norm

dd = Date.(tab.date, Ref("yyyymmdd"))
tab[!, :day] .= getproperty.(dd - Date("2018-03-06") .+ Day(72), :value)


(d, c) = csvread("misc/tabula-2019-Trakhtenbrot-mod.csv", ',', header_exists=true)
tmp = DataFrame(collect(d), Symbol.(c))
tmp[!, :Date] .= string.(tmp.Date)

tab[!, :instr] .= ""
for i in 1:nrow(tab)
    j = findfirst(tmp.Date .== tab[i, :date])
    isnothing(j)  &&  continue
    tab[i, :instr] = tmp[j, :Inst]
end

f = FITS("$(opt.path)/1ES_1927p654_results.fits", "w")
write(f, tab)
close(f)


tab[!, :pt] .= 0
for i in 1:nrow(tab)
    j = findfirst(sort(unique(tab.instr)) .== tab[i, :instr])
    isnothing(j)  &&  continue
    tab[i, :pt] = j
end
cm = StatsBase.countmap(tab.instr)
setindex!.(Ref(cm), 1:length(cm), collect(keys(cm))[sortperm(collect(values(cm)), rev=true)])
tab.pt .= getindex.(Ref(cm), tab.instr)


@gp    tab.contslope tab.Ha    tab.pt "w lp t 'Ha'    pt var ps 2"  ylog=true
@gp :- tab.contslope tab.Hb    tab.pt "w lp t 'Hb'    pt var ps 2"
@gp :- tab.contslope tab.c5100 tab.pt "w lp t 'c5100' pt var ps 2"
@gp :- tab.contslope tab.l5100 tab.pt "w lp t 'l5100' pt var ps 2"
@gp :- tab.contslope tab.c3500 tab.pt "w lp t 'c3500' pt var ps 2"
@gp :- tab.contslope tab.l3500 tab.pt "w lp t 'l3500' pt var ps 2"


@gp    xlog=true ylog=true
@gp :- tab.c5100 tab.Ha tab.pt "w lp t 'c5100' pt var ps 2"
@gp :- tab.l5100 tab.Ha tab.pt "w lp t 'l5100' pt var ps 2"


@gp    tab.day tab.c5100 tab.pt "w lp t 'c5100' pt var ps 2" ylog=true
@gp :- tab.day tab.l5100 tab.pt "w lp t 'l5100' pt var ps 2"
@gp :- tab.day tab.Ha    tab.pt "w lp t 'Ha'    pt var ps 2"
@gp :- tab.day tab.Hb    tab.pt "w lp t 'Hb'    pt var ps 2"

@gp  tab.day tab.contslope tab.pt "w lp t 'SLope' pt var ps 2"


@gp "set grid"
for i in 1:nrow(tab)
    ss = "w l t '$(tab[i, :date])' dt $(tab[i, :pt]) lw 2"
    @gp :- domain(res.multi[i])[:] res.multi[i](:qso_cont) ss
end
save(:prenorm, term="png size 800,600", output="$(opt.path)/evolution.png")


@gp "set grid"
for i in 1:length(res.multi)
    println(res.multi[i][:OIII_5007])
    x  = domain(res.multi[i])[:]
    y  = res.multi[i]()
    y0 = y .- res.multi[i](:OIII_5007)
    ss = "w l t '$(tab[i, :date])' dt $(tab[i, :pt]) lw 2"
    @gp :- x y .- Spline1D(x, y0)(5007.) ss
end
save(:prenorm, term="png size 800,600", output="$(opt.path)/evolution_oiii_norm.png")