text:
"""
`module  JAC.Atomic`  ... a submodel of JAC that contains all methods to set-up and process (simple) atomic and SCF computations as well as
                          atomic cascade computations; it is using JAC, JAC.Radial, JAC.ManyElectron, JAC.Nuclear, JAC.Einstein,
                          JAC.Hfs, JAC.IsotopeShift, JAC.PlasmaShift, JAC.LandeZeeman, JAC.PhotoEmission, JAC.PhotoIonization, JAC.AutoIonization.
"""
module Atomic

    using Interact
    using JAC.Basics, JAC.Radial, JAC.ManyElectron, JAC.Nuclear
    using JAC.Einstein, JAC.Hfs, JAC.IsotopeShift, JAC.PlasmaShift, JAC.LandeZeeman, JAC.AlphaVariation, JAC.FormFactor, JAC.DecayYield,
          JAC.GreenFunction,  JAC.MultipolePolarizibility,
          JAC.PhotoEmission, JAC.PhotoIonization, JAC.PhotoExcitation, JAC.PhotoRecombination,
          JAC.AutoIonization, JAC.Dielectronic, JAC.ImpactExcitation, JAC.CoulombExcitation, JAC.CoulombIonization,
          JAC.PhotoExcitationFluores, JAC.PhotoExcitationAutoion, JAC.RayleighCompton, JAC.MultiPhotonDeExcitation,
          JAC.PhotoIonizationFluores, JAC.PhotoIonizationAutoion, JAC.ImpactExcitationAutoion, JAC.RadiativeAuger,
          JAC.MultiPhotonIonization, JAC.MultiPhotonDoubleIon, JAC.InternalConversion
          ##x , JAC.ElectricDipoleMoment

    """
    `struct  Atomic.CasStep`  ... a type for defining an individual step of a (relativistic) complete active space computation for a specified
                                  set of levels. An instance of this struct provides all information to generate the atomic basis and to perform
                                  the associated SCF and multiplet computations for a selected No. of levels.

        + name             ::String                    ... to assign a name to the given SCF step.
        + NoElectrons      ::Int64                     ... Number of electrons.
        + refConfigs       ::Array{Configuration,1}    ... List of references configurations, at least 1.
        + Sfrom            ::Array{Shell,1}            ... Single-excitations from shells   [sh_1, sh_2, ...]
        + Sto              ::Array{Shell,1}            ... Single-excitations to shells  [sh_1, sh_2, ...]
        + Dfrom            ::Array{Shell,1}            ... Double-excitations from shells   [sh_1, sh_2, ...]
        + Dto              ::Array{Shell,1}            ... Double-excitations to shells  [sh_1, sh_2, ...]
        + Tfrom            ::Array{Shell,1}            ... Triple-excitations from shells   [sh_1, sh_2, ...]
        + Tto              ::Array{Shell,1}            ... Triple-excitations to shells  [sh_1, sh_2, ...]
        + Qfrom            ::Array{Shell,1}            ... Quadrupole-excitations from shells   [sh_1, sh_2, ...]
        + Qto              ::Array{Shell,1}            ... Quadrupole-excitations to shells  [sh_1, sh_2, ...]
        + restrictions     ::Array{String,1}           ... List of Strings to define 'restrictions' to applied to the CSF basis.
        + NoIterations     ::Int64                     ... Number of SCf iterations to be applied in this step of computations.
        + frozenSubshells  ::Array{Subshell,1}         ... List of subshells that are kept 'frozen' in this step.
        + failureHandling  ::Array{String,1}           ... List of Strings to define subsequent steps in case of failure.
    """
    struct  CasStep
        name               ::String
        NoElectrons        ::Int64
        refConfigs         ::Array{Configuration,1}
        Sfrom              ::Array{Shell,1}
        Sto                ::Array{Shell,1}
        Dfrom              ::Array{Shell,1}
        Dto                ::Array{Shell,1}
        Tfrom              ::Array{Shell,1}
        Tto                ::Array{Shell,1}
        Qfrom              ::Array{Shell,1}
        Qto                ::Array{Shell,1}
        restrictions       ::Array{String,1}
        NoIterations       ::Int64
        frozenShells       ::Array{Subshell,1}
        failureHandling    ::Array{String,1}
    end

    """
    `JAC.Atomic.CasStep(name::String, NoElectrons::Int64)`  ... constructor for an 'initial' instance of a variable::Atomic.CasStep with just
                                                                a given name and No. of Electrons.
    """
    function CasStep(name::String, NoElectrons::Int64)
        CasStep(name, NoElectrons, Configuration[], Shell[], Shell[], Shell[], Shell[], Shell[], Shell[], Shell[], Shell[],
                String[], 0, Subshell[], String[])
    end

    """
    `JAC.Atomic.CasStep("interactive"; refine::Atomic.CasStep=..)`  ... constructor to generate a new instance of Atomic.CasStep
                                       interactively by replying to some detailed dialog. This constructor enables one to add, modify or
                                       refine an individual step in the complete active space computations. **Not yet implemented !**
    """
    function CasStep(sa::String; refine::Atomic.CasStep = Atomic.CasStep())
        sa != "interactive"   &&    error("Unsupported keystring = $sa.")

        #= print("Enter the nuclear charge Z:  ");    Zx = Base.parse(Float64, readline(STDIN) )

        yes = yesno("The ...; modify ? [N|y]")
        while  yes
            print("Enter a charge distribution model {point, Fermi, uniform}:  ");    modelx = strip(readline(STDIN))
            !(modelx in [""])   &&   println("Unsupported charge distribution model $modelx ... redo:")
        end =#
    end

    # `Base.show(io::IO, step::Atomic.CasStep)`  ... prepares a proper printout of the (individual step of computations) step::Atomic.CasStep.
    function Base.show(io::IO, step::Atomic.CasStep)
        sa = Base.string(step);                print(io, sa, "\n")
        sa = "Reference configurations: ";     print(io, sa, "\n")
        if  length(step.refConfigs) > 0
           for  conf  in step.refConfigs          print(io, conf )  end;   print(io, "\n")
        end
        if  length(step.Sfrom) > 0
           sa = "Singles from:              { ";      for  sh in step.Sfrom    sa = sa * string(sh) * ", "  end
           sa = sa[1:end-2] *  " }   ... to { ";      for  sh in step.Sto      sa = sa * string(sh) * ", "  end;
           sa = sa[1:end-2] *  " }";            print(io, sa, "\n")
        end
        if  length(step.Dfrom) > 0
           sa = "Doubles from:              { ";      for  sh in step.Dfrom    sa = sa * string(sh) * ", "  end
           sa = sa[1:end-2] *  " }   ... to { ";      for  sh in step.Dto      sa = sa * string(sh) * ", "  end;
           sa = sa[1:end-2] *  " }";            print(io, sa, "\n")
        end
        if  length(step.Tfrom) > 0
           sa = "Triples from:              { ";      for  sh in step.Tfrom    sa = sa * string(sh) * ", "  end
           sa = sa[1:end-2] *  " }   ... to { ";      for  sh in step.Tto      sa = sa * string(sh) * ", "  end;
           sa = sa[1:end-2] *  " }";            print(io, sa, "\n")
        end
        if  length(step.Qfrom) > 0
           sa = "Quadruples from:           { ";      for  sh in step.Qfrom    sa = sa * string(sh) * ", "  end
           sa = sa[1:end-2] *  " }   ... to { ";      for  sh in step.Qto      sa = sa * string(sh) * ", "  end;
           sa = sa[1:end-2] *  " }";            print(io, sa, "\n")
        end
    end

    # `Base.string(step::Atomic.CasStep)`  ... provides a String notation for the variable step::Atomic.CasStep.
    function Base.string(step::Atomic.CasStep)
        sa = "Step: Computational model for $(step.name) with $(step.NoElectrons) electrons and with (step.NoCoreShells) "
        sa = sa * "core shells:"
        return( sa )
    end

    """
    `struct  Atomic.CasComputation`  ... a struct for defining a series of complete active space computations for a specified set of levels.
                                         This type provides all information to generate the atomic basis and to perform the corresponding
                                         SCF and multiplet computations for a selected No. of levels.

        + name             ::String                     ... to assign a name to the given model.
        + previousStep     ::Int64                      ... 0, if no previous steps was yet done in this series.
        + steps            ::Array{Atomic.CasStep,1}    ... List of SCF steps that are to be done in this model computation.
    """
    struct  CasComputation
        name               ::String
        previousStep       ::Int64
        steps              ::Array{Atomic.CasStep,1}
    end

    """
    `JAC.Atomic.CasComputation()`  ... constructor for an 'empty' instance of the a variable::Atomic.CasComputation
    """
    function CasComputation()
        CasComputation("", 0, Array{CasStep,1}[])
    end

    """
    `JAC.Atomic.CasComputation("interactive"; refine::Atomic.CasComputation=..)`  ... constructor to generate a new instance of
                               Atomic.CasComputation interactively by replying to some detailed dialog. This constructor enables one to add,
                               modify or refine the model for a particular complete active space computations. **Not yet implemented !**
    """
    function CasComputation(sa::String; refine::Atomic.CasComputation = Atomic.CasComputation())
        sa != "interactive"   &&    error("Unsupported keystring = $sa.")

        #= print("Enter the nuclear charge Z:  ");    Zx = Base.parse(Float64, readline(STDIN) )

        yes = yesno("The ...; modify ? [N|y]")
        while  yes
            print("Enter a charge distribution model {point, Fermi, uniform}:  ");    modelx = strip(readline(STDIN))
            !(modelx in [""])   &&   println("Unsupported charge distribution model $modelx ... redo:")
        end =#
    end

    # `Base.show(io::IO, model::CasComputation)`  ... prepares a proper printout of the (individual step of computations) model::CasComputation.
    function Base.show(io::IO, model::CasComputation)
        sa = Base.string(model)
        print(io, sa)
    end

    # `Base.string(model::CasComputation)`  ... provides a String notation for the variable model::CasComputation.
    function Base.string(model::CasComputation)
        sa = "Model: $(model.name) includes $(length(model.steps)) steps and has been calculated up to the step $(model.previousStep)."
        return( sa )
    end

    """
    `struct  Atomic.CasSettings`  ... a struct for defining the settings for complete active space computations

        + generateScf       ::Bool               ... True, if a SCF need to be generated, and false if just the start orbitals should
                                                     be applied.
        + startOrbitals     ::String             ... Specify how the start orbitals are obtained ["fromNRorbitals", "fromGrasp", "hydrogenic"].
        + rwfFilename       ::String             ... Filename of orbitals, if taken from Grasp.
        + levels            ::Array{Int64,1}     ... Levels on which the optimization need to be carried out.
        + includeBreit      ::Bool               ... True, if the Breit interaction is included into the SCF generation.
        + maxIterations     ::Int64              ... maximum number of SCF iterations
        + scfAccuracy       ::Float64            ... convergence criterion for the SCF field.
        + iterationSequence ::Array{Subshell,1}  ... Sequence of subshells to be optimized.
    """
    struct  CasSettings
        generateScf         ::Bool
        startOrbitals       ::String
        rwfFilename         ::String
        levels              ::Array{Int64,1}
        includeBreit        ::Bool
        maxIterations       ::Int64
        scfAccuracy         ::Float64
        iterationSequence   ::Array{Subshell,1}
    end

    """
    `JAC.Atomic.CasSettings()`  ... constructor for setting the default values.
    """
    function CasSettings()
        CasSettings(false, "fromNRorbitals", "", Int64[1], false, 24, 1.0e-8, Subshell[] )
    end

    # `Base.show(io::IO, settings::Atomic.CasSettings)`  ... prepares a proper printout of the settings::Atomic.CasSettings.
    function Base.show(io::IO, settings::Atomic.CasSettings)
        println(io, "generateScf:        $(settings.generateScf)  ")
        println(io, "startOrbitals:      $(settings.startOrbitals)  ")
        println(io, "rwfFilename:        $(settings.rwfFilename)  ")
        println(io, "levels:             $(settings.levels)  ")
        println(io, "includeBreit:       $(settings.includeBreit)  ")
        println(io, "maxIterations:      $(settings.maxIterations)  ")
        println(io, "scfAccuracy:        $(settings.scfAccuracy)  ")
        println(io, "iterationSequence:  $(settings.iterationSequence)  ")
    end

    # `Base.string(settings::Atomic.CasSettings)`  ... provides a String notation for the variable settings::Atomic.CasSettings.
    function Base.string(settings::Atomic.CasSettings)
        error("Not yet implemented.")
        sa = "Cas settings: maximum No. of iterations = $(settings.maxIterations), accuracy = (settings.scfAccuracy)"
        return( sa )
    end

    """
    `struct  Computation`  ... defines a type for defining  (the model of simple) atomic computation of a single multiplet,
                               including the SCF and CI as well as level properties and transition property calculations.

        + name                           ::String                        ... A name associated to the computation.
        + nuclearModel                   ::Nuclear.Model          ... Model, charge and parameters of the nucleus.
        + grid                           ::Radial.Grid                   ... The radial grid to be used for the computation.
        ##x + calcLevelProperties            ::Bool                          ... True, if level structures and properties are to be calculated
        + properties                     ::Array{AtomicLevelProperty,1}  ... List of atomic properties to be calculated.
        + configs                        ::Array{Configuration,1}        ... A list of non-relativistic configurations.
        + asfSettings                    ::AsfSettings                   ... Provides the settings for the SCF process and for the CI and QED calculations.
        + initialConfigs                 ::Array{Configuration,1}        ... A list of initial-state configurations for some transition
                                                                             property calculation, such as radiative transition, Auger, etc.
        + initialAsfSettings             ::AsfSettings                   ... Provides the SCF and CI settings for the initial-state multiplet.
        + intermediateConfigs            ::Array{Configuration,1}        ... A list of initial-state configurations.
        + intermediateAsfSettings        ::AsfSettings                   ... Provides the SCF settings for the intermediate-state multiplet.
        + finalConfigs                   ::Array{Configuration,1}        ... A list of final-state configurations.
        + finalAsfSettings               ::AsfSettings                   ... Provides the SCF and CI settings for the final-state multiplet.
        + alphaSettings                  ::AlphaVariation.Settings       ... Settings for alpha-variation parameter calculations.
        + einsteinSettings               ::Einstein.Settings             ... Settings for Einstein coefficient calculations.
        + formSettings                   ::FormFactor.Settings           ... Settings for atomic form factor calculations.
        + greenSettings                  ::GreenFunction.Settings        ... Settings for approximate Green function calculations.
        + hfsSettings                    ::Hfs.Settings                  ... Settings for hyperfine parameter calculations.
        + isotopeSettings                ::IsotopeShift.Settings         ... Settings for isotope shift parameter calculations.
        + plasmaSettings                 ::PlasmaShift.Settings          ... Settings for plasma-shift calculations.
        + polaritySettings               ::MultipolePolarizibility.Settings   ... Settings for plasma-shift calculations.
        + yieldSettings                  ::DecayYield.Settings           ... Settings for fluoresence and Auger yield calculations.
        + zeemanSettings                 ::LandeZeeman.Settings          ... Settings for Lande-Zeeman coefficient calculations.
        + process                        ::JAC.AtomicProcess             ... An (additional) process for which the properties are to be evaluated
                                                                             for the given initial- and final-state configurations.
        + processSettings                ::Union{JAC.PhotoEmission.Settings, JAC.AutoIonization.Settings, JAC.PlasmaShift.AugerSettings,
                                                 JAC.PhotoIonization.Settings, JAC.PlasmaShift.PhotoSettings,
                                                 JAC.PhotoExcitation.Settings, JAC.PhotoExcitationAutoion.Settings, JAC.PhotoRecombination.Settings,
                                                 JAC.ImpactExcitation.Settings, JAC.Dielectronic.Settings, RadiativeAuger.Settings,
                                                 JAC.PairAnnihilation1Photon.Settings, JAC.ImpactExcitationAutoion.Settings,
                                                 JAC.MultiPhotonDeExcitation.Settings, JAC.CoulombExcitation.Settings,
                                                 JAC.CoulombIonization.Settings} ... Provides the settings for the selected process.
    """
    struct  Computation
        name                           ::String
        nuclearModel                   ::Nuclear.Model
        grid                           ::Radial.Grid
        ##x calcLevelProperties            ::Bool
        properties                     ::Array{AtomicLevelProperty,1}
        configs                        ::Array{Configuration,1}
        asfSettings                    ::AsfSettings
        initialConfigs                 ::Array{Configuration,1}
        initialAsfSettings             ::AsfSettings
        intermediateConfigs            ::Array{Configuration,1}
        intermediateAsfSettings        ::AsfSettings
        finalConfigs                   ::Array{Configuration,1}
        finalAsfSettings               ::AsfSettings
        alphaSettings                  ::AlphaVariation.Settings
        einsteinSettings               ::Einstein.Settings
        formSettings                   ::FormFactor.Settings
        greenSettings                  ::GreenFunction.Settings
        hfsSettings                    ::Hfs.Settings
        isotopeSettings                ::IsotopeShift.Settings
        plasmaSettings                 ::PlasmaShift.Settings
        polaritySettings               ::MultipolePolarizibility.Settings
        yieldSettings                  ::DecayYield.Settings
        zeemanSettings                 ::LandeZeeman.Settings
        process                        ::AtomicProcess
        processSettings                ::Union{PhotoEmission.Settings, AutoIonization.Settings, PlasmaShift.AugerSettings,
                                               PhotoIonization.Settings, PlasmaShift.PhotoSettings,
                                               PhotoRecombination.Settings, Dielectronic.Settings, ImpactExcitation.Settings,
                                               CoulombExcitation.Settings, CoulombIonization.Settings,
                                               PhotoExcitation.Settings, PhotoExcitationFluores.Settings,
                                               PhotoExcitationAutoion.Settings, RayleighCompton.Settings,
                                               MultiPhotonDeExcitation.Settings, PhotoIonizationFluores.Settings,
                                               PhotoIonizationAutoion.Settings, ImpactExcitationAutoion.Settings,
                                               RadiativeAuger.Settings, MultiPhotonIonization.Settings,
                                               MultiPhotonDoubleIon.Settings, InternalConversion.Settings} #= ,
                                               #
                                               PairAnnihilation1Photon.Settings } =#
    end

    """
    `JAC.Atomic.Computation()`  ... constructor for an 'empty' instance::Atomic.Computation.
    """
    function Computation()
        Computation("", Nuclear.Model(1.), Radial.Grid(), AtomicLevelProperty[],
                    Configuration[], AsfSettings(),
                    Configuration[], AsfSettings(),
                    Configuration[], AsfSettings(),
                    Configuration[], AsfSettings(),
                    AlphaVariation.Settings(), Einstein.Settings(), FormFactor.Settings(), GreenFunction.Settings(),
                    Hfs.Settings(),
                    IsotopeShift.Settings(), PlasmaShift.Settings(), MultipolePolarizibility.Settings(),
                    DecayYield.Settings(),  LandeZeeman.Settings(),
                    ##x NoAmplitude, ElectricDipoleMoment.Settings(),
                    NoProcess, PhotoEmission.Settings() )
    end

    """
    `JAC.Atomic.Computation(sa::String, nm::Nuclear.Model; grid::Radial.Grid = JAC.Radial.Grid('grid: exponential'), ##x calc::Bool = false,
         properties::Array{AtomicLevelProperty,1} = AtomicLevelProperty[], configs::Array{Configuration,1} = Configuration[],
         asfSettings::AsfSettings = AsfSettings(),
         initialConfigs::Array{Configuration,1} = Configuration[], initialAsfSettings::AsfSettings = AsfSettings(),
         intermediateConfigs::Array{Configuration,1} = Configuration[], intermediateAsfSettings::AsfSettings = AsfSettings(),
         finalConfigs::Array{Configuration,1} = Configuration[], finalAsfSettings::AsfSettings = AsfSettings(),
         alphaSettings::AlphaVariation.Settings = AlphaVariation.Settings(), einsteinSettings::Einstein.Settings = Einstein.Settings(),
         formSettings::FormFactor.Settings = FormFactor.Settings(), greenSettings::GreenFunction.Settings = GreenFunction.Settings(),
         hfsSettings::Hfs.Settings = Hfs.Settings(),
         isotopeSettings::IsotopeShift.Settings = IsotopeShift.Settings(), plasmaSettings::PlasmaShift.Settings = PlasmaShift.Settings(),
         polaritySettings::MultipolePolarizibility.Settings = MultipolePolarizibility.Settings(),
         yieldSettings::DecayYield.Settings = DecayYield.Settings(), zeemanSettings::LandeZeeman.Settings = LandeZeeman.Settings(),
         ##x amplitude::JAC.AtomicAmplitude = JAC.NoAmplitude, amplitudeSettings::Any = true,
         process::JAC.AtomicProcess = JAC.NoProcess, processSettings::Any = true)`   ... constructor for an instance::Atomic.Computation for
         which all requested details are given by proper keyword specification. A few internal checks are  made.
    """
    function Computation(sa::String, nm::Nuclear.Model; grid::Radial.Grid = Radial.Grid("grid: exponential"), ##x calc::Bool = false,
                         properties::Array{AtomicLevelProperty,1} = [Basics.NoProperty], configs::Array{Configuration,1} = Configuration[],
                         asfSettings::AsfSettings = AsfSettings(),
                         initialConfigs::Array{Configuration,1} = Configuration[], initialAsfSettings::AsfSettings = AsfSettings(),
                         intermediateConfigs::Array{Configuration,1} = Configuration[], intermediateAsfSettings::AsfSettings = AsfSettings(),
                         finalConfigs::Array{Configuration,1} = Configuration[], finalAsfSettings::AsfSettings = AsfSettings(),
                         alphaSettings::AlphaVariation.Settings = AlphaVariation.Settings(),
                         einsteinSettings::Einstein.Settings = Einstein.Settings(),
                         formSettings::FormFactor.Settings = FormFactor.Settings(),
                         greenSettings::GreenFunction.Settings = GreenFunction.Settings(),
                         hfsSettings::Hfs.Settings = Hfs.Settings(),
                         isotopeSettings::IsotopeShift.Settings = IsotopeShift.Settings(),
                         plasmaSettings::PlasmaShift.Settings = PlasmaShift.Settings(),
                         polaritySettings::MultipolePolarizibility.Settings = MultipolePolarizibility.Settings(),
                         yieldSettings::DecayYield.Settings = DecayYield.Settings(),
                         zeemanSettings::LandeZeeman.Settings = LandeZeeman.Settings(),
                         process::Basics.AtomicProcess = Basics.NoProcess, processSettings::Any = true)
        if      process == Basics.NoProcess         && typeof(processSettings) == Bool && processSettings             procSettings = PhotoEmission.Settings()
        elseif  process == Basics.Auger             && typeof(processSettings) == Bool && processSettings             procSettings = AutoIonization.Settings()
        elseif  process == Basics.Auger             && typeof(processSettings) == AutoIonization.Settings             procSettings = processSettings
        elseif  process == Basics.AugerInPlasma     && typeof(processSettings) == Bool && processSettings             procSettings = PlasmaShift.AugerSettings()
        elseif  process == Basics.AugerInPlasma     && typeof(processSettings) == PlasmaShift.AugerSettings           procSettings = processSettings
        elseif  process == Basics.Radiative         && typeof(processSettings) == Bool && processSettings             procSettings = PhotoEmission.Settings()
        elseif  process == Basics.Radiative         && typeof(processSettings) == PhotoEmission.Settings              procSettings = processSettings
        elseif  process == Basics.PhotoExc          && typeof(processSettings) == Bool && processSettings             procSettings = PhotoExcitation.Settings()
        elseif  process == Basics.PhotoExc          && typeof(processSettings) == PhotoExcitation.Settings            procSettings = processSettings
        elseif  process == Basics.Photo             && typeof(processSettings) == Bool && processSettings             procSettings = PhotoIonization.Settings()
        elseif  process == Basics.Photo             && typeof(processSettings) == PhotoIonization.Settings            procSettings = processSettings
        elseif  process == Basics.PhotoInPlasma     && typeof(processSettings) == Bool && processSettings             procSettings = PlasmaShift.PhotoSettings()
        elseif  process == Basics.PhotoInPlasma     && typeof(processSettings) == PlasmaShift.PhotoSettings           procSettings = processSettings
        elseif  process == Basics.Rec               && typeof(processSettings) == Bool && processSettings             procSettings = PhotoRecombination.Settings()
        elseif  process == Basics.Rec               && typeof(processSettings) == PhotoRecombination.Settings         procSettings = processSettings
        elseif  process == Basics.Dierec            && typeof(processSettings) == Bool && processSettings             procSettings = Dielectronic.Settings()
        elseif  process == Basics.Dierec            && typeof(processSettings) == Dielectronic.Settings               procSettings = processSettings
        elseif  process == Basics.PhotoExcFluor     && typeof(processSettings) == Bool && processSettings             procSettings = PhotoExcitationFluores.Settings()
        elseif  process == Basics.PhotoExcFluor     && typeof(processSettings) == PhotoExcitationFluores.Settings     procSettings = processSettings
        elseif  process == Basics.PhotoExcAuto      && typeof(processSettings) == Bool && processSettings             procSettings = PhotoExcitationAutoion.Settings()
        elseif  process == Basics.PhotoExcAuto      && typeof(processSettings) == PhotoExcitationAutoion.Settings     procSettings = processSettings
        elseif  process == Basics.PhotoIonFluor     && typeof(processSettings) == Bool && processSettings             procSettings = PhotoIonizationFluores.Settings()
        elseif  process == Basics.PhotoIonFluor     && typeof(processSettings) == PhotoIonizationFluores.Settings     procSettings = processSettings
        elseif  process == Basics.PhotoIonAuto      && typeof(processSettings) == Bool && processSettings             procSettings = PhotoIonizationAutoion.Settings()
        elseif  process == Basics.PhotoIonAuto      && typeof(processSettings) == PhotoIonizationAutoion.Settings     procSettings = processSettings
        elseif  process == Basics.MultiPhotonDE     && typeof(processSettings) == Bool && processSettings             procSettings = MultiPhotonDeExcitation.Settings()
        elseif  process == Basics.MultiPhotonDE     && typeof(processSettings) == MultiPhotonDeExcitation.Settings    procSettings = processSettings
        elseif  process == Basics.Compton           && typeof(processSettings) == Bool && processSettings             procSettings = RayleighCompton.Settings()
        elseif  process == Basics.Compton           && typeof(processSettings) == RayleighCompton.Settings            procSettings = processSettings
        elseif  process == Basics.Eimex             && typeof(processSettings) == Bool && processSettings             procSettings = ImpactExcitation.Settings()
        elseif  process == Basics.Eimex             && typeof(processSettings) == ImpactExcitation.Settings           procSettings = processSettings
        elseif  process == Basics.ImpactExcAuto     && typeof(processSettings) == Bool && processSettings             procSettings = ImpactExcitationAutoion.Settings()
        elseif  process == Basics.ImpactExcAuto     && typeof(processSettings) == ImpactExcitationAutoion.Settings    procSettings = processSettings
        elseif  process == Basics.RAuger            && typeof(processSettings) == Bool && processSettings             procSettings = RadiativeAuger.Settings()
        elseif  process == Basics.RAuger            && typeof(processSettings) == RadiativeAuger.Settings             procSettings = processSettings
        elseif  process == Basics.MultiPI           && typeof(processSettings) == Bool && processSettings             procSettings = MultiPhotonIonization.Settings()
        elseif  process == Basics.MultiPI           && typeof(processSettings) == MultiPhotonIonization.Settings      procSettings = processSettings
        elseif  process == Basics.MultiPDI          && typeof(processSettings) == Bool && processSettings             procSettings = MultiPhotonDoubleIon.Settings()
        elseif  process == Basics.MultiPDI          && typeof(processSettings) == MultiPhotonDoubleIon.Settings       procSettings = processSettings
        elseif  process == Basics.InternalConv      && typeof(processSettings) == Bool && processSettings             procSettings = InternalConversion.Settings()
        elseif  process == Basics.InternalConv      && typeof(processSettings) == InternalConversion.Settings         procSettings = processSettings
        elseif  process == Basics.Coulex            && typeof(processSettings) == Bool && processSettings             procSettings = CoulombExcitation.Settings()
        elseif  process == Basics.Coulex            && typeof(processSettings) == CoulombExcitation.Settings          procSettings = processSettings
        elseif  process == Basics.Coulion           && typeof(processSettings) == Bool && processSettings             procSettings = CoulombIonization.Settings()
        elseif  process == Basics.Coulion           && typeof(processSettings) == CoulombIonization.Settings          procSettings = processSettings
        else    error("The processSettings must fit to the given process or should not occur if NoProcess is to be calculated.")
        end

        Computation(sa, nm, grid, properties, configs, asfSettings,
                    initialConfigs, initialAsfSettings,
                    intermediateConfigs, intermediateAsfSettings,
                    finalConfigs, finalAsfSettings,
                    alphaSettings, einsteinSettings, formSettings, greenSettings, hfsSettings, isotopeSettings, plasmaSettings,
                    polaritySettings, yieldSettings, zeemanSettings,
                    ##x amplitude, ampSettings,
                    process, procSettings)
    end

    # `Base.show(io::IO, computation::Atomic.Computation)`  ... prepares a proper printout of the variable computation::Atomic.Computation.
    function Base.show(io::IO, computation::Atomic.Computation)
        println(io, "name:                           $(computation.name)  ")
        println(io, "nuclearModel:                   $(computation.nuclearModel)  ")
        println(io, "grid:                           $(computation.grid)  ")
        ##x println(io, "calcLevelProperties:            $(computation.calcLevelProperties)  ")
        println(io, "properties:                     $(computation.properties)  ")
        println(io, "configs:                        $(computation.configs)  ")
        println(io, "asfSettings:                      computation.asfSettings  ")
        println(io, "initialConfigs:                 $(computation.initialConfigs)  ")
        println(io, "initialAsfSettings:               computation.initialAsfSettings  ")
        println(io, "intermediateConfigs:            $(computation.intermediateConfigs)  ")
        println(io, "intermediateAsfSettings:          computation.intermediateAsfSettings  ")
        println(io, "finalConfigs:                   $(computation.finalConfigs)  ")
        println(io, "finalAsfSettings:                 computation.finalAsfSettings  ")
        println(io, "alphaSettings:                    computation.alphaSettings  ")
        println(io, "einsteinSettings:                 computation.einsteinSettings  ")
        println(io, "formSettings:                     computation.formSettings  ")
        println(io, "greenSettings:                    computation.greenSettings  ")
        println(io, "hfsSettings:                      computation.hfsSettings  ")
        println(io, "isotopeSettings:                  computation.isotopeSettings  ")
        println(io, "plasmaSettings:                   computation.plasmaSettings  ")
        println(io, "polaritySettings:                 computation.polaritySettings  ")
        println(io, "yieldSettings:                    computation.yieldSettings  ")
        println(io, "zeemanSettings:                   computation.zeemanSettings  ")
        println(io, "process:                        $(computation.process)  ")
        println(io, "processSettings:                  computation.processSettings  ")
    end

    """
    `JAC.Atomic.Computation(gui::Guint; comp::Atomic.Computation=Atomic.Computation())`  ... constructor that is defined by a graphical user interface.
    """
    function Computation(gui::Guint; comp::Atomic.Computation=Atomic.Computation())

        if      gui == Gui
            println("aa")
            output = Atomic.ComputationGui(comp)
        elseif  gui == GuiSettings
            println("bb")
            output = Atomic.ComputationGuiSettings(comp)
        else  error("Unsupported Guint = $gui.")
        end

        return( output )
    end

    """
    `JAC.Atomic.ComputationGui(comp::Atomic.Computation)`  ... constructor that defines the standard part of a computation by a graphical
         user interface apart from the setting and the configurations.
    """
    function ComputationGui(comp::Atomic.Computation)
        nmd = comp.nuclearModel

        title = comp.name;    nmd = comp.nuclearModel;    grid = comp.grid;    calcLevelProperties = comp.calcLevelProperties
        properties          = comp.properties
        configs             = comp.configs;                   asfSettings             = comp.asfSettings
        initialConfigs      = comp.initialConfigs;            initialAsfSettings      = comp.initialAsfSettings
        intermediateConfigs = comp.intermediateConfigs;       intermediateAsfSettings = comp.intermediateAsfSettings
        finalConfigs        = comp.finalConfigs;              finalAsfSettings        = comp.finalAsfSettings
        alphaSettings       = comp.alphaSettings;             einsteinSettings        = comp.einsteinSettings
        formSettings        = comp.formSettings;              greenSettings           = comp.greenSettings
        hfsSettings         = comp.hfsSettings;               isotopeSettings         = comp.isotopeSettings
        plasmaSettings      = comp.plasmaSettings;            polaritySettings        = comp.polaritySettings
        yieldSettings       = comp.yieldSettings;             zeemanSettings          = comp.zeemanSettings
        process             = comp.process;                   processSettings         = comp.processSettings

        # Update the nuclear model
        tn1 = "Update nuclear model: "
        bn1 = slider(1:110, label = "charge", value = nmd.Z)
        bn2 = dropdown([nmd.model, "Fermi", "point", "uniform"])
        bn3 = spinbox(label="mass  ";   value=nmd.mass)
        bn4 = spinbox(label="radius";   value=nmd.radius)
        bn5 = spinbox(label="2*spin";   value=nmd.spinI.num)
        bn6 = spinbox(label="mu    ";   value=nmd.mu)
        bn7 = spinbox(label="Q     ";   value=nmd.Q)
        #
        #
        update = button("Update all")
        ui = vbox( hbox( pad(0em, tn1) ),
                   hbox( pad(1em, bn1) ),
                   hbox( pad(1em, bn2), pad(1em, bn3), pad(1em, bn4) ),
                   hbox( pad(1em, bn5), pad(1em, bn6), pad(1em, bn7) ),
                   hbox( pad(2em, update) )
                  )
        #
        Interact.display(ui)
        output = Interact.@map  (&update;     Computation(title,
                    Nuclear.Model( observe(bn1)[], observe(bn2)[], observe(bn3)[], observe(bn4)[],
                                   AngularJ64( Int64(observe(bn5)[])//2 ), observe(bn6)[], observe(bn7)[] ),
                    grid, calcLevelProperties, properties, configs, asfSettings,
                    initialConfigs, initialAsfSettings,
                    intermediateConfigs, intermediateAsfSettings,
                    finalConfigs, finalAsfSettings,
                    alphaSettings, einsteinSettings, formSettings, greenSettings, hfsSettings, isotopeSettings, plasmaSettings,
                    polaritySettings, yieldSettings, zeemanSettings,
                    process, processSettings) )
        #
        return( output )
    end

    """
    `JAC.Atomic.ComputationGuiSettings(comp::Atomic.Computation)`  ... constructor that defines the settings and configurations of a computation
         by a graphical user interface.
    """
    function ComputationGuiSettings(comp::Atomic.Computation)
        ac = comp;   nmd = ac.nuclearModel

        t1 = "A slider: "
        b1 = slider(1:110, label = "charge", value = nmd.Z)
        update = button("Update")
        ui = vbox( hbox( pad(0em, t1) ),
                   hbox( pad(0em, b1) ),
                   hbox( pad(1em, update) )
                  )
        #
        Interact.display(ui)
        output = Interact.@map  (&update; observe(b1)[] )
        return( output )
    end

end # module


text:
export update!

function update!(n::Node, grads::NamedTuple; η = 0.01, uparams = [])
    if !isnothing(grads.params)
        n.params .+= η * grads.params
    end

    for (i, child) in enumerate(n.children)
        if child isa AbstractLeaf
            n.children[i] = update!(child, grads.children[i]; η = η, uparams = uparams)
        else
            if grads.children[i] === nothing # why is this happening?
                nothing
            else
                update!(child, grads.children[i]; η = η, uparams = uparams)
            end
        end
    end
end

function update!(n::T, grads::NamedTuple; η = 0.01, uparams = []) where {T<:AbstractLeaf}

    v = collect(params(n))
    for (i,k) in enumerate(keys(params(n)))
        if (k ∈ uparams) || isempty(uparams)
            v[i] += grads.params[k] === nothing ? zero(eltype(v)) : η * grads.params[k]
        end
    end

    return T(n.scope, (;zip(keys(params(n)), v)...))
end

ext: jl


text:
#=Given a number N, find its factorial using recursion. It is quite easy
but we want to add recursion algorithms in Julia.=#

## Function

function fact(n)
    if (n == 0)
        return 1
    end
    if (n == 1)
        return 1
    end
    return n * fact(n - 1)
end

## Input

n = readline()
n = parse(Int64, n)

## Calling the Function

fact(n)

#=
Sample Test Case:
Input:
    n=10
Output:
    3628800
Time complexity: O( N )
=#


text:
using SparseArrays
using IterativeSolvers

struct sparsetensor_vec
    # dims is the array of dimensions of the tensor
    dims :: Array{Int64}
    # spv stores the values of the tensor in a sparse vector
    spv :: SparseVector{Float64}
end

struct sparsetensor_mat
    # dims1 is the array of dimensions of the tensor
    # used for the matrix dimension one
    dims1 :: Array{Int64}
    # dims2 is the array of dimensions of the tensor
    # used for the matrix dimension two
    dims2 :: Array{Int64}
    # spm stores the values of the tensor in a sparse matrix
    spm :: SparseMatrixCSC{Float64}
end

function spvec_to_spmat(x :: sparsetensor_vec, lindex)
    # Constructs the sparse matrix x[1:lindex, lindex+1:end]
    res_dims1 = x.dims[1:lindex]
    res_dims2 = x.dims[lindex+1:end]
    m = prod(res_dims1)
    n = prod(res_dims2)
    n_nz = nnz(x.spv)
    z = map(collect, map(a -> divrem(a,m),x.spv.nzind))
    i = zeros(Int, n_nz)
    j = zeros(Int, n_nz)
    for k = 1:n_nz
        i[k] = Int(z[k][2])
        j[k] = Int(z[k][1]) + 1
		if i[k] == 0
			i[k] = m
			j[k] += (-1)
		end
    end
    res_ma = sparse(i, j, x.spv.nzval, m, n)
    return sparsetensor_mat(res_dims1, res_dims2, res_ma)
end

function spmat_to_spvec(x :: sparsetensor_mat)
	# Transforms the sparse matrix into a sparse vector
	res_dim = vcat(x.dims1,x.dims2)
    return sparsetensor_vec(res_dim, reshape(x.spm,prod(res_dim)))
end

function sppermutedims(x :: sparsetensor_vec, perm)
    # permutes the dimensions of the sparse tensor

    # Prepare the transformation to cartesian coordinates
    n_prod = ones(length(perm))
    for k = 1:(length(perm))
        n_prod[1:(end-k)] = n_prod[1:(end-k)] .* x.dims[k]
    end
    # Prepare the transformation to vector coordinates
    n_prod_perm = ones(length(perm))
    for k = 1:(length(perm))
        n_prod_perm[1:(end-k)] = n_prod_perm[1:(end-k)] .* x.dims[perm[k]]
    end
    # Initialize the cartesian coordinates and their permutation
    x_katind = ones(length(x.spv.nzind), length(perm))
    x_katind_perm = ones(length(x.spv.nzind), length(perm))
    # Initialize the vector indices
    res_ind = ones(length(x.spv.nzind))

    for x_ind = 1:length(x.spv.nzind)
        #Calculate the cartesian coordinates
        x_ind_rem = x.spv.nzind[x_ind]
        for k = 1:(length(perm))
            x_katind[x_ind, length(perm) - k + 1], x_ind_rem = divrem(x_ind_rem, n_prod[k])
			if x_ind_rem == 0
				x_ind_rem = n_prod[k]
			else
				x_katind[x_ind, length(perm) - k + 1] += 1
			end
        end
        # Apply the permutation on the cartesian coordinates
        for k = 1:length(perm)
            x_katind_perm[x_ind, k] = x_katind[x_ind, perm[k]]
        end
        # Retransform the cartesian coordinates into vector indices
        for k = 1:length(perm)
            res_ind[x_ind] = res_ind[x_ind] + n_prod_perm[k]*(x_katind_perm[x_ind,length(perm) - k + 1]-1)
        end
    end
	return sparsetensor_vec(x.dims[perm], sparsevec(res_ind, x.spv.nzval))
end

function test_sppermutedims()
	n = 5
	x = randn(n,n,n,n)
	x_spv = sparsetensor_vec([n,n,n,n],sparse(x[:]))
	@test isapprox(Array(sppermutedims(x_spv,[2,3,1,4]).spv),permutedims(x,[2,3,1,4])[:])
end

# Redefine permutedims for a sparsetensor_mat
function Base.permutedims(A::sparsetensor_mat, perm, i_new)
	return permutedims(A.spm, perm, A.dims1, A.dims2, i_new)
end

# Redefine permutedims for the structs SparseVector
function Base.permutedims(A::SparseVector, perm, dims)
	return sppermutedims(sparsetensor_vec(dims, A), perm)
end

# Redefine permutedims for the structs SparseMatrixCSC
function Base.permutedims(A::SparseMatrixCSC, perm, dims1, dims2, i_new)
	return spvec_to_spmat(sppermutedims(spmat_to_spvec(sparsetensor_mat(dims1, dims2, A)), perm), i_new).spm
end

# Define ttv_decomp for sparse input
function ttv_spdecomp(spv :: sparsetensor_vec, index;tol=1e-10)
	# Decomposes a tensor into its tensor train with core matrices at i=index
	# For index < d it doesn't work yet
	dims = spv.dims #dims = [n_1,...,n_d]
	n_max = maximum(dims)
	d = length(dims)
	ttv_vec = Array{Array{Float64}}(undef,d)
	# ttv_ot[i]= -1 if i < index
	# ttv_ot[i] = 0 if i = index
	# ttv_ot[i] = 1 if i > index
	ttv_ot = -ones(d)
	ttv_ot[index] = 0
	if index < d
		ttv_ot[index+1:d] = ones(d-index)
	end
	rks = ones(Int, d+2) # ttv_rks will be rks[2:d+1]
	tensor_curr = convert(SparseMatrixCSC, reshape(spv.spv, dims[1], :))
	# Calculate ttv_vec[i] for i < index
	for i = 1 : (index - 1)
		# Reshape the currently left tensor
		tensor_curr = convert(SparseMatrixCSC, reshape(tensor_curr, Int(rks[i] * dims[i]), :))
		# Perform a qr decomposition
		qro = qr(tensor_curr)
		# Define the i-th rank
		diag_R = abs.(diag(qro.R))
		rks[i+1] = sum(diag_R.>tol)
		Q = qro.Q[invperm(qro.prow),(1:length(diag_R))[diag_R.>tol]]
		R = qro.R[diag_R.>tol,invperm(qro.pcol)]
		# Initialize ttv_vec[i]
		ttv_vec[i] = zeros(dims[i],rks[i],rks[i+1])
		# Fill in the ttv_vec[i]
		for x = 1 : dims[i]
			ttv_vec[i][x, :, :] = Q[(rks[i]*(x-1) + 1):(rks[i]*x), :]
		end
		# Update the currently left tensor
		tensor_curr = R
	end
	# Calculate ttv_vec[i] for i = index
	# Reshape the currently left tensor
	tensor_curr = reshape(tensor_curr, Int(dims[index]*rks[index]),:)
	# Initialize ttv_vec[i]
	ttv_vec[index] = zeros(dims[index], rks[index], rks[index+1])
	# Fill in the ttv_vec[i]
	for x = 1 : dims[index]
		ttv_vec[index][x, :, :] =
		tensor_curr[Int(rks[index]*(x-1) + 1):Int(rks[index]*x), 1:rks[index+1]]
	end
	# Define the return value as a ttvector
	return ttvector(ttv_vec, dims, rks[2:d+1], ttv_ot)
end

function test_ttv_spdecomp()
	n = 5
	x = randn(n,n,n,n)
	x_spv = sparsetensor_vec([n,n,n,n],sparse(x[:]))
	x_sp_tt = ttv_spdecomp(x_spv,4)
	size(ttv_to_tensor(x_sp_tt))
	@test isapprox(ttv_to_tensor(x_sp_tt),x)
end

function tto_spdecomp(spm :: sparsetensor_mat, index)
	# Decomposes a sparse tensor operator into its tensor train
	# with core matrices at i=index
	# For index < d it doesn't work yet
	# The tensor has to be given as element of R^[n_1*...*n_d,n_1*...*n_d]
	d = length(spm.dims1)
	tto_dims = spm.dims1
	n_max = maximum(tto_dims)
	dims_sq = tto_dims.^2
	# The matrix is reshaped, reordered and reshaped
	# into the matrix [(x_1,y_1),((x_2,y_2),...,(x_d,y_d))]
	# and decomposed into its tensor train with core matrices at i = index
	index_sorted = vec(reshape(Transpose(reshape(1:(2*d),:,2)),1,:))
	dims_sorted = zeros(2*d)
    dims_sorted[1:d] = tto_dims
	dims_sorted[(d+1):(2*d)] = tto_dims
	dims_sorted = dims_sorted[index_sorted]
	ttv = ttv_spdecomp(spmat_to_spvec(sparsetensor_mat([dims_sq[1]], dims_sq[2:end], permutedims(spm.spm, index_sorted, spm.dims1, spm.dims2, d))), index)
	# Define the array of ranks [r_0=1,r_1,...,r_d]
	rks = ones(Int, d+1)
	rks[2:(d+1)] = ttv.ttv_rks
	r_max = maximum(rks)
	# Initialize tto_vec
	tto_vec = Array{Array{Float64}}(undef,d)
	# Fill in tto_vec
	for i = 1:d
		# Initialize tto_vec[i]
		tto_vec[i] = zeros(tto_dims[i], tto_dims[i], rks[i], rks[i+1])
		# Fill in tto_vec[i]
		tto_vec[i][:, :, :, :] =
		reshape(ttv.ttv_vec[i][1:dims_sq[i], 1:rks[i], 1:rks[i+1]],
		tto_dims[i], tto_dims[i], rks[i], :)
	end
	# Define the return value as a ttoperator
	return ttoperator(tto_vec, tto_dims, rks[2:(d+1)], ttv.ttv_ot)
end

function Lap_sp(n::Integer, d::Integer)
	# Returns the tensor of the discrete Laplacian in a box [0,1]^d
	# with n equidistant discretization points in each direction
	# as an element of R^[n^d,n^d]
	I = zeros((n^d) * (2 * d + 1))
	J = zeros((n^d) * (2 * d + 1))
	V = zeros((n^d) * (2 * d + 1))

	for x = 0:n^d-1
		x_dig = digits(x, base = n, pad = d)
		I[(2*d + 1) * x + 1] = x + 1
		J[(2*d + 1) * x + 1] = x + 1
		V[(2*d + 1) * x + 1] = 2*d
		for i = 0:d-1
			if x_dig[i + 1] == 0
				I[(2*d + 1) * x + 2*i + 2] = x + 1
				J[(2*d + 1) * x + 2*i + 2] = x + n^i + 1
				V[(2*d + 1) * x + 2*i + 2] = -1
			elseif x_dig[i + 1] == n - 1
				I[(2*d + 1) * x + 2*i + 2] = x + 1
				J[(2*d + 1) * x + 2*i + 2] = x - n^i + 1
				V[(2*d + 1) * x + 2*i + 2] = -1
			else
				I[(2*d + 1) * x + 2*i + 2] = x + 1
				J[(2*d + 1) * x + 2*i + 2] = x + n^i + 1
				V[(2*d + 1) * x + 2*i + 2] = -1
				I[(2*d + 1) * x + 2*i + 3] = x + 1
				J[(2*d + 1) * x + 2*i + 3] = x - n^i + 1
				V[(2*d + 1) * x + 2*i + 3] = -1
			end
		end
	end
	I = I[I .!= 0]
	J = J[J .!= 0]
	V = V[V .!= 0]
	return sparse(I , J, V, n^d, n^d)
end

"""
# Test
n = 25
d = 2
L = Lap(n, d)
x_o = ones(size(L,1))
x_o_tt = ttv_decomp(reshape(x_o, n * ones(Int,d)...), d)
b = L * x_o
x_s = map(round, 10 * randn(size(L,1)))
x_s_tt = ttv_decomp(reshape(x_s, n * ones(Int, d)...),d)
b_tt = ttv_decomp(reshape(b,n*ones(Int,d)...), d)
L_spm = sparsetensor_mat(n*ones(Int, d), n*ones(Int, d), L)
L_tt = tto_spdecomp(L_spm, d)
x_als_tt = als(L_tt, b_tt, x_s_tt, x_o_tt.ttv_rks)
"""

max_stars_repo_name: msdupuy/Tensor-Train-Julia

max_issues_repo_name: msdupuy/Tensor-Train-Julia

max_forks_repo_name: msdupuy/Tensor-Train-Julia


text:
using Documenter, FinEtools, FinEtoolsAcousticsTutorials

makedocs(
	modules = [FinEtoolsAcousticsTutorials],
	doctest = false, clean = true,
	format = Documenter.HTML(prettyurls = false),
	authors = "Petr Krysl",
	sitename = "FinEtoolsAcousticsTutorials.jl",
	pages = Any[
	"Home" => "index.md",
	"Tutorials" => "tutorials/tutorials.md",
	]
	)

deploydocs(
    repo = "github.com/PetrKryslUCSD/FinEtoolsAcousticsTutorials.jl.git",
)


text:
import Statistics

Statistics._var(A::CuArray, corrected::Bool, mean, dims) =
    sum((A .- something(mean, Statistics.mean(A, dims=dims))).^2, dims=dims)/(prod(size(A)[[dims...]])-corrected)

Statistics._var(A::CuArray, corrected::Bool, mean, ::Colon) =
    sum((A .- something(mean, Statistics.mean(A))).^2)/(length(A)-corrected)

Statistics._std(A::CuArray, corrected::Bool, mean, dims) =
    sqrt.(Statistics.var(A; corrected=corrected, mean=mean, dims=dims))

Statistics._std(A::CuArray, corrected::Bool, mean, ::Colon) =
    sqrt.(Statistics.var(A; corrected=corrected, mean=mean, dims=:))


text:
abstract type DepthFormat end

struct Depth{DT} <: DepthFormat
    depth::DT
end

# TODO maybe we should implement this as a 32 bit wide primitive type
# and overload getproperty (getfield on 0.7) to implement depthstencil.depth with masking
# since you almost always want to have depthstencil.depth::Float32
struct DepthStencil{DT, ST} <: DepthFormat
    depth::DT
    stencil::ST
end
#0.7: Base.getproperty(x::DepthStencil, field::Symbol) = field == :depth ? Float32(x.depth) : x.stencil #I actually do want to support v0.6 for now

"""
Float24 storage type for depth
"""
primitive type Float24 <: AbstractFloat 24 end

gl_internal_format(::Type{Depth{Float32}}) = GL_DEPTH_COMPONENT32F
gl_internal_format(::Type{DepthStencil{Float24, N0f8}}) = GL_DEPTH24_STENCIL8

function gl_internal_format(::T) where T
    error("$T doesn't have a valid mapping to an OpenGL internal format enum. Please use DepthStencil/Depth/Color, or overload `gl_internal_format(x::$T)`
    to return the correct OpenGL format enum.
    ")
end

gl_attachment(::Type{<:Depth}) = GL_DEPTH_ATTACHMENT
gl_attachment(::Type{<:DepthStencil}) = GL_DEPTH_STENCIL_ATTACHMENT
function gl_attachment(::T) where T
    error("$T doesn't have a valid mapping to an OpenGL attachment enum. Please use DepthStencil/Depth, or overload `gl_attachment(x::$T)`
    to return the correct OpenGL depth attachment.
    ")
end

const GL_TEXTURE_MAX_ANISOTROPY_EXT = GLenum(0x84FE)

colordim(::Type{T}) where {T} = cardinality(T)
colordim(::Type{T}) where {T <: Real} = 1

#Supported texture modes/dimensions
function texturetype_from_dimensions(ndim::Integer)
    ndim == 1 && return GL_TEXTURE_1D
    ndim == 2 && return GL_TEXTURE_2D
    ndim == 3 && return GL_TEXTURE_3D
    @error "Dimensionality: $(ndim), not supported for OpenGL texture"
end

@generated function textureformat_internal_from_type(::Type{T}) where T
    if T <: Depth{Float32}
        return :GL_DEPTH_COMPONENT32F
    elseif T <: DepthStencil{Float24, N0f8}
        return :GL_DEPTH24_STENCIL8
    end
    glasserteltype(T)
    dim = length(T)
    @assert (dim <= 4 && dim >= 1) "No Textureformat that fits $dim-dimensional eltypes."

    eltyp = eltype(T) #like Float64 in RGBA{Float64}
    sym = "GL_"
    sym *= "RGBA"[1:dim]
    bits = sizeof(eltyp) * 8
    @assert bits <= 32 "$T has too many bits ($bits)"
    sym *= string(bits)

    if eltyp <: AbstractFloat
        sym *= "F"
    elseif eltyp <: FixedPoint
        sym *= eltyp <: Normed ? "" : "_SNORM"
    elseif eltyp <: Signed
        sym *= "I"
    elseif eltyp <: Unsigned
        sym *= "UI"
    end
    glenumsym = Symbol(sym)
    @assert isdefined(ModernGL, glenumsym) "$T doesn't have a proper mapping to an OpenGL format."
    return :($glenumsym)
end

function textureformat_from_type_sym(dim::Integer, isinteger::Bool, order::AbstractString)
    @assert dim <= 4 "no colors with dimension > 4 allowed. Dimension given: $dim"
    sym = "GL_"
    # Handle that colordim == 1 => RED instead of R
    color = dim == 1 ? "RED" : order[1:dim]
    # Handle gray value
    integer = isinteger ? "_INTEGER" : ""
    sym *= color * integer
    return Symbol(sym)
end

textureformat_from_type_sym(::Type{T}) where {T <: Real} = textureformat_from_type_sym(1, T <: Integer, "RED")
textureformat_from_type_sym(::Type{T}) where {T <: AbstractArray} = textureformat_from_type_sym(cardinality(T), eltype(T) <: Integer, "RGBA")

function textureformat_from_type_sym(::Type{T}) where T
    glasserteltype(T)
    typenamestring = string(Base.typename(T).name)
    if typenamestring ∉ ("RGB", "BGR", "RGBA", "ARGB", "ABGR", "RGBA","BGRA")
        typenamestring = "RGBA"
    end
    return textureformat_from_type_sym(cardinality(T), eltype(T) <: Integer, typenamestring)
end

textureformat_from_type_sym(::Type{<:DepthFormat}) = :GL_DEPTH_COMPONENT
@generated function textureformat_from_type(::Type{T}) where T
    sym = textureformat_from_type_sym(T)
    @assert isdefined(ModernGL, sym) "$T doesn't have a proper mapping to an OpenGL format"
    return :($sym)
end

struct TextureParameters{NDim}
    minfilter    ::Symbol
    magfilter    ::Symbol # magnification
    repeat       ::NTuple{NDim, Symbol}
    anisotropic  ::Float32
    swizzle_mask ::Vector{GLenum}
end
function TextureParameters(T, NDim;
        minfilter   = T <: Integer ? :nearest : :linear,
        magfilter   = minfilter, # magnification
        x_repeat    = :clamp_to_edge, #wrap_s
        y_repeat    = x_repeat, #wrap_t
        z_repeat    = x_repeat, #wrap_r
        anisotropic = 1f0
    )
    T <: Integer && (minfilter == :linear || magfilter == :linear) && (@error "Wrong Texture Parameter: Integer texture can't interpolate. Try :nearest")
    repeat = (x_repeat, y_repeat, z_repeat)
    dim = T <: DepthFormat ? 1 : length(T)
    swizzle_mask = if dim == 3 #<: Gray
        GLenum[GL_RED,  GL_GREEN, GL_BLUE, GL_ONE]
    elseif dim == 4 #<: GrayA
        GLenum[GL_RED, GL_GREEN,  GL_BLUE, GL_ALPHA]
    else
        GLenum[]
    end
    TextureParameters(
        minfilter, magfilter, ntuple(i->repeat[i], NDim),
        anisotropic, swizzle_mask
    )
end

map_texture_paramers(s::NTuple{N, Symbol}) where {N} = map(map_texture_paramers, s)

function map_texture_paramers(s::Symbol)

    s == :clamp_to_edge          && return GL_CLAMP_TO_EDGE
    s == :mirrored_repeat        && return GL_MIRRORED_REPEAT
    s == :repeat                 && return GL_REPEAT
    s == :linear                 && return GL_LINEAR
    s == :nearest                && return GL_NEAREST
    s == :nearest_mipmap_nearest && return GL_NEAREST_MIPMAP_NEAREST
    s == :linear_mipmap_nearest  && return GL_LINEAR_MIPMAP_NEAREST
    s == :nearest_mipmap_linear  && return GL_NEAREST_MIPMAP_LINEAR
    s == :linear_mipmap_linear   && return GL_LINEAR_MIPMAP_LINEAR

    @error "$s is not a valid texture parameter"
end

#This is used in the construction of Textures
function set_packing_alignment(a) # at some point we should specialize to array/ptr a
    glPixelStorei(GL_UNPACK_ALIGNMENT,   1)
    glPixelStorei(GL_UNPACK_ROW_LENGTH,  0)
    glPixelStorei(GL_UNPACK_SKIP_PIXELS, 0)
    glPixelStorei(GL_UNPACK_SKIP_ROWS,   0)
end

abstract type OpenglTexture{T, NDIM} <: GPUArray{T, NDIM} end

mutable struct Texture{T, NDIM} <: OpenglTexture{T, NDIM}
    id              ::GLuint
    texturetype     ::GLenum
    pixeltype       ::GLenum
    internalformat  ::GLenum
    format          ::GLenum
    parameters      ::TextureParameters{NDIM}
    size            ::NTuple{NDIM, Int}
    context         ::AbstractContext
    function Texture{T, NDIM}(
            id              ::GLuint,
            texturetype     ::GLenum,
            pixeltype       ::GLenum,
            internalformat  ::GLenum,
            format          ::GLenum,
            parameters      ::TextureParameters{NDIM},
            size            ::NTuple{NDIM, Int}
        )  where {T, NDIM}
        tex = new(
            id,
            texturetype,
            pixeltype,
            internalformat,
            format,
            parameters,
            size,
            current_context()
        )
        finalizer(free!, tex)
        tex
    end
end

function Texture(
        data::Ptr{T}, dims::NTuple{NDim, Int};
        internalformat::GLenum = textureformat_internal_from_type(T),
        texturetype   ::GLenum = texturetype_from_dimensions(NDim),
        format        ::GLenum = textureformat_from_type(T),
        mipmap = false,
        parameters... # rest should be texture parameters
    ) where {T, NDim}
    texparams = TextureParameters(T, NDim; parameters...)
    id = glGenTextures()
    glBindTexture(texturetype, id)
    set_packing_alignment(data)
    if T <: DepthFormat
        numbertype = GL_FLOAT
    else
        numbertype = julia2glenum(eltype(T))
    end
    glTexImage(texturetype, 0, internalformat, dims..., 0, format, numbertype, data)
    mipmap && glGenerateMipmap(texturetype)
    texture = Texture{T, NDim}(
        id, texturetype, numbertype, internalformat, format,
        texparams,
        dims
    )
    set_parameters(texture)
    texture::Texture{T, NDim}
end

"""
Constructor for Array Texture
"""
function Texture(
        data::Vector{Array{T, 2}};
        internalformat::GLenum = textureformat_internal_from_type(T),
        texturetype::GLenum    = GL_TEXTURE_2D_ARRAY,
        format::GLenum         = textureformat_from_type(T),
        parameters...
    ) where T
    glasserteltype(T)
    texparams = TextureParameters(T, 2; parameters...)
    id = glGenTextures()

    glBindTexture(texturetype, id)

    numbertype = julia2glenum(eltype(T))

    layers  = length(data)
    dims    = map(size, data)
    maxdims = foldl((0,0), dims) do v0, x
        a = max(v0[1], x[1])
        b = max(v0[2], x[2])
        (a,b)
    end
    set_packing_alignment(data)
    glTexStorage3D(GL_TEXTURE_2D_ARRAY, 1, internalformat, maxdims..., layers)
    for (layer, texel) in enumerate(data)
        width, height = size(texel)
        glTexSubImage3D(GL_TEXTURE_2D_ARRAY, 0, 0, 0, layer-1, width, height, 1, format, numbertype, texel)
    end

    texture = Texture{T, 2}(
        id, texturetype, numbertype,
        internalformat, format, texparams,
        tuple(maxdims...)
    )
    set_parameters(texture)
    texture
end

"""
Constructor for empty initialization with NULL pointer instead of an array with data.
You just need to pass the wanted color/vector type and the dimensions.
To which values the texture gets initialized is driver dependent
"""
function Texture(::Type{T}, dims::NTuple{N, Int}; kw_args...) where {T, N}
    glasserteltype(T)
    Texture(convert(Ptr{T}, C_NULL), dims; kw_args...)::Texture{T, N}
end

"""
Constructor for a normal array, with color or Abstract Arrays as elements.
So Array{Real, 2} == Texture2D with 1D Colorant dimension
Array{Vec1/2/3/4, 2} == Texture2D with 1/2/3/4D Colorant dimension
Colors from Colors.jl should mostly work as well
"""
function Texture(image::Array{T, NDim}; kw_args...) where {T, NDim}
    glasserteltype(T)
    Texture(pointer(image), size(image); kw_args...)::Texture{T, NDim}
end
Base.size(t::Texture) = t.size
width(t::Texture)     = size(t, 1)
height(t::Texture)    = size(t, 2)
depth(t::Texture)     = size(t, 3)
id(t::Texture)        = t.id

function Base.show(io::IO, t::Texture{T,D}) where {T,D}
    println(io, "Texture$(D)D: ")
    println(io, "                  ID: ", t.id)
    println(io, "                  Size: Dimensions: $(size(t))")
    # println(io, "                Size: ", reduce("Dimensions: ", size(t)) do v0, v1
        # v0*"x"*string(v1)
    # end)
    println(io, "    Julia pixel type: ", T)
    println(io, "   OpenGL pixel type: ", GLENUM(t.pixeltype).name)
    println(io, "              Format: ", GLENUM(t.format).name)
    println(io, "     Internal format: ", GLENUM(t.internalformat).name)
    println(io, "          Parameters: ", t.parameters)
end

is_texturearray(t::Texture) = t.texturetype == GL_TEXTURE_2D_ARRAY
is_texturebuffer(t::Texture) = t.texturetype == GL_TEXTURE_BUFFER

bind(t::Texture) = glBindTexture(t.texturetype, t.id)
bind(t::Texture, id) = glBindTexture(t.texturetype, id)

function resize_nocopy!(t::Texture{T, ND}, newdims::NTuple{ND, Int}) where {T, ND}
    bind(t)
    glTexImage(t.texturetype, 0, t.internalformat, newdims..., 0, t.format, t.pixeltype, C_NULL)
    t.size = newdims
    bind(t, 0)
    t
end

# Resize Texture
function gpu_resize!(t::Texture{T, ND}, newdims::NTuple{ND, Int}) where {T, ND}
    # dangerous code right here...Better write a few tests for this
    newtex   = similar(t, newdims)
    old_size = size(t)
    gpu_setindex!(newtex, t)
    t.size   = newdims
    free(t)
    t.id     = newtex.id
    return t
end
function gpu_setindex!(t::Texture{T, 1}, newvalue::Array{T, 1}, indexes::UnitRange{I}) where {T, I <: Integer}
    glBindTexture(t.texturetype, t.id)
    texsubimage(t, newvalue, indexes)
    glBindTexture(t.texturetype, 0)
end
function gpu_setindex!(t::Texture{T, N}, newvalue::Array{T, N}, indexes::Union{UnitRange,Integer}...) where {T, N}
    glBindTexture(t.texturetype, t.id)
    texsubimage(t, newvalue, indexes...)
    glBindTexture(t.texturetype, 0)
end

function gpu_setindex!(target::Texture{T, 2}, source::Texture{T, 2}, fbo=glGenFramebuffers()) where T
    glBindFramebuffer(GL_FRAMEBUFFER, fbo);
    glFramebufferTexture2D(GL_READ_FRAMEBUFFER, GL_COLOR_ATTACHMENT0,
                           GL_TEXTURE_2D, source.id, 0);
    glFramebufferTexture2D(GL_DRAW_FRAMEBUFFER, GL_COLOR_ATTACHMENT1,
                           GL_TEXTURE_2D, target.id, 0);
    glDrawBuffer(GL_COLOR_ATTACHMENT1);
    w, h = map(minimum, zip(size(target), size(source)))
    glBlitFramebuffer(0, 0, w, h, 0, 0, w, h,
                      GL_COLOR_BUFFER_BIT, GL_NEAREST)
end

# Implementing the GPUArray interface
function gpu_data(t::Texture{T, ND}) where {T, ND}
    result = zeros(T, size(t)...)
    unsafe_copyto!(result, t)
    return result
end

function Base.unsafe_copyto!(dest::Array{T, N}, source::Texture{T, N}) where {T,N}
    bind(source)
    glGetTexImage(source.texturetype, 0, source.format, source.pixeltype, dest)
    bind(source, 0)
    nothing
end

similar(t::Texture{T, NDim}, newdims::Int...) where {T, NDim} = similar(t, newdims)

texsubimage(t::Texture{T, 1}, newvalue::Array{T, 1}, xrange::UnitRange, level=0) where {T} = glTexSubImage1D(
    t.texturetype, level, first(xrange)-1, length(xrange), t.format, t.pixeltype, newvalue
)
function texsubimage(t::Texture{T, 2}, newvalue::Array{T, 2}, xrange::UnitRange, yrange::UnitRange, level=0) where T
    glTexSubImage2D(
        t.texturetype, level,
        first(xrange)-1, first(yrange)-1, length(xrange), length(yrange),
        t.format, t.pixeltype, newvalue
    )
end
texsubimage(t::Texture{T, 3}, newvalue::Array{T, 3}, xrange::UnitRange, yrange::UnitRange, zrange::UnitRange, level=0) where {T} = glTexSubImage3D(
    t.texturetype, level,
    first(xrange)-1, first(yrange)-1, first(zrange)-1, length(xrange), length(yrange), length(zrange),
    t.format, t.pixeltype, newvalue
)

function TextureParameters(t::Texture{T, NDim}; kw_args...) where {T, NDim}
    TextureParameters(T, NDim; kw_args...)
end

function set_parameters(t::Texture{T, N}, params::TextureParameters=t.parameters) where {T, N}
    fnames    = (:minfilter, :magfilter, :repeat)
    data      = Dict([(name, map_texture_paramers(getfield(params, name))) for name in fnames])
    result    = Tuple{GLenum, Any}[]
    push!(result, (GL_TEXTURE_MIN_FILTER, data[:minfilter]))
    push!(result, (GL_TEXTURE_MAG_FILTER, data[:magfilter]))
    push!(result, (GL_TEXTURE_WRAP_S, data[:repeat][1]))
    if !isempty(params.swizzle_mask)
        push!(result, (GL_TEXTURE_SWIZZLE_RGBA, params.swizzle_mask))
    end
    N >= 2 && push!(result, (GL_TEXTURE_WRAP_T, data[:repeat][2]))
    if N >= 3 && !is_texturearray(t) # for texture arrays, third dimension can not be set
        push!(result, (GL_TEXTURE_WRAP_R, data[:repeat][3]))
    end
    # push!(result, (GL_TEXTURE_MAX_ANISOTROPY_EXT, params.anisotropic))
    t.parameters = params
    set_parameters(t, result)
end
function texparameter(t::Texture, key::GLenum, val::GLenum)
    glTexParameteri(t.texturetype, key, val)
end
function texparameter(t::Texture, key::GLenum, val::Vector)
    glTexParameteriv(t.texturetype, key, val)
end
function texparameter(t::Texture, key::GLenum, val::Float32)
    glTexParameterf(t.texturetype, key, val)
end
function set_parameters(t::Texture, parameters::Vector{Tuple{GLenum, Any}})
    bind(t)
    for elem in parameters
        texparameter(t, elem...)
    end
    bind(t, 0)
end

function similar(t::Texture{T, NDim}, newdims::NTuple{NDim, Int}) where {T, NDim}
    Texture(
        Ptr{T}(C_NULL),
        newdims, t.texturetype,
        t.pixeltype,
        t.internalformat,
        t.format,
        t.parameters
    )
end

function free!(x::Texture)
    if !is_current_context(x.context)
        return x
    end
    id = [x.id]
    try
        glDeleteTextures(x.id)
    catch e
        free_handle_error(e)
    end
    return
end

# for bufferSampler, aka Texture Buffer
mutable struct TextureBuffer{T} <: OpenglTexture{T, 1}
    texture::Texture{T, 1}
    buffer::Buffer{T}
end
function TextureBuffer(buffer::Buffer{T}) where T
    glasserteltype(T)
    texture_type = GL_TEXTURE_BUFFER
    id = glGenTextures()
    glBindTexture(texture_type, id)
    internalformat = default_internalcolorformat(T)
    glTexBuffer(texture_type, internalformat, buffer.id)
    tex = Texture{T, 1}(
        id, texture_type, julia2glenum(T), internalformat,
        default_colorformat(T), TextureParameters(T, 1),
        size(buffer)
    )
    TextureBuffer(tex, buffer)
end
function TextureBuffer(buffer::Vector{T}) where T
    glasserteltype(T)
    buff = Buffer(buffer, buffertype = GL_TEXTURE_BUFFER, usage = GL_DYNAMIC_DRAW)
    TextureBuffer(buff)
end
Base.size(t::TextureBuffer) = size(t.buffer)
Base.size(t::TextureBuffer, i::Integer) = size(t.buffer, i)
Base.length(t::TextureBuffer) = length(t.buffer)

# GPUArray interface:
function Base.unsafe_copyto!(a::Vector{T}, readoffset::Int, b::TextureBuffer{T}, writeoffset::Int, len::Int) where T
    copy!(a, readoffset, b.buffer, writeoffset, len)
    glBindTexture(b.texture.texturetype, b.texture.id)
    glTexBuffer(b.texture.texturetype, b.texture.internalformat, b.buffer.id) # update texture
end

function Base.unsafe_copyto!(a::TextureBuffer{T}, readoffset::Int, b::Vector{T}, writeoffset::Int, len::Int) where T
    copy!(a.buffer, readoffset, b, writeoffset, len)
    glBindTexture(a.texture.texturetype, a.texture.id)
    glTexBuffer(a.texture.texturetype, a.texture.internalformat, a.buffer.id) # update texture
end
function Base.unsafe_copyto!(a::TextureBuffer{T}, readoffset::Int, b::TextureBuffer{T}, writeoffset::Int, len::Int) where T
    unsafe_copy!(a.buffer, readoffset, b.buffer, writeoffset, len)

    glBindTexture(a.texture.texturetype, a.texture.id)
    glTexBuffer(a.texture.texturetype, a.texture.internalformat, a.buffer.id) # update texture

    glBindTexture(b.texture.texturetype, btexture..id)
    glTexBuffer(b.texture.texturetype, b.texture.internalformat, b.buffer.id) # update texture
    glBindTexture(t.texture.texturetype, 0)
end
function gpu_setindex!(t::TextureBuffer{T}, newvalue::Vector{T}, indexes::UnitRange{I}) where {T, I <: Integer}
    glBindTexture(t.texture.texturetype, t.texture.id)
    t.buffer[indexes] = newvalue # set buffer indexes
    glTexBuffer(t.texture.texturetype, t.texture.internalformat, t.buffer.id) # update texture
    glBindTexture(t.texture.texturetype, 0)
end

gpu_data(t::TextureBuffer{T}) where {T} = gpu_data(t.buffer)
gpu_getindex(t::TextureBuffer{T}, i::UnitRange{Int64}) where {T} = t.buffer[i]

function similar(t::TextureBuffer{T}, newdims::NTuple{1, Int}) where T
    buff = similar(t.buffer, newdims...)
    return TextureBuffer(buff)
end

# Resize Texture
function gpu_resize!(t::TextureBuffer{T}, newdims::NTuple{1, Int}) where T
    resize!(t.buffer, newdims)
    glBindTexture(t.texture.texturetype, t.texture.id)
    glTexBuffer(t.texture.texturetype, t.texture.internalformat, t.buffer.id) #update data in texture
    t.texture.size  = newdims
    glBindTexture(t.texture.texturetype, 0)
    t
end

next(t::TextureBuffer{T}, state::Tuple{Ptr{T}, Int}) where {T} = next(t.buffer, state)
function done(t::TextureBuffer{T}, state::Tuple{Ptr{T}, Int}) where T
    isdone = done(t.buffer, state)
    if isdone
        glBindTexture(t.texturetype, t.id)
        glTexBuffer(t.texturetype, t.internalformat, t.buffer.id)
        glBindTexture(t.texturetype, 0)
    end
    isdone
end


text:
######### CmdStan program example  ###########

using CmdStan, Test, Statistics

ProjDir = dirname(@__FILE__)
cd(ProjDir) do

  bernoullimodel = "
  data {
    int<lower=1> N;
    int<lower=0,upper=1> y[N];
  }
  parameters {
    real<lower=0,upper=1> theta;
  }
  model {
    theta ~ beta(1,1);
    y ~ bernoulli(theta);
  }
  "

  bernoullidata = Dict("N" => 1, "y" => [0])

  global stanmodel, rc, sim

  stanmodel = Stanmodel(num_samples=1200, thin=2, name="bernoulli",
    output_format=:array, model=bernoullimodel);

  rc, sim, cnames = stan(stanmodel, bernoullidata, ProjDir,
    diagnostics=false, CmdStanDir=CMDSTAN_HOME);

  if rc == 0
    println()
    println("Test round.(mean(theta), 1) ≈ 0.3 rtol=0.1")
    @test round.(mean(sim[:,8,:]), digits=1) ≈ 0.3 rtol=0.1
  end
end # cd


text:
# plot gaussian distribution model

module GaussDistModel
    using Plots
    pyplot()

    function p(z, mu=209.7, dev=23.4)
        return exp(-(z - mu)^2 / (2 * dev)) / sqrt(2 * pi * dev)
    end

    function main(is_test=false)
        zs = Array(190:230)
        ys = [p(z) for z in zs]

        println("Observations: $(zs)")
        println("Probability: $(ys)")

        plot(zs, ys, label="gaussian distribution")

        if is_test == false
            save_path = joinpath(split(@__FILE__, "src")[1], "src/prob_stats/prob_model/gauss_dist_model/gauss_dist_model.png")
            savefig(save_path)
        end
    end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
fname = joinpath(dirname(@__FILE__), "..", "data", "test_Hz19.5-testing.bdf")

a = read_SSR(fname)

b = deepcopy(a)
b.data = rand(size(b.data))

a = extract_epochs(a)
original = a.processing["epochs"]

@test size(a.processing["epochs"]) == (8388,28,6)

a = reduce_epochs(a, 9)
@test size(a.processing["epochs"]) == (8388,9,6)

@test mean(original[:, 1:3, :], 2) == a.processing["epochs"][:, 1, :]
@test mean(original[:, 4:6, :], 2) == a.processing["epochs"][:, 2, :]
@test mean(original[:, 7:9, :], 2) == a.processing["epochs"][:, 3, :]
@test mean(original[:, 25:28, :], 2) == a.processing["epochs"][:, 9, :]

#
# Fake a leadfield
#

a = extract_epochs(a)
b = extract_epochs(b)

x = repmat(collect(-2:2.0), 25)
y = repmat(vec(ones(5) * collect(-2:2.0)'), 5)
z = vec(ones(5*5) * collect(-2:2.0)')
t = [1.0]
H = rand(125, 3, 6)
L = Leadfield(H, x, y, z, channelnames(a))

v = beamformer_lcmv(a, b, L)
v = beamformer_lcmv(a, b, L, subspace = 0.0,  regularisation = 0.0)
v = beamformer_lcmv(a, b, L, subspace = 0.0,  regularisation = 0.001)
v = beamformer_lcmv(a, b, L, subspace = 0.9,  regularisation = 0.001)
v = beamformer_lcmv(a, b, L, subspace = 0.99, regularisation = 0.001)

show(v)

cpsd = cross_spectral_density(a)
@test isa(cpsd, Array{Complex{Float64}, 2})

# Test bilateral correction
v = correct_midline(v)
@test isa(v, EEG.VolumeImage)


text:
"""
Perform EPANET hydraulic simulation (via WNTR) and compute timeseries of flows and heads
"""

function _get_wntr_node_attribute(wntr_simulation, name::AbstractString, attribute::String)
    values = PyCall.getproperty(wntr_simulation.node[attribute], name).values[1:end-1]
    return Float64.(values) # Convert to 64-bit floating point array.
end

function _get_wntr_link_attribute(wntr_simulation, name::AbstractString, attribute::String)
    values = PyCall.getproperty(wntr_simulation.link[attribute], name).values[1:end-1]
    return Float64.(values) # Convert to 64-bit floating point array.
end

# compare nodes
function get_node_dataframe(wm_data, wm_solution, wntr_data, wntr_simulation, node_id)
    num_time_step = length(wm_solution["nw"])            # number of time steps
    time_step = wm_data["time_step"] / 3600.0                        # length per time step (hour)
    node_name = node_id

    head_wntr = _get_wntr_node_attribute(wntr_simulation, node_name, "head")
    pressure_wntr = _get_wntr_node_attribute(wntr_simulation, node_name, "pressure")
    elevation = Array{Float64,1}(undef,num_time_step)
    head_watermodels = Array{Float64,1}(undef,num_time_step)
    pressure_watermodels = Array{Float64,1}(undef,num_time_step)

    for t in 1:num_time_step
        elevation[t] = wm_data["nw"][string(t)]["node"][node_id]["elevation"]
        head_watermodels[t] = wm_solution["nw"][string(t)]["node"][node_id]["h"]
        pressure_watermodels[t] = wm_solution["nw"][string(t)]["node"][node_id]["p"]
    end

    node_df = DataFrames.DataFrame(time = 1:time_step:time_step*num_time_step, elevation = elevation,
        head_wntr = head_wntr, head_watermodels = head_watermodels,
        pressure_wntr = pressure_wntr, pressure_watermodels = pressure_watermodels)

    return node_df
end

# compare tanks
function get_tank_dataframe(wm_data,wm_solution,wntr_data,wntr_simulation, tank_id)
    num_time_step = length(wm_solution["nw"])            # number of time steps
    time_step = wm_data["time_step"]/3600                            # length per time step (hour)
    tank_node_id = string(wm_data["nw"]["1"]["tank"][tank_id]["node"])
    tank_name = tank_node_id
    diameter = wntr_data.nodes._data[tank_name].diameter
    volume_wntr = Array{Float64,1}(undef,num_time_step)
    volume_watermodels = Array{Float64,1}(undef,num_time_step)
    level_wntr = Array{Float64,1}(undef,num_time_step)
    level_watermodels = Array{Float64,1}(undef,num_time_step)

    volume_wntr = _get_wntr_node_attribute(wntr_simulation, tank_name, "pressure") .* (0.25 * pi * diameter^2)
    level_wntr = _get_wntr_node_attribute(wntr_simulation, tank_name, "pressure")

    for t in 1:num_time_step
        volume_watermodels[t] = wm_solution["nw"][string(t)]["tank"][tank_id]["V"]
        level_watermodels[t] = wm_solution["nw"][string(t)]["node"][tank_node_id]["p"]
    end

    tank_df = DataFrames.DataFrame(time = 1:time_step:time_step*num_time_step, volume_wntr = volume_wntr, volume_watermodels = volume_watermodels,
        level_wntr = level_wntr, level_watermodels = level_watermodels)

    return tank_df
end

# compare pipes
function get_pipe_dataframe(wm_data,wm_solution,wntr_data,wntr_simulation,pipe_id)
    num_time_step = length(wm_solution["nw"])            # number of time steps
    time_step = wm_data["time_step"]/3600                            # length per time step (hour)
    pipe_name = "pipe"*pipe_id
    flow_wntr = Array{Float64,1}(undef,num_time_step)
    flow_watermodels = Array{Float64,1}(undef,num_time_step)
    head_loss_wntr = Array{Float64,1}(undef,num_time_step)
    head_loss_watermodels = Array{Float64,1}(undef,num_time_step)

    node_fr = string(wm_data["nw"]["1"]["pipe"][pipe_id]["node_fr"])
    node_to = string(wm_data["nw"]["1"]["pipe"][pipe_id]["node_to"])
    flow_wntr = _get_wntr_link_attribute(wntr_simulation, pipe_name, "flowrate")
    head_start_wntr = _get_wntr_node_attribute(wntr_simulation, node_fr, "head")
    head_end_wntr = _get_wntr_node_attribute(wntr_simulation, node_to, "head")
    head_loss_wntr = head_start_wntr .- head_end_wntr

    for t in 1:num_time_step
        flow_watermodels[t] = wm_solution["nw"][string(t)]["pipe"][pipe_id]["q"]
        head_loss_watermodels[t] = wm_solution["nw"][string(t)]["node"][node_fr]["h"]-wm_solution["nw"][string(t)]["node"][node_to]["h"]
    end

    pipe_df = DataFrames.DataFrame(time = 1:time_step:time_step*num_time_step, flow_wntr = flow_wntr, flow_watermodels = flow_watermodels,
        head_loss_wntr = head_loss_wntr, head_loss_watermodels = head_loss_watermodels)

    return pipe_df
end

# compare short pipes
function get_short_pipe_dataframe(wm_data,wm_solution,wntr_data,wntr_simulation,short_pipe_id)
    num_time_step = length(wm_solution["nw"])            # number of time steps
    time_step = wm_data["time_step"]/3600                            # length per time step (hour)
    short_pipe_name = "short_pipe"*short_pipe_id
    flow_wntr = Array{Float64,1}(undef,num_time_step)
    flow_watermodels = Array{Float64,1}(undef,num_time_step)
    head_loss_wntr = Array{Float64,1}(undef,num_time_step)
    head_loss_watermodels = Array{Float64,1}(undef,num_time_step)

    node_fr = string(wm_data["nw"]["1"]["short_pipe"][short_pipe_id]["node_fr"])
    node_to = string(wm_data["nw"]["1"]["short_pipe"][short_pipe_id]["node_to"])
    flow_wntr = _get_wntr_link_attribute(wntr_simulation, short_pipe_name, "flowrate")
    head_start_wntr = _get_wntr_node_attribute(wntr_simulation, node_fr, "head")
    head_end_wntr = _get_wntr_node_attribute(wntr_simulation, node_to, "head")
    head_loss_wntr = head_start_wntr .- head_end_wntr

    for t in 1:num_time_step
        flow_watermodels[t] = wm_solution["nw"][string(t)]["short_pipe"][short_pipe_id]["q"]
        head_loss_watermodels[t] = wm_solution["nw"][string(t)]["node"][node_fr]["h"]-wm_solution["nw"][string(t)]["node"][node_to]["h"]
    end

    short_pipe_df = DataFrames.DataFrame(time = 1:time_step:time_step*num_time_step, flow_wntr = flow_wntr, flow_watermodels = flow_watermodels,
        head_loss_wntr = head_loss_wntr, head_loss_watermodels = head_loss_watermodels)

    return short_pipe_df
end

# compare valves
function get_valve_dataframe(wm_data,wm_solution,wntr_data,wntr_simulation,valve_id)
    num_time_step = length(wm_solution["nw"])            # number of time steps
    time_step = wm_data["time_step"]/3600                            # length per time step (hour)
    valve_name = "valve"*valve_id
    status_watermodels = Array{Float64,1}(undef,num_time_step)
    flow_wntr = Array{Float64,1}(undef,num_time_step)
    flow_watermodels = Array{Float64,1}(undef,num_time_step)
    head_loss_wntr = Array{Float64,1}(undef,num_time_step)
    head_loss_watermodels = Array{Float64,1}(undef,num_time_step)

    # I presume this is here for development and need not be used now? JJS 5/10/21
    #println("--------- check all valves -----------")
    #println(keys(wm_data["nw"]["1"]["valve"]))

    node_fr = string(wm_data["nw"]["1"]["valve"][valve_id]["node_fr"])
    node_to = string(wm_data["nw"]["1"]["valve"][valve_id]["node_to"])
    status_wntr = _get_wntr_link_attribute(wntr_simulation, valve_name, "status")
    flow_wntr = _get_wntr_link_attribute(wntr_simulation, valve_name, "flowrate")
    head_start_wntr = _get_wntr_node_attribute(wntr_simulation, node_fr, "head")
    head_end_wntr = _get_wntr_node_attribute(wntr_simulation, node_to, "head")
    head_loss_wntr = head_start_wntr .- head_end_wntr

    for t in 1:num_time_step
        status_watermodels[t] = wm_solution["nw"][string(t)]["valve"][valve_id]["status"]
        flow_watermodels[t] = wm_solution["nw"][string(t)]["valve"][valve_id]["q"]
        head_loss_watermodels[t] = wm_solution["nw"][string(t)]["node"][node_fr]["h"]-wm_solution["nw"][string(t)]["node"][node_to]["h"]
    end

    valve_df = DataFrames.DataFrame(time = 1:time_step:time_step*num_time_step, status_wntr = status_wntr, status_watermodels = status_watermodels,
        flow_wntr = flow_wntr, flow_watermodels = flow_watermodels, head_loss_wntr = head_loss_wntr, head_loss_watermodels = head_loss_watermodels)

    return valve_df
end

# compare regulators
# TODO

# compare pumps
function get_pump_dataframe(wm_data,wm_solution,wntr_data,wntr_simulation,pump_id)
    num_time_step = length(wm_solution["nw"])            # number of time steps
    time_step = wm_data["time_step"]/3600                            # length per time step (hour)
    pump_name = "pump"*pump_id
    pump_obj = wntr_data.get_link(pump_name)
    status_watermodels = Array{Float64,1}(undef,num_time_step)
    flow_watermodels = Array{Float64,1}(undef,num_time_step)
    head_gain_watermodels = Array{Float64,1}(undef,num_time_step)
    power_wntr = Array{Float64,1}(undef,num_time_step)               # Watt
    power_watermodels = Array{Float64,1}(undef,num_time_step)        # Watt
    cost_wntr = Array{Float64,1}(undef,num_time_step)                # $
    cost_watermodels = Array{Float64,1}(undef,num_time_step)         # $
    node_fr = string(wm_data["nw"]["1"]["pump"][pump_id]["node_fr"])
    node_to = string(wm_data["nw"]["1"]["pump"][pump_id]["node_to"])

    function compute_pump_power(pump_obj,q,dh,wntr_data)
        if pump_obj.efficiency == nothing
            eff_curve = nothing
        else
            eff_curve = wntr_data.get_curve(pump_obj.efficiency).points
        end
        if eff_curve != nothing
            for j in 1:length(eff_curve)
                if (j == 1) & (q < eff_curve[j][1])
                    eff = eff_curve[j][2]/100 # use the first Q-eff point when q < q_min on the eff curve
                    break
                elseif j == length(eff_curve)
                    eff = eff_curve[j][2]/100 # use the last Q-eff point when q > q_max on the eff curve
                elseif (q >= eff_curve[j][1]) & (q < eff_curve[j+1][1])
                    # linear interpolation of the efficiency curve
                    a1 = (eff_curve[j+1][1]-q)/(eff_curve[j+1][1]-eff_curve[j][1])
                    a2 = (q-eff_curve[j][1])/(eff_curve[j+1][1]-eff_curve[j][1])
                    eff = (a1*eff_curve[j][2]+a2*eff_curve[j+1][2])/100
                    break
                end
            end
        elseif wntr_data.options.energy.global_efficiency != nothing
            eff = wntr_data.options.energy.global_efficiency/100 # e.g. convert 60 to 60% (0.6)
        else
            Memento.info(_LOGGER, "No pump efficiency provided, 75% is used")
            eff = 0.75
        end

        power = _WM._GRAVITY * _WM._DENSITY * dh * q / eff # Watt
        return power
    end

    status_wntr = _get_wntr_link_attribute(wntr_simulation, pump_name, "status")
    flow_wntr = _get_wntr_link_attribute(wntr_simulation, pump_name, "flowrate")
    head_start_wntr = _get_wntr_node_attribute(wntr_simulation, node_fr, "head")
    head_end_wntr = _get_wntr_node_attribute(wntr_simulation, node_to, "head")
    head_gain_wntr = status_wntr .* (head_end_wntr .- head_start_wntr)

    for t in 1:num_time_step
        status_watermodels[t] = wm_solution["nw"][string(t)]["pump"][pump_id]["status"]
        flow_watermodels[t] = wm_solution["nw"][string(t)]["pump"][pump_id]["q"]
        head_gain_watermodels[t] = round(status_watermodels[t])*(wm_solution["nw"][string(t)]["node"][node_to]["h"]-wm_solution["nw"][string(t)]["node"][node_fr]["h"])
        power_wntr[t] = status_wntr[t]*compute_pump_power(pump_obj,flow_wntr[t],head_gain_wntr[t],wntr_data)
        power_watermodels[t] = status_watermodels[t]*compute_pump_power(pump_obj,flow_watermodels[t],head_gain_watermodels[t],wntr_data)
        energy_price = wm_data["nw"][string(t)]["pump"][pump_id]["energy_price"]*3600    # $/Wh
        cost_wntr[t] = power_wntr[t]*time_step*energy_price
        cost_watermodels[t] = power_watermodels[t]*time_step*energy_price
    end

    pump_df = DataFrames.DataFrame(time = 1:time_step:time_step*num_time_step, status_wntr = status_wntr, status_watermodels = status_watermodels,
        flow_wntr = flow_wntr, flow_watermodels = flow_watermodels,
        head_gain_wntr = head_gain_wntr, head_gain_watermodels = head_gain_watermodels,
        power_wntr = power_wntr, power_watermodels = power_watermodels,
        cost_wntr = cost_wntr, cost_watermodels = cost_watermodels)

    return pump_df
end


text:
######### Stan optimize example  ###########

using StanOptimize, Test

bernoulli_model = "
data {
  int<lower=1> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}
";

bernoulli_data = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])

stanmodel = OptimizeModel("optimize_1",  bernoulli_model);

(sample_path, log_path) = stan_sample(stanmodel, data=bernoulli_data);

if sample_path !== Nothing
  optim1, cnames = read_optimize(stanmodel)

  @test optim1["theta"][end] ≈ 0.3 atol=0.001
end

# Same with saved iterations
sm = OptimizeModel("optimize_2", bernoulli_model;
  method = StanOptimize.Optimize(save_iterations = true));

(sample_path, log_path)  = stan_optimize(sm, data=bernoulli_data);

if sample_path !== Nothing
  optim2, cnames = read_optimize(sm)
  @test optim2["theta"][end] ≈ 0.3 atol=0.001
end


text:
include("./Distribution.jl")
using LinearAlgebra

mutable struct Optimisation_problem
    #### This small structure is only to make Newton's method mor readable ####
    δ::Array{Float64, 1} # describes the step
    cost::Float64 # value of the cost function
end

function make_a_step!(dist::Distribution, prb::Optimisation_problem)
    #### modifies β, γ, M -- assumes the step δ is correct ####
    # getting the constants we will need
    n_pop = size(dist.J)[1]
    lmax_half = size(dist.J)[3]

    # let's perform the step
    dist.β -= prb.δ[1]
    dist.γ -= prb.δ[2]
    for l_half = 1:lmax_half
        for a = 1:n_pop
            line = (l_half - 1) * n_pop + a + 2
            dist.M[line - 2] -= prb.δ[line]
        end
    end
end

function update_Newton!(dist::Distribution, prb::Optimisation_problem)
    #### Compute the step of Newton's method, then makes it                                    ####
    #### Assumes that cost and jacobian have already been computed when we enter this function ####

    ldiv!(prb.δ, lu(dist.jacobian), dist.cost) # this computes the step to performs

    make_a_step!(dist, prb) # this makes the step and modifies β, γ, M

    Cost_and_Jacobian!(dist) # this is makes sure the cost and jacobian have been computed when we next call update_Newton!
    prb.cost = norm(dist.cost) # computes how close to the root we are AFTER the step
end

function Root_finder_Newton!(dist::Distribution)
    #### Performs the entire Newton method                                                                    ####
    #### However, it doesn't choose the initial point itself, this should be done by a wrapper to this script ####
    # getting the constants we will need
    n_pop = size(dist.J)[1]
    lmax_half = size(dist.J)[3]

    # preallocating an array we will need
    δ = zeros(lmax_half * n_pop + 2)

    # Let's choose a precision goal and an maximal number of iterations
    ϵ = 10^(-6) # maximum acceptable cost
    n_max = 100 # maximumm acceptable number of iterations

    # Computing the initial cost because the Newton step assumes that when the step begins the cost has already been computed
    Cost_and_Jacobian!(dist)
    cost = norm(dist.cost)

    # Initializing our small Optimization problem with this cost
    prb = Optimisation_problem(δ, cost)

    # Let's perform the Newton method now
    n = 0
    while ((n < n_max) && (prb.cost > ϵ))
        n += 1
        update_Newton!(dist, prb)
    end
end

max_stars_repo_name: NathanMagnan/VRROOMe

max_issues_repo_name: NathanMagnan/VRROOMe

max_forks_repo_name: NathanMagnan/VRROOMe


text:
### A Pluto.jl notebook ###
# v0.12.6

using Markdown
using InteractiveUtils

# ╔═╡ acf67f82-1a85-11eb-2bdf-abc05bfb68eb
begin
	using DataFrames
	using Statistics
	using CSV
end

# ╔═╡ 8061b6ba-1a85-11eb-1a73-130f152eb8d8
VERSION

# ╔═╡ ad304f82-1a85-11eb-1b1e-bba3c3d49c45
isfile.(["Project.toml", "Manifest.toml"])

# ╔═╡ 8849a66c-1adc-11eb-3c30-371f47ba053e
md"""
# Working with text files
"""

# ╔═╡ 8609a440-1ad2-11eb-0696-09f4c5a8fc0f
begin
	url = "https://archive.ics.uci.edu/ml/machine-learning-databases/auto-mpg/auto-mpg.data-original"
	download(url, "auto.txt")
end

# ╔═╡ 8623433e-1ad2-11eb-2589-79e06429f720
readlines("auto.txt")

# ╔═╡ 2b3bd64c-1add-11eb-12a7-dda4b1d43b59
md"""
## Data Cleaning
- Table has no headers
- Last column is tab separated, while earlier columns are separated with varying number of spaces
- Missing values are encoded by "NA"
"""

# ╔═╡ a7f25360-1ac8-11eb-1c1f-552f8c1a9d8a
raw_str = read("auto.txt", String)

# ╔═╡ 96c5bac2-1add-11eb-3852-4d3153fccd80
md"""
Replace all tabs with spaces
"""

# ╔═╡ 95d8036a-1add-11eb-35f4-770cc195ab45
str_no_tab = replace(raw_str, '\t' => ' ')

# ╔═╡ be35efd4-1add-11eb-1652-8bb45ea07ab1
io = IOBuffer(str_no_tab)

# ╔═╡ 6c91e6f4-1add-11eb-30b8-0711ba80a788
md"""
We can think of IO as n in-memory I/O stream. Therefore, we can pass this stream to CSV. While loading our data we are are making the following assumptions.

1. Use spacers as delimeters.
2. Ignore repeated (consecutive) occurrences of the delimiter.
3. We explicitly pass column names via header keyword argument.
4. We specify that missing values are represented using "NA"
"""

# ╔═╡ 0fc39f28-1b01-11eb-172c-c73f51a11353
sin(1.3)

# ╔═╡ 1622caa8-1b01-11eb-2c94-6ddeb98ee29e
1.2 |> sin

# ╔═╡ 6cb070c4-1add-11eb-22c4-91ccd53e92f6
df1 = CSV.File(
	io,
	delim=' ',
	ignorerepeated=true,
	header=[:mpg, :cylinders, :displacement, :horsepower,
            :weight, :acceleration, :year, :origin, :name],
	missingstring="NA") |>
	DataFrame

# ╔═╡ 6ccb47d2-1add-11eb-20f2-a5043f12c160
df1

# ╔═╡ 6ce76124-1add-11eb-29dc-af32921df293
ENV["COLUMNS"], ENV["LINES"] = 200, 15

# ╔═╡ 6d02368e-1add-11eb-252d-bdb11b1c3e41
df1

# ╔═╡ 6d1d1652-1add-11eb-0f14-e541320df058
df_raw = CSV.File("auto.txt", header=[:metrics, :name]) |> DataFrame

# ╔═╡ 6d37f116-1add-11eb-1a4a-fbd3da64ef4c
str_metrics = split.(df_raw.metrics)

# ╔═╡ fc665566-1afc-11eb-1bf5-ebff28a03169
df1_2 = DataFrame(fill(Float64, 8), [:mpg, :cylinders, :displacement, :horsepower, :weight, :acceleration, :year, :origin])

# ╔═╡ fc4dba10-1afc-11eb-32b6-e5eff96c251f
allowmissing!(df1_2, [:mpg, :horsepower])

# ╔═╡ fc3001fa-1afc-11eb-3f59-a779618967c0
for row in str_metrics
	push!(df1_2, [v == "NA" ? missing : parse(Float64, v) for v in row])
end

# ╔═╡ fc1655ca-1afc-11eb-1a6d-2de5383414fd
df1_2

# ╔═╡ a30237b4-1afd-11eb-2a46-4dd62d19b907
df1_2.name = df_raw.name

# ╔═╡ fbfa3032-1afc-11eb-2f21-9f3afc8a66ef
df1_2

# ╔═╡ fbdedf8c-1afc-11eb-33f2-49a07fb38b1e
df1_2.name == df_raw.name

# ╔═╡ fbbeeeaa-1afc-11eb-0cb9-bd01a64ef604
df1_2.name = df_raw[:, :name]

# ╔═╡ fba1329a-1afc-11eb-37a0-55342cda66c9
isequal(df1_2, df1)

# ╔═╡ 6d540338-1add-11eb-0442-67bbdfeba0fb
sum(count(ismissing, col) for col in eachcol(df1))

# ╔═╡ 6d701668-1add-11eb-1494-8ddd04b56334
mapcols(x -> count(ismissing, x), df1)

# ╔═╡ 6ddb83bc-1add-11eb-1e8c-752a464a1ba4
filter(row -> any(ismissing, row), df1)

# ╔═╡ 6df796f6-1add-11eb-2cf8-557ff3745852
df1.brand = first.(split.(df1.name))

# ╔═╡ 6e1b0dca-1add-11eb-1d1a-23b292fe2b3c
df2 = dropmissing(df1)

# ╔═╡ a7d388c2-1ac8-11eb-1bd7-59a0654fd402
df2[df2.brand .== "saab", :]

# ╔═╡ c5824836-1aff-11eb-21cf-2d5e2c4c58b9
filter(:brand => ==("saab"), df2)

# ╔═╡ c565d67e-1aff-11eb-1e6d-2d137f119d62
filter(row -> row.brand == "saab", df2)

# ╔═╡ c54c47a4-1aff-11eb-15d4-97cbec4614d5
CSV.write("auto2.csv", df2)

# ╔═╡ c52fda08-1aff-11eb-1635-7769f5e1fb29
readlines("auto2.csv")

# ╔═╡ c5110f36-1aff-11eb-2712-e56fcfd44535

# ╔═╡ c4f62900-1aff-11eb-319c-6fbed55179ff

# ╔═╡ a79492ea-1ac8-11eb-098c-f12b9d38f313

# ╔═╡ a76f13ea-1ac8-11eb-3800-b1c27e74a466

# ╔═╡ a76f037a-1ac8-11eb-199c-cb332db6cca1

# ╔═╡ a7282752-1ac8-11eb-0f05-83cace66d3b0

# ╔═╡ ab7c9592-1a85-11eb-342f-3dab4c1d54a6

# ╔═╡ Cell order:
# ╠═8061b6ba-1a85-11eb-1a73-130f152eb8d8
# ╠═ad304f82-1a85-11eb-1b1e-bba3c3d49c45
# ╠═acf67f82-1a85-11eb-2bdf-abc05bfb68eb
# ╟─8849a66c-1adc-11eb-3c30-371f47ba053e
# ╠═8609a440-1ad2-11eb-0696-09f4c5a8fc0f
# ╠═8623433e-1ad2-11eb-2589-79e06429f720
# ╟─2b3bd64c-1add-11eb-12a7-dda4b1d43b59
# ╠═a7f25360-1ac8-11eb-1c1f-552f8c1a9d8a
# ╟─96c5bac2-1add-11eb-3852-4d3153fccd80
# ╠═95d8036a-1add-11eb-35f4-770cc195ab45
# ╠═be35efd4-1add-11eb-1652-8bb45ea07ab1
# ╟─6c91e6f4-1add-11eb-30b8-0711ba80a788
# ╠═0fc39f28-1b01-11eb-172c-c73f51a11353
# ╠═1622caa8-1b01-11eb-2c94-6ddeb98ee29e
# ╠═6cb070c4-1add-11eb-22c4-91ccd53e92f6
# ╠═6ccb47d2-1add-11eb-20f2-a5043f12c160
# ╠═6ce76124-1add-11eb-29dc-af32921df293
# ╠═6d02368e-1add-11eb-252d-bdb11b1c3e41
# ╠═6d1d1652-1add-11eb-0f14-e541320df058
# ╠═6d37f116-1add-11eb-1a4a-fbd3da64ef4c
# ╠═fc665566-1afc-11eb-1bf5-ebff28a03169
# ╠═fc4dba10-1afc-11eb-32b6-e5eff96c251f
# ╠═fc3001fa-1afc-11eb-3f59-a779618967c0
# ╠═fc1655ca-1afc-11eb-1a6d-2de5383414fd
# ╠═a30237b4-1afd-11eb-2a46-4dd62d19b907
# ╠═fbfa3032-1afc-11eb-2f21-9f3afc8a66ef
# ╠═fbdedf8c-1afc-11eb-33f2-49a07fb38b1e
# ╠═fbbeeeaa-1afc-11eb-0cb9-bd01a64ef604
# ╠═fba1329a-1afc-11eb-37a0-55342cda66c9
# ╠═6d540338-1add-11eb-0442-67bbdfeba0fb
# ╠═6d701668-1add-11eb-1494-8ddd04b56334
# ╠═6ddb83bc-1add-11eb-1e8c-752a464a1ba4
# ╠═6df796f6-1add-11eb-2cf8-557ff3745852
# ╠═6e1b0dca-1add-11eb-1d1a-23b292fe2b3c
# ╠═a7d388c2-1ac8-11eb-1bd7-59a0654fd402
# ╠═c5824836-1aff-11eb-21cf-2d5e2c4c58b9
# ╠═c565d67e-1aff-11eb-1e6d-2d137f119d62
# ╠═c54c47a4-1aff-11eb-15d4-97cbec4614d5
# ╠═c52fda08-1aff-11eb-1635-7769f5e1fb29
# ╠═c5110f36-1aff-11eb-2712-e56fcfd44535
# ╠═c4f62900-1aff-11eb-319c-6fbed55179ff
# ╠═a79492ea-1ac8-11eb-098c-f12b9d38f313
# ╠═a76f13ea-1ac8-11eb-3800-b1c27e74a466
# ╠═a76f037a-1ac8-11eb-199c-cb332db6cca1
# ╠═a7282752-1ac8-11eb-0f05-83cace66d3b0
# ╠═ab7c9592-1a85-11eb-342f-3dab4c1d54a6

max_stars_repo_name: martin-fabbri/colab-notebooks

max_issues_repo_name: martin-fabbri/colab-notebooks

max_forks_repo_name: martin-fabbri/colab-notebooks


text:
function write_key(file_stream::IOStream, key_name::String, key_value)

	@printf(file_stream, "%s", "\u0009")
	@printf(file_stream, "%s", join(("<", key_name, ">")))

	if typeof(key_value) == Float64
		@printf(file_stream, "%s", key_value)
	elseif typeof(key_value) == Array{Float64, 1}
		for i = 1:length(key_value)-1
			@printf(file_stream, "%s", key_value[i])
			@printf(file_stream, "%s", ",")
		end
		@printf(file_stream, "%s", key_value[end])
	elseif typeof(key_value) == Int64
		@printf(file_stream, "%d", key_value)
	elseif typeof(key_value) == String
		@printf(file_stream, "%s", key_value)
	else
		error("Incompatible data type.")
	end

	@printf(file_stream, "%s", join(("</", key_name, ">")))
	@printf(file_stream, "%s", "\n")

	nothing
end


text:
using MixedModeDebugger
using Test

@testset "MixedModeDebugger.jl" begin
    @testset "$file" for file in ("overdub.jl",)
        include(file)
    end
end


text:
# added structure for Bloch Hamiltonian. Instead of storing  a list of vectors, an in-place function is stored.

"""
    BlochHamiltonian{D, F<:Function}

Stores a tight-binding system in Bloch space.

### Fields
The tight-binsing system consists of:
- 'latbasis::SMatrix{D, D, Float64}': the latbasis for a `D`-dimensional Bravais lattice. The basis vectors are stored as columns of the matrix basis, i. e. a1 = basis[:, 1]
- 'orbitals::Vector{Orbital}': a list of orbitals within the unit cell
- 'hk!::F': The Bloch Hamiltonian: a function that for a given `k` acts on a matrix `mat` and stores there H(k)
"""
struct BlochHamiltonian{D, F}
    basis::SMatrix{D, D, Float64}
    orbitals::Vector{Orbital}
    hk!::F
end

"""
    BlochHamiltonian_noF{D}

Stores a tight-binding system in Bloch space.

### Fields
The tight-binsing system consists of:
- 'latbasis::SMatrix{D, D, Float64}': the latbasis for a `D`-dimensional Bravais lattice. The basis vectors are stored as columns of the matrix basis, i. e. a1 = basis[:, 1]
- 'orbitals::Vector{Orbital}': a list of orbitals within the unit cell
- 'hk!::F': The Bloch Hamiltonian: a function that for a given `k` acts on a matrix `mat` and stores there H(k)
"""
struct BlochHamiltonian_noF{D}
    basis::SMatrix{D, D, Float64}
    orbitals::Vector{Orbital}
    hk!::Function
end