text:
"""
`module  JAC.IsotopeShift`
    ... a submodel of JAC that contains all methods for computing isotope-dependent properties for some level(s).
"""
module IsotopeShift

    using Printf, ..Basics, ..Defaults, ..InteractionStrength, ..ManyElectron, ..Nuclear, ..Radial, ..RadialIntegrals,
                  ..SpinAngular, ..TableStrings

    """
    `struct  IsotopeShift.Settings`
        ... defines a type for the details and parameters of computing isotope-shift M and F parameters.

        + calcNMS                  ::Bool             ... True if mass-shift parameters M_nmn need to be calculated, and false otherwise.
        + calcSMS                  ::Bool             ... True if mass-shift parameters M_sms need to be calculated, and false otherwise.
        + calcF                    ::Bool             ... True if the field-shift parameter need to be calculated, and false otherwise.
        + calcBoson                ::Bool             ... True if the boson-field shift parameter need to be calculated, and false otherwise.
        + printBefore              ::Bool             ... True if a list of selected levels is printed before the actual computations start.
        + bosonMass                ::Float64          ... mass of the scalar boson [e_electron].
        + levelSelection           ::LevelSelection   ... Specifies the selected levels, if any.
    """
    struct Settings
        calcNMS                    ::Bool
        calcSMS                    ::Bool
        calcF                      ::Bool
        calcBoson                  ::Bool
        printBefore                ::Bool
        bosonMass                  ::Float64
        levelSelection             ::LevelSelection
    end

    """
    `IsotopeShift.Settings(; calcNMS::Bool=true,` calcSMS::Bool=false, calcF::Bool=false, calcBoson::Bool=false,
                             printBefore::Bool=true, bosonMass::Float64=0., levelSelection::LevelSelection=LevelSelection())
        ... keyword constructor to overwrite selected value of isoshift computations.
    """
    function Settings(; calcNMS::Bool=true, calcSMS::Bool=false, calcF::Bool=false, calcBoson::Bool=false,
                             printBefore::Bool=true, bosonMass::Float64=0., levelSelection::LevelSelection=LevelSelection())
        Settings(calcNMS, calcSMS, calcF, calcBoson, printBefore, bosonMass, levelSelection)
    end

    # `Base.show(io::IO, settings::IsotopeShift.Settings)`  ... prepares a proper printout of the variable settings::IsotopeShift.Settings.
    function Base.show(io::IO, settings::IsotopeShift.Settings)
        println(io, "calcNMS:                  $(settings.calcNMS)  ")
        println(io, "calcSMS:                  $(settings.calcSMS)  ")
        println(io, "calcF:                    $(settings.calcF)  ")
        println(io, "calcBoson:                $(settings.calcBoson)  ")
        println(io, "printBefore:              $(settings.printBefore)  ")
        println(io, "bosonMass:                $(settings.bosonMass)  ")
        println(io, "levelSelection:           $(settings.levelSelection)  ")
    end

    """
    `struct  IsotopeShift.Outcome`
        ... defines a type to keep the outcome of a isotope-shift computation, such as the K and F parameters as well other
            results.

        + level                     ::Level              ... Atomic level to which the outcome refers to.
        + Knms                      ::Float64            ... K_nms parameter
        + Ksms                      ::Float64            ... K_sms parameter
        + Fme                       ::Float64            ... F parameter from matrix element.
        + Fdensity                  ::Float64            ... F parameter from density.
        + Xboson                    ::Float64            ... X boson-field shift constant.
        + amplitudeKnms             ::Complex{Float64}   ... K_nms amplitude
        + amplitudeKsmsA            ::Complex{Float64}   ... K_sms,A amplitude
        + amplitudeKsmsB            ::Complex{Float64}   ... K_sms,B amplitude
        + amplitudeKsmsC            ::Complex{Float64}   ... K_sms,C amplitude
    """
    struct Outcome
        level                       ::Level
        Knms                        ::Float64
        Ksms                        ::Float64
        Fme                         ::Float64
        Fdensity                    ::Float64
        Xboson                      ::Float64
        amplitudeKnms               ::Complex{Float64}
        amplitudeKsmsA              ::Complex{Float64}
        amplitudeKsmsB              ::Complex{Float64}
        amplitudeKsmsC              ::Complex{Float64}
    end

    """
    `IsotopeShift.Outcome()`  ... constructor for an `empty` instance of Hfs.Outcome for the computation of isotope-shift properties.
    """
    function Outcome()
        Outcome(Level(), 0., 0., 0., 0., 0.,   0., 0., 0., 0.)
    end

    # `Base.show(io::IO, outcome::IsotopeShift.Outcome)`  ... prepares a proper printout of the variable outcome::IsotopeShift.Outcome.
    function Base.show(io::IO, outcome::IsotopeShift.Outcome)
        println(io, "level:                   $(outcome.level)  ")
        println(io, "Knms:                    $(outcome.Knms)  ")
        println(io, "Ksms:                    $(outcome.Ksms)  ")
        println(io, "Fme:                     $(outcome.Fme)  ")
        println(io, "Fdensity:                $(outcome.Fdensity)  ")
        println(io, "Xboson:                  $(outcome.Xboson)  ")
        println(io, "amplitudeKnms:           $(outcome.amplitudeKnms)  ")
        println(io, "amplitudeKsmsA:          $(outcome.amplitudeKsmsA)  ")
        println(io, "amplitudeKsmsB:          $(outcome.amplitudeKsmsB)  ")
        println(io, "amplitudeKsmsC:          $(outcome.amplitudeKsmsC)  ")
    end

    """
    `IsotopeShift.amplitude(kind::String, rLevel::Level, sLevel::Level, nm::Nuclear.Model, grid::Radial.Grid)`
        ... to compute either the  H^(NMS),  H^(SMS,A),  H^(SMS,B) or  H^(SMS,C) normal and specific mass-shift amplitudes
            <alpha_r J_r || H^(A)) || alpha_s J_s>  for a given pair of levels. A value::ComplexF64 is returned.
    """
    function  amplitude(kind::String, rLevel::Level, sLevel::Level, nm::Nuclear.Model, grid::Radial.Grid)
        nr = length(rLevel.basis.csfs);    ns = length(sLevel.basis.csfs);    matrix = zeros(ComplexF64, nr, ns)
        printstyled("Compute mass-shift $(kind[1:9]) matrix of dimension $nr x $ns ... ", color=:light_green)
        #
        if  rLevel.parity != sLevel.parity   return( ComplexF64(0.) )   end
        #
        for  r = 1:nr
            for  s = 1:ns
                me = 0.
                if  rLevel.basis.csfs[r].parity  != rLevel.parity    ||  sLevel.basis.csfs[s].parity  != sLevel.parity  ||
                    rLevel.parity != sLevel.parity    continue
                end
                #
                if      kind == "H^(NMS)   amplitude"
                #------------------------------------
                    ##x wa = compute("angular coefficients: e-e, Ratip2013", rLevel.basis.csfs[r], sLevel.basis.csfs[s])
                    subshellList = rLevel.basis.subshells
                    op = SpinAngular.OneParticleOperator(0, plus, true)
                    wa = SpinAngular.computeCoefficients(op, rLevel.basis.csfs[r], sLevel.basis.csfs[s], subshellList)
                    for  coeff in wa
                        ja = Basics.subshell_2j(rLevel.basis.orbitals[coeff.a].subshell)
                        jb = Basics.subshell_2j(sLevel.basis.orbitals[coeff.b].subshell)
                        tamp  = InteractionStrength.hamiltonian_nms(rLevel.basis.orbitals[coeff.a], sLevel.basis.orbitals[coeff.b], nm, grid)
                        me = me + coeff.T  * sqrt( ja + 1) * tamp
                    end
                #
                elseif  kind == "H^(SMS,A) amplitude"
                #------------------------------------
                    ##x wa = compute("angular coefficients: e-e, Ratip2013", rLevel.basis.csfs[r], sLevel.basis.csfs[s])
                    subshellList = rLevel.basis.subshells
                    op = SpinAngular.TwoParticleOperator(0, plus, true)
                    wa = SpinAngular.computeCoefficients(op, rLevel.basis.csfs[r], sLevel.basis.csfs[s], subshellList)
                    for  coeff in wa
                        if  coeff.nu != 1  continue   end
                        ja = Basics.subshell_2j(rLevel.basis.orbitals[coeff.a].subshell)
                        jb = Basics.subshell_2j(sLevel.basis.orbitals[coeff.b].subshell)
                        tamp  = InteractionStrength.X_smsA(rLevel.basis.orbitals[coeff.a], rLevel.basis.orbitals[coeff.b],
                                                           sLevel.basis.orbitals[coeff.c], sLevel.basis.orbitals[coeff.d], nm, grid)
                        me = me + coeff.V * sqrt( ja + 1) * tamp
                    end
                #
                elseif  kind == "H^(SMS,B) amplitude"
                #------------------------------------
                    ##x wa = compute("angular coefficients: e-e, Ratip2013", rLevel.basis.csfs[r], sLevel.basis.csfs[s])
                    subshellList = rLevel.basis.subshells
                    op = SpinAngular.TwoParticleOperator(0, plus, true)
                    wa = SpinAngular.computeCoefficients(op, rLevel.basis.csfs[r], sLevel.basis.csfs[s], subshellList)
                    for  coeff in wa
                        if  coeff.nu != 1  continue   end
                        ja = Basics.subshell_2j(rLevel.basis.orbitals[coeff.a].subshell)
                        jb = Basics.subshell_2j(sLevel.basis.orbitals[coeff.b].subshell)
                        tamp  = InteractionStrength.X_smsB(rLevel.basis.orbitals[coeff.a], rLevel.basis.orbitals[coeff.b],
                                                           sLevel.basis.orbitals[coeff.c], sLevel.basis.orbitals[coeff.d], nm, grid)
                        me = me + coeff.V * sqrt( ja + 1) * tamp
                    end
                #
                elseif  kind == "H^(SMS,C) amplitude"
                #------------------------------------
                    ##x wa = compute("angular coefficients: e-e, Ratip2013", rLevel.basis.csfs[r], sLevel.basis.csfs[s])
                    subshellList = rLevel.basis.subshells
                    op = SpinAngular.TwoParticleOperator(0, plus, true)
                    wa = SpinAngular.computeCoefficients(op, rLevel.basis.csfs[r], sLevel.basis.csfs[s], subshellList)
                    for  coeff in wa
                        if  coeff.nu != 1  continue   end
                        ja = Basics.subshell_2j(rLevel.basis.orbitals[coeff.a].subshell)
                        jb = Basics.subshell_2j(sLevel.basis.orbitals[coeff.b].subshell)
                        tamp  = InteractionStrength.X_smsC(rLevel.basis.orbitals[coeff.a], rLevel.basis.orbitals[coeff.b],
                                                           sLevel.basis.orbitals[coeff.c], sLevel.basis.orbitals[coeff.d], nm, grid)
                        me = me + coeff.V * sqrt( ja + 1) * tamp
                    end
                #
                else    error("stop a")
                end
                #
                matrix[r,s] = me
            end
        end
        printstyled("done.\n", color=:light_green)
        amplitude = transpose(rLevel.mc) * matrix * sLevel.mc
        #
        return( amplitude )
    end

    """
    `IsotopeShift.amplitude(kind::String, rLevel::Level, sLevel::Level, potential::Array{Float64,1}, grid::Radial.Grid)`
        ... to compute the H^(field-shift) field-shift amplitude  <alpha_r J_r || H^(field-shift) || alpha_s J_s> or the
            H^(boson-field) shift amplitude  <alpha_r J_r || H^(boson-field) || alpha_s J_s>  for a given pair of levels.
            The potential has to provide delta V^(nuc) for the field-shift amplitudes and the effective potential for the
            boson-field shift. A value::ComplexF64 is returned.
    """
    function  amplitude(kind::String, rLevel::Level, sLevel::Level, potential::Array{Float64,1}, grid::Radial.Grid)
        nr = length(rLevel.basis.csfs);    ns = length(sLevel.basis.csfs);    matrix = zeros(ComplexF64, nr, ns)
        printstyled("Compute field-shift $kind matrix of dimension $nr x $ns ... ", color=:light_green)
        #
        if  rLevel.parity != sLevel.parity   return( ComplexF64(0.) )   end
        #
        for  r = 1:nr
            for  s = 1:ns
                me = 0.
                if  rLevel.basis.csfs[r].parity  != rLevel.parity    ||  sLevel.basis.csfs[s].parity  != sLevel.parity  ||
                    rLevel.parity != sLevel.parity    continue
                end
                #
                if      kind == "H^(field-shift) amplitude"
                #------------------------------------------
                    subshellList = rLevel.basis.subshells
                    op = SpinAngular.OneParticleOperator(0, plus, true)
                    wa = SpinAngular.computeCoefficients(op, rLevel.basis.csfs[r], sLevel.basis.csfs[s], subshellList)
                    for  coeff in wa
                        ja = Basics.subshell_2j(rLevel.basis.orbitals[coeff.a].subshell)
                        jb = Basics.subshell_2j(sLevel.basis.orbitals[coeff.b].subshell)
                        tamp  = InteractionStrength.fieldShift(rLevel.basis.orbitals[coeff.a], sLevel.basis.orbitals[coeff.b], potential, grid)
                        me = me + coeff.T  * sqrt( ja + 1) * tamp
                    end
                #
                elseif  kind == "H^(boson-field) amplitude"
                #------------------------------------------
                    subshellList = rLevel.basis.subshells
                    op = SpinAngular.OneParticleOperator(0, plus, true)
                    wa = SpinAngular.computeCoefficients(op, rLevel.basis.csfs[r], sLevel.basis.csfs[s], subshellList)
                    for  coeff in wa
                        ja = Basics.subshell_2j(rLevel.basis.orbitals[coeff.a].subshell)
                        jb = Basics.subshell_2j(sLevel.basis.orbitals[coeff.b].subshell)
                        tamp  = InteractionStrength.bosonShift(rLevel.basis.orbitals[coeff.a], sLevel.basis.orbitals[coeff.b], potential, grid)
                        me = me + coeff.T  * sqrt( ja + 1) * tamp
                    end
                #
                else    error("stop b")
                end
                #
                matrix[r,s] = me
            end
        end
        printstyled("done.\n", color=:light_green)
        amplitude = transpose(rLevel.mc) * matrix * sLevel.mc
        #
        return( amplitude )
    end

    """
    `IsotopeShift.computeAmplitudesProperties(outcome::IsotopeShift.Outcome, nm::Nuclear.Model, grid::Radial.Grid,
                                              settings::IsotopeShift.Settings)
        ... to compute all amplitudes and properties for a given level; an outcome::IsotopeShift.Outcome is returned for
            which the amplitudes and properties are now evaluated explicitly.
    """
    function  computeAmplitudesProperties(outcome::IsotopeShift.Outcome, nm::Nuclear.Model, grid::Radial.Grid, settings::IsotopeShift.Settings)
        amplitudeKnms   = 0.0;   Ksms  = 0.0;    F  = 0.0;    Xboson = 0.0;
        amplitudeKsmsA  = 0.0;   amplitudeKsmsB  = 0.0;  amplitudeKsmsC  = 0.0;  amplitudeF  = 0.0;
        if  settings.calcNMS
            amplitudeKnms  = IsotopeShift.amplitude("H^(NMS)   amplitude", outcome.level, outcome.level, nm, grid)
        end

        if  settings.calcSMS
            amplitudeKsmsA = IsotopeShift.amplitude("H^(SMS,A) amplitude", outcome.level, outcome.level, nm, grid)
            amplitudeKsmsB = IsotopeShift.amplitude("H^(SMS,B) amplitude", outcome.level, outcome.level, nm, grid)
            amplitudeKsmsC = IsotopeShift.amplitude("H^(SMS,C) amplitude", outcome.level, outcome.level, nm, grid)
            Ksms = (amplitudeKsmsA + amplitudeKsmsB + amplitudeKsmsC)
        end

        if  settings.calcF
            nmp       = Nuclear.Model(nm.Z, nm.mass+1.0)
            deltaPot  = Nuclear.nuclearPotential(nm, grid).Zr - Nuclear.nuclearPotential(nmp, grid).Zr
            deltaPot  = deltaPot / (nm.radius^2 - nmp.radius^2)
            lastPoint = Basics.lastPoint(deltaPot, 1.0e-9);   @show lastPoint
            #
            Fme       = IsotopeShift.amplitude("H^(field-shift) amplitude", outcome.level, outcome.level, deltaPot, grid)
            #
            density   = Basics.computeDensity(outcome.level, grid)
            wb        = zeros( lastPoint )
            wb[1:lastPoint] = - density[1:lastPoint] .* deltaPot[1:lastPoint] ./ grid.r[1:lastPoint]
            Fdensity  = RadialIntegrals.V0(wb, lastPoint, grid)
        end

        if  settings.calcBoson
            potential = zeros( grid.NoPoints );   malpha = settings.bosonMass / Defaults.getDefaults("alpha")
            for     i = 1:length(potential)    potential[i] = exp( -malpha*grid.r[i] ) / grid.r[i]    end
            Xboson    = IsotopeShift.amplitude("H^(boson-field) amplitude", outcome.level, outcome.level, potential, grid)
        end

        newOutcome = IsotopeShift.Outcome( outcome.level, amplitudeKnms, Ksms, Fme, Fdensity, Xboson,
                                           amplitudeKnms, amplitudeKsmsA, amplitudeKsmsB, amplitudeKsmsC )
        return( newOutcome )
    end

    """
    `IsotopeShift.computeOutcomes(multiplet::Multiplet, nm::Nuclear.Model, grid::Radial.Grid, settings::IsotopeShift.Settings; output=true)`
        ... to compute (as selected) the isotope-shift M and F parameters for the levels of the given multiplet and as
            specified by the given settings. The results are printed in neat tables to screen but nothing is returned otherwise.
    """
    function computeOutcomes(multiplet::Multiplet, nm::Nuclear.Model, grid::Radial.Grid, settings::IsotopeShift.Settings; output=true)
        println("")
        printstyled("IsotopeShift.computeOutcomes(): The computation of the isotope-shift parameters starts now ... \n", color=:light_green)
        printstyled("-------------------------------------------------------------------------------------------------- \n", color=:light_green)
        #
        ##x @show settings
        outcomes = IsotopeShift.determineOutcomes(multiplet, settings)
        # Display all selected levels before the computations start
        if  settings.printBefore    IsotopeShift.displayOutcomes(outcomes)    end
        # Calculate all amplitudes and requested properties
        newOutcomes = IsotopeShift.Outcome[]
        for  outcome in outcomes
            newOutcome = IsotopeShift.computeAmplitudesProperties(outcome, nm, grid, settings)
            push!( newOutcomes, newOutcome)
        end
        # Print all results to screen
        IsotopeShift.displayResults(stdout, newOutcomes, settings)
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        if  printSummary    IsotopeShift.displayResults(iostream, newOutcomes, settings)   end
        #
        if    output    return( newOutcomes )
        else            return( nothing )
        end
    end

    """
    `IsotopeShift.determineOutcomes(multiplet::Multiplet, settings::IsotopeShift.Settings)`
        ... to determine a list of Outcomes's for the computation of the isotope-shift M and F parameters for the given multiplet.
            It takes into account the particular selections and settings. An Array{IsotopeShift.Outcome,1} is returned. Apart from
            the level specification, all physical properties are set to zero during the initialization process.
    """
    function  determineOutcomes(multiplet::Multiplet, settings::IsotopeShift.Settings)
        outcomes = IsotopeShift.Outcome[]
        @show settings.levelSelection
        for  level  in  multiplet.levels
            if  Basics.selectLevel(level, settings.levelSelection)
                push!( outcomes, IsotopeShift.Outcome(level, 0., 0., 0., 0., 0.,   0., 0., 0., 0.) )
            end
        end
        return( outcomes )
    end

    """
    `IsotopeShift.displayOutcomes(outcomes::Array{IsotopeShift.Outcome,1})`
        ... to display a list of levels that have been selected for the computations. A small neat table of all
            selected levels and their energies is printed but nothing is returned otherwise.
    """
    function  displayOutcomes(outcomes::Array{IsotopeShift.Outcome,1})
        nx = 43
        println(" ")
        println("  Selected IsotopeShift levels:")
        println(" ")
        println("  ", TableStrings.hLine(nx))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(10, "Level"; na=2);                             sb = sb * TableStrings.hBlank(12)
        sa = sa * TableStrings.center(10, "J^P";   na=4);                             sb = sb * TableStrings.hBlank(14)
        sa = sa * TableStrings.center(14, "Energy"; na=4);
        sb = sb * TableStrings.center(14, TableStrings.inUnits("energy"); na=4)
        println(sa);    println(sb);    println("  ", TableStrings.hLine(nx))
        #
        for  outcome in outcomes
            sa  = "  ";    sym = Basics.LevelSymmetry( outcome.level.J, outcome.level.parity)
            sa = sa * TableStrings.center(10, TableStrings.level(outcome.level.index); na=2)
            sa = sa * TableStrings.center(10, string(sym); na=4)
            sa = sa * @sprintf("%.8e", Defaults.convertUnits("energy: from atomic", outcome.level.energy)) * "    "
            println( sa )
        end
        println("  ", TableStrings.hLine(nx))
        #
        return( nothing )
    end

    """
    `IsotopeShift.displayResults(stream::IO, outcomes::Array{Hfs.Outcome,1}, settings::IsotopeShift.Settings)`
        ... to display the energies, M_ms and F-parameters, etc. for the selected levels. A neat table is printed but nothing
            is returned otherwise.
    """
    function  displayResults(stream::IO, outcomes::Array{IsotopeShift.Outcome,1}, settings::IsotopeShift.Settings)
        nx = 133
        println(stream, " ")
        println(stream, "  IsotopeShift parameters and amplitudes:  Just take the difference K_i - K_f for transition lines.")
        println(stream, " ")
        println(stream, "  Boson mass = $(settings.bosonMass)")
        println(stream, " ")
        println(stream, "  ", TableStrings.hLine(nx))
        sa = "  ";   sb = "  ";   sc = "  "
        sa = sa * TableStrings.center(10, "Level"; na=2);                             sb = sb * TableStrings.hBlank(12)
        sa = sa * TableStrings.center(10, "J^P";   na=4);                             sb = sb * TableStrings.hBlank(14)
        sa = sa * TableStrings.center(14, "Energy"; na=5)
        sb = sb * TableStrings.center(14, TableStrings.inUnits("energy"); na=5);  sc = sc * TableStrings.hBlank(45)
        sa = sa * TableStrings.center(11, "K_nms"; na=4)
        sb = sb * TableStrings.center(11, "[GHz u]" ; na=4)
        sc = sc * TableStrings.center(11, "[a.u.]" ; na=4)
        sa = sa * TableStrings.center(11, "K_sms"; na=4)
        sb = sb * TableStrings.center(11, "[GHz u]" ; na=4)
        sc = sc * TableStrings.center(11, "[a.u.]" ; na=4)
        sa = sa * TableStrings.center(11, "K_ms "; na=4)
        sb = sb * TableStrings.center(11, "[GHz u]" ; na=4)
        sc = sc * TableStrings.center(11, "[a.u.]" ; na=4)
        sa = sa * TableStrings.center(11, " F [ME]";   na=4)
        sb = sb * TableStrings.center(11, "[MHz/fm^2]" ; na=4)
        sa = sa * TableStrings.center(11, " F [dens]";   na=5)
        sb = sb * TableStrings.center(11, "[MHz/fm^2]" ; na=5)
        sa = sa * TableStrings.center(11, "X^(boson)";   na=4)
        sb = sb * TableStrings.center(11, "a.u."; na=4)
        println(stream, sa);    println(stream, sb);    ## println(stream, sc);
        println(stream, "  ", TableStrings.hLine(nx))
        #
        for  outcome in outcomes
            sa  = "  ";    sym = LevelSymmetry( outcome.level.J, outcome.level.parity)
            sa = sa * TableStrings.center(10, TableStrings.level(outcome.level.index); na=2)
            sa = sa * TableStrings.center(10, string(sym); na=4)
            sa = sa * @sprintf("%.8e", Defaults.convertUnits("energy: from atomic", outcome.level.energy))            * "    "
            sa = sa * @sprintf("%.5e", outcome.Knms * 3609.4824)                                                      * "    "
            sa = sa * @sprintf("%.5e", outcome.Ksms * 3609.4824)                                                      * "    "
            sa = sa * @sprintf("%.5e", (outcome.Knms + outcome.Ksms) * 3609.4824)                                     * "    "
            sa = sa * @sprintf("%.5e", Defaults.convertUnits("energy: from atomic to Hz", outcome.Fme)      * 1.0e-6) * "    "
            sa = sa * @sprintf("%.5e", Defaults.convertUnits("energy: from atomic to Hz", outcome.Fdensity) * 1.0e-6) * "    "
            sa = sa * @sprintf("%.5e", Defaults.convertUnits("energy: from atomic to Hz", outcome.Xboson) )
            println(stream, sa )
            ## sb = TableStrings.hBlank(47)
        end
        println(stream, "  ", TableStrings.hLine(nx), "\n\n")
        #
        #
        # Now printout the individual amplitudes in an extra table
        nx = 158
        println(stream, "  ", TableStrings.hLine(nx))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(10, "Level"; na=2);                             sb = sb * TableStrings.hBlank(12)
        sa = sa * TableStrings.center(10, "J^P";   na=4);                             sb = sb * TableStrings.hBlank(14)
        sa = sa * TableStrings.center(14, "Energy"; na=4);
        sb = sb * TableStrings.center(14, TableStrings.inUnits("energy"); na=4)
        sa = sa * TableStrings.center(26, "H^(nms) amplitude"    ; na=3);             sb = sb * TableStrings.hBlank(35)
        sa = sa * TableStrings.center(26, "H^(sms,A) amplitude"  ; na=3);             sb = sb * TableStrings.hBlank(35)
        sa = sa * TableStrings.center(26, "H^(sms,B) amplitude"  ; na=3);             sb = sb * TableStrings.hBlank(35)
        sa = sa * TableStrings.center(26, "H^(sms,C) amplitude"  ; na=3);             sb = sb * TableStrings.hBlank(35)
        println(stream, sa);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(nx))
        #
        for  outcome in outcomes
            sa  = "  ";    sym = LevelSymmetry( outcome.level.J, outcome.level.parity)
            sa = sa * TableStrings.center(10, TableStrings.level(outcome.level.index); na=2)
            sa = sa * TableStrings.center(10, string(sym); na=4)
            sa = sa * @sprintf("%.8e", Defaults.convertUnits("energy: from atomic", outcome.level.energy))          * "    "
            sa = sa * @sprintf("%.5e %s %.5e", outcome.amplitudeKnms.re,  " ", outcome.amplitudeKnms.im)  * "    "
            sa = sa * @sprintf("%.5e %s %.5e", outcome.amplitudeKsmsA.re, " ", outcome.amplitudeKsmsA.im) * "    "
            sa = sa * @sprintf("%.5e %s %.5e", outcome.amplitudeKsmsB.re, " ", outcome.amplitudeKsmsB.im) * "    "
            sa = sa * @sprintf("%.5e %s %.5e", outcome.amplitudeKsmsC.re, " ", outcome.amplitudeKsmsC.im) * "    "
            println(stream, sa )
        end
        println(stream, "  ", TableStrings.hLine(nx))
        #
        #
        # Now printout the transition isotope parameters
        nx = 144
        println(stream, " ")
        println(stream, "  IsotopeShift parameters for individual transitions:")
        println(stream, " ")
        println(stream, "  ", TableStrings.hLine(nx))
        sa = "  ";   sb = "  ";   sc = "  "
        sa = sa * TableStrings.center(14, "Transition";    na=2);                        sb = sb * TableStrings.hBlank(16)
        sa = sa * TableStrings.center(14, "I -- J^P -- F"; na=4);                        sb = sb * TableStrings.hBlank(18)
        sa = sa * TableStrings.center(14, "Energy"; na=5)
        sb = sb * TableStrings.center(14, TableStrings.inUnits("energy"); na=5);  sc = sc * TableStrings.hBlank(45)
        sa = sa * TableStrings.center(11, "K_nms"; na=5)
        sb = sb * TableStrings.center(11, "[GHz u]" ; na=5)
        sc = sc * TableStrings.center(11, "[a.u.]" ; na=5)
        sa = sa * TableStrings.center(11, "K_sms"; na=5)
        sb = sb * TableStrings.center(11, "[GHz u]" ; na=5)
        sc = sc * TableStrings.center(11, "[a.u.]" ; na=5)
        sa = sa * TableStrings.center(11, "K_ms "; na=5)
        sb = sb * TableStrings.center(11, "[GHz u]" ; na=5)
        sc = sc * TableStrings.center(11, "[a.u.]" ; na=5)
        sa = sa * TableStrings.center(11, " F [ME]";   na=4)
        sb = sb * TableStrings.center(11, "[MHz/fm^2]" ; na=4)
        sa = sa * TableStrings.center(11, " F [dens]";   na=4)
        sb = sb * TableStrings.center(11, "[MHz/fm^2]" ; na=4)
        sa = sa * TableStrings.center(11, "X^(boson)";   na=4)
        sb = sb * TableStrings.center(11, "a.u."; na=4)
        println(stream, sa);    println(stream, sb);    ## println(stream, sc);
        println(stream, "  ", TableStrings.hLine(nx))
        #
        for  ioutcome in outcomes
            for  foutcome in outcomes
            if ioutcome.level.index >= foutcome.level.index     continue    end
            sa  = "";    isym = LevelSymmetry( ioutcome.level.J, ioutcome.level.parity)
                         fsym = LevelSymmetry( foutcome.level.J, foutcome.level.parity)
            sa = sa * TableStrings.center(14, TableStrings.levels_if(ioutcome.level.index, foutcome.level.index); na=2)
            sa = sa * TableStrings.center(18, string(isym) * "  --   " * string(fsym); na=2)
            en = abs(ioutcome.level.energy - foutcome.level.energy)
            sa = sa * @sprintf("%.8e", Defaults.convertUnits("energy: from atomic", en))                              * "     "
            sa = sa * @sprintf("%.5e", (ioutcome.Knms-foutcome.Knms) * 3609.4824)                                     * "    "
            sa = sa * @sprintf("%.5e", (ioutcome.Ksms-foutcome.Ksms) * 3609.4824)                                     * "     "
            sa = sa * @sprintf("%.5e", (ioutcome.Knms + ioutcome.Ksms - foutcome.Knms - foutcome.Ksms) * 3609.4824)   * "     "
            sa = sa * @sprintf("%.5e", Defaults.convertUnits("energy: from atomic to Hz", ioutcome.Fme-foutcome.Fme)  * 1.0e-6) * "    "
            sa = sa * @sprintf("%.5e", Defaults.convertUnits("energy: from atomic to Hz", ioutcome.Fdensity-foutcome.Fdensity)  * 1.0e-6) * "    "
            sa = sa * @sprintf("%.5e", Defaults.convertUnits("energy: from atomic to Hz", ioutcome.Xboson-foutcome.Xboson))
            println(stream, sa )
            ## sb = TableStrings.hBlank(47)
            end
        end
        println(stream, "  ", TableStrings.hLine(nx), "\n\n")

        return( nothing )
    end

end # module


text:
using Atum
using Atum.Euler

using PGFPlotsX
using StaticArrays: SVector

function sod(law, x⃗)
  FT = eltype(law)
  ρ = x⃗[1] < 1 // 2 ? 1 : 1 // 8
  ρu⃗ = SVector(FT(0))
  p = x⃗[1] < 1 // 2 ? 1 : 1 // 10
  ρe = Euler.energy(law, ρ, ρu⃗, p)
  SVector(ρ, ρu⃗..., ρe)
end

import Atum: boundarystate
boundarystate(law::EulerLaw, n⃗, x⃗, q⁻, _) = sod(law, x⃗)

function run(A, FT, N, K)
  Nq = N + 1

  law = EulerLaw{FT, 1}()

  cell = LobattoCell{FT, A}(Nq)
  v1d = range(FT(0), stop=FT(1), length=K+1)
  grid = brickgrid(cell, (v1d,); periodic=(false,))

  dg = ESDGSEM(; law, cell, grid,
               volume_numericalflux = EntropyConservativeFlux(),
               surface_numericalflux = RusanovFlux())

  cfl = FT(1 // 4)
  dt = cfl * step(v1d) / N / Euler.soundspeed(law, FT(1), FT(1))
  timeend = FT(2 // 10)

  q = sod.(Ref(law), points(grid))

  odesolver = LSRK54(dg, q, dt)
  solve!(q, timeend, odesolver)

  @pgf begin
    ρ, ρu, ρe = components(q)
    p = Euler.pressure.(Ref(law), ρ, ρu, ρe)
    u = ρu ./ ρ
    x = vec(first(components(points(grid))))

    fig = @pgf GroupPlot({group_style= {group_size="2 by 2"}})
    ρ_plot = Plot({no_marks}, Coordinates(x, vec(ρ)))
    u_plot = Plot({no_marks}, Coordinates(x, vec(u)))
    E_plot = Plot({no_marks}, Coordinates(x, vec(ρe)))
    p_plot = Plot({no_marks}, Coordinates(x, vec(p)))

    push!(fig, {}, ρ_plot)
    push!(fig, {}, u_plot)
    push!(fig, {}, E_plot)
    push!(fig, {}, p_plot)

    path = mkpath(joinpath("output", "euler", "sod"))
    pgfsave(joinpath(path, "sod.pdf"), fig)
  end
end

let
  A = Array
  FT = Float64
  N = 4
  K = 32

  errf = run(A, FT, N, K)
end

ext: jl


text:
# Read in a set of trees from a file
# Returns an array of trees
function readtree(filepath::ASCIIString, format="nwk")
	instream = open(expanduser(filepath))
	instring = readall(instream)
	close(instream)
	trees = split(instring, ';')
	trees = [replace(i, r"(\r|\n|\s)", "") for i in trees]
	trees = trees[bool([length(t) > 0 for t in trees])]
	output_trees = Array(Phylogeny, length(trees))
	for i in 1:length(trees)
		if search(trees[i], ":") == 0:-1
			output_trees[i] = cladobuild(trees[i])
		elseif search(trees[i], ":") != 0:-1
			output_trees[i] = treebuild(trees[i])
		end
	end
	return output_trees
end

# Sub function for creation of a Clado structs from newick format strings.
# Used to build Clado structs from newick strings during operation of the
# readtree function.
function cladobuild(tp::ASCIIString)
	function AddInternal(edge, currentNode, node, index, j)
		edge[j, 1] = currentNode
		node += 1
		edge[j, 2] = currentNode = node
		index[node] = j
		j += 1
		return currentNode, node, j
	end
	function AddTerminal(edge, currentNode, tip, index, tipLabel, tpc, k, j)
		edge[j, 1] = currentNode
		edge[j, 2] = tip
		index[tip] = j
		tipLabel[tip] = tpc[k]
		k += 1
		tip += 1
		j += 1
		return currentNode, tip, k, j
	end
	function GoDown(index, currentNode, nodeLabel, nbTip, tpc, k, edge)
		l = index[currentNode]
		nodeLabel[currentNode - nbTip] = tpc[k]
		k += 1
		currentNode = edge[l, 1]
		return k, currentNode
	end
	tp = "$tp;"
	if ismatch(r"^[^\(]+\(", tp)
		cutoff = length(match(r"^[^\(]+\(", tp).match)
		treeName = chop(match(r"^[^\(]+\(", tp).match)
		treeName = treeName[length(treeName)] == ' ' ? treeName[1:length(treeName)-1] : treeName
		tp = tp[cutoff:length(tp)]
	else treeName = ""
	end
	if search(tp, ",") == 0:-1
		edge = Array(Int, 2,2)
		edge[1,1:2] = [2,1]
		edge[2,1:2] = [1,2]
		tp = split(tp, r"[\\(\\):;]")
		edgeLength = tp[3]
		Nnode = 1
		tipLabel = tp[2]
		if tp[4] != ""
			nodeLabel = tp[4]
		end
		phyloobject = Phylo(edge, Nnode, tipLabel, edgeLength, nodeLabel, -1.0)
		return phyloobject
	end
	tsp = split(tp, "")
	tp = replace(tp, r"\s", "")
	tp = replace(tp, ")", ")NA")
	tp = replace(tp, "(", "rem(")
	tpc = split(tp, r"[\\(\\),;]")
	tpc = tpc[1:length(tpc)-1]
	tpc = tpc[tpc .!= "rem"]
	skeleton = tsp[bool([i == "(" || i == ")" || i == "," || i == ";" for i in tsp])]
	nsk = length(skeleton) # Length of the nexus skeleton.
	nbNode = sum(skeleton .== ")")
	nbTip = sum(skeleton .== ",") + 1
	nbEdge = nbNode + nbTip # Number of edges for tree is number of tips and nodes.
	nodeLabel = ["" for i in 1:nbNode]
	tipLabel = ["" for i in 1:nbTip]
	edge = Array(Int, nbEdge,2)
	currentNode = node = nbTip + 1
	edge[nbEdge, 1] = 0
	edge[nbEdge, 2] = node
	index = [0 for i in 1:nbEdge+1]
	index[node] = nbEdge
	j = k = tip = 1
	for i in 2:nsk
		if skeleton[i] == "("
			currentNode, node, j = AddInternal(edge, currentNode, node, index, j)
		end
		if skeleton[i] == "," && skeleton[i-1] != ")"
			currentNode, tip, k, j = AddTerminal(edge, currentNode, tip, index, tipLabel, tpc, k, j)
		end
		if skeleton[i] == ")"
			if skeleton[i - 1] == ","
				currentNode, tip, k, j = AddTerminal(edge, currentNode, tip, index, tipLabel, tpc, k, j)
				k, currentNode = GoDown(index, currentNode, nodeLabel, nbTip, tpc, k, edge)
			end
			if skeleton[i - 1] == ")"
				k, currentNode = GoDown(index, currentNode, nodeLabel, nbTip, tpc, k, edge)
			end
		end
	end
	edge = edge[1:nbEdge-1, 1:2]
        for i in 1:length(nodeLabel)
            nodeLabel[i] = replace(nodeLabel[i],r"^NA","")
        end
	phyloobject = Clado(treeName, edge, nbNode, tipLabel, nodeLabel)
	return phyloobject
end

# Sub function for creation of a Phylo structs from newick format strings.
# Used to build Phylo structs from newick strings during operation of the
# readtree function.
function treebuild(tp::ASCIIString)
	function AddInternal(edge, j, currentNode, node, index)
		edge[j, 1] = currentNode
		node += 1
        edge[j, 2] = currentNode = node
        index[node] = j
        j += 1
        return node, currentNode, j
	end
	function AddTerminal(edge, j, currentNode, tip, index, tpc, k, tipLabel, edgeLength)
		edge[j, 1] = currentNode
        edge[j, 2] = tip
        index[tip] = j
        X = split(tpc[k], ":")
        tipLabel[tip] = X[1]
        edgeLength[j] = X[2]
        k += 1
        tip +=1
        j += 1
        return k, tip, j
	end
	function GoDown(index, currentNode, tpc, k, nodeLabel, nbTip, edgeLength)
		l = index[currentNode]
        X = split(tpc[k], ":")
        nodeLabel[currentNode - nbTip] = X[1]
        edgeLength[l] = X[2]
        k += 1
        currentNode = edge[l, 1]
        return currentNode, k
	end
	tp = "$tp;"
	if ismatch(r"^[^\(]+\(", tp)
		cutoff = length(match(r"^[^\(]+\(", tp).match)
		treeName = chop(match(r"^[^\(]+\(", tp).match)
		treeName = treeName[length(treeName)] == ' ' ? treeName[1:length(treeName)-1] : treeName
		tp = tp[cutoff:length(tp)]
	else treeName = ""
	end
	if search(tp, ",") == 0:-1
		edge = Array(Int, 2,2)
		edge[1,1:2] = [2,1]
		edge[2,1:2] = [1,2]
		tp = split(tp, r"[\\(\\):;]")
		edgeLength = tp[3]
		Nnode = 1
		tipLabel = tp[2]
		if tp[4] != ""
			nodeLabel = tp[4]
		end
		phyloobject = Phylo(edge, Nnode, tipLabel, edgeLength, nodeLabel, -1.0)
		return phyloobject
	end
	tsp = split(tp, "")
	if ismatch(r"\)[^\(:\)\r\n]*;", tp) # Match the root if there is no colon.
		m = match(r"\)[^\(:\)\r\n]*;", tp)
		mstring = m.match[1:length(m.match)-1]
		newString = "$mstring:NA;"
		st1 = tp[1:m.offset-1]
		tp = "$st1$newString"
	end
	tp = replace(tp, ")", ")NA")
	tp = replace(tp, r"\s", "")
	tp = replace(tp, "(", "rem(")
	tpc = split(tp, r"[\\(\\),;]")
	tpc = tpc[1:length(tpc)-1]
	tpc = tpc[tpc .!= "rem"]
    skeleton = tsp[bool([i == "(" || i == ")" || i == "," || i == ";" for i in tsp])]
    nsk = length(skeleton)
    nbNode = sum(skeleton .== ")")
    nbTip = sum(skeleton .== ",") + 1
    nbEdge = nbNode + nbTip
    nodeLabel = ["" for i in 1:nbNode]
	tipLabel = ["" for i in 1:nbTip]
	edgeLength = ["" for i in 1:nbEdge]
    edge = Array(Int, nbEdge, 2)
    currentNode = node = nbTip + 1
    edge[nbEdge, 1] = 0
    edge[nbEdge, 2] = node
    index = [0 for i in 1:(nbEdge + 1)]
    index[node] = nbEdge
    j = k = tip = 1
    for i in 2:nsk
        if skeleton[i] == "("
            node, currentNode, j = AddInternal(edge, j, currentNode, node, index)
        end
        if skeleton[i] == ","
            if skeleton[i - 1] != ")"
                k, tip, j = AddTerminal(edge, j, currentNode, tip, index, tpc, k, tipLabel, edgeLength)
            end
        end
        if skeleton[i] == ")"
            if skeleton[i - 1] == ","
                k, tip, j = AddTerminal(edge, j, currentNode, tip, index, tpc, k, tipLabel, edgeLength)
                currentNode, k = GoDown(index, currentNode, tpc, k, nodeLabel, nbTip, edgeLength)
            end
            if skeleton[i - 1] == ")"
                currentNode, k = GoDown(index, currentNode, tpc, k, nodeLabel, nbTip, edgeLength)
            end
        end
    end
    edge = edge[1:nbEdge-1, 1:2]
    rootEdge = edgeLength[nbEdge]
    if rootEdge != "NA" # Resolve whether there is a rootedge for the tree.
        rootEdge = float64(rootEdge)
    else
        rootEdge = -1.0
    end
    edgeLength = float64([i == "" ? -1.0 : float64(i) for i in edgeLength[1:nbEdge-1]])
    for i in 1:length(nodeLabel)
        nodeLabel[i] = replace(nodeLabel[i],r"^NA","")
    end
    # nodeLabel = [replace(i, r"^NA", "") for i in nodeLabel]
    phyloobject = Phylo(treeName, edge, nbNode, tipLabel, edgeLength, nodeLabel, rootEdge)
    return phyloobject
end

# Function method for writing a single Phylogeny type to file.
function treewrite(tree::Phylogeny,
                   file::ASCIIString = "output.nwk",
                   append::Bool = false,
                   treeNames::Bool = true)
	output = newick(tree, treeNames)
	if append == true
		outstream = open(file, "a")
	else
		outstream = open(file, "w")
	end
	println(outstream, output)
	close(outstream)
end

# Function method for writing multiple trees to file. In order to do this they must be as an array.
function treewrite(tree::Array{Phylogeny},
                   file::ASCIIString = "output.nwk",
                   append::Bool = false,
                   treeNames::Bool = false)
	outputarray = Array(ASCIIString, length(tree))
	for i in 1:length(tree)
		outputarray[i] = newick(tree[i], treeNames)
	end
	if apppend == true
		outstream = open(file, "a")
	else
		outstream = open(file, "w")
	end
	for i in outputarray
		println(outstream, i)
	end
	close(outstream)
end

# function used by the function which write newick trees to check labels.
function checklabels(labels)
	labels = replace(labels, r"^[[:space:]\\(]+", "")
	labels = replace(labels, r"[[:space:]\\)]+$", "")
	labels = replace(labels, r"[[:space:]]", "_")
	labels = replace(labels, r"[,:;]", "")
	labels = replace(labels, r"[\\(\\)]", "-")
	labels = replace(labels, r"_{2,}", "_")
	labels = replace(labels, r"-{2,}", "-")
	return labels
end

# Function used by the functions which write newick trees.
function cp(x, k, STRING)
	STRING[k] = x
	k += 1
	return k, STRING
end

# Returns the internal node which is the root of the tree.
function getroot(children, parents)
	r = parents[[find(in(children, i)) == [1] ? false : true for i in parents]][1]
	return r
end

# This function returns an array of arrays, showing which children
# nodes have.
function getkids(phy::Phylogeny)
	N = length(phy.tipLabel)
	kids = Array(Array{Int64}, N + phy.Nnode)
	for i in 1:N + phy.Nnode
		logic = bool([n == i for n in phy.edge[1:size(phy.edge,1), 1]])
		kids[i] = phy.edge[logic,2]
	end
	return kids
end

# Construct a newick string from a Cladogram
function newick(phy::Clado, name::Bool)
	function addInternal(i, k, STRING, N, nodelab, tiplab, ind)
		k, STRING = cp("(", k, STRING)
		desc = kids[i]
		for j in desc
			if j > N
				STRING, k = addInternal(j, k, STRING, N, nodelab, tiplab)
			else
				k, STRING = addTerminal(ind[j], k, STRING, tiplab, children)
			end
			if j != desc[length(desc)]
				k, STRING = cp(",", k, STRING)
			end
		end
		k, STRING = cp(")", k, STRING)
		if i > N
			k, STRING = cp(nodelab[i - N], k, STRING)
		end
		return STRING, k
	end
	function addTerminal(i, k, STRING, tiplab, children)
		i = i[1]
		k, STRING = cp(tiplab[children[i]], k, STRING)
		return k, STRING
	end
	if name == true
		prefix = phy.name
	else
		prefix = ""
	end
	nodelab = [i != "" ? checklabels(i) : "" for i in phy.nodeLabel]
	tiplab = [i != "" ? checklabels(i) : i for i in phy.tipLabel]
	children = phy.edge[1:size(phy.edge,1), 2]
	parents = phy.edge[1:size(phy.edge,1), 1]
	N = length(phy.tipLabel)
	kids = getkids(phy)
	LS = (4 * N) + 5
	LS = LS + N # if there are nodelabels.
	ind = [findin(children, i) for i in 1:maximum(phy.edge)]
	STRING = ["" for i in 1:LS]
	k = 1
	k, STRING = cp(prefix, k, STRING)
	k, STRING = cp("(", k, STRING)
	desc = kids[getroot(children, parents)]
	for j in desc
		if j > N
			STRING, k = addInternal(j, k, STRING, N, nodelab, tiplab, ind)
		else
			k, STRING = addTerminal(ind[j], k, STRING, tiplab, children)
			if j != desc[length(desc)]
				k, STRING = cp(",", k, STRING)
			end
		end
	end
	k, STRING = cp(")", k, STRING)
	k, STRING = cp(nodelab[1], k, STRING)
	k, STRING = cp(";", k, STRING)
	outstring = ""
	for i in STRING
		outstring = "$outstring$i"
	end
	if name == true && prefix != ""
		namebit = chop(match(r"^[^\(]+\(", outstring).match)
		replace(outstring, namebit, "$namebit ")
	end
	return outstring
end

# Function that creates
function newick(phy::Phylo, name::Bool)
	function addInternal(i, k, STRING, N, nodelab, tiplab, ind)
		k, STRING = cp("(", k, STRING)
		desc = kids[i]
		for j in desc
			if j > N
				STRING, k = addInternal(j, k, STRING, N, nodelab, tiplab)
			else
				k, STRING = addTerminal(ind[j], k, STRING, tiplab, children)
			end
			if j != desc[length(desc)]
				k, STRING = cp(",", k, STRING)
			end
		end
		k, STRING = cp(")", k, STRING)
		if i > N
			k, STRING = cp(nodelab[i - N], k, STRING)
		end
		k, STRING = cp(":", k, STRING)
		edgechar = phy.edgeLength[ind[i]][1]
		k, STRING = cp("$edgechar", k, STRING)
		return STRING, k
	end
	function addTerminal(i, k, STRING, tiplab, children)
		i = i[1]
		k, STRING = cp(tiplab[children[i]], k, STRING)
		k, STRING = cp(":", k, STRING)
		edgechar = phy.edgeLength[i]
		k, STRING = cp("$edgechar", k, STRING)
		return k, STRING
	end
	if name == true
		prefix = phy.name
	else
		prefix = ""
	end
	brl = phy.edgeLength
	nodelab = [i != "" ? checklabels(i) : "" for i in phy.nodeLabel]
	tiplab = [i != "" ? checklabels(i) : i for i in phy.tipLabel]
	children = phy.edge[1:size(phy.edge,1), 2]
	parents = phy.edge[1:size(phy.edge,1), 1]
	N = length(phy.tipLabel)
	kids = getkids(phy)
	LS = (4 * N) + 5
	LS = LS + N
	LS = LS +(4 * N)
	ind = [findin(children, i) for i in 1:maximum(phy.edge)]
	STRING = ["" for i in 1:LS]
	k = 1
	k, STRING = cp(prefix, k, STRING)
	k, STRING = cp("(", k, STRING)
	desc = kids[getroot(children, parents)]
	for j in desc
		if j > N
			STRING, k = addInternal(j, k, STRING, N, nodelab, tiplab, ind)
		else
			k, STRING = addTerminal(ind[j], k, STRING, tiplab, children)
			if j != desc[length(desc)]
				k, STRING = cp(",", k, STRING)
			end
		end
	end
	k, STRING = cp(")", k, STRING)
	k, STRING = cp(nodelab[1], k, STRING)
	k, STRING = cp(";", k, STRING)
	outstring = ""
	for i in STRING
		outstring = "$outstring$i"
	end
	if name == true && prefix != ""
		namebit = chop(match(r"^[^\(]+\(", outstring).match)
		replace(outstring, namebit, "$namebit ")
	end
	return outstring
end

# Macro for making a tree from a string.
macro tr_str(s)
	s = replace(s, r"(\r|\n|\s)", "")
	if search(s, ";") != 0:-1
		s = s[1:length(s)-1]
	end
	if search(s, ":") == 0:-1
		tree = cladobuild(s)
		return tree
	elseif search(s, ":") != 0:-1
		tree = treebuild(s)
		return tree
	end
end


text:
export BrickletIndustrialDualRelayValue
struct BrickletIndustrialDualRelayValue
    channel0::Bool
    channel1::Bool
end

export BrickletIndustrialDualRelayMonoflop
struct BrickletIndustrialDualRelayMonoflop
    value::Bool
    time::Integer
    time_remaining::Integer
end

export BrickletIndustrialDualRelaySPITFPErrorCount
struct BrickletIndustrialDualRelaySPITFPErrorCount
    error_count_ack_checksum::Integer
    error_count_message_checksum::Integer
    error_count_frame::Integer
    error_count_overflow::Integer
end

export BrickletIndustrialDualRelayIdentity
struct BrickletIndustrialDualRelayIdentity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletIndustrialDualRelay
"""
Two relays to switch AC/DC devices
"""
mutable struct BrickletIndustrialDualRelay <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletIndustrialDualRelay(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_industrial_dual_relay")
        deviceInternal = package.BrickletIndustrialDualRelay(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export set_value
"""
    $(SIGNATURES)

Sets the state of the relays, *true* means on and *false* means off.
For example: (true, false) turns relay 0 on and relay 1 off.

If you just want to set one of the relays and don't know the current state
of the other relay, you can get the state with :func:`Get Value` or you
can use :func:`Set Selected Value`.

All running monoflop timers will be aborted if this function is called.
"""
function set_value(device::BrickletIndustrialDualRelay, channel0, channel1)
    device.deviceInternal.set_value(channel0, channel1)
end

export get_value
"""
    $(SIGNATURES)

Returns the state of the relays, *true* means on and *false* means off.
"""
function get_value(device::BrickletIndustrialDualRelay)
    return device.deviceInternal.get_value()
end

export set_monoflop
"""
    $(SIGNATURES)

The first parameter can be 0 or 1 (relay 0 or relay 1). The second parameter
is the desired state of the relay (*true* means on and *false* means off).
The third parameter indicates the time that the relay should hold
the state.

If this function is called with the parameters (1, true, 1500):
Relay 1 will turn on and in 1.5s it will turn off again.

A monoflop can be used as a failsafe mechanism. For example: Lets assume you
have a RS485 bus and a Industrial Dual Relay Bricklet connected to one of the
slave stacks. You can now call this function every second, with a time parameter
of two seconds. The relay will be on all the time. If now the RS485
connection is lost, the relay will turn off in at most two seconds.
"""
function set_monoflop(device::BrickletIndustrialDualRelay, channel, value, time)
    device.deviceInternal.set_monoflop(channel, value, time)
end

export get_monoflop
"""
    $(SIGNATURES)

Returns (for the given relay) the current state and the time as set by
:func:`Set Monoflop` as well as the remaining time until the state flips.

If the timer is not running currently, the remaining time will be returned
as 0.
"""
function get_monoflop(device::BrickletIndustrialDualRelay, channel)
    return device.deviceInternal.get_monoflop(channel)
end

export set_selected_value
"""
    $(SIGNATURES)

Sets the state of the selected relay, *true* means on and *false*
means off.

A running monoflop timer for the selected relay will be aborted if this function
is called.

The other relay remains untouched.
"""
function set_selected_value(device::BrickletIndustrialDualRelay, channel, value)
    device.deviceInternal.set_selected_value(channel, value)
end

export get_spitfp_error_count
"""
    $(SIGNATURES)

Returns the error count for the communication between Brick and Bricklet.

The errors are divided into

* ACK checksum errors,
* message checksum errors,
* framing errors and
* overflow errors.

The errors counts are for errors that occur on the Bricklet side. All
Bricks have a similar function that returns the errors on the Brick side.
"""
function get_spitfp_error_count(device::BrickletIndustrialDualRelay)
    return device.deviceInternal.get_spitfp_error_count()
end

export set_bootloader_mode
"""
    $(SIGNATURES)

Sets the bootloader mode and returns the status after the requested
mode change was instigated.

You can change from bootloader mode to firmware mode and vice versa. A change
from bootloader mode to firmware mode will only take place if the entry function,
device identifier and CRC are present and correct.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_bootloader_mode(device::BrickletIndustrialDualRelay, mode)
    return device.deviceInternal.set_bootloader_mode(mode)
end

export get_bootloader_mode
"""
    $(SIGNATURES)

Returns the current bootloader mode, see :func:`Set Bootloader Mode`.
"""
function get_bootloader_mode(device::BrickletIndustrialDualRelay)
    return device.deviceInternal.get_bootloader_mode()
end

export set_write_firmware_pointer
"""
    $(SIGNATURES)

Sets the firmware pointer for :func:`Write Firmware`. The pointer has
to be increased by chunks of size 64. The data is written to flash
every 4 chunks (which equals to one page of size 256).

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_write_firmware_pointer(device::BrickletIndustrialDualRelay, pointer)
    device.deviceInternal.set_write_firmware_pointer(pointer)
end

export write_firmware
"""
    $(SIGNATURES)

Writes 64 Bytes of firmware at the position as written by
:func:`Set Write Firmware Pointer` before. The firmware is written
to flash every 4 chunks.

You can only write firmware in bootloader mode.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function write_firmware(device::BrickletIndustrialDualRelay, data)
    return device.deviceInternal.write_firmware(data)
end

export set_status_led_config
"""
    $(SIGNATURES)

Sets the status LED configuration. By default the LED shows
communication traffic between Brick and Bricklet, it flickers once
for every 10 received data packets.

You can also turn the LED permanently on/off or show a heartbeat.

If the Bricklet is in bootloader mode, the LED is will show heartbeat by default.
"""
function set_status_led_config(device::BrickletIndustrialDualRelay, config)
    device.deviceInternal.set_status_led_config(config)
end

export get_status_led_config
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Status LED Config`
"""
function get_status_led_config(device::BrickletIndustrialDualRelay)
    return device.deviceInternal.get_status_led_config()
end

export get_chip_temperature
"""
    $(SIGNATURES)

Returns the temperature as measured inside the microcontroller. The
value returned is not the ambient temperature!

The temperature is only proportional to the real temperature and it has bad
accuracy. Practically it is only useful as an indicator for
temperature changes.
"""
function get_chip_temperature(device::BrickletIndustrialDualRelay)
    return device.deviceInternal.get_chip_temperature()
end

export reset
"""
    $(SIGNATURES)

Calling this function will reset the Bricklet. All configurations
will be lost.

After a reset you have to create new device objects,
calling functions on the existing ones will result in
undefined behavior!
"""
function reset(device::BrickletIndustrialDualRelay)
    device.deviceInternal.reset()
end

export write_uid
"""
    $(SIGNATURES)

Writes a new UID into flash. If you want to set a new UID
you have to decode the Base58 encoded UID string into an
integer first.

We recommend that you use Brick Viewer to change the UID.
"""
function write_uid(device::BrickletIndustrialDualRelay, uid)
    device.deviceInternal.write_uid(uid)
end

export read_uid
"""
    $(SIGNATURES)

Returns the current UID as an integer. Encode as
Base58 to get the usual string version.
"""
function read_uid(device::BrickletIndustrialDualRelay)
    return device.deviceInternal.read_uid()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletIndustrialDualRelay)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletIndustrialDualRelay, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
"""
    calculate the mean free path
"""

using Unitful

include("inputting.jl")
include("constants.jl")
include("dviscosity.jl")

_ƒμ = μ

"""
    λg(T, p, M_wt, ..., μ, µ₀, T₀, ...)

Calculate the mean free path based on gas kinetic theory.

Inputs:
T       Temperature             [K]         (default: 298.15)
p       Pressure                [Pa]        (default: 101325)
M_wt    Molar mass              [g/mol]     (default: M_wt_air)
...
μ       Viscosity               [kg/m/s]    (default: μ(T, μ₀, T₀))

Output:
        Mean free path          [m]

Output:
    Dynamic viscosity       [kg/m/s]

# Notes:
* M_wt_air comes from the the constants.jl file.
* Two additional inputs are possible:
    µ₀  Reference viscosity     [kg/m/s]    (default: µ₀_air_stp)
    T₀  Reference temperature   [K]         (default: T₀_air_stp)
    where  µ₀_air_stp and T₀_air_stp come from the constants.jl file.
* μ(T, μ₀, T₀) is defined in dviscosity.jl

# Examples:
```julia-repl
julia> λg()
6.647984982685411e-8 m
julia> using Unitful
julia> uconvert(u"nm", λg(T = 300, p = 1e5))
67.8918099704083 nm
```
"""
function λg(; kwargs...)
    T = ℞T(get(kwargs, :T, 298.15))
    p = ℞p(get(kwargs, :p, 101325))
    M_wt = ℞M_wt(get(kwargs, :M_wt, M_air))
    μ = _ƒμ(; kwargs...)
    μ = ℞μ(get(kwargs, :μ, μ))
    return (2 * μ / p) / (8 * M_wt / (π * R * T)).^ 0.5 |> upreferred
end


text:
function main(r::Robot)
    path = go_to_west_south_corner_and_return_path!(r; go_around_barriers = true)

    x = go_to_border_come_back_and_return_distance!(r, East; go_around_barriers = true)
    marks_to_do = x
    while !isborder(r,Nord) && marks_to_do > 0
        go!(r, East; steps = marks_to_do, go_around_barriers = true, markers = true)
        go!(r, West; steps = marks_to_do, go_around_barriers = true, markers = false)
        move!(r,Nord)
        marks_to_do -= 1
    end
    go!(r, East; steps = marks_to_do, go_around_barriers = true, markers = true)
    go!(r, West; steps = marks_to_do, go_around_barriers = true, markers = false)

    go_to_west_south_corner_and_return_path!(r; go_around_barriers = true)
    go_by_path!(r, path)
end

max_stars_repo_name: Savko-fokin/mirea-progs

max_issues_repo_name: Savko-fokin/mirea-progs

max_forks_repo_name: Savko-fokin/mirea-progs


text:
using Unitful#, UnitfulPlots
using OrdinaryDiffEq, DiffEqBase

const UNITS_TEST_ALGS = [Euler(),Midpoint(),RK4(),SSPRK104(),SSPRK22(),SSPRK33(),
    SSPRK432(),BS3(),BS5(),DP5(),DP5Threaded(),DP8(),Feagin10(),Feagin12(),Feagin14(),
    TanYam7(),Tsit5(),TsitPap8(),Vern6(),Vern7(),Vern8(),Vern9()]

f = (y,p,t) -> 0.5*y / 3.0u"s"
u0 = 1.0u"N"
prob = ODEProblem(f,u0,(0.0u"s",1.0u"s"))

sol =solve(prob,ExplicitRK())

for alg in UNITS_TEST_ALGS
  if !(typeof(alg) <: DP5Threaded)
    @show alg
    sol = solve(prob,alg,dt=1u"s"/10)
  end
end

println("Units for Number pass")

u0 = [1.0u"N" 2.0u"N"
      3.0u"N" 1.0u"N"]
f = (dy,y,p,t) -> (dy .= 0.5.*y ./ 3.0u"s")
prob = ODEProblem(f,u0,(0.0u"s",1.0u"s"))

sol =solve(prob,ExplicitRK())

for alg in UNITS_TEST_ALGS
  @show alg
  sol = solve(prob,alg,dt=1u"s"/10)
end

println("Units for 2D pass")

using Unitful, RecursiveArrayTools, DiffEqBase, OrdinaryDiffEq
using LinearAlgebra
r0 = [1131.340, -2282.343, 6672.423]u"km"
v0 = [-5.64305, 4.30333, 2.42879]u"km/s"
Δt = 86400.0*365u"s"
μ = 398600.4418u"km^3/s^2"
rv0 = ArrayPartition(r0,v0)

function f(dy, y, μ, t)
    r = norm(y.x[1])
    dy.x[1] .= y.x[2]
    dy.x[2] .= -μ .* y.x[1] / r^3
end

prob = ODEProblem(f,rv0,(0.0u"s",1.0u"s"),μ)
sol = solve(prob,Tsit5())

# coordinate: u = [position, momentum]
# parameters: p = [mass, force constanst]
function f_harmonic!(du,u,p,t)
  du[1] = u[2]/p[1]
  du[2] = -p[2]*u[1]
end

mass = 1.0u"kg"
k = 1.0u"N/m"
p = [mass, k]

u0 = [1.0u"m", 0.0u"kg*m/s"] # initial values (position, momentum)
tspan = (0.0u"s", 10.0u"s")
prob = ODEProblem(f_harmonic!, u0, tspan, p)
sol = solve(prob)


text:
############# data preparation ###################

function intended(d::AbstractString)
    m = match(r"\((.+),(.+)\)", d)
    x, y = parse.(Int, m.captures)
    DungBase.Point(x, y)
end

function parsetitle(title, r)
    run = r.data
    return (nest = run.nest,
            feeder = run.feeder,
            fictive_nest = run.fictive_nest,
            track = run.track,
            dropoff = run.dropoff,
            pickup = run.pickup,
            title = title,
            comment = r.metadata.comment,
            nest2feeder = get(r.metadata.setup, :nest2feeder, missing),
            experience = get(r.metadata.setup, :experience, missing),
            pickup_loc = get(r.metadata.setup, :pickup, missing),
            dropoff_loc = get(r.metadata.setup, :dropoff, missing),
            displacement = passmissing(intended)(get(r.metadata.setup, :displacement, missing)))
end

function getdf(data)
    df = DataFrame(parsetitle(k, r) for (k, v) in data for r in v.runs)

    # @. df[!, :displace_direction] = switchdirections(df.displace_direction)
    @. df[!, :group] = getgroup(df.pickup_loc, df.displacement, df.dropoff_loc)
    # @. df[!, :set] = getset(df.transfer, df.group)

    categorical!(df, [:experience, :pickup_loc, :dropoff_loc])
    # levels!(df.group, ["none", "left", "right", "away", "towards", "zero", "back", "far"])

    # df[!, :direction_deviation]  = [angle(r.fictive_nest - r.feeder, turningpoint(r.track) - r.feeder) for r in eachrow(df)]
    # max_direction_deviation = maximum(r.direction_deviation for r in eachrow(df) if r.group ∉ ("far", "zero"))
    # mean_direction_deviation = mean(r.direction_deviation for r in eachrow(df) if r.group ∉ ("far", "zero"))
    # filter!(r -> r.group ≠ "far" || r.direction_deviation < 4mean_direction_deviation, df)

    df[!, :turning_point] .= zero.(df.feeder)
    df[!, :center_of_search] .= zero.(df.feeder)
    for r in eachrow(df)
        trans = createtrans(r.nest, r.dropoff, r.fictive_nest)
        @. r.track.coords = trans(r.track.coords)
        @. r.track.rawcoords.xy .= trans(r.track.rawcoords.xy)
        r.feeder = trans(r.feeder)
        r.fictive_nest = trans(r.fictive_nest)
        r.nest = trans(r.nest)
        r.pickup = trans(r.pickup)
        r.dropoff = trans(r.dropoff)
        Δ = r.displacement - r.fictive_nest
        # r.turning_point = turningpoint(r.track) + Δ
        # r.center_of_search = searchcenter(r.track) + Δ
    end

    groups = levels(df.group)
    nc = length(groups)
    colors = OrderedDict(zip(groups, [colorant"black"; distinguishable_colors(nc - 1, [colorant"white", colorant"black"], dropseed = true)]))

    gdf = groupby(df, :group)
    DataFrames.transform!(gdf, :group => (g -> colors[g[1]]) => :groupcolor)
    DataFrames.transform!(gdf, :groupcolor => getcolor => :color)

    df
end

__x(x) = x == 0 ? "" :
         x < 0 ? "left" :
         "right"
__y(y) = y == 0 ? "" :
         y < 0 ? "away" :
         "towards"
function _f(x, y)
    x = __x(x)
    y = __y(y)
    join(filter(!isempty, [x, y]), " ")
end
_f(xy) = _f(xy...)
_f(::Missing) = ""
getgroup(pickup, displacement, dropoff) = string(pickup, " ", _f(displacement), " ", dropoff)

#=switchdirections(_::Missing) = missing
switchdirections(d) =   d == "left" ? "right" :
                        d == "right" ? "left" :
                        d=#

# getgroup(displace_location::Missing, transfer, displace_direction) = transfer
# getgroup(displace_location, transfer, displace_direction) = displace_location == "nest" ? "zero" : displace_direction
# getset(_::Missing, d) = d == "none" ? "Closed" : "Displacement"
# getset(_, __) = "Transfer"

# intended(::Missing) = missing
# intended(d) = intended(string(d))

_get_rotation_center(nest::Missing, fictive_nest) = fictive_nest
_get_rotation_center(nest, fictive_nest) = nest
function createtrans(nest, dropoff, fictive_nest)
    v = dropoff - _get_rotation_center(nest, fictive_nest)
    α = atan(v[2], v[1])
    rot = LinearMap(Angle2d(-π/2 - α))
    trans = Translation(-_get_rotation_center(nest, fictive_nest))
    passmissing(rot ∘ trans)
end

function highlight(c, i, n)
    h = HSL(c)
    HSL(h.h, h.s, i/(n + 1))
end
function getcolor(g)
    n = length(g)
    [highlight(c, i, n) for (i, c) in enumerate(g)]
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct AutomaticRepairsPolicy <: SwaggerModel
    enabled::Any # spec type: Union{ Nothing, Bool } # spec name: enabled
    gracePeriod::Any # spec type: Union{ Nothing, String } # spec name: gracePeriod

    function AutomaticRepairsPolicy(;enabled=nothing, gracePeriod=nothing)
        o = new()
        validate_property(AutomaticRepairsPolicy, Symbol("enabled"), enabled)
        setfield!(o, Symbol("enabled"), enabled)
        validate_property(AutomaticRepairsPolicy, Symbol("gracePeriod"), gracePeriod)
        setfield!(o, Symbol("gracePeriod"), gracePeriod)
        o
    end
end # type AutomaticRepairsPolicy

const _property_map_AutomaticRepairsPolicy = Dict{Symbol,Symbol}(Symbol("enabled")=>Symbol("enabled"), Symbol("gracePeriod")=>Symbol("gracePeriod"))
const _property_types_AutomaticRepairsPolicy = Dict{Symbol,String}(Symbol("enabled")=>"Bool", Symbol("gracePeriod")=>"String")
Base.propertynames(::Type{ AutomaticRepairsPolicy }) = collect(keys(_property_map_AutomaticRepairsPolicy))
Swagger.property_type(::Type{ AutomaticRepairsPolicy }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_AutomaticRepairsPolicy[name]))}
Swagger.field_name(::Type{ AutomaticRepairsPolicy }, property_name::Symbol) =  _property_map_AutomaticRepairsPolicy[property_name]

function check_required(o::AutomaticRepairsPolicy)
    true
end

function validate_property(::Type{ AutomaticRepairsPolicy }, name::Symbol, val)
end


text:
using GPUifyLoops: @launch, @loop
using Oceananigans.Operators: Δx, Δy, ΔzF
using Oceananigans.Utils: @loop_xy, @loop_xz, @loop_yz, launch_config

#####
##### Algorithm for adding fluxes associated with non-trivial flux boundary conditions.
##### Value and Gradient boundary conditions are handled by filling halos.
#####

"""
    apply_x_bcs!(Gc, arch, grid, args...)

Apply flux boundary conditions to a field `c` by adding the associated flux divergence to
the source term `Gc` at the left and right.
"""
@inline function apply_x_bcs!(Gc, c, arch, args...)
    @launch(device(arch), config=launch_config(Gc.grid, :yz),
            _apply_x_bcs!(Gc.data, Gc.grid, c.boundary_conditions.x.left,
                                            c.boundary_conditions.x.right, args...))
    return nothing
end

"""
    apply_y_bcs!(Gc, arch, grid, args...)

Apply flux boundary conditions to a field `c` by adding the associated flux divergence to
the source term `Gc` at the left and right.
"""
@inline function apply_y_bcs!(Gc, c, arch, args...)
    @launch(device(arch), config=launch_config(Gc.grid, :xz),
            _apply_y_bcs!(Gc.data, Gc.grid, c.boundary_conditions.y.left,
                                            c.boundary_conditions.y.right, args...))
    return nothing
end

"""
    apply_z_bcs!(Gc, arch, grid, args...)

Apply flux boundary conditions to a field `c` by adding the associated flux divergence to
the source term `Gc` at the top and bottom.
"""
@inline function apply_z_bcs!(Gc, c, arch, args...)
    @launch(device(arch), config=launch_config(Gc.grid, :xy),
            _apply_z_bcs!(Gc.data, Gc.grid, c.boundary_conditions.z.bottom,
                                            c.boundary_conditions.z.top, args...))
    return nothing
end

# Avoid some computation / memory accesses for Value, Gradient, Periodic, NormalFlow,
# and zero-flux boundary conditions --- every boundary condition that does *not* prescribe
# a non-trivial flux.
const NotFluxBC = Union{VBC, GBC, PBC, NFBC, ZFBC}

@inline _apply_x_bcs!(Gc, grid, ::NotFluxBC, ::NotFluxBC, args...) = nothing
@inline _apply_y_bcs!(Gc, grid, ::NotFluxBC, ::NotFluxBC, args...) = nothing
@inline _apply_z_bcs!(Gc, grid, ::NotFluxBC, ::NotFluxBC, args...) = nothing

"""
    _apply_x_bcs!(Gc, grid, west_bc, east_bc, args...)

Apply a west and/or east boundary condition to variable `c`.
"""
function _apply_x_bcs!(Gc, grid, west_bc, east_bc, args...)
    @loop_yz j k grid begin
        apply_x_west_bc!(Gc, west_bc, j, k, grid, args...)
        apply_x_east_bc!(Gc, east_bc, j, k, grid, args...)
    end
end

"""
    _apply_y_bcs!(Gc, grid, south_bc, north_bc, args...)

Apply a south and/or north boundary condition to variable `c`.
"""
function _apply_y_bcs!(Gc, grid, south_bc, north_bc, args...)
    @loop_xz i k grid begin
        apply_y_south_bc!(Gc, south_bc, i, k, grid, args...)
        apply_y_north_bc!(Gc, north_bc, i, k, grid, args...)
    end
end

"""
    _apply_z_bcs!(Gc, grid, bottom_bc, top_bc, args...)

Apply a top and/or bottom boundary condition to variable `c`.
"""
function _apply_z_bcs!(Gc, grid, bottom_bc, top_bc, args...)
    @loop_xy i j grid begin
        apply_z_bottom_bc!(Gc, bottom_bc, i, j, grid, args...)
           apply_z_top_bc!(Gc, top_bc,    i, j, grid, args...)
    end
end

# Fall back functions for boundary conditions that are not of type Flux.
@inline apply_x_east_bc!(  args...) = nothing
@inline apply_x_west_bc!(  args...) = nothing
@inline apply_y_north_bc!( args...) = nothing
@inline apply_y_south_bc!( args...) = nothing
@inline apply_z_top_bc!(   args...) = nothing
@inline apply_z_bottom_bc!(args...) = nothing

# Shortcuts for 'zero' flux boundary conditions.
@inline apply_x_east_bc!(  Gc, ::ZFBC, args...) = nothing
@inline apply_x_west_bc!(  Gc, ::ZFBC, args...) = nothing
@inline apply_y_north_bc!( Gc, ::ZFBC, args...) = nothing
@inline apply_y_south_bc!( Gc, ::ZFBC, args...) = nothing
@inline apply_z_top_bc!(   Gc, ::ZFBC, args...) = nothing
@inline apply_z_bottom_bc!(Gc, ::ZFBC, args...) = nothing

"""
    apply_x_west_bc!(Gc, west_flux::BC{<:Flux}, j, k, grid, args...)

Add the flux divergence associated with a west flux boundary condition on `c`.
Note that because

    `tendency = ∂c/∂t = Gc = - ∇ ⋅ flux`

a positive west flux is associated with an *increase* in `Gc` near the west boundary.
If `west_bc.condition` is a function, the function must have the signature

    `west_bc.condition(j, k, grid, boundary_condition_args...)`

The same logic holds for south and bottom boundary conditions in `y`, and `z`, respectively.
"""
@inline apply_x_west_bc!(  Gc,   west_flux::BC{<:Flux}, j, k, grid, args...) = @inbounds Gc[1, j, k] += getbc(west_flux,   j, k, grid, args...) /  Δx(1, j, k, grid)
@inline apply_y_south_bc!( Gc,  south_flux::BC{<:Flux}, i, k, grid, args...) = @inbounds Gc[i, 1, k] += getbc(south_flux,  i, k, grid, args...) /  Δy(i, 1, k, grid)
@inline apply_z_bottom_bc!(Gc, bottom_flux::BC{<:Flux}, i, j, grid, args...) = @inbounds Gc[i, j, 1] += getbc(bottom_flux, i, j, grid, args...) / ΔzF(i, j, 1, grid)

"""
    apply_x_east_bc!(Gc, top_flux::BC{<:Flux}, j, k, grid, args...)

Add the part of flux divergence associated with a east boundary condition on `c`.
Note that because

    `tendency = ∂c/∂t = Gc = - ∇ ⋅ flux`

a positive east flux is associated with a *decrease* in `Gc` near the east boundary.
If `east_bc.condition` is a function, the function must have the signature

    `east_bc.condition(i, j, grid, boundary_condition_args...)`

The same logic holds for north and top boundary conditions in `y`, and `z`, respectively.
"""
@inline apply_x_east_bc!( Gc,  east_flux::BC{<:Flux}, j, k, grid, args...) = @inbounds Gc[grid.Nx, j, k] -= getbc(east_flux,  j, k, grid, args...) /  Δx(grid.Nx, j, k, grid)
@inline apply_y_north_bc!(Gc, north_flux::BC{<:Flux}, i, k, grid, args...) = @inbounds Gc[i, grid.Ny, k] -= getbc(north_flux, i, k, grid, args...) /  Δy(i, grid.Ny, k, grid)
@inline apply_z_top_bc!(  Gc,   top_flux::BC{<:Flux}, i, j, grid, args...) = @inbounds Gc[i, j, grid.Nz] -= getbc(top_flux,   i, j, grid, args...) / ΔzF(i, j, grid.Nz, grid)


text:
module sti
using ValidatedNumerics;

pos = x -> x >= zero(x) ? x : zero(x)
neg = x -> x < zero(x) ? -x : zero(x)

function constituentMatrix{T}(matrix :: Array{T, 2})
    @assert ndims(matrix) == 2
    @assert size(matrix, 1) == size(matrix, 2)

    positivMatrix = map(pos, matrix)
    negativMatrix = map(neg, matrix)
    [
        positivMatrix negativMatrix;
        negativMatrix positivMatrix
    ] :: Array{T, 2}
end

function STI{T}(intervalVector :: Array{Interval{T}})
    @assert ndims(intervalVector) == 1

    function stiReducer(reduced, element)
        reducedSize = size(reduced, 1)
        if (reducedSize == 0)
            return [-element.lo; element.hi]
        end
        halfReducedSize = Int(reducedSize/2)
        vcat(
            slice(reduced, 1:halfReducedSize),
            [-element.lo],
            slice(reduced, (1 + halfReducedSize):reducedSize),
            [element.hi]
        ) # Could use cat instead
    end

    reduce(stiReducer, [], intervalVector) :: Array{T}
end

function reverseSTI{T}(intervalVector :: Array{T})
    @assert ndims(intervalVector) == 1

    vectorSize = size(intervalVector, 1)
    halfVectorSize = Int(vectorSize/2)
    negativLoIntervalPart = slice(intervalVector, 1:halfVectorSize)
    positivHiIntervalPart = slice(intervalVector, (1 + halfVectorSize):vectorSize)

    positivLoIntervalPart = map(-, negativLoIntervalPart)
    zippedIntervalParts = zip(positivLoIntervalPart, positivHiIntervalPart)
    collect(Interval{T}, map(x -> safeInterval(x[1], x[2]), zippedIntervalParts)) :: Array{Interval{T}}
end

function safeInterval(a, b)
    a < b ? @interval(a, b) : @interval(b, a)
end

end


text:
using Documenter, Neurthino

makedocs(;
    modules = [Neurthino],
    authors = "Johannes Schumann, Tamas Gal",
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
        assets = ["assets/logo.ico"],
    ),
    pages=[
        "Introduction" => "index.md",
        "API" => "api.md",
    ],
    repo="https://github.com/KM3NeT/Neurthino.jl/blob/{commit}{path}#L{line}",
    sitename="Neurthino.jl",
)

deploydocs(;
    repo="github.com/KM3NeT/Neurthino.jl",
)


text:
using WebIO
using Test

# NOTE: We use `@eval` a lot here for two reasons.
#   1. Because `struct` definitions are only allowed at the top level.
#   2. To delay execution of the `@register_renderable` macro so that the macro
#      evaluation happens at test time (this avoids invalid invocations from
#      preventing the rest of the test suite from even running and is also
#      necessary to try to catch errors raised by the macro).

@testset "WebIO.@register_renderable" begin

    @testset "@register_renderable for symbol" begin
        MyTypeName = gensym()
        @eval struct $MyTypeName dom::WebIO.Node end
        @eval WebIO.render(x::$MyTypeName) = x.dom
        @eval WebIO.@register_renderable($MyTypeName)

        MyType = @eval $MyTypeName
        @test hasmethod(show, (IO, WebIO.WEBIO_NODE_MIME, MyType))
        myinstance = MyType(node(:p, "Hello, world!"))
        myinstance_json = sprint(show, WebIO.WEBIO_NODE_MIME(), myinstance)
        myinstance_html = sprint(show, MIME("text/html"), myinstance)
        @test occursin("Hello, world!", myinstance_json)
        @test occursin("Hello, world!", myinstance_html)
        @test MyType in WebIO.renderable_types
    end

    @testset "@register_renderable with do-block syntax" begin
        MyTypeName = gensym()
        @eval struct $MyTypeName dom::WebIO.Node end

        MyType = @eval $MyTypeName
        @eval @WebIO.register_renderable($MyType) do x
            return x.dom
        end

        @test hasmethod(show, (IO, WebIO.WEBIO_NODE_MIME, MyType))
        myinstance = MyType(node(:p, "Hello, world!"))
        myinstance_json = sprint(show, WebIO.WEBIO_NODE_MIME(), myinstance)
        myinstance_html = sprint(show, MIME("text/html"), myinstance)
        @test occursin("Hello, world!", myinstance_json)
        @test occursin("Hello, world!", myinstance_html)
    end

    @testset "@register_renderable multiple dispatch" begin
        # This test is meant to make sure that the methods generated refer to
        # the correct types.
        TypeAName, TypeBName = gensym("TypeA"), gensym("TypeB")
        @eval struct $TypeAName end
        @eval struct $TypeBName end
        TypeA, TypeB = (@eval $TypeAName, @eval $TypeBName)
        @eval @WebIO.register_renderable($TypeA) do typea
            return "Type A"
        end
        @eval @WebIO.register_renderable($TypeB) do typeb
            return "Type B"
        end

        @test WebIO.render(TypeA()) == "Type A"
        @test WebIO.render(TypeB()) == "Type B"
    end
end


text:
export K1_1StructuredParams

"""
Type to use the K1.1 formulation with a structured Krylov method, using the package
[`Krylov.jl`](https://github.com/JuliaSmoothOptimizers/Krylov.jl).
This only works for solving Linear Problems.
The outer constructor

    K1_1StructuredParams(; uplo = :L, kmethod = :lsqr, rhs_scale = true,
                         atol0 = 1.0e-4, rtol0 = 1.0e-4,
                         atol_min = 1.0e-10, rtol_min = 1.0e-10,
                         ρ_min = 1e3 * sqrt(eps()), δ_min = 1e4 * sqrt(eps()),
                         mem = 20)

creates a [`RipQP.SolverParams`](@ref) that should be used to create a [`RipQP.InputConfig`](@ref).
The available methods are:
- `:lslq`
- `:lsqr`
- `:lsmr`

"""
mutable struct K1_1StructuredParams <: NormalParams
  uplo::Symbol
  kmethod::Symbol
  rhs_scale::Bool
  atol0::Float64
  rtol0::Float64
  atol_min::Float64
  rtol_min::Float64
  ρ_min::Float64
  δ_min::Float64
  mem::Int
end

function K1_1StructuredParams(;
  uplo::Symbol = :L,
  kmethod::Symbol = :lsqr,
  rhs_scale::Bool = true,
  atol0::T = 1.0e-4,
  rtol0::T = 1.0e-4,
  atol_min::T = 1.0e-10,
  rtol_min::T = 1.0e-10,
  ρ_min::T = 1e3 * sqrt(eps()),
  δ_min::T = 1e4 * sqrt(eps()),
  mem::Int = 20,
) where {T <: Real}
  return K1_1StructuredParams(
    uplo,
    kmethod,
    rhs_scale,
    atol0,
    rtol0,
    atol_min,
    rtol_min,
    ρ_min,
    δ_min,
    mem,
  )
end

mutable struct PreallocatedDataK1_1Structured{T <: Real, S, Ksol <: KrylovSolver} <:
               PreallocatedDataNormalStructured{T, S}
  E::S  # temporary top-left diagonal
  invE::S
  ξ1::S
  ξ2::S # todel
  Δy0::S
  ξ12::S # todel
  rhs_scale::Bool
  regu::Regularization{T}
  KS::Ksol
  atol::T
  rtol::T
  atol_min::T
  rtol_min::T
end

function PreallocatedData(
  sp::K1_1StructuredParams,
  fd::QM_FloatData{T},
  id::QM_IntData,
  itd::IterData{T},
  pt::Point{T},
  iconf::InputConfig{Tconf},
) where {T <: Real, Tconf <: Real}
  @assert iconf.solve_method == :IPF

  # init Regularization values
  E = similar(fd.c, id.nvar)
  if iconf.mode == :mono
    regu =
      Regularization(T(sqrt(eps()) * 1e5), T(sqrt(eps()) * 1e5), T(sp.ρ_min), T(sp.δ_min), :classic)
    E .= T(1.0e0) / 2
  else
    regu = Regularization(
      T(sqrt(eps()) * 1e5),
      T(sqrt(eps()) * 1e5),
      T(sqrt(eps(T)) * 1e0),
      T(sqrt(eps(T)) * 1e0),
      :classic,
    )
    E .= T(1.0e-2)
  end
  if regu.δ_min == zero(T) # gsp for gpmr
    regu.δ = zero(T)
  end
  invE = similar(E)
  invE .= one(T) ./ E

  ξ1 = similar(fd.c, id.nvar)
  ξ2 = similar(fd.c, id.ncon)
  Δy0 = similar(fd.c, id.ncon)
  ξ12 = similar(fd.c, id.nvar)

  KS = init_Ksolver(fd.uplo == :U ? fd.A : fd.A', ξ12, sp)

  return PreallocatedDataK1_1Structured(
    E,
    invE,
    ξ1,
    ξ2,
    Δy0,
    ξ12,
    sp.rhs_scale,
    regu,
    KS,
    T(sp.atol0),
    T(sp.rtol0),
    T(sp.atol_min),
    T(sp.rtol_min),
  )
end

function solver!(
  dd::AbstractVector{T},
  pad::PreallocatedDataK1_1Structured{T},
  dda::DescentDirectionAllocs{T},
  pt::Point{T},
  itd::IterData{T},
  fd::Abstract_QM_FloatData{T},
  id::QM_IntData,
  res::AbstractResiduals{T},
  cnts::Counters,
  T0::DataType,
  step::Symbol,
) where {T <: Real}
  pad.ξ1 .= @views step == :init ? fd.c : dd[1:(id.nvar)]
  pad.ξ2 .= @views (step == :init && all(dd[(id.nvar + 1):end] .== zero(T))) ? one(T) :
         dd[(id.nvar + 1):end]

  if pad.regu.δ == zero(T)
    pad.Δy0 .= zero(T)
  else
    pad.Δy0 .= .-pad.ξ2 ./ pad.regu.δ
  end
  if fd.uplo == :U
    mul!(pad.ξ12, fd.A, pad.Δy0)
  else
    mul!(pad.ξ12, fd.A', pad.Δy0)
  end
  pad.ξ12 .+= pad.ξ1
  if pad.rhs_scale
    ξ12Norm = kscale!(pad.ξ12)
  end
  ksolve!(
    pad.KS,
    fd.uplo == :U ? fd.A : fd.A',
    pad.ξ12,
    Diagonal(pad.invE),
    pad.regu.δ,
    verbose = 0,
    atol = pad.atol,
    rtol = pad.rtol,
  )
  if pad.rhs_scale
    kunscale!(pad.KS.x, ξ12Norm)
  end
  dd[(id.nvar + 1):end] .= pad.KS.x .- pad.Δy0
  if fd.uplo == :U
    @views mul!(pad.ξ1, fd.A, dd[(id.nvar + 1):end], one(T), -one(T))
  else
    @views mul!(pad.ξ1, fd.A', dd[(id.nvar + 1):end], one(T), -one(T))
  end
  dd[1:(id.nvar)] .= pad.ξ1 ./ pad.E
  update_kresiduals_history_K1struct!(
    res,
    fd.uplo == :U ? fd.A' : fd.A,
    pad.E,
    pad.invE, # tmp storage vector
    pad.regu.δ,
    pad.KS.x,
    pad.ξ12,
    :K1_1,
  )

  return 0
end

function update_pad!(
  pad::PreallocatedDataK1_1Structured{T},
  dda::DescentDirectionAllocs{T},
  pt::Point{T},
  itd::IterData{T},
  fd::Abstract_QM_FloatData{T},
  id::QM_IntData,
  res::AbstractResiduals{T},
  cnts::Counters,
  T0::DataType,
) where {T <: Real}
  if cnts.k != 0
    update_regu!(pad.regu)
  end

  if pad.atol > pad.atol_min
    pad.atol /= 10
  end
  if pad.rtol > pad.rtol_min
    pad.rtol /= 10
  end

  pad.E .= pad.regu.ρ
  pad.E[id.ilow] .+= pt.s_l ./ itd.x_m_lvar
  pad.E[id.iupp] .+= pt.s_u ./ itd.uvar_m_x
  pad.invE .= one(T) ./ pad.E

  return 0
end


text:
abstract LinearSolver{T}

factorize!(s::LinearSolver) = error("factorize! not implemented for $(typeof(s))")
solve!(s::LinearSolver) = error("solve! not implemented for $(typeof(s))")


text:
using Test
using TensorNetworkTensors
using TensorOperations: @tensor, tensorcopy
using LinearAlgebra: diag
@testset "TensorNetworkTensors" begin
    include("tensors.jl")
end


text:
function predict(P,x)
    w1, b1, w2 = P
    h = sigm(w1 * x .+ b1)
    return w2 * h
end

#when we have analytical solution available

# function loss_reg(w,x,ygold)
#     sumabs2(ygold - predict(w,x)) / size(ygold,2)
# end

function check_Phi_gradient(P,timepoints,hl_width)
    w,b,v = P
    u0 = 1/2
    t0 = timepoints[1]
    Phi(t) = u0 + (t-t0)*predict(P,t)
    dN_dt(t) = sum([v[i]*w[i]*sig_der(w[i]*t .+ b[i]) for i = 1:hl_width])
    dPhi_dt(t) = predict(P,t)+(t-t0)*dN_dt(t)
    for t in timepoints
        println(:ForwardDiff,gradient(x->Phi(x),t))
        println(:analytical_grad,dPhi_dt(t))
    end
end

function loss_trial(P,timepoints,f,phi,hl_width)
    w,b,v = P
    t0 = timepoints[1]
    dN_dt(t) = sum([v[i]*w[i]*sig_der(w[i]*t .+ b[i]) for i = 1:hl_width])
    #dPhi_dt(t) = predict(P,t)+(t-t0)*dN_dt(t)
    sumabs2([gradient(x->phi(P,x),t) - f(t,phi(P,t)) for t in timepoints][1])

end

lossgradient = grad(loss_trial)

function train(P, prms, timepoints, f, phi, hl_width; maxiters =100)
    #print(size(P),size(g),size(prms))
    for iter=1:maxiters
        #println("epoch no.",epoch)
        for x in timepoints
            g = lossgradient(P,timepoints,f,phi,hl_width)
          #print(size(P),size(g),size(prms))
          update!(P, g, prms)
          #println(P[1][1],P[2][1],P[3][1],P[4][1])
          #println(loss_trial(P,x))
        end

    end
    return P
end

function test(P,timepoints,f,phi,hl_width)
    sumloss = numloss = 0
    for t in timepoints
        sumloss += loss_trial(P,timepoints,f,phi,hl_width)
        numloss += 1
    end
    return sumloss/numloss
end

function init_params(ftype,hl_width;atype=KnetArray{Float32})
    #P = Vector{Vector{Float32}}(4)
    P = Array{Any}(3)
    P[1] = randn(ftype,hl_width,1)
    P[2] = zeros(ftype,hl_width,1)
    P[3] = randn(ftype,1,hl_width)
    #P[4] = zeros(Float32,1,1)
    #P = map(x -> convert(atype, x), P)

    return P
end

function generate_data(low, high, dt; atype=KnetArray{Float32})
    #data = Any[]
    num_points = 1/dt
    x = linspace(low,high,num_points+1)
    #x_ = convert(atype,x)
    # for i=1:num_points
    #     push!(data, (x_[i],y_[i]))
    # end
    return x
end


text:
Base.include(Main, joinpath(@__DIR__, "startup.jl"))

# From Base._start:
@eval Main import Base.MainInclude: eval, include
try
    Base.exec_options(Base.JLOptions())
catch err
    Base.invokelatest(Base.display_error, err, catch_backtrace())
    exit(1)
end
if Base.is_interactive && Base.have_color
    print(Base.color_normal)
end

max_stars_repo_name: tkf/julia-venv

max_issues_repo_name: tkf/julia-venv

max_forks_repo_name: tkf/julia-venv


text:
using FileIO
using Printf

const HORIZONTAL = 1
const VERTICAL   = 2

function thinning_zs_iteration(bim,iter)
    h,w = size(bim)
    marker = zeros(Bool,h,w)
    for i=2:h-1
        for j=2:w-1
            p2 = bim[i-1,j]   & 1
            p3 = bim[i-1,j+1] & 1
            p4 = bim[i,j+1]   & 1
            p5 = bim[i+1,j+1] & 1
            p6 = bim[i+1,j]   & 1
            p7 = bim[i+1,j-1] & 1
            p8 = bim[i,j-1]   & 1
            p9 = bim[i-1,j-1] & 1
            a = (p2 == 0 && p3 == 1)+ (p3 == 0 && p4 == 1)+
                (p4 == 0 && p5 == 1)+ (p5 == 0 && p6 == 1)+
                (p6 == 0 && p7 == 1)+ (p7 == 0 && p8 == 1)+
                (p8 == 0 && p9 == 1)+ (p9 == 0 && p2 == 1)
            b = p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9
            m1 = iter==0 ? (p2*p4*p6) : (p2*p4*p8);
            m2 = iter==0 ? (p4*p6*p8) : (p2*p6*p8);
            if a == 1 && (b >= 2 && b <= 6) && m1 == 0 && m2 == 0
                marker[i,j] = true;
            end
        end
    end
    nim = bim .& (.!marker)
    return nim
end

function thinning_zs(bim)
    while true
        oim = copy(bim)
        bim = thinning_zs_iteration(bim,0)
        bim = thinning_zs_iteration(bim,1)
        dif = sum(bim .⊻ oim)
        if dif == 0
            break
        end
    end
    return bim
end

function merge_impl!(c0,c1,i,sx,isv,mode)
    b0 = Bool(mode >> 1 & 1)
    b1 = Bool(mode >> 0 & 1)
    mj = -1
    md = 4
    p1 = c1[i][b1 ? 1 : end]

    if abs(p1[isv ? 2 : 1] - sx) > 0
        return false
    end
    # find the best match
    for j=1:length(c0)
        p0 = c0[j][b0 ? 1 : end]
        if abs(p0[isv ? 2 : 1] - sx) > 1
            continue
        end
        d = abs(p0[isv ? 1 : 2] - p1[isv ? 1 : 2])
        if d < md
            mj = j
            md = d
        end
    end
    if mj != -1
        j = mj
        if b0 && b1
            splice!(c0[j],1:0,reverse(c1[i]))
        elseif !b0 && b1
            append!(c0[j],c1[i])
        elseif b0 && !b1
            splice!(c0[j],1:0,c1[i])
        else
            append!(c0[j],reverse(c1[i]))
        end
        deleteat!(c1,i)
        return true
    end
    return false
end

function translate!(c,sx,dr)
    x,y = sx-1,0
    if dr == VERTICAL
        y,x = sx-1,0
    end
    for i=1:length(c)
        for j=1:length(c[i])
            c[i][j][1]+=x
            c[i][j][2]+=y
        end
    end
end

function merge_frags!(c0,c1,sx,dr)
    translate!(c1,sx,dr)
    if length(c0) == 0
        append!(c0,c1)
        return
    end
    if length(c1) == 0
        return
    end
    l = length(c1)
    for i=l:-1:1
        if dr == HORIZONTAL
            if merge_impl!(c0,c1,i,sx,false,1) continue end
            if merge_impl!(c0,c1,i,sx,false,3) continue end
            if merge_impl!(c0,c1,i,sx,false,0) continue end
            if merge_impl!(c0,c1,i,sx,false,2) continue end
        elseif dr == VERTICAL
            if merge_impl!(c0,c1,i,sx,true ,1) continue end
            if merge_impl!(c0,c1,i,sx,true ,3) continue end
            if merge_impl!(c0,c1,i,sx,true ,0) continue end
            if merge_impl!(c0,c1,i,sx,true ,2) continue end
        end
    end
    append!(c0,c1)
    return
end

function chunk_to_frags(bim)
    h,w = size(bim)
    frags = []
    on = false
    li = -1
    lj = -1
    for k=1:h+h+w+w-4
        i = 0
        j = 0
        if k <= w
            i = 1; j = k
        elseif k <= w+h-1
            i = k-w+1; j = w
        elseif k <= w+h+w-2
            i = h; j = w-(k-w-h+1)
        else
            i = h-(k-w-h-w+2); j = 1
        end

        if bim[i,j] != 0
            if !on
                on = true
                push!(frags,[[j,i],[w÷2,h÷2]])
            end
        else
            if on # right side of stroke, average to get center of stroke
                on = false
                frags[end][1][1] = (frags[end][1][1]+lj)÷2
                frags[end][1][2] = (frags[end][1][2]+li)÷2
            end
        end
        li = i
        lj = j
    end
    if length(frags) == 2
        frags = [[frags[1][1],frags[2][1]]]
    elseif length(frags) > 2
        ms = 0
        mi = -1
        mj = -1
        # use convolution to find brightest blob
        for i=2:h-1
            for j=2:w-1
                s =
                   (bim[i-1,j-1]) + (bim[i-1,j]) +(bim[i-1,j+1])+
                   (bim[i,j-1]  ) +   (bim[i,j]) +    (bim[i,j+1])+
                   (bim[i+1,j-1]) + (bim[i+1,j]) +  (bim[i+1,j+1]);
                if s > ms
                    mi = i
                    mj = j
                    ms = s
                elseif s == ms && abs(j-(w÷2))+abs(i-(h÷2)) < abs(mj-(w÷2))+abs(mi-(h÷2))
                    mi = i
                    mj = j
                    ms = s
                end
            end
        end
        if mi != -1
            for i=1:length(frags)
                frags[i][2] = [mj,mi]
            end
        end
    end
    return frags
end

function trace_skeleton(bim,csize,maxiter)

    h,w = size(bim)
    if maxiter <= 0
        return []
    end
    if w <= csize && h <= csize
        return chunk_to_frags(bim)
    end
    ms = w+h
    mi = -1
    mj = -1
    if h > csize
        for i=4:h-3
            if bim[i,1]>0 || bim[i-1,1]>0 || bim[i,w]>0 || bim[i-1,w]>0
                continue
            end
            s = 0
            for j=1:w
                s+=bim[i,j]
                s+=bim[i-1,j]
            end
            if s < ms
                ms = s; mi = i
            elseif s == ms && abs(i-h÷2)<abs(mi-h÷2)
                # if there is a draw (very common), we want the seam to be near the middle
                # to balance the divide and conquer tree
                ms = s; mi = i
            end
        end
    end
    if w > csize
        for j=4:w-3
            if bim[1,j]>0 || bim[h,j]>0 || bim[1,j-1]>0 || bim[h,j-1]>0
                continue
            end
            s = 0
            for i=1:h
                s+=bim[i,j]
                s+=bim[i,j-1]
            end
            if s < ms
                ms = s;
                mi = -1
                mj = j
            elseif s == ms && abs(j-w÷2)<abs(mj-w÷2)
                ms = s
                mi = -1
                mj = j
            end
        end
    end
    dr = 0
    sx = 0
    L = []
    R = []

    if h > csize && mi != -1
        L = view(bim, 1:mi-1,:)
        R = view(bim, mi:h,  :)
        dr = VERTICAL;
        sx = mi
    elseif w > csize && mj != -1
        L = view(bim, :, 1:mj-1)
        R = view(bim, :, mj:w)
        dr = HORIZONTAL
        sx = mj
    end
    frags = []
    if dr != 0 && sum(L)>0
        frags = trace_skeleton(L,csize,maxiter-1)
    end
    if dr != 0 && sum(R)>0
        merge_frags!(frags,trace_skeleton(R,csize,maxiter-1),sx,dr)
    end
    if mi == -1 && mj == -1
        return chunk_to_frags(bim)
    end
    return frags
end

function polylines_to_svg(q,w,h)
    svg = @sprintf("<svg xmlns=\"http://www.w3.org/2000/svg\" width=\"%d\" height=\"%d\" fill=\"none\" stroke=\"black\" stroke-width=\"1\">",w,h)
    for i=1:length(q)
        svg *= "<path d=\""
        for j=1:length(q[i])
            svg *= @sprintf("%s%d,%d ",(j==1 ? "M" : "L"),q[i][j][1],q[i][j][2])
        end
        svg *= "\"/>"
    end
    svg *= "</svg>"
    return svg
end

img = load(ARGS[1])
bim = convert(Array{UInt8}, reshape([round(img[i,j].r) for j=1:size(img)[2] for i=1:size(img)[1]], size(img)))

print("<!-- ")
@time bim = thinning_zs(bim)
H,W = size(bim)
@time q = trace_skeleton(bim,10,999)
print(" -->")
print(polylines_to_svg(q,W,H))

# using ColorTypes
# img = convert(Array{ColorTypes.Gray{Bool}}, convert(Array{Bool},bim))
# save("out.jpg",img)

max_stars_repo_name: freder/skeleton-tracing

max_issues_repo_name: freder/skeleton-tracing

max_forks_repo_name: freder/skeleton-tracing


text:
ENV["GKSwstype"] = "100" #src

# # Generalized Linear Models
#
# Generalized Linear Models (GLMs) are a family of statistical models where the prediction
# is based on a function of a linear predictor $x_^T\beta$ where $x$ is a vector of features
# and $\beta$ is the parameter vector.
#
# One of the big assumptions of GLMs is that there is a strictly increasing *link function*
# $g$ that "links" the mean of the distribution $E(Y|X) = \mu$ to the linear predictor
# $x^T\beta$.  In other words, our prediction of $y$ given $x$ is $g^{-1}(x^T\beta)$.

# The two most common GLMs are linear and logistic regression.

# The GLM package
#
#

# ## Linear Regression

# ## Logistic Regression

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VirtualMachineScaleSetNetworkConfigurationProperties <: SwaggerModel
    primary::Any # spec type: Union{ Nothing, Bool } # spec name: primary
    enableAcceleratedNetworking::Any # spec type: Union{ Nothing, Bool } # spec name: enableAcceleratedNetworking
    networkSecurityGroup::Any # spec type: Union{ Nothing, SubResource } # spec name: networkSecurityGroup
    dnsSettings::Any # spec type: Union{ Nothing, VirtualMachineScaleSetNetworkConfigurationDnsSettings } # spec name: dnsSettings
    ipConfigurations::Any # spec type: Union{ Nothing, Vector{VirtualMachineScaleSetIPConfiguration} } # spec name: ipConfigurations
    enableIPForwarding::Any # spec type: Union{ Nothing, Bool } # spec name: enableIPForwarding

    function VirtualMachineScaleSetNetworkConfigurationProperties(;primary=nothing, enableAcceleratedNetworking=nothing, networkSecurityGroup=nothing, dnsSettings=nothing, ipConfigurations=nothing, enableIPForwarding=nothing)
        o = new()
        validate_property(VirtualMachineScaleSetNetworkConfigurationProperties, Symbol("primary"), primary)
        setfield!(o, Symbol("primary"), primary)
        validate_property(VirtualMachineScaleSetNetworkConfigurationProperties, Symbol("enableAcceleratedNetworking"), enableAcceleratedNetworking)
        setfield!(o, Symbol("enableAcceleratedNetworking"), enableAcceleratedNetworking)
        validate_property(VirtualMachineScaleSetNetworkConfigurationProperties, Symbol("networkSecurityGroup"), networkSecurityGroup)
        setfield!(o, Symbol("networkSecurityGroup"), networkSecurityGroup)
        validate_property(VirtualMachineScaleSetNetworkConfigurationProperties, Symbol("dnsSettings"), dnsSettings)
        setfield!(o, Symbol("dnsSettings"), dnsSettings)
        validate_property(VirtualMachineScaleSetNetworkConfigurationProperties, Symbol("ipConfigurations"), ipConfigurations)
        setfield!(o, Symbol("ipConfigurations"), ipConfigurations)
        validate_property(VirtualMachineScaleSetNetworkConfigurationProperties, Symbol("enableIPForwarding"), enableIPForwarding)
        setfield!(o, Symbol("enableIPForwarding"), enableIPForwarding)
        o
    end
end # type VirtualMachineScaleSetNetworkConfigurationProperties

const _property_map_VirtualMachineScaleSetNetworkConfigurationProperties = Dict{Symbol,Symbol}(Symbol("primary")=>Symbol("primary"), Symbol("enableAcceleratedNetworking")=>Symbol("enableAcceleratedNetworking"), Symbol("networkSecurityGroup")=>Symbol("networkSecurityGroup"), Symbol("dnsSettings")=>Symbol("dnsSettings"), Symbol("ipConfigurations")=>Symbol("ipConfigurations"), Symbol("enableIPForwarding")=>Symbol("enableIPForwarding"))
const _property_types_VirtualMachineScaleSetNetworkConfigurationProperties = Dict{Symbol,String}(Symbol("primary")=>"Bool", Symbol("enableAcceleratedNetworking")=>"Bool", Symbol("networkSecurityGroup")=>"SubResource", Symbol("dnsSettings")=>"VirtualMachineScaleSetNetworkConfigurationDnsSettings", Symbol("ipConfigurations")=>"Vector{VirtualMachineScaleSetIPConfiguration}", Symbol("enableIPForwarding")=>"Bool")
Base.propertynames(::Type{ VirtualMachineScaleSetNetworkConfigurationProperties }) = collect(keys(_property_map_VirtualMachineScaleSetNetworkConfigurationProperties))
Swagger.property_type(::Type{ VirtualMachineScaleSetNetworkConfigurationProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VirtualMachineScaleSetNetworkConfigurationProperties[name]))}
Swagger.field_name(::Type{ VirtualMachineScaleSetNetworkConfigurationProperties }, property_name::Symbol) =  _property_map_VirtualMachineScaleSetNetworkConfigurationProperties[property_name]

function check_required(o::VirtualMachineScaleSetNetworkConfigurationProperties)
    (getproperty(o, Symbol("ipConfigurations")) === nothing) && (return false)
    true
end

function validate_property(::Type{ VirtualMachineScaleSetNetworkConfigurationProperties }, name::Symbol, val)
end


text:
function transpose_strings(input)

end


text:
@testset "MembershipTest" begin
    @polyvar x y
    f = x^2 + y^2 + 1

    F = AmoebaFiber2D(f)

    result = membershiptest(F, (0.0, 0.0), PolynomialAmoebas.TorusStartValueGenerator{2}())
    @test result.successfull

    r2 = membershiptest(F, (0.1, 0.1), PolynomialAmoebas.TorusStartValueGenerator{2}(), (result.solution,), MembershipTestOptions())
    @test r2.successfull
    @test r2.tries == 1
end


text:
using Distributed
using Distributions

nb_draws = 100000

function inside_circle(x::Float64, y::Float64)
    output = 0
    if x^2 + y^2 <= 1
	output = 1
    end
    return output
end

function pi_serial(nbPoints::Int64 = 128 * 1000; d=Uniform(-1.0,1.0))

   #draw NbPoints from within the square centered in 0
   #with side length equal to 2
   xDraws = rand(d, nbPoints)
   yDraws = rand(d, nbPoints)
   sumInCircle = 0

   for (xValue, yValue) in zip(xDraws, yDraws)
	sumInCircle+=inside_circle(xValue, yValue)
   end

   return 4*sumInCircle/nbPoints

end

pi_approx = pi_serial(nb_draws);

println("Approximation for Pi $(pi_approx) with $(nb_draws) random draws")
println("True value is $(pi)")

max_stars_repo_name: JulienPascal/MyFirstDockerApp

max_issues_repo_name: JulienPascal/MyFirstDockerApp

max_forks_repo_name: JulienPascal/MyFirstDockerApp