text:
"""
Subroutine dcauchy

This subroutine computes a Cauchy step that satisfies a trust
region constraint and a sufficient decrease condition.

Ths Cauchy step is computed for the quadratic

    q(s) = 0.5*s'*A*s + g'*s

where A is a symmetric matrix in compressed row storage, and
g is a vector. Given a parameter alpha, the Cauchy step is

    s[alpha] = P[x - alpha*g] - x,

with P the projection onto the n-dimensional interval [xl,xu].
The Cauchy step satisfies the trust region constraint and the
sufficient decrease condition

    || s || <= delta,    q(s) <= mu_0*(g'*s),

where mu_0 is a constant in (0,1).

MINPACK-2 Project. March 1999.
Argonne National Laboratory.
Chih-Jen Lin and Jorge J. More'.
"""
function dcauchy(n,x,xl,xu,A,
                         g,delta,alpha,s,wa)
    p5 = 0.5
    one = 1.0

    # Constant that defines sufficient decrease.

    mu0 = 0.01

    # Interpolation and extrapolation factors.

    interpf = 0.1
    extrapf = 10.0

    # Find the minimal and maximal break-point on x - alpha*g.

    dcopy(n,g,1,wa,1)
    dscal(n,-one,wa,1)
    nbrpt,brptmin,brptmax = dbreakpt(n,x,xl,xu,wa)

    # Evaluate the initial alpha and decide if the algorithm
    # must interpolate or extrapolate.

    dgpstep(n,x,xl,xu,-alpha,g,s) # s = P(x - alpha*g) - x
    if dnrm2(n,s,1) > delta
        interp = true
    else
        dssyax(n, A, s, wa)
        gts = ddot(n,g,1,s,1)
        q = p5*ddot(n,s,1,wa,1) + gts
        interp = (q >= mu0*gts)
    end

    # Either interpolate or extrapolate to find a successful step.

    if interp

        # Reduce alpha until a successful step is found.

        search = true
        while search

            # This is a crude interpolation procedure that
            # will be replaced in future versions of the code.

            alpha = interpf*alpha
            dgpstep(n,x,xl,xu,-alpha,g,s)
            if dnrm2(n,s,1) <= delta
                dssyax(n, A, s,wa)
                gts = ddot(n,g,1,s,1)
                q = p5*ddot(n,s,1,wa,1) + gts
                search = (q > mu0*gts)
            end
        end

    else

        # Increase alpha until a successful step is found.

        search = true
        alphas = alpha
        while (search && alpha <= brptmax)

            # This is a crude extrapolation procedure that
            # will be replaced in future versions of the code.

            alpha = extrapf*alpha
            dgpstep(n,x,xl,xu,-alpha,g,s)
            if dnrm2(n,s,1) <= delta
                dssyax(n, A, s, wa)
                gts = ddot(n,g,1,s,1)
                q = p5*ddot(n,s,1,wa,1) + gts
                if q < mu0*gts
                    search = true
                    alphas = alpha
                end
            else
                search = false
            end
        end

        # Recover the last successful step.

        alpha = alphas
        dgpstep(n,x,xl,xu,-alpha,g,s)
    end

    return alpha
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 PublicIPAddressPropertiesFormat2 <: SwaggerModel
    publicIPAllocationMethod::Any # spec type: Union{ Nothing, IPAllocationMethod } # spec name: publicIPAllocationMethod
    publicIPAddressVersion::Any # spec type: Union{ Nothing, IPVersion } # spec name: publicIPAddressVersion
    ipConfiguration::Any # spec type: Union{ Nothing, IPConfiguration } # spec name: ipConfiguration
    dnsSettings::Any # spec type: Union{ Nothing, PublicIPAddressDnsSettings2 } # spec name: dnsSettings
    ddosSettings::Any # spec type: Union{ Nothing, DdosSettings2 } # spec name: ddosSettings
    ipTags::Any # spec type: Union{ Nothing, Vector{IpTag2} } # spec name: ipTags
    ipAddress::Any # spec type: Union{ Nothing, String } # spec name: ipAddress
    publicIPPrefix::Any # spec type: Union{ Nothing, SubResource } # spec name: publicIPPrefix
    idleTimeoutInMinutes::Any # spec type: Union{ Nothing, Int32 } # spec name: idleTimeoutInMinutes
    resourceGuid::Any # spec type: Union{ Nothing, String } # spec name: resourceGuid
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState

    function PublicIPAddressPropertiesFormat2(;publicIPAllocationMethod=nothing, publicIPAddressVersion=nothing, ipConfiguration=nothing, dnsSettings=nothing, ddosSettings=nothing, ipTags=nothing, ipAddress=nothing, publicIPPrefix=nothing, idleTimeoutInMinutes=nothing, resourceGuid=nothing, provisioningState=nothing)
        o = new()
        validate_property(PublicIPAddressPropertiesFormat2, Symbol("publicIPAllocationMethod"), publicIPAllocationMethod)
        setfield!(o, Symbol("publicIPAllocationMethod"), publicIPAllocationMethod)
        validate_property(PublicIPAddressPropertiesFormat2, Symbol("publicIPAddressVersion"), publicIPAddressVersion)
        setfield!(o, Symbol("publicIPAddressVersion"), publicIPAddressVersion)
        validate_property(PublicIPAddressPropertiesFormat2, Symbol("ipConfiguration"), ipConfiguration)
        setfield!(o, Symbol("ipConfiguration"), ipConfiguration)
        validate_property(PublicIPAddressPropertiesFormat2, Symbol("dnsSettings"), dnsSettings)
        setfield!(o, Symbol("dnsSettings"), dnsSettings)
        validate_property(PublicIPAddressPropertiesFormat2, Symbol("ddosSettings"), ddosSettings)
        setfield!(o, Symbol("ddosSettings"), ddosSettings)
        validate_property(PublicIPAddressPropertiesFormat2, Symbol("ipTags"), ipTags)
        setfield!(o, Symbol("ipTags"), ipTags)
        validate_property(PublicIPAddressPropertiesFormat2, Symbol("ipAddress"), ipAddress)
        setfield!(o, Symbol("ipAddress"), ipAddress)
        validate_property(PublicIPAddressPropertiesFormat2, Symbol("publicIPPrefix"), publicIPPrefix)
        setfield!(o, Symbol("publicIPPrefix"), publicIPPrefix)
        validate_property(PublicIPAddressPropertiesFormat2, Symbol("idleTimeoutInMinutes"), idleTimeoutInMinutes)
        setfield!(o, Symbol("idleTimeoutInMinutes"), idleTimeoutInMinutes)
        validate_property(PublicIPAddressPropertiesFormat2, Symbol("resourceGuid"), resourceGuid)
        setfield!(o, Symbol("resourceGuid"), resourceGuid)
        validate_property(PublicIPAddressPropertiesFormat2, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type PublicIPAddressPropertiesFormat2

const _property_map_PublicIPAddressPropertiesFormat2 = Dict{Symbol,Symbol}(Symbol("publicIPAllocationMethod")=>Symbol("publicIPAllocationMethod"), Symbol("publicIPAddressVersion")=>Symbol("publicIPAddressVersion"), Symbol("ipConfiguration")=>Symbol("ipConfiguration"), Symbol("dnsSettings")=>Symbol("dnsSettings"), Symbol("ddosSettings")=>Symbol("ddosSettings"), Symbol("ipTags")=>Symbol("ipTags"), Symbol("ipAddress")=>Symbol("ipAddress"), Symbol("publicIPPrefix")=>Symbol("publicIPPrefix"), Symbol("idleTimeoutInMinutes")=>Symbol("idleTimeoutInMinutes"), Symbol("resourceGuid")=>Symbol("resourceGuid"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_PublicIPAddressPropertiesFormat2 = Dict{Symbol,String}(Symbol("publicIPAllocationMethod")=>"IPAllocationMethod", Symbol("publicIPAddressVersion")=>"IPVersion", Symbol("ipConfiguration")=>"IPConfiguration", Symbol("dnsSettings")=>"PublicIPAddressDnsSettings2", Symbol("ddosSettings")=>"DdosSettings2", Symbol("ipTags")=>"Vector{IpTag2}", Symbol("ipAddress")=>"String", Symbol("publicIPPrefix")=>"SubResource", Symbol("idleTimeoutInMinutes")=>"Int32", Symbol("resourceGuid")=>"String", Symbol("provisioningState")=>"ProvisioningState")
Base.propertynames(::Type{ PublicIPAddressPropertiesFormat2 }) = collect(keys(_property_map_PublicIPAddressPropertiesFormat2))
Swagger.property_type(::Type{ PublicIPAddressPropertiesFormat2 }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_PublicIPAddressPropertiesFormat2[name]))}
Swagger.field_name(::Type{ PublicIPAddressPropertiesFormat2 }, property_name::Symbol) =  _property_map_PublicIPAddressPropertiesFormat2[property_name]

function check_required(o::PublicIPAddressPropertiesFormat2)
    true
end

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

ext: jl


text:
using Test
using TensorFactorizations
using TensorOperations
using Random
using ArgParse
using LinearAlgebra

function parse_pars()
    settings = ArgParseSettings(autofix_names=true)
    @add_arg_table(settings
        , "--iters", arg_type=Int, default=50
        , "--no_test_svd", action=:store_true
        , "--no_test_eig", action=:store_true
        , "--no_test_split", action=:store_true
    )
    pars = parse_args(ARGS, settings; as_symbols=true)
    pars[:test_svd] = !pars[:no_test_svd]
    pars[:test_eig] = !pars[:no_test_eig]
    pars[:test_split] = !pars[:no_test_split]
    return pars
end

""" Generate a random shape. """
function rshape(;n=nothing, chi=nothing, nlow=0, nhigh=5, chilow=1, chihigh=5)
    # n is the number of legs, chi is the bond dimension(s)
    n == nothing && (n = rand(nlow:nhigh))
    shape = (chi == nothing ? rand(chilow:chihigh, (n,))
             : collect(repeat([chi], n)))
    return shape
end

""" Generate a random tensor, either with a given or a random shape. """
function rtensor(;shape=nothing, n=nothing, chi=nothing,
                 nlow=0, nhigh=5, chilow=0, chihigh=6, dtype=nothing)
    if dtype == nothing
        dtype = rand([Int, Float64, Complex{Float64}])
    end
    if shape == nothing
        shape = rshape(n=n, chi=chi, nlow=nlow, nhigh=nhigh, chilow=chilow,
                       chihigh=chihigh)
    end
    A = rand(dtype, (shape...),)
    return A
end

function test_svd()
    shp = rshape(nlow=2)
    n = length(shp)
    perm = randperm(n)
    I = rand(1:n-1)
    J = n - I
    ins = perm[1:I]
    outs = perm[I+1:end]
    A = rtensor(shape=shp)

    # SVD and reconstruct without truncation
    U, S, Vt = tensorsvd(A, ins, outs)
    US = tensorcontract(U, [-1*(1:I); 1], Diagonal(S), [1,-I-1])
    Areco = tensorcontract(US, [-1*(1:I); 1], Vt, [1; -1*(I+1:n)])
    Areco = tensorcopy(Areco, collect(1:ndims(Areco)), invperm(perm))
    @test isapprox(A, Areco)

    # SVD and reconstruct with truncation
    eps = 10.0^(-rand(0:3))
    dim_in = prod(shp[perm][1:I])
    dim_out = prod(shp[perm][I+1:end])
    max_dim = min(dim_in, dim_out)
    max_chi = rand(1:max_dim)
    chis = 1:max_chi
    U, S, Vt, error = tensorsvd(A, ins, outs, eps=eps, chis=chis,
                                return_error=true)
    # This test could fail if there were degenerate singular values. We assume
    # this is not the case because the tensors are random.
    @test (error<eps || length(S) == max_chi)
    US = tensorcontract(U, [-1*(1:I); 1], Diagonal(S), [1,-I-1])
    Areco = tensorcontract(US, [-1*(1:I); 1], Vt, [1; -1*(I+1:n)])
    Areco = tensorcopy(Areco, collect(1:ndims(Areco)), invperm(perm))
    Anorm = norm(A)
    true_error = Anorm > 0 ? norm(Areco - A)/Anorm : norm(Areco)
    @test isapprox(error, true_error; atol=1e-8)
end

function test_eig()
    # Non-hermitian
    I = rand(1:3)
    n = 2*I
    shp = rshape(n=n)
    perm = randperm(n)
    ins = perm[1:I]
    outs = perm[I+1:end]
    # Making sure that the left and right legs match.
    shp = collect(shp)
    for (o,i) in zip(outs, ins)
        shp[o] = shp[i]
    end
    shp = (shp...,)
    A = rtensor(shape=shp)

    # Check that eigenvectors really are eigenvectors
    eps = 10.0^(-rand(0:3))
    dim_in = prod(shp[perm][1:I])
    dim_out = prod(shp[perm][I+1:end])
    max_dim = min(dim_in, dim_out)
    max_chi = rand(1:max_dim)
    chis = 1:max_chi
    E, U, error = tensoreig(A, ins, outs, eps=eps, chis=chis,
                            return_error=true, break_degenerate=true)
    @test (error<eps || length(E) == max_chi)
    A_contract_list = collect(repeat([0], n))
    for (k,i) in enumerate(ins)
        A_contract_list[i] = -k
    end
    for (k,i) in enumerate(outs)
        A_contract_list[i] = k
    end
    AU = tensorcontract(A, A_contract_list, U, [1:I; -(n+1)],
                        [-1*(1:I); -(n+1)])
    UE = tensorcontract(U, [-1*(1:I); 1], Diagonal(E), [1; -(n+1)],
                        [-1*(1:I); -(n+1)])
    @test isapprox(AU, UE)

    # Hermitian
    I = rand(1:3)
    n = 2*I
    shp = rshape(n=n)
    perm = randperm(n)
    ins = perm[1:I]
    outs = perm[I+1:end]
    # Making sure that the left and right legs match.
    shp = collect(shp)
    for (o,i) in zip(outs, ins)
        shp[o] = shp[i]
    end
    shp = (shp...,)
    A = rtensor(shape=shp)
    # Make Hermitian
    hperm = collect(1:n)
    for (o,i) in zip(outs, ins)
        hperm[o], hperm[i] = hperm[i], hperm[o]
    end
    At = conj!(tensorcopy(A, collect(1:ndims(A)), hperm))
    A = (A + At)/2

    # Decompose and reconstruct without truncation
    E, U = tensoreig(A, ins, outs, hermitian=true)
    UE = tensorcontract(U, [-1*(1:I); 1], Diagonal(E), [1,-I-1])
    Areco = tensorcontract(UE, [-1*(1:I); 1], conj(U), [-1*(I+1:n); 1])
    Areco = tensorcopy(Areco, collect(1:ndims(Areco)), invperm(perm))

    # Decompose and reconstruct with truncation
    eps = 10.0^(-rand(0:3))
    dim_in = prod(shp[perm][1:I])
    dim_out = prod(shp[perm][I+1:end])
    max_dim = min(dim_in, dim_out)
    max_chi = rand(1:max_dim)
    chis = 1:max_chi
    E, U, error = tensoreig(A, ins, outs, eps=eps, chis=chis,
                            return_error=true, hermitian=true)
    # This test could fail if there were degenerate singular values. We assume
    # this is not the case because the tensors are random.
    @test (error<eps || length(E) == max_chi)
    UE = tensorcontract(U, [-1*(1:I); 1], Diagonal(E), [1,-I-1])
    Areco = tensorcontract(UE, [-1*(1:I); 1], conj(U), [-1*(I+1:n); 1])
    Areco = tensorcopy(Areco, collect(1:ndims(Areco)), invperm(perm))
    Anorm = norm(A)
    true_error = Anorm > 0 ? norm(Areco - A)/Anorm : norm(Areco)
    @test isapprox(error, true_error; atol=1e-8)
end

function test_split()
    # Non-hermitian, SVD
    shp = rshape(nlow=2)
    n = length(shp)
    perm = randperm(n)
    I = rand(1:n-1)
    J = n - I
    ins = perm[1:I]
    outs = perm[I+1:end]
    A = rtensor(shape=shp)

    B1, B2 = tensorsplit(A, ins, outs)
    U, S, Vt = tensorsvd(A, ins, outs)
    Ssqrt = Diagonal(sqrt.(S))
    B1reco = tensorcontract(U, [-1*(1:I); 1], Ssqrt, [1, -n-1])
    B2reco = tensorcontract(Ssqrt, [-1,1], Vt, [1; -1*(2:J+1)])
    @test isapprox(B1reco, B1)
    @test isapprox(B2reco, B2)

    # Hermitian, eigendecomposition
    I = rand(1:3)
    n = 2*I
    shp = rshape(n=n)
    perm = randperm(n)
    ins = perm[1:I]
    outs = perm[I+1:end]
    # Making sure that the left and right legs match.
    shp = collect(shp)
    for (o,i) in zip(outs, ins)
        shp[o] = shp[i]
    end
    shp = (shp...,)
    A = rtensor(shape=shp)
    # Make Hermitian
    hperm = collect(1:n)
    for (o,i) in zip(outs, ins)
        hperm[o], hperm[i] = hperm[i], hperm[o]
    end
    At = conj!(tensorcopy(A, collect(1:ndims(A)), hperm))
    A = (A + At)/2

    B1, B2 = tensorsplit(A, ins, outs, hermitian=true)
    E, U = tensoreig(A, ins, outs, hermitian=true)
    Esqrt = Diagonal(sqrt.(complex.(E)))
    B1reco = tensorcontract(U, [-1*(1:I); 1], Esqrt, [1, -n-1])
    B2reco = tensorcontract(Esqrt, [-1,1], conj(U), [-1*(2:I+1); 1])
    @test isapprox(B1reco, B1)
    @test isapprox(B2reco, B2)
end

function main()
    pars = parse_pars()
    if pars[:test_svd]
        println("Testing SVD.")
        for iter_num in 1:pars[:iters]
            test_svd()
        end
        println("Done testing SVD.")
    end
    if pars[:test_eig]
        println("Testing eig.")
        for iter_num in 1:pars[:iters]
            test_eig()
        end
        println("Done testing eig.")
    end
    if pars[:test_split]
        println("Testing split.")
        for iter_num in 1:pars[:iters]
            test_split()
        end
        println("Done testing split.")
    end
end

main()


text:
module MPSToolkit

using Logging
using LinearAlgebra
import ITensors.HDF5
import ITensors.Strided
using Dates

using ITensors
export setmaxdim!, setmindim!, setcutoff!, setnoise!, Sweeps, productMPS, randomMPS

include("Observables.jl")
include("Models/Model.jl")
include("Models/SpinHalf.jl")
include("Models/TransverseIsing1D.jl")
include("DMRG.jl")
include("TEBD.jl")
include("Util.jl")

if Base.VERSION >= v"1.4.2"
    include("Precompile/precompile_MPSToolkit.jl")
    _precompile_()
end

end


text:
using JDFManager
using Test

@testset "JDFManager.jl" begin
    # Write your tests here.
end


text:
"""
# StationXML

Implement the FDSN StationXML specification to describe seismic (and other)
stations recording regularly-sampled information such as ground velocity.

Full specification can be found at:

    https://www.fdsn.org/xml/station/fdsn-station-1.1.xsd

## Reading FDSN StationXML data

Two unexported functions are available for use in creating `FDSNStationXML` objects:

- `StationXML.read(filename)`: Read from a file on disk.
- `StationXML.readstring(string)`: Parse from a `String`.

For instance (using an example StationXML file supplied with this module):

```julia
julia> using StationXML

julia> sxml = StationXML.read(joinpath(dirname(pathof(StationXML)), "..", "data", "JSA.xml"))
StationXML.FDSNStationXML
  source: String "IRIS-DMC"
  sender: String "IRIS-DMC"
  module_name: String "IRIS WEB SERVICE: fdsnws-station | version: 1.1.36"
  module_uri: String "http://service.iris.edu/fdsnws/station/1/query?network=GB&station=JSA&level=response&format=xml&nodata=204"
  created: Dates.DateTime
  network: Array{StationXML.Network}((1,))
  schema_version: String "1.0"

```

## Accessing fields

You should access the fields of `FDSNStationXML` objects directly.  These match the
StationXML specification directly, and can also be listed for each of the types with
the usual `fieldnames(::DataType)` function.  (E.g., `fieldnames(Channel)`.)

To find out how many stations are in each of the networks returned in your request
XML, and what the network code is, you can do:

```julia
julia> for net in sxml.network
           println(net.code, ": ", net.total_number_stations, " stations")
       end
29
```

To get a vector of the station codes in the one network (GB) returned in our request:

```julia
julia> gb = sxml.network[1];

julia> stas = [sta.code for sta in gb.station]
1-element Array{String,1}:
 "JSA"

```

You can also access individual networks, stations and channels using functions with
these names.  For example, using the [SeisRequests](https://github.com/anowacki/SeisRequests.jl)
package to get all the broadband, high-gain channels stations in the GB network:

```julia
julia> using SeisRequests, Dates

julia> xml = get_request(
                 FDSNStation(starttime=DateTime(2012), endtime=now(), network="GB",
                             station="*", location="--", channel="BH?",
                             level="channel")).body |> String;

julia> [s.code for s in stations(sxml)]
28-element Array{String,1}:
 "BIGH"
 "CCA1"
 "CLGH"
 "CWF"
 "DRUM"
 "DYA"
 "EDI"
 "EDMD"
 "ELSH"
 "ESK"
 "FOEL"
 "GAL1"
 "HMNX"
 "HPK"
 "HTL"
 "IOMK"
 "JSA"
 "KESW"
 "KPL"
 "LBWR"
 "LMK"
 "LRW"
 "MCH1"
 "SOFL"
 "STNC"
 "SWN1"
 "WACR"
 "WLF1"

```

## Accessor functions

You can easily construct vectors of all the networks, stations and channels in the StationXML
using the following accessor functions:

- networks(stationxml)
- stations(stationxml_or_network)
- channels(stationxml_or_network_or_station)

Note that `station`, for instance, accepts either a `Network` or a whole `FDSNStationXML`
object, whilst either of those or a `Station` can be given to `channels`.

```julia
julia> stations(gb)
28-element Array{StationXML.Station,1}:
 StationXML.Station(missing, StationXML.Comment[], "BIGH", 2009-12-15T00:00:00, 2599-12-31T23:59:59, StationXML.RestrictedStatus
  value: String "open"
...
```

The `channel_codes` function returns a list of all of the channel codes within
a `FDSNStationXML` document or a `Network`.

## Dot-access to arrays of objects

The module defines `getproperty` methods for conveniently accessing the fields of each member
of arrays of `Network`s, `Station`s and `Channel`s.  So our previous example of finding all
the station codes could actually have been done like this:

```julia
julia> stations(sxml).code
28-element Array{String,1}:
 "BIGH"
 "CCA1"
 "CLGH"
 "CWF"
...
```

We can equally access any other field of the items this way:

```julia
julia> channels(sxml).longitude
84-element Array{Float64,1}:
 -3.9087
 -3.9087
 -3.9087
 -5.227299
 -5.227299
 -5.227299
 -6.110599
 -6.110599
 -6.110599
...
```

## Preamble to the specification

#### FDSN StationXML (www.fdsn.org/xml/station)

The purpose of this schema is to define an XML representation of the most important
and commonly used structures of SEED 2.4 metadata.

The goal is to allow mapping between SEED 2.4 dataless SEED volumes and this schema with as
little transformation or loss of information as possible while at the same time simplifying
station metadata representation when possible.  Also, content and clarification has been added
where lacking in the SEED standard.

When definitions and usage are underdefined the SEED manual should be referred to for
clarification.  SEED specifiation: http://www.fdsn.org/publications.htm

Another goal is to create a base schema that can be extended to represent similar data types.

#### Versioning for FDSN StationXML:

The `version` attribute of the schema definition identifies the version of the schema.  This
version is not enforced when validating documents.

The required `schemaVersion` attribute of the root element identifies the version of the schema
that the document is compatible with.  Validation only requires that a value is present but
not that it matches the schema used for validation.

The `targetNamespace` of the document identifies the major version of the schema and document,
version 1.x of the schema uses a target namespace of "http://www.fdsn.org/xml/station/1".
All minor versions of a will be backwards compatible with previous minor releases.  For
example, all 1.x schemas are backwards compatible with and will validate documents for 1.0.
Major changes to the schema that would break backwards compabibility will increment the major
version number, e.g. 2.0, and the namespace, e.g. "http://www.fdsn.org/xml/station/2".

This combination of attributes and `targetNamespace`s allows the schema and documents to be
versioned and allows the schema to be updated with backward compatible changes (e.g. 1.2)
and still validate documents created for previous major versions of the schema (e.g. 1.0).

"""
module StationXML

using Dates

using Mixers: @pour
using Parameters: @with_kw
import EzXML
import DocStringExtensions

export
  FDSNStationXML,
  channel_codes,
  channels,
  networks,
  stations,
  xmldoc

const M{T} = Union{T,Missing}

"The schema version used for writing files"
const DEFAULT_SCHEMA_VERSION = "1.1"

include("compat.jl")
include("util.jl")
include("base_types.jl")
include("derived_types.jl")
include("types.jl")
include("deprecations.jl")
include("accessors.jl")
include("io.jl")
include("merge.jl")

include("precompile.jl")

end # module


text:
## =======================================================================================
## Remove before v0.11
## =======================================================================================

# see JuliaRobotics/Caesar.jl#811
@deprecate plotKDE(w...;kw...) plotBelief(w...;kw...)

## =======================================================================================
## Remove before v0.10
## =======================================================================================

# function getColorsByLength(len::Int=7)
#   len > 99 ? error("Don't have enough colors, 100 is the max.") : nothing
#   COLORS = String["red";"green";"blue";"magenta";"yellow";"deepskyblue"]
#   if len > 6
#     scale = len -7 + 2 # + 2 is artificial and avoids gray100==white
#     scale = 100/scale
#     for i in 7:len
#       push!(COLORS, "gray$(floor(Int,(i-7)*scale))")
#     end
#   end
#   return COLORS[1:len]
# end

export togglePrtStbLines
export stbPrtLineLayers!
# export drawHorDens # from KDEPlotting
# export plotMCMC, plotPose2DMC!, mcmcPose2D!

# for some reason we still need msgPlots of right size in the global for these functions to work.
# precall drawTreeUpwardMsgs or drawFrontalDens to make this work properly TODO
function vstackedDensities(msgPlots)
  #msgPlots = f(fg, bt) # drawTreeUpwardMsgs
  evalstr = ""
  for i in 1:length(msgPlots)
      evalstr = string(evalstr, ",msgPlots[$(i)]")
  end
  eval(parse(string("vstack(",evalstr[2:end],")")))
end

global DISABLESTBPRTLINES = false

function togglePrtStbLines()
  global DISABLESTBPRTLINES
  DISABLESTBPRTLINES = !DISABLESTBPRTLINES
end

## TODO -- you were here with port starboard lines
function stbPrtLineLayers!(pl, Xpp, Ypp, Thpp; l::Real=5.0)
  if DISABLESTBPRTLINES
    return nothing
  end

  lnstpr = [0.0;l;0.0]
  lnstpg = [0.0;-l;0.0]

  Rd  =SE2(lnstpr)
  Gr = SE2(lnstpg)

  for i in 1:length(Xpp)
    lnstt = [Xpp[i];Ypp[i];Thpp[i]]
    Ps = SE2(lnstt)
    lnr = se2vee(Ps*Rd)
    lng = se2vee(Ps*Gr)
    xsr = [Xpp[i];lnr[1]]
    ysr = [Ypp[i];lnr[2]]
    xsg = [Xpp[i];lng[1]]
    ysg = [Ypp[i];lng[2]]

    push!(pl.layers, layer(x=xsr, y=ysr, Geom.path(), Gadfly.Theme(default_color=colorant"red", line_width=1.5pt))[1] )
    push!(pl.layers, layer(x=xsg, y=ysg, Geom.path(), Gadfly.Theme(default_color=colorant"green", line_width=1.5pt))[1] )
  end
  nothing
end

#


text:
include("../src/NativeInterface/options_structure.jl")

@testset "Triangle.NativeInterface.TriangulateOptions" begin
    @testset "getTriangulateStringOptions()" begin
        @testset "Default" begin
            opt = TriangulateOptions()
            @test getTriangulateStringOptions(opt) == "BPIQ"
        end

        @testset "All disable" begin
            opt = TriangulateOptions()
            opt.nobound = false
            opt.nopolywritten = false
            opt.noiterationnum = false
            opt.quiet = false
            @test getTriangulateStringOptions(opt) == ""
        end

        @testset "All enable" begin
            opt = TriangulateOptions()
            opt.pslg = true
            opt.regionattrib = true
            opt.convex = true
            opt.jettison = true
            opt.firstnumberiszero = true
            opt.edgesout = true
            opt.voronoi = true
            opt.neighbors = true
            opt.nobound = true
            opt.nopolywritten = true
            opt.nonodewritten = true
            opt.noelewritten = true
            opt.noiterationnum = true
            opt.noholes = true
            opt.noexactaritmetic = true
            opt.order = true
            opt.orderHow = 1
            opt.dwyer = true
            opt.quiet = true
            opt.verbose = true
            @test getTriangulateStringOptions(opt) == "pAcjzevnBPNEIOXo1lQV"
        end
    end
end


text:
# This file is a part of PicoScopes.jl, licensed under the MIT License (MIT).

import Compat.Test
Test.@testset "Package PicoScopes" begin
    # write your own tests here
end

max_stars_repo_name: mppmu/PicoScopes

max_issues_repo_name: mppmu/PicoScopes

max_forks_repo_name: mppmu/PicoScopes


text:
cd(@__DIR__)
using Pkg
# you may need to use Pkg.instantiate() the first time you run the file
Pkg.activate("")
using Plots
include("functions.jl")

# current chain
x = [.8,.5]
# random chain
z = [1.5,2.5]
# difference of current and random z
# forms dotted line
d = x - z

# random chains r1 and r2
r1 = [0,-1.5]
r2 = [-.3,-2.1]

# project r1 onto line spanned by d = x - z
proj1 = project(r1, d)
# project r2 onto line spanned by d = x - z
proj2 = project(r2, d)

# difference scalar set to 1 for simplicity
γ = 1.0

# difference in projections
dproj1_2 = proj1 - proj2

# proposal x + γ * (proj1 - proj2)
proposal = x + γ * dproj1_2

scatter(
    [x[1]],
    [x[2]],
    xlims = (-3.0, 3.0),
    ylims = (-3.0, 3.0),
    color = :black,
    linewidth = 2,
    label = "",
    framestyle = :origin,
)
annotate!(x[1]-.1,x[2]+.2, text("X", :black, 12))

scatter!(
    [z[1]],
    [z[2]],
    xlims = (-3.0, 3.0),
    ylims = (-3.0, 3.0),
    color = :black,
    linewidth = 2,
    label = "",
    framestyle = :origin,
)
annotate!(z[1]-.1,z[2] +.1, text("Z", :black, 12))

xvals = -30:.1:30
β = get_coefs(x, z)
yvals = make_line(β..., xvals)
plot!(
    xvals,
    yvals,
    color = :black,
    linestyle = :dash,
    linewidth = 2,
    label = "",
    framestyle = :origin,
)

scatter!(
    [r1[1]],
    [r1[2]],
    color = :black,
    linewidth = 2,
    label = "",
    framestyle = :origin,
)
annotate!(r1[1]-.1,r1[2] +.2, text("R1", :black, 12))

scatter!(
    [r2[1]],
    [r2[2]],
    color = :black,
    linewidth = 2,
    label = "",
    framestyle = :origin,
)
annotate!(r2[1]-.1,r2[2] +.2, text("R2", :black, 12))

scatter!(
    [proposal[1]],
    [proposal[2]],
    color = :red,
    linewidth = 2,
    label = "",
    framestyle = :origin,
)
annotate!(proposal[1] - .1,proposal[2] +.1, text("X*", :black, 12))


text:
"""
    DomainDecompositions

Provides a domain decomposition into
 - `grid_mean`
 - `environment`
 - `updrafts`
"""
module DomainDecompositions

export DomainDecomposition
export gridmean, environment, updraft
export over_sub_domains, allcombinations
export var_names, var_string, var_suffix

include("types.jl")
include("subsets.jl")
include("indexes.jl")
include("variable_mapper.jl")

"""
    DomainDecomposition(vars, domain_set::DomainSet) where {T}

Return a domain decomposition, given a tuple of tuples of variable
names and the number of their subdomains.
"""
function DomainDecomposition(vars, domain_set::DomainSet) where {T}
    n_subdomains = sum(domain_set)
    n_vars = sum([sum(domain_set, dss) for (v, dss) in vars])
    var_mapper, dss_per_var, var_names = get_var_mapper(vars, domain_set)
    idx = DomainIdx(domain_set)
    sd_unmapped = get_sd_unmapped(vars, idx, domain_set)

    idx_ss_per_var =
        Dict(ϕ => DomainIdx(domain_set, dss_per_var[ϕ]) for ϕ in var_names)
    a_map = Dict(ϕ => get_sv_a_map(idx, idx_ss_per_var[ϕ]) for ϕ in var_names)

    return DomainDecomposition{typeof(var_mapper)}(
        idx,
        var_mapper,
        idx_ss_per_var,
        a_map,
        sd_unmapped,
        domain_set,
        dss_per_var,
    )
end

gridmean(domain_decomp::DomainDecomposition) = gridmean(domain_decomp.idx)
environment(domain_decomp::DomainDecomposition) = environment(domain_decomp.idx)
updraft(domain_decomp::DomainDecomposition) = updraft(domain_decomp.idx)
alldomains(domain_decomp::DomainDecomposition) = alldomains(domain_decomp.idx)

"""
    over_sub_domains(domain_decomp::DomainDecomposition, ϕ::Symbol)

Get list of indexes over all subdomains for variable `ϕ`.
"""
function over_sub_domains(domain_decomp::DomainDecomposition, ϕ::Symbol)
    return [x for x in domain_decomp.sd_unmapped[ϕ] if !(x == 0)]
end

"""
    var_names(domain_decomp::DomainDecomposition)

Get tuple of variable names.
"""
var_names(domain_decomp::DomainDecomposition{VM}) where {VM} = fieldnames(VM)

function Base.show(io::IO, domain_decomp::DomainDecomposition)
    println(io, "\n----------------------- Domain Decomposition")
    println(io, "DomainIdx")
    println(io, domain_decomp.idx)
    println(io, "var_mapper")
    println(io, domain_decomp.var_mapper)
    println(io, "\n-----------------------")
end

DomainIdx(domain_decomp::DomainDecomposition) = domain_decomp.idx

allcombinations(domain_decomp::DomainDecomposition) =
    allcombinations(DomainIdx(domain_decomp))

end


text:
#Playing around with the mean polynomial function

using Gadfly, GaussianProcesses

#Training data
x=[-4.0,-3.0,-1.0,0.0,2.0];

# Create quadratic function
y = 0.5*x.^2 + 2.0*x + 0.5*rand(5)

#Test data
xpred = collect(-5.0:0.1:5.0);

#Specify mean and covariance function
beta = [0.5 2.0]
mPoly = MeanPoly(beta)
kern = SE(0.0,0.0)
gp = GP(x,y,mPoly,kern)  #Fit the GP

plot(gp)   #Plot the GP


text:
module DipTransform
using SparseArrays
using LinearAlgebra
using jInv.Mesh
using jInv.Utils
using jInv.InverseSolve
using ShapeReconstructionPaLS.Utils
using ParamLevelSet
using MAT
using Distributed

import jInv.ForwardShare.getData
import jInv.ForwardShare.getSensTMatVec
import jInv.ForwardShare.getSensMatVec
import jInv.ForwardShare.ForwardProbType

## The rotation planes are:
# first rotation (theta) is for planes x1 and x2
# second rotation (phi) is for planes x1 and x3
# Sampling is done by summing planes of x1 and x2 for all values of x3.

export DipParam, getDipParam
mutable struct DipParam <: ForwardProbType
    Mesh      				:: RegularMesh
	ndipsAll				:: Int64
	workerSubIdxs
    theta_phi_rad			:: Array{Float64,2}
	b						:: Array{Float64,2}
	samplingBinning			:: Int64
	method					:: String
	Jacobian				:: SparseMatrixCSC{Float32,Int32}
	S						:: SparseMatrixCSC{Float32,Int32}
end

function getDipParamInternal(Mesh::RegularMesh,ndipsAll::Int64,workerSubIdxs,theta_phi_rad::Array{Float64,2},b,samplingBinning::Int64,method)
	return DipParam(Mesh, ndipsAll, workerSubIdxs, theta_phi_rad,b, samplingBinning, method, spzeros(Float32,Int32,0,0), spzeros(Float32,Int32,0,0));
end

function getDipParam(Mesh::RegularMesh,theta_phi_rad::Array{Float64,2},b,samplingBinning::Int64,numWorkers::Int64,method = MATBased)
	## This function does use the parallel mechanism of jInv (i.e., returns a RemoteChannel), even if numWorkers=1.
	if numWorkers > nworkers()
		numWorkers = nworkers();
	end
	SourcesSubInd = Array{Array{Int64,1}}(undef,numWorkers);
	ActualWorkers = workers();
	if numWorkers < nworkers()
		ActualWorkers = ActualWorkers[1:numWorkers];
	end
	pFor   = Array{RemoteChannel}(undef,numWorkers)
	i = 1; nextidx() = (idx=i; i+=1; idx)

	ndips  = size(theta_phi_rad,1);
	# send out jobs
	@sync begin
		for p = ActualWorkers
			@async begin
				while true
					idx = nextidx()
					if idx > numWorkers
						break
					end
					I_p = getIndicesOfKthWorker(numWorkers,idx,ndips);
					# println("Sending ",collect(I_p)," to worker ", p);
					# find src and rec on mesh
					pFor[idx]  = initRemoteChannel(getDipParamInternal,p, Mesh,ndips,I_p,theta_phi_rad[I_p,:],b[I_p,:],samplingBinning,method);
					wait(pFor[idx]);
				end
			end
		end
	end
	return pFor # Array of Remote Refs
end

import jInv.Utils.clear!
function clear!(pFor::DipParam)
	pFor.SampleMat = spzeros(0);
	clear!(pFor.Mesh);
	return pFor;
end

include("generateDipMatrix.jl");
include("prepareSyntheticDipData.jl");
include("getData.jl");
include("sensitivityFuncs.jl");
end

max_stars_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_issues_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_forks_repo_name: BGUCompSci/ShapeReconstructionPaLS


text:
const ASSETS_SUMMARY_PERIOD = 1994:2015
const RETURNS_SUMMARY_PERIOD = 2000:2015
const BY_CHARITY_FIELD = :charitytype
const BY_CATEGORY_FIELD = :categoryname

function assetsByYearDF(data::NCCSData; summaryPeriod::UnitRange = ASSETS_SUMMARY_PERIOD,
  focalField::Symbol = :adjnetassets)::DataFrame

  local outdf::DataFrame = DataFrame(fisyr = collect(summaryPeriod)) #this will hold the summary stats
  local subdf::SubDataFrame =  view(data.df,  (y::Int-> y ∈ summaryPeriod).(data.df[:fisyr]),
    [:ein, :fisyr, :charitytype, focalField])
  local nyears::Int = length(summaryPeriod) # of years

  local N::Vector{MInt} = Vector{MInt}(undef, nyears)
  local total::Vector{MFloat64} = Vector{MFloat64}(undef, nyears)
  local mean::Vector{MFloat64} = Vector{MFloat64}(undef, nyears)
  local stddev::Vector{MFloat64} = Vector{MFloat64}(undef, nyears)
  local skew::Vector{MFloat64} = Vector{MFloat64}(undef, nyears)
  local kurt::Vector{MFloat64} = Vector{MFloat64}(undef, nyears)
  local p5::Vector{MFloat64} = Vector{MFloat64}(undef, nyears)
  local p25::Vector{MFloat64} = Vector{MFloat64}(undef, nyears)
  local p50::Vector{MFloat64} = Vector{MFloat64}(undef, nyears)
  local p75::Vector{MFloat64} = Vector{MFloat64}(undef, nyears)
  local p95::Vector{MFloat64} = Vector{MFloat64}(undef, nyears)
  local fracFoundation::Vector{MFloat64} = Vector{MFloat64}(undef, nyears)
  local wealth80::Vector{MFloat64} = Vector{MFloat64}(undef, nyears)
  local wealth95::Vector{MFloat64} = Vector{MFloat64}(undef, nyears)

  #construct the info for each year

  #make an index so we can use the high performance by function
  local outdfIndex::Dict = Dict(outdf[i, :fisyr]=>i for i::Int ∈ 1:nyears)

  #this format is much faster since it does not imply recombination
  subdfByfisyr::GroupedDataFrame = groupby(subdf, :fisyr)
  for ssubdf::SubDataFrame in subdfByfisyr
  #by(subdf, :fisyr) do ssubdf
    r::Int = outdfIndex[ssubdf[1, :fisyr]]
    focal::Vector = ssubdf[:, focalField]
    nfirms::Int = length(focal)
    #println("Vector Length nfirms: $nfirms")
    focal = processmissing(focal) #handle missing values

    #println("Vector Length: $(length(focal))")
    N[r] = length(focal)
    total[r] = sum(focal)
    mean[r] = StatsBase.mean(focal)
    stddev[r] = StatsBase.std(focal)
    skew[r] = StatsBase.skewness(focal)
    kurt[r] = StatsBase.kurtosis(focal)
    p5[r] = StatsBase.quantile(focal, .05)
    p25[r] = StatsBase.quantile(focal, .25)
    p50[r] = StatsBase.quantile(focal, 0.5)
    p75[r] = StatsBase.quantile(focal, 0.75)
    p95[r] = StatsBase.quantile(focal, 0.95)
    fracFoundation[r] = sum(ssubdf[:charitytype] .== :pf) / N[r]

    #now get a couple distributional values
    cumv::Vector{Float64} = similar(focal) #will hold the running total
    cumv[1] = focal[1]
    for j::Int ∈ 2:length(focal)
      cumv[j] = cumv[j-1] + focal[j]
    end
    cumv .= cumv ./ cumv[end] #normalize

    wealth80[r] = sum(cumv .≤ 0.8) / N[r]
    wealth95[r] = sum(cumv .≤ 0.95) / N[r]
  end

  #now put the data into the dataframe
  outdf[:N] = N
  outdf[:total] = total
  outdf[:mean] = mean
  outdf[:stddev] = stddev
  outdf[:skew] = skew
  outdf[:kurt] = kurt
  outdf[:p5] = p5
  outdf[:p25] = p25
  outdf[:p50] = p50
  outdf[:p75] = p75
  outdf[:p95] = p95
  outdf[:fracFoundation] = fracFoundation
  outdf[:wealth80] = wealth80
  outdf[:wealth95] = wealth95

  return outdf
end

function assetsByYear(outdf::DataFrame; summaryPeriod::UnitRange = ASSETS_SUMMARY_PERIOD,
  focalField::Symbol = :adjnetassets, decimals::Int = DECIMALS)::String

  columnSymbols::Vector{Symbol} = setdiff(names(outdf), [:fisyr])

  numColumns = length(columnSymbols) #number of columns excluding the year
  columnIndex::Dict = Dict(columnSymbols[i]=>i for i ∈ 1:numColumns)

  #include all columns in the names
  columnNames::Vector{Vector{String}} = [(string).(columnSymbols)]
  descRowNames::Vector{String} = (string).(outdf[:fisyr])

  numRows = size(outdf,1)

  #need to print the descriptive rows
  descContent::Vector{Vector{String}} = ((i::Int)->Vector{String}(undef, numColumns)).(1:length(descRowNames))

  for r::Int ∈ 1:numRows
    haskey(columnIndex, :N) && (
      descContent[r][columnIndex[:N]] = "$(outdf[r, :N])")
    haskey(columnIndex, :total) && (
      descContent[r][columnIndex[:total]] = num2Str(outdf[r, :total], 0, scaleFactor=10^-6))
    haskey(columnIndex, :mean) && (
      descContent[r][columnIndex[:mean]] = num2Str(outdf[r, :mean], decimals, scaleFactor=10^-6))
    haskey(columnIndex, :stddev) && (
      descContent[r][columnIndex[:stddev]] = num2Str(outdf[r, :stddev], decimals, scaleFactor=10^-6))
    haskey(columnIndex, :skew) && (
      descContent[r][columnIndex[:skew]] = num2Str(outdf[r, :skew], decimals, scaleFactor=1.0))
    haskey(columnIndex, :kurt) && (
      descContent[r][columnIndex[:kurt]] = num2Str(outdf[r, :kurt], decimals, scaleFactor=1.0))
    haskey(columnIndex, :p5) && (
      descContent[r][columnIndex[:p5]] = num2Str(outdf[r, :p5], decimals, scaleFactor=10^-6))
    haskey(columnIndex, :p25) && (
      descContent[r][columnIndex[:p25]] = num2Str(outdf[r, :p25], decimals, scaleFactor=10^-6))
    haskey(columnIndex, :p50) && (
      descContent[r][columnIndex[:p50]] = num2Str(outdf[r, :p50], decimals, scaleFactor=10^-6))
    haskey(columnIndex, :p75) && (
      descContent[r][columnIndex[:p75]] = num2Str(outdf[r, :p75], decimals, scaleFactor=10^-6))
    haskey(columnIndex, :p95) && (
      descContent[r][columnIndex[:p95]] = num2Str(outdf[r, :p95], decimals, scaleFactor=10^-6))
    haskey(columnIndex, :fracFoundation) && (
      descContent[r][columnIndex[:fracFoundation]] = num2Str(outdf[r, :fracFoundation], 2))
    haskey(columnIndex, :wealth80) && (
      descContent[r][columnIndex[:wealth80]] = num2Str(outdf[r, :wealth80], 2))
    haskey(columnIndex, :wealth95) && (
      descContent[r][columnIndex[:wealth95]] = num2Str(outdf[r, :wealth95], 2))
  end

  assetsByYearTable::String = texTable(titleCaption="Descriptive Data by Year",
      caption="""See Tex File""", #caption
      colNames= columnNames, #colNames
      descRowNames=descRowNames, #descRowNames
      descContent = descContent, #descContent
      nakedTable=true
    )

    return assetsByYearTable

end

function assetsSummaryDF(data::NCCSData; summaryPeriod::UnitRange = ASSETS_SUMMARY_PERIOD,
  focalField::Symbol = :adjnetassets, decimals::Int = DECIMALS,
  byCharityField::Symbol = BY_CHARITY_FIELD,
  byCategoryField::Symbol = BY_CATEGORY_FIELD)::DataFrame

  local subdf::SubDataFrame =  view(data.df,  (y::Int-> y ∈ summaryPeriod).(data.df[:fisyr]),
    [:ein; :fisyr; :fullsample; byCharityField; byCategoryField; focalField])
  local byCharityValues::Vector{Symbol} = unique(processmissing(subdf[byCharityField]))
  local byCategoryValues::Vector{Symbol} = unique(processmissing(subdf[byCategoryField]))
  local rowFields::Vector{Symbol} = [byCharityValues; byCategoryValues; :fullsample]
  local nrows::Int = length(rowFields) # of rows we will have

  local N::Vector{MInt} = Vector{MInt}(undef, nrows)
  local minyear::Vector{MInt} = Vector{MInt}(undef, nrows)
  local maxyear::Vector{MInt} = Vector{MInt}(undef, nrows)
  local total::Vector{MFloat64} = Vector{MFloat64}(undef, nrows)
  local mean::Vector{MFloat64} = Vector{MFloat64}(undef, nrows)
  local stddev::Vector{MFloat64} = Vector{MFloat64}(undef, nrows)
  local skew::Vector{MFloat64} = Vector{MFloat64}(undef, nrows)
  local kurt::Vector{MFloat64} = Vector{MFloat64}(undef, nrows)
  local p5::Vector{MFloat64} = Vector{MFloat64}(undef, nrows)
  local p25::Vector{MFloat64} = Vector{MFloat64}(undef, nrows)
  local p50::Vector{MFloat64} = Vector{MFloat64}(undef, nrows)
  local p75::Vector{MFloat64} = Vector{MFloat64}(undef, nrows)
  local p95::Vector{MFloat64} = Vector{MFloat64}(undef, nrows)

  #performance optimized way to split up the dataframe
  local outdfIndex::Dict = Dict(rowFields[i]=>i for i::Int ∈ 1:nrows)
  local rowNames::Vector{Symbol} =
    (s::Symbol->haskey(charitytype, s) ? charitytype[s] : s).(rowFields) #clarify some of the names
  local outdf::DataFrame = DataFrame(rownames = rowNames) #this will hold the summary stats

  #for each field we are considering
  #println("keys: $(collect(keys(outdfIndex)))")
  for f::Symbol ∈ [byCharityField; byCategoryField; :fullsample]
    #println("field $f")
    subdfByf::GroupedDataFrame = groupby(subdf, f)
    for ssubdf::SubDataFrame ∈ subdfByf
    #by(subdf, f) do ssubdf::SubDataFrame
      if !ismissing(ssubdf[1, f])
        r::Int = outdfIndex[ssubdf[1, f]]

        focal::Vector = ssubdf[:, focalField]
        nfirms::Int = length(focal)
        focal = processmissing(focal) #handle missing values

        N[r] = length(focal)
        minyear[r] = minimum(processmissing(ssubdf[:, :fisyr]))
        maxyear[r] = maximum(processmissing(ssubdf[:, :fisyr]))
        total[r] = StatsBase.sum(focal)
        mean[r] = StatsBase.mean(focal)
        stddev[r] = StatsBase.std(focal)
        skew[r] = StatsBase.skewness(focal)
        kurt[r] = StatsBase.kurtosis(focal)
        p5[r] = StatsBase.quantile(focal, .05)
        p25[r] = StatsBase.quantile(focal, .25)
        p50[r] = StatsBase.quantile(focal, 0.5)
        p75[r] = StatsBase.quantile(focal, 0.75)
        p95[r] = StatsBase.quantile(focal, 0.95)
      end
    end
  end

  #now put the data into the dataframe
  outdf[:N] = N
  outdf[:minyear] = minyear
  outdf[:maxyear] = maxyear
  outdf[:total] = total
  outdf[:mean] = mean
  outdf[:stddev] = stddev
  outdf[:skew] = skew
  outdf[:kurt] = kurt
  outdf[:p5] = p5
  outdf[:p25] = p25
  outdf[:p50] = p50
  outdf[:p75] = p75
  outdf[:p95] = p95

  return outdf
end

function assetsSummary(outdf::DataFrame; summaryPeriod::UnitRange = ASSETS_SUMMARY_PERIOD,
  focalField::Symbol = :adjnetassets, decimals::Int = DECIMALS,
  byCharityField::Symbol = BY_CHARITY_FIELD,
  byCategoryField::Symbol = BY_CATEGORY_FIELD)::String

  columnSymbols::Vector{Symbol} = setdiff(names(outdf), [:rownames])

  numColumns = length(columnSymbols) #number of columns excluding the year
  columnIndex::Dict = Dict(columnSymbols[i]=>i for i ∈ 1:numColumns)

  #include all columns in the names
  columnNames::Vector{Vector{String}} = [(string).(columnSymbols)]
  descRowNames::Vector{String} = (string).(outdf[:rownames])

  #NOTE: a couple of hacks to improve readability
  descRowNames[length(keys(charitytype))+1] = "\\midrule $(descRowNames[length(keys(charitytype))+1])"
  descRowNames[end] = "\\midrule $(descRowNames[end])"

  numRows = size(outdf,1)

  #need to print the descriptive rows
  descContent::Vector{Vector{String}} = (
    (i::Int)->Vector{String}(undef, numColumns)).(1:length(descRowNames))

  for r::Int ∈ 1:numRows
    haskey(columnIndex, :N) && (
      descContent[r][columnIndex[:N]] = "$(outdf[r, :N])")
    haskey(columnIndex, :minyear) && (
      descContent[r][columnIndex[:minyear]] = "$(outdf[r, :minyear])")
    haskey(columnIndex, :maxyear) && (
      descContent[r][columnIndex[:maxyear]] = "$(outdf[r, :maxyear])")
    haskey(columnIndex, :total) && (
      descContent[r][columnIndex[:total]] = num2Str(outdf[r, :total], 0, scaleFactor=10^-6))
    haskey(columnIndex, :mean) && (
      descContent[r][columnIndex[:mean]] = num2Str(outdf[r, :mean], decimals, scaleFactor=10^-6))
    haskey(columnIndex, :stddev) && (
      descContent[r][columnIndex[:stddev]] = num2Str(outdf[r, :stddev], decimals, scaleFactor=10^-6))
    haskey(columnIndex, :skew) && (
      descContent[r][columnIndex[:skew]] = num2Str(outdf[r, :skew], decimals, scaleFactor=1.0))
    haskey(columnIndex, :kurt) && (
      descContent[r][columnIndex[:kurt]] = num2Str(outdf[r, :kurt], decimals, scaleFactor=1.0))
    haskey(columnIndex, :p5) && (
      descContent[r][columnIndex[:p5]] = num2Str(outdf[r, :p5], decimals, scaleFactor=10^-6))
    haskey(columnIndex, :p25) && (
      descContent[r][columnIndex[:p25]] = num2Str(outdf[r, :p25], decimals, scaleFactor=10^-6))
    haskey(columnIndex, :p50) && (
      descContent[r][columnIndex[:p50]] = num2Str(outdf[r, :p50], decimals, scaleFactor=10^-6))
    haskey(columnIndex, :p75) && (
      descContent[r][columnIndex[:p75]] = num2Str(outdf[r, :p75], decimals, scaleFactor=10^-6))
    haskey(columnIndex, :p95) && (
      descContent[r][columnIndex[:p95]] = num2Str(outdf[r, :p95], decimals, scaleFactor=10^-6))
  end

  assetsByYearTable::String = texTable(titleCaption="Descriptive Data by Charity Type and Category",
      caption="""See Tex File""", #caption
      colNames= columnNames, #colNames
      descRowNames=descRowNames, #descRowNames
      descContent = descContent, #descContent
      nakedTable=true
    )

    return assetsByYearTable
end

function returnsByYearDF(data::NCCSData; summaryPeriod::UnitRange = RETURNS_SUMMARY_PERIOD,
  focalField::Symbol = :return, assetsField::Symbol = :adjnetassets)::DataFrame

  local outdf::DataFrame = DataFrame(fisyr = collect(summaryPeriod)) #this will hold the summary stats
  local subdf::SubDataFrame =  view(data.df,  (y::Int-> y ∈ summaryPeriod).(data.df[:fisyr]),
    [:ein, :fisyr, :charitytype, focalField, assetsField])
  local nyears::Int = length(summaryPeriod) # of years

  local N::Vector{MInt} = Vector{MInt}(undef, nyears)
  local equalmean::Vector{MFloat64} = Vector{MFloat64}(undef, nyears)
  local valuemean::Vector{MFloat64} = Vector{MFloat64}(undef, nyears)
  local stddev::Vector{MFloat64} = Vector{MFloat64}(undef, nyears)
  local skew::Vector{MFloat64} = Vector{MFloat64}(undef, nyears)
  local kurt::Vector{MFloat64} = Vector{MFloat64}(undef, nyears)
  local p5::Vector{MFloat64} = Vector{MFloat64}(undef, nyears)
  local p25::Vector{MFloat64} = Vector{MFloat64}(undef, nyears)
  local p50::Vector{MFloat64} = Vector{MFloat64}(undef, nyears)
  local p75::Vector{MFloat64} = Vector{MFloat64}(undef, nyears)
  local p95::Vector{MFloat64} = Vector{MFloat64}(undef, nyears)

  #construct the info for each year
  local outdfIndex::Dict = Dict(outdf[i, :fisyr]=>i for i::Int ∈ 1:nyears)

  subdfByfisyr::GroupedDataFrame = groupby(subdf, :fisyr)
  for ssubdf::SubDataFrame ∈ subdfByfisyr
    r::Int = outdfIndex[ssubdf[1, :fisyr]]
    focal::Vector = ssubdf[:, focalField]
    nfirms::Int = length(focal)
    #println("Vector Length nfirms: $nfirms")
    focal = processmissing(focal) #handle missing values

    #println("Vector Length: $(length(focal))")
    N[r] = length(focal)
    equalmean[r] = StatsBase.mean(focal)
    valuemean[r] = #compute the value-weighted average
      sum(processmissing(ssubdf[:, focalField] .* ssubdf[:, assetsField])) /
      sum(processmissing(ssubdf[(!ismissing).(ssubdf[:, focalField]), assetsField]))
    stddev[r] = StatsBase.std(focal)
    skew[r] = StatsBase.skewness(focal)
    kurt[r] = StatsBase.kurtosis(focal)
    p5[r] = StatsBase.quantile(focal, .05)
    p25[r] = StatsBase.quantile(focal, .25)
    p50[r] = StatsBase.quantile(focal, 0.5)
    p75[r] = StatsBase.quantile(focal, 0.75)
    p95[r] = StatsBase.quantile(focal, 0.95)
  end

  #now put the data into the dataframe
  outdf[:N] = N
  outdf[:equalmean] = equalmean
  outdf[:valuemean] = valuemean
  outdf[:stddev] = stddev
  outdf[:skew] = skew
  outdf[:kurt] = kurt
  outdf[:p5] = p5
  outdf[:p25] = p25
  outdf[:p50] = p50
  outdf[:p75] = p75
  outdf[:p95] = p95

  return outdf
end

function returnsByYear(outdf::DataFrame; summaryPeriod::UnitRange = RETURNS_SUMMARY_PERIOD,
  focalField::Symbol = :return, decimals::Int = DECIMALS_RETURN)::String

  columnSymbols::Vector{Symbol} = setdiff(names(outdf), [:fisyr])

  numColumns = length(columnSymbols) #number of columns excluding the year
  columnIndex::Dict = Dict(columnSymbols[i]=>i for i ∈ 1:numColumns)

  #include all columns in the names
  columnNames::Vector{Vector{String}} = [(string).(columnSymbols)]
  descRowNames::Vector{String} = (string).(outdf[:fisyr])

  numRows = size(outdf,1)

  #need to print the descriptive rows
  descContent::Vector{Vector{String}} = ((i::Int)->Vector{String}(undef, numColumns)).(1:length(descRowNames))

  for r::Int ∈ 1:numRows
    haskey(columnIndex, :N) && (
      descContent[r][columnIndex[:N]] = "$(outdf[r, :N])")
    haskey(columnIndex, :equalmean) && (
      descContent[r][columnIndex[:equalmean]] = num2Str(outdf[r, :equalmean], decimals, scaleFactor=100.0))
    haskey(columnIndex, :valuemean) && (
      descContent[r][columnIndex[:valuemean]] = num2Str(outdf[r, :valuemean], decimals, scaleFactor=100.0))
    haskey(columnIndex, :stddev) && (
      descContent[r][columnIndex[:stddev]] = num2Str(outdf[r, :stddev], decimals, scaleFactor=100.0))
    haskey(columnIndex, :skew) && (
      descContent[r][columnIndex[:skew]] = num2Str(outdf[r, :skew], decimals, scaleFactor=1.0))
    haskey(columnIndex, :kurt) && (
      descContent[r][columnIndex[:kurt]] = num2Str(outdf[r, :kurt], decimals, scaleFactor=1.0))
    haskey(columnIndex, :p5) && (
      descContent[r][columnIndex[:p5]] = num2Str(outdf[r, :p5], decimals, scaleFactor=100.0))
    haskey(columnIndex, :p25) && (
      descContent[r][columnIndex[:p25]] = num2Str(outdf[r, :p25], decimals, scaleFactor=100.0))
    haskey(columnIndex, :p50) && (
      descContent[r][columnIndex[:p50]] = num2Str(outdf[r, :p50], decimals, scaleFactor=100.0))
    haskey(columnIndex, :p75) && (
      descContent[r][columnIndex[:p75]] = num2Str(outdf[r, :p75], decimals, scaleFactor=100.0))
    haskey(columnIndex, :p95) && (
      descContent[r][columnIndex[:p95]] = num2Str(outdf[r, :p95], decimals, scaleFactor=100.0))
  end

  returnByYearTable::String = texTable(titleCaption="Descriptive Return Data by Year",
      caption="""See Tex File""", #caption
      colNames= columnNames, #colNames
      descRowNames=descRowNames, #descRowNames
      descContent = descContent, #descContent
      nakedTable=true
    )

    return returnByYearTable
end

function returnsSummaryDF(data::NCCSData; summaryPeriod::UnitRange = RETURNS_SUMMARY_PERIOD,
  focalField::Symbol = :return, decimals::Int = DECIMALS_RETURN,
  byCharityField::Symbol = BY_CHARITY_FIELD,
  byCategoryField::Symbol = BY_CATEGORY_FIELD, assetsField::Symbol = :adjnetassets)::DataFrame

  local subdf::SubDataFrame =  view(data.df,  (y::Int-> y ∈ summaryPeriod).(data.df[:fisyr]),
    [:ein; :fisyr; :fullsample; byCharityField; byCategoryField; focalField; assetsField])
  local byCharityValues::Vector{Symbol} = unique(processmissing(subdf[byCharityField]))
  local byCategoryValues::Vector{Symbol} = unique(processmissing(subdf[byCategoryField]))
  local rowFields::Vector{Symbol} = [byCharityValues; byCategoryValues; :fullsample]
  local nrows::Int = length(rowFields) # of rows we will have

  local N::Vector{MInt} = Vector{MInt}(undef, nrows)
  local minyear::Vector{MInt} = Vector{MInt}(undef, nrows)
  local maxyear::Vector{MInt} = Vector{MInt}(undef, nrows)
  local equalmean::Vector{MFloat64} = Vector{MFloat64}(undef, nrows)
  local valuemean::Vector{MFloat64} = Vector{MFloat64}(undef, nrows)
  local stddev::Vector{MFloat64} = Vector{MFloat64}(undef, nrows)
  local skew::Vector{MFloat64} = Vector{MFloat64}(undef, nrows)
  local kurt::Vector{MFloat64} = Vector{MFloat64}(undef, nrows)
  local p5::Vector{MFloat64} = Vector{MFloat64}(undef, nrows)
  local p25::Vector{MFloat64} = Vector{MFloat64}(undef, nrows)
  local p50::Vector{MFloat64} = Vector{MFloat64}(undef, nrows)
  local p75::Vector{MFloat64} = Vector{MFloat64}(undef, nrows)
  local p95::Vector{MFloat64} = Vector{MFloat64}(undef, nrows)

  #performance optimized way to split up the dataframe
  local outdfIndex::Dict = Dict(rowFields[i]=>i for i::Int ∈ 1:nrows)
  local rowNames::Vector{Symbol} =
    (s::Symbol->haskey(charitytype, s) ? charitytype[s] : s).(rowFields) #clarify some of the names
  local outdf::DataFrame = DataFrame(rownames = rowNames) #this will hold the summary stats

  #for each field we are considering
  for f::Symbol ∈ [byCharityField; byCategoryField; :fullsample]
    subdfByf::GroupedDataFrame = groupby(subdf, f)
    for ssubdf::SubDataFrame ∈ subdfByf
    #by(subdf, f) do ssubdf::SubDataFrame
      if !ismissing(ssubdf[1, f])
        r::Int = outdfIndex[ssubdf[1, f]]

        focal::Vector = ssubdf[:, focalField]
        nfirms::Int = length(focal)
        focal = processmissing(focal) #handle missing values

        N[r] = length(focal)
        minyear[r] = minimum(processmissing(ssubdf[:fisyr]))
        maxyear[r] = maximum(processmissing(ssubdf[:fisyr]))
        equalmean[r] = StatsBase.mean(focal)
        valuemean[r] = #compute the value-weighted average
          sum(processmissing(ssubdf[:, focalField] .* ssubdf[:, assetsField])) /
          sum(processmissing(ssubdf[(!ismissing).(ssubdf[:, focalField]), assetsField]))
        stddev[r] = StatsBase.std(focal)
        skew[r] = StatsBase.skewness(focal)
        kurt[r] = StatsBase.kurtosis(focal)
        p5[r] = StatsBase.quantile(focal, .05)
        p25[r] = StatsBase.quantile(focal, .25)
        p50[r] = StatsBase.quantile(focal, 0.5)
        p75[r] = StatsBase.quantile(focal, 0.75)
        p95[r] = StatsBase.quantile(focal, 0.95)
      end
    end
  end

  #now put the data into the dataframe
  outdf[:N] = N
  outdf[:minyear] = minyear
  outdf[:maxyear] = maxyear
  outdf[:equalmean] = equalmean
  outdf[:valuemean] = valuemean
  outdf[:stddev] = stddev
  outdf[:skew] = skew
  outdf[:kurt] = kurt
  outdf[:p5] = p5
  outdf[:p25] = p25
  outdf[:p50] = p50
  outdf[:p75] = p75
  outdf[:p95] = p95

  return outdf
end

function returnsSummary(outdf::DataFrame; summaryPeriod::UnitRange = RETURNS_SUMMARY_PERIOD,
  focalField::Symbol = :return, decimals::Int = DECIMALS_RETURN,
  byCharityField::Symbol = BY_CHARITY_FIELD,
  byCategoryField::Symbol = BY_CATEGORY_FIELD)::String

  columnSymbols::Vector{Symbol} = setdiff(names(outdf), [:rownames])

  numColumns = length(columnSymbols) #number of columns excluding the year
  columnIndex::Dict = Dict(columnSymbols[i]=>i for i ∈ 1:numColumns)

  #include all columns in the names
  columnNames::Vector{Vector{String}} = [(string).(columnSymbols)]
  descRowNames::Vector{String} = (string).(outdf[:rownames])

  #NOTE: a couple of hacks to improve readability
  descRowNames[length(keys(charitytype))+1] = "\\midrule $(descRowNames[length(keys(charitytype))+1])"
  descRowNames[end] = "\\midrule $(descRowNames[end])"

  numRows = size(outdf,1)

  #need to print the descriptive rows
  descContent::Vector{Vector{String}} = ((i::Int)->Vector{String}(undef, numColumns)).(1:length(descRowNames))

  for r::Int ∈ 1:numRows
    haskey(columnIndex, :N) && (
      descContent[r][columnIndex[:N]] = "$(outdf[r, :N])")
    haskey(columnIndex, :minyear) && (
      descContent[r][columnIndex[:minyear]] = "$(outdf[r, :minyear])")
    haskey(columnIndex, :maxyear) && (
      descContent[r][columnIndex[:maxyear]] = "$(outdf[r, :maxyear])")
    haskey(columnIndex, :equalmean) && (
      descContent[r][columnIndex[:equalmean]] = num2Str(outdf[r, :equalmean], decimals, scaleFactor=100.0))
    haskey(columnIndex, :valuemean) && (
      descContent[r][columnIndex[:valuemean]] = num2Str(outdf[r, :valuemean], decimals, scaleFactor=100.0))
    haskey(columnIndex, :stddev) && (
      descContent[r][columnIndex[:stddev]] = num2Str(outdf[r, :stddev], decimals, scaleFactor=100.0))
    haskey(columnIndex, :skew) && (
      descContent[r][columnIndex[:skew]] = num2Str(outdf[r, :skew], decimals, scaleFactor=1.0))
    haskey(columnIndex, :kurt) && (
      descContent[r][columnIndex[:kurt]] = num2Str(outdf[r, :kurt], decimals, scaleFactor=1.0))
    haskey(columnIndex, :p5) && (
      descContent[r][columnIndex[:p5]] = num2Str(outdf[r, :p5], decimals, scaleFactor=100.0))
    haskey(columnIndex, :p25) && (
      descContent[r][columnIndex[:p25]] = num2Str(outdf[r, :p25], decimals, scaleFactor=100.0))
    haskey(columnIndex, :p50) && (
      descContent[r][columnIndex[:p50]] = num2Str(outdf[r, :p50], decimals, scaleFactor=100.0))
    haskey(columnIndex, :p75) && (
      descContent[r][columnIndex[:p75]] = num2Str(outdf[r, :p75], decimals, scaleFactor=100.0))
    haskey(columnIndex, :p95) && (
      descContent[r][columnIndex[:p95]] = num2Str(outdf[r, :p95], decimals, scaleFactor=100.0))
  end

  assetsByYearTable::String = texTable(titleCaption="Descriptive Return Data by Charity Type and Category",
      caption="""See Tex File""", #caption
      colNames= columnNames, #colNames
      descRowNames=descRowNames, #descRowNames
      descContent = descContent, #descContent
      nakedTable=true
    )

    return assetsByYearTable
end

function getAssetsTables(data::NCCSData; summaryPeriod::UnitRange = ASSETS_SUMMARY_PERIOD,
  focalField::Symbol = :adjnetassets, sorted::Bool=false, outputPath::String = TABLE_PATH,
  headerName::String = HEADER_NAME, footerName::String = FOOTER_NAME)::Nothing

  #sort if we need to
  (!sorted) && sort!(data.df, [:fisyr, :adjnetassets], (:fisyr, order(:adjnetassets, rev=true)))

  #we get a full table and a shorter table for the decks
  #start with getting the wealth info by year
  outassetsyeardf::DataFrame = assetsByYearDF(data, focalField=focalField, summaryPeriod=summaryPeriod)
  assetsByYearTable::String = assetsByYear(outassetsyeardf, focalField=focalField, summaryPeriod=summaryPeriod)
  assetsByYearTableShort::String = assetsByYear(
    outassetsyeardf[[:fisyr, :N, :total, :mean, :stddev, :p5, :p50, :p95, :wealth80, :wealth95]],
    focalField=focalField, summaryPeriod=summaryPeriod)

  #now get the wealth info by category
  outassetssummarydf::DataFrame = assetsSummaryDF(data, focalField=focalField, summaryPeriod=summaryPeriod)
  assetsSummaryTable::String = assetsSummary(outassetssummarydf, focalField=focalField, summaryPeriod=summaryPeriod)
  assetsSummaryTableShort::String = assetsSummary(
    outassetssummarydf[[:rownames, :N, :minyear, :maxyear, :total, :mean, :stddev, :p5, :p50, :p95]],
    focalField=focalField, summaryPeriod=summaryPeriod)

  dataReligious::NCCSData=NCCSData(
    data.df[(s::MSymbol->((!ismissing(s)) && s==:ReligionRelated)).(data.df[:categoryname]),:])
  outreligousassetsdf::DataFrame = assetsByYearDF(dataReligious, focalField=focalField, summaryPeriod=summaryPeriod)
  religiousAssetsByYearTable::String =
    assetsByYear(outreligousassetsdf, focalField=focalField, summaryPeriod=summaryPeriod)

  #write the tables
  writeNakedTable(assetsByYearTable, path=outputPath,  outName = "assetsTableByYear.tex")
  writeNakedTable(assetsByYearTableShort, path=outputPath,  outName = "assetsTableByYearShort.tex")
  writeNakedTable(assetsSummaryTable, path=outputPath,  outName = "assetsSummaryTable.tex")
  writeNakedTable(assetsSummaryTableShort, path=outputPath,  outName = "assetsSummaryTableShort.tex")
  writeNakedTable(religiousAssetsByYearTable, path=outputPath,
    outName = "assetsReligiousTable.tex")

  return nothing
end

function getReturnsTables(data::NCCSData; summaryPeriod::UnitRange = RETURNS_SUMMARY_PERIOD,
  focalField::Symbol = :return, sorted::Bool=false, outputPath::String = TABLE_PATH,
  headerName::String = HEADER_NAME, footerName::String = FOOTER_NAME)::Nothing

  #sort if we need to
  (!sorted) && sort!(data.df, [:fisyr, :adjnetassets], (:fisyr, order(:adjnetassets, rev=true)))

  outreturnsyeardf::DataFrame = returnsByYearDF(data, focalField=focalField, summaryPeriod=summaryPeriod)
  returnsByYearTable::String = returnsByYear(outreturnsyeardf, focalField=focalField, summaryPeriod=summaryPeriod)
  returnsByYearTableShort::String = returnsByYear(
    outreturnsyeardf[[:fisyr, :N, :equalmean, :valuemean, :stddev, :skew, :p5, :p50, :p95]],
    focalField=focalField, summaryPeriod=summaryPeriod)

  outreturnssummarydf::DataFrame = returnsSummaryDF(data, focalField=focalField, summaryPeriod=summaryPeriod)
  returnsSummaryTable::String = returnsSummary(outreturnssummarydf, focalField=focalField, summaryPeriod=summaryPeriod)
  returnsSummaryTableShort::String = returnsSummary(
    outreturnssummarydf[[:rownames, :N, :equalmean, :valuemean, :p5, :p50, :p95]],
    focalField=focalField, summaryPeriod=summaryPeriod)

  #write the tables
  writeNakedTable(returnsByYearTable, path=outputPath,  outName = "returnsTableByYear.tex")
  writeNakedTable(returnsByYearTableShort, path=outputPath,  outName = "returnsTableByYearShort.tex")
  writeNakedTable(returnsSummaryTable, path=outputPath,  outName = "returnsSummaryTable.tex")
  writeNakedTable(returnsSummaryTableShort, path=outputPath,  outName = "returnsSummaryTableShort.tex")

  return nothing
end

#NOTE: Consider slicing by type with collors being the years

function makeDescriptiveTables(data::NCCSData)::Nothing

  sort!(data.df, (:fisyr, order(:adjnetassets, rev=true)))

  data.df[:fullsample] = :fullsample #this is a hack for analyzing a full sample

  getReturnsTables(data, sorted=true)
  getAssetsTables(data, sorted=true)

  deletecols!(data.df, :fullsample)

  return nothing
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
__precompile__()
module DiscreteRanges

using ArgCheck: @argcheck
using Dates: Date, days, Day

import Base:
    show, isequal, ==, hash, convert,
    in, isempty, issubset, union, intersect,
    minimum, maximum, extrema,
    length, size, IndexStyle, getindex, start, next, done, IteratorSize,
    IteratorEltype

export DiscreteRange, ..

"""
    DiscreteRange(left::T, right::T) where T
    DiscreteRange(singleton)

An object representing all values `x::T` such that `left ≤ x ≤ right`. Supports
`∈`, `⊆`, indexing, iteration, and (immutable) array interfaces. For integers,
it is equivalent to `UnitRange{T}`. With a single argument, `left == right`.

## Definition for custom types

The package defines methods for `<: Integer` and `Date`. To extend this, define
the following: [`isdiscrete`](@ref), [`discrete_gap`](@ref),
[`discrete_next`](@ref), [`Base.isless`](@ref).

The methods for [`Base.==`](@ref), [`Base.isequal`](@ref) and
[`Base.hash`](@ref) only work for `DiscreteRange{T}` when they are defined for
`T`. The unit tests of the
"""
struct DiscreteRange{T} <: AbstractVector{T}
    left::T
    right::T
    function DiscreteRange(left::T, right::T) where T
        @argcheck isdiscrete(T) "Only discrete types are supported. Define `isdiscrete($T)`."
        new{T}(left, right)
    end
end

convert(::Type{DiscreteRange{T}}, D::DiscreteRange{T}) where T = D

convert(::Type{DiscreteRange{T}}, D::DiscreteRange) where T =
    DiscreteRange(T(D.left), T(D.right))

..(left, right) = DiscreteRange(left, right)

"""
    isdiscrete(T)

Creation of `DiscreteRange{T}` objects is only allowed when this returns `true`.

Needs to be defined for valid `DiscreteRange` type parameters.
"""
isdiscrete(::Type) = false
isdiscrete(::Type{<:Integer}) = true
isdiscrete(::Type{Date}) = true

"""
    discrete_gap(x::T, y::T)

Return the difference `y - x` as an *integer*.

Needs to be defined for valid `DiscreteRange` type parameters.
"""
discrete_gap(x::Integer, y::Integer) = x - y
discrete_gap(x::Date, y::Date) = days(x - y)

"""
    discrete_next(x::T, Δ::Integer)

Return `y::T` such that `discrete_gap(y, x) == Δ`.

Needs to be defined for valid `DiscreteRange` type parameters.
"""
discrete_next(x) = discrete_next(x, 1)
discrete_next(x::T, Δ) where {T <: Integer} = x + T(Δ)
discrete_next(x::Date, Δ) = x + Day(Δ)

DiscreteRange(x) = DiscreteRange(x, x) # single value

isempty(D::DiscreteRange) = D.left > D.right

function show(io::IO, D::DiscreteRange{T}) where T
    if isempty(D)
        print(io, "empty DiscreteRange{$(T)}")
    else
        print(io, D.left, "..", D.right)
    end
end

show(io::IO, ::MIME"text/plain", D::DiscreteRange) =
    show(io, D) # fallback to prevent AbstractArray printing

in(x::T, D::DiscreteRange{T}) where T = D.left ≤ x ≤ D.right

in(x, D::DiscreteRange{T}) where T = T(x) ∈ D

isequal(A::DiscreteRange, B::DiscreteRange) =
    (isequal(A.left, B.left) & isequal(A.right, B.right)) | (isempty(A) & isempty(B))

==(A::DiscreteRange, B::DiscreteRange) =
    (A.left == B.left) & (A.right == B.right) | (isempty(A) & isempty(B))

const _DISCRETERANGE_HASH = UInt === UInt64 ? 0x3c467b76e5939730 : 0xff9d6b2d

hash(A::DiscreteRange, h::UInt) = hash(A.left, hash(A.right, hash(_DISCRETERANGE_HASH, h)))

function minimum(D::DiscreteRange)
    @argcheck !isempty(D)
    D.left
end

function maximum(D::DiscreteRange)
    @argcheck !isempty(D)
    D.right
end

function extrema(D::DiscreteRange)
    @argcheck !isempty(D)
    D.left, D.right
end

intersect(A::DiscreteRange, B::DiscreteRange) =
    DiscreteRange(max(A.left, B.left), min(A.right, B.right))

function union(A::DiscreteRange, B::DiscreteRange)
    @argcheck(discrete_gap(max(A.left, B.left), min(A.right, B.right)) ≤ 1,
              "Disjoint discrete ranges.")
    DiscreteRange(min(A.left, B.left), max(A.right, B.right))
end

issubset(A::DiscreteRange{T}, B::DiscreteRange{T}) where {T} =
    isempty(A) || (B.left ≤ A.left ≤ A.right ≤ B.right)

issubset(A::DiscreteRange, B::DiscreteRange) = issubset(oftype(B, A), B)

length(D::DiscreteRange) = isempty(D) ? 0 : Int(discrete_gap(D.right, D.left)) + 1

size(D::DiscreteRange) = (length(D),)

IndexStyle(::Type{<:DiscreteRange}) = IndexLinear()

function getindex(D::DiscreteRange, i)
    i ≥ 1 || throw(BoundsError(D, i))
    elt = discrete_next(D.left, i-1)
    elt ≤ D.right || throw(BoundsError(D, i))
    elt
end

start(D::DiscreteRange) = D.left

next(D::DiscreteRange, state) = state, discrete_next(state)

done(D::DiscreteRange, state) = state > D.right

IteratorSize(::Type{DiscreteRange{T}}) where T = Base.HasLength()

IteratorEltype(::Type{DiscreteRange{T}}) where T = Base.HasEltype()

convert(::Type{R}, D::DiscreteRange) where {R <: StepRange} =
    R(D.left, discrete_next(D.left) - D.left, D.right)

end # module


text:
using Plots
using LaTeXStrings
t = 0:.001:.999999
N = 1 ./(exp.(1 .- t) .- t)
plot(t, N,label="", yaxis="ideal population (log scale)",xaxis=L"\varepsilon")
yaxis!(:log)
savefig("zeroMatters.jpg")


text:
# # Forest fire model

# ![](forest.gif)

# The forest fire model is defined as a cellular automaton on a grid.
# A position can be empty or occupied by a tree which is ok, burning or burnt.
# We implement a slightly different ruleset to that of
# [Drossel and Schwabl (1992)](https://en.wikipedia.org/wiki/Forest-fire_model),
# so that our implementation can be compared with other ABM frameworks
#
# 1. A burning position turns into a burnt position
# 1. A tree will burn if at least one neighbor is burning

# The forest has an innate `density`, which is the proportion of trees initialized as
# `green`, however all trees that reside at `x=1` on the grid are `burning`.
# The model is also available from the `Models` module as [`Models.forest_fire`](@ref).

# ## Defining the core structures

# We start by defining the agent type
using Agents, Random, Plots
gr() # hide

mutable struct Tree <: AbstractAgent
    id::Int
    pos::Dims{2}
    status::Symbol  #:green, :burning, :burnt
end
nothing # hide

# The agent type `Tree` has three fields: `id` and `pos`, which have to be there for any agent,
# and a `status` field that we introduce for this specific model.
# The `status` field will be `:green` when the tree is ok, `:burning` when on fire,
# and finally `:burnt`.

# We then make a setup function that initializes the model.
function forest_fire(; density = 0.7, griddims = (100, 100))
    space = GridSpace(griddims; periodic = false, metric = :euclidean)
    forest = AgentBasedModel(Tree, space)
    ## create and add trees to each position with a probability
    ## determined by the `density`.
    for position in positions(forest)
        if rand() < density
            ## Set the trees at position x=1 on fire
            state = position[1] == 1 ? :burning : :green
            add_agent!(position, forest, state)
        end
    end
    return forest
end

forest = forest_fire()

# ## Defining the step!
# Because of the way the forest fire model is defined, we only need a
# stepping function for the agents

function tree_step!(tree, forest)
    ## The current tree is burning
    if tree.status == :burning
        ## Find all green neighbors and set them on fire
        for neighbor in nearby_agents(tree, forest)
            if neighbor.status == :green
                neighbor.status = :burning
            end
        end
        tree.status = :burnt
    end
end
nothing # hide

# ## Running the model

step!(forest, tree_step!, 1)
count(t->t.status == :burnt, allagents(forest))

#

step!(forest, tree_step!, 10)
count(t->t.status == :burnt, allagents(forest))

# Now we can do some data collection as well using an aggregate function `percentage`:

Random.seed!(2)
forest = forest_fire(griddims = (20, 20))
burnt_percentage(m) = count(t->t.status == :burnt, allagents(m)) / length(positions(m))
mdata = [burnt_percentage]

_, data = run!(forest, tree_step!, 10; mdata)
data

# Now let's plot the model. We use green for unburnt trees, red for burning and a
# dark red for burnt.
forest = forest_fire()
step!(forest, tree_step!, 1)

function treecolor(a)
    color = :green
    if a.status == :burning
        color = :red
    elseif a.status == :burnt
        color = :darkred
    end
    color
end

plotabm(forest; ac = treecolor, ms = 5)

# or animate it
cd(@__DIR__) #src
forest = forest_fire(density = 0.6)
anim = @animate for i in 0:10
    i > 0 && step!(forest, tree_step!, 5)
    p1 = plotabm(forest; ac = treecolor, ms = 5, msw = 0)
    title!(p1, "step $(i)")
end

gif(anim, "forest.gif", fps = 2)


text:
# Run line by line in the REPL, or ctrl + enter in VSCode
# First code section
using MechanicalUnits
m_air = 1000kg; c_p = 1.00kJ/(kg*K)
@import_expand ~W   # Watt = Joule / Second is not exported by default. Several: (u1, u2,..)
Q_cp(T1, T2) = m_air * c_p * (T2 - T1) |> (kW*h)
Q_cp(20°C, 985°C)
dm |> upreferred
preferunits(m) # No effect, since upreferred was called once this session
m_s = [30kg/m 28.8lb/ft]
l_s = 93ft*[3 4]m/s
m_s.*l_s .|> (kg*m)
E=206GPa; h_y = 100mm; b = 30mm; I = 1/12 * b * h_y^3
L = 2m; F=100kg*g |> N
F*L^3/(3E*I) |> upreferred
l_wire = 20m
k(d) = E * 0.691 * π/4 * d^2 / l_wire |> N/mm
k.([5 6 8]mm)
δ(d)= F / k(d) |> mm
δ.([5, 6, 8]mm)
d = 6mm
dimension(d)
1d |> s
@import_expand ~V ~W ~A  G
sqrt(1G²)
[1V*12.0A 2W 1kg*g*1m/2s]*30minute .|> kJ
ω = 50*2π*rad/s
t = (0:0.006:0.02)s
u = 220V*exp.(im∙(ω∙t))
u*1.5A .|> J

# Second code section
import MechanicalUnits: @import_expand, ∙
@import_expand ~m dyn     # ~ : also import SI prefixes for metre
(1.0cm², 2.0mm∙m, 3.0dm⁴/m² ) .|> mm²
typeof(dyn)
1dyn |> μm

# Third code section
strinp = "2 [s]\t11364.56982421875 [N]\t-44553.50244140625 [N]\t-26.586366176605225 [N]\t0.0[N mm]\t0.0[N mm]\t0.0[N mm]\t1561.00350618362 [mm]\t-6072.3729133606 [mm]\t2825.15907287598 [mm]";

time, Fx, Fy, Fz, Mx, My, Mz, px, py, pz = parse.(Quantity{Float64}, split(strinp, '\t'))

# Fourth code section
strain = 10.6μm/m
strain |> upreferred
strain *m/μm
strain |> NoUnits

max_stars_repo_name: hustf/MechanicalUnits

max_issues_repo_name: hustf/MechanicalUnits


text:
abstract InterpType
type InterpForward <: InterpType; end
type InterpBackward <: InterpType; end
type InterpNearest <: InterpType; end
type InterpLinear <: InterpType; end
type InterpQuadratic <: InterpType; end
type InterpCubic <: InterpType; end


text:
"""
# module TypeGraph

Taken from [Google Groups](https://groups.google.com/forum/#!msg/julia-users/ECkQa8vAkko/6tqVeWrexr4J).
A short julia program that generates the graph of all subtypes of a given type, in png format by default.
Uses `graphviz`.

# Examples

```jldoctest
julia> makegraph(Number, "result.dot")
=> Result here: https://imgur.com/a/CkoOA#0
julia> makegraph(Any, "alltypes.dot", "svg")
=> Result here: https://imgur.com/U5vTDNp
```
"""
module TypeGraph

using InteractiveUtils

export makegraph

isalnum(c) = isletter(c) || isnumeric(c)

function typeit(fn, typ)
    seen = Dict{Type,Bool}()
    function apply(fn, tp)
        if !get(seen, tp, false)
            seen[tp] = true
            for child in subtypes(tp)
                fn(tp, child)
                apply(fn, child)
            end
        end
    end
    apply(fn, typ)
    collect(keys(seen))
end

function edge(typ1, typ2)
    let convert = typ->map(x->if isalnum(x) x else '_' end, string(typ))
        c1 = convert(typ1)
        c2 = convert(typ2)
        if c2 == "Graph" ; c2 = "G;\nG [label=\"Graph\"]" ; end
        "$(c1) -> $(c2);\n"
    end
end

const graphviz_formats = ["png","gif","svg","jpg","svgz","pcl", "mif", "pcl", "dia", "ps","fig", "imap","cmapx"]

function makegraph(typ, file, format = "png")
    if !in(format, graphviz_formats)
        println("Warning: the output format may not be supported by graphviz.\nHere are the default formats:\n$(join(graphviz_formats, " "))")
    end
    open(file, "w") do f
        write(f, "digraph{\n")                # oriented graph
        write(f, "rankdir=LR\n")              # layout of the graph from left to right
        typeit((x, y)->write(f, edge(x, y)), typ)
        write(f, "}")
    end
    run(`dot -T$(format) -O $file`)        # graphviz needed here.
end

end # module


text:
function state_to_msg(multicopter::Multicopter, x)
    # Msg = pyimport("geometry_msgs.msg")
    Msg = pyimport("fsim_interfaces.msg")  # TODO
    @unpack p, v, R, ω = x
    quat = dcm_to_quat(SMatrix{3, 3}(R))
    # pose
    msg = Msg.PoseTwist()
    msg.pose.position.x, msg.pose.position.y, msg.pose.position.z = p
    msg.pose.orientation.x, msg.pose.orientation.y, msg.pose.orientation.z, msg.pose.orientation.w = quat.q1, quat.q2, quat.q3, quat.q0  # be careful for the order
    # twist
    msg.twist.linear.x, msg.twist.linear.y, msg.twist.linear.z = v
    msg.twist.angular.x, msg.twist.angular.y, msg.twist.angular.z = ω
    msg
end

function msg_to_state(multicopter::Multicopter, msg)
    Msg = pyimport("fsim_interfaces.msg")
    p = [msg.pose.position.x, msg.pose.position.y, msg.pose.position.z]
    v = [msg.twist.linear.x, msg.twist.linear.y, msg.twist.linear.z]
    quat = Quaternion(msg.pose.orientation.w, msg.pose.orientation.x, msg.pose.orientation.y, msg.pose.orientation.z)  # be careful for the order
    R = quat_to_dcm(quat)
    ω = [msg.twist.angular.x, msg.twist.angular.y, msg.twist.angular.z]
    state = State(multicopter)(p, v, R, ω)
end


text:
using Revise
using Plots
using Pkg
# Pkg.activate(".")
using RadiativeTransfer
using RadiativeTransfer.Architectures
using RadiativeTransfer.Absorption
using RadiativeTransfer.Scattering
using RadiativeTransfer.vSmartMOM
using RadiativeTransfer.SolarModel
using InstrumentOperator
using Interpolations
using Polynomials
using ForwardDiff
using Distributions
using NCDatasets

## Atmospheric Radiative Transfer

# Load parameters from file
parameters = vSmartMOM.parameters_from_yaml("test/test_parameters/O2Parameters.yaml")
#parameters.architecture = CPU()
FT = Float64

# Load OCO Data:
# File names:
L1File = "/net/fluo/data1/group/oco2/L1bSc/oco2_L1bScND_26780a_190715_B10003r_200429212407.h5"
metFile = "/net/fluo/data1/group/oco2/L2Met/oco2_L2MetND_26780a_190715_B10003r_200429212406.h5"
dictFile = "/home/cfranken/code/gitHub/InstrumentOperator.jl/json/oco2.yaml"

# Load L1 file (could just use filenames here as well)
oco = InstrumentOperator.load_L1(dictFile,L1File, metFile);

# oco2_L1bScND_18688a_180105_B8100r_180206190633

oco_file = "/net/fluo/data1/group/oco2/L1bSc/oco2_L1bScGL_15258a_170515_B10003r_200214061601.h5"
oco_met_file = "/net/fluo/data1/group/oco2/L2Met/oco2_L2MetGL_26777a_190715_B10003r_200429213029.h5"
ils_file = "test/prototyping/ils_oco2.json"

# Pick some bands as tuple (or just one)
bands = (1,);
# Indices within that band:
indices = (92:885,);
# Geo Index (footprint,sounding):
GeoInd = [5,5000];

# Get data for that sounding:
oco_sounding = InstrumentOperator.getMeasurement(oco, bands, indices, GeoInd)

fp = 5
iOrbit = 5000
band = 1

sza_ = 32.4436
vza_ = [0.072]

function conv_spectra_local(m::VariableKernelInstrument, ν, spectrum; stride=1)
    # FT = eltype(m.ν_out)
    # Define grid where to perform convolution:

    # Padding at both sides required:
    off = ceil(Int, size(m.kernel, 1) / 2)
    ind = off:stride:(length(ν) - off)

    # knots where convolution will be applied to
    knots = view(ν, ind)
    te = LinearInterpolation(m.ν_out, Float32.(m.ind_out))
    spec_out = zeros(Real, length(knots));
    for i in eachindex(knots)
        # Simple first, nearest neighbor ILS
        ind_fraction = round(Int, te(knots[i]));
        kernel = view(m.kernel, :, ind_fraction)
        for j in eachindex(kernel)
            spec_out[i] += kernel[j] * spectrum[ind[i] + j]
        end
    end
    # Change this later to only perform conv around output grid!
    fin = LinearInterpolation(ν[ind], spec_out; extrapolation_bc=Interpolations.Flat())
    return fin(m.ν_out)
end;

# Runner is used to set AD fields as duals
function runner!(y, x, parameters=parameters, oco_file=oco_file,
                                              oco_met_file=oco_met_file,
                                              ils_file=ils_file, fp=fp, iOrbit=iOrbit,
                                              sza_=sza_, vza_=vza_)

    # Set parameters fields as the dual numbers
    parameters.brdf = [vSmartMOM.LambertianSurfaceLegendre([x[1],x[3],x[4]])]

    parameters.scattering_params.rt_aerosols[1].τ_ref = exp(x[2]);
    parameters.scattering_params.rt_aerosols[1].p₀    = 70000.0; #x[4]
    # parameters.scattering_params.rt_aerosols[1].aerosol.size_distribution = LogNormal(log(x[3]), log(x[4]), check_args=false)

    # parameters.scattering_params.rt_aerosols[1].aerosol.nᵣ = x[5];
    # parameters.scattering_params.rt_aerosols[1].aerosol.nᵢ = x[6];

    # parameters.scattering_params.rt_aerosols[1].p₀ = x[7];
    # parameters.scattering_params.rt_aerosols[1].σp = x[8];
    ocoData = Dataset(oco_file)
    SoundingGeometry = ocoData.group["SoundingGeometry"]
    ϕ = SoundingGeometry["sounding_polarization_angle"][fp,iOrbit]
    m₁ = 0.5
    m₂ = cosd(2ϕ)/2
    m₃ = sind(2ϕ)/2
    # Stokes:

    # Set profiles properly
    met = Dataset(oco_met_file);
    T_met = met.group["Meteorology"]["temperature_profile_met"][:,fp,iOrbit];
    ak = met.group["MeteorologyDiagnostics"]["ak"][:,fp,iOrbit];
    bk = met.group["MeteorologyDiagnostics"]["bk"][:,fp,iOrbit];
    ak = [0.000000e+00, 4.804826e-02, 6.593752e+00, 1.313480e+01, 1.961311e+01, 2.609201e+01,
3.257081e+01, 3.898201e+01, 4.533901e+01, 5.169611e+01, 5.805321e+01, 6.436264e+01,
7.062198e+01, 7.883422e+01, 8.909992e+01, 9.936521e+01, 1.091817e+02, 1.189586e+02,
1.286959e+02, 1.429100e+02, 1.562600e+02, 1.696090e+02, 1.816190e+02, 1.930970e+02,
2.032590e+02, 2.121500e+02, 2.187760e+02, 2.238980e+02, 2.243630e+02, 2.168650e+02,
2.011920e+02, 1.769300e+02, 1.503930e+02, 1.278370e+02, 1.086630e+02, 9.236572e+01,
7.851231e+01, 6.660341e+01, 5.638791e+01, 4.764391e+01, 4.017541e+01, 3.381001e+01,
2.836781e+01, 2.373041e+01, 1.979160e+01, 1.645710e+01, 1.364340e+01, 1.127690e+01,
9.292942e+00, 7.619842e+00, 6.216801e+00, 5.046801e+00, 4.076571e+00, 3.276431e+00,
2.620211e+00, 2.084970e+00, 1.650790e+00, 1.300510e+00, 1.019440e+00, 7.951341e-01,
6.167791e-01, 4.758061e-01, 3.650411e-01, 2.785261e-01, 2.113490e-01, 1.594950e-01,
1.197030e-01, 8.934502e-02, 6.600001e-02, 4.758501e-02, 3.270000e-02, 2.000000e-02,
1.000000e-02][end:-1:1]

bk = [1.000000e+00, 9.849520e-01, 9.634060e-01, 9.418650e-01, 9.203870e-01, 8.989080e-01,
8.774290e-01, 8.560180e-01, 8.346609e-01, 8.133039e-01, 7.919469e-01, 7.706375e-01,
7.493782e-01, 7.211660e-01, 6.858999e-01, 6.506349e-01, 6.158184e-01, 5.810415e-01,
5.463042e-01, 4.945902e-01, 4.437402e-01, 3.928911e-01, 3.433811e-01, 2.944031e-01,
2.467411e-01, 2.003501e-01, 1.562241e-01, 1.136021e-01, 6.372006e-02, 2.801004e-02,
6.960025e-03, 8.175413e-09, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00,
0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00,
0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00,
0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00,
0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00,
0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00,
0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00,
0.000000e+00][end:-1:1]

    p_surf = met.group["Meteorology"]["surface_pressure_met"][fp,iOrbit];
    p_half = (ak + bk * p_surf);
    #p_half = vcat(p_half, p_surf);
    q = met.group["Meteorology"]["specific_humidity_profile_met"][:,fp,iOrbit];
    parameters.p = p_half / 100
    parameters.q = q # zeros(size(q))
    parameters.T = T_met
    #parameters.T[end-30:end] .-= 10
    parameters.sza = sza_
    parameters.vza = vza_

    @show parameters.sza, parameters.vza

    # @show parameters.sza
    # plot!(parameters.p)
    # return

    model = model_from_parameters(parameters);
    @show sum(model.τ_rayl[1] )
    # Run the model to obtain reflectance matrix
    R = vSmartMOM.rt_run(model, i_band=1);
    RR = m₁*R[1,1,:] + m₂*R[1,2,:] + m₃*R[1,3,:]
    @show size(R)
    # Produce black-body in wavenumber range
    T = 5777
    λ_grid = reverse(1e4 ./ parameters.spec_bands[1]) #collect(757:0.01:777)
    black_body = planck_spectrum_wl(T, λ_grid) * 2.1629e-05 * pi
    black_body = SolarModel.watts_to_photons(λ_grid, black_body)

    # Get solar transmittance spectrum
    #solar_transmission = solar_transmission_from_file("/home/rjeyaram/RadiativeTransfer/src/solar_merged_20160127_600_26316_100.out", parameters.spec_bands[1])
    solar_transmission = solar_transmission_from_file("/home/rjeyaram/RadiativeTransfer/src/SolarModel/solar.out", parameters.spec_bands[1])
    # Get outgoing solar radiation
    sun_out = reverse(solar_transmission) .* black_body

    # Apply Earth reflectance matrix
    earth_out = sun_out .* reverse(RR[:])

    # Set up InstrumentOperator
    ils_Δ, ils_in, dispersion = InstrumentOperator.read_ils_table(oco_file, ils_file);

    # Define model grid:
    res = 0.001

    # Just consider the ILS within ± 0.35nm
    grid_x = -0.35e-3:res*1e-3:0.35e-3

    extended_dims = [fp,band] # Footprint, band

    # Re-interpolate I from ν_grid to new grid/resolution
    interp_I = LinearInterpolation(λ_grid, earth_out);
    wl = 757.5:res:771.0
    I_wl = interp_I(wl/1000)

    # Pixels to be used
    ind_out = collect(0:1015);

    # Eventual grid of OCO-2 for Band 1, FP 5:
    dispPoly = Polynomial(view(dispersion, :, extended_dims...))
    ν = Float32.(dispPoly.(1:1016))

    # Prepare ILS table:
    ils_pixel   = InstrumentOperator.prepare_ils_table(grid_x, ils_in, ils_Δ,extended_dims)
    oco2_kernel = VariableKernelInstrument(ils_pixel, ν, ind_out)

    # Convolve input spectrum with variable kernel
    I_conv = conv_spectra_local(oco2_kernel, wl./1e3, I_wl)

    y[:] = I_conv[:]

end

directory = "/net/fluo/data1/group/oco2/L1bSc/"
files_list = filter(x->endswith(x, ".h5") && startswith(x, "oco2_L1bScND_26"), readdir(directory))

met_files_list = filter(x->endswith(x, ".h5"), readdir("/net/fluo/data1/group/oco2/L2Met/"))

I_convs_all = zeros(1016, length(files_list))
oco2_Abands_all = zeros(1016, length(files_list))

albedos = zeros(length(files_list))
cloud_levels = zeros(length(files_list))

optimal_x = zeros(8, length(files_list))

#for i in 1:length(files_list)
i = 1
ind = 92:885
@show i

oco_file = files_list[i]

    # try
ocoData = Dataset(directory * oco_file)
oco_met_file = "/net/fluo/data1/group/oco2/L2Met/" * filter(x->contains(x, oco_file[14:19]), met_files_list)[1]
oco_file = directory * oco_file
met = Dataset(oco_met_file)
cloud_levels[i] = met.group["Meteorology"]["cloud_liquid_water_path_met"][fp,iOrbit]

## Getting an OCO spectrum for fun to fit:
o2AbandSpectra = ocoData.group["SoundingMeasurements"]["radiance_o2"]
SoundingGeometry = ocoData.group["SoundingGeometry"]
ϕ = SoundingGeometry["sounding_polarization_angle"][fp,iOrbit]
m₁ = 0.5
m₂ = cosd(2ϕ)/2
m₃ = sind(2ϕ)/2
vza = SoundingGeometry["sounding_zenith"]
sza = SoundingGeometry["sounding_solar_zenith"]
lat = SoundingGeometry["sounding_latitude"]
lon = SoundingGeometry["sounding_longitude"]
o2AbandSpectra = ocoData.group["SoundingMeasurements"]["radiance_o2"]
oco2_Aband = o2AbandSpectra[:,fp,iOrbit]

println("$(round(Float64(lat[fp,iOrbit]), digits=2)) $(round(Float64(lon[fp,iOrbit]), digits=4))")
# continue
sza_ = sza[fp,iOrbit]
vza_ = [vza[fp,iOrbit]]

# State vector
x = FT[0.2377,
        -1, -0.00244, 0]#,

# Run FW model:
# @time runner(x);
I_conv = zeros(1016)
y = oco2_Aband[ind];
for i=1:3
    dfdx = ForwardDiff.jacobian(runner!, I_conv, x);
    Fx = I_conv[ind];
    K = dfdx[ind,:];
    dx = K \ (y-Fx)
    x += dx
    @show x
end
runner!(I_conv,x)
Fx = I_conv[ind];

# plot fits
ils_Δ, ils_in, dispersion = InstrumentOperator.read_ils_table(oco_file, ils_file);
extended_dims = [fp,band]; # Footprint, band
# Eventual grid of OCO-2 for Band 1, FP 5:
dispPoly = Polynomial(view(dispersion, :, extended_dims...));
ν = Float32.(dispPoly.(1:1016));
plot(ν[ind]*1e3,  y/1e20, label="Meas")
plot!(ν[ind]*1e3, Fx/1e20, label="Mod")
plot!(ν[ind]*1e3, (y-Fx)/1e20, label="Meas-mod")

function sampleOCO_Orbit(fp,iOrbit,l1,met)

end


text:
#!/usr/bin/env julia6
using ReadGlobal

function makeTsgs!(r::Array{T,N},u::Array{T,N},v::Array{T,N}) where {T<:Number,N}
    Threads.@threads for i::Int in 1:length(u)
        @inbounds r[i] = r[i] - u[i]*v[i]
    end
    return r
end

function makeTsum!(r::Array{T,N},u::Array{T,N},v::Array{T,N}) where {T<:Number,N}
    Threads.@threads for i::Int in 1:length(u)
        @inbounds r[i] = r[i]+u[i]+v[i]
    end
    return r
end

function makeTtraceless!(r::Array{T,N},u::Array{T,N}) where {T<:Number,N}
    third::T = 1/3
    Threads.@threads for i::Int in 1:length(u)
        @inbounds r[i] = r[i] - third*u[i]
    end
    return r
end

function makeT()

    N, Fil = getnfilter()

    nx, ny, nz, lx, ly, lz = getdimsize()

    aux1 = readfield("$(Fil)u1u1_N$N",nx,ny,nz)
    aux2 = readfield("$(Fil)u1_N$N",nx,ny,nz)
    write("$(Fil)Tr11_N$N",makeTsgs!(aux1,aux2,aux2))

    read!("$(Fil)u1u2_N$N",aux1)
    aux3 = readfield("$(Fil)u2_N$N",nx,ny,nz)
    write("$(Fil)T12_N$N",makeTsgs!(aux1,aux2,aux3))

    read!("$(Fil)u1u3_N$N",aux1)
    read!("$(Fil)u3_N$N",aux3)
    write("$(Fil)T13_N$N",makeTsgs!(aux1,aux2,aux3))

    read!("$(Fil)u2u2_N$N",aux1)
    read!("$(Fil)u2_N$N",aux2)
    write("$(Fil)Tr22_N$N",makeTsgs!(aux1,aux2,aux2))

    read!("$(Fil)u2u3_N$N",aux1)
    write("$(Fil)T23_N$N",makeTsgs!(aux1,aux2,aux3))

    read!("$(Fil)u3u3_N$N",aux1)
    write("$(Fil)Tr33_N$N",makeTsgs!(aux1,aux3,aux3))

    read!("$(Fil)Tr11_N$N",aux2)
    read!("$(Fil)Tr22_N$N",aux3)

    write("2k",makeTsum!(aux1,aux2,aux3))

    write("$(Fil)T11_N$N",makeTtraceless!(aux2,aux1))
    write("$(Fil)T22_N$N",makeTtraceless!(aux3,aux1))

    read!("$(Fil)Tr33_N$N",aux3)
    write("$(Fil)T33_N$N",makeTtraceless!(aux3,aux1))

    return 0
end

makeT()


text:
using MibS_jll

function _build_single_model(model::BilevelModel, check_MIPMIP::Bool = false)
    upper = JuMP.backend(model.upper)
    lower = JuMP.backend(model.lower)
    lower_to_upper =
        Dict(JuMP.index(v) => JuMP.index(k) for (k, v) in model.link)
    lower_only = Dict(
        JuMP.index(k) => JuMP.index(v) for (k, v) in model.lower_to_upper_link
    )
    return _build_single_model(upper, lower, lower_to_upper, lower_only, check_MIPMIP)
end

function _build_single_model(
    upper::MOI.ModelLike,
    lower::MOI.ModelLike,
    lower_to_upper_link::Dict{MOI.VariableIndex,MOI.VariableIndex},
    lower_only::Dict{MOI.VariableIndex,MOI.VariableIndex},
    check_MIPMIP::Bool = false
)
    model = MOI.FileFormats.MPS.Model()
    upper_to_model_link = MOI.copy_to(model, upper)
    lower_variables = [upper_to_model_link[k] for k in values(lower_only)]
    lower_constraints = Vector{MOI.ConstraintIndex}()
    for (F, S) in MOI.get(lower, MOI.ListOfConstraints())
        for ci in MOI.get(lower, MOI.ListOfConstraintIndices{F,S}())
            lower_f = MOI.get(lower, MOI.ConstraintFunction(), ci)
            lower_s = MOI.get(lower, MOI.ConstraintSet(), ci)
            lower_f = MOI.Utilities.map_indices(lower_f) do x
                return upper_to_model_link[lower_to_upper_link[x]]
            end
            new_ci = MOI.add_constraint(model, lower_f, lower_s)
            if F == MOI.ScalarAffineFunction{Float64}
                push!(lower_constraints, new_ci)
            end
        end
    end
    lower_objective = MOI.get(
        lower,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
    )
    lower_objective = MOI.Utilities.map_indices(lower_objective) do x
        return upper_to_model_link[lower_to_upper_link[x]]
    end

    # Testing if the model is MIP-MIP or not.
    if check_MIPMIP
        int_var = MOI.get(model, MOI.NumberOfConstraints{MOI.SingleVariable, MOI.Integer}())
        int_var = int_var + MOI.get(model, MOI.NumberOfConstraints{MOI.SingleVariable, MOI.ZeroOne}())
        all_var = MOI.get(model, MOI.NumberOfVariables())
        if int_var != all_var
            throw("Currently MibS works on only MIP-MIP problems and the input model is not MIP-MIP!!")
        end
    end

    lower_sense = MOI.get(lower, MOI.ObjectiveSense())
    return model, lower_variables, lower_objective, lower_constraints, lower_sense
end

function _index_to_row_link(model::MOI.FileFormats.MPS.Model)
    i = 0
    dict = Dict{MOI.ConstraintIndex, Int}()
    for (S, _) in MOI.FileFormats.MPS.SET_TYPES
        for ci in MOI.get(
            model,
            MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{Float64}, S}(),
        )
            dict[ci] = i
            i += 1
        end
    end
    return dict
end

function _index_to_column_link(model::MOI.FileFormats.MPS.Model)
    variables = MOI.get(model, MOI.ListOfVariableIndices())
    return Dict{MOI.VariableIndex,Int}(
        x => i - 1 for (i, x) in MOI.enumerate(variables)
    )
end

function _write_auxillary_file(
    new_model::MOI.FileFormats.MPS.Model,
    lower_variables::Vector{MOI.VariableIndex},
    lower_objective::MOI.ScalarAffineFunction,
    lower_constraints::Vector{MOI.ConstraintIndex},
    lower_sense::MOI.OptimizationSense,
    aux_filename::String
)
    rows = _index_to_row_link(new_model)
    cols = _index_to_column_link(new_model)
    obj_coefficients = Dict{MOI.VariableIndex,Float64}(
        x => 0.0 for x in lower_variables
    )
    for term in lower_objective.terms
        if haskey(obj_coefficients, term.variable_index)
            obj_coefficients[term.variable_index] += term.coefficient
        end
    end
    open(aux_filename, "w") do io
        println(io, "N $(length(lower_variables))")
        println(io, "M $(length(lower_constraints))")
        for x in lower_variables
            println(io, "LC $(cols[x])")
        end
        for y in lower_constraints
            println(io, "LR $(rows[y])")
        end
        for x in lower_variables
            println(io, "LO $(obj_coefficients[x])")
        end
        println(io, "OS ", lower_sense == MOI.MAX_SENSE ? -1 : 1)
    end
    return
end

function _call_mibs(mps_filename, aux_filename)
    io = IOBuffer()
    MibS_jll.mibs() do exe
        run(
            pipeline(
                `$(exe) -Alps_instance $(mps_filename) -MibS_auxiliaryInfoFile $(aux_filename)`,
                stdout = io,
            )
        )
    end
    seekstart(io)
    return read(io, String)
end

function _parse_output(
    output::String,
    new_model::MOI.FileFormats.MPS.Model,
    lower_variables::Vector{MOI.VariableIndex}
    )
    lines = split(output, '\n')
    found_status = false
    objective_value = NaN

    upper = Dict{Int,Float64}()
    lower = Dict{Int,Float64}()

    all_var = MOI.get(new_model, MOI.ListOfVariableIndices())

    CntU = 0
    CntD = 0

    Dict_Lower_Name = Dict()
    Dict_Lower_Value = Dict()
    Dict_Upper_Name = Dict()
    Dict_Upper_Value = Dict()
    Dict_Upper_IndexToModel = Dict()
    Dict_Lower_IndexToModel = Dict()
    Dict_All = Dict()

    for (x, y) in MOI.enumerate(all_var)
        nameofvar = MOI.get(new_model, MOI.VariableName(), y)
        if y in lower_variables
            Dict_Lower_Name[CntD] = nameofvar
            Dict_Lower_Value[nameofvar] = 0
            Dict_Lower_IndexToModel[CntD] = y
            CntD = CntD + 1
        else
            Dict_Upper_Name[CntU] = nameofvar
            Dict_Upper_Value[nameofvar] = 0
            Dict_Upper_IndexToModel[CntU] = y
            CntU = CntU + 1
        end
        Dict_All[y] = 0
    end

    for line in lines
        if !found_status
            if occursin("Optimal solution", line)
                found_status = true
            end
            continue
        end
        m = match(r"([xy])\[([0-9]+)\] \= (.+)", line)
        if m === nothing
            m = match(r"Cost \= (.+)", line)
            if m !== nothing
                objective_value = parse(Float64, m[1])
            end
            continue
        end

        column = parse(Int, m[2])
        value = parse(Float64, m[3])

        if m[1] == "x"
            upper[column] = value
            nameofvar = Dict_Upper_Name[column]
            Dict_Upper_Value[nameofvar] = value
            indexofvar = Dict_Upper_IndexToModel[column]
            Dict_All[indexofvar] = value
        else
            lower[column] = value
            nameofvar = Dict_Lower_Name[column]
            Dict_Lower_Value[nameofvar] = value
            indexofvar = Dict_Lower_IndexToModel[column]
            Dict_All[indexofvar] = value
        end
    end

    return (
        status = found_status,
        objective = objective_value,
        nonzero_upper = upper,
        nonzero_lower = lower,
        all_upper = Dict_Upper_Value,
        all_lower = Dict_Lower_Value,
        all_var = Dict_All
    )
end

"""
    solve_with_MibS(model::BilevelModel; silent::Bool = true)

## Inputs
* `model::BilevelModel`:  the model to optimize
* `silent::Bool = true`: controls the verbosity of the solver. If `silent`, nothing is printed. Set to `false` to display the MibS output.
## Outputs
This function returns a `NamedTuple` with fields:
* `status::Bool`: `true` if the problem is feasible and has an optimal solution. `false` otherwise.
* `objective::Float64`: objective value (cost) of the upper problem
* `nonzero_upper::Dict{Int, Float64}`: it returns `Dict{index => value}`, in which the `index` refers to the index of upper variables with non zero values and the index starts from `0`. Here, the order of the variables is based on their order of appearance in the MPS file.
* `nonzero_lower::Dict{Int, Float64}`: it has the same structure as `nonzero_upper`, but it represents the index of non-zero variables in the lower problem.
* `all_upper::Dict{String, Float64}`: it returns `Dict{name => value}` which contains all upper variables values (zero and non-zero). For recalling the variables, you need to use the same name as you used to define the variables, e.g., for `@variable(Upper(model), y, Int)`, we need to use `all_upper["y"]` to get the value of the variable `y`.
* `all_lower::Dict{String, Float64}`: it has the same structure as the `all_upper` but is defined for lower variables.
* `all_var::Dict{MOI.VariableIndex, Float64}`: it contains information on all variables (upper and lower) in the format of `MOI.VariableIndex` and their output values.

!!! warning
    Currently, `MibS` is designed to solve MIP-MIP problems only. Thus, if you define LP-MIP, MIP-LP, or LP-LP, it will throw an error.
"""
function solve_with_MibS(model::BilevelModel; silent::Bool = true)
    mktempdir() do path
        mps_filename = joinpath(path, "model.mps")
        aux_filename = joinpath(path, "model.aux")
        new_model, variables, objective, constraints, sense  =
            _build_single_model(model, true)
        MOI.write_to_file(new_model, mps_filename)
        _write_auxillary_file(
            new_model,
            variables,
            objective,
            constraints,
            sense,
            aux_filename,
        )
        output = _call_mibs(mps_filename, aux_filename)
        if !silent
            println(output)
        end
        return _parse_output(output, new_model, variables)
    end
end