text:
#using Ipopt: optimize!
using Base: Int64
using XLSX: length
#=using Pkg
Pkg.add("Tables")
Pkg.add("DataFrames")
Pkg.add("XLSX")
Pkg.add("ExcelReaders")
Pkg.add("JuMP")
Pkg.add("Ipopt")
=#
using XLSX, ExcelReaders, DataFrames, Tables, JuMP, Ipopt;
IOSource = XLSX.readdata("IO.xlsx", "io-table-5!A1:DV130");

#filepath cross system compatability code
if Sys.KERNEL === :NT || kern === :Windows
    pathmark = "\\"
else
    pathmark = "/"
end

#indexing vectors for initial data import groups
intermediaryTotalsCol = findall(x -> occursin("T4", x), string.(IOSource[3,:]));
intermediaryTotalsRow = findall(x -> occursin("T1", x), string.(IOSource[:,1]));
finalTotalsCol = findall(x -> occursin("T6", x), string.(IOSource[3,:]));
finalTotalsRow = findall(x -> occursin("Australian Production", x), string.(IOSource[:,2]));
finalDemandCol = findall(x -> occursin('Q', x), string.(IOSource[3,:]));
factorRow = findall(x -> occursin('P', x), string.(IOSource[:,1]));
IOSourceCol = vcat(intermediaryTotalsCol, finalDemandCol, finalTotalsCol);
IOSourceRow = vcat(intermediaryTotalsRow, factorRow, finalTotalsRow);
#initialising IO
IO = zeros(length(IOSourceRow), length(IOSourceCol));
#import numerical data into IO
IO[1:length(IOSourceRow), 1:length(IOSourceCol)] = IOSource[IOSourceRow, IOSourceCol];
#creating vectors of titles for IO
IOCodeRow = IOSource[IOSourceRow, 1];
IOCodeCol = IOSource[3, IOSourceCol];
IONameRow = IOSource[IOSourceRow, 2];
IONameCol = IOSource[2, IOSourceCol];

#code to sum public and private entities into one collumn
investment = findall(x -> occursin("Capital Formation", x), IONameCol);
IO[:, investment[1]]=sum(eachcol(IO[:, investment[1:2]]));
IO = IO[:,Not(investment[2])];
#alter title vectors accordingly (include Q in total investment collumn in IOcode)
IOCodeCol[investment[1]] = "Q3+Q4";
IOCodeCol = IOCodeCol[Not(investment[2])];
IONameCol[investment[1]] = "Private and Public Gross Fixed Capital Formation";
IONameCol = IONameCol[Not(investment[2])];
#creating a dictionary for the index of each collumn and row in IO by IOCode
IOColDict = Dict(IOCodeCol .=> [1:1:8;]);
IORowDict = Dict(IOCodeRow .=> [1:1:8;]);
IOCapForm = findall(x -> occursin("Capital Formation", x), IONameCol);
IOChangeInv = findall(x -> occursin("Changes in Inventories", x), IONameCol);

#importing relevant ASNA data for table 5
ASNAHouseCap = ExcelReaders.readxl("ASNAData"*pathmark*"5204039_Household_Capital_Account.xls", "Data1!A1:T71");
ASNANonFinCap = ExcelReaders.readxl("ASNAData"*pathmark*"5204018_NonFin_Corp_Capital_Account.xls", "Data1!A1:T71");
ASNAFinCap = ExcelReaders.readxl("ASNAData"*pathmark*"5204026_Fin_Corp_Capital_Account.xls", "Data1!A1:S71");
ASNAGovCap = ExcelReaders.readxl("ASNAData"*pathmark*"5204032_GenGov_Capital_Account.xls", "Data1!A1:AV71");
ASNAYearRow = findall(x -> occursin("2019", x), string.(ASNAHouseCap[:,1]));

#table 5
#creating table 5a - allocation of investment expenditure (broken into subsections for dict referencing purposes)
#subsection a is fixed capital expenditure
table5aNameCol = ["Households", "Non-Financial Corporations", "Financial Corporations", "General Government", "Total"];
table5aNameRow = ["Domestic Commodities", "Imported Commodities, complementary", "Imported Commodities, competing",
"Taxes less subsidies on products", "Other taxes less subsidies on investment", "Total indirect taxes",
"Total fixed capital expenditure"];
table5aRowDict = Dict(table5aNameRow .=> [1:1:length(table5aNameRow);]);
table5aColDict = Dict(table5aNameCol .=> [1:1:length(table5aNameCol);]);
table5a = zeros(length(table5aNameRow), length(table5aNameCol));

#filling in totals collumn from corresponding IO data
table5a[table5aRowDict["Domestic Commodities"], table5aColDict["Total"]] = sum(IO[IORowDict["T1"],IOCapForm]);
table5a[table5aRowDict["Imported Commodities, complementary"], table5aColDict["Total"]] = sum(IO[IORowDict["P5"],IOCapForm]);
table5a[table5aRowDict["Imported Commodities, competing"], table5aColDict["Total"]] = sum(IO[IORowDict["P6"],IOCapForm]);
table5a[table5aRowDict["Taxes less subsidies on products"], table5aColDict["Total"]] = sum(IO[IORowDict["P3"],IOCapForm]);
table5a[table5aRowDict["Other taxes less subsidies on investment"], table5aColDict["Total"]] = sum(IO[IORowDict["P4"],IOCapForm]);
table5aTaxes = findall(x -> occursin("taxes", lowercase(x)), table5aNameRow);
table5aTaxes = table5aTaxes[Not(3)];
table5a[table5aRowDict["Total indirect taxes"], table5aColDict["Total"]] = sum(table5a[table5aTaxes,table5aColDict["Total"]]);
table5a[table5aRowDict["Total fixed capital expenditure"], table5aColDict["Total"]] = sum(table5a[Not(table5aRowDict["Total indirect taxes"]),table5aColDict["Total"]]);

#creating index variables for the measurements that we want
ASNAHouseCapTotCapForm = findall(x -> occursin("Gross fixed capital formation ;", x), string.(ASNAHouseCap[1,:]));
ASNANonFinCapTotCapForm = findall(x -> occursin("Gross fixed capital formation ;", x), string.(ASNANonFinCap[1,:]));
ASNAFinCapTotCapForm = findall(x -> occursin("Gross fixed capital formation ;", x), string.(ASNAFinCap[1,:]));
ASNAGenGovCapTotCapForm = findall(x -> occursin("General government ;  Gross fixed capital formation ;", x), string.(ASNAGovCap[1,:]));

#filling in totals row from ASNA Data
table5a[table5aRowDict["Total fixed capital expenditure"], table5aColDict["Households"]] = first(ASNAHouseCap[ASNAYearRow, ASNAHouseCapTotCapForm]);
table5a[table5aRowDict["Total fixed capital expenditure"], table5aColDict["Non-Financial Corporations"]] = first(ASNANonFinCap[ASNAYearRow, ASNANonFinCapTotCapForm]);
table5a[table5aRowDict["Total fixed capital expenditure"], table5aColDict["Financial Corporations"]] = first(ASNAFinCap[ASNAYearRow, ASNAFinCapTotCapForm]);
table5a[table5aRowDict["Total fixed capital expenditure"], table5aColDict["General Government"]] = first(ASNAGovCap[ASNAYearRow, ASNAGenGovCapTotCapForm]);

#filling in non-total values
for ring in [1:1:length(table5aColDict)-1;];
    table5a[table5aRowDict["Domestic Commodities"],ring] = table5a[table5aRowDict["Total fixed capital expenditure"],ring]*IO[IORowDict["T1"],IOCapForm[1]] / IO[IORowDict[missing],IOCapForm[1]];
    table5a[table5aRowDict["Imported Commodities, complementary"],ring] = table5a[table5aRowDict["Total fixed capital expenditure"],ring]*IO[IORowDict["P5"],IOCapForm[1]] / IO[IORowDict[missing],IOCapForm[1]];
    table5a[table5aRowDict["Imported Commodities, competing"],ring] = table5a[table5aRowDict["Total fixed capital expenditure"],ring]*IO[IORowDict["P6"],IOCapForm[1]] / IO[IORowDict[missing],IOCapForm[1]];
    table5a[table5aRowDict["Taxes less subsidies on products"],ring] = table5a[table5aRowDict["Total fixed capital expenditure"],ring]*IO[IORowDict["P3"],IOCapForm[1]] / IO[IORowDict[missing],IOCapForm[1]];
    table5a[table5aRowDict["Other taxes less subsidies on investment"],ring] = table5a[table5aRowDict["Total fixed capital expenditure"],ring]*IO[IORowDict["P4"],IOCapForm[1]] / IO[IORowDict[missing],IOCapForm[1]];
    table5a[table5aRowDict["Total indirect taxes"], ring] = sum(table5a[table5aTaxes, ring]);
end

#creating table 5b - allocation of investment expenditure (broken into subsections for dict referencing purposes)
#subsection b is fixed capital expenditure
table5bNameCol = ["Households", "Non-Financial Corporations", "Financial Corporations", "General Government", "Total"];
table5bNameRow = ["Domestic Commodities", "Imported Commodities, complementary", "Imported Commodities, competing",
"Taxes less subsidies on products", "Total change in inventories"];
table5bRowDict = Dict(table5bNameRow .=> [1:1:length(table5bNameRow);]);
table5bColDict = Dict(table5bNameCol .=> [1:1:length(table5bNameCol);]);
table5b = zeros(length(table5bNameRow), length(table5bNameCol));

#filling in totals collumn from corresponding IO data
table5b[table5bRowDict["Domestic Commodities"], table5bColDict["Total"]] = sum(IO[IORowDict["T1"],IOChangeInv]);
table5b[table5bRowDict["Imported Commodities, complementary"], table5bColDict["Total"]] = sum(IO[IORowDict["P5"],IOChangeInv]);
table5b[table5bRowDict["Imported Commodities, competing"], table5bColDict["Total"]] = sum(IO[IORowDict["P6"],IOChangeInv]);
table5b[table5bRowDict["Taxes less subsidies on products"], table5bColDict["Total"]] = sum(IO[IORowDict["P3"],IOChangeInv]);
table5b[table5bRowDict["Total change in inventories"], table5bColDict["Total"]] = sum(table5b[:,table5bColDict["Total"]]);

#creating index variables for the measurements that we want
ASNAHouseCapChangeInv = findall(x -> occursin("Changes in inventories ;", x), string.(ASNAHouseCap[1,:]));
ASNANonFinCapChangeInv = findall(x -> occursin("Changes in inventories ;", x), string.(ASNANonFinCap[1,:]));
ASNAFinCapChangeInv = findall(x -> occursin("Changes in inventories ;", x), string.(ASNAFinCap[1,:]));
ASNAGenGovCapChangeInv = findall(x -> occursin("General government ;  Changes in inventories ;", x), string.(ASNAGovCap[1,:]));

#filling in totals row from ASNA Data
table5b[table5bRowDict["Total change in inventories"], table5bColDict["Households"]] = first(ASNAHouseCap[ASNAYearRow, ASNAHouseCapChangeInv]);
table5b[table5bRowDict["Total change in inventories"], table5bColDict["Non-Financial Corporations"]] = first(ASNANonFinCap[ASNAYearRow, ASNANonFinCapChangeInv]);
table5b[table5bRowDict["Total change in inventories"], table5bColDict["Financial Corporations"]] = first(ASNAFinCap[ASNAYearRow, ASNAFinCapChangeInv]);
table5b[table5bRowDict["Total change in inventories"], table5bColDict["General Government"]] = first(ASNAGovCap[ASNAYearRow, ASNAGenGovCapChangeInv]);

#calculate non-total values with lagrangian optimisation
table5bScalingFact = abs(minimum(table5b)) * 2;
mod5b = Model(Ipopt.Optimizer);
@variable(mod5b, x[1:(length(table5bNameRow)-1), 1:(length(table5bNameCol)-1)]);
@NLobjective(mod5b, Min, sum((x[i,j] - table5bScalingFact) ^ 2 for i in 1:(length(table5bNameRow)-1), j in 1:(length(table5bNameCol)-1)));
for i in 1:(length(table5bNameRow)-1);
    @constraint(mod5b, sum(x[:,i]) == table5b[table5bRowDict["Total change in inventories"],i]+table5bScalingFact);
end;
for i in 1:(length(table5bNameCol)-1);
    @constraint(mod5b, sum(x[i,:]) == table5b[i,table5bColDict["Total"]]+table5bScalingFact);
end;
optimize!(mod5b);

#plug back into table 5b
table5b[1:(length(table5bNameRow)-1),1:(length(table5bNameCol)-1)]=value.(x).-table5bScalingFact/4;

#creating table 5c - allocation of investment expenditure (broken into subsections for dict referencing purposes)
#subsection c is totals
table5cNameCol = ["Households", "Non-Financial Corporations", "Financial Corporations", "General Government", "Total"];
table5cNameRow = ["Domestic Commodities", "Imported Commodities", "Taxes less subsidies on products", "Other taxes less subsidies on investment", "Total investment expenditure"];
table5cRowDict = Dict(table5cNameRow .=> [1:1:length(table5cNameRow);]);
table5cColDict = Dict(table5cNameCol .=> [1:1:length(table5cNameCol);]);
table5c = zeros(length(table5cNameRow), length(table5cNameCol));

#do totals calcuations to get all values in 5c
table5c[table5cRowDict["Domestic Commodities"],:] = (table5a[table5aRowDict["Domestic Commodities"],:] +
table5b[table5bRowDict["Domestic Commodities"],:]);
table5c[table5cRowDict["Imported Commodities"],:] = sum(eachrow(table5a[[table5aRowDict["Imported Commodities, competing"],table5aRowDict["Imported Commodities, complementary"]],:] +
table5b[[table5bRowDict["Imported Commodities, competing"],table5bRowDict["Imported Commodities, complementary"]],:]));
table5c[table5cRowDict["Taxes less subsidies on products"],:] = table5a[table5aRowDict["Taxes less subsidies on products"],:];
table5c[table5cRowDict["Other taxes less subsidies on investment"],:] = (table5a[table5aRowDict["Other taxes less subsidies on investment"],:] +
table5b[table5bRowDict["Taxes less subsidies on products"],:]);
table5c[table5cRowDict["Total investment expenditure"],:] = (table5a[table5aRowDict["Total fixed capital expenditure"],:] +
table5b[table5bRowDict["Total change in inventories"],:]);

#table 6
#importing relevant ASNA data
ASNAHouseInc = ExcelReaders.readxl("ASNAData"*pathmark*"5204036_Household_Income_Account.xls", "Data1!A1:AN71");
ASNANonFinInc = ExcelReaders.readxl("ASNAData"*pathmark*"5204017_NonFin_Corp_Income_Account.xls", "Data1!A1:AE71");
ASNAFinInc = ExcelReaders.readxl("ASNAData"*pathmark*"5204025_Fin_Corp_Income_Account.xls", "Data1!A1:AD71");
ASNAGovInc = ExcelReaders.readxl("ASNAData"*pathmark*"5204030_GenGov_Income_Account.xls", "Data1!A1:DA71");
ASNAExtInc = ExcelReaders.readxl("ASNAData"*pathmark*"5204043_External_Accounts.xls", "Data1!A1:AI71");
#initialising table
table6Name = ["Households", "Non-Financial Corporations", "Financial Corporations", "General Government", "External", "Total"];
table6Dict = Dict(table6Name .=> [1:1:length(table6Name);]);
table6 = zeros(length(table6Name),length(table6Name));
#allocating total collumn and row data
table6[table6Dict["Total"],table6Dict["Households"]] = (first(ASNAHouseInc[ASNAYearRow,findall(x -> occursin("receivable - Interest", x), string.(ASNAHouseInc[1,:]))])
+first(ASNAHouseInc[ASNAYearRow,findall(x -> occursin("receivable - Imputed interest", x), string.(ASNAHouseInc[1,:]))]));
table6[table6Dict["Total"],table6Dict["Non-Financial Corporations"]] = (first(ASNANonFinInc[ASNAYearRow,findall(x -> occursin("receivable - Interest", x), string.(ASNANonFinInc[1,:]))])
+first(ASNANonFinInc[ASNAYearRow,findall(x -> occursin("Property income attributed to insurance policyholders", x), string.(ASNANonFinInc[1,:]))]));
table6[table6Dict["Total"],table6Dict["Financial Corporations"]] = first(ASNAFinInc[ASNAYearRow,findall(x -> occursin("receivable - Interest", x), string.(ASNAFinInc[1,:]))]);
table6[table6Dict["Total"],table6Dict["General Government"]] = first(ASNAGovInc[ASNAYearRow,findall(x -> occursin("General government ;  Property income receivable - Interest ;", x), string.(ASNAGovInc[1,:]))]);
table6[table6Dict["Total"],table6Dict["External"]] = first(ASNAExtInc[ASNAYearRow,findall(x -> occursin("receivable - Interest", x), string.(ASNAExtInc[1,:]))]);

table6[table6Dict["Households"],table6Dict["Total"]] = sum(ASNAHouseInc[ASNAYearRow,findall(x -> occursin("Property income payable - Interest", x), string.(ASNAHouseInc[1,:]))]);
table6[table6Dict["Non-Financial Corporations"],table6Dict["Total"]] = first(ASNANonFinInc[ASNAYearRow,findall(x -> occursin("Property income payable - Interest", x), string.(ASNANonFinInc[1,:]))]);
table6[table6Dict["Financial Corporations"],table6Dict["Total"]] = (first(ASNAFinInc[ASNAYearRow,findall(x -> occursin("Property income payable - Interest", x), string.(ASNAFinInc[1,:]))])
+first(ASNAFinInc[ASNAYearRow,findall(x -> occursin("payable - Property income attributed to insurance policy holders", x), string.(ASNAFinInc[1,:]))]));
table6[table6Dict["General Government"],table6Dict["Total"]] = first(ASNAGovInc[ASNAYearRow,findall(x -> occursin("General government ;  Property income payable - Total interest ;", x), string.(ASNAGovInc[1,:]))]);
table6[table6Dict["External"],table6Dict["Total"]] = first(ASNAExtInc[ASNAYearRow,findall(x -> occursin("Property income payable - Interest", x), string.(ASNAExtInc[1,:]))]);

if 0.98*sum(table6[:,length(table6Name)])<sum(table6[length(table6Name),:])<1.02*sum(table6[:,length(table6Name)])
else
    error("Large discrepancy in row and collumn total sums in table 6")
end

#=solve missing values with Ipopt
mod6 = Model(Ipopt.Optimizer);
@variable(mod6, x[1:(length(table6Name)-1), 1:(length(table6Name)-1)]>=0);
@NLobjective(mod6, Min, sum((x[i,j]) ^ 2 for i in 1:(length(table6Name)-1), j in 1:(length(table6Name)-1)));
for i in 1:(length(table6Name)-1);
    @constraint(mod6, sum(x[:,i]) == table6[table6Dict["Total"],i]-sum(table6[1:(length(table6Name)-1),i]));
end;
for i in 1:(length(table6Name)-1);
    @constraint(mod6, sum(x[i,:]) == table6[i,table6Dict["Total"]]-sum(table6[i,1:(length(table6Name)-1)]));
end;
for i in 1:(length(table6Name)-1);
    @constraint(mod6, x[i,table6Dict["General Government"]] == 0);
end;
@constraint(mod6, x[table6Dict["Households"],table6Dict["Households"]] == 0);
@constraint(mod6, x[table6Dict["External"],table6Dict["Households"]] == 0);
@constraint(mod6, x[table6Dict["Households"],table6Dict["External"]] == 0);
@constraint(mod6, x[table6Dict["External"],table6Dict["External"]] == 0);
optimize!(mod6);

table6[1:(length(table6Name)-1), 1:(length(table6Name)-1)] = table6[1:(length(table6Name)-1), 1:(length(table6Name)-1)] + value.(x);
=#

#solve for missing values with scaling
table6Step3 = zeros(length(table6Name),length(table6Name));
table6Step3Row = [table6Dict["Non-Financial Corporations"],table6Dict["Financial Corporations"],table6Dict["General Government"]];
table6Step3Col = [table6Dict["Households"],table6Dict["External"]];
for i in table6Step3Col;
    for ring in table6Step3Row;
        table6Step3[ring,i] = (table6[table6Dict["Total"],i]-sum(table6[1:(length(table6Name)-1),i]))*(
            table6[ring,table6Dict["Total"]]-sum(table6[ring,1:(length(table6Name)-1)]))/sum(table6[
            table6Step3Row,table6Dict["Total"]]-sum(eachcol(table6[table6Step3Row,1:(length(table6Name)-1)])));
    end
end
table6 = table6+table6Step3;

table6Step4 = zeros(length(table6Name),length(table6Name));
table6Step4Row = [1:1:(length(table6Name)-1);];
table6Step4Col = [table6Dict["Non-Financial Corporations"],table6Dict["Financial Corporations"]];
for i in table6Step4Col;
    for ring in table6Step4Row;
        table6Step4[ring,i] = (table6[table6Dict["Total"],i]-sum(table6[1:(length(table6Name)-1),i]))*(
            table6[ring,table6Dict["Total"]]-sum(table6[ring,1:(length(table6Name)-1)]))/sum(table6[
            table6Step4Row,table6Dict["Total"]]-sum(eachcol(table6[table6Step4Row,1:(length(table6Name)-1)])));
    end
end
table6 = table6+table6Step4;

#=convert dataframe to dictionary
function increment!( d::Dict{S, T}, k::S, i::T) where {T<:Real, S<:Any}
    if haskey(d, k)
        d[k] += i
    else
        d[k] = i
    end
end
increment!(d::Dict{S, T}, k::S ) where {T<:Real, S<:Any} = increment!( d, k, one(T))

function df2dict( df::DataFrame, key_col::Symbol, val_col::Symbol=:null)
    keytype = typeof(df[1,key_col])
    if val_col == :null
        valtype = Int
    else
        valtype = typeof(df[1,val_col])
    end
    D = Dict{keytype, valtype}()
    for i=1:size(df,1)
        if !ismissing(df[i,key_col])
            if val_col == :null
                increment!( D, df[i,key_col] )
            elseif valtype <: Real
                increment!( D, df[i,key_col], df[i,val_col] )
            else
                if haskey(D, df[i,key_col])
                    @warn("non-unique entry: $(df[i,key_col])")
                else
                    D[df[i,key_col]] = df[i,val_col]
                end
            end
        end
    end
    return D
end
df[!, "IOCode"]=IOcode
insertcols!(df, 2, :name => vector)
D = df2dict(df, :IOCode, :x3)
=#


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

@doc raw"""CronJobStatus represents the current state of a cron job.

    IoK8sApiBatchV2alpha1CronJobStatus(;
        active=nothing,
        lastScheduleTime=nothing,
    )

    - active::Vector{IoK8sApiCoreV1ObjectReference} : A list of pointers to currently running jobs.
    - lastScheduleTime::IoK8sApimachineryPkgApisMetaV1Time : Information when was the last time the job was successfully scheduled.
"""
mutable struct IoK8sApiBatchV2alpha1CronJobStatus <: SwaggerModel
    active::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1ObjectReference} } # spec name: active
    lastScheduleTime::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1Time } # spec name: lastScheduleTime

    function IoK8sApiBatchV2alpha1CronJobStatus(;active=nothing, lastScheduleTime=nothing)
        o = new()
        validate_property(IoK8sApiBatchV2alpha1CronJobStatus, Symbol("active"), active)
        setfield!(o, Symbol("active"), active)
        validate_property(IoK8sApiBatchV2alpha1CronJobStatus, Symbol("lastScheduleTime"), lastScheduleTime)
        setfield!(o, Symbol("lastScheduleTime"), lastScheduleTime)
        o
    end
end # type IoK8sApiBatchV2alpha1CronJobStatus

const _property_map_IoK8sApiBatchV2alpha1CronJobStatus = Dict{Symbol,Symbol}(Symbol("active")=>Symbol("active"), Symbol("lastScheduleTime")=>Symbol("lastScheduleTime"))
const _property_types_IoK8sApiBatchV2alpha1CronJobStatus = Dict{Symbol,String}(Symbol("active")=>"Vector{IoK8sApiCoreV1ObjectReference}", Symbol("lastScheduleTime")=>"IoK8sApimachineryPkgApisMetaV1Time")
Base.propertynames(::Type{ IoK8sApiBatchV2alpha1CronJobStatus }) = collect(keys(_property_map_IoK8sApiBatchV2alpha1CronJobStatus))
Swagger.property_type(::Type{ IoK8sApiBatchV2alpha1CronJobStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiBatchV2alpha1CronJobStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiBatchV2alpha1CronJobStatus }, property_name::Symbol) =  _property_map_IoK8sApiBatchV2alpha1CronJobStatus[property_name]

function check_required(o::IoK8sApiBatchV2alpha1CronJobStatus)
    true
end

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

ext: jl


text:
#=
CVaRO
=========
    CVaR Optimization Model
Author: Azamat Berdyshev
Date: 15/03/2016
=#

type CVaRO{R<:Real, S<:AbstractString} <: AbstractModel
    sense::Symbol
    vars::Vector
    _objective::Expr
    _default_constraints::Vector{Expr}
    constraints::Dict{Symbol,Expr}
    assets::AssetsCollection{R,S}
    objVal::R
    weights::Vector{R}
    status::Symbol

    # Inner constructor
    function CVaRO(assets::AssetsCollection{R, S},
                   losses::Matrix{R},
                   constraints::Dict{Symbol,Expr},
                   α::R,
                   short_sale::Bool)

        # n - number of assets
        # k - number of Monte-Carlo Scenarios
        # α - CVaR level
        # losses - Monte-Carlo sample of losses for every asset in every scenario

        (n,k) = size(losses)
        vars = Union{Expr, Symbol}[]

        # if no short sale => add corresponding constraint
        if short_sale
            push!(vars, :(w[1:$n]))
        else
            push!(vars, :(w[1:$n] >= 0))
        end
        push!(vars, :(y[1:$k] >= 0))
        push!(vars, :q)

        _objective = :(q + sum(y)/($k*(1-$α)))

        _default_constraints = [:(($(losses[:,i])'w)[1] - q - y[$i] <= 0) for i in 1:k]
        push!(_default_constraints, :((ones($n)'w)[1] == 1))

        # _default_constraints = [:(dot($(losses[:,i]), w) - q - y[$i] <= 0) for i in 1:k]

        # _default_constraints = [:($losses*w - q*ones($k) - y .≤ 0),
                                # :(dot(ones($n),w) == 1)]

        new(:Min,
            vars,
            _objective,
            _default_constraints,
            constraints,
            assets,
            NaN,
            fill(NaN,n),
            :Unsolved)
    end
end

# Outer constructor
CVaRO{R<:Real, S<:AbstractString}(
            assets::AssetsCollection{R, S},
            losses::Matrix{R},
            constraints=Dict{Symbol,Expr}()::Dict{Symbol,Expr};
            alpha=.95::R,
            short_sale=false::Bool) = CVaRO{R, S}(assets,
                                                  losses,
                                                  constraints,
                                                  alpha,
                                                  short_sale)


text:
#   Debug.AST:
# ==============
# Extensions to the julia AST format shared by decorate, instrument, and graft

module AST
using Compat
import Base.haskey, Base.show, Base.isequal, Base.promote_rule, Base.==

export Env, LocalEnv, NoEnv, child
export State, SimpleState, Def, Lhs, Rhs, Nonsyntax
export empty_symbol
export Loc, Plain, ExValue, Location
export Node, ExNode, PLeaf, SymNode, LocNode
export is_emittable, is_evaluable
export parentof, valueof, envof, exof
export Event, Enter, Leave

# ---- Env: analysis-time scope -----------------------------------------------

abstract Env
type NoEnv    <: Env; end
type LocalEnv <: Env
    source
    parent::Env
    defined::Set{Symbol}
    assigned::Set{Symbol}
end
child(source, env::Env) = LocalEnv(source, env, Set{Symbol}(), Set{Symbol}())

function show(io::IO, env::LocalEnv)
    print(io, "LocalEnv(,$(env.parent),$(env.defined),$(env.assigned))")
end

haskey(env::NoEnv,    sym::Symbol) = false
haskey(env::LocalEnv, sym::Symbol) = (sym in env.defined) || haskey(env.parent,sym)

add_defined( ::Env, ::Symbol) = nothing
add_assigned(::Env, ::Symbol) = nothing
add_defined( env::LocalEnv, sym::Symbol) = push!(env.defined,  sym)
add_assigned(env::LocalEnv, sym::Symbol) = push!(env.assigned, sym)

# ---- State: Context for a Node ----------------------------------------------

abstract State
promote_rule{S<:State,T<:State}(::Type{S},::Type{T}) = State

abstract SimpleState <: State
type Def <: SimpleState;  env::Env;  end  # definition, e.g. inside local
type Lhs <: SimpleState;  env::Env;  end  # e.g. to the left of =
type Rhs <: SimpleState;  env::Env;  end  # plain evaluation

type Nonsyntax <: State; end

envof(s::SimpleState) = s.env

# ---- Node: decorated AST node format ----------------------------------------

const empty_symbol = symbol("")

type Loc;   ex; line::Int; file::Symbol;  end
Loc(ex, line)       = Loc(ex, line, empty_symbol)

parent_type = VERSION >= v"0.4.0-dev" ? :(Union{Node,Void}) : :(Union(Node,Nothing))
@eval type Node{T}
    value::T
    parent::$parent_type
    introduces_scope::Bool
    state::State
    loc::Loc

    function adopt_args!{T}(node::Node{T})
        if T <: ExValue
            for arg in node.value.args
                set_parent(arg, node)
            end
        end
        node
    end

    Node(value::T)                 = adopt_args!(new(value,nothing,false))
    Node(value::T, s::State)       = adopt_args!(new(value,nothing,false,s))
    Node(value::T,s::State,l::Loc) = adopt_args!(new(value,nothing,false,s,l))
end
Node{T}(value::T, args...) = Node{T}(value, args...)

isequal(x::Node, y::Node)  = isequal(x.value, y.value)
==(x::Node, y::Node)       = isequal(x.value, y.value)
show(io::IO, ex::Node) = (print(io,"Node("); show(io,ex.value); print(io,")"))

function set_parent(ex::Node, parent::Node)
    if parentof(ex) === nothing; ex.parent = parent
    else; error("$ex already has a parent!")
    end
end

## Types for Node.value ##

type Plain; ex; end
type Location;  end
type ExValue
    head::Symbol
    args::Vector{Node}

    ExValue(head::Symbol, args) = new(head, Node[args...])
end

typealias ExNode  Node{ExValue}
typealias PLeaf   Node{Plain}
typealias SymNode Node{Symbol}
typealias LocNode Node{Location}

# override for node types that should not be emitted
is_emittable(ex) = true

is_evaluable(ex) = is_emittable(ex)
is_evaluable(::Node{Location}) = false
is_evaluable(::LocNode) = false

## Accessors for Node:s ##

parentof(node::Node) = node.parent
valueof( node::Node) = node.value
envof(   node::Node) = envof(node.state) # will only work State:s with envof

exof(node::SymNode) = valueof(node)
exof(node::LocNode) = node.loc.ex
exof(node::Node)    = exof(valueof(node))
exof(value::@compat(Union{Plain, Loc})) = value.ex

# ---- Events -----------------------------------------------------------------

abstract Event{T<:Node}

type Enter{T<:Node} <: Event{T};  node::Node;  end
type Leave{T<:Node} <: Event{T};  node::Node;  end

Enter{T<:Node}(node::T) = Enter{T}(node)
Leave{T<:Node}(node::T) = Leave{T}(node)

end # module


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM, Test

function calc_model(mesh_name; with_volume_load=false)
    meshfile = @__DIR__()*"/testdata/3d_block.med"
    mesh = aster_read_mesh(meshfile, mesh_name)

    block = Problem(Elasticity, "BLOCK", 3)
    block.properties.finite_strain = false
    block.properties.geometric_stiffness = false

    block_elements = create_elements(mesh, "BLOCK")
    update!(block_elements, "youngs modulus", 288.0)
    update!(block_elements, "poissons ratio", 1/3)
    with_volume_load && update!(block_elements, "displacement load 3", 576.0)
    add_elements!(block, block_elements)

    traction_elements = create_elements(mesh, "LOAD")
    update!(traction_elements, "displacement traction force 3", 288.0)
    add_elements!(block, traction_elements)

    bc = Problem(Dirichlet, "symmetry boundary condition", 3, "displacement")
    symyz = create_elements(mesh, "SYMYZ")
    symxz = create_elements(mesh, "SYMXZ")
    symxy = create_elements(mesh, "SYMXY")
    update!(symyz, "displacement 1", 0.0)
    update!(symxz, "displacement 2", 0.0)
    update!(symxy, "displacement 3", 0.0)
    add_elements!(bc, symyz, symxz, symxy)

    analysis = Analysis(Linear)
    add_problems!(analysis, block, bc)
    run!(analysis)

    u = block("displacement", 0.0)
    max_u = maximum(maximum(ui for ui in values(u)))
    return max_u
end

# 3d block HEX8
max_u = calc_model("BLOCK_HEX8"; with_volume_load=true)
@test isapprox(max_u, 2.0)

# 3d block TET4
max_u = calc_model("BLOCK_TET4"; with_volume_load=false)
@test isapprox(max_u, 1.0)

# 3d block TET10
max_u = calc_model("BLOCK_TET10"; with_volume_load=false)
@test isapprox(max_u, 1.0)


text:
export MWFarField3D

"""
Operator to compute the far field of a current distribution. In particular, given the current distribution ``j`` this operator allows for the computation of

```math
A j = n × ∫_Γ e^{γ ̂x ⋅ y} dy
```

where ``̂x`` is the unit vector in the direction of observation. Note that the assembly routing expects the observation directions to be normalised by the caller.
"""
type MWFarField3D{K}
  gamma::K
end

quaddata(op::MWFarField3D,rs,els) = quadpoints(rs,els,(3,))
quadrule(op::MWFarField3D,refspace,p,y,q,el,qdata) = qdata[1,q]

kernelvals(op::MWFarField3D,y,p) = exp(op.gamma*dot(y,cartesian(p)))
function integrand(op::MWFarField3D,krn,y,f,p)
    y × (krn * f[1])
end


text:
using Base: invokelatest

dummy() = return

function build_codeinfo(ir::IR)
  ir = copy(ir)
  ci = code_lowered(dummy, Tuple{})[1]
  ci.inlineable = true
  for arg in arguments(ir)
    push!(ci.slotnames, Symbol(""))
    push!(ci.slotflags, 0)
  end
  argument!(ir, at = 1)
  update!(ci, ir)
end

function func(m::Module, ir::IR)
  @eval @generated function $(gensym())($([Symbol(:arg, i) for i = 1:length(arguments(ir))]...))
    return build_codeinfo($ir)
  end
end

func(ir::IR) = func(Main, ir)

evalir(m::Module, ir::IR, args...) = invokelatest(func(m, ir), args...)
evalir(ir::IR, args...) = evalir(Main, ir, args...)


text:
using CSV
using DataFrames
using Ipopt
using JuMP

n_steps = 121 # Setting number of ODE discretisation steps

# Data
println("Reading measurement data...")
data_path = "/media/sf_DPhil_Project/Project07_Parameter Fitting/df_software/sbml2julia/examples/Vinod_FEBS2015/measurementData_Vinod_FEBS2015.tsv"
df = CSV.read(data_path)
insert!(df, 1, (1:length(df[:,1])), :id)

df_by_o = groupby(df, :observableId)
df_by_c = groupby(df, :simulationConditionId)

# Setting variables
t_sim = range(0, stop=maximum(df[!, :time]), length=n_steps)
t_exp = Dict()
m_exp = Dict()
t_sim_to_exp = Dict()
deduplicated_time_idx = Dict()

cond2idx = Dict()
for i_cond in 1:length(keys(df_by_c))
    cond_name = values(keys(df_by_c)[i_cond])[1]
    cond2idx[cond_name] = i_cond
end

for i_obs in 1:length(keys(df_by_o))
    obs_name = values(keys(df_by_o)[i_obs])[1]
    t_exp[obs_name] = Dict()
    m_exp[obs_name] = Dict()
    t_sim_to_exp[obs_name] = Dict()
    deduplicated_time_idx[obs_name] = Dict()
    df_by_o_c = groupby(DataFrame(df_by_o[i_obs]), :simulationConditionId)
    for i_cond in 1:length(keys(df_by_o_c))
        cond_name = values(keys(df_by_o_c)[i_cond])[1]
        cond_idx = cond2idx[cond_name]
        t_exp[obs_name][cond_idx] = df_by_o_c[i_cond][!, :time]
        m_exp[obs_name][cond_idx] = df_by_o_c[i_cond][!, :measurement]

        t = df_by_o_c[i_cond][!, :time]
        tmp = []
        for i in 1:length(t)
            idx = argmin(abs.(t[i] .- t_sim))
            append!(tmp, idx)
        end
        t_sim_to_exp[obs_name][cond_idx] = tmp

        deduplicated_time_idx[obs_name][cond_idx] = []
        idx = 0
        prev = "a"
        for t in df_by_o_c[i_cond][!, :time]
            if t != prev
                idx = idx + 1
            end
            push!(deduplicated_time_idx[obs_name][cond_idx], idx)
            prev = t
        end
    end
end

obs2conds = Dict()
for obs_dict in t_exp
    obs_name = obs_dict[1]
    obs2conds[obs_name] = []
    for cond_dict in obs_dict[2]
        cond_idx = cond_dict[1]
        push!(obs2conds[obs_name], cond_idx)
    end
end

results = Dict()
results["objective_value"] = Dict()
results["parameters"] = Dict()
results["species"] = Dict()
results["observables"] = Dict()

j_to_cond_par = [1, 2, 3, 4]
cond_without_preequ = [1, 2, 3, 4]
cond_with_preequ = []
preequ = []

i_start = 1
    m = Model(with_optimizer(Ipopt.Optimizer))

    set_optimizer_attribute(m,"linear_solver","MA57")

    # Define global parameters
    println("Defining global parameters...")
    @variable(m, 0.025 <= kDpEnsa <= 0.1, start=0.025+(0.1-(0.025))*rand(Float64))
    @variable(m, 0.5 <= kPhGw <= 2.0, start=0.5+(2.0-(0.5))*rand(Float64))
    @variable(m, 0.125 <= kDpGw1 <= 0.5, start=0.125+(0.5-(0.125))*rand(Float64))
    @variable(m, 5.0 <= kDpGw2 <= 20.0, start=5.0+(20.0-(5.0))*rand(Float64))
    @variable(m, 0.005 <= kWee1 <= 0.02, start=0.005+(0.02-(0.005))*rand(Float64))
    @variable(m, 0.495 <= kWee2 <= 1.98, start=0.495+(1.98-(0.495))*rand(Float64))
    @variable(m, 0.5 <= kPhWee <= 2.0, start=0.5+(2.0-(0.5))*rand(Float64))
    @variable(m, 5.0 <= kDpWee <= 20.0, start=5.0+(20.0-(5.0))*rand(Float64))
    @variable(m, 0.05 <= kCdc25_1 <= 0.2, start=0.05+(0.2-(0.05))*rand(Float64))
    @variable(m, 0.45 <= kCdc25_2 <= 1.8, start=0.45+(1.8-(0.45))*rand(Float64))
    @variable(m, 0.5 <= kPhCdc25 <= 2.0, start=0.5+(2.0-(0.5))*rand(Float64))
    @variable(m, 5.0 <= kDpCdc25 <= 20.0, start=5.0+(20.0-(5.0))*rand(Float64))
    @variable(m, 0.0034 <= kDipEB55 <= 0.0136, start=0.0034+(0.0136-(0.0034))*rand(Float64))
    @variable(m, 28.5 <= kAspEB55 <= 114.0, start=28.5+(114.0-(28.5))*rand(Float64))
    @variable(m, 1.0 <= fCb <= 4.0, start=1.0+(4.0-(1.0))*rand(Float64))
    @variable(m, 0.05 <= jiWee <= 0.2, start=0.05+(0.2-(0.05))*rand(Float64))
    @variable(m, 0.5 <= fB55_wt <= 2.0, start=0.5+(2.0-(0.5))*rand(Float64))
    @variable(m, 0.45 <= fB55_iWee <= 1.8, start=0.45+(1.8-(0.45))*rand(Float64))
    @variable(m, 0.55 <= fB55_Cb_low <= 2.2, start=0.55+(2.2-(0.55))*rand(Float64))
    @variable(m, 0.5 <= fB55_pGw_weak <= 2.0, start=0.5+(2.0-(0.5))*rand(Float64))
    @variable(m, 0.05 <= kPhEnsa_wt <= 0.2, start=0.05+(0.2-(0.05))*rand(Float64))
    @variable(m, 0.05 <= kPhEnsa_iWee <= 0.2, start=0.05+(0.2-(0.05))*rand(Float64))
    @variable(m, 0.05 <= kPhEnsa_Cb_low <= 0.2, start=0.05+(0.2-(0.05))*rand(Float64))
    @variable(m, 0.045 <= kPhEnsa_pGw_weak <= 0.18, start=0.045+(0.18-(0.045))*rand(Float64))

    # Define condition-specific parameters
    println("Defining condition-specific parameters...")
    @variable(m, iWee_0[1:4])
    @constraint(m, iWee_0[1] == 0.0)
    @constraint(m, iWee_0[2] == 0.1)
    @constraint(m, iWee_0[3] == 0.0)
    @constraint(m, iWee_0[4] == 0.0)

    @variable(m, Cb_0[1:4])
    @constraint(m, Cb_0[1] == 0.8)
    @constraint(m, Cb_0[2] == 0.8)
    @constraint(m, Cb_0[3] == 0.75)
    @constraint(m, Cb_0[4] == 0.8)

    @variable(m, fB55[1:4])
    @constraint(m, fB55[1] == fB55_wt)
    @constraint(m, fB55[2] == fB55_iWee)
    @constraint(m, fB55[3] == fB55_Cb_low)
    @constraint(m, fB55[4] == fB55_pGw_weak)

    @variable(m, kPhEnsa[1:4])
    @constraint(m, kPhEnsa[1] == kPhEnsa_wt)
    @constraint(m, kPhEnsa[2] == kPhEnsa_iWee)
    @constraint(m, kPhEnsa[3] == kPhEnsa_Cb_low)
    @constraint(m, kPhEnsa[4] == kPhEnsa_pGw_weak)

    # Define overrides
    # Define observable overrides
    println("Defining observableParameter overrides...")

    # Define noise overrides
    println("Defining noiseParameter overrides...")

    # Model compartments
    println("Defining compartments...")
    @variable(m, compartment == 1.0, start=1.0)

    # Model species
    println("Defining species...")
    @variable(m, 0 <= Cb[j in 1:4, k in 1:(length(t_sim)+1)] <= 1.1)
    @variable(m, 0 <= pCb[j in 1:4, k in 1:(length(t_sim)+1)] <= 1.1)
    @variable(m, 0 <= Wee[j in 1:4, k in 1:(length(t_sim)+1)] <= 1.1)
    @variable(m, 0 <= pWee[j in 1:4, k in 1:(length(t_sim)+1)] <= 1.1)
    @variable(m, 0 <= Cdc25[j in 1:4, k in 1:(length(t_sim)+1)] <= 1.1)
    @variable(m, 0 <= pCdc25[j in 1:4, k in 1:(length(t_sim)+1)] <= 1.1)
    @variable(m, 0 <= Gw[j in 1:4, k in 1:(length(t_sim)+1)] <= 1.1)
    @variable(m, 0 <= pGw[j in 1:4, k in 1:(length(t_sim)+1)] <= 1.1)
    @variable(m, 0 <= Ensa[j in 1:4, k in 1:(length(t_sim)+1)] <= 1.1)
    @variable(m, 0 <= pEnsa[j in 1:4, k in 1:(length(t_sim)+1)] <= 1.1)
    @variable(m, 0 <= pEB55[j in 1:4, k in 1:(length(t_sim)+1)] <= 1.1)
    @variable(m, 0 <= B55[j in 1:4, k in 1:(length(t_sim)+1)] <= 1.1)
    @variable(m, iWee[j in 1:4, k in 1:(length(t_sim)+1)])

    # Model initial assignments
    println("Defining initial assignments...")
    @constraint(m, [j in 1:4], Cb[cond_without_preequ[j],1] == Cb_0[cond_without_preequ[j]])
    @constraint(m, [j in 1:4], iWee[cond_without_preequ[j],1] == iWee_0[cond_without_preequ[j]])

    # Model ODEs
    println("Defining ODEs...")
    println("Cb")
    @NLconstraint(m, [j in 1:4, k in 1:length(t_sim)-1],
        Cb[j, k+1] == Cb[j, k] + ( -1.0*( compartment * kWee2 * Cb[j, k+1] * 0.5 * (Wee[j, k+1] - iWee[j] - jiWee + sqrt(^(-Wee[j, k+1] + iWee[j] + jiWee, 2) + 4 * jiWee * Wee[j, k+1])) ) -1.0*( compartment * kWee1 * Cb[j, k+1] ) +1.0*( compartment * kCdc25_1 * pCb[j, k+1] ) +1.0*( compartment * kCdc25_2 * pCb[j, k+1] * pCdc25[j, k+1] )     ) * ( t_sim[k+1] - t_sim[k] ) )
    println("pCb")
    @NLconstraint(m, [j in 1:4, k in 1:length(t_sim)-1],
        pCb[j, k+1] == pCb[j, k] + ( +1.0*( compartment * kWee2 * Cb[j, k+1] * 0.5 * (Wee[j, k+1] - iWee[j] - jiWee + sqrt(^(-Wee[j, k+1] + iWee[j] + jiWee, 2) + 4 * jiWee * Wee[j, k+1])) ) +1.0*( compartment * kWee1 * Cb[j, k+1] ) -1.0*( compartment * kCdc25_1 * pCb[j, k+1] ) -1.0*( compartment * kCdc25_2 * pCb[j, k+1] * pCdc25[j, k+1] )     ) * ( t_sim[k+1] - t_sim[k] ) )
    println("Wee")
    @NLconstraint(m, [j in 1:4, k in 1:length(t_sim)-1],
        Wee[j, k+1] == Wee[j, k] + ( -1.0*( compartment * kPhWee * Cb[j, k+1] * Wee[j, k+1] ) +1.0*( compartment * kDpWee * pWee[j, k+1] * B55[j, k+1] )     ) * ( t_sim[k+1] - t_sim[k] ) )
    println("pWee")
    @NLconstraint(m, [j in 1:4, k in 1:length(t_sim)-1],
        pWee[j, k+1] == pWee[j, k] + ( +1.0*( compartment * kPhWee * Cb[j, k+1] * Wee[j, k+1] ) -1.0*( compartment * kDpWee * pWee[j, k+1] * B55[j, k+1] )     ) * ( t_sim[k+1] - t_sim[k] ) )
    println("Cdc25")
    @NLconstraint(m, [j in 1:4, k in 1:length(t_sim)-1],
        Cdc25[j, k+1] == Cdc25[j, k] + ( -1.0*( compartment * kPhCdc25 * Cb[j, k+1] * Cdc25[j, k+1] ) +1.0*( compartment * kDpCdc25 * pCdc25[j, k+1] * B55[j, k+1] )     ) * ( t_sim[k+1] - t_sim[k] ) )
    println("pCdc25")
    @NLconstraint(m, [j in 1:4, k in 1:length(t_sim)-1],
        pCdc25[j, k+1] == pCdc25[j, k] + ( +1.0*( compartment * kPhCdc25 * Cb[j, k+1] * Cdc25[j, k+1] ) -1.0*( compartment * kDpCdc25 * pCdc25[j, k+1] * B55[j, k+1] )     ) * ( t_sim[k+1] - t_sim[k] ) )
    println("Gw")
    @NLconstraint(m, [j in 1:4, k in 1:length(t_sim)-1],
        Gw[j, k+1] == Gw[j, k] + ( -1.0*( compartment * kPhGw * Gw[j, k+1] * Cb[j, k+1] ) +1.0*( compartment * kDpGw1 * pGw[j, k+1] ) +1.0*( compartment * kDpGw2 * pGw[j, k+1] * B55[j, k+1] )     ) * ( t_sim[k+1] - t_sim[k] ) )
    println("pGw")
    @NLconstraint(m, [j in 1:4, k in 1:length(t_sim)-1],
        pGw[j, k+1] == pGw[j, k] + ( +1.0*( compartment * kPhGw * Gw[j, k+1] * Cb[j, k+1] ) -1.0*( compartment * kDpGw1 * pGw[j, k+1] ) -1.0*( compartment * kDpGw2 * pGw[j, k+1] * B55[j, k+1] )     ) * ( t_sim[k+1] - t_sim[k] ) )
    println("Ensa")
    @NLconstraint(m, [j in 1:4, k in 1:length(t_sim)-1],
        Ensa[j, k+1] == Ensa[j, k] + ( -1.0*( compartment * kPhEnsa[j_to_cond_par[j]] * Ensa[j, k+1] * pGw[j, k+1] ) +1.0*( compartment * kDpEnsa * pEB55[j, k+1] )     ) * ( t_sim[k+1] - t_sim[k] ) )
    println("pEnsa")
    @NLconstraint(m, [j in 1:4, k in 1:length(t_sim)-1],
        pEnsa[j, k+1] == pEnsa[j, k] + ( +1.0*( compartment * kPhEnsa[j_to_cond_par[j]] * Ensa[j, k+1] * pGw[j, k+1] ) -1.0*( compartment * kAspEB55 * B55[j, k+1] * pEnsa[j, k+1] ) +1.0*( compartment * kDipEB55 * pEB55[j, k+1] )     ) * ( t_sim[k+1] - t_sim[k] ) )
    println("pEB55")
    @NLconstraint(m, [j in 1:4, k in 1:length(t_sim)-1],
        pEB55[j, k+1] == pEB55[j, k] + ( +1.0*( compartment * kAspEB55 * B55[j, k+1] * pEnsa[j, k+1] ) -1.0*( compartment * kDipEB55 * pEB55[j, k+1] ) -1.0*( compartment * kDpEnsa * pEB55[j, k+1] )     ) * ( t_sim[k+1] - t_sim[k] ) )
    println("B55")
    @NLconstraint(m, [j in 1:4, k in 1:length(t_sim)-1],
        B55[j, k+1] == B55[j, k] + ( -1.0*( compartment * kAspEB55 * B55[j, k+1] * pEnsa[j, k+1] ) +1.0*( compartment * kDipEB55 * pEB55[j, k+1] ) +1.0*( compartment * kDpEnsa * pEB55[j, k+1] )     ) * ( t_sim[k+1] - t_sim[k] ) )
    @constraint(m, [j in 1:4, k in 1:length(t_sim)], iWee[j, k] == iWee_0[j])

    # Pre-equilibration constraints
    @constraint(m, [j in 1:4], Cb[j, length(t_sim)+1] == 0) # Dummy preequilibration
    @constraint(m, [j in 1:4], pCb[j, length(t_sim)+1] == 0) # Dummy preequilibration
    @constraint(m, [j in 1:4], Wee[j, length(t_sim)+1] == 0) # Dummy preequilibration
    @constraint(m, [j in 1:4], pWee[j, length(t_sim)+1] == 0) # Dummy preequilibration
    @constraint(m, [j in 1:4], Cdc25[j, length(t_sim)+1] == 0) # Dummy preequilibration
    @constraint(m, [j in 1:4], pCdc25[j, length(t_sim)+1] == 0) # Dummy preequilibration
    @constraint(m, [j in 1:4], Gw[j, length(t_sim)+1] == 0) # Dummy preequilibration
    @constraint(m, [j in 1:4], pGw[j, length(t_sim)+1] == 0) # Dummy preequilibration
    @constraint(m, [j in 1:4], Ensa[j, length(t_sim)+1] == 0) # Dummy preequilibration
    @constraint(m, [j in 1:4], pEnsa[j, length(t_sim)+1] == 0) # Dummy preequilibration
    @constraint(m, [j in 1:4], pEB55[j, length(t_sim)+1] == 0) # Dummy preequilibration
    @constraint(m, [j in 1:4], B55[j, length(t_sim)+1] == 0) # Dummy preequilibration
    @constraint(m, [j in 1:4], iWee[j, length(t_sim)+1] == 0) # Dummy preequilibration

    # Define observables
    println("Defining observables...")
    @variable(m, -1.2171211 <= obs_Cb[j in obs2conds["obs_Cb"], k in 1:length(t_exp["obs_Cb"][j])] <= 3.0085605500000003, start=1.)
    @NLconstraint(m, [j in obs2conds["obs_Cb"], k in 1:length(t_exp["obs_Cb"][j])], obs_Cb[j, k] == fCb*Cb[j, t_sim_to_exp["obs_Cb"][j][k]])
    @variable(m, -0.5130194399999999 <= obs_Gw[j in obs2conds["obs_Gw"], k in 1:length(t_exp["obs_Gw"][j])] <= 1.75650972, start=1.)
    @NLconstraint(m, [j in obs2conds["obs_Gw"], k in 1:length(t_exp["obs_Gw"][j])], obs_Gw[j, k] == 2*Gw[j, t_sim_to_exp["obs_Gw"][j][k]]/2)
    @variable(m, -0.585472911 <= obs_pEnsa[j in obs2conds["obs_pEnsa"], k in 1:length(t_exp["obs_pEnsa"][j])] <= 1.170945822, start=1.)
    @NLconstraint(m, [j in obs2conds["obs_pEnsa"], k in 1:length(t_exp["obs_pEnsa"][j])], obs_pEnsa[j, k] == 1+pEnsa[j, t_sim_to_exp["obs_pEnsa"][j][k]]-1)
    @variable(m, -0.994651293 <= obs_pWee[j in obs2conds["obs_pWee"], k in 1:length(t_exp["obs_pWee"][j])] <= 1.989302586, start=1.)
    @NLconstraint(m, [j in obs2conds["obs_pWee"], k in 1:length(t_exp["obs_pWee"][j])], obs_pWee[j, k] == pWee[j, t_sim_to_exp["obs_pWee"][j][k]])
    @variable(m, -0.989302586 <= obs_Cdc25[j in obs2conds["obs_Cdc25"], k in 1:length(t_exp["obs_Cdc25"][j])] <= 1.994651293, start=1.)
    @NLconstraint(m, [j in obs2conds["obs_Cdc25"], k in 1:length(t_exp["obs_Cdc25"][j])], obs_Cdc25[j, k] == Cdc25[j, t_sim_to_exp["obs_Cdc25"][j][k]])
    @variable(m, -0.670096254 <= obs_Ensa[j in obs2conds["obs_Ensa"], k in 1:length(t_exp["obs_Ensa"][j])] <= 1.8350481269999999, start=1.)
    @NLconstraint(m, [j in obs2conds["obs_Ensa"], k in 1:length(t_exp["obs_Ensa"][j])], obs_Ensa[j, k] == Ensa[j, t_sim_to_exp["obs_Ensa"][j][k]])
    @variable(m, -0.75650972 <= obs_pGw[j in obs2conds["obs_pGw"], k in 1:length(t_exp["obs_pGw"][j])] <= 1.51301944, start=1.)
    @NLconstraint(m, [j in obs2conds["obs_pGw"], k in 1:length(t_exp["obs_pGw"][j])], obs_pGw[j, k] == pGw[j, t_sim_to_exp["obs_pGw"][j][k]])
    @variable(m, -0.249575216 <= obs_pEB55[j in obs2conds["obs_pEB55"], k in 1:length(t_exp["obs_pEB55"][j])] <= 0.499150432, start=1.)
    @NLconstraint(m, [j in obs2conds["obs_pEB55"], k in 1:length(t_exp["obs_pEB55"][j])], obs_pEB55[j, k] == pEB55[j, t_sim_to_exp["obs_pEB55"][j][k]])
    @variable(m, -0.989302586 <= obs_Wee[j in obs2conds["obs_Wee"], k in 1:length(t_exp["obs_Wee"][j])] <= 1.994651293, start=1.)
    @NLconstraint(m, [j in obs2conds["obs_Wee"], k in 1:length(t_exp["obs_Wee"][j])], obs_Wee[j, k] == Wee[j, t_sim_to_exp["obs_Wee"][j][k]])
    @variable(m, -0.994651293 <= obs_pCdc25[j in obs2conds["obs_pCdc25"], k in 1:length(t_exp["obs_pCdc25"][j])] <= 1.989302586, start=1.)
    @NLconstraint(m, [j in obs2conds["obs_pCdc25"], k in 1:length(t_exp["obs_pCdc25"][j])], obs_pCdc25[j, k] == pCdc25[j, t_sim_to_exp["obs_pCdc25"][j][k]])
    @variable(m, -0.274235388 <= obs_B55[j in obs2conds["obs_B55"], k in 1:length(t_exp["obs_B55"][j])] <= 0.549617694, start=1.)
    @NLconstraint(m, [j in obs2conds["obs_B55"], k in 1:length(t_exp["obs_B55"][j])], obs_B55[j, k] == B55[j, t_sim_to_exp["obs_B55"][j][k]]*fB55[j])
    @variable(m, -0.695490749 <= obs_pCb[j in obs2conds["obs_pCb"], k in 1:length(t_exp["obs_pCb"][j])] <= 1.390981498, start=1.)
    @NLconstraint(m, [j in obs2conds["obs_pCb"], k in 1:length(t_exp["obs_pCb"][j])], obs_pCb[j, k] == pCb[j, t_sim_to_exp["obs_pCb"][j][k]])

    # Define objective
    println("Defining objective...")
    @NLobjective(m, Min, sum(0.5 * log(2*pi*(( 1 ))^2) + 0.5*((obs_Cb[j, deduplicated_time_idx["obs_Cb"][j][k]]-m_exp["obs_Cb"][j][k])/(( 1 )))^2 for j in obs2conds["obs_Cb"] for k in 1:length(t_exp["obs_Cb"][j]))
        + sum(0.5 * log(2*pi*(( 1 ))^2) + 0.5*((obs_Gw[j, deduplicated_time_idx["obs_Gw"][j][k]]-m_exp["obs_Gw"][j][k])/(( 1 )))^2 for j in obs2conds["obs_Gw"] for k in 1:length(t_exp["obs_Gw"][j]))
        + sum(0.5 * log(2*pi*(( 1 ))^2) + 0.5*((obs_pEnsa[j, deduplicated_time_idx["obs_pEnsa"][j][k]]-m_exp["obs_pEnsa"][j][k])/(( 1 )))^2 for j in obs2conds["obs_pEnsa"] for k in 1:length(t_exp["obs_pEnsa"][j]))
        + sum(0.5 * log(2*pi*(( 1 ))^2) + 0.5*((obs_pWee[j, deduplicated_time_idx["obs_pWee"][j][k]]-m_exp["obs_pWee"][j][k])/(( 1 )))^2 for j in obs2conds["obs_pWee"] for k in 1:length(t_exp["obs_pWee"][j]))
        + sum(0.5 * log(2*pi*(( 1 ))^2) + 0.5*((obs_Cdc25[j, deduplicated_time_idx["obs_Cdc25"][j][k]]-m_exp["obs_Cdc25"][j][k])/(( 1 )))^2 for j in obs2conds["obs_Cdc25"] for k in 1:length(t_exp["obs_Cdc25"][j]))
        + sum(0.5 * log(2*pi*(( 1 ))^2) + 0.5*((obs_Ensa[j, deduplicated_time_idx["obs_Ensa"][j][k]]-m_exp["obs_Ensa"][j][k])/(( 1 )))^2 for j in obs2conds["obs_Ensa"] for k in 1:length(t_exp["obs_Ensa"][j]))
        + sum(0.5 * log(2*pi*(( 1 ))^2) + 0.5*((obs_pGw[j, deduplicated_time_idx["obs_pGw"][j][k]]-m_exp["obs_pGw"][j][k])/(( 1 )))^2 for j in obs2conds["obs_pGw"] for k in 1:length(t_exp["obs_pGw"][j]))
        + sum(0.5 * log(2*pi*(( 1 ))^2) + 0.5*((obs_pEB55[j, deduplicated_time_idx["obs_pEB55"][j][k]]-m_exp["obs_pEB55"][j][k])/(( 1 )))^2 for j in obs2conds["obs_pEB55"] for k in 1:length(t_exp["obs_pEB55"][j]))
        + sum(0.5 * log(2*pi*(( 1 ))^2) + 0.5*((obs_Wee[j, deduplicated_time_idx["obs_Wee"][j][k]]-m_exp["obs_Wee"][j][k])/(( 1 )))^2 for j in obs2conds["obs_Wee"] for k in 1:length(t_exp["obs_Wee"][j]))
        + sum(0.5 * log(2*pi*(( 1 ))^2) + 0.5*((obs_pCdc25[j, deduplicated_time_idx["obs_pCdc25"][j][k]]-m_exp["obs_pCdc25"][j][k])/(( 1 )))^2 for j in obs2conds["obs_pCdc25"] for k in 1:length(t_exp["obs_pCdc25"][j]))
        + sum(0.5 * log(2*pi*(( 1 ))^2) + 0.5*((obs_B55[j, deduplicated_time_idx["obs_B55"][j][k]]-m_exp["obs_B55"][j][k])/(( 1 )))^2 for j in obs2conds["obs_B55"] for k in 1:length(t_exp["obs_B55"][j]))
        + sum(0.5 * log(2*pi*(( 1 ))^2) + 0.5*((obs_pCb[j, deduplicated_time_idx["obs_pCb"][j][k]]-m_exp["obs_pCb"][j][k])/(( 1 )))^2 for j in obs2conds["obs_pCb"] for k in 1:length(t_exp["obs_pCb"][j]))
        )

    println("Optimizing:")
    optimize!(m)

    println("Transfering results to Python...")
    parameter_names = [kDpEnsa, kPhGw, kDpGw1, kDpGw2, kWee1, kWee2, kPhWee, kDpWee, kCdc25_1, kCdc25_2, kPhCdc25, kDpCdc25, kDipEB55, kAspEB55, fCb, jiWee, fB55_wt, fB55_iWee, fB55_Cb_low, fB55_pGw_weak, kPhEnsa_wt, kPhEnsa_iWee, kPhEnsa_Cb_low, kPhEnsa_pGw_weak]
    parameter_values = Dict()
    for p in parameter_names
        if occursin("[", string(p))
            parameter_values[split(string(p[1]), "[")[1]] = JuMP.value.(p)
        else
            parameter_values[string(p)] = JuMP.value.(p)
        end
    end

    species_names = ["Cb", "pCb", "Wee", "pWee", "Cdc25", "pCdc25", "Gw", "pGw", "Ensa", "pEnsa", "pEB55", "B55", "iWee", ]
    species_values = Dict(string(spec)=>Dict(cond_idx=>[value(eval(Symbol(spec))[cond_idx, val_idx]) for val_idx=1:length(t_sim)]
        for cond_idx in 1:length(keys(df_by_c))) for spec in species_names)

    observable_values = Dict(obs=>Dict(cond_idx=>[value(eval(Symbol(obs))[cond_idx, val_idx]) for val_idx=1:length(list)]
        for (cond_idx, list) in dict) for (obs, dict) in m_exp)

    objective_val = objective_value(m)

    results["objective_value"][string(i_start)] = objective_val
    results["parameters"][string(i_start)] = parameter_values
    results["species"][string(i_start)] = species_values
    results["observables"][string(i_start)] = observable_values

results

max_stars_repo_name: paulflang/DisFit

max_issues_repo_name: paulflang/DisFit

max_forks_repo_name: paulflang/DisFit


text:
using Test

include("../src/matrixcore.jl")
using .matrixcore

@testset "test_matrixcore" begin

    @testset "StoMtoS" begin
        # Running random tests to check that MtoS and StoM are each others inverse
        for _ in 1:100
            S = rand(ComplexF64, 2, 2)
            @test S ≈ (StoM ∘ MtoS)(S)
        end
    end

    @testset "test_cascades" begin
        # Getting random transmittance and reflectances, testing two mirror systems
        # Saleh & Teich 3.ed eq.7.1-8
        for _ in 1:100
            t12 = rand(ComplexF64)
            t21 = rand(ComplexF64)

            r12 = rand(ComplexF64)
            r21 = rand(ComplexF64)

            t23 = rand(ComplexF64)
            t32 = rand(ComplexF64)

            r23 = rand(ComplexF64)
            r32 = rand(ComplexF64)

            ϕ = rand() * 2π

            mirror1 = [
                t12 r21;
                r12 t21
            ]
            air = [ℯ^(-ϕ * im) 0; 0 ℯ^(-ϕ * im)]
            mirror2 = [
                t23 r32;
                r23 t32
            ]

            machine = CascadeScattering([mirror1, air, mirror2])

            @test machine[1, 1] ≈ (t12 * t23 * ℯ^(-ϕ * im)) / (1 - (r21 * r23 * ℯ^(-ϕ * 2im)))
            @test machine[2, 1] ≈ r12 + (t12 * t21 * r23 * ℯ^(-ϕ * 2im)) / (1 - (r21 * r23 * ℯ^(-ϕ * 2im)))
        end
    end

end


text:
#******************************************************************************#
#                                                                              #
#                                                         :::      ::::::::    #
#    run_command.jl                                     :+:      :+:    :+:    #
#                                                     +:+ +:+         +:+      #
#    By: mngomane <mngomane@student.42.fr>          +#+  +:+       +#+         #
#                                                 +#+#+#+#+#+   +#+            #
#    Created: 2015/02/12 07:32:54 by mngomane          #+#    #+#              #
#    Updated: 2015/02/12 07:32:54 by mngomane         ###   ########.fr        #
#                                                                              #
#******************************************************************************#

function	run_command(cmd, arg)
	try
		run(`$cmd $arg`)
	end
end

function	tests_state(data)
	state = true
	data = split(data, '\n', false)
	for c in data
		try
			if c[end-3:end] != "[OK]"
				println("reached")
				state = false
			end
		end
	end
	return state
end

function	check_return()
	originalSTDOUT = STDOUT
	(outRead, outWrite) = redirect_stdout()
	run(`./libftasm`)
	close(outWrite)
	data = readavailable(outRead)
	redirect_stdout(originalSTDOUT)
	if tests_state(data) == true
		println(valid_string("OK"))
	else
		println(error_string("KO"))
	end
end

max_stars_repo_name: mngomane/projects

max_issues_repo_name: mngomane/projects

max_forks_repo_name: mngomane/projects


text:
fatalerrors = false
anyerrors = false
quiet = false

my_tests = ["dataframes.jl",
            "eachrow.jl",
            "grouping.jl",
            "subset.jl",
            "function_compilation.jl",
            "chaining.jl",
            "linqmacro.jl",
            "deprecated.jl",
            "byrow.jl",
            "astable.jl"]

println("Running tests:")

for my_test in my_tests
    try
        include(my_test)
        println("\t\033[1m\033[32mPASSED\033[0m: $(my_test)")
    catch e
        global anyerrors = true
        println("\t\033[1m\033[31mFAILED\033[0m: $(my_test)")
        if fatalerrors
            rethrow(e)
        elseif !quiet
            showerror(stdout, e, backtrace())
            println()
        end
    end
end

if anyerrors
    throw("Tests failed")
end


text:
using SparseArrays, LinearAlgebra, Arpack
include("io/io_helper.jl")

function diagonalize(ham::Hamiltonian, param::Dict{Any,Any})
    """ Performs the diagonalization as specified in the paramter dictionary.
        Done either by iterated Arnoldi method with ARPACK (when full_diag is
        set to `false`) or by full diagonalization.

        Notes:
         -  Julia has a matrix-independent iteration number of 300 per default,
            this is potentially problematic. Needs to be adjusted in some cases.
    """
    sparse_ham = sparse(ham.row, ham.col, ham.data, ham.n_fock, ham.n_fock)

    if ham.n_fock > 1
        @info "Starting to diagonalize the Hamiltonian." n_fock=Int(ham.n_fock)
        if get(param, "full_diag", false)
            # Full diagonalization with built-in method.
            time = @elapsed mem = @allocated F = eigen(Matrix(sparse_ham))
            @info "Done computing the full spectrum." time=time allocated=MemoryTag(mem) spectrum=F.values
            return F.values, F.vectors
        else
            # Arnoldi method with ARPACK.
            which_map = Dict([
                ("SA", :SR)
            ])
            # Here, we set the default values that are otherwise overwritten
            # if specified in the parameter dictionary.
            time = @elapsed mem = @allocated ev, est = eigs(
                sparse_ham,
                nev=get(param, "n_eigenvalues", 10),
                which=which_map[get(param, "ev_type", "SA")],
                tol=get(param, "tol", 1e-15),
                maxiter=get(param, "max_iter", 10*ham.n_fock)
            )
            @info "Done computing the lower spectrum." time=time allocated=MemoryTag(mem) spectrum=ev
            return ev, est
        end
    end
    @info "Trivial 1x1 matrix to diagonalize."
    return Matrix(sparse_ham)[1,1], [1]
end


text:
for (Func,Len) in ((:(Base.sum),:complexlength),(:linesum,:arclength))
    @eval begin
        function $Func(G::Function,u::Fun{<:JacobiWeight{<:Chebyshev}},z)
           d,α,β,n=domain(u),u.space.α,u.space.β,2ncoefficients(u)
           vals,t = ichebyshevtransform(pad(u.coefficients,n)),points(d,n)
           sp,p = space(u),plan_chebyshevtransform(complex(vals))
           map(z->$Func(Fun(sp,p*(G.(z,t).*vals))),z)
       end

       # TODO: remove the following hack
        function $Func(G::Function,u::Fun{<:JacobiWeight{<:Chebyshev,<:IntervalOrSegment}},z)
            d,α,β,n=domain(u),u.space.α,u.space.β,2ncoefficients(u)
            vals,t = ichebyshevtransform(pad(u.coefficients,n)),points(d,n)
            if α == β == -0.5
                return 0.5*$Len(d)*map(z->mean(G.(z,t).*vals),z)*π
            else
                sp,p = space(u),plan_chebyshevtransform(complex(vals))
                map(z->$Func(Fun(sp,p*(G.(z,t).*vals))),z)
            end
        end
    end
end

function logkernel(G::Function,u::Fun{<:JacobiWeight{<:Chebyshev}},z)
    sp,n=space(u),2ncoefficients(u)
    vals,t = ichebyshevtransform(pad(u.coefficients,n)),points(sp,n)
    p = plan_chebyshevtransform(complex(vals))
    return map(z->logkernel(Fun(sp,p*(G.(z,t).*vals)),z),z)
end

function cauchy(G::Function,u::Fun{<:JacobiWeight{<:Chebyshev}},z)
    sp,n=space(u),2ncoefficients(u)
    vals,t = ichebyshevtransform(pad(u.coefficients,n)),points(sp,n)
    p = plan_chebyshevtransform(vals)#complex(vals))
    return map(z->cauchy(Fun(sp,p*(G.(z,t).*vals)),z),z)
end

for TYP in (:Fourier,:Laurent)
    @eval begin
        function Base.sum(G::Function,u::Fun{<:$TYP{<:PeriodicSegment}},z)
            d,n=domain(u),2ncoefficients(u)
            vals,t = values(pad(u,n)),points(d,n)

            return map(z->mean(G.(z,t).*vals),z)*arclength(d)
        end

        function Base.sum(G::Function,u::Fun{<:$TYP{<:Circle}},z)
            d,n=domain(u),2ncoefficients(u)
            vals,t = values(pad(u,n)),points(d,n)
            #TODO: shouldn't this depend on which circle?
            return map(z->mean(G.(z,t).*vals.*t),z)*2π*im
        end

        function linesum(G::Function,u::Fun{<:$TYP{<:PeriodicSegment}},z)
            d,n=domain(u),2ncoefficients(u)
            vals,t = values(pad(u,n)),points(d,n)
            map(z->mean(G.(z,t).*vals),z)*arclength(d)
        end
    end
end

function logkernel(G::Function,sp::Fourier{DD,RR},u,z) where {DD,RR}
    n=2length(u)
    vals,t = values(pad(Fun(sp,u),n)),points(sp,n)
    p = plan_transform(sp,vals)
    return map(z->logkernel(sp,transform(sp,G.(z,t).*vals,p),z),z)
end
logkernel(G::Function,sp::Laurent{DD,RR},u,z) where {DD,RR} =
    logkernel(G,Fun(Fun(sp,u),Fourier),z)

for Func in (:(Base.sum),:linesum,:logkernel,:cauchy)
    @eval begin
        $Func(G::Function,u::Vector{F},z) where {F<:Fun} = mapreduce(u->$Func(G,u,z),+,u)
        $Func(G::Function,sp::PiecewiseSpace,u,z) = $Func(G,vec(Fun(sp,u)),z)
        $Func(G::Function,u::Fun,z) = $Func(G,space(u),coefficients(u),z)
        $Func(G::Function,sp::JacobiWeight{PS},f,z) where {PS<:PolynomialSpace} =
            $Func(G,Fun(Fun(sp,f),JacobiWeight(sp.β,sp.α,Chebyshev(domain(sp)))),z)
    end
end

function logkernel(G::Function,sp::Space{LS,RR},u,z) where {LS,RR<:Arc}
    n=2length(u)
    vals,t = itransform(sp,pad(u,n)),points(sp,n)
    p = plan_transform(sp,complex(vals))
    return map(z->logkernel(sp,transform(sp,G.(z,t).*vals,p),z),z)
end

function linesum(G::Function,sp::Space{LS,RR},u,z) where {LS,RR<:Arc}
    n=2length(u)
    vals,t = itransform(sp,pad(u,n)),points(sp,n)
    p = plan_transform(sp,complex(vals))
    return map(z->linesum(sp,transform(sp,G.(z,t).*vals,p),z),z)
end


text:
import MathProgBase
export can_solve_mip, can_solve_socp, can_solve_sdp, can_solve_exp
export set_default_solver, get_default_solver

function set_default_solver(solver::MathProgBase.AbstractMathProgSolver)
  global DEFAULT_SOLVER
  DEFAULT_SOLVER = solver
end

function get_default_solver()
  if DEFAULT_SOLVER == nothing
    error("The default solver is set to `nothing`
         You must have at least one solver installed to use Convex.
         You can install a solver such as SCS by running:
         Pkg.add(\"SCS\").
         You will have to restart Julia after that.")
  end
  return DEFAULT_SOLVER
end

solvers = [("SCS", "SCSSolver"), ("ECOS", "ECOSSolver"), ("Gurobi", "GurobiSolver"), ("Mosek", "MosekSolver"),
          ("GLPKMathProgInterface", "GLPKSolverMIP")]

for (dir, solver) in solvers
  if isdir(Pkg.dir(dir)) && DEFAULT_SOLVER == nothing
    eval(parse("using "*dir))
    eval(parse("set_default_solver("*solver*"())"))
  end
end

if get_default_solver() == nothing
  packages = ""
  for (dir, solver) in solvers
    packages = packages*dir*" | "
  end
  warn("***********************************************************************************************
       You don't have any of
       "*packages*" installed.
       You must have at least one of these solvers. You can install a solver such as SCS by running:
       Pkg.add(\"SCS\")
       You will have to restart Julia after that.
       ***********************************************************************************************")
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

# Collection of sources required to build libpng
sources = [
    "http://www.imagemagick.org/download/delegates/libpng-1.6.31.tar.gz" =>
    "042e8701abc737a72f1393d0d8fca8ca86146460fe94c6c396261868f489c2de",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd libpng-1.6.31/
mkdir build
cd build/
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain ..
make -j${ncore}
make install

if [[ ${target} == *-w64-mingw* ]]; then
    mkdir ${WORKSPACE}/tmp
    mv $prefix/lib/libpng.dll.a $prefix/bin # fix broken symlink
    cp -r -L $prefix/* ${WORKSPACE}/tmp
    rm -r $prefix
    mv ${WORKSPACE}/tmp $prefix
fi

exit

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libpng16", :libpng)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/bicycle1885/ZlibBuilder/releases/download/v1.0.2/build_Zlib.v1.2.11.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, "libpng", v"1.6.31", sources, script, platforms, products, dependencies)


text:
"""
    jmst(g::JGraph, mst_func)

Function used to animate the appearance of a minimum spannin tree
on a `JGraph`.

# Arguments
- `g::JGraph` the `JGraph` to use for the animation.
- `mst_func` a function that takes an `AbstractGraph` as an arguments and
return a list of `edge`s as defined in `LightGraphs`. One can use `kruskal_mst`
and `prim_mst` defined in `LightGraphs`.
"""

function jmst(g::JGraph, mst_func)
    mst = mst_func(g.data.graph)
    Jnodes = jnodes(g)
    Jedges = jedges(g)
    step_length = g.data.frames[end] ÷ length(mst)
    color_anim =
        Animation([0.0, 1.0], [Lab(colorant"black"), Lab(colorant"red")], [sineio()])
    acted_nodes = []
    for (idx, e) in enumerate(mst)
        from = (idx - 1) * step_length + 1
        to = idx * step_length

        !(e.src in acted_nodes) &&
            act!(Jnodes[e.src], Action(from:to, color_anim, sethue()))
        push!(acted_nodes, e.src)

        !(e.dst in acted_nodes) &&
            act!(Jnodes[e.dst], Action(from:to, color_anim, sethue()))
        push!(acted_nodes, e.dst)

        if haskey(Jedges, e.src => e.dst)
            act!(Jedges[e.src => e.dst], Action(from:to, color_anim, sethue()))
        else
            act!(Jedges[e.dst => e.src], Action(from:to, color_anim, sethue()))
        end
    end
end


text:
using Oceananigans.Operators: ℑzᵃᵃᶠ
using Oceananigans.Coriolis

####
#### Element-wise forcing and right-hand-side calculations
####

@inline ρu_slow_source_term(i, j, k, grid, coriolis, closure, ρ, ρũ, K̃) =
    (- x_f_cross_U(i, j, k, grid, coriolis, ρũ)
     + ∂ⱼτ₁ⱼ(i, j, k, grid, closure, ρ, ρũ, K̃))

@inline ρv_slow_source_term(i, j, k, grid, coriolis, closure, ρ, ρũ, K̃) =
    (- y_f_cross_U(i, j, k, grid, coriolis, ρũ)
     + ∂ⱼτ₂ⱼ(i, j, k, grid, closure, ρ, ρũ, K̃))

@inline ρw_slow_source_term(i, j, k, grid, coriolis, closure, ρ, ρũ, K̃) =
    (- z_f_cross_U(i, j, k, grid, coriolis, ρũ)
     + ∂ⱼτ₃ⱼ(i, j, k, grid, closure, ρ, ρũ, K̃))

@inline ρc_slow_source_term(i, j, k, grid, closure, tracer_index, ρ, ρc, K̃) =
    ∂ⱼDᶜⱼ(i, j, k, grid, closure, tracer_index, ρ, ρc, K̃)

@inline ρt_slow_source_term(i, j, k, grid, closure, tvar::Energy, gases, gravity, ρ, ρũ, ρc̃, K̃) =
    ∂ⱼDᵖⱼ(i, j, k, grid, closure, 1, diagnose_p_over_ρ, tvar, gases, gravity, ρ, ρũ, ρc̃, K̃)

@inline function ρt_slow_source_term(i, j, k, grid, closure, tvar::Entropy, gases, gravity, ρ, ρũ, ρc̃, K̃)
    @inbounds begin
        Ṡ = 0.0
        for gas_index = 1:length(gases)
            tracer_index = gas_index + 1
            ρc = ρc̃[tracer_index]
            Ṡ += ∂ⱼtᶜDᶜⱼ(i, j, k, grid, closure, diagnose_ρs, tracer_index, tvar, gases, gravity, ρ, ρũ, ρc̃, ρc)
        end
        T = diagnose_temperature(i, j, k, grid, tvar, gases, gravity, ρ, ρũ, ρc̃)
        Ṡ += Q_dissipation(i, j, k, grid, closure, ρ, ρũ) / T
        return Ṡ
    end
end

@inline function ρu_fast_source_term(i, j, k, grid, tvar, gases, gravity, advection_scheme, ρ, ρũ, ρc̃, FU)
    @inbounds begin
        return (- div_ρuũ(i, j, k, grid, advection_scheme, ρ, ρũ)
                - ∂p∂x(i, j, k, grid, tvar, gases, gravity, ρ, ρũ, ρc̃)
                + FU[i, j, k])
    end
end

@inline function ρv_fast_source_term(i, j, k, grid, tvar, gases, gravity, advection_scheme, ρ, ρũ, ρc̃, FV)
    @inbounds begin
        return (- div_ρvũ(i, j, k, grid, advection_scheme, ρ, ρũ)
                - ∂p∂y(i, j, k, grid, tvar, gases, gravity, ρ, ρũ, ρc̃)
                + FV[i, j, k])
    end
end

@inline function ρw_fast_source_term(i, j, k, grid, tvar, gases, gravity, advection_scheme, ρ, ρũ, ρc̃, FW)
    @inbounds begin
        return (- div_ρwũ(i, j, k, grid, advection_scheme, ρ, ρũ)
                - ∂p∂z(i, j, k, grid, tvar, gases, gravity, ρ, ρũ, ρc̃)
                - gravity * ℑzᵃᵃᶠ(i, j, k, grid, ρ)
                + FW[i, j, k])
    end
end

@inline function ρc_fast_source_term(i, j, k, grid, advection_scheme, ρ, ρũ, ρc, FC)
    @inbounds begin
        return -div_ρUc(i, j, k, grid, advection_scheme, ρ, ρũ, ρc) + FC[i, j, k]
    end
end

@inline ρt_fast_source_term(i, j, k, grid::AbstractGrid{T}, tvar::Entropy, gases, gravity, ρ, ρũ, ρc̃) where T = zero(T)

@inline ρt_fast_source_term(i, j, k, grid, tvar::Energy, gases, gravity, ρ, ρũ, ρc̃) =
    -∂ⱼpuⱼ(i, j, k, grid, diagnose_pressure, tvar, gases, gravity, ρ, ρũ, ρc̃)


text:
module VarianceComponentTest

import Distributions

export
       vctestnullsim,
       vctest

include("vctestnullsim.jl")
include("vctest.jl")

end


text:
function int2str(y,vocab)
    y = vec(y)
    ysos = findnext(w->!isequal(w,vocab.eos), y, 1)
    ysos == nothing && return ""
    yeos = something(findnext(isequal(vocab.eos), y, ysos), 1+length(y))
    join(vocab.i2w[y[ysos:yeos-1]], " ")
end

max_stars_repo_name: denizyuret/research

max_issues_repo_name: denizyuret/research

max_forks_repo_name: denizyuret/research


text:
module Exchange

using CompEcon

include("states.jl")
include("util.jl")
using .EDSUtil

# ---------- #
# Agent type #
# ---------- #

immutable Agent
    ρ::Float64
    β::Float64
    α::Float64
    η::Float64
    ν::Float64
    ω::Float64
    σ::Float64
    δ::Float64
end

_unpack(a::Agent) = (a.ρ, a.β, a.α, a.η, a.ν, a.ω, a.σ, a.δ)

# cobb douglass (σ = 0)
# _agg(a::Agent, x, y) = x.^(1-a.ω) .* y.^(a.ω)

# CES
function _agg(agent::Agent, x, y)
    ω, σ = agent.ω, agent.σ
    ((1-ω)*x.^σ + ω*y.^σ).^(1/σ)
end

# ---------- #
# Model type #
# ---------- #

"""
Stores the parameters for the following model:

- 2 Agents with EZ time aggregators:
```tex
U_j = V[c_j, μ_t(U_j')] =[(1-β)c_j^ρ+ β μ_t(U'_j)^ρ]^{1/ρ}
```
- Certainty equivalents:
```tex
μ_t(x) = (E_t[x^α])^{1/α}
```
- Production in coungry j:
```tex
y_j = f(k_j, z_j) = [(1-η)k_j^ν + η z_j^{ν}]^{1/ν}
```
- Aggregators of composite goods
```tex
c_j + i_j = h(a_j, b_j) = [(1-ω)a_j^{σ} + ω b_j^{σ}]^{1/σ}
```
- Resource constarints
```tex
a_j + b_j = y_j = f(k_j, z_j)
```
- Capital law of motion
```tex
k'_j = (1-δ) k_j + i_j = (1-δ) k_j + h(a_j, b_j) - c_j
```
- Shocks
```tex
log z_{1,t+1} = log g + γ log z_{1,t} + (1-γ) log z_{2,t} + ζ_1 w_{1,t+1}
log z_{2,t+1} = log g + γ log z_{2,t} + (1-γ) log z_{1,t} + ζ_2 w_{2,t+1}
w_{1,t+1}, w_{2,t+1} ∼ N(0, 1)
```

Where the parameters (ρ, β, α, η, ν, ω, σ, δ, v) are all potentially different
for each agent (notice that γ is crucially the same across agents)
"""
type BCFL22C
    agent1::Agent
    agent2::Agent
    γ::Float64
    ϵ::Matrix{Float64}
    Π::Vector{Float64}
    ζ1::Float64
    ζ2::Float64
end

# TODO: check the sigma parameters
function BCFL22C(;ρ1::Real=-1.0 , α1::Real=-9.0, β1::Real=0.999,  # EZ 1
                  η1::Real=1.0, ν1::Real=0.1,                     # production 1
                  σ1::Real=0.6, ω1::Real=0.1,                     # composite 1
                  δ1::Real=0.025, ζ1::Real=0.001,                 # other 1
                  ρ2::Real=-1.0 , α2::Real=-9.0, β2::Real=0.999,  # EZ 2
                  η2::Real=1.0, ν2::Real=0.1,                     # production 2
                  σ2::Real=0.6, ω2::Real=0.1,                     # composite 2
                  δ2::Real=0.025, ζ2::Real=0.001,                 # other 2
                  γ::Real=0.9, nϵ::Int=4)                         # exog
    agent1 = Agent(ρ1, β1, α1, η1, ν1, ω1, σ1, δ1)
    agent2 = Agent(ρ2, β2, α2, η2, ν2, ω2, σ2, δ2)
    ϵ, Π = EDSUtil.qnwgh(nϵ, 2)
    BCFL22C(agent1, agent2, γ, ϵ, Π, ζ1, ζ2)
end

Base.writemime(io::IO, ::MIME"text/plain", m::BCFL22C) =
    println(io, "BCFL model with 2 goods and 2 agents")

@inline function exog_step(m::BCFL22C, lzbar, ϵ1, ϵ2)
    lzbarp = (2*m.γ - 1)*lzbar .+ (m.ζ1*ϵ1 - m.ζ2*ϵ2)
    lgp = (1-m.γ)*lzbar .+ m.ζ1*ϵ1
    lzbarp, lgp
end

function simulate_exog(m::BCFL22C, capT::Int=10000, seed::Int=1)
    # set random seed, then hand of to real method
    srand(seed)
    simulate_exog(m, randn(capT), randn(capT))
end

# This is `simulate_exog` called in IRF
function simulate_exog(m::BCFL22C, ϵ1::Vector, ϵ2::Vector)

    capT = length(ϵ1)
    capT == length(ϵ2) || error("ϵ1 and ϵ2 should be the same length")

    # Allocate Space
    lzbar = zeros(capT)
    lg = zeros(capT)

    @inbounds for t=2:capT
        lzbar[t], lg[t] = exog_step(m, lzbar[t-1], ϵ1[t], ϵ2[t])
    end

    return lzbar, lg
end

# ------------------------------ #
# Routines to compute allocation #
# ------------------------------ #
"""
If we combine all 4 FOC and allow all params to be different, we get

```
((1-ω1)*(1-ω2))/(ω1*ω2) = (1-sa)^(σ2-1)/(sa^(σ1-1)) * (sb)^(σ1-1)/((1-sb)^(σ2-1))
```

Given a guess for sa = a1, we can solve the above for b2.

Given a1, b2, zbar=z2/z1 we have

```
a2 = 1 - a1
b1 = zbar - b2
```

The above follows from the fact that in the exchange economy scaled by z1 we
have y1 = 1 and y2 = zbar

Then we can compute c1, c2

Given c1, c2 a1, a2 and stuff we can compute the residual

```
♠ - (c2/a2)^(σ2-1.0) * (c1/a1)^(1.0-σ1) * (1.0-ω1)/ω2
```

or

```
♠ - (c2/b2)^(σ2-1.0) * (c1/b1)^(1.0-σ1) * ω1/(1.0-ω2)
```

The use of either residual is accceptable. Given how we solved for b2, they are
equivalent. However, when a1 is very small, the b2 can be very close to zbar,
making it numerically unstable to divide by b2. So, given that we directly
control a1 and a2, we prefer the first equation above
"""
function a1_resid{T<:Real}(m::BCFL22C, st::TimeTState{T}, a1)
    # simplify notation
    ρ1, β1, α1, η1, ν1, ω1, σ1, δ1 = _unpack(m.agent1)
    ρ2, β2, α2, η2, ν2, ω2, σ2, δ2 = _unpack(m.agent2)
    l♠t, lzbar = st.l♠, st.lzbar

    a2 = 1.0 - a1

    zbar = exp(lzbar)

    if abs(σ1 - σ2) < 1e-14  # sigmas are the same, used closed form sb
        y = ((1-ω1)*(1-ω2)/(ω1*ω2))^(1/(σ1-1))*(a1/a2)
        b2 = zbar / (1 + y)
    else
        error("implement solver for sb")
    end

    # get rest of allocation
    b1 = zbar - b2
    c1 = ((1-ω1)*a1^σ1 + ω1*b1^σ1)^(1/σ1)
    c2 = ((1-ω2)*b2^σ2 + ω2*a2^σ2)^(1/σ2)

    # compute RHS of residual for the a_1 FOC
    rhs = (c2/a2)^(σ2-1.0) * (c1/a1)^(1.0-σ1) * (1.0-ω1)/ω2
    # rhs2 = (c2/b2)^(σ2-1.0) * (c1/b1)^(1.0-σ1) * ω1/(1.0-ω2)
    # if abs(rhs - rhs2) > 1e-12
    #     warn("rhs are different for statei l♠t = $(l♠t), zbar=$zbar, a1=$a1")
    #     @show rhs, rhs2
    # end

    # return the residual as first argument and allocation as second
    log(rhs) - l♠t, (a1, a2, b1, b2, c1, c2)
end

function get_allocation{T<:Real}(m::BCFL22C, st::TimeTState{T})
    a1t = brent(foo->a1_resid(m, st, foo)[1], 1e-15, 1-1e-15)

    # evaluate once more to get the allocation at the optimal a1
    a1t, a2t, b1t, b2t, c1t, c2t = a1_resid(m, st, a1t)[2]
end

get_allocation{T<:Real}(m::BCFL22C, fst::FullState{T}) =
    get_allocation(m, TimeTState(fst.l♠, fst.lzbar))

# ------------------- #
# PolicyFunction type #
# ------------------- #
immutable PolicyFunction{D}  # type parameter is complete polynomial degree
    coefs::Vector{Float64}
end

evaluate{T<:Real}(pf::PolicyFunction, fst::FullState{T}) =
    evaluate(pf, Matrix(fst))[1]  # extract first (only) element

evaluate{T<:Real}(pf::PolicyFunction, fst::FullState{Vector{T}}) =
    evaluate(pf, Matrix(fst))

function evaluate{D}(pf::PolicyFunction{D}, state::Vector)
    X = complete_polynomial(reshape(state, 1, length(state)), D)
    y = X*pf.coefs
    return y[1]  # return only element of our output
end

function evaluate{D}(pf::PolicyFunction{D}, state::Matrix)
    X = complete_polynomial(state, D)
    y = X*pf.coefs
    return y
end

function evaluate!{D}(pf::PolicyFunction{D}, state::Matrix, X_buf::Matrix,
                      y_buf::Vector)
    complete_polynomial!(state, D, X_buf)  # fill in X_buf
    A_mul_B!(y_buf, X_buf, pf.coefs)       # fill in y_buf = X_buf * coefs
    y_buf                                  # return y_buf
end

# ----------------- #
# Updating the grid #
# ----------------- #
"""
Given scalars for ♠ and zbar, return vectors of ♠' and zbar' for every
ϵ' in m.ϵ
"""
function get_next_grid{T<:Real}(m::BCFL22C, pf::PolicyFunction,
                                state::TimeTState{T})
    # unpack
    ϵ1, ϵ2 = m.ϵ[:, 1], m.ϵ[:, 2]
    lzbarp, lgp = exog_step(m, state.lzbar, ϵ1, ϵ2)

    # build state and use pf to step ♠ forward
    fst = FullState(state, lzbarp, lgp)
    l♠p = evaluate(pf, fst)

    return l♠p, lzbarp, lgp
end

"""
Given vectors for ♠ and zbar, return matrices of ♠' and zbar' for every
ϵ' in m.ϵ.

The layout of the matrix is such that

```
♠[i, j] = ♠'(♠_i, ϵ_j),
```

which means

```
E[♠'] = matrix_from_this_func * m.Π
```
"""
function get_next_grid{T<:Real}(m::BCFL22C, pf::PolicyFunction,
                                states::TimeTState{Vector{T}})
    # extract dimensions
    N = length(states)
    J = size(m.ϵ, 1)

    # allocte memory. Create transpose. Fill in columns, transpose before return
    l♠p = Array(Float64, J, N)
    lzbarp = similar(l♠p)
    lgp = similar(l♠p)

    # fill in row by row by calling one state function above
    for i=1:N
        l♠p[:, i], lzbarp[:, i], lgp[:, i] = get_next_grid(m, pf, states[i])
    end

    # transpose before returning
    return l♠p', lzbarp', lgp'

end

# ------------------- #
# ValueFunctions type #
# ------------------- #
immutable ValueFunctions{D}
    coefs::Matrix{Float64}  # first column U_1, second col U_2
end

# return a Vector of [U_1(l♠, zbar), U_2(l♠, lzbar)]. Reshape into vector
value_funcs{D,T<:Real}(vfs::ValueFunctions{D}, st::TimeTState{T}) =
    reshape(complete_polynomial(Matrix(st), D) * vfs.coefs, 2)

value_funcs{D,T<:Real}(vfs::ValueFunctions{D}, st::TimeTState{Vector{T}}) =
    complete_polynomial(Matrix(st), D) * vfs.coefs

value_funcs(vfs::ValueFunctions, fst::FullState) =
    value_funcs(vfs, TimeTState(fst))

"""
This method assumes you are passing the values of (l♠, zbar) in period t.
It will construct the implied (l♠', zbar', g') for all integration nodes in
m.ϵ and then hand off to the function below
"""
function certainty_equiv{T<:Real}(m::BCFL22C, pf::PolicyFunction,
                                  vfs::ValueFunctions, st::TimeTState{T})
    l♠p, lzbarp, lgp = get_next_grid(m, pf, st)
    certainty_equiv(m, vfs, l♠p, lzbarp, lgp)
end

"""
This method assumes you are passing (l♠', zbar', g') for all the integration
nodes in m.ϵ
"""
function certainty_equiv(m::BCFL22C, vfs::ValueFunctions, l♠p::Vector,
                         lzbarp::Vector, lgp::Vector)
    J = length(m.Π)
    if J != length(l♠p) || J != length(lzbarp) || J != length(lgp)
        error("l♠p, lzbarp and lgp should have $J elements")
    end

    # Compute value functions and evaluate expectations that form certainty
    # equiv
    VF = value_funcs(vfs, TimeTState(l♠p, lzbarp))
    μ1 = dot(m.Π, (exp(lgp).*VF[:, 1]).^(m.agent1.α))^(1.0/m.agent1.α)
    μ2 = dot(m.Π, (exp(lgp).*VF[:, 2]).^(m.agent2.α))^(1.0/m.agent2.α)
    μ1, μ2
end

# TODO: add a "vectorized" version of certainty_equiv that evaluates
#       at mulitple time t states at once. Will help in the euler equation func

# ----------------------- #
# VFI for post simulation #
# ----------------------- #
function vfi_from_simulation(m::BCFL22C, pf::PolicyFunction, grid::Matrix{Float64};
                             deg::Int=3, tol=1e-8, maxit::Int=5000)
    # unpack
    ρ1, β1, α1, η1, ν1, ω1, σ1, δ1 = _unpack(m.agent1)
    ρ2, β2, α2, η2, ν2, ω2, σ2, δ2 = _unpack(m.agent2)
    st = TimeTState(grid[:, 1], grid[:, 2])
    Ngrid = length(st)

    # get consumption on the grid
    cs = [get_allocation(m, st_t)[5:6] for st_t in st]
    c1 = map(x->getindex(x, 1), cs)
    c2 = map(x->getindex(x, 2), cs)

    # construct _constant_ basis matrix
    basis_mat = complete_polynomial(grid, deg)

    # get coefficients for degree deg complete poly in ♠ and zbar for the
    # initial guess that U_i = c_i
    coefs = basis_mat \ [c1 c2]
    vfs = ValueFunctions{deg}(coefs)

    # advance grid one time step using integration nodes m.ϵ
    l♠p, lzbarp, lgp = get_next_grid(m, pf, st)

    # iteration manager stuff + allocate for value functions
    err = 10.0
    it = 0
    U1_old, U2_old = copy(c1), copy(c2)
    U1, U2 = similar(c1), similar(c2)

    # TODO: this slices rows (yuck)
    while err > tol && it < maxit
        it += 1
        for i=1:Ngrid
            # compute certainty equivalents
            μ1, μ2 = certainty_equiv(m, vfs, l♠p[i, :][:], lzbarp[i, :][:],
                                     lgp[i, :][:])

            # apply backward step to get U_{i,t}
            U1[i] = ((1-β1)*c1[i]^ρ1 + β1*μ1^ρ1)^(1/ρ1)
            U2[i] = ((1-β2)*c2[i]^ρ2 + β2*μ2^ρ2)^(1/ρ2)
        end

        # compute residual
        err = max(maxabs(U1 - U1_old), maxabs(U2 - U2_old))

        # update coefficient vector
        coefs = basis_mat \ [U1 U2]
        vfs = ValueFunctions{deg}(coefs)

        # update cache of previous U_i
        copy!(U1_old, U1)
        copy!(U2_old, U2)
    end

    vfs
end

# ------------------------- #
# Main simulation algorithm #
# ------------------------- #

#=
now construct LHS of regression using the euler equation:

♠' = ♠ * beta_1/beta_2 * g'^(α1-α2) * (c_1'/c_1)^(ρ_1-1) *
     (c_2'/c_2)^(1-ρ_2) * (U_1'/μ_1)^(α_1 - ρ_1) *(U_2'/μ_2)^(ρ_2 - α_2)
=#
function eval_euler_eqn!(m::BCFL22C, pf::PolicyFunction, vfs::ValueFunctions,
                         fsts::FullState, c1::Vector, c2::Vector, LHS::Vector)
    ρ1, β1, α1, η1, ν1, ω1, σ1, δ1 = _unpack(m.agent1)
    ρ2, β2, α2, η2, ν2, ω2, σ2, δ2 = _unpack(m.agent2)

    for t=1:length(c1)-1  # can only go to capT-1 b/c we need c[t+1]
        # get μ_t based on current state
        μ1, μ2 = certainty_equiv(m, pf, vfs, TimeTState(fsts[t]))

        # get U_{t+1} from simulation
        U1p, U2p = value_funcs(vfs, TimeTState(fsts[t+1]))

        # now package up the LHS of the Euler eqn, which is ♠'
        LHS[t] = (exp(fsts.l♠[t])*(β2/β1) *
                  exp(fsts.lgp[t])^(α2-α1) *
                  (c1[t+1]/c1[t])^(1.0 - ρ1) *
                  (c2[t+1]/c2[t])^(ρ2 - 1.0) *
                  (U1p/μ1)^(ρ1 - α1) *
                  (U2p/μ2)^(α2 - ρ2))
    end

    # take logs of LHS so it is log♠'
    map!(log, LHS)
end

"""
Simulate ♠ forward and solve for c1, c2 along the way. Need the cs so I can
evaluate the euler equation later
"""
function do_simulation!(m::BCFL22C, pf::PolicyFunction, fsts::FullState,
                        c1::Vector, c2::Vector)
    for t=1:length(fsts) - 1
        # step l♠ forward one period
        fst = fsts[t]
        l♠p = evaluate(pf, fst)
        fsts.l♠[t+1] = l♠p

        # Solve for the optimal allocation at this ♠, zbart and store
        # consumption in vectors
        c1[t], c2[t] = get_allocation(m, fst)[5:6]
    end
end

function linear_coefs(m::BCFL22C, lzbar::Vector{Float64}=simulate_exog(m)[1],
                      lg::Vector{Float64}=simulate_exog(m)[2],
                      pf::PolicyFunction{1}=PolicyFunction{1}([0.0, 0.95, -0.95, 0.0, 0.0]);
                      maxiter::Int=10_000, grid_skip::Int=5,
                      verbose::Bool=true, print_skip::Int=1,
                      ξ::Float64=0.05)
    # Make sure we are working with the endowment economy
    if abs(m.agent1.η-1.0) > 1e-14 || abs(m.agent2.η - 1.0) > 1e-14
        msg = "Solution only implemented for endowment economy (η1=η2=1.0)"
        error(msg)
    end

    # unpack
    ρ1, β1, α1, η1, ν1, ω1, σ1, δ1 = _unpack(m.agent1)
    ρ2, β2, α2, η2, ν2, ω2, σ2, δ2 = _unpack(m.agent2)
    capT = length(lzbar) - 1

    # unpack exog
    γ = m.γ
    ϵ1 = m.ϵ[:, 1]
    ϵ2 = m.ϵ[:, 2]
    Π = m.Π

    # Iteration management stuff
    tol = 1e-4*ξ
    dist = 100.0
    it = 0

    # ♠ is ratio multipliers: φ_2/φ_1. Starts at 1.0, so logs at 0
    l♠ = zeros(capT+1)
    LHS = Array(Float64, capT-1)  # LHS of regression from euler equation

    l♠_old = fill(100.0, capT+1)

    c1 = Array(Float64, capT)
    c2 = Array(Float64, capT)

    # costruct state vector.
    # NOTE: I pad the last elements of time t+1 exog states with a NaN because
    #       we need the full capT+1 length time series of the time t states
    #       in order to fill in the allocation during the simulation.
    fsts = FullState(l♠, lzbar, [lzbar[2:capT+1]; NaN], [lg[2:capT+1]; NaN])

    # core algorithm
    while dist > tol && it < maxiter
        start_time = time()

        # do simulation. Updates ♠, c1, c2
        do_simulation!(m, pf, fsts, c1, c2)

        # build EDS grid on ♠ and lzbar terms in regression matrix
        if mod(it, grid_skip) == 0 || it == 0
            verbose && info("Updating grid")
            grid = get_eds_grid(Matrix(fsts)[1:capT-1, 1:2])
        end

        # do VFI to get approximation to value functions/certainty equiv
        vfs = vfi_from_simulation(m, pf, grid; deg=3, tol=1e-8, maxit=5000)

        # update LHS
        eval_euler_eqn!(m, pf, vfs, fsts, c1, c2, LHS)

        # compute dist, copy cache
        dist = mean(abs(1.0 - exp(fsts.l♠ - l♠_old)))
        copy!(l♠_old, fsts.l♠)

        # update policy function
        X = complete_polynomial(Matrix(fsts)[1:capT-1, :], 1)
        κ_hat = X \ LHS
        κ = ξ*κ_hat + (1-ξ)*pf.coefs
        pf = PolicyFunction{1}(κ)

        it += 1
        if verbose && mod(it, print_skip) == 0
            tot_time = time() - start_time
            @printf "Iteration %i, dist %2.5e, time %5.5e\n" it dist tot_time
        end
    end

    fsts, pf

end

function main()
    m = BCFL22C()
    lzbar, lg = simulate_exog(m)
    κ = [-4.594267427095823e-7,0.9999976514132554,1.8822879898974099,-1.8874374369713394,2.9809481869877624e-7]
    pf = PolicyFunction{1}(κ)
    ξ = 0.05
    fsts, pf = linear_coefs(m, lzbar, lg, pf, maxiter=500)
    m, fsts, pf
end

include("analysis.jl")

end  # module

#=
function foobar(m::BCFL22C, a1::Real, zbar::Real)
    a2 = 1.0 - a1

    if abs(σ1 - σ2) < 1e-14  # sigmas are the same, used closed form sb
        y = ((1-ω1)*(1-ω2)/(ω1*ω2))^(1/(σ1-1))*(a1/a2)
        b2 = zbar / (1 + y)

        # cobb douglass
        # x = (1-ω1)*(1-ω2)/(ω1*ω2)*(a2/a1)
        # b2 = zbar*x/(1+x)
    else
        error("implement solver for sb")
    end

    # get rest of allocation
    b1 = zbar - b2
    c1 = ((1-ω1)*a1^σ1 + ω1*b1^σ1)^(1/σ1)
    c2 = ((1-ω2)*b2^σ2 + ω2*a2^σ2)^(1/σ2)

    ♠ = (c2/a2)^(1.0-σ2) * (c1/a1)^(σ1-1.0) * ω2/(1-ω1)
    ♠2 = (c2/b2)^(1.0-σ2) * (c1/b1)^(σ1-1.0) * (1-ω2)/ω1

    @show ♠, ♠2
    return ♠

end

function foobar(m::BCFL22C, a1_grid::Vector, zbar_grid::Vector)
    grid = gridmake(a1_grid, zbar_grid)
end
=#

max_stars_repo_name: NYUEcon/GrowthModels

max_issues_repo_name: NYUEcon/GrowthModels

max_forks_repo_name: NYUEcon/GrowthModels


text:
using DependencyTrees: deps, leftdeps, rightdeps, leftmostdep, rightmostdep

@testset "Graphs" begin

    sent = [
        ("Economic", "ATT", 2),
        ("news", "SBJ", 3),
        ("had", "PRED", 0),
        ("little", "ATT", 5),
        ("effect", "OBJ", 3),
        ("on", "ATT", 5),
        ("financial", "ATT", 8),
        ("markets", "PC", 6),
        (".", "PU", 3),
    ]

    graph = deptree(sent) do ((form, deprel, head))
        DependencyTrees.Token(form, head, deprel)
    end
    @test graph == deepcopy(graph)

    @test [t.form for t in graph.tokens] == [first(t) for t in sent]
    for (i, t) in enumerate(sent)
        @test graph.tokens[i].form == first(t)
    end

    @test DependencyTrees.to_conllu(graph) |> strip == """
1	Economic	_	_	_	_	2	ATT	_	_
2	news	_	_	_	_	3	SBJ	_	_
3	had	_	_	_	_	0	PRED	_	_
4	little	_	_	_	_	5	ATT	_	_
5	effect	_	_	_	_	3	OBJ	_	_
6	on	_	_	_	_	5	ATT	_	_
7	financial	_	_	_	_	8	ATT	_	_
8	markets	_	_	_	_	6	PC	_	_
9	.	_	_	_	_	3	PU	_	_
""" |> strip

    @test showstr(graph) ==
        """
        ┌────────────── ROOT
        │           ┌─► Economic
        │        ┌─►└── news
        └─►┌──┌──└───── had
           │  │     ┌─► little
           │  └─►┌──└── effect
           │  ┌──└────► on
           │  │     ┌─► financial
           │  └────►└── markets
           └──────────► .
        """ |> strip

    @test length(graph) == length(sent) == 9
    for (i, (word, tag, id_)) in enumerate(sent)
        sent_deps = filter(x -> x[3] == i, sent)
        deps_ = deps(graph, i)
        @test length(deps_) == length(sent_deps)
        @test Set(graph.tokens[id].form for id in deps_) == Set([d[1] for d in sent_deps])
    end

    @test leftdeps(graph, 0) == []
    @test rightdeps(graph, 0) == [3]
    @test leftdeps(graph, 1) == []
    @test rightdeps(graph, 1) == []
    @test leftdeps(graph, 2) == [1]
    @test rightdeps(graph, 2) == []
    @test leftdeps(graph, 3) == [2]
    @test rightdeps(graph, 3) == [5,9]
    @test leftdeps(graph, 4) == []
    @test rightdeps(graph, 4) == []
    @test leftdeps(graph, 5) == [4]
    @test rightdeps(graph, 5) == [6]
    @test leftdeps(graph, 6) == []
    @test rightdeps(graph, 6) == [8]
    @test leftdeps(graph, 7) == []
    @test rightdeps(graph, 7) == []
    @test leftdeps(graph, 8) == [7]
    @test rightdeps(graph, 8) == []
    @test leftdeps(graph, 9) == []
    @test rightdeps(graph, 9) == []

    @test leftmostdep(graph, 1) == -1
    @test rightmostdep(graph, 1) == -1
    @test leftmostdep(graph, 2) == 1
    @test rightmostdep(graph, 2) == -1
    @test leftmostdep(graph, 3) == 2
    @test rightmostdep(graph, 3) == 9
    @test leftmostdep(graph, 4) == -1
    @test rightmostdep(graph, 4) == -1
    @test leftmostdep(graph, 5) == 4
    @test rightmostdep(graph, 5) == 6
    @test leftmostdep(graph, 6) == -1
    @test rightmostdep(graph, 6) == 8
    @test leftmostdep(graph, 7) == -1
    @test rightmostdep(graph, 7) == -1
    @test leftmostdep(graph, 8) == 7
    @test rightmostdep(graph, 8) == -1
    @test leftmostdep(graph, 9) == -1
    @test rightmostdep(graph, 9) == -1

    sent = [
        ("Pierre", "NNP", 2),
        ("Vinken", "NNP", 8),
        (",", ",", 2),
        ("61", "CD", 5),
        ("years", "NNS", 6),
        ("old", "JJ", 2),
        (",", ",", 2),
        ("will", "MD", 0),
        ("join", "VB", 8),
        ("the", "DT", 11),
        ("board", "NN", 9),
        ("as", "IN", 9),
        ("a", "DT", 15),
        ("nonexecutive", "JJ", 15),
        ("director", "NN", 12),
        ("Nov.", "NNP", 9),
        ("29", "CD", 16),
        (".", ".", 8)
    ]

    graph = deptree(sent) do ((form, label, head))
        DependencyTrees.Token(form, head, label)
    end

    @test length(graph) == length(sent) == 18
    # @test isroot(graph[0])
    for i in 1:length(sent)
        sent_deps = filter(x -> x[3] == i, sent)
        deps_ = deps(graph, i)
        @test length(deps_) == length(sent_deps)
        @test Set(graph[id].form for id in deps_) == Set([d[1] for d in sent_deps])
    end
    # @test isprojective(graph)

    @testset "Projectivity" begin

        # mcdonald & pereira 2005
        # https://www.seas.upenn.edu/~strctlrn/bib/PDF/nonprojectiveHLT-EMNLP2005.pdf
        sent = [
            ("john", 2),      # 1
            ("saw", 0),       # 2
            ("a", 4),         # 3
            ("dog", 2),       # 4
            ("yesterday", 2), # 5
            ("which", 7),     # 6
            ("was", 4),       # 7
            ("a", 9),         # 8
            ("yorkshire", 10),# 9
            ("terrier", 7)    # 10
        ]
        # graph = DependencyTree([Token(fm, hd, id=i) for (i, (fm, hd)) in enumerate(sent)], 2)
        graph = deptree(t -> DependencyTrees.Token(t...), sent)
        # @test !isprojective(graph)

        @test DependencyTrees.to_conllu(graph) |> strip == """
1	john	_	_	_	_	2	_	_	_
2	saw	_	_	_	_	0	_	_	_
3	a	_	_	_	_	4	_	_	_
4	dog	_	_	_	_	2	_	_	_
5	yesterday	_	_	_	_	2	_	_	_
6	which	_	_	_	_	7	_	_	_
7	was	_	_	_	_	4	_	_	_
8	a	_	_	_	_	9	_	_	_
9	yorkshire	_	_	_	_	10	_	_	_
10	terrier	_	_	_	_	7	_	_	_
""" |> strip

        @test showstr(graph) |> strip ==
            """
               ┌─────────── ROOT
               │        ┌─► john
               └─►┌──┌──└── saw
                  │  │  ┌─► a
               ┌──│  └─►└── dog
               │  └───────► yesterday
               │        ┌─► which
            ┌──└───────►└── was
            │           ┌─► a
            │        ┌─►└── yorkshire
            └───────►└───── terrier
            """ |> strip

        # jurafsky & martin, speech & language processing (3ed)
        sent = [
            ("jetblue", "nsubj", 2), # 1
            ("canceled", "root", 0), # 2
            ("our", "det", 4),       # 3
            ("flight", "dobj", 2),   # 4
            ("this", "det", 2),      # 5
            ("morning", "nmod", 7),  # 6
            ("which", "case", 4),    # 7
            ("was", "mod", 9),       # 8
            ("already", "adv", 10),  # 9
            ("late", "mod", 7)       # 10
        ]
        graph = deptree(sent) do ((form, label, head))
            DependencyTrees.Token(form, head, label)
        end
        @test graph[1].label == "nsubj"
    end

    @testset "Errors" begin

        using DependencyTrees: NonProjectiveGraphError

        sent = [
            ("john", 2),      # 1
            ("saw", 0),       # 2
            ("a", 4),         # 3
            ("dog", 2),       # 4
            ("yesterday", 2), # 5
            ("which", 7),     # 6
            ("was", 4),       # 7
            ("a", 9),         # 8
            ("yorkshire", 10),# 9
            ("terrier", 7)    # 10
        ]
        tree = deptree(x -> Token(x...), sent)

        @test !is_projective(tree)

        for system in (ArcEager(), ArcHybrid(), ArcStandard())
            oracle = Oracle(system, static_oracle)
            sequence = oracle(tree)
            @test sequence isa DependencyTrees.TransitionParsing.UnparsableTree
        end
    end

end


text:
module TestNonFoldable

using FGenerators
using FLoops
using Test

struct OrganPipe  # not inheriting `Foldable`
    n::Int
end

@fgenerator(foldable::OrganPipe) do
    n = foldable.n
    @yieldfrom 1:n
    @yieldfrom n-1:-1:1
end

@testset "n=$n" for (n, desired) in [
    (2, [1, 2, 1]),
    (3, [1, 2, 3, 2, 1]),
    (4, [1:4; 3:-1:1]),
    (5, [1:5; 4:-1:1]),
]
    actual = []
    @floop for x in OrganPipe(n)
        push!(actual, x)
    end
    @test actual == desired
end

end  # module