Recipe_single.jl

function fit(source::QSO{TRecipe}; id=1) where TRecipe <: q1927p654
    elapsed = time()
    mzer = GFit.cmpfit()
    mzer.config.ftol = mzer.config.gtol = mzer.config.xtol = 1.e-6

    # Initialize components and guess initial values
    println(source.log, "\nFit continuum components...")
    λ = source.domain[id][:]
    model = Model(source.domain[id], :Continuum => Reducer(sum, [:qso_cont]),
                  :qso_cont => QSFit.qso_cont_component(TRecipe))

    if source.options[:instr_broadening]
        GFit.set_instr_response!(model[1], (l, f) -> instrumental_broadening(l, f, source.spectra[id].resolution))
    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
    if source.options[:use_host_template]  &&
        (minimum(λ) .< 5500 .< maximum(λ))
        add!(model, :Continuum => Reducer(sum, [:qso_cont, :galaxy]),
             :galaxy => QSFit.hostgalaxy(source.options[:host_template]))
        model[:galaxy].norm.val = Spline1D(λ, source.data[id].val, k=1, bc="error")(5500.)
    end

    # Balmer continuum and pseudo-continuum
    if source.options[:use_balmer]
        tmp = [:qso_cont, :balmer]
        (:galaxy in keys(model))  &&  push!(tmp, :galaxy)
        add!(model, :Continuum => Reducer(sum, tmp),
             :balmer => QSFit.balmercont(0.1, 0.5))
        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) do m
            m[:balmer].norm *= m[:qso_cont].norm
        end
    end

    bestfit = fit!(model, source.data, minimizer=mzer);  show(source.log, bestfit)

    # QSO continuum renormalization
    freeze(model, :qso_cont)
    c = model[:qso_cont]
    if c.norm.val > 0
        println(source.log, "Cont. norm. (before): ", c.norm.val)
        initialnorm = 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(source.log, "Cont. norm. (after) : ", c.norm.val)
    else
        println(source.log, "Skipping cont. renormalization")
    end
    freeze(model, :qso_cont)
    (:galaxy in keys(model))  &&  freeze(model, :galaxy)
    (:balmer in keys(model))  &&  freeze(model, :balmer)
    evaluate!(model)

    # Fit iron templates
    println(source.log, "\nFit iron templates...")
    iron_components = Vector{Symbol}()
    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
            add!(model, :ironuv => comp)
            model[:ironuv].norm.val = 0.5
            push!(iron_components, :ironuv)
        else
            println(source.log, "Ignoring ironuv component (threshold: $threshold)")
        end
    end

    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))
            add!(model,
                 :ironoptbr => comp,
                 :ironoptna => QSFit.ironopt_narrow(fwhm))
            model[:ironoptbr].norm.val = 0.5
            model[:ironoptna].norm.val = 0.0
            freeze(model, :ironoptna)  # will be freed during last run
            push!(iron_components, :ironoptbr, :ironoptna)
        else
            println(source.log, "Ignoring ironopt component (threshold: $threshold)")
        end
    end
    if length(iron_components) > 0
        add!(model, :Iron => Reducer(sum, iron_components))
        add!(model, :main => Reducer(sum, [:Continuum, :Iron]))
        evaluate!(model)
        bestfit = fit!(model, source.data, minimizer=mzer); show(source.log, bestfit)
    else
        add!(model, :Iron => Reducer(() -> [0.], Symbol[]))
        add!(model, :main => Reducer(sum, [: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]))
    line_groups = unique(collect(values(source.line_names[id])))
    println(source.log, "\nFit known emission lines...")
    resid = source.data[id].val - model()  # will be used to guess line normalization
    add!(model, source.line_comps[id])
    for (group, lnames) in QSFit.invert_dictionary(source.line_names[id])
        add!(model, group  => Reducer(sum, lnames))
    end
    add!(model, :main => Reducer(sum, [: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
    for cname in line_names
        model[cname].norm_integrated = source.options[:norm_integrated]
    end

    # Guess values
    evaluate!(model)
    y = source.data[id].val - model()
    for cname in line_names
        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(source.log, "Line $cname is not a Gaussian profile: Can't take spectral resolution into account")
            end
        end
    end

    # Patch parameters
    if  haskey(model, :OIII_4959)  &&
        haskey(model, :OIII_5007)
        model[:OIII_4959].voff.fixed = true
        patch!(model) do m
            m[:OIII_4959].voff = m[:OIII_5007].voff
        end
    end
    if  haskey(model, :NII_6549)  &&
        haskey(model, :NII_6583)
        model[:NII_6549].voff.fixed = true
        patch!(model) do m
            m[:NII_6549].voff = m[:NII_6583].voff
        end
    end
    if  haskey(model, :OIII_5007_bw)  &&
        haskey(model, :OIII_5007)
        patch!(model) do m
            m[:OIII_5007_bw].voff += m[:OIII_5007].voff
            m[:OIII_5007_bw].fwhm += m[:OIII_5007].fwhm
        end
    end

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

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

    # Avoid division by zero
    model[:br_Ha].norm.low = 1.e-10

    model[:na_Hg].norm.fixed = 1
    model[:na_Hg].fwhm.fixed = 1
    model[:na_Hg].voff.fixed = 1
    model[:bb_Hg].norm.fixed = 1
    model[:bb_Hg].fwhm.fixed = 1
    model[:bb_Hg].voff.fixed = 1
    model[:br_Hg].fwhm.fixed = 1
    model[:br_Hg].voff.fixed = 1
    patch!(model) do m
        m[:na_Hg].norm = m[:br_Hg].norm * (m[:na_Ha].norm / m[:br_Ha].norm)
        m[:na_Hg].voff = m[:na_Ha].voff
        m[:na_Hg].fwhm = m[:na_Ha].fwhm

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

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

    model[:na_Hb].norm.fixed = 1
    model[:na_Hb].fwhm.fixed = 1
    model[:na_Hb].voff.fixed = 1
    model[:bb_Hb].norm.fixed = 1
    model[:bb_Hb].fwhm.fixed = 1
    model[:bb_Hb].voff.fixed = 1
    model[:br_Hb].fwhm.fixed = 1
    model[:br_Hb].voff.fixed = 1
    patch!(model) do m
        m[:na_Hb].norm = m[:br_Hb].norm * (m[:na_Ha].norm / m[:br_Ha].norm)
        m[:na_Hb].voff = m[:na_Ha].voff
        m[:na_Hb].fwhm = m[:na_Ha].fwhm

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

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

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

    # Add unknown lines
    println(source.log, "\nFit unknown emission lines...")
    if source.options[:n_unk] > 0
        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
        add!(model, :UnkLines => Reducer(sum, collect(keys(tmp))), tmp)
        add!(model, :main => Reducer(sum, [: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, minimizer=mzer); show(source.log, bestfit)
            freeze(model, cname)
        end
    end
    evaluate!(model)

    # Last run with all parameters free
    println(source.log, "\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 line_names
        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, minimizer=mzer); show(source.log, bestfit)

    # 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(source.log, "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(source.log, "Disabling $cname (unc. / norm. > 3)")
        end
    end
    if needs_fitting
        println(source.log, "\nRe-run fit...")
        bestfit = fit!(model, source.data, minimizer=mzer); show(source.log, bestfit)
    end

    println(source.log, "\nFinal model and bestfit:")
    show(source.log, model)
    println(source.log)
    show(source.log, bestfit)

    elapsed = time() - elapsed
    println(source.log, "\nElapsed time: $elapsed s")
    QSFit.close_log(source)

    return (model, bestfit)
end