Recipe_single.jl

function minimizer(source::QSO{T}) where T <: DefaultRecipe
    mzer = GFit.cmpfit()
    mzer.Δfitstat_theshold = 1.e-5
    return mzer
end


function add_qso_continuum!(source::QSO{T}, model::Model; id=1) where T <: DefaultRecipe
    λ = domain(model)[:]

    comp = QSFit.powerlaw(3000)
    comp.x0.val = median(λ)
    comp.norm.val = Spline1D(λ, source.data[id].val, k=1, bc="error")(comp.x0.val)
    comp.alpha.val  = -1.5
    comp.alpha.low  = -3
    comp.alpha.high =  1
    
    model[:qso_cont] = comp
    model[:Continuum] = SumReducer([:qso_cont])
end


function add_host_galaxy!(source::QSO{T}, model::Model; id=1) where T <: DefaultRecipe
    λ = domain(model)[:]
    if source.options[:use_host_template]  &&
        (minimum(λ) .< 5500 .< maximum(λ))
        model[:galaxy] = QSFit.hostgalaxy(source.options[:host_template])
        model[:Continuum] = SumReducer([:qso_cont, :galaxy])

        # Split total flux between continuum and host galaxy
        vv = Spline1D(λ, source.data[id].val, k=1, bc="error")(5500.)
        model[:galaxy].norm.val  = 1/2 * vv
        model[:qso_cont].x0.val *= 1/2 * vv / Spline1D(λ, model(:qso_cont), k=1, bc="error")(5500.)
    end
end


function add_balmer_cont!(source::QSO{T}, model::Model; id=1) where T <: DefaultRecipe
    if source.options[:use_balmer]
        tmp = [:qso_cont, :balmer]
        (:galaxy in keys(model))  &&  push!(tmp, :galaxy)
        model[:balmer] = QSFit.balmercont(0.1, 0.5)
        model[:Continuum] = SumReducer(tmp)
        c = model[:balmer]
        c.norm.val  = 0.1
        c.norm.fixed = false
        c.norm.high = 0.5
        c.ratio.val = 0.5
        c.ratio.fixed = false
        c.ratio.low  = 0.1
        c.ratio.high = 1
        @patch! model[:balmer].norm *= model[:qso_cont].norm
    end
end


function renorm_cont!(source::QSO{T}, model::Model; id=1) where T <: DefaultRecipe
    freeze(model, :qso_cont)
    c = model[:qso_cont]
    initialnorm = c.norm.val
    if c.norm.val > 0
        println(logio(source), "Cont. norm. (before): ", c.norm.val)
        while true
            residuals = (model() - source.data[id].val) ./ source.data[id].unc
            (count(residuals .< 0) / length(residuals) > 0.9)  &&  break
            (c.norm.val < initialnorm / 5)  &&  break # give up
            c.norm.val *= 0.99
            evaluate!(model)
        end
        println(logio(source), "Cont. norm. (after) : ", c.norm.val)
    else
        println(logio(source), "Skipping cont. renormalization")
    end
end


function add_iron_uv!(source::QSO{T}, model::Model; id=1) where T <: DefaultRecipe
    λ = domain(model)[:]
    if source.options[:use_ironuv]
        fwhm = 3000.
        source.options[:instr_broadening]  ||  (fwhm = sqrt(fwhm^2 + (source.spectra[id].resolution * 2.355)^2))
        comp = QSFit.ironuv(fwhm)
        (_1, _2, coverage) = spectral_coverage(λ, source.spectra[id].resolution, comp)
        threshold = get(source.options[:min_spectral_coverage], :ironuv, source.options[:min_spectral_coverage][:default])
        if coverage >= threshold
            model[:ironuv] = comp
            model[:ironuv].norm.val = 0.5
            push!(model[:Iron].list, :ironuv)
        else
            println(logio(source), "Ignoring ironuv component (threshold: $threshold)")
        end
    end
end


function add_iron_opt!(source::QSO{T}, model::Model; id=1) where T <: DefaultRecipe
    λ = domain(model)[:]
    if source.options[:use_ironopt]
        fwhm = 3000.
        source.options[:instr_broadening]  ||  (fwhm = sqrt(fwhm^2 + (source.spectra[id].resolution * 2.355)^2))
        comp = QSFit.ironopt_broad(fwhm)
        (_1, _2, coverage) = spectral_coverage(λ, source.spectra[id].resolution, comp)
        threshold = get(source.options[:min_spectral_coverage], :ironopt, source.options[:min_spectral_coverage][:default])
        if coverage >= threshold
            fwhm = 500.
            source.options[:instr_broadening]  ||  (fwhm = sqrt(fwhm^2 + (source.spectra[id].resolution * 2.355)^2))
            model[:ironoptbr] = comp
            model[:ironoptna] = QSFit.ironopt_narrow(fwhm)
            model[:ironoptbr].norm.val = 0.1  # TODO: guess a sensible value
            model[:ironoptna].norm.val = 0.0
            freeze(model, :ironoptna)  # will be freed during last run
            push!(model[:Iron].list, :ironoptbr)
            push!(model[:Iron].list, :ironoptna)
        else
            println(logio(source), "Ignoring ironopt component (threshold: $threshold)")
        end
    end
end


function add_emission_lines!(source::QSO{T}, model::Model; id=1) where T <: DefaultRecipe
    line_groups = unique(collect(values(source.line_names[id])))
    for (cname, comp) in source.line_comps[id]
        comp.norm_integrated = source.options[:norm_integrated]
        model[cname] = comp
    end
    for (group, lnames) in QSFit.invert_dictionary(source.line_names[id])
        model[group] = SumReducer(lnames)
    end
    model[:main] = SumReducer([:Continuum, :Iron, line_groups...])

    if haskey(model, :MgII_2798)
        model[:MgII_2798].voff.low  = -1000
        model[:MgII_2798].voff.high =  1000
    end
    if haskey(model, :OIII_5007_bw)
        model[:OIII_5007_bw].fwhm.val  = 500
        model[:OIII_5007_bw].fwhm.low  = 1e2
        model[:OIII_5007_bw].fwhm.high = 1e3
        model[:OIII_5007_bw].voff.low  = 0
        model[:OIII_5007_bw].voff.high = 2e3
    end
    evaluate!(model)
end


function guess_emission_lines_values!(source::QSO{T}, model::Model; id=1) where T <: DefaultRecipe
    resid = source.data[id].val - model(:Continuum) .- model(:Iron)
    for cname in keys(source.line_names[id])
        c = model[cname]
        resid_at_line = Spline1D(λ, resid, k=1, bc="nearest")(c.center.val)
        c.norm.val *= abs(resid_at_line) / maximum(model(cname))

        # If instrumental broadening is not used and the line profile
        # is a Gaussian one take spectral resolution into account.
        # This is significantly faster than convolving with an
        # instrument response but has some limitations:
        # - works only with Gaussian profiles;
        # - all components must be additive (i.e. no absorptions)
        # - further narrow components (besides known emission lines)
        #   will not be corrected for instrumental resolution
        if !source.options[:instr_broadening]
            if isa(c, QSFit.SpecLineGauss)
                c.spec_res_kms = source.spectra[id].resolution
            else
                println(logio(source), "Line $cname is not a Gaussian profile: Can't take spectral resolution into account")
            end
        end
    end
    evaluate!(model)
end

    
function add_patch_functs!(source::QSO{T}, model::Model; id=1) where T <: DefaultRecipe
    # Patch parameters
    @patch! begin
         # model[:OIII_4959].norm = model[:OIII_5007].norm / 3
         model[:OIII_4959].voff = model[:OIII_5007].voff
     end
    @patch! begin
        model[:OIII_5007_bw].voff += model[:OIII_5007].voff
        model[:OIII_5007_bw].fwhm += model[:OIII_5007].fwhm
    end
    @patch! begin
        # model[:OI_6300].norm = model[:OI_6364].norm / 3
        model[:OI_6300].voff = model[:OI_6364].voff
    end
    @patch! begin
        # model[:NII_6549].norm = model[:NII_6583].norm / 3
        model[:NII_6549].voff = model[:NII_6583].voff
    end
    @patch! begin
        # model[:SII_6716].norm = model[:SII_6731].norm / 1.5
        model[:SII_6716].voff = model[:SII_6731].voff
    end

    @patch! model[:na_Hb].voff = model[:na_Ha].voff

    # The following are required to avoid degeneracy with iron
    # template
    @patch! begin
        model[:Hg].voff = model[:br_Hb].voff
        model[:Hg].fwhm = model[:br_Hb].fwhm
    end
    @patch! begin
        model[:br_Hg].voff = model[:br_Hb].voff
        model[:br_Hg].fwhm = model[:br_Hb].fwhm
    end
    @patch! begin
        model[:na_Hg].voff = model[:na_Hb].voff
        model[:na_Hg].fwhm = model[:na_Hb].fwhm
    end

    # Ensure luminosity at peak of the broad base component is
    # smaller than the associated broad component:
    if  haskey(model, :br_Hb)  &&
        haskey(model, :bb_Hb)
        model[:bb_Hb].norm.high = 1
        model[:bb_Hb].norm.val  = 0.5
        @patch! model[:bb_Hb].norm *= model[:br_Hb].norm / model[:br_Hb].fwhm * model[:bb_Hb].fwhm
    end
    if  haskey(model, :br_Ha)  &&
        haskey(model, :bb_Ha)
        model[:bb_Ha].norm.high = 1
        model[:bb_Ha].norm.val  = 0.5
        @patch! model[:bb_Ha].norm *= model[:br_Ha].norm / model[:br_Ha].fwhm * model[:bb_Ha].fwhm
    end
end

    
function add_unknown_lines!(source::QSO{T}, model::Model; id=1) where T <: DefaultRecipe
    (source.options[:n_unk] > 0)  &&  (return nothing)
    tmp = OrderedDict{Symbol, GFit.AbstractComponent}()
    for j in 1:source.options[:n_unk]
        tmp[Symbol(:unk, j)] = line_component(TRecipe, QSFit.UnkLine(5e3))
        tmp[Symbol(:unk, j)].norm_integrated = source.options[:norm_integrated]
    end
    for (cname, comp) in tmp
        model[cname] = comp
    end
    model[:UnkLines] = SumReducer(collect(keys(tmp)))
    model[:main] = SumReducer([:Continuum, :Iron, line_groups..., :UnkLines])
    evaluate!(model)
    for j in 1:source.options[:n_unk]
        freeze(model, Symbol(:unk, j))
    end
    evaluate!(model)

    # Set "unknown" line center wavelength where there is a maximum in
    # the fit residuals, and re-run a fit.
    λunk = Vector{Float64}()
    while true
        (length(λunk) >= source.options[:n_unk])  &&  break
        evaluate!(model)
        Δ = (source.data[id].val - model()) ./ source.data[id].unc

        # Avoid considering again the same region (within 1A) TODO: within resolution
        for l in λunk
            Δ[findall(abs.(l .- λ) .< 1)] .= 0.
        end

        # Avoidance regions
        for rr in source.options[:unk_avoid]
            Δ[findall(rr[1] .< λ .< rr[2])] .= 0.
        end

        # Do not add lines close to from the edges since these may
        # affect qso_cont fitting
        Δ[findall((λ .< minimum(λ)*1.02)  .|
                  (λ .> maximum(λ)*0.98))] .= 0.
        iadd = argmax(Δ)
        (Δ[iadd] <= 0)  &&  break  # No residual is greater than 0, skip further residuals....
        push!(λunk, λ[iadd])

        cname = Symbol(:unk, length(λunk))
        model[cname].norm.val = 1.
        model[cname].center.val  = λ[iadd]
        model[cname].center.low  = λ[iadd] - λ[iadd]/10. # allow to shift 10%
        model[cname].center.high = λ[iadd] + λ[iadd]/10.

        thaw(model, cname)
        bestfit = fit!(model, source.data[id], minimizer=mzer); show(logio(source), bestfit)
        freeze(model, cname)
    end
    evaluate!(model)
end


function neglect_weak_features!(source::QSO{T}, model::Model; id=1) where T <: DefaultRecipe
    # Disable "unknown" lines whose normalization uncertainty is larger
    # than X times the normalization
    needs_fitting = false
    for ii in 1:source.options[:n_unk]
        cname = Symbol(:unk, ii)
        isfixed(model, cname)  &&  continue
        if bestfit[cname].norm.val == 0.
            freeze(model, cname)
            needs_fitting = true
            println(logio(source), "Disabling $cname (norm. = 0)")
        elseif bestfit[cname].norm.unc / bestfit[cname].norm.val > 3
            model[cname].norm.val = 0.
            freeze(model, cname)
            needs_fitting = true
            println(logio(source), "Disabling $cname (unc. / norm. > 3)")
        end
    end
    return needs_fitting
end


function fit(source::QSO{TRecipe}; id=1) where TRecipe <: q1927p654
    elapsed = time()
    mzer = default_minimizer(TRecipe)
    model = Model(source.domain[id])

    # TODO if source.options[:instr_broadening]
    # TODO     GFit.set_instr_response!(model[1], (l, f) -> instrumental_broadening(l, f, source.spectra[id].resolution))
    # TODO end

    println(logio(source), "\nFit continuum components...")
    add_qso_continuum!(source, model)
    add_host_galaxy!(source, model)
    add_balmer_cont!(source, model)
    bestfit = fit!(model, source.data[id], minimizer=mzer);  show(logio(source), bestfit)
    renorm_cont!(source, model)
    freeze(model, :qso_cont)
    haskey(model, :galaxy)  &&  freeze(model, :galaxy)
    haskey(model, :balmer)  &&  freeze(model, :balmer)
    evaluate!(model)

    println(logio(source), "\nFit iron templates...")
    model[:Iron] = SumReducer([])
    model[:main] = SumReducer([:Continuum, :Iron])
    add_iron_uv!( source, model)
    add_iron_opt!(source, model)
    evaluate!(model)

    if length(model[:Iron].list) > 0
        bestfit = fit!(model, source.data[id], minimizer=mzer); show(logio(source), bestfit)
        haskey(model, :ironuv   )  &&  freeze(model, :ironuv)
        haskey(model, :ironoptbr)  &&  freeze(model, :ironoptbr)
        haskey(model, :ironoptna)  &&  freeze(model, :ironoptna)
    end
    evaluate!(model)

    println(logio(source), "\nFit known emission lines...")
    add_emission_lines!(source, model)
    guess_emission_lines_values!(source, model)
    add_patch_functs!(source, model)


    # Force Hg and Hb lines to have the same shape as Ha
    model[:br_Ha].norm.low = 1.e-10  # avoid division by zero

    @patch! begin
        model[:na_Hg].norm = model[:br_Hg].norm * (model[:na_Ha].norm / model[:br_Ha].norm)
        model[:na_Hg].voff = model[:na_Ha].voff
        model[:na_Hg].fwhm = model[:na_Ha].fwhm

        model[:bb_Hg].norm = model[:br_Hg].norm * (model[:bb_Ha].norm / model[:br_Ha].norm)
        model[:bb_Hg].voff = model[:bb_Ha].voff
        model[:bb_Hg].fwhm = model[:bb_Ha].fwhm

        model[:br_Hg].voff = model[:br_Ha].voff
        model[:br_Hg].fwhm = model[:br_Ha].fwhm
    end

    @patch! begin
        model[:na_Hb].norm = model[:br_Hb].norm * (model[:na_Ha].norm / model[:br_Ha].norm)
        model[:na_Hb].voff = model[:na_Ha].voff
        model[:na_Hb].fwhm = model[:na_Ha].fwhm

        model[:bb_Hb].norm = model[:br_Hb].norm * (model[:bb_Ha].norm / model[:br_Ha].norm)
        model[:bb_Hb].voff = model[:bb_Ha].voff
        model[:bb_Hb].fwhm = model[:bb_Ha].fwhm

        model[:br_Hb].voff = model[:br_Ha].voff
        model[:br_Hb].fwhm = model[:br_Ha].fwhm
    end

    bestfit = fit!(model, source.data[id], minimizer=mzer); show(logio(source), bestfit)
    for lname in line_names
        freeze(model, lname)
    end

    println(logio(source), "\nFit unknown emission lines...")
    add_unknown_lines!(source, model)

    println(logio(source), "\nLast run with all parameters free...")
    thaw(model, :qso_cont)
    (:galaxy in keys(model))     &&  thaw(model, :galaxy)
    (:balmer in keys(model))     &&  thaw(model, :balmer)
    (:ironuv    in keys(model))  &&  thaw(model, :ironuv)
    (:ironoptbr in keys(model))  &&  thaw(model, :ironoptbr)
    (:ironoptna in keys(model))  &&  thaw(model, :ironoptna)
    for lname in keys(source.line_names[id])
        thaw(model, lname)
    end
    for j in 1:source.options[:n_unk]
        cname = Symbol(:unk, j)
        if model[cname].norm.val > 0
            thaw(model, cname)
        else
            freeze(model, cname)
        end
    end
    bestfit = fit!(model, source.data[id], minimizer=mzer)

    if neglect_weak_features!(source, model)
        println(logio(source), "\nRe-run fit...")
        bestfit = fit!(model, source.data[id], minimizer=mzer)
    end

    println(logio(source))
    show(logio(source), bestfit)

    out = QSFit.QSFitResults(source, model, bestfit)
    elapsed = time() - elapsed
    println(logio(source), "\nElapsed time: $elapsed s")
    QSFit.close_logio(source)
    return out
end