Updated

function fit(source::QSO{TRecipe}; id=1) where TRecipe <: q1927p654

    elapsed = time()

function minimizer(source::QSO{T}) where T <: DefaultRecipe

    mzer = GFit.cmpfit()

    mzer.Δfitstat_theshold = 1.e-5

    return mzer

end

    # Initialize components and guess initial values

    println(logio(source), "\nFit continuum components...")

    λ = source.domain[id][:]

    model = Model(source.domain[id],

                  :qso_cont => QSFit.qso_cont_component(TRecipe),

                  :Continuum => SumReducer([:qso_cont]))

    # 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

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

    c = model[:qso_cont]

    c.x0.val = median(λ)

    c.norm.val = Spline1D(λ, source.data[id].val, k=1, bc="error")(c.x0.val)

    # Host galaxy template

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[: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

    # Balmer continuum and pseudo-continuum

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)

        c.ratio.high = 1

        @patch! model[:balmer].norm *= model[:qso_cont].norm

    end

end

    bestfit = fit!(model, source.data[id], minimizer=mzer);  show(logio(source), bestfit)

    # QSO continuum renormalization

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

    else

        println(logio(source), "Skipping cont. renormalization")

    end

    freeze(model, :qso_cont)

    (:galaxy in keys(model))  &&  freeze(model, :galaxy)

    (:balmer in keys(model))  &&  freeze(model, :balmer)

    evaluate!(model)

end

    # Fit iron templates

    println(logio(source), "\nFit iron templates...")

    iron_components = Vector{Symbol}()

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))

        if coverage >= threshold

            model[:ironuv] = comp

            model[:ironuv].norm.val = 0.5

            push!(iron_components, :ironuv)

            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))

            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!(iron_components, :ironoptbr, :ironoptna)

            push!(model[:Iron].list, :ironoptbr)

            push!(model[:Iron].list, :ironoptna)

        else

            println(logio(source), "Ignoring ironopt component (threshold: $threshold)")

        end

    end

    if length(iron_components) > 0

        model[:Iron] = SumReducer(iron_components)

        model[:main] = SumReducer([:Continuum, :Iron])

        evaluate!(model)

        bestfit = fit!(model, source.data[id], minimizer=mzer); show(logio(source), bestfit)

    else

        model[:Iron] = @expr m -> [0.]

        model[:main] = SumReducer([:Continuum, :Iron])

end

    (:ironuv    in keys(model))  &&  freeze(model, :ironuv)

    (:ironoptbr in keys(model))  &&  freeze(model, :ironoptbr)

    (:ironoptna in keys(model))  &&  freeze(model, :ironoptna)

    evaluate!(model)

    # Add emission lines

    line_names = collect(keys(source.line_names[id]))

function add_emission_lines!(source::QSO{T}, model::Model; id=1) where T <: DefaultRecipe

    line_groups = unique(collect(values(source.line_names[id])))

    println(logio(source), "\nFit known emission lines...")

    resid = source.data[id].val - model()  # will be used to guess line normalization

    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[:OIII_5007_bw].voff.low  = 0

        model[:OIII_5007_bw].voff.high = 2e3

    end

    for cname in line_names

        model[cname].norm_integrated = source.options[:norm_integrated]

    evaluate!(model)

end

    # Guess values

    evaluate!(model)

    for cname in line_names

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))

            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[:bb_Ha].norm.val  = 0.5

        @patch! model[:bb_Ha].norm *= model[:br_Ha].norm / model[:br_Ha].fwhm * model[:bb_Ha].fwhm

    end

    # 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

    # Add unknown lines

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

    if source.options[:n_unk] > 0

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))

        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)

    # ----------------------------------------------------------------

    # Last run with all parameters free

    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)

    (:ironuv    in keys(model))  &&  thaw(model, :ironuv)

    (:ironoptbr in keys(model))  &&  thaw(model, :ironoptbr)

    (:ironoptna in keys(model))  &&  thaw(model, :ironoptna)

    for lname in line_names

    for lname in keys(source.line_names[id])

        thaw(model, lname)

    end

    for j in 1:source.options[:n_unk]

    end

    bestfit = fit!(model, source.data[id], minimizer=mzer)

    # Disable "unknown" lines whose normalization uncertainty is larger

    # than 3 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

    if needs_fitting

    if neglect_weak_features!(source, model)

        println(logio(source), "\nRe-run fit...")

        bestfit = fit!(model, source.data[id], minimizer=mzer)

    end