text:
cdir = pwd(); ii = findlast("MDP", cdir); MDPpath = cdir[1:ii[end]] * "/";
include(MDPpath * "Installation/setup.jl");

# periodic boundary conditions
pbc = [1; 1; 1]

# path to the database
datapath = "/Users/ngoccuongnguyen/Dropbox (MIT)/Research/Software/MDP/fitting/snap/Ta/JSON/";
dataformat = "json"
fileextension = "json"
atomspecies = ["Ta"];

# folders in the datapath
folders = ["Displaced_A15", "Displaced_BCC", "Displaced_FCC", "Elastic_BCC",
            "Elastic_FCC", "GSF_110", "GSF_112", "Liquid", "Surface",
            "Volume_A15", "Volume_BCC", "Volume_FCC"];

# weights for energies, forces, and stresses in the linear fit
weightinner=[100            1               1.00E-08
            100             1               1.00E-08
            100             1               1.00E-08
            1.00E-08        1.00E-08        0.0001
            1.00E-09        1.00E-09        1.00E-09
            100             1               1.00E-08
            100             1               1.00E-08
            4.67E+02        1               1.00E-08
            100             1               1.00E-08
            1.00E+00        1.00E-09        1.00E-09
            1.00E+00        1.00E-09        1.00E-09
            1.00E+00        1.00E-09        1.00E-09]

# randomly selecting the configurations in the database
randomize = false;

# use all the data
percentage = 100.0;

# translate atom positions
translationvector = nothing

# rotate atom positions
rotationmatrix = nothing

# transpose lattice vectors because the lattice ordering is python-based
transposelattice = true

# training data
for i = 1:length(folders)
    traindata[i] = adddata(datapath * folders[i], dataformat, fileextension,
                percentage, randomize, atomspecies, weightinner[i,:], translationvector,
                rotationmatrix, transposelattice)
end

rcut = 4.67637   # cut-off radiuss

# single descriptors to catpure energy of isolated pure elements
descriptors[1] = PotentialStruct("single", "nonbonded", [1], "single", [0.0], rcut);
include(descriptors[1].potentialfunction * ".jl");  # include the descriptors file

# SNAP parameters
twojmax = 6
ntypes = 1
rcutfac = rcut;
rfac0 = 0.99363;
rmin0 = 0;
bzeroflag = 0;
quadraticflag = 0;
switchflag = 1;
chemflag = 0;
bnormflag = 0;
wselfallflag = 0;

snaparam = [ntypes, twojmax, rcutfac, rfac0, rmin0, bzeroflag, switchflag, quadraticflag, chemflag, bnormflag, wselfallflag]
elemradius = 0.5*ones(Float64, ntypes)
elemweight = ones(Float64, ntypes)

# SNAP descriptors
descriptors[2] = initsna(snaparam, elemradius, elemweight)

# ZBL reference potential
potentials[1] = addpotential("pair", "nonbonded", [1, 1], "zbl", [4.0, 4.8, 73], 4.8);
include(potentials[1].potentialfunction * ".jl");  # include the potential file

# loss function style must be: "energy", "force", "energyforce", "energyforcestress"
lossfunc = "energyforcestress"
# use least-square method
method = "lsq"
# divide energies by the number of atoms
normalizeenergy = true
# divide virials by the box volume
normalizestress = true
# optimization parameters
optim = setoptim(lossfunc, method, normalizeenergy, normalizestress)

# linear fit to compute SNAP coefficients
coeff,~ = linearfit2(traindata, descriptors, potentials, optim, pbc)

print("SNAP Coeffficients: "), show(stdout, "text/plain", coeff)

# compute unweighted MAE, RMSE, RSQ errors
energyerrors, forceerrors, stresserrors = validate2(traindata, descriptors, potentials, optim, coeff, pbc)

printerrors(folders, energyerrors, "Energy Errors")
printerrors(folders, forceerrors, "Force Errors")
printerrors(folders, stresserrors, "Stress Errors")

# ----------------------------------------------------------------------------------------
# Energy Errors |          MAE          |          RMSE          |          RSQ          |
# ----------------------------------------------------------------------------------------
# ALL           |   0.11265174691267    |    0.37977565433503    |    0.99419038677079   |
# Displaced_A15 |   0.00235792938939    |    0.00246855977428    |    -1.2156363735745   |
# Displaced_BCC |   0.00188371436433    |    0.00273064614363    |    0.99587522319726   |
# Displaced_FCC |   0.00061741229597    |    0.00075009168170    |    0.97160166434659   |
# Elastic_BCC   |   0.00599579098947    |    0.00604563796454    |    0.99230919870267   |
# Elastic_FCC   |   0.00354421105172    |    0.00445843637691    |    0.99333159772842   |
# GSF_110       |   0.00652425814607    |    0.00883023587747    |    0.77951696510583   |
# GSF_112       |   0.00868251914489    |    0.00961453139000    |    0.87318809800944   |
# Liquid        |   0.00566811969259    |    0.00655584081006    |    0.99601358415720   |
# Surface       |   0.01352360826034    |    0.01571242406414    |    0.98816795126889   |
# Volume_A15    |   0.24578877837149    |    0.34490495333075    |    0.99781261498874   |
# Volume_BCC    |   0.32857920069640    |    0.51271284190945    |    0.99853838108547   |
# Volume_FCC    |   0.81209055282782    |    1.18123504852115    |    0.98443646179057   |
# ----------------------------------------------------------------------------------------

# ----------------------------------------------------------------------------------------
# Force Errors  |          MAE          |          RMSE          |          RSQ          |
# ----------------------------------------------------------------------------------------
# ALL           |   0.08359239990607    |    0.16311584659554    |    0.98512903562289   |
# Displaced_A15 |   0.09806643270309    |    0.12820169392009    |    0.94957955563905   |
# Displaced_BCC |   0.14339589323270    |    0.18486360986945    |    0.98517811287921   |
# Displaced_FCC |   0.05740637316446    |    0.07166197131771    |    0.98061307246834   |
# Elastic_BCC   |   0.06439647995478    |    0.07899854689009    |    0.41298668776769   |
# Elastic_FCC   |   0.05065258662499    |    0.06490610225971    |    0.63719969979068   |
# GSF_110       |   0.02651057125534    |    0.05418215537551    |    0.99889872481223   |
# GSF_112       |   0.05698006876321    |    0.10025780716990    |    0.99737013744671   |
# Liquid        |   0.35379720579374    |    0.48364455518089    |    0.94779367837188   |
# Surface       |   0.04760235280733    |    0.10895743153042    |    0.99671811972583   |
# Volume_A15    |   2.79473115136733    |    6.25088426150716    |    0.69215396831169   |
# Volume_BCC    |   2.88326079591402    |    6.06857960431047    |    -4.7705594018581   |
# Volume_FCC    |   1.39818002339307    |    2.88854746833441    |    0.32360003124122   |
# ----------------------------------------------------------------------------------------

# ----------------------------------------------------------------------------------------
# Stress Errors |          MAE          |          RMSE          |          RSQ          |
# ----------------------------------------------------------------------------------------
# ALL           |   68331.2692226880    |    381756.218018856    |    0.93961426214254   |
# Displaced_A15 |   19553.8208834055    |    26783.3681374323    |    0.99755940235506   |
# Displaced_BCC |   14933.0136738189    |    20558.4972892687    |    0.99870935725821   |
# Displaced_FCC |   35323.5274494020    |    49684.5017070339    |    0.99113396279263   |
# Elastic_BCC   |   170.405533015481    |    268.553751189739    |    0.99999977908789   |
# Elastic_FCC   |   48132.1295201740    |    66784.2466913797    |    0.98419798844541   |
# GSF_110       |   1623.00572832491    |    2322.56048399292    |    0.99993839179953   |
# GSF_112       |   2068.90488220449    |    3116.39773656099    |    0.99987644651027   |
# Liquid        |   42920.1434033482    |    58562.0314191833    |    0.98628360736784   |
# Surface       |   2309.28255361132    |    4277.50368379242    |    0.99977822976790   |
# Volume_A15    |   132939.655632498    |    357432.429798387    |    0.95261560122466   |
# Volume_BCC    |   271466.858168366    |    770712.743079258    |    0.97424571846356   |
# Volume_FCC    |   304211.704271920    |    1.07923030094437    |    0.79235217127926   |
# ----------------------------------------------------------------------------------------

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
__precompile__()

module DiffEqPDEBase

using DiffEqBase

end # module

ext: jl


text:
using SuperResolution
using Test

@testset "SuperResolution" begin
    include("data.jl")
end


text:
@defcomp n2oforcing begin
    c_N2Oconcentration = Parameter(index=[time], unit="ppbv")
    c_CH4concentration = Parameter(index=[time], unit="ppbv")
    f0_N2Obaseforcing = Parameter(unit="W/m2", default=0.17)
    fslope_N2Oforcingslope = Parameter(unit="W/m2", default=0.12)
    c0_baseN2Oconc = Parameter(unit="ppbv", default=324.0)
    c0_baseCH4conc = Parameter(unit="ppbv", default=1803.0)
    f_N2Oforcing = Variable(index=[time], unit="W/m2")
    over_baseoverlap = Variable(unit="W/m2")

    function run_timestep(p, v, d, t)

        # from p.16 in Hope 2009
        if is_first(t)
            # calculate baseline forcing overlap in first time period
            v.over_baseoverlap = -0.47 * log(1 + 2.0e-5 * (p.c0_baseN2Oconc * p.c0_baseCH4conc)^0.75 + 5.3e-15 * p.c0_baseCH4conc * (p.c0_baseCH4conc * p.c0_baseN2Oconc)^1.52)
        end

        over = -0.47 * log(1 + 2.0e-5 * (p.c0_baseCH4conc * p.c_N2Oconcentration[t])^0.75 + 5.3e-15 * p.c0_baseCH4conc * (p.c0_baseCH4conc * p.c_N2Oconcentration[t])^1.52)
        v.f_N2Oforcing[t] = p.f0_N2Obaseforcing + p.fslope_N2Oforcingslope * (sqrt(p.c_N2Oconcentration[t]) - sqrt(p.c0_baseN2Oconc)) + over - v.over_baseoverlap
    end
end


text:
using StationXML, Test
using Dates: DateTime
import EzXML

@testset "Parsing" begin
    @testset "Dates" begin
        let
            # Time zone and too much precision
            @test StationXML.readstring("""
                <?xml version="1.0" encoding="UTF-8"?>
                <FDSNStationXML xmlns="http://www.fdsn.org/xml/station/1" schemaVersion="1.0">
                  <Source>StationXML.jl</Source>
                  <Sender>Test</Sender>
                  <Created>2019-04-01T09:58:21.123456789+01:00</Created>
                </FDSNStationXML>
                """).created == DateTime(2019, 04, 1, 10, 58, 21, 123)
            # Time zone ahead of UTC
            @test StationXML.readstring("""
                <?xml version="1.0" encoding="UTF-8"?>
                <FDSNStationXML xmlns="http://www.fdsn.org/xml/station/1" schemaVersion="1.0">
                  <Source>StationXML.jl</Source>
                  <Sender>Test</Sender>
                  <Created>2000-01-01T12:34:00.123456789-12:34</Created>
                </FDSNStationXML>
                """).created == DateTime(2000, 1, 1, 0, 0, 0, 123)
            # UTC specified
            @test StationXML.readstring("""
                <?xml version="1.0" encoding="UTF-8"?>
                <FDSNStationXML xmlns="http://www.fdsn.org/xml/station/1" schemaVersion="1.0">
                  <Source>StationXML.jl</Source>
                  <Sender>Test</Sender>
                  <Created>2019-04-01T09:58:21Z</Created>
                </FDSNStationXML>
                """).created == DateTime(2019, 04, 1, 9, 58, 21)
            # Invalid format throws a helpful error
            for dtstring in ("2000", "2000-01", "2000-01-01", "2000-01-01T",
                    "2000-01-01T00", "2000-01-01T00:00")
                node = EzXML.parsexml("<element>$dtstring</element>").root
                @test_throws ArgumentError StationXML.parse_node(DateTime, node)
            end
            # Dates in attributes
            @test StationXML.readstring("""
                <?xml version='1.0' encoding='UTF-8'?>
                <FDSNStationXML xmlns="http://www.fdsn.org/xml/station/1" schemaVersion="1.1">
                  <Source>SeisComP3</Source>
                  <Sender>ODC</Sender>
                  <Module/>
                  <ModuleURI/>
                  <Created>2020-04-28T21:17:08.868852Z</Created>
                  <Network code="NL" startDate="1993-01-01T00:00:00.987654Z" restrictedStatus="open">
                    <Description>Netherlands Seismic and Acoustic Network</Description>
                  </Network>
                </FDSNStationXML>
                """).network[1].start_date == DateTime(1993, 1, 1, 0, 0, 0, 987)
        end
    end

    @testset "Required fields" begin
        @testset "Missing source" begin
            @test_throws UndefVarError StationXML.readstring("""
                <?xml version="1.0" encoding="UTF-8"?>
                <FDSNStationXML xmlns="http://www.fdsn.org/xml/station/1" schemaVersion="1.1">
                  <Created>2000-01-01T12:34:00.123456789-12:34</Created>
                </FDSNStationXML>
                """)
        end
        @testset "Missing network code" begin
            @test_throws UndefVarError StationXML.readstring("""
                <?xml version="1.0" encoding="UTF-8"?>
                <FDSNStationXML xmlns="http://www.fdsn.org/xml/station/1" schemaVersion="1.1">
                  <Created>2000-01-01T12:34:00.123456789-12:34</Created>
                  <Source>Test</Source>
                  <Network>
                  </Network>
                </FDSNStationXML>
                """)
        end
        @testset "Missing station latitude" begin
            @test_throws UndefVarError StationXML.readstring("""
                <?xml version="1.0" encoding="UTF-8"?>
                <FDSNStationXML xmlns="http://www.fdsn.org/xml/station/1" schemaVersion="1.1">
                  <Created>2000-01-01T12:34:00.123456789-12:34</Created>
                  <Source>Test</Source>
                  <Network code="XX">
                    <Station code="YYY">
                      <Longitude>0</Longitude>
                      <Elevation>0</Elevation>
                      <Site>
                        <Name>Example site</Name>
                      </Site>
                    </Station>
                  </Network>
                </FDSNStationXML>
                """)
        end
    end

    # Changes from v1.0 to v1.1
    @testset "Version 1.1" begin
        let sxml = StationXML.readstring("""
                <FDSNStationXML xmlns="http://www.fdsn.org/xml/station/1"
                  xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
                  xmlns:iris="http://www.fdsn.org/xml/station/1/iris"
                  xsi:schemaLocation="http://www.fdsn.org/xml/station/1 http://www.fdsn.org/xml/station/fdsn-station-1.1.xsd"
                  schemaVersion="1.1">
                 <Source>IRIS-DMC</Source>
                 <Sender>IRIS-DMC</Sender>
                 <Module>IRIS WEB SERVICE: fdsnws-station | version: 1.1.45</Module>
                 <ModuleURI>http://service.iris.edu/fdsnws/station/1/query?nodata=204</ModuleURI>
                 <Created>2020-05-05T19:13:31</Created>
                 <Network code="1B" startDate="2014-01-01T00:00:00" endDate="2014-12-31T23:59:59"
                   restrictedStatus="open" sourceID="mailto:net@example.com">
                  <Description>Sweetwater Array (1B)</Description>
                  <TotalNumberStations>2268</TotalNumberStations>
                  <Identifier type="some type">Identifying text</Identifier>
                  <Identifier type="another type">More text</Identifier>
                  <SelectedNumberStations>2268</SelectedNumberStations>
                  <Station code="5R536" startDate="2014-01-01T00:00:00"
                    endDate="2014-12-12T23:59:59" restrictedStatus="open"
                    sourceID="https://example.com"
                    iris:alternateNetworkCodes=".UNRESTRICTED">
                   <DataAvailability>
                    <Extent start="2000-01-01T00:00:00.000999Z" end="3000-01-01T00:00:00">
                    </Extent>
                    <Span start="1990-01-01T00:00:00" end="1991-01-01T00:00:00" numberSegments="2"
                     maximumTimeTear="0.1"/>
                   </DataAvailability>
                   <Identifier>ID</Identifier>
                   <Latitude>32.748402</Latitude>
                   <Longitude>-100.535698</Longitude>
                   <Elevation>634.7</Elevation>
                   <Site>
                    <Name>207536</Name>
                   </Site>
                   <TotalNumberChannels>1</TotalNumberChannels>
                   <SelectedNumberChannels>1</SelectedNumberChannels>
                   <Channel code="SXZ" locationCode="00">
                    <Identifier type="A">B</Identifier>
                    <Latitude>1</Latitude>
                    <Longitude>2</Longitude>
                    <Elevation>3</Elevation>
                    <Depth>4</Depth>
                   </Channel>
                  </Station>
                 </Network>
                </FDSNStationXML>
                """)
            @test sxml.schema_version == "1.1"
            @testset "Add Identifier" begin
                @test length(sxml.network[1].identifier) == 2
                @test sxml.network[1].identifier[1].value == "Identifying text"
                @test sxml.network[1].identifier[1].type == "some type"
                @test sxml.network[1].identifier[2].value == "More text"
                @test sxml.network[1].identifier[2].type == "another type"
                @test length(sxml.network[1].station[1].identifier) == 1
                @test sxml.network[1].station[1].identifier[1].value == "ID"
                @test sxml.network[1].station[1].identifier[1].type === missing
                @test sxml.network[1].station[1].channel[1].identifier[1].value == "B"
                @test sxml.network[1].station[1].channel[1].identifier[1].type == "A"
            end
            @testset "Optional CreationDate" begin
                @test sxml.network[1].station[1].creation_date === missing
            end
            @testset "Add sourceID" begin
                @test sxml.network[1].source_id == "mailto:net@example.com"
                @test sxml.network[1].station[1].source_id == "https://example.com"
                @test sxml.network[1].station[1].channel[1].source_id === missing
            end
            @testset "Add DataAvailability" begin
                @test sxml.network[1].station[1].data_availability.extent.start == DateTime(2000)
                @test sxml.network[1].station[1].data_availability.extent.end_ == DateTime(3000)
                @test sxml.network[1].station[1].data_availability.span[1].start == DateTime(1990)
                @test sxml.network[1].station[1].data_availability.span[1].end_ == DateTime(1991)
                @test sxml.network[1].station[1].data_availability.span[1].number_segments == 2
                @test sxml.network[1].station[1].data_availability.span[1].maximum_time_tear == 0.1
            end
        end
    end

    @testset "Pressure channel" begin
        let node = EzXML.root(EzXML.parsexml("""
                <Channel code="LDO" locationCode="00" startDate="2008-10-10T00:00:00" restrictedStatus="open">
                 <Latitude>10.00207</Latitude>
                 <Longitude>-84.111389</Longitude>
                 <Elevation>1186</Elevation>
                 <Depth>0</Depth>
                 <Azimuth>0</Azimuth>
                 <Dip>0</Dip>
                 <Type>CONTINUOUS</Type>
                 <Type>WEATHER</Type>
                 <SampleRate>1E00</SampleRate>
                 <ClockDrift>1E-04</ClockDrift>
                 <Sensor>
                  <Description>Microbaro VAISALA</Description>
                 </Sensor>
                 <Response>
                 <InstrumentSensitivity>
                   <InputUnits>
                     <Name>PA</Name>
                     <Description>PRESSURE in Pascals</Description>
                   </InputUnits>
                   <OutputUnits>
                     <Name>COUNTS</Name>
                     <Description>DIGITAL UNIT in Counts</Description>
                   </OutputUnits>
                 </InstrumentSensitivity>
                 </Response>
                </Channel>
                """))
            channel = StationXML.parse_node(StationXML.Channel, node)
            @test channel.response.instrument_sensitivity.value === missing
            @test channel.response.instrument_sensitivity.input_units.name == "PA"
        end
    end
end


text:
using Ghost
using Ghost: Tape, Input, mkcall, Variable, V

struct ONNXCtx
    name2var::Dict{String, Variable}
    backends::Vector{Symbol}
    exec::Bool
end

ONNXCtx(backends; exec=true) = ONNXCtx(Dict(), backends, exec)
ONNXCtx(;exec=true) = ONNXCtx(Dict(), [:ONNX], exec)

# TODO: implement rebind_context!()

"""
    getindex(tape::Tape{ONNXCtx}, onnx_name::String)

Get operation on the tape using the name in ONNX graph
"""
Base.getindex(tape::Tape{ONNXCtx}, onnx_name::String) =
    tape[tape.c.name2var[onnx_name]]

###############################################################################
#                               Operations                                    #
###############################################################################

mrev(x) = x
mrev(x::AbstractVector) = reverse(x)
prev(x) = x
prev(x::AbstractVector) = reshape(permutedims(reverse(reshape(x, length(x) ÷ 2,:);dims=1)),:)

# mrev = maybe reverse. prev = rearrange padding, e.g. (1,2,1,2) => (2,2,1,1) or (1,2,3,1,2,3) => (3,3,2,2,1,1)
_akpsd(params) = get(params, :activation, identity), mrev(get(params, :kernel_shape, 1)), prev(get(params, :pads, 0)), mrev(get(params, :strides, 1)), mrev(get(params, :dilations, 1))
akpsd(params) = a2t.(_akpsd(params))
a2t(x) = x
a2t(a::AbstractArray) = Tuple(a)

conv_attr_onnx2tape(attrs) = Dict(
    :stride => mrev(get(attrs, :strides, 1)),
    :pad => prev(get(attrs, :pads, 0)),
    :dilation => mrev(get(attrs, :dilations, 1)),
    :groups => get(attrs, :group, 1),
    # kenrnel_shape => mrev(get(params, :kernel_shape, 1)) -- not used in NNlib.conv
)

"""
    push_call!(tape::Tape{ONNXCtx}, fn, args...; kwargs)

Shortcut for `push!(tape, mkcall(fn, args..))` also handling
keyword arguments and respecting `ONNXCtx.exec` setting.
"""
function push_call!(tape::Tape{ONNXCtx}, fn, args...; kwargs...)
    kwargs = NamedTuple(kwargs)
    if !isempty(kwargs)
        args = (kwargs, fn, args...)
        fn = Core.kwfunc(fn)
    end
    op = tape.c.exec ? mkcall(fn, args...) : mkcall(fn, args...; val=nothing)
    return push!(tape, op)
end

# A few constants to keep function signatures concise
struct OpConfig{BE, Op} end
const VarVec = Vector{Ghost.Variable}
const AttrDict = Dict{Symbol, Any}

function load_node!(tape::Tape, nd::NodeProto, backend::Symbol)
    args = [tape.c.name2var[name] for name in nd.input]
    attrs = convert(Dict{Symbol, Any}, Dict(nd.attribute))
    conf = OpConfig{backend, Symbol(nd.op_type)}()
    out = load_node!(tape, conf, args, attrs)
    if out isa Tuple
        for i=1:length(nd.output)
            tape.c.name2var[nd.output[i]] = out[i]
        end
    else
        tape.c.name2var[nd.output[1]] = out
    end
end

load_node!(::Tape, ::OpConfig{BE, Op}, ::VarVec, ::AttrDict) where {BE, Op} = missing

# function load_node!(tape::Tape, ::OpConfig{:ONNX, :Conv}, args::VarVec, attrs::AttrDict)
#     _,_,p,s,d = akpsd(attrs)
#     kw = (stride = s, pad = p, dilation = d, groups = get(attrs, "group", 1))
#     return push_call!(tape, conv, args...; kw...)
# end

# function load_node!(tape::Tape, ::OpConfig{:ONNX, :Gemm}, args::VarVec, attrs::AttrDict)
#     kw = Dict(
#         :tA => get(attrs, :transA, 0),
#         :tB => get(attrs, :transB, 0),
#         :α => get(attrs, :alpha, 1),
#         :β => get(attrs, :beta, 0)
#     )
#     return push_call!(tape, onnx_gemm, args...; kw...)
# end

# function load_node!(tape::Tape, ::OpConfig{:ONNX, :Flatten}, args::VarVec, attrs::AttrDict)
#     return push_call!(tape, onnx_flatten, args...; attrs...)
# end

function load_node!(tape::Tape, ::OpConfig{:ONNX, :Add}, args::VarVec, attrs::AttrDict)
    return push_call!(tape, add, args...)
end

# function load_node!(tape::Tape, ::OpConfig{:ONNX, :Mul}, args::VarVec, attrs::AttrDict)
#     return push_call!(tape, mul, args...)
# end

# function load_node!(tape::Tape, ::OpConfig{:ONNX, :Relu}, args::VarVec, attrs::AttrDict)
#     return push_call!(tape, relu, args[1])
# end

# function load_node!(tape::Tape, ::OpConfig{:ONNX, :MaxPool}, args::VarVec, attrs::AttrDict)
#     _,k,p,s,_ = akpsd(attrs)
#     return push_call!(tape, maxpool, args[1], k; pad=p, stride=s)
# end

# function load_node!(tape::Tape, ::OpConfig{:ONNX, :BatchNormalization}, args::VarVec, attrs::AttrDict)
#     ϵ = get(attrs, :epsilon, 1f-5)
#     momentum = get(attrs, :momentum, 9f-1)
#     training_mode = Bool(get(attrs, :training_mode, 0))
#     res = push_call!(tape, batch_norm, args..., ϵ, momentum, training_mode)
#     if training_mode
#         y = push_call!(tape, getfield, 1)
#         μ_new = push_call!(tape, getfield, 2)
#         σ²_new = push_call!(tape, getfield, 3)
#         return y, μ_new, σ²_new
#     else
#         return res
#     end
# end

# function load_node!(tape::Tape, ::OpConfig{:ONNX, :GlobalAveragePool}, args::VarVec, attrs::AttrDict)
#     return push_call!(tape, global_average_pool, args...)
# end

###############################################################################
#                                    API                                      #
###############################################################################

"""
    load(io::IO, model_args...; backends=[:ONNX], exec::Bool=true)
    load(filename::String, model_args...; backends=[:ONNX], exec::Bool=true)

Load an ONNX model as a Ghost.Tape. The way a particular ONNX node is deserialized is
controlled by methods of [load_node!](@ref) dispatched by backend and node's op_type.

`backends` can be used to customize the (in-order) choices for which backend libraries to use.

`exec` parameter instructs the loader to execute every operation pushed onto the tape just after
the push, making the debugging easier. Default is `true`.

See also: [save!](@ref)
"""
function load(io::IO, model_args...; backends=[:ONNX], exec::Bool=true)
    onnx_model = readproto(io, ModelProto());
    g = onnx_model.graph;
    tape = Tape(ONNXCtx(backends; exec=exec))
    # create map of initializers
    init_vals = Dict{String, Any}()
    for init in g.initializer
        # TODO: consider non-array inputs
        init_vals[init.name] = array(init)
    end
    # load inputs
    arg_idx = 1
    for inp in g.input
        val = get(init_vals, inp.name, missing)
        if val === missing && exec == true
            val = model_args[arg_idx]
            arg_idx += 1
        end
        v = push!(tape, Input(val))
        tape.c.name2var[inp.name] = v
    end
    # load nodes
    for nd in g.node
        success = false
        for backend in tape.c.backends
          if !ismissing(load_node!(tape, nd, backend))
            success = true
            @debug "Loaded $(nd.op_type) using backend $(backend)"
            break
          end
        end
        success || error("Couldn't load node for $(nd.op_type), " *
                         "tried the following backends: $(tape.c.backends)")
    end
    tape.result = Ghost.bound(tape, V(length(tape)))
    return tape
end

function load(filename::String, model_args...; backends=[:ONNX], exec::Bool=true)
    return open(filename) do io
        load(io, model_args...; backends=backends, exec=exec)
    end
end


text:
@doc_str """
    eigenvector_centrality(g)

Compute the eigenvector centrality for the graph `g`.

Eigenvector centrality computes the centrality for a node based on the
centrality of its neighbors. The eigenvector centrality for node `i` is

    ``\mathbf{Ax} = \lambda \mathbf{x}``

where `A` is the adjacency matrix of the graph `g` with eigenvalue `\lambda`.
By virtue of the Perron–Frobenius theorem, there is a unique and positive
solution if `\lambda` is the largest eigenvalue associated with the
eigenvector of the adjacency matrix `A`.

### References

- Phillip Bonacich: Power and Centrality: A Family of Measures.
    American Journal of Sociology 92(5):1170–1182, 1986
    http://www.leonidzhukov.net/hse/2014/socialnetworks/papers/Bonacich-Centrality.pdf
- Mark E. J. Newman: Networks: An Introduction.
       Oxford University Press, USA, 2010, pp. 169.
"""
eigenvector_centrality(g::AbstractGraph) = abs.(vec(eigs(adjacency_matrix(g), nev=1)[2]))::Vector{Float64}

max_stars_repo_name: AsileBcd/JohnsonCycles

max_issues_repo_name: AsileBcd/JohnsonCycles

max_forks_repo_name: AsileBcd/JohnsonCycles


text:
export outertangents, scale, scalearea

function Circle(p1::Point, p2::Point, p3::Point)
    circlethrough(p1, p2, p3)
end

function circlethrough(p1::Point, p2::Point, p3::Point)
    l1 = Line(p1, p2)
    c1 = fraction(l1, 0.5)
    l2 = Line(p2, p3)
    c2 = fraction(l2, 0.5)

    dl1 = vector(l1)
    dl2 = vector(l2)
    cl1 = Line(c1, c1 + Point(dl1.y, -dl1.x))
    cl2 = Line(c2, c2 + Point(dl2.y, -dl2.x))

    center = intersection(cl1, cl2)
    radius = from_to(p1, center) |> magnitude
    Circle(center, radius)
end

function Circle(center::Point, p1::Point)
    radius = from_to(center, p1) |> magnitude
    Circle(center, radius)
end

needed_attributes(::Type{Circle}) = (Visible, Fill, Stroke, Linewidth, Linestyle)

function intersection(c::Circle, l::Line)
    # algorithm works for circle at (0, 0)
    frm_moved = l.from - c.center
    to_moved = l.to - c.center
    dx, dy = from_to(frm_moved, to_moved).xy
    dr = sqrt(dx ^ 2 + dy ^ 2)
    D = cross(frm_moved, to_moved)

    delta = (c.radius ^ 2) * (dr ^ 2) - (D ^ 2)

    if abs(delta) < 1e-10
        # tangent
        x = D * dy / (dr ^ 2)
        y = -D * dx / (dr ^ 2)
        return Point(x, y) + c.center
    elseif delta < 0
        return nothing
    else
        xplusminus = sign(dy) * dx * sqrt(delta)
        yplusminus = abs(dy) * sqrt(delta)

        x1 = (D * dy + xplusminus) / (dr ^ 2)
        x2 = (D * dy - xplusminus) / (dr ^ 2)
        y1 = (-D * dx + yplusminus) / (dr ^ 2)
        y2 = (-D * dx - yplusminus) / (dr ^ 2)
        return (Point(x1, y1) + c.center, Point(x2, y2) + c.center)
    end
end

function tangentpoints(c::Circle, through::Point)
    p_to_center = from_to(through, c.center)
    a = asin(c.radius / magnitude(p_to_center))
    b = atan(p_to_center.y, p_to_center.x)
    t1 = b - a
    p_tangent_1 = Point(sin(t1), -cos(t1)) * c.radius + c.center
    t2 = b + a
    p_tangent_2 = Point(-sin(t2), cos(t2)) * c.radius + c.center
    return p_tangent_1, p_tangent_2
end

function tangents(c::Circle, through::Point)
    p_tangent_1, p_tangent_2 = tangentpoints(c, through)
    return Line(through, p_tangent_1), Line(through, p_tangent_2)
end

function outertangents(c1::Circle, c2::Circle)
    big_circle, small_circle = c1.radius > c2.radius ? (c1, c2) : (c2, c1)
    radius_difference = big_circle.radius - small_circle.radius
    gamma = -atan(big_circle.center.y - small_circle.center.y, big_circle.center.x - small_circle.center.x)

    beta = asin(radius_difference / magnitude(from_to(small_circle.center, big_circle.center)))
    alpha1 = gamma - beta
    alpha2 = gamma + beta

    small_tangent_point_1 = small_circle.center + Point(
        small_circle.radius * cos(pi / 2 - alpha1),
        small_circle.radius * sin(pi / 2 - alpha1),
    )

    big_tangent_point_1 = big_circle.center + Point(
        big_circle.radius * cos(pi / 2 - alpha1),
        big_circle.radius * sin(pi / 2 - alpha1),
    )

    small_tangent_point_2 = small_circle.center + Point(
        small_circle.radius * cos(-pi / 2 - alpha2),
        small_circle.radius * sin(-pi / 2 - alpha2),
    )

    big_tangent_point_2 = big_circle.center + Point(
        big_circle.radius * cos(-pi / 2 - alpha2),
        big_circle.radius * sin(-pi / 2 - alpha2),
    )

    return [Line(small_tangent_point_1, big_tangent_point_1), Line(small_tangent_point_2, big_tangent_point_2)]
end

function Point(c::Circle, angle::Angle)
    c.center + Point(cos(angle), sin(angle)) * c.radius
end

function Point(c::Circle, angle::Angle, fraction::Real)
    c.center + Point(cos(angle), sin(angle)) * c.radius * fraction
end

function closestto(c::Circle, p::Point)
    normalize(from_to(c.center, p)) * c.radius
end

area(c::Circle) = pi * (c.radius ^ 2)
circumference(c::Circle) = 2pi * r
scale(c::Circle, s::Real) = Circle(c.center, c.radius * s)
scalearea(c::Circle, s::Real) = Circle(c.center, c.radius * sqrt(s))
grow(c::Circle, radius::Real) = Circle(c.center, c.radius + radius)


text:
# This file is part of the IntervalArithmetic.jl package; MIT licensed

"""
    in(x, a)
    ∈(x, a)

Checks if the number `x` is a member of the interval `a`, treated as a set.
Corresponds to `isMember` in the ITF-1788 Standard.
"""
function in(x::T, a::Interval) where T<:Real
    isinf(x) && return false
    a.lo <= x <= a.hi
end

"""
    issubset(a,b)
    ⊆(a,b)

Checks if all the points of the interval `a` are within the interval `b`.
"""
function ⊆(a::Interval, b::Interval)
    isempty(a) && return true
    b.lo ≤ a.lo && a.hi ≤ b.hi
end

function ⊂(a::Interval, b::Interval)
    a == b && return false
    a ⊆ b
end

function ⊇(a::Interval, b::Interval)
    b ⊆ a
end

function ⊃(a::Interval, b::Interval)
    b ⊂ a
end

# isinterior
function isinterior(a::Interval, b::Interval)
    isempty(a) && return true
    islessprime(b.lo, a.lo) && islessprime(a.hi, b.hi)
end
const ⪽ = isinterior  # \subsetdot

# Disjoint:
function isdisjoint(a::Interval, b::Interval)
    (isempty(a) || isempty(b)) && return true
    islessprime(b.hi, a.lo) || islessprime(a.hi, b.lo)
end

# Intersection
"""
    intersect(a, b)
    ∩(a,b)

Returns the intersection of the intervals `a` and `b`, considered as
(extended) sets of real numbers. That is, the set that contains
the points common in `a` and `b`.
"""
function intersect(a::Interval{T}, b::Interval{T}) where T
    isdisjoint(a,b) && return emptyinterval(T)

    Interval(max(a.lo, b.lo), min(a.hi, b.hi))
end
# Specific promotion rule for intersect:
intersect(a::Interval{T}, b::Interval{S}) where {T,S} =
    intersect(promote(a, b)...)

## Hull
"""
    hull(a, b)

Returns the "interval hull" of the intervals `a` and `b`, considered as
(extended) sets of real numbers, i.e. the smallest interval that contains
all of `a` and `b`.
"""
hull(a::Interval, b::Interval) = Interval(min(a.lo, b.lo), max(a.hi, b.hi))
#
# hull{T,S}(a::Interval{T}, b::Interval{S}) = hull(promote(a, b)...)

"""
    union(a, b)
    ∪(a,b)

Returns the union (convex hull) of the intervals `a` and `b`; it is equivalent
to `hull(a,b)`.
"""
union(a::Interval, b::Interval) = hull(a, b)
#
# union(a::Interval, b::Interval) = union(promote(a, b)...)

"""
    setdiff(x::Interval, y::Interval)

Calculate the set difference `x ∖ y`, i.e. the set of values
that are inside the interval `x` but not inside `y`.

Returns an array of intervals.
The array may:

- be empty if `x ⊆ y`;
- contain a single interval, if `y` overlaps `x`
- contain two intervals, if `y` is strictly contained within `x`.
"""
function setdiff(x::Interval, y::Interval)
    intersection = x ∩ y

    isempty(intersection) && return [x]
    intersection == x && return typeof(x)[]  # x is subset of y; setdiff is empty

    x.lo == intersection.lo && return [Interval(intersection.hi, x.hi)]
    x.hi == intersection.hi && return [Interval(x.lo, intersection.lo)]

    return [Interval(x.lo, y.lo),
            Interval(y.hi, x.hi)]

end


text:
module Forall

typealias Callable Union{Function, DataType}

export forall

"""
    forall(f, iter)

Iterate over `iter`, applying `f` to each element. The call

    forall(iter) do x
        f(x)
    end

is equivalent to

    for x in iter
        f(x)
    end

This also corresponds to `(map(f, iter); nothing)`.
"""
function forall(f::Callable, iter)
    for x in iter
        f(x)
    end
end

end


text:
import CoverageTools

const isevaldef = CoverageTools.isevaldef
const isfuncexpr = CoverageTools.isfuncexpr
const function_body_lines = CoverageTools.function_body_lines
const function_body_lines! = CoverageTools.function_body_lines!


text:
module Molecules

using Unitful
using PeriodicTable
using LinearAlgebra
using Formatting

import PhysicalConstants.CODATA2018: a_0

export Atom

bohr_to_angstrom = convert(Float64, a_0 / 1u"Å")

struct Atom{I,F}
    Z::I
    mass::F
    xyz::Vector{F}
end

# For checking if two atoms are the same
function Base.:(==)(A1::Atom, A2::Atom)
    A1.Z    == A2.Z    &&
    A1.mass == A2.mass &&
    A1.xyz  == A2.xyz
end

# Alias for a vector of atoms
Molecule = Vector{A} where A <: Atom

include("Parse.jl")
include("Transformations.jl")
#include("ZMAT.jl")
include("Properties.jl")
# include("RotationalEnergy.jl")
include("PointGroup/PointGroup.jl")
# include("Symmetry/Symmetry.jl")
end # module


text:
@inline fieldindex(::T, field::Val) where {T} = fieldindex(T, field)
@generated function fieldindex(::Type{T}, ::Val{field}) where {T,field}
    field isa Symbol || return :(error("`field` must be a symbol; given: $field"))
    return findfirst(n -> n === field, fieldnames(T))
end

@inline function fieldpointer(obj, field::Val)
    offset = fieldoffset(typeof(obj), fieldindex(obj, field))
    return Ptr{UInt}(pointer_from_objref(obj) + offset)
end

@inline atomic_modifyfield!(obj, field::Val, op, x) =
    atomic_modifyfield!(obj, field, op, x, seq_cst)
@inline function atomic_modifyfield!(obj, field::Val, op::OP, x, order) where {OP}
    FieldType = fieldtype(typeof(obj), fieldindex(obj, field))
    fptr = Ptr{FieldType}(fieldpointer(obj, field))
    v = convert(FieldType, x)
    GC.@preserve obj begin
        old = UnsafeAtomics.modify!(fptr, op, v, order)
    end
    return old
end


text:
using Documenter, Pretalx

makedocs(;
    modules=[Pretalx],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/oxinabox/Pretalx.jl/blob/{commit}{path}#L{line}",
    sitename="Pretalx.jl",
    authors="Lyndon White",
    assets=String[],
)

deploydocs(;
    repo="github.com/oxinabox/Pretalx.jl",
)


text:
#Code written by Skipper af Rosenborg and Esben Bager
##### MEAN-VARIANCE OPTIMIZATION
#This is not supposed to produce great results, since mean and covariance
#is supposed to be very exact before good weights are found

#Load packages and support functions
using JuMP
using Gurobi
using StatsBase
using DataFrames
using CSV
include("SupportFunction.jl")
include("DataLoad.jl")
include("CEOSupportFunctions.jl")
println("Leeeeroooy Jenkins")
#Esben's path
#path = "$(homedir())/Documents/GitHub/Thesis/Data/IndexData"

#Set path for dataloading
#path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Thesis/Data/IndexDataDiff/"
path = "/zhome/9f/d/88706/SpecialeCode/Thesis/Data/IndexDataDiff/"

#Set parameters
trainingSize = 240
possibilities = 5
industries = ["NoDur", "Durbl", "Manuf", "Enrgy", "HiTec", "Telcm", "Shops", "Hlth", "Utils", "Other"]
industriesTotal = length(industries)

#Define the input and transformations
modelMatrix = zeros(industriesTotal, possibilities)
noDurModel = [1 0 0 1 1]
modelMatrix[1, :] = noDurModel
for i=2:industriesTotal
    modelMatrix[i, :] = noDurModel
end

#Create and fill X and Y arrays
XArrays = Array{Array{Float64, 2}}(industriesTotal)
YArrays = Array{Array{Float64, 2}}(industriesTotal)
XArrays, YArrays = generateXandYs(industries, modelMatrix)

#Set
nGammas = 5
standY = YArrays[1]
lambdaValues = [0.1, 0.5, 1, 2, 10]
nRows = size(standY)[1]
amountOfModels = nGammas^3

#Initialization of parameters
w1N = repeat([1/10], outer = 10) #1/N weights
gamma = 2.4 #risk aversion
validationPeriod = 0 # Number of periods to do parameter training
PMatrix = zeros(nRows-trainingSize-2, amountOfModels)

bestModelAmount = nGammas^3
bestModelConfigs = zeros(bestModelAmount, 4)
bestModelIndexes = zeros(bestModelAmount)

return1NMatrix = zeros(nRows-trainingSize)
returnCEOMatrix = zeros(nRows-trainingSize, bestModelAmount)
returnPerfectMatrix = zeros(nRows-trainingSize)

weightsPerfect = zeros(nRows-trainingSize, 10)
weightsCEO     = zeros(nRows-trainingSize, 10, bestModelAmount)
expectedReturnMatrix = zeros(nRows-trainingSize, 10)
forecastErrors = zeros(nRows-trainingSize, 10)

#path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results/CEO/RFR/VIXTimeTA2.4/"
path = "/zhome/9f/d/88706/SpecialeCode/Results/CEORFRSearch/"
weightsPerfect = Array{Float64}(CSV.read(path*"weightsPerfect.csv",header=false, datarow=1, nullable=false))
returnPerfectMatrix = Array{Float64}(CSV.read(path*"returnPerfectMatrix.csv", header=false, datarow=1, nullable=false))

counter = 1
for l2 = 1:nGammas
    for l3 = 1:nGammas
        for l4 = 1:nGammas
            bestModelConfigs[counter,:] = [1 lambdaValues[l2] lambdaValues[l3] lambdaValues[l4]]
            counter += 1
        end
    end
end

futureLambdaValues = zeros(3,5)
futureLambdaValues[1,:] = lambdaValues
futureLambdaValues[2,:] = lambdaValues
futureLambdaValues[3,:] = lambdaValues

PMatrix = zeros(nRows-trainingSize-2, amountOfModels)
println("Lambda 2 ", futureLambdaValues[1,:])
println("Lambda 3 ", futureLambdaValues[2,:])
println("Lambda 4 ", futureLambdaValues[3,:])
for j = 1:parse(Int64,"4")#ARGS[2])
    PMatrix = zeros(nRows-trainingSize-2, amountOfModels)
    println(j)
    for i = 1:50
        #path = "/zhome/9f/d/88706/SpecialeCode/Results/CEORFRSearch/"
        path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results/CEO/CEOSearch/"
        PMatrixTemp = Array{Float64}(CSV.read(path*"Iter"*string(j-1)*"/"*string(i-1)*"_returnPvalueOutcome1to200.csv",header=false, datarow=1, nullable=false, types=fill(Float64,252)))

        PMatrix[i,:] += PMatrixTemp[i,3+125:end]
    end
    bestIndex = indmax(mean(PMatrix[1:50,:],1))
    bestConfig = bestModelConfigs[bestIndex,:]

    for k = 1:3
        bestIndx = find(x -> bestConfig[k+1] == x,futureLambdaValues[k,:])[1]

        best = futureLambdaValues[k,bestIndx]

        if bestIndx == 1
            min = futureLambdaValues[k,1]*0
            max = futureLambdaValues[k,2]
        elseif bestIndx == 5
            min = futureLambdaValues[k,4]
            max = futureLambdaValues[k,5]*2
        else
            min = futureLambdaValues[k,bestIndx-1]
            max = futureLambdaValues[k,bestIndx+1]
        end

        futureLambdaValues[k,1] = best-(best-min)*0.8
        futureLambdaValues[k,2] = best-(best-min)*0.2
        futureLambdaValues[k,3] = best
        futureLambdaValues[k,4] = best+(max-best)*0.2
        futureLambdaValues[k,5] = best+(max-best)*0.8
    end
    println("Lambda 2 ", futureLambdaValues[1,:])
    println("Lambda 3 ", futureLambdaValues[2,:])
    println("Lambda 4 ", futureLambdaValues[3,:])

    counter = 1
    for l2 = 1:nGammas
        for l3 = 1:nGammas
            for l4 = 1:nGammas
                bestModelConfigs[counter,:] = [1 futureLambdaValues[1,l2] futureLambdaValues[2,l3] futureLambdaValues[3,l4]]
                counter += 1
            end
        end
    end
end

println(bestModelConfigs[1:10,:])

for i = 1:bestModelAmount
	bestModelIndexes[i]=i
end

PMatrix = zeros(nRows-trainingSize-2, amountOfModels)
inputArg1 = 0 #83 breaks the trainingSize
inputArg1 = parse(Int64,ARGS[1])
println("Starting CEO Validation loop")
t=1+inputArg1
println("Time $t/50")
trainingXArrays, trainingYArrays, validationXRows, validationY, OOSXArrays, OOSYArrays, OOSRow, OOSY = createDataSplits(XArrays, YArrays, t, trainingSize)

for m = 1:bestModelAmount
    expectedReturns = zeros(10)

    betaArray, U = @time(runCEO(trainingXArrays, trainingYArrays, bestModelConfigs[m, :], gamma))

    expectedReturns[1:10] = generateExpectedReturns(betaArray, trainingXArrays, trainingYArrays, validationXRows)

    if m == 1
        expectedReturnMatrix[t, 1:10] = (exp.(expectedReturns[1:10])-1)*100

    end

    valY = zeros(10)
    for i = 1:10
        valY[i] = validationY[i][1]
    end

    trainX = trainingXArrays[1][:,1:10]
    Sigma =  cov(trainX)
    F = lufact(Sigma)
    U = F[:U]  #Cholesky factorization of Sigma

    return1N, returnCEO, wStar, forecastRow = performMVOptimization(expectedReturns, U, gamma, valY, valY)
    weightsCEO[t, 1:10, Int64(m)] = wStar
    return1NMatrix[t]      = return1N
    returnCEOMatrix[t, Int64(bestModelIndexes[m])] = returnCEO
    returnPerfect = returnPerfectMatrix[t]
    if m == 1
        forecastErrors[t, 1:10]  = forecastRow-expectedReturnMatrix[t,:]
    end
    println("1N returns is $return1N, returnPerfect is $returnPerfect and returnCEO is $returnCEO")
    PMatrix[t, Int64(bestModelIndexes[m])] = calculatePvalue(return1N, returnPerfect, returnCEO)
end

#Output data
combinedPortfolios = hcat(returnPerfectMatrix[1:nRows-trainingSize-2, 1], return1NMatrix[1:nRows-trainingSize-2, 1],
    returnCEOMatrix[1:nRows-trainingSize-2, Array{Int64}(bestModelIndexes)], PMatrix[1:nRows-trainingSize-2, Array{Int64}(bestModelIndexes)])
#path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results/CEO/"
path = "/zhome/9f/d/88706/SpecialeCode/Results/CEORFRSearch/Iter"*ARGS[2]*"/"
writedlm(path*string(inputArg1)*"_returnPvalueOutcome1to200.csv", combinedPortfolios, ",")
for i = 1:bestModelAmount
    writedlm(path*"wightsCEO_Model"*string(i)*"_"*string(inputArg1)*".csv",weightsCEO[:,:,i], ",")
end

combinedPortfolios2 = hcat(expectedReturnMatrix, forecastErrors)
writedlm(path*string(inputArg1)*"_forecastsAndErrors1to200.csv", combinedPortfolios2, ",")
println("Finished Everything")

max_stars_repo_name: skipperafrosenborg/Thesis

max_issues_repo_name: skipperafrosenborg/Thesis

max_forks_repo_name: skipperafrosenborg/Thesis


text:
using Optim
export dmdobjective, projBl2

# This file is available under the terms of the MIT License

@doc """
Functions for evaluating the objective and gradient.
"""

function prox_null(x)
    return
end

function dmdobjective(params;updatep=true)
    if updatep
        return aBFunc(params)
    else
        return BFunc(params)
    end
end

function BFunc(params)
    T = eltype(params.X);
    # update residual
    BLAS.gemm!('N', 'N', T(1.0), params.P, params.B, T(0.0), params.R);
    params.R .-= params.X
    #
    return params.lossFunc(view(params.R, :, params.ikeep))
end

function bFunc(params, id)
    T = eltype(params.X);
    # update residual
    BLAS.gemv!('N', T(1.0), params.P, params.b[id], T(0.0), params.r[id]);
    params.r[id] .-= params.x[id]
    #
    return params.lossFunc(params.r[id])
end

function bGrad(gbr, params, id)
    T = eltype(params.X);
    Tr = typeof(real(params.X[1]));
    # update residual
    BLAS.gemv!('N', T(1.0), params.P, params.b[id], T(0.0), params.r[id]);
    params.r[id] .-= params.x[id];
    #
    params.lossGrad(params.r[id]);
    #
    gb = vr2vc(gbr);
    BLAS.gemv!('T', T(1.0), params.P, params.r[id], T(0.0), gb);
    gbr[1:2:end] .*= Tr(2.0);
    gbr[2:2:end] .*= Tr(-2.0);
end

function bFuncGrad(gbr, params, id)
    T = eltype(params.X);
    Tr = typeof(real(params.X[1]));
    # update residual
    BLAS.gemv!('N', T(1.0), params.P, params.b[id], T(0.0), params.r[id]);
    params.r[id] .-= params.x[id];
    #

    f = params.lossFunc(params.r[id])

    params.lossGrad(params.r[id]);
    #
    gb = vr2vc(gbr);
    BLAS.gemv!('T', T(1.0), params.P, params.r[id], T(0.0), gb);
    gbr[1:2:end] .*= Tr(2.0);
    gbr[2:2:end] .*= Tr(-2.0);

    return f
end

function bfg!(ft, gbr, br, params, id)
    copyto!(params.br[id], br);
    if ft != nothing && gbr != nothing
        return bFuncGrad(gbr, params, id);
    end
    if ft == nothing && gbr != nothing
        bGrad(gbr, params, id)
        return
    end
    if ft != nothing && gbr == nothing
        return bFunc(params, id)
    end
end

function projb(params, id)

    # define the function and gradient interface for optim
    #

    function obfg!(f, g, x)
        return bfg!(f, g, x, params, id)
    end

    optimfg! = Optim.only_fg!(obfg!)
    res = optimize(optimfg!, params.br[id], params.inner_solver,
                   params.inner_opts);

    copyto!(params.br[id], res.minimizer);

end

function projB(params)
    for id = 1:params.n
        projb(params, id);
    end
end

function abFunc(params, id)
    # update P
    update_P!(params);
    # partially minimize over b
    if params.inner_directl2
        projbl2(params,id);
    else
        projb(params,id);
    end
    #
    return bFunc(params, id)
end

function aBFunc(params)
    # update P
    update_P!(params);
    # partially minimize over B
    if params.inner_directl2
        projBl2(params);
    else
        projB(params);
    end
    #
    return BFunc(params)
end

function abGrad(gar, params, id)
    # update P
    update_P!(params);
    #
    T = eltype(params.X);
    Tr = typeof(real(params.X[1]));
    # partially minimize over b

    if params.inner_directl2
        projbl2(params,id);
    else
        projb(params,id);
    end
    # update residual
    BLAS.gemv!('N', T(1.0), params.P, params.b[id], T(0.0), params.r[id]);
    params.r[id] .-= params.x[id];
    #
    params.lossGrad(params.r[id]);
    #
    ga = vr2vc(gar);
    #
    params.r[id] .*= params.t;
    BLAS.gemv!('T', T(1.0), params.P, params.r[id], T(0.0), ga);
    ga .*= params.b[id];
    #
    gar[1:2:end] .*= Tr(2.0);
    gar[2:2:end] .*= Tr(-2.0);
end

function aBGrad(gar, params)
    # update P
    update_P!(params);
    #
    T = eltype(params.X);
    Tr = typeof(real(params.X[1]));
    # partially minimize over B

    if params.inner_directl2
        projBl2(params);
    else
        projB(params);
    end

    # trim
    trimstandard(params)

    # update residual
    BLAS.gemm!('N', 'N', T(1.0), params.P, params.B, T(0.0), params.R);
    params.R .-= params.X
    #
    params.lossGrad(params.R);
    #
    ga = vr2vc(gar);
    #
    params.R .*= params.t;
    BLAS.gemm!('T', 'N', T(1.0), params.P, params.R, T(0.0), params.tM);
    params.tM .*= params.B;
    sum!(ga, view(params.tM, :, params.ikeep));
    #
    gar[1:2:end] .*= Tr(2.0);
    gar[2:2:end] .*= Tr(-2.0);
end

function aBFuncGrad(gar, params)
    # update P
    update_P!(params);
    #
    T = eltype(params.X);
    Tr = typeof(real(params.X[1]));
    # partially minimize over B
    if params.inner_directl2
        projBl2(params);
    else
        projB(params);
    end

    # trim
    trimstandard(params)

    # update residual
    BLAS.gemm!('N', 'N', T(1.0), params.P, params.B, T(0.0), params.R);
    params.R .-= params.X
    #

    f = params.lossFunc(view(params.R, :, params.ikeep))

    params.lossGrad(params.R);
    #
    ga = vr2vc(gar);
    #
    params.R .*= params.t;
    BLAS.gemm!('T', 'N', T(1.0), params.P, params.R, T(0.0), params.tM);
    params.tM .*= params.B;
    sum!(ga, view(params.tM, :, params.ikeep));
    #
    gar[1:2:end] .*= Tr(2.0);
    gar[2:2:end] .*= Tr(-2.0);

    return f
end

#####################################################
# L2 Penalty

function updateqr!(params)
    if params.a != params.qra
        qrfact = params.qrfact;

        # grab pointers
        fact = qrfact.factors
        jpvt = qrfact.jpvt
        tau = qrfact.τ

        # prep and call factorization
        copyto!(fact, params.P)
        jpvt .= 0 # do not forget!
        LinearAlgebra.LAPACK.geqp3!(fact, jpvt, tau)
        copyto!(params.qra, params.a)
    end
end

function projbl2(params, id)
    updateqr!(params)
    copyto!(params.qrtemp[id], params.x[id])
    ldiv!(params.qrfact, params.qrtemp[id])
    k = params.k
    copyto!(params.b[id], view(params.qrtemp[id], 1:k))
end

function projBl2(params)
    if any(isinf, params.P)
        params.B .= 0.0
    else
        updateqr!(params)
        copyto!(params.qrTemp, params.X)
        ldiv!(params.qrfact, params.qrTemp)
        k = params.k
        copyto!(params.B, view(params.qrTemp, 1:k, :))
    end
end

################################################
# trimming

function trimstandard(params)
    if params.n == params.nkeep
        return
    end
    T = eltype(params.X)
    # update residual
    BLAS.gemm!('N', 'N', T(1.0), params.P, params.B, T(0.0), params.R);
    params.R .-= params.X
    for id = 1:params.n
        params.tcnorm[id] = params.lossFunc(params.r[id])
    end
    sortperm!(params.idsort, params.tcnorm, rev=false)
    copyto!(params.ikeep, view(params.idsort, 1:params.nkeep))
    sort!(params.ikeep)
end


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler equations
gamma = 1.4
equations = CompressibleEulerEquations2D(gamma)

initial_condition = initial_condition_density_wave

surface_flux = flux_central
polydeg = 5
solver = DGSEM(polydeg, surface_flux)

coordinates_min = (-1, -1)
coordinates_max = ( 1,  1)
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=2,
                n_cells_max=30_000)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 2.0)
ode = semidiscretize(semi, tspan)

stepsize_callback = StepsizeCallback(cfl=1.6)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=stepsize_callback);


text:
# trendy
@testset "Trendlines" begin
    Random.seed!(SEED)
    @testset "Array" begin
        x = cumsum(randn(N))
        tmp = resistance(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = support(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = minima(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
        tmp = maxima(x)
        @test size(tmp, 1) == N
        @test size(tmp, 2) == 1
        @test sum(isnan.(tmp)) != N
    end
end


text:
module RiemannHilbert
using Base, ApproxFun, SingularIntegralEquations, DualNumbers, LinearAlgebra,
        SpecialFunctions, FillArrays, DomainSets, FastTransforms, SparseArrays

import DomainSets: UnionDomain, TypedEndpointsInterval

import FastTransforms: ichebyshevtransform!

import SingularIntegralEquations: stieltjesforward, stieltjesbackward, undirected, Directed, stieltjesmoment!, JacobiQ, istieltjes, ComplexPlane, ℂ,
                                    mxa_₂F₁, _₂F₁general, directed_mxa_₂F₁, directed_₂F₁general

import ApproxFunBase: mobius, pieces, npieces, piece, BlockInterlacer, interlacer, pieces_npoints,
                    ArraySpace, tocanonical, components_npoints, ScalarFun, VectorFun, MatrixFun,
                    dimension, evaluate, prectype, cfstype, Space, SumSpace, spacescompatible,
                    pieces
import ApproxFunOrthogonalPolynomials: PolynomialSpace, recA, recB, recC, IntervalOrSegmentDomain, IntervalOrSegment

import Base: values, convert, getindex, setindex!, *, +, -, ==, <, <=, >, |, !, !=, eltype,
                >=, /, ^, \, ∪, size, reindex, tail, broadcast, broadcast!,
                isinf, in

# we need to import all special functions to use Calculus.symbolic_derivatives_1arg
# we can't do importall Base as we replace some Base definitions
import Base: sinpi, cospi, exp,
                asinh, acosh,atanh,
                sin, cos, sinh, cosh,
                exp2, exp10, log2, log10,
                tan, tanh, csc, asin, acsc, sec, acos, asec,
                cot, atan, acot, sinh, csch, asinh, acsch,
                sech, acosh, asech, tanh, coth, atanh, acoth,
                expm1, log1p, sinc, cosc,
                abs, sign, log, expm1, tan, abs2, sqrt, angle, max, min, cbrt, log,
                atan, acos, asin, inv, real, imag, abs, conj

import LinearAlgebra: conj, transpose

import SpecialFunctions: airy, besselh, erfcx, dawson, erf, erfi,
                airyai, airybi, airyaiprime, airybiprime,
                hankelh1, hankelh2, besselj, bessely, besseli, besselk,
                besselkx, hankelh1x, hankelh2x, lfact,
                erfinv, erfcinv, erfc, beta, lbeta,
                eta, zeta, gamma,  lgamma, polygamma, invdigamma, digamma, trigamma

import DualNumbers: Dual, realpart, epsilon, dual
import FillArrays: AbstractFill

export cauchymatrix, rhmatrix, rhsolve, ℂ, istieltjes, KdV

include("LogNumber.jl")

function component_indices(it::BlockInterlacer, N::Int, kr::UnitRange)
    ret = Vector{Int}()
    ind = 1
    k_end = last(kr)
    for (M,j) in it
        N == M && j > k_end && return ret
        N == M && j ∈ kr && push!(ret, ind)
        ind += 1
    end
    ret
end

function component_indices(it::BlockInterlacer{NTuple{N,<:AbstractFill{Bool}}}, k::Int, kr::AbstractUnitRange) where N
    b = length(it.blocks)
    k + (first(kr)-1)*b:b:k + (last(kr)-1)*b
end

component_indices(sp::Space, k...) = component_indices(interlacer(sp), k...)

# # function fpstieltjes(f::Fun,z::Dual)
# #     x = mobius(domain(f),z)
# #     if !isinf(mobius(domain(f),Inf))
# #         error("Not implemented")
# #     end
# #     cfs = coefficients(f,Chebyshev)
# #     if realpart(x) ≈ 1
# #         c = -(log(dualpart(x))-log(2)) * sum(cfs)
# #         r = 0.0
# #         for k=2:2:length(cfs)-1
# #             r += 1/(k-1)
# #             c += -r*4*cfs[k+1]
# #         end
# #         r = 1.0
# #         for k=1:2:length(cfs)-1
# #             r += 1/(k-2)
# #             c += -(r+1/(2k))*4*cfs[k+1]
# #         end
# #         c
# #     elseif realpart(x) ≈ -1
# #         v = -(log(-dualpart(x))-log(2))
# #         if !isempty(cfs)
# #             c = -v*cfs[1]
# #         end
# #         r = 0.0
# #         for k=2:2:length(cfs)-1
# #             r += 1/(k-1)
# #             c += r*4*cfs[k+1]
# #             c += -v*cfs[k+1]
# #         end
# #         r = 1.0
# #         for k=1:2:length(cfs)-1
# #             r += 1/(k-2)
# #             c += -(r+1/(2k))*4*cfs[k+1]
# #             c += v*cfs[k+1]
# #         end
# #         c
# #     else
# #         error("Not implemented")
# #     end
# # end
# #
# # fpcauchy(x...) = fpstieltjes(x...)/(-2π*im)
#
#
#
# function stieltjesmatrix(space,pts::Vector,s::Bool)
#     n=length(pts)
#     C=Array(ComplexF64,n,n)
#     for k=1:n
#          C[k,:] = stieltjesforward(s,space,n,pts[k])
#     end
#     C
# end
#
# function stieltjesmatrix(space,pts::Vector)
#     n=length(pts)
#     C=zeros(ComplexF64,n,n)
#     for k=1:n
#         cfs = stieltjesbackward(space,pts[k])
#         C[k,1:min(length(cfs),n)] = cfs
#     end
#
#     C
# end

# stieltjesmatrix(space,n::Integer,s::Bool)=stieltjesmatrix(space,points(space,n),s)
# stieltjesmatrix(space,space2,n::Integer)=stieltjesmatrix(space,points(space2,n))

orientedleftendpoint(d::IntervalOrSegment) = RiemannDual(leftendpoint(d), sign(d))
orientedrightendpoint(d::IntervalOrSegment) = RiemannDual(rightendpoint(d), -sign(d))

# use 2nd kind to include endpoints
collocationpoints(d::IntervalOrSegmentDomain, m::Int) = points(d, m; kind=2)
collocationpoints(d::UnionDomain, ms::AbstractVector{Int}) = vcat(collocationpoints.(pieces(d), ms)...)
collocationpoints(d::UnionDomain, m::Int) = collocationpoints(d, pieces_npoints(d,m))

collocationpoints(sp::Space, m) = collocationpoints(domain(sp), m)

collocationvalues(f::ScalarFun, n) = f.(collocationpoints(space(f), n))
collocationvalues(f::Fun{<:Chebyshev}, n) = ichebyshevtransform!(coefficients(pad(f,n)); kind=2)
function collocationvalues(f::VectorFun, n)
    m = n÷size(f,1)
    mapreduce(f̃ -> collocationvalues(f̃,m), vcat, f)
end
function collocationvalues(f::MatrixFun, n)
    M = size(f,2)
    ret = Array{cfstype(f)}(undef, n, M)
    for J=1:M
        ret[:,J] = collocationvalues(f[:,J], n)
    end
    ret
end

collocationvalues(f::Fun{<:PiecewiseSpace}, n) = vcat(collocationvalues.(components(f), pieces_npoints(domain(f),n))...)

function evaluationmatrix!(E, sp::PolynomialSpace, x)
    x .= real(tocanonical.(Ref(sp), x))

    E[:,1] .= 1
    E[:,2] .= (recA(Float64,sp,0) .* x .+ recB(Float64,sp,0)) .* view(E,:,1)
    for j = 3:size(E,2)
        E[:,j] .= (recA(Float64,sp,j-2) .* x .+ recB(Float64,sp,j-2)) .* view(E,:,j-1) .- recC(Float64,sp,j-2).*view(E,:,j-2)
    end
    E
end

evaluationmatrix!(E, sp::PolynomialSpace) =
    evaluationmatrix!(E, sp, collocationpoints(sp, size(E,1)))

evaluationmatrix(sp::PolynomialSpace, x, n) =
    evaluationmatrix!(Array{Float64}(undef, length(x), n), sp,copy(x))

function evaluationmatrix!(C, sp::PiecewiseSpace, ns::AbstractVector{Int}, ms::AbstractVector{Int})
    N, M = length(ns), length(ms)
    @assert N == M == npieces(sp)
    n, m = sum(ns), sum(ms)
    @assert size(C) == (n,m)

    C .= 0

    for J = 1:M
        jr = component_indices(sp, J, 1:ms[J])
        k_start = sum(view(ns,1:J-1))+1
        kr = k_start:k_start+ns[J]-1
        evaluationmatrix!(view(C, kr, jr), component(sp, J))
    end

    C
end

function evaluationmatrix!(C, sp::ArraySpace, ns::AbstractVector{Int}, ms::AbstractVector{Int})
    @assert length(ns) == length(ms) == length(sp)
    N = length(ns)

    n, m = sum(ns), sum(ms)
    @assert size(C) == (n,m)

    C .= 0

    for J = 1:N
        jr = component_indices(sp, J, 1:ms[J]) ∩ (1:m)
        k_start = sum(view(ns,1:J-1))+1
        kr = k_start:k_start+ns[J]-1
        evaluationmatrix!(view(C, kr, jr), sp[J])
    end

    C
end

evaluationmatrix!(C, sp::PiecewiseSpace) =
    evaluationmatrix!(C, sp, pieces_npoints(sp, size(C,1)), pieces_npoints(sp, size(C,2)))

evaluationmatrix!(C, sp::ArraySpace) =
    evaluationmatrix!(C, sp, components_npoints(sp, size(C,1)), components_npoints(sp, size(C,2)))

evaluationmatrix(sp::Space, n::Int) = evaluationmatrix!(Array{Float64}(undef,n,n), sp)

fprightstieltjesmoment!(V, sp) = stieltjesmoment!(V, sp, Directed{false}(orientedrightendpoint(domain(sp))), finitepart)
fpleftstieltjesmoment!(V, sp) = stieltjesmoment!(V, sp, Directed{false}(orientedleftendpoint(domain(sp))), finitepart)
fprightstieltjesmoment!(V, sp, d) = stieltjesmoment!(V, sp, orientedrightendpoint(d), finitepart)
fpleftstieltjesmoment!(V, sp, d) = stieltjesmoment!(V, sp, orientedleftendpoint(d), finitepart)

function fpstieltjesmatrix!(C, sp, d)
    m, n = size(C)
    pts = collocationpoints(d, m)
    if d == domain(sp)
        fprightstieltjesmoment!(view(C,1,:), sp)
        for k=2:m-1
            stieltjesmoment!(view(C,k,:), sp, Directed{false}(pts[k]))
        end
        fpleftstieltjesmoment!(view(C,m,:), sp)
    elseif leftendpoint(d) ∈ domain(sp) && rightendpoint(d) ∈ domain(sp)
        fprightstieltjesmoment!(view(C,1,:), sp, d)
        for k=2:m-1
            stieltjesmoment!(view(C,k,:), sp, pts[k])
        end
        fpleftstieltjesmoment!(view(C,m,:), sp, d)
    elseif leftendpoint(d) ∈ domain(sp)
        for k=1:m-1
            stieltjesmoment!(view(C,k,:), sp, pts[k])
        end
        fpleftstieltjesmoment!(view(C,m,:), sp, d)
    elseif rightendpoint(d) ∈ domain(sp)
        fprightstieltjesmoment!(view(C,1,:), sp, d)
        for k=2:m
            stieltjesmoment!(view(C,k,:), sp, pts[k])
        end
    else
        for k=1:m
            stieltjesmoment!(view(C,k,:), sp, pts[k])
        end
    end
    C
end

fpstieltjesmatrix!(C, sp) = fpstieltjesmatrix!(C, sp, domain(sp))

fpstieltjesmatrix(sp::Space, d::Domain, n::Int, m::Int) =
    fpstieltjesmatrix!(Array{ComplexF64}(undef, n, m), sp, d)

fpstieltjesmatrix(sp::Space, n::Int, m::Int) =
    fpstieltjesmatrix!(Array{ComplexF64}(undef, n, m), sp, domain(sp))

# we group points together by piece
function fpstieltjesmatrix!(C, sp::PiecewiseSpace, ns::AbstractVector{Int}, ms::AbstractVector{Int})
    N, M = length(ns), length(ms)
    @assert N == M == npieces(sp)
    n, m = sum(ns), sum(ms)
    @assert size(C) == (n,m)

    for J = 1:M
        jr = component_indices(sp, J, 1:ms[J])
        k_start = 1
        for K = 1:N
            k_end = k_start + ns[K] - 1
            kr = k_start:k_end
            fpstieltjesmatrix!(view(C, kr, jr), component(sp, J),  domain(component(sp, K)))
            k_start = k_end+1
        end
    end

    C
end

fpstieltjesmatrix(sp::PiecewiseSpace, ns::AbstractVector{Int}, ms::AbstractVector{Int}) =
    fpstieltjesmatrix!(Array{ComplexF64}(undef, sum(ns), sum(ms)), sp, ns, ms)

fpstieltjesmatrix!(C, sp::PiecewiseSpace) = fpstieltjesmatrix!(C, sp, pieces_npoints(sp, size(C,1)), pieces_npoints(sp, size(C,2)))
fpstieltjesmatrix(sp::PiecewiseSpace, n::Int, m::Int) = fpstieltjesmatrix(sp, pieces_npoints(sp, n), pieces_npoints(sp, m))

# we group indices together by piece
function fpstieltjesmatrix(sp::ArraySpace, ns::AbstractArray{Int}, ms::AbstractArray{Int})
    @assert size(ns) == size(ms) == size(sp)
    N = length(ns)

    n, m = sum(ns), sum(ms)
    C = zeros(ComplexF64, n, m)

    for J = 1:N
        jr = component_indices(sp, J, 1:ms[J]) ∩ (1:m)
        k_start = sum(view(ns,1:J-1))+1
        kr = k_start:k_start+ns[J]-1
        fpstieltjesmatrix!(view(C, kr, jr), sp[J])
    end

    C
end

fpstieltjesmatrix(sp::ArraySpace, n::Int, m::Int) =
    fpstieltjesmatrix(sp, reshape(pieces_npoints(sp, n), size(sp)), reshape(pieces_npoints(sp, m), size(sp)))

cauchymatrix(x...) = stieltjesmatrix(x...)/(-2π*im)
function fpcauchymatrix(x...)
    C = fpstieltjesmatrix(x...)
    C ./= (-2π*im)
    C
end

## riemannhilbert
function multiplicationmatrix(G, n)
    N, M = size(G)
    @assert N == M
    sp = space(G)
    ret = spzeros(cfstype(G), n, n)
    m = n ÷ N
    pts = collocationpoints(sp, m)
    for K=1:N,J=1:M
        kr = (K-1)*m .+ (1:m)
        jr = (J-1)*m .+ (1:m)
        V = view(ret, kr, jr)
        view(V, diagind(V)) .= collocationvalues(G[K,J],m)
    end
    ret
end

function rhmatrix(g::ScalarFun, n)
    sp = rhspace(g)
    C₋ = fpcauchymatrix(sp, n, n)
    g_v = collocationvalues(g-1, n)
    E = evaluationmatrix(sp, n)
    C₋ .= g_v .* C₋
    E .- C₋
end

function rhmatrix(g::MatrixFun, n)
    sp = vector_rhspace(g)
    C₋ = fpcauchymatrix(sp, n, n)
    G = multiplicationmatrix(g-I, n)
    E = evaluationmatrix(sp, n)
    E .- G*C₋
end

function rh_sie_solve(G::MatrixFun, n)
    sp = vector_rhspace(G)
    cfs = rhmatrix(G, n) \ (collocationvalues(G-I, n))
    U = hcat([Fun(sp, cfs[:,J]) for J=1:size(G,2)]...)
end

struct RHProblem{GTyp,CM,RHMTyp}
    G::GTyp
    C₋::CM
    RP::RHMTyp
end

# function RHProblem(G)
#     RHProblem(G,

scalar_rhspace(d::AbstractInterval) = Legendre(d)
scalar_rhspace(d::UnionDomain) = PiecewiseSpace(Legendre.(components(d)))
array_rhspace(sz, d::Domain) = ArraySpace(scalar_rhspace(d), sz)
vector_rhspace(sz1, d::Domain) = ArraySpace(scalar_rhspace(d), sz1)
vector_rhspace(f::Fun) = vector_rhspace(size(f,1), domain(f))

rhspace(g::Fun{<:ArraySpace}) = array_rhspace(size(g), domain(g))
rhspace(g::Fun) = scalar_rhspace(domain(g))

rhsolve(g::ScalarFun, n) = 1+cauchy(Fun(rhspace(g), rhmatrix(g, n) \ (collocationvalues(g-1, n))))
function rhsolve(G::MatrixFun, n)
    U = rh_sie_solve(G, n)
    I+cauchy(U)
end

## AffineSpace

struct AffineSpace{DD,RR} <: Space{DD,RR}
    domain::DD
end

AffineSpace(d::Domain) = AffineSpace{typeof(d),prectype(d)}(d)
spacescompatible(::AffineSpace, ::AffineSpace) = true

dimension(::AffineSpace) = 2

function evaluate(v::AbstractVector{T}, s::AffineSpace, x::V) where {T,V}
    @assert length(v) ≤ 2
    (isempty(v) || x ∉ domain(s)) && return zero(promote_type(T,V))
    length(v) == 1 && return v[1] + zero(x)
    v[1] + v[2]*x
end

Fun(::typeof(identity), S::AffineSpace) = Fun(S, [0.0,1.0])
Fun(::typeof(identity), S::ComplexPlane) = Fun(Space(S), [0.0,1.0])

Space(d::ComplexPlane) = AffineSpace(d)

*(φ::Fun, z::Fun{<:AffineSpace}) = z*φ

function *(z::Fun{<:AffineSpace}, φ::Fun{<:JacobiQ})
    a = coefficient(z,1)
    b = coefficient(z,2)
    u = istieltjes(φ)
    x = Fun(domain(u))
    b*sum(u)+stieltjes(a*u + b*x*u)
end

*(z::Fun{<:AffineSpace}, φ::Fun{<:ConstantSpace}) = Fun(space(z),Number(φ)*coefficients(z))
*(z::Fun{<:AffineSpace}, Φ::Fun{<:SumSpace}) = mapreduce(f -> z*f, +, components(Φ))
*(z::Fun{<:AffineSpace}, Φ::Fun{<:ArraySpace}) = Fun(z.*Array(Φ))

include("KdV.jl")

function unorientedangles(ds, z₀)
    ret = Vector{Float64}()
    for d in ds
        if z₀ ≈ leftendpoint(d)
            push!(ret, angle(rightendpoint(d)-z₀))
        elseif z₀ ≈ rightendpoint(d)
            push!(ret, angle(leftendpoint(d)-z₀))
        else
            throw(ArgumentError("Must contain"))
        end
    end
    ret
end

function productcondition(G, z₀)
    Gs = filter(g -> z₀ ∈ domain(g),  pieces(G))
    p =  sortperm( unorientedangles(domain.(Gs), z₀))

    g₀ = Matrix{ComplexF64}(I, size(G))
    for g in Gs[p]
        if leftendpoint(domain(g)) ≈ z₀
            g₀ = g₀ * first(g)
        else
            g₀ = g₀ * inv(last(g))
        end
    end
    g₀
end

function productcondition(G)
    error("Implement")
end

end #module