text:
####################################
#
# Extension of SLiDE model to include
#       dynamics and perfect foresight
#
####################################

using SLiDE
using CSV
using JuMP
using Complementarity
using DataFrames

############
# LOAD DATA
############

#SLiDE data needs to be built or point to pre-existing build directory
#can pass a name (d, set) = build_data("name_of_build_directory")
!(@isdefined(d_in) && @isdefined(set_in)) && ((d_in, set_in) = build("state_model"))
d = copy(d_in)
set = copy(set_in)

#Define range of years in time horizon
# !!!! or load from CSV or YAML
bmkyr = 2016
endyr = 2018
years = bmkyr:endyr
#years = [2017, 2016, 2019, 2018]

#Define years and associated flags
years = ensurearray(sort(years))
yrl = years[length(years)]
yrf = years[1]
islast = Dict(years[k] => (years[k] == yrl ? 1 : 0) for k in keys(years))
isfirst = Dict(years[k] => (years[k] == yrf ? 1 : 0) for k in keys(years))
yrint = Dict(years[k+1] => years[k+1]-years[k] for k in 1:(length(years)-1))

#Load slide data and time horizon to produce model data and appropriate time-indexed subsets
(sld, set, idx) = model_input(years, d, set)

########## Model ##########
cge = MCPModel();

##############
# PARAMETERS
##############

#benchmark values
#benchmark values
@NLparameter(cge, ys0[r in set[:r], s in set[:s], g in set[:g]] == sld[:ys0][r,s,g]); # Sectoral supply
@NLparameter(cge, id0[r in set[:r], s in set[:s], g in set[:g]] == sld[:id0][r,s,g]); # Intermediate demand
@NLparameter(cge, ld0[r in set[:r], s in set[:s]] == sld[:ld0][r,s]); # Labor demand
@NLparameter(cge, kd0[r in set[:r], s in set[:s]] == sld[:kd0][r,s]); # Capital demand
@NLparameter(cge, ty0[r in set[:r], s in set[:s]] == sld[:ty0][r,s]); # Production tax (benchmark)
@NLparameter(cge, m0[r in set[:r], g in set[:g]] == sld[:m0][r,g]); # Imports
@NLparameter(cge, x0[r in set[:r], g in set[:g]] == sld[:x0][r,g]); # Exports of goods and services
@NLparameter(cge, rx0[r in set[:r], g in set[:g]] == sld[:rx0][r,g]); # Re-exports of goods and services
@NLparameter(cge, md0[r in set[:r], m in set[:m], g in set[:g]] == sld[:md0][r,m,g]); # Total margin demand
@NLparameter(cge, nm0[r in set[:r], g in set[:g], m in set[:m]] == sld[:nm0][r,g,m]); # Margin demand from national market
@NLparameter(cge, dm0[r in set[:r], g in set[:g], m in set[:m]] == sld[:dm0][r,g,m]); # Margin supply from local market
@NLparameter(cge, s0[r in set[:r], g in set[:g]] == sld[:s0][r,g]); # Aggregate supply
@NLparameter(cge, a0[r in set[:r], g in set[:g]] == sld[:a0][r,g]); # Armington supply
@NLparameter(cge, ta0[r in set[:r], g in set[:g]] == sld[:ta0][r,g]); # Tax net of subsidy rate on intermediate demand
@NLparameter(cge, tm0[r in set[:r], g in set[:g]] == sld[:tm0][r,g]); # Import tariff
@NLparameter(cge, cd0[r in set[:r], g in set[:g]] == sld[:cd0][r,g]); # Final demand
@NLparameter(cge, c0[r in set[:r]] == sld[:c0][r]); # Aggregate final demand
@NLparameter(cge, yh0[r in set[:r], g in set[:g]] == sld[:yh0][r,g]); #Household production
@NLparameter(cge, bopdef0[r in set[:r]] == sld[:bopdef0][r]); #Balance of payments
@NLparameter(cge, hhadj[r in set[:r]] == sld[:hhadj][r]); # Household adjustment
@NLparameter(cge, g0[r in set[:r], g in set[:g]] == sld[:g0][r,g]); # Government demand
@NLparameter(cge, i0[r in set[:r], g in set[:g]] == sld[:i0][r,g]); # Investment demand
@NLparameter(cge, xn0[r in set[:r], g in set[:g]] == sld[:xn0][r,g]); # Regional supply to national market
@NLparameter(cge, xd0[r in set[:r], g in set[:g]] == sld[:xd0][r,g]); # Regional supply to local market
@NLparameter(cge, dd0[r in set[:r], g in set[:g]] == sld[:dd0][r,g]); # Regional demand to local market
@NLparameter(cge, nd0[r in set[:r], g in set[:g]] == sld[:nd0][r,g]); # Regional demand to national market

#counterfactual taxes
@NLparameter(cge, ty[r in set[:r], s in set[:s]] == sld[:ty0][r,s]); # Output tax
@NLparameter(cge, ta[r in set[:r], g in set[:g]] == sld[:ta0][r,g]); # Armington tax
@NLparameter(cge, tm[r in set[:r], g in set[:g]] == sld[:tm0][r,g]); # Import tariff

# benchmark value share parameters
@NLparameter(cge, alpha_kl[r in set[:r], s in set[:s]] == ensurefinite(value(ld0[r,s]) / (value(ld0[r,s]) + value(kd0[r,s]))));
@NLparameter(cge, alpha_x[r in set[:r], g in set[:g]] == ensurefinite((value(x0[r,g]) - value(rx0[r,g])) / value(s0[r,g])));
@NLparameter(cge, alpha_d[r in set[:r], g in set[:g]] == ensurefinite(value(xd0[r,g]) / value(s0[r,g])));
@NLparameter(cge, alpha_n[r in set[:r], g in set[:g]] == ensurefinite(value(xn0[r,g]) / value(s0[r,g])));
@NLparameter(cge, theta_n[r in set[:r], g in set[:g]] == ensurefinite(value(nd0[r,g]) / (value(nd0[r,g]) - value(dd0[r,g]))));
@NLparameter(cge, theta_m[r in set[:r], g in set[:g]] == ensurefinite((1+value(tm0[r,g])) * value(m0[r,g]) / (value(nd0[r,g]) + value(dd0[r,g]) + (1 + value(tm0[r,g])) * value(m0[r,g]))));

#Substitution and transformation elasticities
# !!!! Currently assigned to global variable in definitions.jl
@NLparameter(cge, es_va[r in set[:r], s in set[:s]] == SUB_ELAST[:va]); # value-added nest - substitution elasticity
@NLparameter(cge, es_y[r in set[:r], s in set[:s]]  == SUB_ELAST[:y]); # Top-level Y nest (VA,M) - substitution elasticity
@NLparameter(cge, es_m[r in set[:r], s in set[:s]]  == SUB_ELAST[:m]); # Materials nest - substitution elasticity
@NLparameter(cge, et_x[r in set[:r], g in set[:g]]    == TRANS_ELAST[:x]); # Disposition, distribute regional supply to local, national, export - transformation elasticity
@NLparameter(cge, es_a[r in set[:r], g in set[:g]]    == SUB_ELAST[:a]); # Top-level A nest for aggregate demand (set[:m], set[:g]) - substitution elasticity
@NLparameter(cge, es_mar[r in set[:r], g in set[:g]]  == SUB_ELAST[:mar]); # Margin supply - substitution elasticity
@NLparameter(cge, es_d[r in set[:r], g in set[:g]]    == SUB_ELAST[:d]); # Domestic demand aggregation nest (intranational) - substitution elasticity
@NLparameter(cge, es_f[r in set[:r], g in set[:g]]    == SUB_ELAST[:f]); # Domestic and foreign demand aggregation nest (international) - substitution elasticity

# Temporal/Dynamic modifications
@NLparameter(cge, ir == 0.05); # Interest rate
@NLparameter(cge, gr == 0.01); # Growth rate
@NLparameter(cge, dr == 0.02); # Depreciation rate

@NLparameter(cge, pvm[yr in years] == (1/(1+value(ir)))^(yr-yrf)); # Reference price path - Present value multiplier
@NLparameter(cge, qvm[yr in years] == (1 + value(gr))^(yr-yrf)); # Reference quantity path multiplier

@NLparameter(cge, rk0 == value(ir) + value(dr)); # Initial return to capital

# !!!! Not used, will likely be used in updated model
# and intertemporal choice
# Initial capital earnings vk0
@NLparameter(cge, vk0[r in set[:r], s in set[:s]] == value(kd0[r,s]));
# Initial capital stock k0
@NLparameter(cge, k0[r in set[:r], s in set[:s]] == value(vk0[r,s])/value(rk0));
# Initial investment inv0
@NLparameter(cge, inv0[r in set[:r], s in set[:s]] == (value(dr)+value(gr))*value(k0[r,s]));
# !!!!

##################
# -- VARIABLES --
##################

# Set lower bound - small value that acts as a lower limit to variable values
# default is zero
lo = MODEL_LOWER_BOUND

# !!!! Define variables separately or in blocks
# sectors
@variable(cge, Y[(yr,r,s) in set[:Y]] >= lo, start = value(qvm[yr])); # Output
@variable(cge, X[(yr,r,g) in set[:X]] >= lo, start = value(qvm[yr])); # Exports
@variable(cge, A[(yr,r,g) in set[:A]] >= lo, start = value(qvm[yr])); # Armington
@variable(cge, C[yr in years, r in set[:r]] >= lo, start = value(qvm[yr])); # Consumption
@variable(cge, MS[yr in years, r in set[:r], m in set[:m]] >= lo, start = value(qvm[yr])); # Margin supply

@variable(cge, K[(yr,r,s) in set[:PK]] >= lo, start = value(qvm[yr])*value(kd0[r,s])); # Capital
@variable(cge, I[(yr,r,s) in set[:PK]] >= lo, start = value(qvm[yr])*(value(dr)+value(gr))*value(kd0[r,s])); # Investment
@variable(cge, TK[(r,s) in set[:PKT]] >= lo, start = value(kd0[r,s]) * value(qvm[yrl])*(1+value(gr))); # Terminal Capital

# !!!! Future updates to capital stock/investment
# @variable(cge, K[(yr,r,s) in set[:PK]] >= lo, start = value(qvm[yr])*value(k0[r,s])); # Capital
# @variable(cge, I[(yr,r,s) in set[:PK]] >= lo, start = value(qvm[yr])*value(inv0[r,s])); # Investment
# @variable(cge, TK[(r,s) in set[:PKT]] >= lo, start = value(k0[r,s]) * value(qvm[yr]) * (1 + value(gr))); # Terminal Capital

# commodities
@variable(cge, PA[(yr,r,g) in set[:PA]] >= lo, start = value(pvm[yr])); # Regional market (input)
@variable(cge, PY[(yr,r,g) in set[:PY]] >= lo, start = value(pvm[yr])); # Regional market (output)
@variable(cge, PD[(yr,r,g) in set[:PD]] >= lo, start = value(pvm[yr])); # Local market price
@variable(cge, PN[yr in years, g in set[:g]] >= lo, start = value(pvm[yr])); # National market
@variable(cge, PL[yr in years, r in set[:r]] >= lo, start = value(pvm[yr])); # Wage rate
@variable(cge, PM[yr in years, r in set[:r], m in set[:m]] >= lo, start = value(pvm[yr])); # Margin price
@variable(cge, PC[yr in years, r in set[:r]] >= lo, start = value(pvm[yr])); # Consumer price index
@variable(cge, PFX[yr in years] >= lo, start = value(pvm[yr])); # Foreign exchange

@variable(cge, PK[(yr,r,s) in set[:PK]] >= lo, start = value(pvm[yr]) * (1 + value(ir))); # Price of capital
@variable(cge, RK[(yr,r,s) in set[:PK]] >= lo, start = value(pvm[yr])*value(rk0)); # Capital rental rate
@variable(cge, PKT[(r,s) in set[:PKT]] >= lo, start = value(pvm[yrl])); # Terminal capital price
#@variable(cge, PKT[(r,s) in set[:PKT]] >= lo, start = start_value(PK[yrl,r,s])/(1+value(ir)); # Terminal capital price

# consumer
# !!!! Updates needed for Intertemporal Consumption
@variable(cge,RA[yr in years, r in set[:r]]>=lo,start = value(qvm[yr])*value(pvm[yr])*value(c0[r])); #Representative Agent

###############################
# -- PLACEHOLDER VARIABLES --
###############################

#cobb-douglas function for value added (VA)
@NLexpression(cge,CVA[yr in years,r in set[:r],s in set[:s]],
  PL[yr,r]^alpha_kl[r,s] * ((haskey(RK.lookup[1], (yr,r,s)) ? RK[(yr,r,s)] : 1.0)/rk0)^ (1-alpha_kl[r,s]) );

#demand for labor in VA
@NLexpression(cge,AL[yr in years, r in set[:r], s in set[:s]], ld0[r,s] * CVA[yr,r,s] / PL[yr,r] );

#demand for capital in VA
@NLexpression(cge,AK[yr in years, r in set[:r],s in set[:s]],
  kd0[r,s] * CVA[yr,r,s] / ((haskey(RK.lookup[1], (yr,r,s)) ? RK[(yr,r,s)] : 1.0)/rk0));

#CES function for output demand - including
# exports (absent of re-exports) times the price for foreign exchange,
# region's supply to national market times the national market price
# regional supply to local market times domestic price

@NLexpression(cge,RX[yr in years,r in set[:r],g in set[:g]],
  (alpha_x[r,g]*PFX[yr]^(1 + et_x[r,g])+alpha_n[r,g]*PN[yr,g]^(1 + et_x[r,g])+alpha_d[r,g]*(haskey(PD.lookup[1], (yr,r,g)) ? PD[(yr,r,g)] : 1.0)^(1 + et_x[r,g]))^(1/(1 + et_x[r,g])) );

#demand for exports via demand function
@NLexpression(cge,AX[yr in years,r in set[:r],g in set[:g]], (x0[r,g] - rx0[r,g])*(PFX[yr]/RX[yr,r,g])^et_x[r,g] );

#demand for contribution to national market
@NLexpression(cge,AN[yr in years,r in set[:r],g in set[:g]], xn0[r,g]*(PN[yr,g]/(RX[yr,r,g]))^et_x[r,g] );

#demand for regionals supply to local market
@NLexpression(cge,AD[yr in years,r in set[:r],g in set[:g]],
  xd0[r,g] * ((haskey(PD.lookup[1], (yr,r,g)) ? PD[(yr,r,g)] : 1.0) / (RX[yr,r,g]))^et_x[r,g] );

# CES function for tradeoff between national and domestic market
@NLexpression(cge,CDN[yr in years,r in set[:r],g in set[:g]],
(theta_n[r,g]*PN[yr,g]^(1-es_d[r,g])+(1-theta_n[r,g])*(haskey(PD.lookup[1], (yr,r,g)) ? PD[(yr,r,g)] : 1.0)^(1-es_d[r,g]))^(1/(1-es_d[r,g])) );

# CES function for tradeoff between domestic consumption and foreign exports
# recall tm in the import tariff thus tm / tm0 is the relative change in import tariff rates
@NLexpression(cge,CDM[yr in years,r in set[:r],g in set[:g]],
((1-theta_m[r,g])*CDN[yr,r,g]^(1-es_f[r,g])+theta_m[r,g]*(PFX[yr]*(1+tm[r,g])/(1+tm0[r,g]))^(1-es_f[r,g]))^(1/(1-es_f[r,g])) );

# set[:r] demand from the national market <- note nesting of CDN in CDM
@NLexpression(cge,DN[yr in years,r in set[:r],g in set[:g]],
  nd0[r,g]*(CDN[yr,r,g]/PN[yr,g])^es_d[r,g]*(CDM[yr,r,g]/CDN[yr,r,g])^es_f[r,g] );

# region demand from local market <- note nesting of CDN in CDM
@NLexpression(cge,DD[yr in years,r in set[:r],g in set[:g]],
  dd0[r,g]*(CDN[yr,r,g]/(haskey(PD.lookup[1], (yr,r,g)) ? PD[(yr,r,g)] : 1.0))^es_d[r,g]*(CDM[yr,r,g]/CDN[yr,r,g])^es_f[r,g] );

# import demand
@NLexpression(cge,MD[yr in years,r in set[:r],g in set[:g]],
  m0[r,g]*(CDM[yr,r,g]*(1+tm[r,g])/(PFX[yr]*(1+tm0[r,g])))^es_f[r,g] );

# final demand
@NLexpression(cge,CD[yr in years,r in set[:r],g in set[:g]],
  cd0[r,g]*PC[yr,r] / (haskey(PA.lookup[1], (yr,r,g)) ? PA[(yr,r,g)] : 1.0) );

###############################
# -- Zero Profit Conditions --
###############################

@mapping(cge,profit_y[(yr,r,s) in set[:Y]],
# cost of intermediate demand
    sum((haskey(PA.lookup[1], (yr,r,g)) ? PA[(yr,r,g)] : 1.0) * id0[r,g,s] for g in set[:g])
# cost of labor inputs
    + PL[yr,r] * AL[yr,r,s]
# cost of capital inputs
    + ((haskey(RK.lookup[1], (yr,r,s)) ? RK[(yr,r,s)] : 1.0)/rk0)* AK[yr,r,s]
    -
# revenue from sectoral supply (take note of r/s/g indices on ys0)
    sum((haskey(PY.lookup[1], (yr,r,g)) ? PY[(yr,r,g)] : 1.0)  * ys0[r,s,g] for g in set[:g]) * (1-ty[r,s])
);

@mapping(cge,profit_x[(yr,r,g) in set[:X]],
# output 'cost' from aggregate supply
    (haskey(PY.lookup[1], (yr,r,g)) ? PY[(yr,r,g)] : 1.0) * s0[r,g]
    - (
# revenues from foreign exchange
        PFX[yr] * AX[yr,r,g]
# revenues from national market
        + PN[yr,g] * AN[yr,r,g]
# revenues from domestic market
        + (haskey(PD.lookup[1], (yr,r,g)) ? PD[(yr,r,g)] : 1.0) * AD[yr,r,g]
    )
);

@mapping(cge,profit_a[(yr,r,g) in set[:A]],
# costs from national market
    PN[yr,g] * DN[yr,r,g]
# costs from domestic market
    + (haskey(PD.lookup[1], (yr,r,g)) ? PD[(yr,r,g)] : 1.0) * DD[yr,r,g]
# costs from imports, including import tariff
    + PFX[yr] * (1+tm[r,g]) * MD[yr,r,g]
# costs of margin demand
    + sum(PM[yr,r,m] * md0[r,m,g] for m in set[:m])
    - (
# revenues from regional market based on armington supply
        (haskey(PA.lookup[1], (yr,r,g)) ? PA[(yr,r,g)] : 1.0) * (1-ta[r,g]) * a0[r,g]
# revenues from re-exports
        + PFX[yr] * rx0[r,g]
    )
);

@mapping(cge, profit_c[yr in years,r in set[:r]],
# costs of inputs - computed as final demand times regional market prices
    sum((haskey(PA.lookup[1], (yr,r,g)) ? PA[(yr,r,g)] : 1.0) * CD[yr,r,g] for g in set[:g])
    -
# revenues/benefit computed as CPI * reference consumption
    PC[yr,r] * c0[r]
);

@mapping(cge,profit_k[(yr,r,s) in set[:PK]],
    (haskey(PK.lookup[1], (yr,r,s)) ? PK[(yr,r,s)] : 1.0)
    - (
        (haskey(RK.lookup[1], (yr,r,s)) ? RK[(yr,r,s)] : 1.0)
        + (1-dr) * (yr!=yrl ? (haskey(PK.lookup[1], (yr+1,r,s)) ? PK[(yr+1,r,s)] : 1.0) : (haskey(PKT.lookup[1], (r,s)) ? PKT[(r,s)] : 1.0))
    )
);

@mapping(cge,profit_i[(yr,r,s) in set[:PK]],
    (haskey(PY.lookup[1], (yr,r,s)) ? PY[(yr,r,s)] : 1.0)
    -
    (yr!=yrl ? (haskey(PK.lookup[1], (yr+1,r,s)) ? PK[(yr+1,r,s)] : 1.0) : (haskey(PKT.lookup[1], (r,s)) ? PKT[(r,s)] : 1.0))
);

@mapping(cge,profit_ms[yr in years,r in set[:r],m in set[:m]],
# provision of set[:m] to national market
    sum(PN[yr,gm]   * nm0[r,gm,m] for gm in set[:gm])
# provision of set[:m] to domestic market
    + sum((haskey(PD.lookup[1], (yr,r,gm)) ? PD[(yr,r,gm)] : 1.0) * dm0[r,gm,m] for gm in set[:gm])
    -
# total margin demand
    PM[yr,r,m] * sum(md0[r,m,gm] for gm in set[:gm])
);

###################################
# -- Market Clearing Conditions --
###################################

@mapping(cge,market_pa[(yr,r,g) in set[:PA]],
# absorption or supply
    (haskey(A.lookup[1], (yr,r,g)) ? A[(yr,r,g)] : 1.0) * a0[r,g]
    - (
# government demand (exogenous)
        g0[r,g]*qvm[yr]
# demand for investment (exogenous)
        + i0[r,g]*qvm[yr]
# final demand
        + C[yr,r] * CD[yr,r,g]
# intermediate demand
        + sum((haskey(Y.lookup[1], (yr,r,s)) ? Y[(yr,r,s)] : 1.0) * id0[r,g,s] for s in set[:s] if (yr,r,s) in set[:Y])
    )
)

@mapping(cge,market_py[(yr,r,g) in set[:PY]],
# sectoral supply
    sum((haskey(Y.lookup[1], (yr,r,s)) ? Y[(yr,r,s)] : 1.0) *ys0[r,s,g] for s in set[:s])
# household production (exogenous)
    + yh0[r,g]*qvm[yr]
    -
# aggregate supply (akin to market demand)
    (haskey(X.lookup[1], (yr,r,g)) ? X[(yr,r,g)] : 1.0) * s0[r,g]
);

@mapping(cge,market_pd[(yr,r,g) in set[:PD]],
# aggregate supply
    (haskey(X.lookup[1], (yr,r,g)) ? X[(yr,r,g)] : 1.0)  * AD[yr,r,g]
    - (
# demand for local market
        (haskey(A.lookup[1], (yr,r,g)) ? A[(yr,r,g)] : 1.0) * DD[yr,r,g]
# margin supply from local market
        + sum(MS[yr,r,m] * dm0[r,g,m] for m in set[:m] if (g in set[:gm] ) )
    )
);

@mapping(cge,market_pn[yr in years,g in set[:g]],
# supply to the national market
    sum((haskey(X.lookup[1], (yr,r,g)) ? X[(yr,r,g)] : 1.0)  * AN[yr,r,g] for r in set[:r])
    - (
# demand from the national market
        sum((haskey(A.lookup[1], (yr,r,g)) ? A[(yr,r,g)] : 1.0) * DN[yr,r,g] for r in set[:r])
# market supply to the national market
        + sum(MS[yr,r,m] * nm0[r,g,m] for r in set[:r] for m in set[:m] if (g in set[:gm]) )
    )
);

@mapping(cge,market_pl[yr in years,r in set[:r]],
# supply of labor
    sum(ld0[r,s]*qvm[yr] for s in set[:s])
    -
# demand for labor in all set[:s]
    sum((haskey(Y.lookup[1], (yr,r,s)) ? Y[(yr,r,s)] : 1.0) * AL[yr,r,s] for s in set[:s])
);

@mapping(cge,market_pk[(yr,r,s) in set[:PK]],
# if first year, initial capital
# else investment plus previous year's decayed capital
    (yr==yrf ? kd0[r,s] : (haskey(I.lookup[1], (yr-1,r,s)) ? I[(yr-1,r,s)] : 0.0))
    +(1-dr) * (yr>yrf ? (haskey(K.lookup[1], (yr-1,r,s)) ? K[(yr-1,r,s)] : 0.0) : 0.0)
    -
# capital in year [yr]
    (haskey(K.lookup[1], (yr,r,s)) ? K[(yr,r,s)] : 0.0)
);

@mapping(cge,market_rk[(yr,r,s) in set[:PK]],
    (haskey(K.lookup[1], (yr,r,s)) ? K[(yr,r,s)] : 0.0)
    -
    (haskey(Y.lookup[1], (yr,r,s)) ? Y[(yr,r,s)] : 1.0) * AK[yr,r,s]
);

@mapping(cge,market_pkt[(r,s) in set[:PKT]],
    (1-dr)*(haskey(K.lookup[1], (yrl,r,s)) ? K[(yrl,r,s)] : 0.0)
    + (haskey(I.lookup[1], (yrl,r,s)) ? I[(yrl,r,s)] : 0.0)
    -
    (haskey(TK.lookup[1], (r,s)) ? TK[(r,s)] : 0.0)
);

@mapping(cge,termk[(r,s) in set[:PKT]],
    (haskey(I.lookup[1], (yrl,r,s)) ? I[(yrl,r,s)] : 0.0)
    / ((haskey(I.lookup[1], (yrl-1,r,s)) ? I[(yrl-1,r,s)] : 0.0))
    -
    (haskey(Y.lookup[1], (yrl,r,s)) ? Y[(yrl,r,s)] : 1.0)
    / (haskey(Y.lookup[1], (yrl-1,r,s)) ? Y[(yrl-1,r,s)] : 1.0)
);

@mapping(cge,market_pm[yr in years,r in set[:r], m in set[:m]],
# margin supply
    MS[yr,r,m] * sum(md0[r,m,gm] for gm in set[:gm])
    -
# margin demand
    sum((haskey(A.lookup[1], (yr,r,g)) ? A[(yr,r,g)] : 1.0) * md0[r,m,g] for g in set[:g])
);

@mapping(cge,market_pfx[yr in years],
# balance of payments (exogenous)
    sum(bopdef0[r] for r in set[:r])*qvm[yr]
# supply of exports
    + sum((haskey(X.lookup[1], (yr,r,g)) ? X[(yr,r,g)] : 1.0)  * AX[yr,r,g] for r in set[:r] for g in set[:g])
# supply of re-exports
    + sum((haskey(A.lookup[1], (yr,r,g)) ? A[(yr,r,g)] : 1.0) * rx0[r,g] for r in set[:r] for g in set[:g] if (yr,r,g) in set[:A])
    -
# import demand
    sum((haskey(A.lookup[1], (yr,r,g)) ? A[(yr,r,g)] : 1.0) * MD[yr,r,g] for r in set[:r] for g in set[:g] if (yr,r,g) in set[:A])
);

@mapping(cge,market_pc[yr in years,r in set[:r]],
# a period's final demand
    C[yr,r] * c0[r]
    -
# consumption / utiltiy
    RA[yr,r] / PC[yr,r]
);

# !!!! Updates needed for Intertemporal Consumption
@mapping(cge,income_ra[yr in years,r in set[:r]],
# consumption/utility
    RA[yr,r]
    -
    (
# labor income
        PL[yr,r] * sum(ld0[r,s] for s in set[:s])*qvm[yr]
# provision of household supply
        + sum( (haskey(PY.lookup[1], (yr,r,g)) ? PY[(yr,r,g)] : 1.0) * yh0[r,g]*qvm[yr] for g in set[:g])
# revenue or costs of foreign exchange including household adjustment
        + PFX[yr] * (bopdef0[r] + hhadj[r])*qvm[yr]
# government and investment provision
        - sum((haskey(PA.lookup[1], (yr,r,g)) ? PA[(yr,r,g)] : 1.0) * (g0[r,g] + i0[r,g])*qvm[yr] for g in set[:g])
# import taxes - assumes lumpsum recycling
        + sum((haskey(A.lookup[1], (yr,r,g)) ? A[(yr,r,g)] : 1.0) * MD[yr,r,g] * PFX[yr] * tm[r,g] for g in set[:g] if (yr,r,g) in set[:A])
# taxes on intermediate demand - assumes lumpsum recycling
        + sum((haskey(A.lookup[1], (yr,r,g)) ? A[(yr,r,g)] : 1.0) * a0[r,g]*(haskey(PA.lookup[1], (yr,r,g)) ? PA[(yr,r,g)] : 1.0)*ta[r,g] for g in set[:g] if (yr,r,g) in set[:A])
# production taxes - assumes lumpsum recycling
        + sum( pvm[yr]*(haskey(Y.lookup[1], (yr,r,s)) ? Y[(yr,r,s)] : 1.0) * ys0[r,s,g] * ty[r,s] for s in set[:s], g in set[:g])
# capital income
        + (1-islast[yr]) * sum((haskey(PK.lookup[1], (yr,r,s)) ? PK[(yr,r,s)] : 1.0) * (haskey(K.lookup[1], (yr,r,s)) ? K[(yr,r,s)] : 0.0) for s in set[:s]) / (1+ir)
        + (islast[yr]) * sum((haskey(PKT.lookup[1], (r,s)) ? PKT[(r,s)] : 1.0) * (haskey(TK.lookup[1], (r,s)) ? TK[(r,s)] : 0.0) for s in set[:s]) / (1+gr)
    )
);

####################################
# -- Complementarity Conditions --
####################################

# define complementarity conditions
# note the pattern of ZPC -> primal variable  &  MCC -> dual variable (price)
@complementarity(cge,profit_y,Y);
@complementarity(cge,profit_x,X);
@complementarity(cge,profit_a,A);
@complementarity(cge,profit_c,C);
@complementarity(cge,profit_ms,MS);
@complementarity(cge,market_pa,PA);
@complementarity(cge,market_py,PY);
@complementarity(cge,market_pd,PD);
@complementarity(cge,market_pn,PN);
@complementarity(cge,market_pl,PL);
@complementarity(cge,market_pm,PM);
@complementarity(cge,market_pc,PC);
@complementarity(cge,market_pfx,PFX);
@complementarity(cge,income_ra,RA);
@complementarity(cge,market_pk,PK);
@complementarity(cge,market_rk,RK);
@complementarity(cge,market_pkt,PKT)
@complementarity(cge,termk,TK)
@complementarity(cge,profit_k,K)
@complementarity(cge,profit_i,I)

####################
# -- Model Solve --
####################

#set up the options for the path solver
#PATHSolver.options(convergence_tolerance=1e-6, output=:yes, time_limit=3600)
PATHSolver.options(convergence_tolerance=1e-6, output=:yes, time_limit=3600, cumulative_iteration_limit=0)

# solve the model
status = solveMCP(cge)

max_stars_repo_name: NREL/SLiDE

max_issues_repo_name: NREL/SLiDE

max_forks_repo_name: NREL/SLiDE


text:
using Base.Test, DataFrames, EventHistory

function test_show(x)
    io = IOBuffer()
    show(io, x)
end

d = DataFrame(start=[1,2,5,2,1,7,3,4,8,8],
              stop=[2,3,6,7,8,9,9,9,14,17],
              event=[1,1,1,1,1,1,1,0,0,0],
              x=[1,0,0,1,0,1,1,1,0,0]);
d[:S] = EventHistory.Event([:start,:stop,:event],d);
mm = EventHistory.phreg(S~x,d)

@test_approx_eq_eps EventHistory.coef(mm) [-0.02110521] 1e-8
@test_approx_eq_eps mm.invhess 0.6323063 1e-6

@test mm.grad[1]<1e-6

p = EventHistory.predict(mm,X=[0; 1], time=[6]; order=true);
@test_approx_eq_eps p[2:3] [0.3522254,0.3599854] 1e-6

ext: jl


text:
mutable struct SystemState
    i::Int64
    θ::Vector{Float64}
    σ2::Matrix{Float64}
    γs::Vector{Float64}
    γmins::Vector{Float64}
    ns::Vector{Int64}
    θs::Vector{Float64}

    SystemState(; i::Int64 = 0, θ::Vector{Float64} = Vector{Float64}(), σ2::Matrix{Float64} = Matrix{Float64}(undef, 0, 0), γs::Vector{Float64} = Vector{Float64}(), δmin_history::Vector{Float64} = Vector{Float64}(), ns::Vector{Float64} = Vector{Float64}(), θs::Vector{Float64} = Vector{Float64}()) = new(i, θ, σ2, γs, δmin_history, ns, θs)
end

function Base.show(io::IO, s::SystemState)
    if s.i == 0
        print("\n\t\tIterations: ", s.i, ",\n\t\tCurrent θ: ", s.θ, ",\n\t\tVariance ", diag(s.σ2), ",\n\t\tCost Value: Inf,\n\t\tMin Cost Value: Inf.")
    else
        print("\n\t\tIterations: ", s.i, ",\n\t\tCurrent θ: ", s.θ, ",\n\t\tVariance ", diag(s.σ2), ",\n\t\tCost Value: ", s.γs[end], ",\n\t\tMin Cost Value: ", s.γmins[end], ",\n\t\tIteration samples: ", s.ns[end], ",\n\t\tCumulative samples: ", sum(s.ns), ".")
    end
end

max_stars_repo_name: pzuliani/SPICE

max_issues_repo_name: pzuliani/SPICE

max_forks_repo_name: pzuliani/SPICE


text:
# test packing functions

using RoME
using Distributions
using TransformUtils
using KernelDensityEstimate

using Base.Test

N = 100
fg = initfg()

initCov = diagm([0.03;0.03;0.001])
odoCov = diagm([3.0;3.0;0.01])

# Some starting position
v1 = addNode!(fg, :x1, zeros(3,1), diagm([1.0;1.0;0.1]), N=N)
ipp = PriorPose2(zeros(3,1), initCov, [1.0])
f1  = addFactor!(fg,[v1], ipp)

# and a second pose
v2 = addNode!(fg, :x2, vectoarr2([50.0;0.0;pi/2]), diagm([1.0;1.0;0.05]), N=N)
ppc = Pose2Pose2(vectoarr2([50.0;0.0;pi/2]), odoCov, [1.0])
f2 = addFactor!(fg, [v1;v2], ppc)

println("test conversions of PriorPose2")
dd = convert(PackedPriorPose2, ipp)
upd = convert(RoME.PriorPose2, dd)

@test RoME.compare(ipp, upd)

packeddata = FNDencode(IncrementalInference.PackedFunctionNodeData{RoME.PackedPriorPose2}, getData(f1))
unpackeddata = FNDdecode(IncrementalInference.FunctionNodeData{GenericWrapParam{RoME.PriorPose2}}, packeddata)

@test IncrementalInference.compare(getData(f1), unpackeddata)

packedv1data = VNDencoder(IncrementalInference.PackedVariableNodeData, getData(v1))
upv1data = VNDdecoder(IncrementalInference.VariableNodeData, packedv1data)

@test IncrementalInference.compare(getData(v1), upv1data)

println("test conversions of Pose2Pose2")

dd = convert(PackedPose2Pose2, ppc)
upd = convert(RoME.Pose2Pose2, dd)

# TODO -- fix ambiguity in compare function
@test RoME.compare(ppc, upd)

packeddata = FNDencode(IncrementalInference.PackedFunctionNodeData{RoME.PackedPose2Pose2}, getData(f2))
unpackeddata = FNDdecode(IncrementalInference.FunctionNodeData{GenericWrapParam{RoME.Pose2Pose2}}, packeddata)

# TODO -- fix ambibuity in compare function
@test IncrementalInference.compare(getData(f2), unpackeddata)

println("test conversions of Pose2DPoint2DBearingRange")

# and a second pose
v3 = addNode!(fg, :l1, vectoarr2([50.0,50.0]), diagm([1.0;1.0]), N=N)
# ppc = Pose2DPoint2DBearingRange([50.0;0.0;pi/2], 0.01*eye(2), [1.0])
ppbr = Pose2DPoint2DBearingRange{Normal{Float64}, Normal{Float64}}(
              Normal(0.0, 0.005 ),
              Normal(50, 0.5) )
f3 = addFactor!(fg, [:x2;:l1], ppbr)

dd = convert(PackedPose2DPoint2DBearingRange, ppbr)
upd = convert(
        RoME.Pose2DPoint2DBearingRange{Normal{Float64},Normal{Float64}},
        dd
      )

packeddata = FNDencode(IncrementalInference.PackedFunctionNodeData{RoME.PackedPose2DPoint2DBearingRange}, getData(f3))
unpackeddata = FNDdecode(IncrementalInference.FunctionNodeData{GenericWrapParam{RoME.Pose2DPoint2DBearingRange{Normal{Float64},Normal{Float64}}}}, packeddata)

@test ppbr.bearing.μ == unpackeddata.fnc.usrfnc!.bearing.μ
@test ppbr.bearing.σ == unpackeddata.fnc.usrfnc!.bearing.σ

@test ppbr.range.μ == unpackeddata.fnc.usrfnc!.range.μ
@test ppbr.range.σ == unpackeddata.fnc.usrfnc!.range.σ

println("test conversions of PriorPose3")

# parameters
N = 300
initCov = eye(6)
[initCov[i,i] = 0.01 for i in 4:6];
odoCov = deepcopy(initCov)

# start new factor graph
fg = initfg()

v1 = addNode!(fg,:x1,  0.1*randn(6,N))
ipp = PriorPose3(MvNormal(zeros(6), initCov) )
f1  = addFactor!(fg,[v1], ipp)

dd = convert(PackedPriorPose3, ipp)
upd = convert(RoME.PriorPose3, dd)

# @test TransformUtils.compare(ipp.Zi, upd.Zi)
@test norm(ipp.Zi.μ - upd.Zi.μ) < 1e-10
@test norm(ipp.Zi.Σ.mat - upd.Zi.Σ.mat) < 1e-8

packeddata = FNDencode(IncrementalInference.PackedFunctionNodeData{RoME.PackedPriorPose3}, getData(f1))
unpackeddata = FNDdecode(IncrementalInference.FunctionNodeData{GenericWrapParam{RoME.PriorPose3}}, packeddata)

# TODO -- fix ambibuity in compare function
@test IncrementalInference.compare(getData(f1), unpackeddata)

packedv1data = VNDencoder(IncrementalInference.PackedVariableNodeData, getData(v1))
upv1data = VNDdecoder(IncrementalInference.VariableNodeData, packedv1data)

@test IncrementalInference.compare(getData(v1), upv1data)

println("test conversions of Pose3Pose3")

pp3 = Pose3Pose3( MvNormal([25.0;0;0;0;0;0], odoCov) )
v2, f2 = addOdoFG!(fg, pp3 )

dd = convert(PackedPose3Pose3, pp3)
upd = convert(RoME.Pose3Pose3, dd)

@test norm(pp3.Zij.μ - upd.Zij.μ) < 1e-10
@test norm(pp3.Zij.Σ.mat - upd.Zij.Σ.mat) < 1e-8

packeddata = FNDencode(IncrementalInference.PackedFunctionNodeData{RoME.PackedPose3Pose3}, getData(f2))
unpackeddata = FNDdecode(IncrementalInference.FunctionNodeData{GenericWrapParam{RoME.Pose3Pose3}}, packeddata)

# TODO -- fix ambibuity in compare function
@test IncrementalInference.compare(getData(f2), unpackeddata)

println("test conversions of Pose3Pose3NH")

odo = SE3(randn(3), convert(Quaternion, so3(0.1*randn(3)) ))
odoc = Pose3Pose3NH( MvNormal(veeEuler(odo),odoCov), [0.5;0.5])

f3 = addFactor!(fg,[:x1;:x2],odoc)

dd = convert(PackedPose3Pose3NH, odoc)
upd = convert(RoME.Pose3Pose3NH, dd)

@test norm(odoc.Zij.μ - upd.Zij.μ) < 1e-8
@test norm(odoc.Zij.Σ.mat - upd.Zij.Σ.mat) < 1e-8
@test norm(odoc.nullhypothesis.p - upd.nullhypothesis.p) < 1e-8

packeddata = FNDencode(IncrementalInference.PackedFunctionNodeData{RoME.PackedPose3Pose3NH}, getData(f3))
unpackeddata = FNDdecode(IncrementalInference.FunctionNodeData{GenericWrapParam{RoME.Pose3Pose3NH}}, packeddata)

# TODO -- fix ambibuity in compare function
@test IncrementalInference.compare(getData(f3), unpackeddata)

println("test conversions of PartialPriorRollPitchZ")

prpz = PartialPriorRollPitchZ(MvNormal([0.0;0.5],0.1*eye(2)),Normal(3.0,0.5))

pprpz = convert(PackedPartialPriorRollPitchZ, prpz)
unp = convert(PartialPriorRollPitchZ, pprpz)

@test RoME.compare(prpz, unp)

println("test conversions of PartialPose3XYYaw")

xyy = PartialPose3XYYaw(MvNormal([1.0;2.0;0.5],0.1*eye(3)))

pxyy = convert(PackedPartialPose3XYYaw, xyy)
unp = convert(PartialPose3XYYaw, pxyy)

@test RoME.compare(xyy, unp)

println("test conversions of PartialPose3XYYawNH")

xyy = PartialPose3XYYawNH(MvNormal([1.0;2.0;0.5],0.1*eye(3)), [0.6;0.4])

pxyy = convert(PackedPartialPose3XYYawNH, xyy)
unp = convert(PartialPose3XYYawNH, pxyy)

@test RoME.compare(xyy, unp)

println("test PriorPoint2DensityNH")

prpt2 = PriorPoint2DensityNH(kde!(randn(2,100)),[0.25;0.75]  )

pprpt2 = convert(PackedPriorPoint2DensityNH, prpt2)
uprpt2 = convert(PriorPoint2DensityNH, pprpt2)

@test norm(getPoints(prpt2.belief)-getPoints(uprpt2.belief)) < 1e-8

@test norm(prpt2.nullhypothesis.p-uprpt2.nullhypothesis.p) < 1e-8

#


text:
module PLists

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

include("plist_lexer.jl")
include("xml_lexer.jl")
include("parser.jl")
include("plist_parser.jl")
include("plist_writer.jl")
include("xml_dom.jl")
include("xml_parser.jl")
include("xml_show.jl")
include("xml_plist_parser.jl")

end


text:
module TestRDA
    using Base.Test
    using DataFrames
    using Compat

    # R code generating test .rdas
        # df = data.frame(num=c(1.1, 2.2))
        # save(df, file='minimal.rda')
        # save(df, file='minimal_ascii.rda', ascii=TRUE)

        # df['int'] = c(1L, 2L)
        # df['logi'] = c(TRUE, FALSE)
        # df['chr'] = c('ab', 'c')
        # df['factor'] = factor(df$chr)
        # df['cplx'] = complex( real=c(1.1,0.0), imaginary=c(0.5,1.0) )
        # #utf=c('Ж', '∰')) R handles it, read_rda doesn't.
        # save(df, file='types.rda')
        # save(df, file='types_ascii.rda', ascii=TRUE)

        # df[2, ] = NA
        # df[3, ] = df[2, ]
        # df[3,'num'] = NaN
        # df[,'cplx'] = complex( real=c(1.1,1,NaN), imaginary=c(NA,NaN,0) )
        # save(df, file='NAs.rda')
        # save(df, file='NAs_ascii.rda', ascii=TRUE)

        # names(df) = c('end', '!', '1', '%_B*\tC*', NA, 'x')
        # save(df, file='names.rda')
        # save(df, file='names_ascii.rda', ascii=TRUE)

    testdir = dirname(@__FILE__)

    df = DataFrame(num = [1.1, 2.2])
    @test isequal(DataFrame(read_rda("$testdir/data/RDA/minimal.rda")["df"]), df)
    @test isequal(DataFrame(open(read_rda,"$testdir/data/RDA/minimal_ascii.rda")["df"]), df)
    @test isequal(read_rda("$testdir/data/RDA/minimal.rda",convertdataframes=true)["df"], df)

    df[:int] = Int32[1, 2]
    df[:logi] = [true, false]
    df[:chr] = ["ab", "c"]
    df[:factor] = pool(df[:chr])
    df[:cplx] = Complex128[1.1+0.5im, 1.0im]
    @test isequal(DataFrame(read_rda("$testdir/data/RDA/types.rda")["df"]), df)
    @test isequal(DataFrame(read_rda("$testdir/data/RDA/types_ascii.rda")["df"]), df)

    df[2, :] = NA
    append!(df, df[2, :])
    df[3, :num] = NaN
    df[:, :cplx] = @data [NA, @compat(Complex128(1,NaN)), NaN]
    @test isequal(DataFrame(read_rda("$testdir/data/RDA/NAs.rda")["df"]), df)
    # ASCII format saves NaN as NA
    df[3, :num] = NA
    df[:, :cplx] = @data [NA, NA, NA]
    @test isequal(DataFrame(read_rda("$testdir/data/RDA/NAs_ascii.rda")["df"]), df)

    rda_names = names(DataFrame(read_rda("$testdir/data/RDA/names.rda")["df"]))
    expected_names = [:_end, :x!, :x1, :_B_C_, :x, :x_1]
    @test rda_names == expected_names
    rda_names = names(DataFrame(read_rda("$testdir/data/RDA/names_ascii.rda")["df"]))
    @test rda_names == [:_end, :x!, :x1, :_B_C_, :x, :x_1]

end


text:
@testset "WeighedMean fit!" begin
    test_fit(WeightedMean{Float64}, x, w, mean, (x, w) -> mean(x, weights(w)))
    test_fit(WeightedMean{Float32}, x, w, mean, (x, w) -> mean(x, weights(w)))

    test_fit(WeightedMean{Float64}, xmis, wmis, mean,
             (x, w) -> sum(skipmissing(x .* w)) / sum(skipmissing(w)))

    m = mean(x, weights(w))
    s = sum(broadcast(*, w, x))

    mval = mean(fit!(WeightedMean(), x, w))
    mzip = mean(fit!(WeightedMean(), zip(x, w)))
    sval = sum(fit!(WeightedMean(), x, w))
    szip = sum(fit!(WeightedMean(), zip(x, w)))

    @test m ≈ mzip
    @test m ≈ mval
    @test s ≈ szip
    @test s ≈ sval

    m2val = mean(fit!(WeightedMean(Float32), x, w))
    m2zip = mean(fit!(WeightedMean(Float32), zip(x, w)))

    @test m ≈ m2zip
    @test m ≈ m2val
    @test typeof(m2val) == Float32
    @test typeof(m2zip) == Float32
end

@testset "WeighedMean merge!" begin
    m = mean(x, weights(w))
    wm = fit!(WeightedMean(), x, w) |> mean

    om = map(x, w) do xi, wi
        fit!(WeightedMean(), xi, wi)
    end;

    rm = reduce(merge!, om) |> mean
    rm2 = merge!(
        fit!(WeightedMean(), x[1:end ÷ 2], w[1:end ÷ 2]),
        fit!(WeightedMean(), x[end ÷ 2 + 1:end], w[end ÷ 2 + 1:end])
    ) |> mean

    rm_32 = reduce(merge!, om, init = WeightedMean(Float32)) |> mean
    rm2_32 = merge!(
        fit!(WeightedMean(Float32), x[1:end ÷ 2], w[1:end ÷ 2]),
        fit!(WeightedMean(), x[end ÷ 2 + 1:end], w[end ÷ 2 + 1:end])
    ) |> mean

    @test rm ≈ m
    @test rm ≈ wm
    @test rm2 ≈ m
    @test rm_32 ≈ m
    @test rm2_32 ≈ m

    @test typeof(rm) == Float64
    @test typeof(rm2) == Float64
    @test typeof(rm_32) == Float32
    @test typeof(rm2_32) == Float32
end

@testset "WeightedMean constructor" begin
    @test WeightedMean{Float64}() == WeightedMean()
    @test WeightedMean{Float32}() == WeightedMean(Float32)
    @test WeightedMean() == WeightedMean(0.0, 0.0, 0)
    @test WeightedMean() == WeightedMean{Float64}(0, 0, 0)
end


text:
using ..Fields: Field, PrimeField, BinaryField

abstract type EllipticCurve end
abstract type AbstractPoint end

struct AffinePoint{E <: EllipticCurve, T <: Field} <: AbstractPoint
    x::T
    y::T

    AffinePoint{E, F}(x::F, y::F) where {E <: EllipticCurve, F <: Field} = new{E, F}(x, y)
    AffinePoint{E}(x::F, y::F) where {E <: EllipticCurve, F <: Field} = new{E, F}(x, y)
end

AffinePoint{E, F}(x, y) where {E <: EllipticCurve, F <: Field} = AffinePoint{E, F}(F <| x, F <| y)

eq(::Type{AffinePoint{EQ, F}}) where {EQ <: EllipticCurve, F <: Field} = EQ
field(::Type{AffinePoint{EQ, F}}) where {EQ <: EllipticCurve, F <: Field} = F

### More through subtyping could be in here
Base.convert(::Type{P}, x::Tuple{BigInt, BigInt}) where P <: AffinePoint = P(x...)
Base.convert(::Type{P}, x::Tuple{BitVector, BitVector}) where P <: AffinePoint = P(x...)

gx(p::AffinePoint) = p.x
gy(p::AffinePoint) = p.y

modulus(::Type{AffinePoint{<:EllipticCurve, F}}) where F <: PrimeField = modulus(F)

Base.zero(::Type{AffinePoint{E, F}}) where {E <: EllipticCurve, F <: Field} = AffinePoint{E , F}(zero(F), zero(F))
Base.zero(x::P) where P <: AffinePoint = zero(P)

Base.:(==)(x::P, y::P) where P <: AffinePoint = x.x == y.x && x.y == y.y

Base.:-(u::P) where P <: AffinePoint{<:EllipticCurve, <: PrimeField} = P(u.x, -u.y)
Base.:-(u::P) where P <: AffinePoint{<:EllipticCurve, <: BinaryField} = P(u.x, u.x + u.y)

Base.:-(u::P, v::P) where P <: AffinePoint = u + (-v)

Base.isless(x::P, y::P) where P <: AbstractPoint = gx(x) == gx(y) ? gx(x) < gx(y) : gy(x) < gy(y)

function Base.:*(P::AffinePoint, k::Integer)

    h = tobits(3*k)
    r = length(h) - 1
    e = tobits(k, length(h))

    R = P

    for i in 2:r #+1

        R = double(R) #R + R

        if h[i] == 1 && e[i] == 0
            R = R + P
        elseif h[i] == 0 && e[i] == 1
            R = R - P
        end
    end

    return R
end

Base.:*(k::Integer, P::AffinePoint) = P * k

function Base.isvalid(::Type{P}) where P <: AffinePoint

    g = generator(x)

    if oncurve(g) == false
        return false
    end

    if !(g * order(P) == zero(P))
        return false
    end

    return true
end

Base.isvalid(x::AffinePoint) = oncurve(x)

### Definition of some elliptic curves and coresponding operations

struct Weierstrass{a, b} <: EllipticCurve end # May assume that a, b are

a(::Type{Weierstrass{A, B}}) where {A, B} = A
b(::Type{Weierstrass{A, B}}) where {A, B} = B

a(::Type{AffinePoint{W, F}}) where {W <: Weierstrass, F <: Field} = F <| a(W)
b(::Type{AffinePoint{W, F}}) where {W <: Weierstrass, F <: Field} = F <| b(W)

function Base.:+(u::AffinePoint{E, F}, v::AffinePoint{E, F}) where {E <: Weierstrass, F <: Field}

    if u.x == v.x && u.y + v.y == zero(F) # I could put this case into ECPoint as assertion
        return zero(AffinePoint{E, F})
    end

    @assert u.x != v.x

    λ = (v.y - u.y)/(v.x - u.x)

    x₃ = λ^2 - u.x - v.x
    y₃ = λ * (u.x - x₃) - u.y

    return AffinePoint{E, F}(x₃, y₃)
end

function double(u::P) where P <: AffinePoint{<:Weierstrass, <:PrimeField}

    (; x, y) = u

    λ = (3 * x^2 + a(P))/ (2 * y)

    x₃ = λ^2 - 2x
    y₃ = λ*(x - x₃) - y

    return P(x₃, y₃)
end

function oncurve(u::P) where P <: AffinePoint{<:Weierstrass, <:Field}
    (; x, y) = u
    return y^2 == x^3 + a(P)*x + b(P)
end

struct BinaryCurve{a, b} <: EllipticCurve end

a(::Type{BinaryCurve{A, B}}) where {A, B} = A
b(::Type{BinaryCurve{A, B}}) where {A, B} = B

a(::Type{AffinePoint{W, F}}) where {W <: BinaryCurve, F <: Field} = F <| a(W)
b(::Type{AffinePoint{W, F}}) where {W <: BinaryCurve, F <: Field} = F <| b(W)

function Base.:+(u::AffinePoint{E, F}, v::AffinePoint{E, F}) where {E<:BinaryCurve, F<:BinaryField}

    if u.x == v.x && u.y + v.y + u.x == zero(F)
        return zero(AffinePoint{E, F})
    end

    @assert u.x != v.x

    _a = a(AffinePoint{E, F})

    λ = (v.y + u.y)/(v.x + u.x)

    x₃ = λ^2 + λ + u.x + v.x + _a
    y₃ = λ * (u.x + x₃) + x₃ + u.y

    return AffinePoint{E, F}(x₃, y₃)
end

function double(u::P) where P <: AffinePoint{<:BinaryCurve, <:BinaryField}

    (; x, y) = u

    λ = x + y/x

    x₃ = λ^2 + λ + a(P)
    y₃ = x^2 + (λ + one(λ))*x₃

    return P(x₃, y₃)
end

function oncurve(u::P) where P <: AffinePoint{<:BinaryCurve, <:Field}
    (; x, y) = u

    return y^2 + x*y == x^3 + a(P)*x^2 + b(P)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 54f89e42-7eb6-11eb-0a18-79732e59d627
begin
	using PlutoUI
	using NativeSVG
end

# ╔═╡ 21cefe9c-7eb3-11eb-2469-2b7b3d1eccc3
md"""# Arrays

This chapter presents one of Julia’s most useful built-in types, the array. You will also learn about objects and what can happen when you have more than one name for the same object."""

# ╔═╡ 62637476-7eb6-11eb-0cb4-01f126c2feb8
md"""## An Array Is a Sequence

Like a string, an *array* is a sequence of values. In a string, the values are characters; in an array, they can be any type. The values in an array are called *elements* or sometimes *items*.

There are several ways to create a new array; the simplest is to enclose the elements in square brackets (`[ ]`):

```julia
[10, 20, 30, 40]
["crunchy frog", "ram bladder", "lark vomit"]
```

The first example is an array of four integers. The second is an array of three strings. The elements of an array don’t have to be the same type. The following array contains a string, a float, an integer, and another array:

```julia
["spam", 2.0, 5, [10, 20]]
```

An array within another array is *nested*.

An array that contains no elements is called an empty array; you can create one with empty brackets, `[]`.

As you might expect, you can assign array values to variables:

```julia
julia> cheeses = ["Cheddar", "Edam", "Gouda"];

julia> numbers = [42, 123];

julia> empty = [];

julia> print(cheeses, " ", numbers, " ", empty)
["Cheddar", "Edam", "Gouda"] [42, 123] Any[]
```

The function `typeof` can be used to find out the type of the array:

```julia
julia> typeof(cheeses)
Array{String,1}
julia> typeof(numbers)
Array{Int64,1}
julia> typeof(empty)
Array{Any,1}
```

The number indicates the dimensions (we’ll talk more about this in “Arrays”). The array `empty` contains values of type `Any`; that is, it can hold values of all types.
"""

# ╔═╡ daca2c98-7eb6-11eb-1c09-8f1e08b31315
md"""## Arrays Are Mutable

The syntax for accessing the elements of an array is the same as for accessing the characters of a string—using the bracket operator. The expression inside the brackets specifies the index. Remember that the indices start at 1:

```julia
julia> cheeses[1]
"Cheddar"
```

Unlike strings, arrays are *mutable*; that is, their values can be changed. When the bracket operator appears on the left side of an assignment, it identifies the element of the array that will be assigned:

```julia
julia> numbers[2] = 5
5
julia> print(numbers)
[42, 5]
```

The second element of `numbers`, which used to be `123`, is now `5`.
Figure 10-1 shows the state diagrams for `cheeses`, `numbers`, and `empty`.
"""

# ╔═╡ 1b70a1d2-7eb7-11eb-2ab9-6739d60fcf89
Drawing(width=720, height=230) do
	defs() do
        marker(id="arrow", markerWidth="10", markerHeight="10", refX="0", refY="3", orient="auto", markerUnits="strokeWidth") do
      		path(d="M0,0 L0,6 L9,3 z", fill="black")
		end
	end
	text(x=300, y=55, font_family="JuliaMono, monospace", text_anchor="end", font_size="0.85rem", font_weight=600) do
		str("cheeses ->")
	end
    rect(x=310, y=10, width=160, height=85, fill="rgb(242, 242, 242)", stroke="black")
	text(x=320, y=30, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("1")
	end
	line(x1=335, y1=25, x2=355, y2=25, stroke="black", marker_end="url(#arrow)")
	text(x=370, y=30, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("\"Cheddar\"")
	end
	text(x=320, y=55, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("2")
	end
	line(x1=335, y1=50, x2=355, y2=50, stroke="black", marker_end="url(#arrow)")
	text(x=370, y=55, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("\"Edam\"")
	end
	text(x=320, y=80, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("3")
	end
	line(x1=335, y1=75, x2=355, y2=75, stroke="black", marker_end="url(#arrow)")
	text(x=370, y=80, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("\"Gouda\"")
	end
	text(x=300, y=150, font_family="JuliaMono, monospace", text_anchor="end", font_size="0.85rem", font_weight=600) do
		str("numbers ->")
	end
    rect(x=310, y=105, width=160, height=85, fill="rgb(242, 242, 242)", stroke="black")
	text(x=320, y=125, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("1")
	end
	line(x1=335, y1=120, x2=355, y2=120, stroke="black", marker_end="url(#arrow)")
	text(x=370, y=125, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("42")
	end
	text(x=320, y=165, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("2")
	end
	text(x=370, y=150, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("123")
	end
	line(x1=335, y1=158, x2=355, y2=150, stroke="black", stroke_dasharray="5,5", marker_end="url(#arrow)")
	text(x=370, y=180, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("5")
	end
	line(x1=335, y1=162, x2=355, y2=170, stroke="black", marker_end="url(#arrow)")
	text(x=300, y=215, font_family="JuliaMono, monospace", text_anchor="end", font_size="0.85rem", font_weight=600) do
		str("empty ->")
	end
    rect(x=310, y=200, width=160, height=20, fill="rgb(242, 242, 242)", stroke="black")
end

# ╔═╡ c41c5b10-7ebb-11eb-06f6-57f0ccbbf9ac
md"*Figure 10-1. State diagram.*"

# ╔═╡ d3894a90-7ebb-11eb-0541-133459689a2d
md"""Arrays are represented by boxes with the elements of the array inside them. `cheeses` refers to an array with three elements indexed: 1, 2, and 3. `numbers` contains two elements; the diagram shows that the value of the second element has been reassigned from `123` to `5`. `empty` refers to an array with no elements.

Array indices work the same way as string indices (but without UTF-8 caveats):

- Any integer expression can be used as an index.
- If you try to read or write an element that does not exist, you get a BoundsError.
- The keyword end points to the last index of the array.

The `∈` operator also works on arrays:

```julia
julia> "Edam" ∈ cheeses
true
julia> "Brie" in cheeses
false
```
"""

# ╔═╡ 1628711e-7ebc-11eb-2cec-f71c338de835
md"""## Traversing an Array

The most common way to traverse the elements of an array is with a `for` loop. The syntax is the same as for strings:

```julia
for cheese in cheeses
	println(cheese)
end
```

This works well if you only need to read the elements of the array. But if you want to write or update the elements, you need the indices. One technique is to use the built-in function `eachindex`:

```julia
for i in eachindex(numbers)
	numbers[i] = numbers[i] * 2
end
```

This loop traverses the array and updates each element. Each time through the loop `i` gets the index of the next element. The assignment statement in the body uses `i` to read the old value of the element and to assign the new value.

A `for` loop over an empty array never runs the body:

```julia
for x in []
	println("This can never happen.")
end
```

Although an array can contain another array, the nested array still counts as a single element. The length of this array is 4:

```julia
["spam", 1, ["Brie", "Roquefort", "Camembert"], [1, 2, 3]]
```
"""

# ╔═╡ 9c7423e4-7ebc-11eb-391d-b3388978ea47
md"""## Array Slices

The slice operator (`[n:m]`) also works on arrays:

```julia
julia> t = ['a', 'b', 'c', 'd', 'e', 'f'];

julia> print(t[1:3])
['a', 'b', 'c']
julia> print(t[3:end])
['c', 'd', 'e', 'f']
```

Using the slice operator without arguments makes a copy of the whole array:

```julia
julia> print(t[:])
['a', 'b', 'c', 'd', 'e', 'f']
```

Since arrays are mutable, it is often useful to make a copy before performing operations that modify arrays.

A slice operator on the left side of an assignment can update multiple elements:

```julia
julia> t[2:3] = ['x', 'y'];

julia> print(t)
['a', 'x', 'y', 'd', 'e', 'f']
```
"""

# ╔═╡ eadf298e-7ebc-11eb-26e8-65a4f6640fbd
md"""## Array Library

Julia provides functions that operate on arrays. For example, `push!` adds a new element to the end of an array:

```julia
julia> t = ['a', 'b', 'c'];

julia> push!(t, 'd');

julia> print(t) ['a', 'b', 'c', 'd']
```

`append!` adds the elements of the second array to the end of the first:

```julia
julia> t1 = ['a', 'b', 'c'];

julia> t2 = ['d', 'e'];

julia> append!(t1, t2);

julia> print(t1)
['a', 'b', 'c', 'd', 'e']
```

This example leaves `t2` unmodified.

`sort!` arranges the elements of the array from low to high:

```julia
julia> t = ['d', 'c', 'e', 'b', 'a'];

julia> sort!(t);

julia> print(t)
['a', 'b', 'c', 'd', 'e']
```

`sort` returns a copy of the elements of the array in order:

```julia
julia> t1 = ['d', 'c', 'e', 'b', 'a'];

julia> t2 = sort(t1);

julia> print(t1)
['d', 'c', 'e', 'b', 'a']
julia> print(t2)
['a', 'b', 'c', 'd', 'e']
```

!!! tip
    As a style convention in Julia, `!` is appended to names of functions that modify their arguments.
"""

# ╔═╡ 5b124fc4-7ebd-11eb-1a06-c9ea65c0953c
md"""## Map, Filter, and Reduce

To add up all the numbers in an array, you can use a loop like this:

```julia
function addall(t)
	total = 0
	for x in t
		total += x
	end
	total
end
```

`total` is initialized to 0. Each time through the loop, `x` gets one element from the array. The `+=` operator provides a short way to update a variable. This *augmented assignment statement*:

```julia
total += x
```

is equivalent to:

```julia
total = total + x
```

As the loop runs, `total` accumulates the sum of the elements; a variable used this way is sometimes called an *accumulator*.

Adding up the elements of an array is such a common operation that Julia provides it
as a built-in function, `sum`:

```julia
julia>t=[1,2,3,4];

julia> sum(t)
10
```

An operation like this that combines a sequence of elements into a single value is sometimes called a *reduce operation*.

Often you want to traverse one array while building another. For example, the following function takes an array of strings and returns a new array that contains capitalized strings:

```julia
function capitalizeall(t)
	res=[]
	for s in t
		push!(res, uppercase(s))
	end
	res
end
```

`res` is initialized with an empty array; each time through the loop, we append the next element. So, `res` is another kind of accumulator.

An operation like `capitalizeall` is sometimes called a *map* because it “maps” a function (in this case `uppercase`) onto each of the elements in a sequence.

Another common operation is to select some of the elements from an array and return a subarray. For example, the following function takes an array of strings and returns an array that contains only the uppercase strings:

```julia
function onlyupper(t)
	res=[]
	for s int
		if s == uppercase(s)
			push!(res, s)
		end
	end
	res
end
```

An operation like `onlyupper` is called a *filter* because it selects some of the elements and filters out the others.

Most common array operations can be expressed as a combination of map, filter, and reduce.
"""

# ╔═╡ 8b9a193c-7ebe-11eb-3fc3-b76b75cd0186
md"""## Dot Syntax

For every binary operator, like `^`, there is a corresponding *dot operator*, like `.^`, that is automatically defined to perform the operation element-by-element on arrays. For example, `[1, 2, 3] ^ 3` is not defined, but `[1, 2, 3] .^ 3` is defined as computing the elementwise result `[1^3, 2^3, 3^3]`:

```julia
julia> print([1, 2, 3] .^ 3)
[1, 8, 27]
```

Any Julia function `f` can be applied elementwise to any array with the *dot syntax*. For example, to capitalize an array of strings, we don’t need an explicit loop:

```julia
julia> t = uppercase.(["abc", "def", "ghi"]);

julia> print(t)
["ABC", "DEF", "GHI"]
```

This is an elegant way to create a map. The function `capitalizeall` can be imple‐ mented by a one-liner:

```julia
capitalizeall(t) = uppercase.(t)
```
"""

# ╔═╡ e1966318-7ebe-11eb-37a1-775a6c3621aa
md"""## Deleting (Inserting) Elements

There are several ways to delete elements from an array. If you know the index of the element you want, you can use `splice!`:

```julia
julia> t = ['a', 'b', 'c'];

julia> splice!(t, 2)
'b': ASCII/Unicode U+0062 (category Ll: Letter, lowercase)
julia> print(t)
['a', 'c']
```

`splice!` modifies the array and returns the element that was removed.

`pop!` deletes and returns the last element:

```julia
julia> t = ['a', 'b', 'c'];

julia> pop!(t)
'c': ASCII/Unicode U+0063 (category Ll: Letter, lowercase)
julia> print(t)
['a', 'b']
```

`popfirst!` deletes and returns the first element:

```julia
julia> t = ['a', 'b', 'c'];

julia> popfirst!(t)
'a': ASCII/Unicode U+0061 (category Ll: Letter, lowercase)
julia> print(t)
['b', 'c']
```

The functions `pushfirst!` and `push!` insert an element at the beginning and end, respectively, of the array.

If you don’t need the removed value, you can use the function `deleteat!`:

```julia
julia> t = ['a', 'b', 'c'];

julia> print(deleteat!(t, 2))
['a', 'c']
```

The function `insert!` inserts an element at a given index:

```julia
julia> t = ['a', 'b', 'c'];

julia> print(insert!(t, 2, 'x'))
['a', 'x', 'b', 'c']
```
"""

# ╔═╡ 53f60b52-7ebf-11eb-140a-75c0c412b628
md"""## Arrays and Strings

A string is a sequence of characters and an array is a sequence of values, but an array of characters is not the same as a string. To convert from a string to an array of characters, you can use the function `collect`:

```julia
julia> t = collect("spam");

julia> print(t)
['s', 'p', 'a', 'm']
```

The collect function breaks a string or another sequence into individual elements.

If you want to break a string into words, you can use the `split` function:

```julia
julia> t = split("pining for the fjords");

julia> print(t)
SubString{String}["pining", "for", "the", "fjords"]
```

An *optional argument* called a *delimiter* specifies which characters to use as word boundaries (we talked briefly about optional arguments in “Exercise 8-7”). The following example uses a hyphen as a delimiter:

```julia
julia> t = split("spam-spam-spam", '-');

julia> print(t)
SubString{String}["spam", "spam", "spam"]
```

`join` is the inverse of `split`. It takes an array of strings and concatenates the elements:

```julia
julia> t = ["pining", "for", "the", "fjords"];

julia> s = join(t, ' ')
"pining for the fjords"
```

In this case the delimiter is a space character. To concatenate strings without spaces, you don’t specify a delimiter.
"""

# ╔═╡ e08fed44-7ebf-11eb-1c6a-55ee8e49968e
md"""## Objects and Values

An *object* is something a variable can refer to. Until now, you might have used “object” and “value” interchangeably.

If we run these assignment statements:

```julia
a = "banana"
b = "banana"
```

we know that `a` and `b` both refer to a string, but we don’t know whether they refer to the same string. There are two possible states, as shown in Figure 10-2.
"""

# ╔═╡ 04e15ba6-7ec0-11eb-3c7d-211e84ba32c0
Drawing(width=720, height=70) do
	defs() do
        marker(id="arrow", markerWidth="10", markerHeight="10", refX="0", refY="3", orient="auto", markerUnits="strokeWidth") do
      		path(d="M0,0 L0,6 L9,3 z", fill="black")
		end
	end
    rect(x=210, y=10, width=140, height=50, fill="rgb(242, 242, 242)", stroke="black")
	rect(x=360, y=10, width=140, height=50, fill="rgb(242, 242, 242)", stroke="black")
	text(x=220, y=30, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("a")
	end
	line(x1=235, y1=25, x2=255, y2=25, stroke="black", marker_end="url(#arrow)")
	text(x=270, y=30, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("\"banana\"")
	end
	text(x=220, y=50, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("b")
	end
	line(x1=235, y1=45, x2=255, y2=45, stroke="black", marker_end="url(#arrow)")
	text(x=270, y=50, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("\"banana\"")
	end
	text(x=370, y=30, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("a")
	end
	line(x1=385, y1=25, x2=405, y2=30, stroke="black", marker_end="url(#arrow)")
	text(x=420, y=40, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("\"banana\"")
	end
	text(x=370, y=50, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("b")
	end
	line(x1=385, y1=45, x2=405, y2=40, stroke="black", marker_end="url(#arrow)")
end

# ╔═╡ a2a6face-7ec2-11eb-313a-73f3916a14ad
md"*Figure 10-2. State diagrams.*"

# ╔═╡ b218e204-7ec2-11eb-3321-3f0082f5e74d
md"""In one case, `a` and `b` refer to two different objects that have the same value. In the second case, they refer to the same object.

To check whether two variables refer to the same object, you can use the `≡` (**`\equiv TAB`**) or `===` operator:

```julia
julia> a = "banana";

julia> b = "banana";

julia> a ≡ b
true
```

In this example, Julia only created one string object, and both `a` and `b` refer to it. But when you create two arrays, you get two objects:

```julia
julia> a = [1, 2, 3];

julia> b = [1, 2, 3];

julia> a ≡ b false
```

So the state diagram looks like Figure 10-3.
"""

# ╔═╡ 91741fb0-7ec3-11eb-1c35-03f99160c75d
Drawing(width=720, height=70) do
	defs() do
        marker(id="arrow", markerWidth="10", markerHeight="10", refX="0", refY="3", orient="auto", markerUnits="strokeWidth") do
      		path(d="M0,0 L0,6 L9,3 z", fill="black")
		end
	end
    rect(x=290, y=10, width=140, height=50, fill="rgb(242, 242, 242)", stroke="black")
	text(x=300, y=30, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("a")
	end
	line(x1=315, y1=25, x2=335, y2=25, stroke="black", marker_end="url(#arrow)")
	text(x=350, y=30, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("[1, 2, 3]")
	end
	text(x=300, y=50, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("b")
	end
	line(x1=315, y1=45, x2=335, y2=45, stroke="black", marker_end="url(#arrow)")
	text(x=350, y=50, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("[1, 2, 3]")
	end
end

# ╔═╡ fec78276-7ec3-11eb-2876-ebd0e582766f
md"*Figure 10-3. State diagram.*"

# ╔═╡ 0d8047dc-7ec4-11eb-2b0c-91f0a3c01be7
md"""In this case we would say that the two arrays are equivalent, because they have the same elements, but not *identical*, because they are not the same object. If two objects are identical, they are also equivalent, but if they are equivalent, they are not necessarily identical.

To be precise, an object has a value. If you evaluate `[1, 2, 3]`, you get an array object whose value is a sequence of integers. If another array has the same elements, we say it has the same value, but it is not the same object."""

# ╔═╡ 32e15b56-7ec4-11eb-11d1-d576f5f7e355
md"""## Aliasing

If `a` refers to an object and you assign `b = a`, then both variables refer to the same object:

```julia
julia> a = [1, 2, 3];

julia> b = a;

julia> b ≡ a
true
```

The state diagram looks like Figure 10-4.
"""

# ╔═╡ 4ff0e338-7ec4-11eb-15ff-a17321526a29
Drawing(width=720, height=70) do
	defs() do
        marker(id="arrow", markerWidth="10", markerHeight="10", refX="0", refY="3", orient="auto", markerUnits="strokeWidth") do
      		path(d="M0,0 L0,6 L9,3 z", fill="black")
		end
	end
    rect(x=290, y=10, width=140, height=50, fill="rgb(242, 242, 242)", stroke="black")
	text(x=300, y=30, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("a")
	end
	line(x1=315, y1=25, x2=335, y2=30, stroke="black", marker_end="url(#arrow)")
	text(x=350, y=40, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("[1, 2, 3]")
	end
	text(x=300, y=50, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("b")
	end
	line(x1=315, y1=45, x2=335, y2=40, stroke="black", marker_end="url(#arrow)")
end

# ╔═╡ 85204788-7ec4-11eb-2337-c7f0728c2c5b
md"""The association of a variable with an object is called a *reference*. In this example, there are two references to the same object.

An object with more than one reference has more than one name, so we say that the object is *aliased*.

If the aliased object is mutable, changes made with one alias affect the other:

```julia
julia> b[1] = 42
42
julia> print(a)
[42, 2, 3]
```

!!! danger
    Although this behavior can be useful, it is error-prone. In general, it is safer to avoid aliasing when you are working with mutable objects.

For immutable objects like strings, aliasing is not as much of a problem. In this example:

```julia
a = "banana"
b = "banana"
```

it almost never makes a difference whether `a` and `b` refer to the same string or not.
"""

# ╔═╡ d9c3898a-7ec4-11eb-16de-a547447e325e
md"""## Array Arguments

When you pass an array to a function, the function gets a reference to the array. If the function modifies the array, the caller sees the change. For example, `deletehead!` removes the first element from an array:

```julia
deletehead!(t) = popfirst!(t)
```

Here’s how it is used:

```julia
julia> letters = ['a', 'b', 'c'];

julia> deletehead!(letters);

julia> print(letters)
['b', 'c']
```

The parameter `t` and the variable `letters` are aliases for the same object. The stack diagram looks like Figure 10-5.
"""

# ╔═╡ 11e91384-7ec5-11eb-325d-859e008c7cb1
Drawing(width=720, height=90) do
	defs() do
        marker(id="arrow", markerWidth="10", markerHeight="10", refX="0", refY="3", orient="auto", markerUnits="strokeWidth") do
      		path(d="M0,0 L0,6 L9,3 z", fill="black")
		end
	end
	text(x=140, y=30, font_family="JuliaMono, monospace", text_anchor="end", font_size="0.85rem", font_weight=600) do
		str("Main")
	end
    rect(x=170, y=10, width=200, height=30, fill="rgb(242, 242, 242)", stroke="black")
	text(x=180, y=30, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("letters")
	end
	line(x1=250, y1=25, x2=390, y2=25, stroke="black", marker_end="url(#arrow)")
	text(x=140, y=70, font_family="JuliaMono, monospace", text_anchor="end", font_size="0.85rem", font_weight=600) do
		str("deletehead!")
	end
    rect(x=170, y=50, width=200, height=30, fill="rgb(242, 242, 242)", stroke="black")
	text(x=180, y=70, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("t")
	end
	line(x1=200, y1=65, x2=390, y2=65, stroke="black", marker_end="url(#arrow)")
	rect(x=400, y=10, width=90, height=70, fill="rgb(242, 242, 242)", stroke="black")
	text(x=410, y=30, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("1")
	end
	line(x1=425, y1=25, x2=440, y2=25, stroke="black", marker_end="url(#arrow)")
	text(x=455, y=30, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("'a'")
	end
	text(x=410, y=50, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("2")
	end
	line(x1=425, y1=45, x2=440, y2=45, stroke="black", marker_end="url(#arrow)")
	text(x=455, y=50, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("'b'")
	end
	text(x=410, y=70, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("3")
	end
	line(x1=425, y1=65, x2=440, y2=65, stroke="black", marker_end="url(#arrow)")
	text(x=455, y=70, font_family="JuliaMono, monospace", font_size="0.85rem", font_weight=600) do
		str("'c'")
	end
end

# ╔═╡ e4ab98a4-7ec6-11eb-2293-2f901f40ffb0
md"*Figure 10-5. Stack diagram.*"

# ╔═╡ f68d4f86-7ec6-11eb-167b-cbf4bac179d6
md"""Since the array is shared by two frames, I drew it between them.

It is important to distinguish between operations that modify arrays and operations that create new arrays. For example, push! modifies an array, but vcat creates a new array.

Here’s an example using `push!`:

```julia
julia> t1 = [1, 2];

julia> t2 = push!(t1, 3);

julia> print(t1)
[1, 2, 3]
```

`t2` is an alias of `t1`.

Here’s an example using `vcat`:

```julia
julia> t3 = vcat(t1, [4]);

julia> print(t1)
[1, 2, 3]
julia> print(t3)
[1, 2, 3, 4]
```

The result of `vcat` is a new array, and the original array is unchanged.

This difference is important when you write functions that are supposed to modify
arrays.

For example, this function *does not* delete the head of an array:

```julia
function baddeletehead(t)
	t = t[2:end]               # WRONG!
end
```

The slice operator creates a new array and the assignment makes `t` refer to it, but that doesn’t affect the caller:

```julia
julia> t4 = baddeletehead(t3);

julia> print(t3)
[1, 2, 3, 4]
julia> print(t4)
[2, 3, 4]
```

At the beginning of `baddeletehead`, `t` and `t3` refer to the same array. At the end, `t` refers to a new array, but `t3` still refers to the original unmodified array.

An alternative is to write a function that creates and returns a new array. For example, `tail` returns all but the first element of an array:

```julia
tail(t) = t[2:end]
```

This function leaves the original array unmodified. Here’s how it is used:

```julia
julia> letters = ['a', 'b', 'c'];

julia> rest = tail(letters);

julia> print(rest)
['b', 'c']
```
"""

# ╔═╡ 96e036b0-7ec7-11eb-0308-1d445f8d21fa
md"""## Debugging

Careless use of arrays (and other mutable objects) can lead to long hours of debugging. Here are some common pitfalls and ways to avoid them:

* Most array functions modify the argument. This is the opposite of the string functions, which return a new string and leave the original alone.

  If you are used to writing string code like this:

  ```julia
  new_word = strip(word)
  ```

  It is tempting to write array code like this:

  ```julia
  t2 = sort!(t1)
  ```

  Because `sort!` returns the modified original array `t1`, `t2` is an alias of `t1`.

  !!! tip
      Before using array functions and operators, you should read the documentation carefully and then test them in interactive mode.

* Pick an idiom and stick with it.

  Part of the problem with arrays is that there are too many ways to do things. For example, to remove an element from an array, you can use `pop!`, `popfirst!`, `delete_at`, or even a slice assignment. To add an element, you can use `push!`, `pushfirst!`, `insert!`, or `vcat`. Assuming that `t` is an array and `x` is an array element, these are correct:

  ```julia
  insert!(t, 4, x)
  push!(t, x)
  append!(t, [x])
  ```

  And these are wrong:

  ```julia
  insert!(t, 4, [x])     # WRONG!
  push!(t, [x])          # WRONG!
  vcat(t, [x])           # WRONG!
  ```

* Make copies to avoid aliasing.

  If you want to use a function like sort! that modifies the argument, but you need to keep the original array as well, you can make a copy:

  ```julia
  julia> t = [3, 1, 2];

  julia> t2 = t[:]; # t2 = copy(t)

  julia> sort!(t2);

  julia> print(t)
  [3, 1, 2]
  julia> print(t2)
  [1, 2, 3]
  ```

  In this example you could also use the built-in function `sort`, which returns a new sorted array and leaves the original alone:

  ```julia
  julia> t2 = sort(t);

  julia> println(t)
  [3, 1, 2]
  julia> println(t2)
  [1, 2, 3]
  ```
"""

# ╔═╡ eb6f80ea-7ec8-11eb-1b60-916222e958db
md"""## Glossary

*array*:
A sequence of values.

*element*:
One of the values in an array (or other sequence); also called items.

*nested array*:
An array that is an element of another array.

*mutable*:
The property of a value that can be modified.

*augmented assignment*:
A statement that updates the value of a variable using an operator like `=.`

*accumulator*:
A variable used in a loop to add up or accumulate a result.

*reduce operation*:
A processing pattern that traverses a sequence and accumulates the elements into a single result.

*map*:
A processing pattern that traverses a sequence and performs an operation on each element.

*filter*:
A processing pattern that traverses a sequence and selects the elements that sat‐ isfy some criterion.

*dot operator*:
A binary operator that is applied element-by-element to arrays.

*dot syntax*:
Syntax used to apply a function elementwise to any array.

*optional argument*:
An argument that is not required.

*delimiter*:
A character or string used to indicate where a string should be split.

*object*:
Something a variable can refer to. An object has a type and a value.

*equivalent*:
Having the same value.

*identical*:
Being the same object (which implies equivalence).

*reference*:
The association between a variable and its value.

*aliasing*:
A circumstance where two or more variables refer to the same object.
"""

# ╔═╡ 4ac0178a-7ec9-11eb-2e69-51752445551a
md"""## Exercises

#### Exercise 10-1

Write a function called `nestedsum` that takes an array of arrays of integers and adds up the elements from all of the nested arrays. For example:

```julia
julia> t = [[1, 2], [3], [4, 5, 6]];

julia> nestedsum(t)
21
```
"""

# ╔═╡ 6c3c929e-7ec9-11eb-00e3-87c9c28c60bc
md"""#### Exercise 10-2

Write a function called `cumulsum` that takes an array of numbers and returns the cumulative sum; that is, a new array where the ith element is the sum of the first ``i`` elements from the original array. For example:

```julia
julia> t = [1, 2, 3];

julia> print(cumulsum(t))
Any[1, 3, 6]
```
"""

# ╔═╡ a81cc482-7ec9-11eb-018d-f933e7d86b01
md"""#### Exercise 10-3

Write a function called `interior` that takes an array and returns a new array that contains all but the first and last elements. For example:

```julia
julia>t=[1,2,3,4];

julia> print(interior(t))
[2, 3]
```
"""

# ╔═╡ c1a1c43e-7ec9-11eb-3454-6d078f7a70de
md"""#### Exercise 10-4

Write a function called `interior!` that takes an array, modifies it by removing the first and last elements, and returns `nothing`. For example:

```julia
julia> t=[1,2,3,4];

julia> interior!(t)

julia> print(t)
[2, 3]
```
"""

# ╔═╡ dd4164ce-7ec9-11eb-0504-e3e6966b861b
md"""#### Exercise 10-5

Write a function called `issort` that takes an array as a parameter and returns `true` if the array is sorted in ascending order and `false` otherwise. For example:

```julia
julia> issort([1, 2, 2])
true
julia> issort(['b', 'a'])
false
```
"""

# ╔═╡ 0393e21e-7eca-11eb-1d46-a394c6c31f3f
md"""#### Exercise 10-6

Two words are anagrams if you can rearrange the letters from one to spell the other. Write a function called `isanagram` that takes two strings and returns `true` if they are anagrams.
"""

# ╔═╡ 18bd9022-7eca-11eb-0f83-43fc924505a1
md"""#### Exercise 10-7

Write a function called `hasduplicates` that takes an array and returns `true` if there is any element that appears more than once. It should not modify the original array.
"""

# ╔═╡ 33230758-7eca-11eb-34b6-6b89f1c8afd4
md"""#### Exercise 10-8

This exercise pertains to the so-called Birthday Paradox.

If there are 23 students in your class, what are the chances that 2 of you have the same birthday? You can estimate this probability by generating random samples of 23 birthdays and checking for matches.

!!! tip
    You can generate random birthdays with `rand(1:365)`.
"""

# ╔═╡ 4db45f7c-7eca-11eb-0d84-697b0f34ba19
md"""#### Exercise 10-9

Write two versions of a function that reads the file *words.txt* and builds an array with one element per word, one using `push!` and the other using the idiom `t = [t..., x]`. Which one takes longer to run? Why?
"""

# ╔═╡ 644ba7a4-7eca-11eb-005c-450b665506f8
md"""#### Exercise 10-10

To check whether a word is in the word array you just built, you could use the `∈` operator, but it would be slow because it searches through the words in order.

Because the words are in alphabetical order, we can speed things up with a bisection search (also known as a binary search), which is similar to what you do when you look a word up in the dictionary. You start in the middle and check to see whether the word you are looking for comes before the word in the middle of the array. If so, you search the first half of the array the same way. Otherwise, you search the second half.

Either way, you cut the remaining search space in half. If the word array has 113,809 words, it will take about 17 steps to find the word or conclude that it’s not there.

Write a function called `inbisect` that takes a sorted array and a target value and returns `true` if the word is in the array and `false` if it’s not.
"""

# ╔═╡ 9b0559d6-7eca-11eb-0e90-b35daf751f26
md"""#### Exercise 10-11

Two words are a “reverse pair” if each is the reverse of the other. Write a function `reversepairs` that finds all the reverse pairs in the word array.
"""

# ╔═╡ b70a0e20-7eca-11eb-2181-a3b160b3c369
md"""#### Exercise 10-12

Two words “interlock” if taking alternating letters from each forms a new word. For example, “shoe” and “cold” interlock to form “schooled.”

1. Write a program that finds all pairs of words that interlock.

   !!! tip
       Don’t enumerate all pairs!

   *Credit*: This exercise is inspired by an example at [http://puzzlers.org](http://puzzlers.org).

2. Can you find any words that are three-way interlocked (i.e., every third letter forms a word, starting from the first, second, or third letter)?
"""

# ╔═╡ Cell order:
# ╟─54f89e42-7eb6-11eb-0a18-79732e59d627
# ╟─21cefe9c-7eb3-11eb-2469-2b7b3d1eccc3
# ╟─62637476-7eb6-11eb-0cb4-01f126c2feb8
# ╟─daca2c98-7eb6-11eb-1c09-8f1e08b31315
# ╟─1b70a1d2-7eb7-11eb-2ab9-6739d60fcf89
# ╟─c41c5b10-7ebb-11eb-06f6-57f0ccbbf9ac
# ╟─d3894a90-7ebb-11eb-0541-133459689a2d
# ╟─1628711e-7ebc-11eb-2cec-f71c338de835
# ╟─9c7423e4-7ebc-11eb-391d-b3388978ea47
# ╟─eadf298e-7ebc-11eb-26e8-65a4f6640fbd
# ╟─5b124fc4-7ebd-11eb-1a06-c9ea65c0953c
# ╟─8b9a193c-7ebe-11eb-3fc3-b76b75cd0186
# ╟─e1966318-7ebe-11eb-37a1-775a6c3621aa
# ╟─53f60b52-7ebf-11eb-140a-75c0c412b628
# ╟─e08fed44-7ebf-11eb-1c6a-55ee8e49968e
# ╟─04e15ba6-7ec0-11eb-3c7d-211e84ba32c0
# ╟─a2a6face-7ec2-11eb-313a-73f3916a14ad
# ╟─b218e204-7ec2-11eb-3321-3f0082f5e74d
# ╟─91741fb0-7ec3-11eb-1c35-03f99160c75d
# ╟─fec78276-7ec3-11eb-2876-ebd0e582766f
# ╟─0d8047dc-7ec4-11eb-2b0c-91f0a3c01be7
# ╟─32e15b56-7ec4-11eb-11d1-d576f5f7e355
# ╟─4ff0e338-7ec4-11eb-15ff-a17321526a29
# ╟─85204788-7ec4-11eb-2337-c7f0728c2c5b
# ╟─d9c3898a-7ec4-11eb-16de-a547447e325e
# ╟─11e91384-7ec5-11eb-325d-859e008c7cb1
# ╟─e4ab98a4-7ec6-11eb-2293-2f901f40ffb0
# ╟─f68d4f86-7ec6-11eb-167b-cbf4bac179d6
# ╟─96e036b0-7ec7-11eb-0308-1d445f8d21fa
# ╟─eb6f80ea-7ec8-11eb-1b60-916222e958db
# ╟─4ac0178a-7ec9-11eb-2e69-51752445551a
# ╟─6c3c929e-7ec9-11eb-00e3-87c9c28c60bc
# ╟─a81cc482-7ec9-11eb-018d-f933e7d86b01
# ╟─c1a1c43e-7ec9-11eb-3454-6d078f7a70de
# ╟─dd4164ce-7ec9-11eb-0504-e3e6966b861b
# ╟─0393e21e-7eca-11eb-1d46-a394c6c31f3f
# ╟─18bd9022-7eca-11eb-0f83-43fc924505a1
# ╟─33230758-7eca-11eb-34b6-6b89f1c8afd4
# ╟─4db45f7c-7eca-11eb-0d84-697b0f34ba19
# ╟─644ba7a4-7eca-11eb-005c-450b665506f8
# ╟─9b0559d6-7eca-11eb-0e90-b35daf751f26
# ╟─b70a0e20-7eca-11eb-2181-a3b160b3c369


text:
"""
    struct ProductSpace{S<:ElementarySpace, N} <: CompositeSpace{S}

A `ProductSpace` is a tensor product space of `N` vector spaces of type
`S<:ElementarySpace`. Only tensor products between [`ElementarySpace`](@ref) objects of the
same type are allowed.
"""
struct ProductSpace{S<:ElementarySpace, N} <: CompositeSpace{S}
    spaces::NTuple{N, S}
end
ProductSpace(spaces::Vararg{S,N}) where {S<:ElementarySpace, N} =
    ProductSpace{S,N}(spaces)
ProductSpace{S,N}(spaces::Vararg{S,N}) where {S<:ElementarySpace, N} =
    ProductSpace{S,N}(spaces)

# Corresponding methods
#-----------------------
"""
    dims(::ProductSpace{S,N}) -> Dims{N} = NTuple{N,Int}

Return the dimensions of the spaces in the tensor product space as a tuple of integers.
"""
dims(P::ProductSpace) = map(dim, P.spaces)
dim(P::ProductSpace, n::Int) = dim(P.spaces[n])
dim(P::ProductSpace) = prod(dims(P))

Base.axes(P::ProductSpace) = map(axes, P.spaces)
Base.axes(P::ProductSpace, n::Int) = axes(P.spaces[n])

dual(P::ProductSpace{<:ElementarySpace,0}) = P
dual(P::ProductSpace) = ProductSpace(map(dual, reverse(P.spaces)))

# Base.conj(P::ProductSpace) = ProductSpace(map(conj, P.spaces))

function Base.show(io::IO, P::ProductSpace{S}) where {S<:ElementarySpace}
    spaces = P.spaces
    if length(spaces) == 0
        print(io,"ProductSpace{", S, ",0}")
    end
    if length(spaces) == 1
        print(io,"ProductSpace")
    end
    print(io,"(")
    for i in 1:length(spaces)
        i==1 || print(io," ⊗ ")
        show(io, spaces[i])
    end
    print(io,")")
end

# more specific methods
sectors(P::ProductSpace) = _sectors(P, sectortype(P))
_sectors(P::ProductSpace{<:ElementarySpace, N}, ::Type{Trivial}) where {N} =
    (ntuple(n->Trivial(), StaticLength{N}()),) # speed up sectors for ungraded spaces
_sectors(P::ProductSpace{<:ElementarySpace, N}, ::Type{<:Sector}) where {N} =
    product(map(sectors, P.spaces)...)

hassector(V::ProductSpace{<:ElementarySpace,N}, s::NTuple{N}) where {N} =
    reduce(&, map(hassector, V.spaces, s); init = true)

dims(P::ProductSpace{<:ElementarySpace, N}, sector::NTuple{N,<:Sector}) where {N} =
    map(dim, P.spaces, sector)
dim(P::ProductSpace{<:ElementarySpace, N}, sector::NTuple{N,<:Sector}) where {N} =
    prod(dims(P, sector))

Base.axes(P::ProductSpace{<:ElementarySpace,N}, sectors::NTuple{N,<:Sector}) where {N} =
    map(axes, P.spaces, sectors)

"""
    dims(::ProductSpace{S,N}) -> Dims{N} = NTuple{N,Int}

Return the dimensions of the spaces in the tensor product space as a tuple of integers.
"""
function blocksectors(P::ProductSpace{S,N}) where {S,N}
    G = sectortype(S)
    if G == Trivial
        return (Trivial(),)
    end
    bs = Vector{G}()
    if N == 0
        push!(bs, one(G))
    elseif N == 1
        for s in sectors(P)
            push!(bs, first(s))
        end
    else
        for s in sectors(P)
            for c in ⊗(s...)
                if !(c in bs)
                    push!(bs, c)
                end
            end
        end
        # return foldl(union!, Set{G}(), (⊗(s...) for s in sectors(P)))
    end
    return bs
end
function blockdim(P::ProductSpace, c::Sector)
    sectortype(P) == typeof(c) || throw(SectorMismatch())
    d = 0
    for s in sectors(P)
        ds = dim(P, s)
        d += length(fusiontrees(s, c))*ds
    end
    return d
end

Base.:(==)(P1::ProductSpace, P2::ProductSpace) = (P1.spaces == P2.spaces)

# Default construction from product of spaces
#---------------------------------------------
⊗(V1::S, V2::S) where {S<:ElementarySpace}= ProductSpace((V1, V2))
⊗(P1::ProductSpace{S}, V2::S) where {S<:ElementarySpace} =
    ProductSpace(tuple(P1.spaces..., V2))
⊗(V1::S, P2::ProductSpace{S}) where {S<:ElementarySpace} =
    ProductSpace(tuple(V1, P2.spaces...))
⊗(P1::ProductSpace{S}, P2::ProductSpace{S}) where {S<:ElementarySpace} =
    ProductSpace(tuple(P1.spaces..., P2.spaces...))
⊗(P::ProductSpace{S,0}, ::ProductSpace{S,0}) where {S<:ElementarySpace} = P
⊗(P::ProductSpace{S}, ::ProductSpace{S,0}) where {S<:ElementarySpace} = P
⊗(::ProductSpace{S,0}, P::ProductSpace{S}) where {S<:ElementarySpace} = P
⊗(V::ElementarySpace) = ProductSpace((V,))
⊗(P::ProductSpace) = P

# unit element with respect to the monoidal structure of taking tensor products
"""
    one(::S) where {S<:ElementarySpace} -> ProductSpace{S,0}
    one(::ProductSpace{S}) where {S<:ElementarySpace} -> ProductSpace{S,0}

Return a tensor product of zero spaces of type `S`, i.e. this is the unit object under the
tensor product operation, such that `V ⊗ one(V) == V`.
"""
Base.one(V::VectorSpace) = one(typeof(V))
Base.one(::Type{<:ProductSpace{S}}) where {S<:ElementarySpace} = ProductSpace{S,0}(())
Base.one(::Type{S}) where {S<:ElementarySpace} = ProductSpace{S,0}(())

Base.convert(::Type{<:ProductSpace}, V::ElementarySpace) = ProductSpace((V,))
Base.literal_pow(::typeof(^), V::ElementarySpace, p::Val{N}) where N =
    ProductSpace{typeof(V), N}(ntuple(n->V, p))
Base.convert(::Type{S}, P::ProductSpace{S,0}) where {S<:ElementarySpace} = oneunit(S)
Base.convert(::Type{S}, P::ProductSpace{S}) where {S<:ElementarySpace} = fuse(P.spaces...)
fuse(P::ProductSpace{S,0}) where {S<:ElementarySpace} = oneunit(S)
fuse(P::ProductSpace{S}) where {S<:ElementarySpace} = fuse(P.spaces...)

# Functionality for extracting and iterating over spaces
#--------------------------------------------------------
Base.length(P::ProductSpace) = length(P.spaces)
Base.getindex(P::ProductSpace, n::Integer) = P.spaces[n]
Base.getindex(P::ProductSpace{S}, I::NTuple{N,Integer}) where {S<:ElementarySpace,N} =
    ProductSpace{S,N}(TupleTools.getindices(P.spaces, I))

Base.iterate(P::ProductSpace) = Base.iterate(P.spaces)
Base.iterate(P::ProductSpace, s) = Base.iterate(P.spaces, s)

Base.eltype(P::ProductSpace{S}) where {S<:ElementarySpace} = S

Base.IteratorEltype(P::ProductSpace) = Base.IteratorEltype(P.spaces)
Base.IteratorSize(P::ProductSpace) = Base.IteratorSize(P.spaces)


text:
using EVQueues, Plots, ProgressMeter

lambda=120.0;
mu=1.0;
gamma=0.5;
#C=80.0;
C=60.0;

Tfinal=60.0;

frames = 24*30;

snaps = collect(range(0.01,stop=Tfinal,length=frames));

sim = ev_edf(lambda,mu,gamma,Tfinal,C,snapshots=snaps)
compute_statistics!(sim)

prog=Progress(length(snaps), dt=1, desc="Creando animacion... ");

azulcito = RGB(0.0,140/256,240/256);

rojito = RGB(240/256,0.0,140/256);

anim = @animate for i=1:length(snaps)

    p=stateplot(sim.snapshots[i].charging, markercolor=[azulcito rojito], xlims=(0,3/mu), ylims=(0,3*(1/mu+1/gamma)))

    xlabel!(p,"Carga remanente");
    ylabel!(p,"Tiempo remanente");
    title!(p,"Población EV - "*sim.parameters["Policy"])

    next!(prog);
end

gif(anim, "/tmp/"*sim.parameters["Policy"]*".gif", fps = 24)


text:
using Rcl

rcl.init()

nb = rcl.NodeBase("minimal_publisher")

# ret = Rcl.internal.rcl_node_fini(nb.rcl_node_handle)
# @show ret
# rcl.finalize(nb)
# finalize(nb)
finalize(nb)
# GC.gc()
# nb.rcl_node_handle[] = C_NULL
# ret = Rcl.internal.rcl_node_fini(nb.rcl_node_handle)
# @show ret

nb = nothing

rcl.shutdown()
# nb = nothing
GC.gc()

sleep(1.0)


text:
module Clutter

using ..Ahorn, Maple
using Random

const validClutterNames = String["yellowBlocks", "redBlocks", "greenBlocks"]
const clutterDisplayNames = Dict{String, String}(
    "yellowBlocks" => "Boxes",
    "redBlocks" => "Laundry",
    "greenBlocks" => "Books"
)
const clutterFunctions = Dict{String, Type}(
    "yellowBlocks" => Maple.YellowBlock,
    "redBlocks" => Maple.RedBlock,
    "greenBlocks" => Maple.GreenBlock
)
const simpleName = Dict{String, String}(
    "yellowBlocks" => "Yellow",
    "redBlocks" => "Red",
    "greenBlocks" => "Green"
)
const clutterBlockUnion = Union{
    Maple.YellowBlock,
    Maple.RedBlock,
    Maple.GreenBlock
}

const placements = Ahorn.PlacementDict(
    "Clutter Cabinet" => Ahorn.EntityPlacement(
        Maple.ClutterCabinet
    ),
)

for (raw, name) in clutterDisplayNames
    placements["Clutter Block ($name)"] = Ahorn.EntityPlacement(
        clutterFunctions[raw],
        "rectangle"
    )

    placements["Clutter Switch ($name)"] = Ahorn.EntityPlacement(
        Maple.ColorSwitch,
        "point",
        Dict{String, Any}(
            "type" => simpleName[raw]
        )
    )

    placements["Clutter Door ($name)"] = Ahorn.EntityPlacement(
        Maple.ClutterDoor,
        "point",
        Dict{String, Any}(
            "type" => simpleName[raw]
        )
    )
end

textureName(entity::clutterBlockUnion, i::Integer) = "objects/resortclutter/$(lowercase(simpleName[entity.name]))_$(lpad(i, 2, "0"))"

function getTextures(entity::clutterBlockUnion)
    i = 0
    res = Ahorn.Sprite[]

    atlas = Ahorn.getAtlas("Gameplay")

    while haskey(atlas, textureName(entity, i))
        push!(res, Ahorn.getSprite(textureName(entity, i), "Gameplay"))
        i += 1
    end

    return res
end

function getEntityRng(entity::clutterBlockUnion, dr::Ahorn.DrawableRoom)
    x = floor(Int, get(entity.data, "x", 0) / 8)
    y = floor(Int, get(entity.data, "y", 0) / 8)

    width = floor(Int, get(entity.data, "width", 32) / 8)
    height = floor(Int, get(entity.data, "height", 32) / 8)

    randStates = dr.bgTileStates.rands
    rands = get(randStates, (y:y + height - 1, x:x + width - 1), 1)

    # Rands are stored as 1:4, we need 0:3 here
    seed = abs(foldl((n, m) -> n << 2 | m, rands .- 1))

    return MersenneTwister(seed)
end

function renderClutterBlock(ctx::Ahorn.Cairo.CairoContext, entity::clutterBlockUnion, dr::Ahorn.DrawableRoom)
    entityX = Int(get(entity.data, "x", 0))
    entityY = Int(get(entity.data, "y", 0))

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    tw = ceil(Int, width / 8)
    th = ceil(Int, height / 8)

    needsDrawing = fill(true, (th, tw))

    rng = getEntityRng(entity, dr)
    textures = getTextures(entity)

    for y in 1:th, x in 1:tw
        if needsDrawing[y, x]
            choices = copy(textures)
            shuffle!(rng, choices)

            for choice in choices
                w, h = floor(Int, choice.width / 8), floor(Int, choice.height / 8)
                if all(get(needsDrawing, (y:y + h - 1, x:x + w - 1), false))
                    needsDrawing[y:y + h - 1, x:x + w - 1] .= false
                    Ahorn.drawImage(ctx, choice, entityX + x * 8 - 8, entityY + y * 8 - 8)

                    break
                end
            end
        end
    end
end

Ahorn.editingOptions(entity::Maple.ClutterDoor) = Dict{String, Any}(
    "type" => Maple.clutter_block_colors
)
Ahorn.editingOptions(entity::Maple.ColorSwitch) = Dict{String, Any}(
    "type" => Maple.clutter_block_colors
)

Ahorn.minimumSize(entity::clutterBlockUnion) = 8, 8
Ahorn.resizable(entity::clutterBlockUnion) = true, true

Ahorn.selection(entity::clutterBlockUnion) = Ahorn.getEntityRectangle(entity)

function Ahorn.selection(entity::Maple.ColorSwitch)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 4, y - 2, 40, 18)
end

function Ahorn.selection(entity::Maple.ClutterCabinet)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x, y, 16, 16)
end

function Ahorn.selection(entity::Maple.ClutterDoor)
    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 24))
    height = Int(get(entity.data, "height", 24))

    return Ahorn.Rectangle(x, y, width, height)
end

clutterDoorColor = (74, 71, 135, 0.6) ./ (255.0, 255.0, 255.0, 1.0)

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::clutterBlockUnion, room::Maple.Room)
    dr = Ahorn.getDrawableRoom(Ahorn.loadedState.map, room)

    renderClutterBlock(ctx, entity, dr)
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::ColorSwitch, room::Maple.Room)
    x, y = Ahorn.position(entity)

    variant = lowercase(get(entity.data, "type", "red"))
    sprite = Ahorn.getSprite("objects/resortclutter/icon_$variant", "Gameplay")

    Ahorn.drawImage(ctx, "objects/resortclutter/clutter_button00", x - 4, y - 2)
    Ahorn.drawImage(ctx, sprite, x + (32 - sprite.width) / 2, y + (16 - sprite.height) / 2)
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::ClutterCabinet, room::Maple.Room)
    x, y = Ahorn.position(entity)

    Ahorn.drawImage(ctx, "objects/resortclutter/cabinet00", x, y)
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::ClutterDoor, room::Maple.Room)
    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 24))
    height = Int(get(entity.data, "height", 24))

    variant = lowercase(get(entity.data, "type", "red"))
    sprite = Ahorn.getSprite("objects/resortclutter/icon_$variant", "Gameplay")

    Ahorn.drawRectangle(ctx, x, y, width, height, clutterDoorColor, (1.0, 1.0, 1.0, 8.0))
    Ahorn.drawImage(ctx, sprite, x + width / 2 - (sprite.width) / 2, y + height / 2 - (sprite.height) / 2)
end

end

max_stars_repo_name: asdfsdfadfasfasd/Ahorn

max_issues_repo_name: asdfsdfadfasfasd/Ahorn

max_forks_repo_name: asdfsdfadfasfasd/Ahorn


text:
function propagate_emit!(
    state,
    model,
    timestep::Timestep = RandomTimestep(Exponential(0.121))
)
    δ = get_step(timestep)
    propagate_emit!(state, model, δ)
end

function propagate_emit!(
    state::BinomialState{<:Array},
    model::AbstractBinomialModel,
    δ::Number
)
    if maximum(size(state.n)) > 1
        error("Emission not supported for non-scalar states")
        return nothing
    else
        propagate!(state, model, δ)
        q = model.q[1]
        σ = model.σ[1]
        k = state.k[1,1]
        EPSP = rand(Normal(q*k, σ))
        return BinomialObservation(EPSP, δ)
    end
end


text:
"""
$(TYPEDEF)

The supertype for Lattices in `D`-dimensional space. Type `T` is used to represent coordinates. The boundary conditions can be either periodic (`PB=true`) or free (`PB=false`).
"""
abstract type AbstractLattice{D, T, PB} <: AbstractNodeCollection{D, T} end

"""
$(TYPEDEF)
$(TYPEDFIELDS)

This type describes a regular arrangment of nodes in `D`-dimensional space with boundary condition `PB` (`PB=true` for periodic boundary conditions, `PB=false` for free periodic boundary conditions). A general lattice consists of identicall cells (combinations of nodes) arranged as a Bravais lattice.
"""
struct RegularLattice{D, T, PB, CT, L<:Union{Symbol,Nothing}} <: AbstractLattice{D, T, PB}
    """
    The number of basis cells along each of the basis directions.
    """
    lattice_dims::NTuple{D, Int}
    """
    Coordinates of the primitive vectors of the underlying Bravais lattice. `primitive_vecs[:, i]` gives the ``i``-th primitive vector.
    """
    primitive_vecs::SMatrix{D,D, T}
    """
    Repeated basis cell of the lattice.
    """
    basis_cell::CT
    """
    Label of the Lattice.
    """
    label::L
    """
    Total number of cells.
    """
    num_of_cells::Int
    """
    Total number of nodes.
    """
    num_of_nodes::Int
    """
    Cartesian index of the central cell (in the case of even length among of dimensions it is an approximation).
    """
    central_cell::CartesianIndex{D}
    """
    $(SIGNATURES)

    Constructs a Regular Lattice enforcing the condition that `typeof(basis_cell)<:AbstractCell{D,T}`.

    #Arguments
    - `lattice_dims::NTuple{D,Int}`: specifies the number of basis cells along each of the basis directions.
    - `primitive_vecs::SMatrix{D,D,T}`: stores as columns the primitive vectors of the underlying Bravais lattice.
    - `basis_cell::AbstractCell{D,T}`: specifies the repeated basis cell of the lattice; if the basis cell does not have any special structure, an instance of `TrivialCell` should be used.
    - `label::Union{Symbol,Nothing}`: the label of the lattice.
    - `PB::Bool`: specifies if periodic or free boundary conditions are applied; default value is `true`.
    """
	function RegularLattice(lattice_dims::NTuple{D, Int}, primitive_vecs::SMatrix{D,D,T}, basis_cell::AbstractCell{D,T}, label::Union{Symbol,Nothing}, PB::Bool) where {D,T}
        @assert D>=1
        num_of_cells = prod(lattice_dims)
        central_cell = CartesianIndex(div.(lattice_dims,2).+1)
        new{D, T, PB, typeof(basis_cell), typeof(label)}(lattice_dims, primitive_vecs, basis_cell, label, num_of_cells, length(basis_cell)*num_of_cells, central_cell)
    end
end

### Outer constructors.

function RegularLattice(lattice_dims::NTuple{D, Int}, primitive_vecs::SMatrix{D,D,T1}, basis_cell::AbstractCell{D,T2}, label, PB) where {D,T1,T2}
    T = promote_type(T1,T2)
    return RegularLattice(lattice_dims, SMatrix{D,D,T}(primitive_vecs), switch_coord_type(basis_cell, T), label, PB)
end

"""
`RegularLattice(lattice_dims::NTuple{D, Int}, a::T=1; label=:cubic, periodic=true)`

Constructs hypercubic lattice with lattice parameter `a` and trivial unit cell.
"""
function RegularLattice(lattice_dims::NTuple{D, Int}, a::T=1; periodic=true, label=:cubic) where {D, T<:Number}
    primitive_vecs = one(SMatrix{D,D,Int})*a
    basis_cell = TrivialCell{D,T}()
    return RegularLattice(lattice_dims, primitive_vecs, basis_cell, label, periodic)
end
function RegularLattice(lattice_dims::NTuple{D,Int}, primitive_vecs::SMatrix{D,D,T}; label=:simple, periodic = true) where {D,T}
    basis_cell = TrivialCell{D,T}()
    return RegularLattice(lattice_dims, primitive_vecs, basis_cell, label, periodic)
end
"""
`RegularLattice(lattice_dims::NTuple{D, Int}, primitive_vecs::SMatrix{D,D}, basis_cell::AbstractCell{D}; label=nothing, periodic=true)`

Convenient constructor which allows to specify the label and boundary condition as keyword arguments.
"""
RegularLattice(lattice_dims::NTuple{D, Int}, primitive_vecs::SMatrix{D,D}, basis_cell::AbstractCell{D}; label=nothing, periodic=true) where {D} =
    RegularLattice(lattice_dims, primitive_vecs, basis_cell, label, periodic)

"""
$(TYPEDSIGNATURES)

Return the number of lattice separate groups.
"""
@inline num_of_groups(lattice::RegularLattice) = num_of_groups(lattice.basis_cell)

Base.@propagate_inbounds group_size(lattice::RegularLattice, ig) = lattice.num_of_cells*group_size(lattice.basis_cell, ig)

Base.length(lattice::RegularLattice) = lattice.num_of_nodes

##
## Indexing interface.

"""
`getindex(lattice::RegularLattice{D}, Index...)`

Returns the coordinate of the node with index `Index`. The configurations of indices for different types of unit cell are the following:
- `TrivialCell`: `I::CartesianIndex{D}`.
- `HomogeneousCell`: `I::CartesianIndex{D}, ic::Int`.
- `InhomogeneousCell`: `I::CartesianIndex{D}, ic::Int, ig::Int`.
When iterating over lattice, the index at the left is iterated faster. For example, in the case of `InhomogeneousCell` I is iterated first, then `ic`, then `ig`.

In the case of periodic lattice, if `CartesianIndex` is outside of the lattice dims, it is simply translated back inside.
"""

Base.@propagate_inbounds function Base.getindex(lattice::RegularLattice{D,T}, I::CartesianIndex{D}, Ic...) where {D,T}
    @boundscheck check_lattice_index(lattice, I)
	@boundscheck check_cell_index(lattice.basis_cell, Ic...)
    @inbounds lattice.primitive_vecs*SVector{D}(I.I) + lattice.basis_cell[Ic...]
end
Base.@propagate_inbounds function Base.getindex(lattice::RegularLattice{D,T,false,<:TrivialCell}, I::CartesianIndex{D}) where {D,T}
    @boundscheck check_lattice_index(lattice, I)
    @inbounds lattice.primitive_vecs*SVector{D}(I.I)
end
Base.@propagate_inbounds function Base.getindex(lattice::RegularLattice{D,T,true}, I::CartesianIndex{D}, Ic...) where {D,T}
    @boundscheck check_cell_index(lattice.basis_cell, Ic...)
    @inbounds lattice.primitive_vecs*SVector{D}(mod1.(I.I, lattice.lattice_dims)) + lattice.basis_cell[Ic...]
end
Base.@propagate_inbounds function Base.getindex(lattice::RegularLattice{D,T,true,<:TrivialCell}, I::CartesianIndex{D}) where {D,T}
    @inbounds lattice.primitive_vecs*SVector{D}(mod1.(I.I, lattice.lattice_dims))
end

@inline check_lattice_index(lattice::RegularLattice{D}, I::CartesianIndex{D}) where {D} = _check_cartesian_index(true, I.I, lattice.lattice_dims) || throw(BoundsError(lattice, I))
@inline _check_cartesian_index(b, i, stop) = _check_cartesian_index(b & (1 <= i[1] <= stop[1]), Base.tail(i), Base.tail(stop))
@inline _check_cartesian_index(b, i::Tuple, ::Tuple{}) = false
@inline _check_cartesian_index(b, i::Tuple{}, stop::Tuple) = false
@inline _check_cartesian_index(b, i::Tuple{}, stop::Tuple{}) = b

## Relative coordinate

Base.@propagate_inbounds function relative_coordinate(lattice::RegularLattice{D,T,false}, I1::TI, I2::TI) where {D,T, TI<:Union{CartesianIndex{D}, Tuple{CartesianIndex{D},Vararg{Int}}}}
    return lattice[I1] - lattice[I2]
end
"""
$(TYPEDSIGNATURES)

For periodic lattice, the "shortest" relative coordinate is calculated instead.

For, that the following heuristic is used. Cartesian indices of the two nodes are shifted by the same amount, so that the cartesian index of the second node corresponds to the central cell of the lattice. Then, the cartesian index of the first node is translated back inside the lattice. The relative coordinate is computed using the resulting indices.

This heuristic guarantees that
`relative_coordinate(lattice, I1, I2) == -relative_coordinate(lattice, I2, I1)`.
"""
Base.@propagate_inbounds function relative_coordinate(lattice::RegularLattice{D,T,true,<:TrivialCell}, I1::TI, I2::TI) where {D,T, TI<:CartesianIndex{D}}
    return lattice[I1 + lattice.central_cell - I2] - lattice[lattice.central_cell]
end
Base.@propagate_inbounds function relative_coordinate(lattice::RegularLattice{D,T,true}, I1::TI, I2::TI) where {D,T, TI<:Tuple{CartesianIndex{D},Vararg{Int}}}
    return lattice[first(I1) + lattice.central_cell - first(I2), Base.tail(I1)...] - lattice[lattice.central_cell, Base.tail(I2)...]
end

#Base.@propagate_inbounds function translate_indices(lattice::RegularLattice{D,T,PB,<:TrivialCell}, Is) where {D,T,PB}
#    @boundscheck map(x -> check_lattice_index(lattice, x), Is)
#    return Is
#end
#Base.@propagate_inbounds function translate_indices(lattice::RegularLattice{D,T,PB,<:HomogeneousCell}, Is) where {D,T,PB}
#    @boundscheck map(x -> check_lattice_index(lattice, x), Is)
#    return Base.product(Is, Base.OneTo(length(lattice.unit_cell)))
#end
#Base.@propagate_inbounds function translate_indices(lattice::RegularLattice{D,T,PB,<:InhomogeneousCell}, Is, ig::Int) where {D,T,PB}
#    @boundscheck map(x -> check_lattice_index(lattice, x), Is)
#    return Base.product(Is, Base.OneTo(group_size(lattice.unit_cell, ig)), ig)
#end

#@inline to_default(lattice::RegularLattice{D,T,PB,<:TrivialCell}, x::Int, y::Int, z::Int) = CartesianIndex