text:
"""
`module  JAC.ManyElectron`
    ... a submodel of JAC that contains all structs and methods for many-electron CSF, basis and levels.
"""
module ManyElectron

    using  ..Basics, ..Defaults,  ..Radial
    export Configuration, ConfigurationR, AsfSettings, CsfR, Basis, Level, Multiplet,
           AbstractQedModel, NoneQed, QedPetersburg, QedSydney

    """
    `abstract type ManyElectron.AbstractQedModel`
        ... defines an abstract and a number of singleton types for dealing with QED estimates in the many-electron
            computations.

      + struct QedPetersburg
        ... to estimate the QED corrections due to the model Hamiltonian by Shabaev and coworkers (2013).

      + struct QedSydney
        ... to estimate the QED corrections by means of the radiative potential by Flambaum and Ginges .
    """
    abstract type  AbstractQedModel                      end
    struct         QedPetersburg  <:  AbstractQedModel   end
    struct         QedSydney      <:  AbstractQedModel   end
    struct         NoneQed        <:  AbstractQedModel   end

    """
    `struct  ManyElectron.Configuration`
        ... defines a struct for a non-relativistic electron configuration that is fully speficied by its shell notations
            (such as "1s", "2s", "2p", ....) and the corresponding occupation numbers (>= 0). An electron configuration
            is independent of any order of the shells, and a zero occupation is assumed for all shells that do not appear as a
            key in the shell (dictionary). Therefore, the number of keys do not allow any conclusion about the underlying orbital
            space of any considered computation that include more than just a single configuration.

        + shells         ::Dict{Shell,Int64}   ... Dictionary that maps shells to their occupation.
        + NoElectrons    ::Int64               ... No. of electrons.
    """
    struct  Configuration
        shells           ::Dict{Shell,Int64}
        NoElectrons      ::Int64
    end

    """
    `ManyElectron.Configuration(sa::String)`
        ... constructor for a given configuration string, such as "[He]", "[Ne]", "[Ne] 3s 3p^6"  or "1s 2p^6 3s^2 3p".
    """
    function Configuration(sa::String)
        sax = strip(sa)
        shellOccList = split(sax, " ")
        ##x println("shellOccList = $shellOccList")
        NoElectrons = 0;    shellDict    = Dict{Shell,Int64}()
        for shocc in shellOccList
            if       shocc == " "    continue
            elseif   shocc  in [ "[He]", "[Ne]", "[Ar]", "[Kr]", "[Xe]" ]
                shellDict = Base.merge( shellDict, Dict( Shell("1s") => 2));                                       NoElectrons = NoElectrons + 2
                shocc == "[He]"  &&  continue
                shellDict = Base.merge( shellDict, Dict(Shell("2s") => 2,  Shell("2p") => 6));                     NoElectrons = NoElectrons + 8
                shocc == "[Ne]"  &&  continue
                shellDict = Base.merge( shellDict, Dict(Shell("3s") => 2,  Shell("3p") => 6));                     NoElectrons = NoElectrons + 8
                shocc == "[Ar]"  &&  continue
                shellDict = Base.merge( shellDict, Dict(Shell("3d") => 10, Shell("4s") => 2, Shell("4p") => 6));   NoElectrons = NoElectrons + 18
                shocc == "[Kr]"  &&  continue
                shellDict = Base.merge( shellDict, Dict(Shell("5s") => 10, Shell("5s") => 2, Shell("5p") => 6));   NoElectrons = NoElectrons + 18
                shocc == "[Xe]"  &&  continue
            else
                ia = findall(in("^"), shocc);
                if     length(ia) == 0     shellDict = Base.merge( shellDict, Dict( Shell(shocc) => 1));           NoElectrons = NoElectrons + 1
                else   occ = parse(Int64, shocc[ia[1]+1:end] );    sh = Shell( shocc[1:ia[1]-1] )
                       shellDict = Base.merge( shellDict, Dict( sh => occ));                                       NoElectrons = NoElectrons + occ
                end
            end
        end
        return( Configuration(shellDict, NoElectrons) )
    end

    # `Base.show(io::IO, conf::Configuration)`  ... prepares a proper printout of the non-relativistic configuration conf::Configuration.
    function Base.show(io::IO, conf::Configuration)
        sa = "Configuration: " * string(conf)
        print(io, sa)
    end

    # `Base.string(conf::Configuration)`  ... provides a String notation for the variable conf::Configuration.
    function Base.string(conf::Configuration)
        wa = keys(conf.shells);   va = values(conf.shells)
        wb = Defaults.getDefaults("ordered shell list: non-relativistic", 7)

        sa = ""
        for  k  in  wb
            if   k  in  wa
              occ = conf.shells[k]
              sa = sa * string(k) * "^$occ "
            end
        end

        return( sa )
    end

    """
    `Base.:(==)(confa::Configuration, confb::Configuration)`
        ... compares (recursively) two non-relativistic configurations and return true if all subfields are equal,
            and false otherwise
    """
    function  Base.:(==)(confa::Configuration, confb::Configuration)
        if   confa.NoElectrons  != confb.NoElectrons    return( false )    end
        wk = keys(confa.shells)
        for  k in wk
            if  !haskey(confb.shells, k)  ||   confa.shells[k] != confb.shells[k]    return( false )    end
        end
        return( true )
    end

    """
    `struct  ManyElectron.ConfigurationR`
        ... defines a type for a relativistic electron configuration that is fully speficied by its subshell notations
            (such as "1s_1/2", "2s_1/2", "2p_3/2", ....) and the corresponding occupation numbers (>= 0). An electron
            configuration is independent of any order of the subshells, and a zero occupation is assumed for all shells that
            do not appears as a key. Therefore, the number of keys do not allow any conclusion about the underlying orbital
            space of any considered computation that include more than a single configuration.

        + subshells      ::Dict{Subshell,Int64}  ... Dictionary that maps subshells to their occupation.
        + NoElectrons    ::Int64                 ... No. of electrons.
    """
    struct  ConfigurationR
        subshells        ::Dict{Subshell,Int64}
        NoElectrons      ::Int64
    end

    # `Base.show(io::IO, conf::ConfigurationR)`  ... prepares a proper printout of the relativistic configuration conf::ConfigurationR.
    function Base.show(io::IO, conf::ConfigurationR)
        sa = "Configuration: " * string(conf)
        print(io, sa)
    end

    # `Base.string(conf::ConfigurationR)`  ... provides a String notation for the variable conf::ConfigurationR.
    function Base.string(conf::ConfigurationR)
        wa = keys(conf.subshells);   va = values(conf.subshells)
        wb = Defaults.getDefaults("ordered subshell list: relativistic", 7)

        sa = ""
        for  k  in  wb
           if   k  in  wa
               occ = conf.subshells[k]
               sa = sa * string(k) * "^$occ "
           end
        end

        return( sa )
    end

    """
    `Base.:(==)(confa::ConfigurationR, confb::ConfigurationR)`
        ... compares (recursively) two relativistic configurations and return true if all subfields are equal,
            and false otherwise
    """
    function  Base.:(==)(confa::ConfigurationR, confb::ConfigurationR)
        if   confa.NoElectrons  != confb.NoElectrons    return( false )    end
        wk = keys(confa.subshells)
        for  k in wk
            if  !haskey(confb.subshells, k)  ||   confa.subshells[k] != confb.subshells[k]    return( false )    end
        end
        return( true )
    end

    """
    `struct  ManyElectron.LSjjSettings`  ... defines a type for the details and parameters for performing jj-LS expansions.

        + makeIt        ::Bool                         ... True, if the jj-LS expansion is to be made and false otherwise.
        + minWeight     ::Float64                      ... minimum weight with which a (relativistic) CSF much contribute to
                                                           (at least one) selected level.
        + printWeight   ::Float64                      ... minimum weight of a nonrelativistic CSF to be printed out in the
                                                           final expansion.
        + selectLines   ::Bool                         ... True, if lines are selected individually for the computations.
        + selectedLines ::Array{Tuple{Int64,Int64},1}  ... List of lines, given by tupels (inital-level, final-level).
    """
    struct LSjjSettings
        makeIt          ::Bool
        minWeight       ::Float64
        printWeight     ::Float64
        selectLines     ::Bool
        selectedLines   ::Array{Tuple{Int64,Int64},1}
    end

    """
    `LSjjSettings(makeIt::Bool)`  ... constructor for the default values of jj-LS transformations.
    """
    function LSjjSettings(makeIt::Bool)
        LSjjSettings(makeIt, 0.05, 0.1, false, Tuple{Int64,Int64}[])
    end

    # `Base.show(io::IO, settings::LSjjSettings)`  ... prepares a proper printout of the variable settings::LSjjSettings.
    function Base.show(io::IO, settings::LSjjSettings)
        println(io, "makeIt:            $(settings.makeIt)  ")
        println(io, "minWeight:         $(settings.minWeight)  ")
        println(io, "printWeight:       $(settings.printWeight)  ")
        println(io, "selectLines:       $(settings.selectLines)  ")
        println(io, "selectedLines:     $(settings.selectedLines)  ")
    end

    """
    `struct  ManyElectron.AsfSettings`
        ... a struct for defining the settings for the atomic state functions, i.e. the self-consistent-field (SCF)
            and CI computations

        + generateScf          ::Bool               ... True, if a SCF need to be generated, and false otherwise
                                                        (frozen orbitals).
        + breitScf             ::Bool               ... True, if Breit interaction is to be included into the SCF
                                                        computations.
        + methodScf            ::String             ... Specify the SCF method: ["AL", "OL", "EOL", "meanDFS", "meanHS"].
        + startScf             ::String             ... Specify how the start orbitals are obtained
                                                        ["fromNRorbitals", "fromGrasp", "hydrogenic"].
        + orbitalFileScf       ::String             ... Filename of orbitals, if taken from Grasp.
        + levelsScf            ::Array{Int64,1}     ... Levels on which the optimization need to be carried out.
        + maxIterationsScf     ::Int64              ... maximum number of SCF iterations
        + accuracyScf          ::Float64            ... convergence criterion for the SCF field.
        + shellSequenceScf     ::Array{Subshell,1}  ... Sequence of subshells to be optimized.

    	+ coulombCI            ::Bool               ... logical flag to include Coulomb interactions.
    	+ breitCI              ::Bool               ... logical flag to include Breit interactions.
    	+ qedModel             ::AbstractQedModel   ... specifies the applied model for estimating QED corrections;
    	                                                {NoneQed(), QedPetersburg(), QedSydney()}
    	+ methodCI             ::String             ... method for diagonalizing the matrix.
    	+ jjLS                 ::LSjjSettings       ... settings to enforce and control a jj-LS transformation of atomic level.
    	+ selectLevelsCI       ::Bool               ... true, if specific level (number)s have been selected.
    	+ selectedLevelsCI     ::Array{Int64,1}     ... Level number that have been selected.
    	+ selectSymmetriesCI   ::Bool               ... true, if specific level symmetries have been selected.
    	+ selectedSymmetriesCI ::Array{LevelSymmetry,1}    ... Level symmetries that have been selected.
    """
    struct  AsfSettings
        generateScf            ::Bool
        breitScf               ::Bool
        methodScf              ::String
        startScf               ::String
        orbitalFileScf         ::String
        levelsScf              ::Array{Int64,1}
        maxIterationsScf       ::Int64
        accuracyScf            ::Float64
        shellSequenceScf       ::Array{Subshell,1}
        #
    	coulombCI		       ::Bool
    	breitCI		           ::Bool
    	qedModel               ::AbstractQedModel
    	methodCI               ::String
    	jjLS                   ::LSjjSettings
    	selectLevelsCI         ::Bool
    	selectedLevelsCI       ::Array{Int64,1}
    	selectSymmetriesCI     ::Bool
    	selectedSymmetriesCI   ::Array{LevelSymmetry,1}
     end

    """
    `ManyElectron.AsfSettings()`  ... constructor for setting the default values.
    """
    function AsfSettings()
    	AsfSettings(false, false, "meanDFS", "hydrogenic", "", Int64[1], 24, 1.0e-6, Subshell[],
    	            true, false, NoneQed(), "eigen", LSjjSettings(false), false, Int64[], false, LevelSymmetry[] )
    end

    # `Base.show(io::IO, settings::AsfSettings)`  ... prepares a proper printout of the settings::AsfSettings.
    function Base.show(io::IO, settings::AsfSettings)
    	  println(io, "generateScf:          $(settings.generateScf)  ")
    	  println(io, "breitScf:             $(settings.breitScf)  ")
    	  println(io, "methodScf:            $(settings.methodScf)  ")
    	  println(io, "startScf:             $(settings.startScf)  ")
    	  println(io, "orbitalFileScf:       $(settings.orbitalFileScf)  ")
    	  println(io, "levelsScf:            $(settings.levelsScf)  ")
    	  println(io, "maxIterationsScf:     $(settings.maxIterationsScf)  ")
    	  println(io, "accuracyScf:          $(settings.accuracyScf)  ")
    	  println(io, "shellSequenceScf:     $(settings.shellSequenceScf)  ")
    	  #
    	  println(io, "coulombCI:            $(settings.coulombCI)  ")
    	  println(io, "breitCI:              $(settings.breitCI)  ")
    	  println(io, "qedModel :            $(settings.qedModel)  ")
    	  println(io, "methodCI:             $(settings.methodCI)  ")
    	  println(io, "jjLS :                $(settings.jjLS)  ")
    	  println(io, "selectLevelsCI:       $(settings.selectLevelsCI)  ")
    	  println(io, "selectedLevelsCI:     $(settings.selectedLevelsCI)  ")
    	  println(io, "selectSymmetriesCI:   $(settings.selectSymmetriesCI)  ")
    	  println(io, "selectedSymmetriesCI: $(settings.selectedSymmetriesCI)  ")

    end

    # `Base.string(settings::AsfSettings)`  ... provides a String notation for the variable settings::AsfSettings.
    function Base.string(settings::AsfSettings)
    	  error("Not yet implemented.")
    	  sa = "Asf settings: maximum No. of iterations = $(settings.maxIterationsScf), accuracy = (settings.accuracyScf)"
    	  return( sa )
    end

    """
    `struct  ManyElectron.CsfR`
        ... defines a type for a relativistic configuration state function (CSF) in terms of its orbitals (sequence of orbitals),
            their occupation as well as the angular momenta, seniority and the coupling of the corresponding antisymmetric subshell states.

        + useStandardSubshells ::Bool                 ... Determines whether the subshell list has to be obtained from some outer structure
                                                          given by the subshells field below.
        + J                    ::AngularJ64           ... Total angular momentum J.
        + parity               ::Parity               ... Total parity.
        + occupation           ::Array{Int64,1}       ... occupation of the orbitals with regard to the specified orbital list.
        + seniority            ::Array{Int64,1}       ... list of seniority values of the antisymmetric subshell states
        + subshellJ            ::Array{AngularJ64,1}  ... list of J-values of the antisymmetric subshell states.
        + subshellX            ::Array{AngularJ64,1}  ... intermediate X-values in the coupling of the antisymmetric subshell states.
        + subshells            ::Array{Subshell,1}    ... Explicitly given subshell list if useStandardSubshells == false.
    """
    struct  CsfR
        useStandardSubshells   ::Bool
        J		               ::AngularJ64
        parity  	           ::Parity
        occupation	           ::Array{Int64,1}
        seniority	           ::Array{Int64,1}
        subshellJ	           ::Array{AngularJ64,1}
        subshellX	           ::Array{AngularJ64,1}
        subshells	           ::Array{Subshell,1}
    end

    """
    `ManyElectron.CsfR(J::AngularJ64, parity::Parity)`
        ... simple constructor for given J and parity, and where standardOrbitals is set to false.
    """
    function CsfR(J::AngularJ64, parity::String)
    	  CsfR(false, J, parity, Int64[], Int64[], AngularJ64[], AngularJ64[], Subshell[] )
    end

    """
    `ManyElectron.CsfR(sa::String)`
        ... simple constructor for given closed-shell occupation (such as "[Ne]"), and where standardOrbitals is set to true.
    """
    function CsfR(sa::String)
    	  wa = subshellsFromClosedShellConfiguration(sa)

    	  occupation = Int64[];	  subshellJ = AngularJ64[];	subshellX  = AngularJ64[]
    	  subshells  = Subshell[];  seniority = Int64[]

        for  sb  in  wa
    	      #x wb = Basics.SubshellQuantumNumbers(string(sb));   kappa = wb[2];   occ = wb[4] + 1
    	      occ = Basics.subshell_2j(sb) + 1
    	      occupation = push!(occupation, occ)
    	      seniority  = push!(seniority,    0)
    	      subshellJ  = push!(subshellJ, AngularJ64(0) )
    	      subshellX  = push!(subshellX, AngularJ64(0) )
    	      subshells  = push!(subshells, sb )
    	  end
    	  J = AngularJ64(0);    parity = plus

    	  CsfR(false, J, parity, occupation, seniority, subshellJ, subshellX, subshells)
    end

    # `Base.show(io::IO, csf::CsfR)`  ... prepares a proper printout of the variable orbital::CsfR.
    function Base.show(io::IO, csf::CsfR)
    	  print(io, string(csf) )
    end

    # `Base.string(csf::CsfR)`  ... provides a String notation for csf::CsfR.
    function Base.string(csf::CsfR)
    	  sa = "\n   CSF: "
        if    csf.useStandardSubshells

    	      for  i in 1:length(csf.occupation)
    		    sa = sa * Basics.subshellStateString(string(Defaults.GBL_STANDARD_SUBSHELL_LIST[i]), csf.occupation[i], csf.seniority[i],
                                                  csf.subshellJ[i], csf.subshellX[i])
    		    sa = sa * ", "
    	      end
    	      sa = sa * ": J=" * string(csf.J) * string(csf.parity)

            else
    	      for  i in 1:length(csf.subshells)
    		    sa = sa * Basics.subshellStateString(string(csf.subshells[i]), csf.occupation[i], csf.seniority[i],
                                                  csf.subshellJ[i], csf.subshellX[i])
    		    sa = sa * ", "
    	      end
    	      sa = sa * ": J=" * string(csf.J) * string(csf.parity)
    	  end
    	  return( sa )
    end

    """
    `Base.:(==)(csfa::CsfR, csfb::CsfR)`
        ... compares (recursively) two relativistic CSFs and return true if all subfields are equal,
            and false otherwise
    """
    function  Base.:(==)(csfa::CsfR, csfb::CsfR)
        # Stop the comparison with an error message if the CSF are not defined w.r.t a common subshell list.
        if   !(csfa.useStandardSubshells)   ||   !(csfb.useStandardSubshells)	## ||   length(csfa.occupation) != length(csfb.occupation)
            error("To determine the equivalence of two CSF requires that both are defined on a common subshell list.")
        end

        if  length(csfa.occupation) != length(csfb.occupation)      return( false )    end
        if  csfa.J  !=  csfb.J  ||  csfa.parity  !=  csfb.parity	return( false )    end
    	if  csfa.occupation      !=  csfb.occupation			    return( false )    end
        if  csfa.seniority       !=  csfb.seniority			        return( false )    end
        if  csfa.subshellJ       !=  csfb.subshellJ			        return( false )    end
        if  csfa.subshellX       !=  csfb.subshellX			        return( false )    end

    	  return( true )
    end

    """
    `ManyElectron.CsfRGrasp92(subshells::Array{Subshell,1}, coreSubshells::Array{Subshell,1}, sa::String, sb::String, sc::String)`
        ... to construct a CsfR from 3 Grasp-Strings and the given list of subshells and core subshells.
    """
    function CsfRGrasp92(subshells::Array{Subshell,1}, coreSubshells::Array{Subshell,1}, sa::String, sb::String, sc::String)
    	occupation = Int64[];	 seniority = Int64[];	 subshellJ = AngularJ64[];    subshellX = AngularJ64[];    subshellsx = Subshell[]

    	# Define the occupation and quantum numbers of the closed shells
    	for i in 1:length(coreSubshells)
    	    push!(occupation, Basics.subshell_2j( coreSubshells[i] ) +1);
    	    push!(seniority, 0);     push!(subshellJ, AngularJ64(0));	 push!(subshellX, AngularJ64(0))
    	end

    	scc = sc * "	  ";   scx = " "
    	i  = -1;   ishell = length(coreSubshells)
    	while true
    	    i = i + 1;    sax = sa[9i+1:9i+9];     sbx = sb[9i+1:9i+9];    scx = strip( scc[9i+6:9i+14] )

    	    sh = strip( sax[1:5] );    sh = Basics.subshellGrasp(sh);    occ = parse( sax[7:8] )
    	    search(sbx, ',') > 0    &&    error("stop a: missing seniority; sb = $sb")   # Include seniority more properly if it occurs
    	    subJx = parse( sbx )
    	    if      typeof(subJx) == Void     subJ = AngularJ64(0)
    	    elseif  typeof(subJx) == Int64    subJ = AngularJ64(subJx)
    	    elseif  typeof(subJx) == Expr     subJ = AngularJ64( subJx.args[2]// subJx.args[3] )
    	    else    error("stop b; type = $(typeof(subJx))")
    	    end

    	    if   length(scx) > 0  &&  scx[end:end] in ["+", "-"]    scx = scx[1:end-1]   end
    	    subXx = parse( scx )
    	    if      ishell == 0  &&  typeof(subXx)     == Void                                       subX = subJ
    	    elseif  ishell >= 1  &&  subshellX[ishell] == AngularJ64(0)  && typeof(subXx) == Void    subX = subJ
    	    elseif  ishell >= 1  &&  typeof(subXx)     == Void                                       subX = subshellX[ishell]
    	    elseif  typeof(subXx) == Void                     subX = subshellX[ishell-1]
    	    elseif  typeof(subXx) == Int64                    subX = AngularJ64(subXx)
    	    elseif  typeof(subXx) == Expr                     subX = AngularJ64( subXx.args[2]// subXx.args[3] )
    	    else    error("stop c")
    	    end

    	    # Fill the arrays due to the pre-defined sequence of subshells
    	    while  true
    		ishell = ishell + 1
    		if  subshells[ishell] == sh
    		    wa = Basics.provide("subshell states: antisymmetric, seniority", sh, occ);    nu = -1
    		    for  a in wa
    			if   AngularJ64( a.Jsub2//2 ) == subJ	 nu = a.nu;    break   end
    		    end
    		    nu < 0    &&    error("stop d")
    		    push!(occupation, occ);    push!(seniority, nu);	 push!(subshellJ, subJ);    push!(subshellX, subX)
    		    break
    		else
    		    push!(occupation, 0);      push!(seniority, 0);	 push!(subshellJ, AngularJ64(0));    push!(subshellX, subshellX[ishell-1])
    		end
    	    end

    	    ishell > length(subshells)    &&	error("stop e")

    	    # Terminate if all subshells are read
    	    if  9i + 10 >  length(sa)	 break    end
    	end

      while  ishell < length(subshells)
    	    ishell = ishell + 1
          push!(occupation, 0);    push!(seniority, 0);	 push!(subshellJ, AngularJ64(0));    push!(subshellX, subshellX[ishell-1])
      end

    	# Fill the remaining subshells

    	J = subshellX[end];    scx = strip(sc)
    	if	  scx[end:end] == "+"  parity = Parity("+")
    	elseif  scx[end:end] == "-"  parity = Parity("-")
    	else	  error("stop f")
    	end

    	wa = CsfR(true, J, parity, occupation, seniority, subshellJ, subshellX, subshellsx)
    end

    """
    `ManyElectron.CsfRTransformToStandardSubshells(csf::CsfR, subshells::Array{Subshell,1})`
        ... constructor to re-define a CSF with regard to the given list of subshells. It incorporates 'empty' subshell state
            and their coupling (if necessary) but does not support a re-arrangement of the subshell order. This constructor
            is used to 'combine' individual CsfR into a common list (such as in a Basis) and, thus, use StandardSubshells = true
            is set. This method terminates with an error message if the given subshell list requires a re-coupling of the CSF.
    """
    function CsfRTransformToStandardSubshells(csf::CsfR, subshells::Array{Subshell,1})
    	J = csf.J;    parity = csf.parity;
    	occupation = Int64[];	 seniority = Int64[];	  subshellJ = AngularJ64[];	subshellX = AngularJ64[]

    	ish = 0
    	for  subsh in subshells
    	    push!(subshells, subsh)
    	    found = false
    	    for  i in 1:length(csf.subshells)
    		if     subsh == csf.subshells[i]  &&   i >   ish    i = ish;  found = true;   break
    		elseif subsh == csf.subshells[i]  &&   i <=  ish    error("CSF requires recoupling; this is not supported.")
    		end
    	    end
    	    if  found
    		push!(occupation, csf.occupation[ish]);    push!(seniority, csf.seniority[ish])
    		push!(subshellJ,  csf.subshellJ[ish]);     push!(subshellX, csf.subshellX[ish])
    	    else
    		push!(occupation, 0);	 push!(seniority, 0);	 push!(subshellJ,  0)
    		if    ish == 0    push!(subshellX, AngularJ64(0) )
    		else		  push!(subshellX, subshellX[end])
    		end
    	    end
    	end

    	CsfR(true, J, parity, occupation, seniority, subshellJ, subshellX, Subshell[] )
    end

    """
    `struct  ManyElectron.Basis`
        ... defines a type for a relativistic atomic basis, including the (full) specification of the configuration space
            and radial orbitals.

        + isDefined      ::Bool                      ... Determines whether this atomic basis 'points' to some a well-defined basis
                                                         or to an empty instance, otherwise.
        + NoElectrons    ::Int64                     ... No. of electrons.
        + subshells      ::Array{Subshell,1}         ... Explicitly given subshell list for this basis.
        + csfs           ::Array{CsfR,1}             ... List of CSF.
        + coreSubshells  ::Array{Subshell,1}         ... List of (one-electron) core orbitals.
        + orbitals       ::Dict{Subshell, Orbital}   ... Dictionary of (one-electron) orbitals.
    """
    struct  Basis
        isDefined	       ::Bool
    	  NoElectrons	 ::Int64
        subshells	       ::Array{Subshell,1}
        csfs             ::Array{CsfR,1}
    	  coreSubshells	 ::Array{Subshell,1}
    	  orbitals	       ::Dict{Subshell, Orbital}
    end

    """
    `ManyElectron.Basis()`  ... constructor for an 'empty' instance::Basis with isDefined == false.
    """
    function Basis()
    	  Basis(false, 0, Subshell[], CsfR[], Subshell[], Dict{Subshell, Orbital}() )
    end

    """
    `ManyElectron.Basis("from Grasp2013", cslFilename::String, rwfFilename::String)`
        ... to construct an instance::Basis from the Grasp92/Grasp2013 .csl and .rwf files.
    """
    function Basis(sa::String, cslFilename::String, rwfFilename::String)
        if  sa == "from Grasp2013"
           isDefined = true
           #
           println("ManyElectron.Basis-aa: warning ... The standard grid is set; make an interactive request if not appropriate.")
           grid     = Radial.Grid("grid: exponential")
           basis    = Basics.read("CSF list: Grasp92", cslFilename)
           orbitals = Basics.readOrbitalFileGrasp92(rwfFilename, grid)
        else  error("stop a")
        end

    	  Basis( isDefined, basis.NoElectrons, basis.subshells, basis.csfs, basis.coreSubshells, orbitals)
    end

    """
    `Base.show(io::IO, basis::Basis; full=false)`  ... prepares a proper printout of the variable basis::Basis.
    """
    function Base.show(io::IO, basis::Basis)
        print(io, string(basis), "\n")
    	  if  basis.isDefined
            println(io, "isDefined:          $(basis.isDefined)  ")
            println(io, "NoElectrons:        $(basis.NoElectrons)  ")
            println(io, "subshells:          $(basis.subshells)  ")
            println(io, "csfs:               $(basis.csfs)  ")
            println(io, "coreSubshells:      $(basis.coreSubshells)  ")
            println(io, "orbitals:           $(basis.orbitals)  ")
    	  end
    end

    # `Base.string(basis::Basis)`  ... provides a String notation for variable basis::Basis.
    function Base.string(basis::Basis)
    	if  basis.isDefined
    	    sa = "Atomic basis:  $(length(basis.csfs)) CSF defined for $(basis.NoElectrons) elecrons."
    	else
    	    sa = "Atomic basis:  Undefined."
    	end
    	return(sa)
    end

    """
    `struct  ManyElectron.Level`
        ... defines a type for an atomic level in terms of its quantum number, energy and with regard to a relativistic basis.

        + J            ::AngularJ64       ... Total angular momentum J.
        + M            ::AngularM64       ... Total projection M, only defined if a particular sublevel is referred to.
        + parity       ::Parity           ... Parity of the level.
        + index        ::Int64            ... Index of this level in its original multiplet, or 0.
        + energy       ::Float64          ... energy
        + relativeOcc  ::Float64          ... Relative occupation of this level, if involved in the evolution of a cascade.
        + hasStateRep  ::Bool             ... Determines whether this level 'point' to a physical representation of an state
                                              (i.e. an basis and corresponding mixing coefficient vector), or just to empty instances, otherwise.
        + basis        ::Basis            ... relativistic atomic basis
        + mc           ::Vector{Float64}  ... Vector of mixing coefficients w.r.t basis.
    """
    struct  Level
        J		     ::AngularJ64
        M		     ::AngularM64
    	  parity  	     ::Parity
    	  index  	     ::Int64
    	  energy  	     ::Float64
    	  relativeOcc    ::Float64
    	  hasStateRep    ::Bool
    	  basis	     ::Basis
    	  mc		     ::Vector{Float64}
    end

    """
    `ManyElectron.Level()`  ... constructor an empty level without useful data.
    """
    function Level()
        Level( AngularJ64(0),  AngularM64(0), Parity("+"), 0, 0., 0., false, Basis(), Vector{Float64}[] )
    end

    """
    `ManyElectron.Level(J::AngularJ64, parity::Parity, energy::Float64, relativeOcc::Float64)`
        ... constructor an atomic level without a state representation: hasStateRep == false.
    """
    function Level(J::AngularJ64, parity::Parity, energy::Float64, relativeOcc::Float64)
        Level(J, J, parity, 0, energy, relativeOcc, false, Basis(), Vector{Float64}[] )
    end

    # `Base.show(io::IO, level::Level)`  ... prepares a proper printout of the variable level::Level.
    function Base.show(io::IO, level::Level)
        println(io, "Level: J = $(level.J), M = $(level.M), parity = $(level.parity), index = $(level.index) ")
        println(io, "energy:         $(level.energy)  ")
        println(io, "relativeOcc:    $(level.relativeOcc)  ")
        println(io, "hasStateRep:    $(level.hasStateRep)  ")
        println(io, "basis:           (level.basis)  ")
        println(io, "mc:             $(level.mc)  ")
    end

    """
    `struct  ManyElectron.Multiplet`
        ... defines a type for an ordered list of atomic levels; has only the default constructor.

    	  + name  	 ::String	         ... A name associated to the multiplet.
    	  + levels	 ::Array{Level,1}   ... A list of levels (pointers).
    """
    struct  Multiplet
        name       ::String
    	  levels  	 ::Array{Level,1}
    end

    """
     `ManyElectron.Multiplet()`  ... constructor for providing an 'empty' instance of this struct.
    """
    function Multiplet()
        Multiplet("", Level[] )
    end

    """
    `ManyElectron.Multiplet("from Ratip2012", cslFilename::String, rwfFilename::String, mixFilename::String)`
        ... to construct an instance::Multiplet from the Grasp92/Ratip2012 .csl, .rwf and .mix files.  Some consistency
            checks are made and the method terminates with an error message if the files don't fit together.

    + `("from Grasp2018", cslFilename::String, rwfFilename::String, mixFilename::String)`
        ... to do the same but for files from Grasp2018.
    """
    function Multiplet(sa::String, cslFilename::String, rwfFilename::String, mixFilename::String)
        if  sa == "from Ratip2012"
           name      = "from "*cslFilename[1:end-4]
           basis     = ManyElectron.Basis("from Grasp2013", cslFilename, rwfFilename)
           multiplet = Basics.readMixFileRelci(mixFilename, basis)
        else  error("Unrecognized keystrings:  $sa")
        end

        return( multiplet )
    end

    # `Base.show(io::IO, multiplet::Multiplet)`  ... prepares a proper printout of the variable multiplet::Multiplet.
    function Base.show(io::IO, multiplet::Multiplet)
        println(io, "name:        $(multiplet.name)  ")
        println(io, "levels:      $(multiplet.levels)  ")
    end

    """
    `struct  ManyElectron.MultipletSettings`
        ... a struct for defining the atomic interactions to be incorporated into the representation of a multiplet.

    	+ Coulomb		           ::Bool 		      ... logical flag to include Coulomb interactions.
    	+ Breit 		           ::Bool 		      ... logical flag to include Breit interactions.
    	+ QED			           ::Bool 		      ... logical flag to include QED interactions.
    	+ diagonalizationMethod    ::String		      ... method for diagonalizing the matrix.
    	+ selectLevels  	       ::Bool 		      ... true, if specific level (number)s have been selected for computation.
    	+ selectedLevels	       ::Array{Int64,1}	  ... Level number that have been selected.
    	+ selectSymmetries	       ::Bool 		      ... true, if specific level symmetries have been selected for computation.
    	+ selectededSymmetries     ::Array{LevelSymmetry,1}   ... Level symmetries that have been selected.

    	### `Methods for diagonalization`

    	* `eigval`	       The internal Julia method for diagonalizing a quadratic matrix
    	*
    """
    struct  MultipletSettings
        Coulomb 		         ::Bool
        Breit			         ::Bool
        QED			             ::Bool
    	diagonalizationMethod	 ::String
        selectLevels		     ::Bool
        selectedLevels  	     ::Array{Int64,1}
        selectSymmetries	     ::Bool
        selectedSymmetries	     ::Array{LevelSymmetry,1}
    end

end # module


text:
function searchsortedrange(a, x)
    x0, x1 = extrema(x)
    i0 = searchsortedfirst(a, x0)
    i1 = searchsortedlast(a, x1)
    return i0:i1
end

# Only apply `f` in place within the shapes provided when they are all inbounds
function broadcast_within!(f, args...)
    idxs_within = map(args) do (array, idxs)
        return map(searchsortedrange, idxs, axes(array))
    end
    selected_idxs = map(intersect, idxs_within...)
    views = map(args) do (array, idxs)
        return view(array, map(getindex, idxs, selected_idxs)...)
    end
    broadcast!(f, views...)
    (target, _), _... = args
    return target
end

add_within!(t, args...) = broadcast_within!(+, t, t, args...)

"""
    slice!(buffer::AbstractVector, parent::AbstractArray, shape::Tuple{Vararg{AbstractRange}})

Non-allocating contiguous view of `parent` data at indices `shape`, which uses
memory reserved in `buffer`. Indices outside the axes of `parent` are set to `0`.
"""
function slice!(buffer::AbstractVector, parent::AbstractArray, shape::Tuple{Vararg{AbstractRange}})
    # maybe consider `slice = unsafe_wrap(Array, pointer(buffer), map(length, shape))` for performance
    trimmed = view(buffer, Base.OneTo(prod(length, shape)))
    slice = reshape(trimmed, map(length, shape))
    fill!(slice, 0)
    return broadcast_within!(identity, slice => axes(slice), parent => shape)
end

function slice(parent, shape)
    slice = similar(parent, map(length, shape))
    fill!(slice, 0)
    return broadcast_within!(identity, slice => axes(slice), parent => shape)
end

struct Slices{T, U, V}
    parent::T
    shapes::U
    buffer::V
end

function Slices(parent, shapes)
    l = maximum(shape -> prod(length, shape), shapes)
    buffer = similar(parent, l)
    return Slices(parent, shapes, buffer)
end

Base.getindex(s::Slices, i) = slice!(s.buffer, s.parent, s.shapes[i])

function selectchannels(y::AbstractArray{T, N}, rg::AbstractRange) where {T, N}
    shape = ntuple(n -> n == N - 1 ? rg : Colon(), N)
    return view(y, shape...)
end

function expandchannels(y₀::AbstractArray{T, N}, rg::AbstractRange) where {T, N}
    if axes(y₀, N - 1) ⊈ rg
        msg = "input has incompatible channels: $(axes(y₀, N - 1)) ⊈ $(rg)"
        throw(ArgumentError(msg))
    end
    shape = ntuple(n -> n == N - 1 ? rg : axes(y₀, n), N)
    return slice(y₀, shape)
end

ext: jl


text:
"""

Create a `name`_summary.csv file.

$(SIGNATURES)

# Extended help

### Required arguments
```julia
* `model::SampleModel             : SampleModel
```

### Optional arguments
```julia
* `printsummary=false             : Display summary
```

After completion a ..._summary.csv file has been created.
This file can be read as a DataFrame in by `df = read(summary(model))`

"""
function stan_summary(
  model::T,
  printsummary=false) where {T <: CmdStanModels}

  #local csvfile
  n_chains = get_n_chains(model)

  samplefiles = String[]
  for i in 1:get_n_chains(model)
    push!(samplefiles, "$(model.output_base)_chain_$(i).csv")
  end
  try
    pstring = joinpath("$(model.cmdstan_home)", "bin", "stansummary")
    if Sys.iswindows()
      pstring = pstring * ".exe"
    end
    csvfile = "$(model.output_base)_summary.csv"
    isfile(csvfile) && rm(csvfile)
    cmd = `$(pstring) -c $(csvfile) $(par(samplefiles))`
    outb = IOBuffer()
    run(pipeline(cmd, stdout=outb));
    if printsummary
      cmd = `$(pstring) $(par(samplefiles))`
      resfile = open(cmd; read=true);
      print(read(resfile, String))
    end
  catch e
    println(e)
  end

  return

end


text:
"""
    diagview(M) = view(M, diagind(M))

Like `diag(M)` but makes a view.
"""
diagview(A::AbstractMatrix) = view(A, diagind(A))

diagview(A::LinearAlgebra.Diagonal) = A.diag

"""
    rvec(x) = [x]
    rvec(A) = vec(A)

Really-vec... extends `LinearAlgebra.vec` to work on scalars too.
"""
rvec(x::Number) = [x]
rvec(A) = LinearAlgebra.vec(A)

"""
    mul!(Z,A,B)

Exactly `LinearAlgebra.mul!` except that it can write into a zero-array.
"""
mul!(Z,A,B) = LinearAlgebra.mul!(Z,A,B)
mul!(Z::AbstractArray{T,0}, A,B) where {T} = copyto!(Z, A * B)

"""
    star(x,y,...)

Used for multiplying axes now, not sizes.
"""
star(x, y) = Base.OneTo(length(x) * length(y))
star(x,y,zs...) = star(star(x,y), zs...)


text:
import LibGit2
using CBindingGen
using Libdl

# Set base path
deps_dir = @__DIR__

# Clone/update repo
krome_repo = "https://bitbucket.org/tgrassi/krome.git"
clone_path = joinpath(deps_dir, "krome")
build_dir = joinpath(clone_path, "build")
build_dir_backup = joinpath(clone_path, "build.orig")
if isdir(clone_path)
  @info "Restoring build directory..."
  rm(build_dir, recursive=true)
  mv(build_dir_backup, build_dir)

  @info "Updating repository in $clone_path..."
  repo = LibGit2.GitRepo(clone_path)
  LibGit2.fetch(repo)
  LibGit2.merge!(repo, fastforward=true)

  @info "Creating a backup of the build directory..."
  cp(build_dir, build_dir_backup)
else
  @info "Cloning repository $krome_repo to $clone_path..."
  LibGit2.clone(krome_repo, clone_path)

  @info "Creating a backup of the build directory..."
  cp(build_dir, build_dir_backup)
end

# Call krome to set up a build
reactions_target = joinpath(deps_dir, "reactions_verbatim.dat")
krome_default_args = ["-interfacePy", "-interfaceC", "-unsafe",
                      "-verbatimFilename=$reactions_target"]
krome_custom_args = split(get(ENV, "JULIA_KROME_CUSTOM_ARGS", "-test=hello"), ";")
krome_args = vcat(krome_default_args, krome_custom_args)
python3_exec = get(ENV, "JULIA_KROME_PYTHON3_EXEC", "python3")

cd(clone_path) do
  setup_cmd = `$python3_exec krome $krome_args`
  # Set locale if unset
  if !haskey(ENV, "LC_ALL")
    @info "Setting environment variable LC_ALL to 'C.UTF-8' to avoid Python errors"
    ENV["LC_ALL"] = "C.UTF-8"
  end
  @info "Running command '$setup_cmd' in '$clone_path'..."
  run(setup_cmd)
end

# Call make
if lowercase(get(ENV, "JULIA_KROME_OPTIMIZED_BUILD", "true")) == "true"
  make_target = "csharedlib_opt"
else
  make_target = "csharedlib"
end

cd(build_dir) do
  make_cmd = `make $make_target`
  @info "Running command '$make_cmd'..."
  run(make_cmd)
end

# Generate C bindings
@info "Generate C bindings..."

# Manually set header files to consider
hdrs = ["krome.h", "krome_user.h"]

# Build list of arguments for Clang
include_args = String[]
include_directories = [normpath(build_dir)]
for dir in include_directories
  append!(include_args, ("-I", dir))
end

# Convert symbols in header
cvts = convert_headers(hdrs, args=include_args) do cursor
  header = CodeLocation(cursor).file
  name   = string(cursor)

  # only wrap the KROME headers
  dirname, filename = splitdir(header)
  if !(filename in hdrs)
    return false
  end

  return true
end

# Write generated C bindings to file
@info "Write generated C bindings to file..."
krome_library = joinpath(deps_dir, "libkrome." * Libdl.dlext)
const bindings_filename = joinpath(deps_dir, "libkrome.jl")
open(bindings_filename, "w+") do io
  generate(io, krome_library => cvts)
end

# Final procedure: copy products (library + reactions) to `deps/` such that in case of a build
# problem, the previous files remain intact
# Copy KROME library from build directory to deps/
library_source = joinpath(build_dir, "libkrome." * Libdl.dlext)
library_target = krome_library
@info "Copying '$library_source' to '$library_target'..."
cp(library_source, library_target, force=true)

# Copy reactions_verbatim.dat from build directory to deps/
reactions_source = joinpath(build_dir, "reactions_verbatim.dat")
@info "Copying '$reactions_source' to '$reactions_target'..."
cp(reactions_source, reactions_target, force=true)


text:
#######################################
#
# sfp file
#
#######################################

"""
Read sfp file containing sensor locations

#### Input
* `fname`: Path for the sfp file

#### Output
* `elec`: Electrodes object
"""
function read_sfp(fname::AbstractString; coordinate = Talairach)

    @info("Reading dat file = $fname")

    # Create an empty electrode set
    elecs = Electrode[]

    # Read file and match to expected file format
    file = read(fname, String)
    regexp = r"(\S+)\s+(\S+)\s+(\S+)\s+(\S+)"
    m = collect((m.match for m in eachmatch(regexp, file)))

    # Convert label to ascii and remove '
    for idx = 1:length(m)
        local_matches = match(regexp, m[idx])
        push!(
            elecs,
            Electrode(
                replace(local_matches[1], "'" => ""),
                coordinate(
                    parse(Float64, local_matches[2]),
                    parse(Float64, local_matches[3]),
                    parse(Float64, local_matches[4]),
                ),
                Dict(),
            ),
        )
    end

    @debug("Imported $(length(elecs)) electrodes")

    return elecs
end


text:
### Julia OpenStreetMap Package ###
### MIT License                 ###
### Copyright 2014              ###

### Default Speed Limits in Kilometers Per Hour ###
const SPEED_ROADS_URBAN = Compat.@Dict(
    1 => 95,    # Motorway
    2 => 72,    # Trunk
    3 => 48,    # Primary
    4 => 32,    # Secondary
    5 => 22,    # Tertiary
    6 => 12,    # Residential/Unclassified
    7 => 8,     # Service
    8 => 5)     # Living street

const SPEED_ROADS_RURAL = Compat.@Dict(
    1 => 110,
    2 => 90,
    3 => 80,
    4 => 72,
    5 => 55,
    6 => 40,
    7 => 15,
    8 => 15)


text:
@testset "Test bisection method for strictly increasing function" begin
    f(x) = tanh(3.0 + 2.0* x)

    a = -3.0
    b = 3.0
    ϵx = 1e-14
    ϵf = 1e-14
    x = a + 0.5*(b-a)

    ya = f(a)
    yb = f(b)

    while b - a > ϵx || abs(f(x)) > ϵf
        x, fx, a, ya, b, yb = bisection(x, f(x), a , ya, b, yb)
    end

    @test (b - a) <= ϵx
    @test abs(f(x)) <= ϵf
end

@testset "Test bisection, Newton, hybrid method for wavy globally increasing function" begin
    f(x) = 3*x - exp(0.5-x) + cos(3.5π*x)
    g(x) = ForwardDiff.derivative(f, x)

    xroot = 0.4144911998972488

    a = -5.0
    b = 7.0
    ϵx = 1e-6
    ϵf = 1e-6
    # Start from midpoint
    x = -4.9#a + 0.5*(b-a)
    xold = x
    fx = f(x)
    gx = g(x)

    ya = f(a)
    yb = f(b)

    # newton method
    while abs(fx) > ϵf
        fx = f(x)
        gx = g(x)
        x -= fx/gx
    end

    xrootnewton = x
    @test isapprox(xroot, xrootnewton, atol = ϵx)
    @test isapprox(f(xrootnewton), 0.0, atol = ϵx)

    # bisection method

    a = -5.0
    b = 7.0
    ϵx = 1e-4
    ϵf = 1e-4
    # Start from midpoint
    x = -4.9#a + 0.5*(b-a)
    xold = x
    fx = f(x)
    gx = g(x)

    ya = f(a)
    yb = f(b)

    # bisection method
    while b - a > ϵx
        x, fx, a, ya, b, yb = bisection(x, f(x), a , ya, b, yb)
    end

    xrootbisection = x
    @test isapprox(xroot, xrootbisection, atol = ϵx)
    @test isapprox(a, b, atol = ϵx)

    # hybrid method Newton + bisection

    a = -5.0
    b = 7.0
    ϵx = 1e-4
    ϵf = 1e-4

    # Start with a bad initial guees
    x = -4.9

    xroothybrid = hybridsolver(f, g, x, a, b)

    @test isapprox(xroot, xroothybrid, atol = ϵx)
end

@testset "Test hybridinverse function" begin

    Ne = 200
    Nx = 3

    idx = [0 0 0; 0 0 1; 0 0 2; 0 1 0; 0 0 1; 0 1 1; 0 2 1; 0 3 1; 1 3 1]

    Nψ = size(idx,1)

    coeff = randn(Nψ)

    B = MultiBasis(CstProHermiteBasis(6), Nx)
    f = ExpandedFunction(B, idx, coeff)
    R = IntegratedFunction(f)

    X = randn(Nx, Ne) .* randn(Nx, Ne)

    F = evaluate(R, X)

    Xmodified = deepcopy(X)

    Xmodified[end,:] .+= cos.(2*π*randn(Ne)) .* exp.(-randn(Ne).^2/2)

    @test norm(Xmodified[end,:]-X[end,:])>1e-6

    S = Storage(f, Xmodified)

    hybridinverse!(Xmodified, F, R, S)

    @test norm(Xmodified[end,:] - X[end,:])<2e-8

    @test norm(Xmodified[1:end-1,:] - X[1:end-1,:])<2e-8
end

@testset "Test hybridinverse HermiteMapComponent" begin

    Ne = 200
    Nx = 3

    idx = [0 0 0; 0 0 1; 0 0 2; 0 1 0; 0 0 1; 0 1 1; 0 2 1; 0 3 1; 1 3 1]

    Nψ = size(idx,1)

    coeff = randn(Nψ)

    C = HermiteMapComponent(20, Nx, idx, coeff)

    X = randn(Nx, Ne) .* randn(Nx, Ne)

    F = evaluate(C.I, X)

    Xmodified = deepcopy(X)

    Xmodified[end,:] .+= cos.(2*π*randn(Ne)) .* exp.(-randn(Ne).^2/2)

    @test norm(Xmodified[end,:]-X[end,:])>1e-6

    S = Storage(C.I.f, Xmodified)

    hybridinverse!(Xmodified, F, C, S)

    @test norm(Xmodified[end,:] - X[end,:])<1e-5

     # Check that the first entries are not modified
    @test norm(Xmodified[1:end-1,:] - X[1:end-1,:])<1e-5

    # Verify the inversion
    @test norm(evaluate(C.I, X) - evaluate(C.I, Xmodified))<1e-5
end


text:
module CYCLOPS_v6_2a_CircularStats_U

export Circular_Mean
export Fischer_Circular_CorrelationMeasures
export Jammalamadka_Circular_CorrelationMeasures
export Circular_Error
export Circular_Error_List

using MultivariateStats
using StatsBase

############################################################
function Circular_Mean(phases::Array{Float64,1})
  sinterm=sum(sin.(phases))
  costerm=sum(cos.(phases))
  atan2(sinterm,costerm)
end

############################################################
function Fischer_Circular_Correlations(rphases,sphases)
	n1=length(rphases)
	n2=length(sphases)

	num=n1

	rphases=mod.(rphases+2*pi,2*pi)
	sphases=mod.(sphases+2*pi,2*pi)

	numtot=0.
	d1tot=0.
	d2tot=0.
	for i in 1:num
		for j in (i+1):num
			numeratorterm=sin(sphases[i]-sphases[j])*sin(rphases[i]-rphases[j])
			denomterm1=(sin(sphases[i]-sphases[j]))^2
			denomterm2=(sin(rphases[i]-rphases[j]))^2
			numtot=numtot+numeratorterm
			d1tot=d1tot+denomterm1
			d2tot=d2tot+denomterm2
		end
	end
	fischercor=numtot/(sqrt(d1tot)*sqrt(d2tot))
	fischercor
end
############################################################
function Jammalamadka_Circular_Correlations(rphases,sphases)
	numtot=0.
	d1tot=0.
	d2tot=0.

	rbar=mod.(2*pi+Circular_Mean(rphases),2*pi)
	sbar=mod.(2*pi+Circular_Mean(sphases),2*pi)

	numtot=sum(sin.(rphases-rbar) .* sin.(sphases-sbar))
	d1tot=sqrt(sum( sin.(rphases-rbar) .^ 2))
	d2tot=sqrt(sum( sin.(sphases-sbar) .^ 2))
	Jammalamadka=numtot/(d1tot*d2tot)
	Jammalamadka
end

############################################################
function Circular_Rank_Phases(rphases)
	number=length(rphases)
	rphases=mod.(rphases+2*pi,2*pi)
	rranks=tiedrank(rphases)
	rrankphases=rranks*2*pi/number
	rrankphases
end

############################################################
function Jammalamadka_Rank_Circular_Correlations(rphases,sphases)

	rphases=Circular_Rank_Phases(rphases)
	sphases=Circular_Rank_Phases(sphases)

	r_minus_s_bar=mod(atan2(sum(sin.(rphases-sphases)),sum(cos.(rphases-sphases))),2*pi)
	r_plus_s_bar=mod(atan2(sum(sin.(rphases+sphases)),sum(cos.(rphases+sphases))),2*pi)

	Ntot=length(rphases)

	term1=cos.(rphases-sphases-r_minus_s_bar)
	term2=cos.(rphases+sphases-r_plus_s_bar)

	Jammalamadka=1/Ntot*(sum(term1))-1/Ntot*(sum(term2))
	Jammalamadka
end
############################################################
##########################################################

function FindComponentAngles(angle_sum,angle_diff)
    rang=(angle_sum + angle_diff)/2
    sang=(angle_sum - angle_diff)/2
	[rang,sang]
end
##########################################################
function Jammalamadka_Uniform_Circular_Correlations(rphases,sphases)
# This is a modification of he Jammalamadka Circular Correlation
# Described in Topics in Circular Statistics
# It is required beacuse the circular average is not well defined with circular uniform data
# This measure should only be used when 1 or both of the data sets being compared are uniform

	rphases=mod.(rphases,2*pi)
	sphases=mod.(sphases,2*pi)

	r_minus_s_bar=mod(atan2(sum(sin.(rphases-sphases)),sum(cos.(rphases-sphases))),2*pi)
	r_plus_s_bar=mod(atan2(sum(sin.(rphases+sphases)),sum(cos.(rphases+sphases))),2*pi)

	bars=FindComponentAngles(r_plus_s_bar,r_minus_s_bar)
	rbar=bars[1]
	sbar=bars[2]

	numtot=sum(sin.(rphases-rbar) .* sin.(sphases-sbar))
	d1tot=sqrt(sum( sin.(rphases-rbar) .^ 2))
	d2tot=sqrt(sum( sin.(sphases-sbar) .^ 2))

	Jammalamadka=numtot/(d1tot*d2tot)
	Jammalamadka
end

############################################################

function Fischer_Circular_CorrelationMeasures(rphases,sphases)

	rrankphases=Circular_Rank_Phases(rphases)
	srankphases=Circular_Rank_Phases(sphases)

	F =Fischer_Circular_Correlations(rphases,sphases)

	FR =Fischer_Circular_Correlations(rrankphases,srankphases)

	[F,FR]
end

############################################################

function Jammalamadka_Circular_CorrelationMeasures(rphases,sphases)

	J =Jammalamadka_Circular_Correlations(rphases,sphases)

	JU =Jammalamadka_Uniform_Circular_Correlations(rphases,sphases)

	JR =Jammalamadka_Rank_Circular_Correlations(rphases,sphases)

	[J, JU, JR]
end
############################################################

function Circular_Error(truth,estimate)
	truth=mod(2*pi+truth,2*pi)
	estimate=mod(2*pi+estimate,2*pi)

	diff1=mod((truth-estimate),2*pi)
	diff2=diff1-(2*pi)

	if (abs(diff1)<=abs(diff2))
		diff=diff1
	end

	if (abs(diff2)<abs(diff1))
		diff=diff2
	end

	e=diff
	e
end
############################################################

function Circular_Error_List(true_list,estimate_list)
	n1=length(true_list)
	n2=length(estimate_list)
	if (n1 != n2) print("Error Warning") end

	error_list=zeros(n1)
	for count in 1:n1
		error_list[count]=Circular_Error(true_list[count],estimate_list[count])
	end
	error_list
end

############################################################
end

max_stars_repo_name: AdithyaVenkateshMohan/CYCLOPS

max_issues_repo_name: AdithyaVenkateshMohan/CYCLOPS

max_forks_repo_name: AdithyaVenkateshMohan/CYCLOPS


text:
import Base: size, show, summary, similar

abstract type PointData{N,T} <: AbstractVector{T} end

"""
    ScalarData <: PointData

A wrapper for a one-dimensional array of scalar-valued data. The resulting wrapper
can be indexed in the same way as the array itself.

# Constructors
- `ScalarData(d[,dtype=Float64])` constructs a wrapper for the one-dimensional array of data `d`
- `ScalarData(n::Int)` constructs a wrapper for an array of zeros of length `n`.
- `ScalarData(x::PointData)` constructs a wrapper for an array of zeros of the
   same length as that wrapped by `x`.
- `ScalarData(n::Int,dtype=ComplexF64)` constructs a wrapper for complex-valued data.
- `ScalarData(x::PointData,dtype=ComplexF64)` constructs a wrapper for an array of
   complex zeros of the same length as that wrapped by `x`.

# Example

```jldoctest
julia> f = ScalarData(10);

julia> f[5] = 1.0;

julia> f
10 points of scalar-valued Float64 data
10-element Array{Float64,1}:
 0.0
 0.0
 0.0
 0.0
 1.0
 0.0
 0.0
 0.0
 0.0
 0.0
```
"""
struct ScalarData{N,T} <: PointData{N,T}
    data::Vector{T}
end

@wraparray ScalarData data 1

"""
    VectorData <: GridData

A wrapper for a one-dimensional array of two-component vector-valued data. The
resulting wrapper can be indexed as though the first component and second
component are stacked on top of each other.

# Constructors
- `VectorData(u,v)` constructs a wrapper for the vector components data `u` and `v`.
- `VectorData(n::Int)` constructs a wrapper with zeros of length `n` for both components.
- `VectorData(x::GridData)` constructs a wrapper for zero components of the
   same length as that wrapped by `x`.
- `VectorData(n::Int,dtype=ComplexF64)` constructs a wrapper with complex-valued zeros
   of length `n` for both components.

# Example

```jldoctest
julia> f = VectorData(10,dtype=ComplexF64);

julia> f.v[1:5] = 1:5;

julia> f
10 points of vector-valued Complex{Float64} data
10×2 Array{Complex{Float64},2}:
 0.0+0.0im  1.0+0.0im
 0.0+0.0im  2.0+0.0im
 0.0+0.0im  3.0+0.0im
 0.0+0.0im  4.0+0.0im
 0.0+0.0im  5.0+0.0im
 0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im

julia> f[7] = 1im; f[18] = 0.2;

julia> f
10 points of vector-valued Complex{Float64} data
10×2 Array{Complex{Float64},2}:
 0.0+0.0im  1.0+0.0im
 0.0+0.0im  2.0+0.0im
 0.0+0.0im  3.0+0.0im
 0.0+0.0im  4.0+0.0im
 0.0+0.0im  5.0+0.0im
 0.0+0.0im  0.0+0.0im
 0.0+1.0im  0.0+0.0im
 0.0+0.0im  0.2+0.0im
 0.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im
```
"""
struct VectorData{N,T} <: PointData{N,T}
    u::ScalarData{N,T}
    v::ScalarData{N,T}
end

"""
    TensorData

A wrapper for a one-dimensional array of 2x2 tensor-valued data, with fields
`dudx`, `dudy`, `dvdx`, `dvdy`. The resulting wrapper can be indexed as though these four components are stacked
on top of each other.

# Constructors
- `TensorData(dudx,dudy,dvdx,dvdy)` constructs a wrapper for the tensor components data.
- `TensorData(n::Int)` constructs a wrapper with zeros of length `n` for all components.
- `TensorData(x::ScalarData/VectorData/TensorData)` constructs a wrapper for zero components of the
   same length as that wrapped by `x`.

# Example

"""
struct TensorData{N,T} <: PointData{N,T}
    dudx::ScalarData{N,T}
    dudy::ScalarData{N,T}
    dvdx::ScalarData{N,T}
    dvdy::ScalarData{N,T}
end

function ScalarData(data::Vector{T}) where {T <: Number}
  ScalarData{length(data),T}(data)
end

function VectorData(u::Vector{T},v::Vector{T}) where {T <: Number}
  @assert length(u) == length(v)
  VectorData{length(u),T}(ScalarData(u),ScalarData(v))
end

function TensorData(dudx::Vector{T},dudy::Vector{T},
                    dvdx::Vector{T},dvdy::Vector{T}) where {T <: Number}
  @assert length(dudx) == length(dudy) == length(dvdx) == length(dvdy)
  TensorData{length(dudx),T}(ScalarData(dudx),ScalarData(dudy),
                           ScalarData(dvdx),ScalarData(dvdy))
end

ScalarData(n::Int;dtype=Float64) = ScalarData(zeros(dtype,n))

VectorData(x::Tuple{Vector{T},Vector{T}}) where {T <: Number} = VectorData(x[1],x[2])
VectorData(n::Int;dtype=Float64) = VectorData(zeros(dtype,n),zeros(dtype,n))

TensorData(x::NTuple{4,Vector{T}}) where {T <: Number} = TensorData(x[1],x[2],x[3],x[4])
TensorData(n::Int;dtype=Float64) = TensorData(zeros(dtype,n),zeros(dtype,n),zeros(dtype,n),zeros(dtype,n))

for f in (:ScalarData,:VectorData,:TensorData)
  @eval (::Type{$f{N,T}})() where {N,T} = $f(N,dtype=T)
  @eval $f(::PointData{N,T};dtype=T) where {N,T} = $f(N,dtype=dtype)
  @eval similar(::$f{N,T};element_type=T) where {N,T} = $f(N,dtype=element_type)
end

"""
    size(A::VectorData,d::Int) -> Int

Return twice the number of vector data points if `d` is 1 (the sum of the length of the `u` and
`v` vectors) and 1 if `d` is 2. This is consistent with the interpretation of VectorData as
a stacked pair of columns, corresponding to the `u` and `v` components, respectively.
"""
Base.size(::VectorData{N,T},d::Int) where {N,T} = d == 1 ? 2*N : 1

"""
    size(A::VectorData) -> Tuple

Return a tuple of the number of vector data points by the number of dimensions.
"""
Base.size(A::VectorData) = (size(A,1),)
@propagate_inbounds Base.getindex(A::VectorData{N,T},i::Int) where {N,T} =
   i > N ? A.v[i-N] : A.u[i]
@propagate_inbounds Base.setindex!(A::VectorData{N,T}, v, i::Int) where {N,T} =
   i > N ? A.v[i-N] = convert(T, v) : A.u[i] = convert(T, v)

"""
    size(A::TensorData,d::Int) -> Int

Return four times the number of tensor data points if `d` is 1 (the sum of the length of the four components)
and 1 if `d` is 2. This is consistent with the interpretation of TensorData as
a stacked set of columns.
"""
Base.size(::TensorData{N,T},d::Int) where {N,T} = d == 1 ? 4*N : 1

"""
    size(A::TensorData) -> Tuple

Return a tuple of the number of tensor data points by the number of dimensions.
"""
Base.size(A::TensorData) = (size(A,1),)

@propagate_inbounds Base.getindex(A::TensorData{N,T},i::Int) where {N,T} =
   i > N ? (i > 2*N ? (i > 3*N ? A.dvdy[i-3*N] : A.dvdx[i-2*N]) : A.dudy[i-N] ) : A.dudx[i]
@propagate_inbounds Base.setindex!(A::TensorData{N,T}, v, i::Int) where {N,T} =
   i > N ? (i > 2*N ? (i > 3*N ? A.dvdy[i-3*N] = convert(T, v) : A.dvdx[i-2*N] = convert(T, v) ) : A.dudy[i-N] = convert(T, v) ) : A.dudx[i] = convert(T, v)

function show(io::IO, m::MIME"text/plain", pts::ScalarData{N,T}) where {N,T}
  println(io,"$N points of scalar-valued $T data")
  show(io,m,pts.data)
end

function show(io::IO, m::MIME"text/plain", pts::VectorData{N,T}) where {N,T}
  println(io,"$N points of vector-valued $T data")
  show(io,m,hcat(pts.u,pts.v))
end

function show(io::IO, m::MIME"text/plain", pts::TensorData{N,T}) where {N,T}
  println(io,"$N points of tensor-valued $T data dudx, dudy, dvdx, dvdy")
  show(io,m,hcat(pts.dudx,pts.dudy,pts.dvdx,pts.dvdy))
end

include("basicpointoperations.jl")


text:
import HTTP

"""Pluto Events interface

Use this interface to hook up functionality into the Pluto world.
Events are guaranteed to be run at least every time something interesting happens,
but keep in mind that Pluto may be run the events multiple times "logically". For instance,
the FileSaveEvent may be triggered whenever pluto wants to make sure the file is saved,
which may be more often than the file is actually changed. Deduplicate on your own if
you care about this.

Define your function to handle the events using multiple dispatch:

First assign a handler for all the types you will not use, using the supertype:

```julia-repl
julia> function myfn(a::PlutoEvent)
            nothing
        end
```

And then create a special function for each event you want to handle specially

```julia-repl
julia> function myfn(a::FileSaveEvent)
            HTTP.post("https://my.service.com/count_saves")
        end
```

Finally, pass the listener to Pluto's configurations with a keyword argument

```julia-repl
julia> Pluto.run(; on_event = myfn)
```
"""
abstract type PlutoEvent end

function try_event_call(session, event::PlutoEvent)
    return try
        session.options.server.on_event(event)
    catch e
        @warn "Couldn't run event listener" event exception=(e, catch_backtrace())
        nothing
    end
end

# Triggered when the web server gets started
struct ServerStartEvent <: PlutoEvent
    address::String
    port::UInt16
end

# Triggered when a notebook is saved
struct FileSaveEvent <: PlutoEvent
    notebook::Notebook
    file_contents::String
    path::String
end

FileSaveEvent(notebook::Notebook) = begin
    file_contents = sprint() do io
        save_notebook(io, notebook)
    end
    FileSaveEvent(notebook, file_contents, notebook.path)
end

# Triggered when the local code has changed (user typed something),
# but the code hasn't run yet.
struct FileEditEvent <: PlutoEvent
    notebook::Notebook
    file_contents::String
    path::String
end

FileEditEvent(notebook::Notebook) = begin
    file_contents = sprint() do io
        # TODO: https://github.com/fonsp/Pluto.jl/pull/1729: serialize_temp flag
        # to only get local changes; the workspace edit of the notebook!
        save_notebook(io, notebook #=; serialize_temp=true =#)
    end
    FileEditEvent(notebook, file_contents, notebook.path)
end

# Triggered when we open a new notebook
struct NewNotebookEvent <: PlutoEvent
end

# Triggered when a user opens a notebook
struct OpenNotebookEvent <: PlutoEvent
    notebook::Notebook
end

# This will be fired ONLY if URL params don't match anything else.
# Useful if you want to create a file in a custom way,
# before opening the notebook
# Should return a redirect to /edit?id={notebook_id}
# Note: use Pluto.open - try_launch_notebook_response will return a fitting response.
struct CustomLaunchEvent <: PlutoEvent
    params::Dict{Any, Any}
    request::HTTP.Request
    try_launch_notebook_response::Function
end

# Triggered when a notebook has shut down.
struct ShutdownNotebookEvent <: PlutoEvent
    notebook::Notebook
end


text:
function run_commands(commands::Commands)
    Base.run(commands.cmd_1_before_script)
    Base.run(commands.cmd_2_script)
    Base.run(commands.cmd_3_after_script)
    return nothing
end

function generate_julia_pkg_commands(cfg::Config,
                                     temporary_environment::TemporaryEnvironment,
                                     proj_or_mnfst::ProjectFileOrManifestFile)
    julia_command = Base.julia_cmd()

    env_without_project = temporary_environment.env_without_project
    env_with_project    = temporary_environment.env_with_project

    before_script = _get_before_script(cfg, proj_or_mnfst)::String
    script        = _get_script(       cfg, proj_or_mnfst)::String
    after_script  = _get_after_script( cfg, proj_or_mnfst)::String

    cmd_1_before_script = setenv(`$(julia_command) -e $(before_script)`, env_without_project)
    cmd_2_script        = setenv(`$(julia_command) -e $(script)`,        env_with_project)
    cmd_3_after_script  = setenv(`$(julia_command) -e $(after_script)`,  env_with_project)

    commands = Commands(;
        cmd_1_before_script,
        cmd_2_script,
        cmd_3_after_script,
    )

    return commands
end

for field_name in [:before_script, :script, :after_script]
    fn = Symbol(:_get_, field_name)
    @eval begin
        function $(fn)(cfg::Config, proj_or_mnfst::ProjectFileOrManifestFile)::String
            value_from_proj_or_mnfst = proj_or_mnfst.$(field_name)
            if value_from_proj_or_mnfst !== nothing
                return value_from_proj_or_mnfst
            end
            value_from_cfg = cfg.$(field_name)
            return value_from_cfg
        end
    end
end


text:
#=
    This file includes all the funtions needed to simulate the DES
 =#

 mutable struct Options
     mode::String
     firstyear::Bool
     Options(; mode = "multithreads", firstyear = false) = new(mode, firstyear)
 end

# Simulate function
function simulate!(mg::Microgrid,
                   controller::AbstractController,
                   designer::AbstractDesigner,
                   ω_simu::Scenarios;
                   options::Options = Options())

    # Parameters
    ns = mg.parameters.ns

    if options.mode == "serial"
        # We simulate over the horizon for all the scenarios
        @showprogress for s in 1:ns
            simulate!(s, mg, controller, designer, ω_simu, options)
        end
    elseif options.mode == "multicores"
        # Init
        s = 0
        # We simulate over the horizon for all the scenarios in parallel
        @sync begin
            for p in workers()
                @async begin
                    while true
                        # Increment scenario
                        s += 1
                        # Break if ns is reached...
                        if s > ns
                            break
                        end
                        # Execute the function
                        remotecall_fetch(simulate!, p, s, mg, controller, designer, ω_simu, options)
                    end
                end
            end
        end
    elseif options.mode == "distributed"
        # We simulate over the horizon for all the scenarios in parallel using the distributed macro
        @sync @distributed for s in 1:ns
            simulate!(s, mg, controller, designer, ω_simu, options)
        end
    elseif options.mode == "multithreads"
        # We simulate over the horizon for all the scenarios in parallel using the distributed macro
        Threads.@threads for s in 1:ns
            simulate!(s, mg, controller, designer, ω_simu, options)
        end
    else
        println("Unknown mode... Please chose between 'serial', 'multicores', 'distributed' or 'multithreads'.")
    end
end
function simulate!(s::Int64,
                   mg::Microgrid,
                   controller::AbstractController,
                   designer::AbstractDesigner,
                   ω_simu::Scenarios,
                   options::Options)

    # Parameters
    nh = mg.parameters.nh
    ny = mg.parameters.ny

    # We simulate over the horizon for a single scenario
    for y in 1:ny
        simulate!(y, s, mg, controller, designer, ω_simu, options)
    end
end
function simulate!(y::Int64,
                   s::Int64,
                   mg::Microgrid,
                   controller::AbstractController,
                   designer::AbstractDesigner,
                   ω_simu::Scenarios,
                   options::Options)

    # Parameters
    nh = mg.parameters.nh

    if y == 1 && !options.firstyear
        # Update investment informations
        update_investment_informations!(y, s, mg, ω_simu)

        # Compute investment decision variables
        compute_investment_decisions!(y, s, mg, designer)

        # Compute investment dynamics
        compute_investment_dynamics!(y, s, mg, designer)
    else
        for h in 1:nh
            simulate!(h, y, s, mg, controller, designer, ω_simu, options)
        end

        # Update investment informations
        update_investment_informations!(y, s, mg, ω_simu)

        # Compute investment decision variables
        compute_investment_decisions!(y, s, mg, designer)

        # Compute investment dynamics
        compute_investment_dynamics!(y, s, mg, designer)
    end
end
function simulate!(h::Int64,
                   y::Int64,
                   s::Int64,
                   mg::Microgrid,
                   controller::AbstractController,
                   designer::AbstractDesigner,
                   ω_simu::Scenarios,
                   options::Options)

    # Update operation informations
    update_operation_informations!(h, y, s, mg, ω_simu)

    # Compute operation decision variables
    compute_operation_decisions!(h, y, s, mg, controller)

    # Compute operation dynamics
    compute_operation_dynamics!(h, y, s, mg, controller)

    # Power balance constraint checked for each node
    compute_power_balances!(h, y, s, mg)
end


text:
using Random: seed!

struct Bandit{A<:Integer, R<:AbstractFloat}
    actionspace::ActionSpace{A}
    rwdparams::Matrix{R}
    optsol::R
    function Bandit{A, R}(n=10, seed=1234) where {A<:Integer, R<:AbstractFloat}
        seed!(seed)
        rwdparams = randn(R, n, 2)
        new(ActionSpace(A(n)), rwdparams, maximum(rwdparams[:, 1]))
    end
end

Bandit() = Bandit{Int64, Float64}()

function step!(env::Bandit, action::Integer)
    reward = env.rwdparams[action, 1] + env.rwdparams[action, 2] * randn(eltype(env.rwdparams))
    regret = abs(env.optsol - env.rwdparams[action, 1])
    return reward, regret
end

reset!(::Bandit) = nothing


text:
# ==============================================================================
# Motor
#
# Part of the SLICOT benchmarks
#
# system type: continuous LTI system
# state dimension: 8
# input dimension: 2
# ==============================================================================
using SparseArrays

function state_matrix_motor()
    # system matrix
    I = [1, 2, 2, 3, 3, 3, 3, 4, 5, 6, 6, 7, 7, 7, 7, 8]
    J = [2, 3, 2, 1, 2, 3, 4, 1, 6, 7, 6, 5, 6, 7, 8, 5]
    V = [1, 8487.2, -1.0865, -2592.1, -21.119, -698.91, -141399.0, 1.0, 1.0,
         8487.2, -1.0865, -2592.1, -21.119, -698.91, -141399.0, 1.0]
    A = sparse(I, J, V)
end

function initial_state_motor()
    # initial set:
    # xᵢ ∈ [0.002, 0.0025] for i = 1
    # xᵢ ∈ [0.001, 0.0015] for i = 5
    # xᵢ = 0 otherwise
    X0 = Hyperrectangle(low=[0.002, 0.0, 0.0, 0.0, 0.001, 0.0, 0.0, 0.0],
                        high=[0.0025, 0.0, 0.0, 0.0, 0.0015, 0.0, 0.0, 0.0])
end

function motor_homog()
    A = state_matrix_motor()
    X0 = initial_state_motor()

    prob = @ivp(x' = Ax, x(0) ∈ X0) # continuous LTI system
    tspan = (0.0, 20.0)
    return prob, tspan
end

function motor()
    A = state_matrix_motor()
    B = sparse([4, 8], [1, 2], [-1.0, -1.0]) # input matrix
    X = Universe(8) # state domain: unconstrained
    U = Hyperrectangle([0.23, 0.3], [0.07, 0.1]) # input domain
    X0 = initial_state_motor()

    prob = @ivp(x' = Ax + Bu, x ∈ X, u ∈ U, x(0) ∈ X0) # continuous LTI system
    tspan = (0.0, 20.0)
    return prob, tspan
end


text:
module Render
using Base: get
export render, @format, format

@nospecialize
struct Discard end
const discard = Discard()

function render(expr, pair :: Pair{Symbol, Any})
    render(expr, Dict(pair))
end

function render(expr, config :: Dict{Symbol, Any})
    function visit(expr :: Expr)
        hd = expr.head
        tl = filter(x -> x !== discard, map(visit, expr.args))
        Expr(hd, tl...)
    end

    function visit(sym :: Symbol)
        get(config, sym) do
            sym
        end
    end

    function visit(:: LineNumberNode)
        discard
    end

    function visit(a)
        a
    end

    visit(expr)
end

function format(args, template)
    function dispatch(arg :: Symbol)
        Expr(:call, :(=>), QuoteNode(arg), arg)
    end

    function dispatch(arg :: Pair)
        Expr(:call, :(=>), QuoteNode(arg[1]), arg[2])
    end

    function dispatch(arg :: Expr)
        @assert arg.head == :(=)
        sym = arg.args[1]
        @assert sym isa Symbol "$sym"
        value = arg.args[2]
        Expr(:call, :(=>), QuoteNode(sym), value)
    end
    function  dispatch(_)
        throw("Unknown argtype")
    end

    constlist = map(dispatch, args.args)
    constlist = Expr(:vect, constlist...)
    config = Expr(:call, Dict{Symbol, Any}, constlist)

    wrap = @static VERSION < v"1.1.0" ? (x -> x) : (x -> Expr(:call, merge, :(Base.@locals), x))
    Expr(:call, render, template, wrap(config))
end

function format(template)
    Expr(:call, render, template, :(Base.@locals))
end

macro format(args, template)
    esc(format(args, template))
end

macro format(template)
    esc(format(template))
end

for deprecated_sym in [:format, :render, Symbol("@format")]
   @eval Base.@deprecate_moved $deprecated_sym "ASTRender.jl"
end
@specialize

end


text:
#Generative algorithms

import Base: integer_partitions

export cool_lex, ncpartitions

# Lists the partitions of the number n, the order is consistent with GAP
function integer_partitions(n::Integer)
    if n < 0
        throw(DomainError())
    elseif n == 0
        return Vector{Int}[]
    elseif n == 1
        return Vector{Int}[[1]]
    end

    list = Vector{Int}[]

    for p in integer_partitions(n-1)
        push!(list, [p, 1])
        if length(p) == 1 || p[end] < p[end-1]
            push!(list, [p[1:end-1], p[end]+1])
        end
    end

    list
end

# Produces (n,k)-combinations in cool-lex order
#
#Implements the cool-lex algorithm to generate (n,k)-combinations
#@article{Ruskey:2009fk,
#	Author = {Frank Ruskey and Aaron Williams},
#	Doi = {10.1016/j.disc.2007.11.048},
#	Journal = {Discrete Mathematics},
#	Month = {September},
#	Number = {17},
#	Pages = {5305-5320},
#	Title = {The coolest way to generate combinations},
#	Url = {http://www.sciencedirect.com/science/article/pii/S0012365X07009570},
#	Volume = {309},
#	Year = {2009}}
function cool_lex(n::Integer, t::Integer)
  s = n-t
  if n > 64 error("Not implemented for n > 64") end
  if t < 1 error("Need t>1") end
  R0::Int64=R1::Int64=0
  R2::Int64 = 1 << (s+t)
  R3::Int64 = (1 << t) - 1
  while R3 & R2 == 0
    produce(visit(R3))
    R0 = R3 & (R3 + 1)
    R1 = R0 $ (R0 - 1)
    R0 = R1 + 1
    R1 &= R3
    R0 = max((R0 & R3) - 1, 0)
    R3 += R1 - R0
  end
end

#Converts an integer word X into a subset
#If X & 2^k == 1, then k is in the subset
function visit(X::Integer)
  subset = Int[]
  n::Int=1
  while X != 0
    if X & 1 == 1 push!(subset, n) end
    X >>= 1
    n += 1
  end
  subset
end

#Produces noncrossing partitions of length n
ncpartitions(n::Integer)=ncpart(1,n,n,{})
function ncpart(a::Integer, b::Integer, nn::Integer,
    x::Array{Any,1})
  n=b-a+1
  for k=1:n
    for root in @task cool_lex(n, k)
      root += a-1
      #Abort if construction is out of lex order
      if length(x)>0 && x[end] > root return end
      #Produce if we've filled all the holes
      sofar = {x..., root}
      ssofaru = sort(union(sofar...))
      if length(ssofaru)==nn && ssofaru==[1:nn]
        produce(sofar)
        return
      end
      #otherwise patch all remaining holes
      blob = [ssofaru; nn+1]
      for l=1:length(blob)-1
        ap, bp = blob[l]+1, blob[l+1]-1
        if ap <= bp ncpart(ap, bp, nn, sofar) end
      end
    end
  end
end


text:
# white.jl
# White and Breusch-Pagan tests for heteroskedasticity
#
# Copyright (C) 2020   Paul Söderlind
#
# Permission is hereby granted, free of charge, to any person obtaining
# a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including
# without limitation the rights to use, copy, modify, merge, publish,
# distribute, sublicense, and/or sell copies of the Software, and to
# permit persons to whom the Software is furnished to do so, subject to
# the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
# LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
# OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
# WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

export WhiteTest, BreuschPaganTest

struct WhiteTest <: HypothesisTest
    dof::Int                 #degrees of freedom in test
    lm::Float64              #Lagrange Multiplier test statistic, distributed as Chisq(dof)
    type::Symbol
end

"""
    WhiteTest(X, e; type = :White)

Compute White's (or Breusch-Pagan's) test for heteroskedasticity.

`X` is a matrix of regressors and `e` is the vector of residuals from the original model.
The keyword argument `type` is either `:linear` for the Breusch-Pagan/Koenker test,
`:linear_and_squares` for White's test with linear and squared terms only (no cross-products), or
`:White` (the default) for the full White's test (linear, squared and cross-product terms).
`X` should include a constant and at least one more regressor, with observations in rows and
regressors in columns. In some applications, `X` is a subset of the regressors in
the original model, or just the fitted values. This saves degrees of freedom and may
give a more powerful test. The `lm` (Lagrange multiplier) test statistic is T*R2
where R2 is from the regression of `e^2` on the terms mentioned above.
Under the null hypothesis it is distributed as `Chisq(dof)` where `dof` is the number
of independent terms (not counting the constant), so the null is rejected when the
test statistic is large enough.

Implements: [`pvalue`](@ref)

# References
  * H. White, (1980): A heteroskedasticity-consistent covariance matrix estimator and a direct test for heteroskedasticity, Econometrica, 48, 817-838.
  * T.S. Breusch & A.R. Pagan (1979), A simple test for heteroscedasticity and random coefficient variation, Econometrica, 47, 1287-1294
  * R. Koenker (1981), A note on studentizing a test for heteroscedasticity, Journal of Econometrics, 17, 107-112

# External links
  * [White's test on Wikipedia](https://en.wikipedia.org/wiki/White_test)
"""
function WhiteTest(X::AbstractMatrix{<:Real}, e::AbstractVector{<:Real}; type = :White)
    (n,K) = (size(X,1),size(X,2))    #nobs,nvars
    n == length(e) || throw(DimensionMismatch("inputs must have the same length"))

    K >= 2 || throw(ArgumentError("X must have >= 2 columns"))

    intercept_col = false
    for i = 1:K
        col           = view(X,:,i)
        intercept_col = first(col) != 0 && all(==(first(col)), col)
        intercept_col && break
    end
    intercept_col || throw(ArgumentError("one of the colums of X must be a non-zero constant"))

    if type == :linear                    #Breush-Pagan/Koenker
        z = X
    elseif type == :linear_and_squares    #White with linear and squares
        z = [X X.^2]
    else                                      #White with linear, squares and cross-products
        z = fill(NaN,n,round(Int,K*(K+1)/2))
        vv = 1
        @views for i = 1:K, j = i:K
            z[:,vv] .= X[:,i] .* X[:,j]       #eg. x1*x1, x1*x2, x2*x2
            vv      += 1
        end
    end

    dof = rank(z) - 1                        #number of independent regressors in z
    e2  = e.^2
    b   = z\e2
    res = e2 - z*b
    R2  = 1 - var(res)/var(e2)
    lm  = n*R2
    return WhiteTest(dof, lm, type)
end

"""
    BreuschPaganTest(X, e)

Compute Breusch-Pagan's test for heteroskedasticity.

`X` is a matrix of regressors from the original model and `e` the vector of residuals.
This is equivalent to `WhiteTest(X, e, type = :linear)`.
See [`WhiteTest`](@ref) for further details.
"""
BreuschPaganTest(X, e) = WhiteTest(X, e, type = :linear)

testname(t::WhiteTest) = "White's (or Breusch-Pagan's) test for heteroskedasticity"

population_param_of_interest(t::WhiteTest) = ("T*R2", 0, t.lm)
default_tail(test::WhiteTest)              = :right

function show_params(io::IO, t::WhiteTest, ident = "")
    println(io, ident, "T*R^2 statistic:        ", t.lm)
    println(io, ident, "degrees of freedom:     ", t.dof)
    println(io, ident, "type:                   ", t.type)
end

StatsBase.dof(t::WhiteTest) = t.dof
pvalue(t::WhiteTest) = pvalue(Chisq(t.dof), t.lm, tail=:right)


text:
# Author: Youngjun Kim, youngjun@stanford.edu
# Date: 06/16/2014

module Sensor

export
    AbstractSensor,

    initialize,
    update,

    updateSensor,

    SimpleTCASSensor,
    SimpleTCASSRState,
    SimpleTCASSensorInput,
    SimpleTCASSensorOutput,

    ACASXSensor,
    ACASXSensorState,
    ACASXSensorInput,
    ACASXSensorOutput

import CommonInterfaces: initialize, update

using AbstractSensorImpl

using SimpleTCASSensorImpl

using ACASXSensorImpl

end


text:
using Base.Test

@testset "Scaling" begin
  include("qpTestProblems.jl")
  include("socp_TestProblems.jl")
end


text:
function thread_blocks_heuristic(len::Integer)
    # TODO better threads default
    threads = min(len, 256)
    blocks = ceil(Int, len / threads)
    (blocks,), (threads,)
end


text:
# Integration test combining AtomsBase, NBodySimulator, and InteratomicPotentials

@testset "Integration: AB & NBS & IP" begin
    # TODO
end


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

mutable struct ComputeLongRunningOperationProperties <: SwaggerModel
    output::Any # spec type: Union{ Nothing, Any } # spec name: output

    function ComputeLongRunningOperationProperties(;output=nothing)
        o = new()
        validate_property(ComputeLongRunningOperationProperties, Symbol("output"), output)
        setfield!(o, Symbol("output"), output)
        o
    end
end # type ComputeLongRunningOperationProperties

const _property_map_ComputeLongRunningOperationProperties = Dict{Symbol,Symbol}(Symbol("output")=>Symbol("output"))
const _property_types_ComputeLongRunningOperationProperties = Dict{Symbol,String}(Symbol("output")=>"Any")
Base.propertynames(::Type{ ComputeLongRunningOperationProperties }) = collect(keys(_property_map_ComputeLongRunningOperationProperties))
Swagger.property_type(::Type{ ComputeLongRunningOperationProperties }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_ComputeLongRunningOperationProperties[name]))}
Swagger.field_name(::Type{ ComputeLongRunningOperationProperties }, property_name::Symbol) =  _property_map_ComputeLongRunningOperationProperties[property_name]

function check_required(o::ComputeLongRunningOperationProperties)
    true
end

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


text:
function pretty_print(constraints::Constraints)
    @printf("==========================================\n")
    @printf("Constraints\n")
    @printf("* Lower Bounds:\n")
    for c in 'a':'z'
        if constraints.lower_bound[c] > 0
            @printf(" * %s ≥ %d\n", c, constraints.lower_bound[c])
        end
    end
    @printf("\n")

    @printf("* Exact Matches:\n")
    @printf(" * ")
    for i in 1:5
        if constraints.is_exact_match[i]
            print(constraints.match[i])
        else
            print('·')
        end
    end
    @printf("\n\n")

    @printf("* Global Ban:\n")
    for c in 'a':'z'
        if constraints.global_ban[c]
            @printf(" * %s = 0\n", c)
        end
    end
    @printf("\n")

    @printf("* Local Ban:\n")
    for i in 1:5
        for c in 'a':'z'
            if constraints.local_ban[i][c]
                @printf(" * word[%d] ≠ %s\n", i, c)
            end
        end
    end
    @printf("==========================================\n\n")
end


text:
using Plots

pa = readdlm("pa")
pa = reshape(pa', length(pa))
i = readdlm("i")
i = reshape(i', length(i))
q = readdlm("q")
q = reshape(q', length(q))
u = readdlm("u")
u = reshape(u', length(u))
v = readdlm("v")
v = reshape(v', length(v))

plot(pa,i)

max_stars_repo_name: KajWiik/PolarimetryCalibration

max_issues_repo_name: KajWiik/PolarimetryCalibration

max_forks_repo_name: KajWiik/PolarimetryCalibration


text:
function ordinal(n::Integer)
    n < 0 && throw(DomainError())
    suffixes = ("st", "nd", "rd")
    u = n % 10
    t = n ÷ 10 % 10
    if u > 3 || u == 0 || t == 1
        suf = "th"
    else
        suf = suffixes[u]
    end
    return string(n, suf)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using ODEInterfaceDiffEq, DiffEqProblemLibrary, DiffEqBase
using Test

@time @testset "Algorithms" begin include("algorithm_tests.jl") end
@time @testset "Saving" begin include("saving_tests.jl") end
@time @testset "Mass Matrix" begin include("mass_matrix_tests.jl") end
@time @testset "Jacobian Tests" begin include("jac_tests.jl") end
@time @testset "Callback Tests" begin include("callbacks.jl") end


text:
function gc2gc!(dst_ceps::AbstractVector, dst_γ,
                src_ceps::AbstractVector, src_γ)
    fill!(dst_ceps, zero(eltype(dst_ceps)))
    dst_order = length(dst_ceps)-1
    m1 = length(src_ceps)-1
    dst_ceps[1] = src_ceps[1]

    for m=2:dst_order+1
        ss1, ss2 = 0.0, 0.0
        min = (m1 < m-1) ? m1 : m - 2

        for k=2:min+1
            @inbounds cc = src_ceps[k] * dst_ceps[m-k+1]
            ss2 += (k-1) * cc
            ss1 += (m-k) * cc
        end

        if m <= m1+1
            @inbounds dst_ceps[m] = src_ceps[m] + (dst_γ*ss2 - src_γ*ss1)/(m-1)
        else
            @inbounds dst_ceps[m] = (dst_γ*ss2 - src_γ*ss1)/(m-1)
        end
    end

    dst_ceps
end

function gc2gc(src_ceps::AbstractVector, src_γ=0.0,
               dst_order=length(src_ceps)-1, dst_γ=0.0)
    gc2gc!(Array(eltype(src_ceps), dst_order+1), dst_γ, src_ceps, src_γ)
end

function gc2gc{T<:MelGeneralizedCepstrum}(state::SpectralParamState{T},
                                          dst_order=param_order(paramdef(state)),
                                          dst_γ=glog_gamma(paramdef(state)))
    assert_not_ready_to_filt(state)

    def = paramdef(state)
    order = param_order(def)
    src_α = allpass_alpha(def)
    src_γ = glog_gamma(def)
    data = rawdata(state)
    normalized = (dst_γ == zero(dst_γ))
    newdef = MelGeneralizedCepstrum(order, src_α, dst_γ)
    SpectralParamState(newdef, gc2gc(data, src_γ, dst_order, dst_γ),
                       has_loggain(state), normalized)
end


text:
AbsoluteRectangle{N,T}(mini::Vec{N,T}, maxi::Vec{N,T}) = HyperRectangle{N,T}(mini, maxi-mini)

@compat (::Type{AABB})(a) = AABB{Float32}(a)
@compat function (B::Type{AABB{T}}){T}(a::Pyramid)
    w,h = a.width/T(2), a.length
    m = Vec{3,T}(a.middle)
    B(m-Vec{3,T}(w,w,0), m+Vec{3,T}(w, w, h))
end
@compat (B::Type{AABB{T}}){T}(a::Cube) = B(origin(a), widths(a))
@compat (B::Type{AABB{T}}){T}(a::AbstractMesh) = B(vertices(a))
@compat (B::Type{AABB{T}}){T}(a::NativeMesh) = B(gpu_data(a.data[:vertices]))

@compat function (B::Type{AABB{T}}){T}(
        positions, scale, rotation,
        primitive::AABB{T}
    )

    ti = TransformationIterator(positions, scale, rotation)
    B(ti, primitive)
end
@compat function (B::Type{AABB{T}}){T}(instances::Instances)
    ti = TransformationIterator(instances)
    B(ti, B(instances.primitive))
end

function transform(translation, scale, rotation, points)
    _max = Vec3f0(typemin(Float32))
    _min = Vec3f0(typemax(Float32))
    for p in points
        x = scale.*Vec3f0(p)
        x = Vec3f0(rotation*Vec4f0(x, 1f0))
        x = Vec3f0(translation)+x
        _min = min(_min, x)
        _max = max(_max, x)
    end
    AABB{Float32}(_min, _max-_min)
end

@compat function (B::Type{AABB{T}}){T}(
      ti::TransformationIterator, primitive::AABB{T}
    )
    state = start(ti)
    if done(ti, state)
        return primitive
    end
    trans_scale_rot, state = next(ti, state)
    points = decompose(Point3f0, primitive)
    bb = transform(trans_scale_rot..., points)
    while !done(ti, state)
        trans_scale_rot, state = next(ti, state)
        translatet_bb = transform(trans_scale_rot..., points)
        bb = union(bb, translatet_bb)
    end
    bb
end


text:
using Test

using MARY_fRG.Fermi
using MARY_fRG.FlowEquation
using MARY_fRG.Drawers

@testset "切分费米面" begin
    qs = common_square_lattice(0.20)
    Γ4 = TFGamma4(qs, 8.0, 16, 4)
    plt = draw_lines(Γ4.model.brillouin.edges)
    draw_polygon!(plt, Γ4.ltris; pcidx=Γ4.lpats)
    @savepng plt "brlu-20"
    newΓ4 = TFGamma4_reweight_ltris_mt(Γ4, 0.0, 1.0, 4)
    plt = draw_lines(newΓ4.model.brillouin.edges)
    draw_polygon!(plt, newΓ4.ltris; pcidx=newΓ4.lpats)
    @savepng plt "brlu-201"
end


text:
# Load Julia packages (libraries) needed

using StatisticalRethinking
#gr(size=(600,600));

# ### snippet 2.3

# Define a grid

N = 20
p_grid = range( 0 , stop=1 , length=N )

# Define the (uniform) prior

prior = ones( 20 );

# Compute likelihood at each value in grid

likelihood = [pdf(Binomial(9, p), 6) for p in p_grid]
likelihood[1:5]

# Compute product of likelihood and prior

unstd_posterior = likelihood .* prior;

# Standardize the posterior, so it sums to 1

posterior = unstd_posterior  ./ sum(unstd_posterior);

# ### snippet 2.4

p1 = plot( p_grid , posterior ,
    xlab="probability of water" , ylab="posterior probability",
    lab = "interpolated", title="20 points" )
p2 = scatter!( p1, p_grid , posterior, lab="computed" )

# ### snippet 2.5

prior1 = [p < 0.5 ? 0 : 1 for p in p_grid]
prior2 = [exp( -5*abs( p - 0.5 ) ) for p in p_grid]

p3 = plot(p_grid, prior1,
  xlab="probability of water" , ylab="posterior probability",
  lab = "semi_uniform", title="Other priors" )
scatter!(p3, p_grid, prior1, lab = "semi_uniform grid point")
plot!(p3, p_grid, prior2,  lab = "double_exponential" )
scatter!(p3, p_grid, prior2,  lab = "double_exponential grid point" )

# End of `02/clip-03-05.jl`


text:
const MessyData = Union{AbstractMatrix, DataFrame}

mutable struct GGLRM{T <: MessyData} <: AbstractGLRM
  A::T                         # The data table
  losses::Array{Loss,1}        # array of loss functions
  rx::Regularizer              # Regularizer to apply to each row of X
  ry::Regularizer              # Array of regularizers to be applied to each column of Y
  k::Int                       # Desired rank
  observed_features::ObsArray  # for each example, an array telling which features were observed
  observed_examples::ObsArray  # for each feature, an array telling in which examples the feature was observed
  X::AbstractArray{Float64,2}  # Representation of data in low-rank space. A ≈ X'Y
  Y::AbstractArray{Float64,2}  # Representation of features in low-rank space. A ≈ X'Y
end

function no_multidim_edges(ry::AbstractGraphReg, losses::Array)
  for edge in edges(ry.idxgraph.graph)
    i = edge.src; j = edge.dst
    if (embedding_dim(losses[i]) > 1) || (embedding_dim(losses[j]) > 1)
      error("Graph regularizer cannot have any edges into or out of a multidimensional loss")
    end
  end
end

function GGLRM(A::MessyData, losses::Array, rx::Regularizer, ry::Regularizer, k::Int;
          X::Matrix = randn(k, size(A,1)), Y::Matrix = randn(k,embedding_dim(losses)),
          obs = nothing,                                    # [(i₁,j₁), (i₂,j₂), ... (iₒ,jₒ)]
          observed_features = fill(1:size(A,2), size(A,1)), # [1:n, 1:n, ... 1:n] m times
          observed_examples = fill(1:size(A,1), size(A,2)), # [1:m, 1:m, ... 1:m] n times)# [(i₁,j₁), (i₂,j₂), ... (iₒ,jₒ)]
          offset::Bool=false, scale::Bool=false,
          sparse_na::Bool=true)

  n, d = size(A)
  if isa(rx, AbstractGraphReg)
    #A graph regularizer on X should be able to multiply with X
    if (n,n) != size(matrix(rx))
      error("Graph regularizer on X must have as many nodes as there are rows in X")
    end
  end

  if isa(ry, AbstractGraphReg)
    #Same goes for Y
    if (d,d) != size(matrix(ry))
      error("Graph regularizer on Y must have as many nodes as there are columns in A")
    end
    no_multidim_edges(ry, losses)
    #ry = GraphQuadReg(embed_graph(ry.idxgraph, get_yidxs(losses)), ry.scale, ry.quadamt)
    ry = embed(ry, get_yidxs(losses))
  end

  if size(X) != (k, size(A,1))
    error("X must have the same number of columns as there are rows in A. This is the transpose of the monograph's description but makes for efficient memory access")
  end

  if size(Y)!=(k,sum(map(embedding_dim, losses)))
    error("Y must be of size (k,d) where d is the sum of the embedding dimensions of all the losses. \n(1 for real-valued losses, and the number of categories for categorical losses).")
  end

  # Determine observed entries of data
  if obs==nothing && sparse_na && isa(A,SparseMatrixCSC)
    obs = [zip(findn(A)...)...] # observed indices (list of tuples)
  end
  if obs==nothing # if no specified array of tuples, use what was explicitly passed in or the defaults (all)
  # println("no obs given, using observed_features and observed_examples")
    glrm = GGLRM(A,losses,rx,ry,k, observed_features, observed_examples, X,Y)
  else
    glrm = GGLRM(A,losses,rx,ry,k, sort_observations(obs,size(A)...)..., X,Y)
  end

  if (offset) add_offset!(glrm) end
  if (scale) equilibrate_variance!(glrm) end

  glrm
end

function GGLRM(A::MessyData, loss::Loss, rx::Regularizer, ry::Regularizer, k::Int; kwargs...)
  losses = convert(Array{Loss, 1}, [copy(loss) for i in 1:size(A,2)])
  GGLRM(A, losses, rx, ry, k; kwargs...)
end