text:
module GMT

using Printf, Dates, Statistics, Pkg

struct CTRLstruct
	limits::Vector{Float64}
	proj_linear::Vector{Bool}		# To know if images sent to GMT need Pad
	callable::Array{Symbol}			# Modules that can be called inside other modules
	pocket_call::Vector{Any}		# To temporarily store data needed by modules sub-calls
	pocket_B::Vector{String}		# To temporarily store opt_B grid and fill color to be reworked in psclip
	gmt_mem_bag::Vector{Ptr{Cvoid}}	# To temporarily store a GMT owned memory to be freed in gmt()
end

struct CTRLstruct2
	first::Vector{Bool}				# Signal that we are starting a new plot (used to set params)
	points::Vector{Bool}			# If maps are using points as coordinates
	fname::Vector{String}			# Store the full name of PS being constructed
end

# Function to change data of GMT.jl and hence force a rebuild in next Julia session
force_precompile() = Sys.iswindows() ? run(`cmd /c copy /b "$(pathof(GMT))" +,, "$(pathof(GMT))"`) : run(`touch '$(pathof(GMT))'`)

# Need to know what GMT version is available or if none at all to warn users on how to install GMT.
function get_GMTver()
	out = v"0.0"
	try						# First try to find an existing GMT installation (RECOMENDED WAY)
		(get(ENV, "FORCE_INSTALL_GMT", "") != "") && error("Forcing an automatic GMT install")
		ver = readlines(`gmt --version`)[1]
		out = ((ind = findfirst('_', ver)) === nothing) ? VersionNumber(ver) : VersionNumber(ver[1:ind-1])
		return out, false, "", "", "", ""
	catch err1;		println(err1)		# If not, install GMT
		ENV["INSTALL_GMT"] = "1"
		try
			depfile = joinpath(dirname(@__FILE__),"..","deps","deps.jl")	# File with shared lib names
			if isfile(depfile)
				include(depfile)		# This loads the shared libs names
				if (Sys.iswindows() && !isfile(_GMT_bindir * "\\gmt.exe"))		# If GMT was removed but depfile still exists
					Pkg.build("GMT");	include(depfile)
				end
			else
				Pkg.build("GMT");		include(depfile)
			end
			ver = readlines(`$(joinpath("$(_GMT_bindir)", "gmt")) --version`)[1]
			out = ((ind = findfirst('_', ver)) === nothing) ? VersionNumber(ver) : VersionNumber(ver[1:ind-1])
			return out, true, _libgmt, _libgdal, _libproj, _GMT_bindir
		catch err2;		println(err2)
			return out, false, "", "", "", ""
		end
	end
end
_GMTver, GMTbyConda, _libgmt, _libgdal, _libproj, _GMT_bindir = get_GMTver()

if (!GMTbyConda)		# In the other case (the non-existing ELSE branch) lib names already known at this point.
	_libgmt = haskey(ENV, "GMT_LIBRARY") ? ENV["GMT_LIBRARY"] : string(chop(read(`gmt --show-library`, String)))
	@static Sys.iswindows() ?
		(Sys.WORD_SIZE == 64 ? (_libgdal = "gdal_w64.dll") : (_libgdal = "gdal_w32.dll")) : (
			Sys.isapple() ? (_libgdal = string(split(readlines(pipeline(`otool -L $(_libgmt)`, `grep libgdal`))[1])[1])) : (
				Sys.isunix() ? (_libgdal = string(split(readlines(pipeline(`ldd $(_libgmt)`, `grep libgdal`))[1])[3])) :
				error("Don't know how to install this package in this OS.")
			)
		)
	@static Sys.iswindows() ?
		(Sys.WORD_SIZE == 64 ? (_libproj = "proj_w64.dll") : (_libproj = "proj_w32.dll")) : (
			Sys.isapple() ? (_libproj = string(split(readlines(pipeline(`otool -L $(_libgdal)`, `grep libproj`))[1])[1])) : (
				Sys.isunix() ? (_libproj = string(split(readlines(pipeline(`ldd $(_libgdal)`, `grep libproj`))[1])[3])) :
				error("Don't know how to use PROJ4 in this OS.")
			)
		)
end
const GMTver, libgmt, libgdal, libproj, GMT_bindir = _GMTver, _libgmt, _libgdal, _libproj, _GMT_bindir

const global G_API = [C_NULL]
const global img_mem_layout = [""]			# "TCP"	 For Images.jl. The default is "TRBa"
const global grd_mem_layout = [""]			# "BRP" is the default for GMT PS images.
const global current_view   = [""]			# To store the current viewpoint (-p)
const global multi_col   = Vector{Bool}(undef, 1);multi_col[1] = false	# To allow plottig multiple columns at once (init to false)
const global IamModern   = Vector{Bool}(undef, 1);IamModern[1] = false		# To know if we are in modern mode
const global FirstModern = Vector{Bool}(undef, 1);FirstModern[1] = false	# To know
const global IamSubplot  = Vector{Bool}(undef, 1);IamSubplot[1]  = false	# To know if we are in subplot mode
const global usedConfPar = Vector{Bool}(undef, 1);usedConfPar[1] = false	# Hacky solution for the session's memory trouble
const global ThemeIsOn   = Vector{Bool}(undef, 1);ThemeIsOn[1] = false		# To know if we have an active plot theme
const global convert_syntax = Vector{Bool}(undef, 1);convert_syntax[1] = false	# To only convert to hard core GMT syntax (like Vd=2)
const global show_kwargs = Vector{Bool}(undef, 1);show_kwargs[1] = false	# To just print the kwargs of a option call)
const global FMT = ["png"]
const global box_str = [""]
const def_fig_size  = "14c/9.5c"                   # Default fig size for plot like programs. Approx 16/11
const def_fig_axes_bak     = " -Baf -BWSen"        # Default fig axes for plot like programs
const def_fig_axes3_bak    = " -Baf -Bza"          # 		"" but for 3D views
const global def_fig_axes  = [def_fig_axes_bak]    # This one may be be changed by theme()
const global def_fig_axes3 = [def_fig_axes3_bak]   #		""
const global CTRL = CTRLstruct(zeros(6), [true], [:clip, :coast, :colorbar, :basemap, :logo, :text, :arrows, :lines, :scatter, :scatter3, :plot, :plot3, :hlines, :vlines], [nothing], ["",""], [C_NULL])
const global CTRLshapes = CTRLstruct2([true], [true], [""])
const prj4WGS84 = "+proj=longlat +datum=WGS84 +units=m +no_defs"	# This is used in many places
const CPTaliases = [:C :color :cmap :colormap :colorscale]
const global VMs = Union{Nothing, Vector{Symbol}, Matrix{Symbol}}
const global VMr = Union{Vector{<:Real}, Matrix{<:Real}}
# GItype = Union{GMTgrid, GMTimage} and GDtype = Union{GMTdataset, Vector{GMTdataset}} are edeclared in gmt_main
#const global unused_opts = [()]					# To track consumed options
#const global unused_subopts = [()]					# To track consumed options in sub-options

if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@optlevel"))
	@eval Base.Experimental.@optlevel 1
end

export
	GMTgrid, GMTimage, GMTdataset, GMTver, FMT, gmt, libgdal,
	arrows, arrows!, bar, bar!, bar3, bar3!, hlines, hlines!, lines, lines!, legend, legend!, vlines, vlines!,
	basemap, basemap!, blockmean, blockmedian, blockmode, clip, clip!, coast, coast!, colorbar, colorbar!,
	colorscale, colorscale!, contour, contour!, contourf, contourf!, events, filter1d, fitcircle, gmt2kml,
	gmtconnect, gmtconvert, gmtinfo, gmtmath, gmtregress, gmtread, gmtselect, gmtset, gmtsimplify, gmtspatial,
	gmtvector, gmtwrite, gmtwhich, grd2cpt, grd2kml, grd2xyz, grdblend, grdclip, grdcontour, grdcontour!, grdcut,
	grdedit, grdfft, grdfill, grdfilter, grdgradient, grdhisteq, grdimage, grdimage!, grdinfo, grdinterpolate,
	grdlandmask, grdmath, grdmask, grdpaste, grdproject, grdsample, grdtrack, grdtrend, grdvector, grdvector!,
	grdview, grdview!, grdvolume, greenspline, histogram, histogram!, image, image!, image_alpha!, image_cpt!,
	imshow, ind2rgb, isnodata, kml2gmt, logo, logo!, makecpt, mask, mask!, mapproject, movie, nearneighbor, plot, plot!,
	plot3, plot3!, plot3d, plot3d!, plotyy, project, pscontour, pscontour!, psconvert, psbasemap, psbasemap!,
	psclip, psclip!, pscoast, pscoast!, psevents, pshistogram, pshistogram!,
	psimage, psimage!, pslegend, pslegend!, psmask, psmask!, psrose, psrose!, psscale, psscale!, pssolar, pssolar!,
	psternary, psternary!, pstext, pstext!, pswiggle, pswiggle!, psxy, psxy!, psxyz, psxyz!, regress, resetGMT, rose,
	rose!, sample1d, scatter, scatter!, scatter3, scatter3!, solar, solar!, spectrum1d, sphdistance, sphinterpolate,
	sphtriangulate, surface, ternary, ternary!, text, text!, text_record, trend1d, trend2d, triangulate, gmtsplit,
	decorated, vector_attrib, wiggle, wiggle!, xyz2grd, gmtbegin, gmtend, gmthelp, subplot, gmtfig, inset, showfig,
	earthtide, gmtgravmag3d, pscoupe, pscoupe!, coupe, coupe!, psmeca, psmeca!, meca, meca!, psvelo, psvelo!, velo, velo!,

	mbimport, mbgetdata, mbsvplist, mblevitus,

	blendimg!, lonlat2xy, xy2lonlat, mat2ds, mat2grid, mat2img, slicecube, linspace, logspace, fields, tic, toc, theme,
	tern2cart, geodetic2enu, cpt4dcw, gd2gmt, gmt2gd, gdalread, gdalshade, gdalwrite, gadm,

	magic, rescale, stackgrids, delrows!,

	getband, getdriver, getlayer, getproj, getgeom, getgeotransform, toPROJ4, toWKT, importPROJ4,
	importWKT, importEPSG, gdalinfo, gdalwarp, gdaldem, gdaltranslate, gdalgrid, gdalvectortranslate, ogr2ogr,
	gdalrasterize, gdalbuildvrt, readraster, setgeotransform!, setproj!, destroy,
	delaunay, dither, buffer, centroid, intersection, intersects, polyunion, fromWKT,
	convexhull, difference, symdifference, distance, geomarea, pointalongline, polygonize, simplify,
	wkbUnknown, wkbPoint, wkbLineString, wkbPolygon, wkbMultiPoint, wkbMultiLineString, wkbMultiPolygon,
	wkbGeometryCollection, set_config_option,

	buffergeo, circgeo, epsg2proj, epsg2wkt, geod, invgeod, loxodrome, loxodrome_direct, loxodrome_inverse,
	orthodrome, proj2wkt, wkt2proj,

	colorzones!, rasterzones!, crop, doy2date, date2doy, yeardecimal, median, mean, std, nanmean, nanstd,

	append2fig, regiongeog, wmsinfo, wmstest, wmsread

include("common_docs.jl")
include("libgmt_h.jl")
(GMTver >= v"6") && include("libgmt.jl")
include("gmt_main.jl")
include("utils_types.jl")
include("grd_operations.jl")
include("common_options.jl")
include("gmtbegin.jl")
include("blendimg.jl")
include("blocks.jl")
include("contourf.jl")
include("filter1d.jl")
include("fitcircle.jl")
include("gadm.jl")
include("gmt2kml.jl")
include("gmtconnect.jl")
include("gmtconvert.jl")
include("gmtinfo.jl")
include("gmtlogo.jl")
include("gmtreadwrite.jl")
include("gmtset.jl")
include("gmtselect.jl")
include("gmtsimplify.jl")
include("gmtspatial.jl")
include("gmtregress.jl")
include("gmtvector.jl")
include("gmtwich.jl")
include("grd2cpt.jl")
include("grd2kml.jl")
include("grd2xyz.jl")
include("grdblend.jl")
include("grdclip.jl")
include("grdcut.jl")
include("grdedit.jl")
include("grdcontour.jl")
include("grdfft.jl")
include("grdfill.jl")
include("grdfilter.jl")
include("grdgmtmath.jl")
include("grdhisteq.jl")
include("grdimage.jl")
include("grdinfo.jl")
include("grdinterpolate.jl")
include("grdgradient.jl")
include("grdlandmask.jl")
include("grdmask.jl")
include("grdpaste.jl")
include("grdproject.jl")
include("grdsample.jl")
include("grdtrack.jl")
include("grdtrend.jl")
include("grdvector.jl")
include("grdview.jl")
include("grdvolume.jl")
include("greenspline.jl")
include("kml2gmt.jl")
include("loxodromics.jl")
include("makecpt.jl")
include("mapproject.jl")
include("movie.jl")
include("nearneighbor.jl")
include("plot.jl")
include("project.jl")
include("psbasemap.jl")
include("psclip.jl")
include("pscoast.jl")
include("pscontour.jl")
include("psconvert.jl")
include("pshistogram.jl")
include("psimage.jl")
include("pslegend.jl")
include("psmask.jl")
include("psscale.jl")
include("psrose.jl")
include("pssolar.jl")
include("pstext.jl")
include("psxy.jl")
include("pswiggle.jl")
include("rasterpolygonfuns.jl")
include("sample1d.jl")
include("spectrum1d.jl")
include("sphdistance.jl")
include("sphinterpolate.jl")
include("sphtriangulate.jl")
include("splitxyz.jl")
include("surface.jl")
include("subplot.jl")
include("themes.jl")
include("triangulate.jl")
include("trend1d.jl")
include("trend2d.jl")
include("xyz2grd.jl")
include("utils_project.jl")
include("choropleth_utils.jl")
include("webmapserver.jl")
include("seis/psmeca.jl")
include("geodesy/psvelo.jl")
include("geodesy/earthtide.jl")
include("MB/mbimport.jl")
include("MB/mbgetdata.jl")
include("MB/mbsvplist.jl")
include("MB/mblevitus.jl")
(GMTver > v"6.1.1") && include("potential/gmtgravmag3d.jl")
include("drawing.jl")

if (GMTver >= v"6")			# Needed to cheat the autoregister autobot
	include("get_enums.jl")
	include("gdal.jl")
	include("gdal_utils.jl")
	include("proj_utils.jl")
	using GMT.Gdal
end
include("imshow.jl")		# Include later because one method depends on knowing about GDAL

const global current_cpt = [GMTcpt()]		# To store the current palette
const global legend_type = [legend_bag()]	# To store Legends info

function __init__(test::Bool=false)
	if (GMTver == v"0.0" || test)
		println("\n\nYou don't seem to have GMT installed and the automatic installation also failed.\nYou will have to do it yourself.")
		t = Sys.iswindows() ? "Download and install the official version at (the '..._win64.exe':" *
		                      "\n\t\t https://github.com/GenericMappingTools/gmt/releases" : (
		                      Sys.isapple() ? "Install GMT with Homebrew: brew install gmt ghostscript ffmpeg" :
		                      "https://github.com/GenericMappingTools/gmt/blob/master/INSTALL.md#linux")
		println(t)
		return
	end

	if !isfile(libgmt) || ( !Sys.iswindows() && (!isfile(libgdal) || !isfile(libproj)) )
		println("\nDetected a previously working GMT.jl version but something has broken meanwhile.\n" *
		"(like updating your GMT instalation). Run this command in REPL and restart Julia.\n\n\t\tGMT.force_precompile()\n")
		return
	end

	clear_sessions(3600)		# Delete stray sessions dirs older than 1 hour
	G_API[1] = GMT_Create_Session("GMT", 2, GMT_SESSION_BITFLAGS)
	theme_modern()				# Set the MODERN theme
	haskey(ENV, "JULIA_GMT_IMGFORMAT") && (FMT[1] = ENV["JULIA_GMT_IMGFORMAT"])
	f = joinpath(readlines(`$(joinpath("$(GMT_bindir)", "gmt")) --show-userdir`)[1], "theme_jl.txt")
	(isfile(f)) && (theme(readline(f));	ThemeIsOn[1] = false)	# False because we don't want it reset in showfig()
	gmtlib_setparameter(G_API[1], "COLOR_NAN", "255")			# Stop those ugly grays
	#(GMTver == v"6.4") && gmtlib_setparameter(G_API[1], "MAP_EMBELLISHMENT", "auto")
end

include("precompile_GMT_i.jl")
_precompile_()

"""
GMT manipulating geographic and Cartesian data sets (including filtering, trend fitting, gridding, projecting, etc.)
and producing high quality illustrations.

Full modules list and docs in terse GMT style at $(GMTdoc)

Documentation for GMT.jl at https://www.generic-mapping-tools.org/GMT.jl/latest/
"""
GMT

end # module


text:
using VehicleModels, Plots, DataFrames
using PrettyPlots
# Name For Resutls
results_dir = string("testing");

dfs = Vector{DataFrame}(2) # create am empty DataFrame

# read in sample data
s_data = readtable("settings_data.csv")
obs_data = readtable("obstacle_data.csv")
dfs[1] = readtable("STATES.csv") # according to optimization
t_data = float(dfs[1][:t]);
SR_data = float(dfs[1][:SR]);
Jx_data = float(dfs[1][:Jx]);
t0=dfs[1][:t][1];tf=dfs[1][:t][end];
x0 = [dfs[1][:X][1],                 # 1. X, position
      dfs[1][:Y][1],                 # 2. Y, position
      dfs[1][:V][1],                 # 3. V, Lateral Speed
      dfs[1][:r][1],                 # 4. r, Yaw Rate
      dfs[1][:SA][1],                # 5. SA, Steering Angle
      dfs[1][:PSI][1],               # 6. PSI, Yaw angle
      dfs[1][:U][1],                 # 7. U, Longitudinal Speed
      dfs[1][:Ax][1]];               # 8. Ax, Longitudinal Acceleration

pa = Vpara(); # initialize parameter set
sol = Three_DOF(pa,x0,t_data,SR_data,Jx_data,t0,tf)

# extract data
t0=sol.t[1]; tf=sol.t[end]; pts = length(dfs[1][:V]);  #does not have to be same size, just doing this so the animation looks OK
x = [sol(t)[1] for t in linspace(t0,tf,pts)];
y = [sol(t)[2] for t in linspace(t0,tf,pts)];
v = [sol(t)[3] for t in linspace(t0,tf,pts)];
r = [sol(t)[4] for t in linspace(t0,tf,pts)];
sa = [sol(t)[5] for t in linspace(t0,tf,pts)];
psi = [sol(t)[6] for t in linspace(t0,tf,pts)];
u = [sol(t)[7] for t in linspace(t0,tf,pts)];
ax = [sol(t)[8] for t in linspace(t0,tf,pts)];

# build data for plots
t = convert(Vector, linspace(t0,tf,pts));
sp_SR = Linear_Spline(t_data, SR_data);
sp_Jx = Linear_Spline(t_data, Jx_data);
#sr  = sp_SR(t); # using Interpolations.jl
#jx  = sp_Jx(t); # using Interpolations.jl
sr  = sp_SR[t];
jx  = sp_Jx[t];

# put the data into a DataFrame
dfs[2]=x_to_DataFrame(t,x,y,sa,ax,psi,u,v,r,jx,sr) # according to ODE equations (RK4)

# plot the results
@unpack_Vpara pa  # unpack the parameters for the plots

all_plots(dfs,["opt.","RK4"],results_dir,obs_data,s_data,pa)

# simulate
panim_fun(dfs,["opt.","RK4"], results_dir, obs_data,s_data,pa)

ext: jl


text:
#writes one CSV file for how an ideal observer would update the v matrix
#another CSV file for the ground truth V

using MetaGen
using JSON
using Pipe: @pipe

include("helper_function.jl")
include("scripts/useful_functions.jl")

path = "../../scratch_work_07_16_21/09_18/shuffle_0/"

#could equally use input dictionary
#dict = @pipe "output.json" |> open |> read |> String |> JSON.parse
dict = @pipe (path * "output.json") |> open |> read |> String |> JSON.parse

num_videos = 50
num_frames = 20
threshold = 0.0
top_n = 5

params = Video_Params(n_possible_objects = 5)
receptive_fields = make_receptive_fields(params)
objects_observed, camera_trajectories = make_observations_office(dict, receptive_fields, num_videos, num_frames, threshold, top_n)

###############################################################################
#Set up the output file
outfile = string(path * "ideal_v_matrix.csv")
file = open(outfile, "w")

#return
print(file, "video_number&")
for i = 1:params.n_possible_objects
	print(file, "ideal_fa_", string(i), "&")
	print(file, "ideal_m_", string(i), "&")
end
print(file, "ground_truth_world_states")
print(file, "\n")

function print_helper(file, v::Int64, avg_v::Matrix{Float64}, gt_objects::Vector{Any})
	print(file, v, "&")
	n_rows,_ = size(avg_v)
    for i = 1:n_rows
        print(file, avg_v[i,1], "&")
        print(file, avg_v[i,2], "&")
    end
	print(file, gt_objects)
	print(file, "\n")
end

function print_ideal_v(file, dict::Array{Any}, params::Video_Params,
	camera_trajectories::Matrix{Camera_Params},
	objects_observed::Matrix{Array{Detection2D}})

    num_videos, num_frames = size(camera_trajectories)

	real_detections = Matrix{Array{Detection2D}}(undef, num_videos, num_frames)

    #could make these zeros
    alphas = fill(1, (params.n_possible_objects,2))
    betas = fill(1, (params.n_possible_objects,2))
    avg_v = fill(0.0, (params.n_possible_objects,2))

	gt_objects = get_ground_truth(dict, num_videos) #3D objects
    for v = 1:num_videos
		#println("v ", v)
		#println("gt_objects ", gt_objects[v])
		old_beta = deepcopy(betas)
		old_alpha = deepcopy(alphas)
        for f = 1:num_frames
            paramses = fill(params, length(gt_objects[v]))
            camera_paramses = fill(camera_trajectories[v, f], length(gt_objects[v]))
            gt_objects_2D = map(render, paramses, camera_paramses, gt_objects[v])
            gt_objects_2D = Array{Detection2D}(gt_objects_2D)

			#real_detections = filter!(within_frame, gt_objects_2D)

			#observations

            #inferences = convert(Vector{Int64}, dict[v]["views"][f]["inferences"]["labels"])
            #println("inferences ", inferences)
            #alphas, betas = update_alpha_beta(alphas, betas, inferences, gt_objects_2D)
			alphas, betas = update_alpha_beta(false, alphas, betas, objects_observed[v,f], gt_objects_2D)
        end
		change_alpha = alphas - old_alpha
		change_beta = betas - old_beta
		miss_rate = (change_alpha[:,2]) ./ ((change_alpha[:,2]) .+ (change_beta[:,2])) #don't need -1 since only looking at changes
		println(v)
		println(miss_rate)

        #println("alphas fa ", alphas[:,1])
        #println("betas fa ", betas[:,1])
		#println("alphas miss ", alphas[:,2])
        #println("betas miss ", betas[:,2])
        avg_v[:,1] = alphas[:,1] ./ (betas[:,1].^2)#fa rate. mean of gamma distribution
        avg_v[:,2] = alphas[:,2] ./ (alphas[:,2] .+ betas[:,2])#miss rates. mean of beta distribution
        print_helper(file, v, avg_v, gt_objects[v])
    end
    return alphas, betas
end

alphas, betas = print_ideal_v(file, dict, params, camera_trajectories, objects_observed)
close(file)

################################################################################
#ground truth
file = open(string(path * "ground_truth_V.csv"), "w")

print(file, "video_number&")
for i = 1:params.n_possible_objects
	print(file, "gt_fa_", string(i), "&")
	print(file, "gt_m_", string(i), "&")
end
print(file, "\n")

#for ground truth, would just subtract 1 from the alphas and betas matrices
gt_v = fill(0.0, (params.n_possible_objects,2))
gt_v[:,1] = (alphas[:,1] .- 1) ./ ((betas[:,1] .- 1).^2)
gt_v[:,2] = (alphas[:,2] .- 1) ./ ((alphas[:,2] .- 1) .+ (betas[:,2] .- 1))

gt_objects = get_ground_truth(dict, num_videos)
for v=1:num_videos
	print_helper(file, v, gt_v, gt_objects[v])
end

close(file)

################################################################################
#add ground truth per frame to dictionary
out = write_gt_to_dict(dict, camera_trajectories, gt_objects)

open(path * "output_with_gt.json","w") do f
	JSON.print(f,out)
end

max_stars_repo_name: zhangir-azerbayev/MetaGen

max_issues_repo_name: zhangir-azerbayev/MetaGen

max_forks_repo_name: zhangir-azerbayev/MetaGen


text:
module Router # Bukdu

using ..Bukdu.Naming
using ..Bukdu.Deps
using ..Bukdu: handle_request

"""
    Router.call

To get the object of the action without going through HTTP server.
Results are named tuple `(got= , resp= , route= )`.
"""
function call
end

"""
    Router.call(verb, path::String, headers=[], body=UInt8[])::NamedTuple{(:got, :resp, :route)}
"""
function call(verb, path::String, headers=[], body=UInt8[])::NamedTuple{(:got, :resp, :route)}
    method = Naming.verb_name(verb)
    req = Deps.Request(method, path, headers, body)
    call(req)
end

"""
    Router.call(req::Deps.Request)::NamedTuple{(:got, :resp, :route)}
"""
function call(req::Deps.Request)::NamedTuple{(:got, :resp, :route)}
    handle_request(req)
end

module Helpers # Bukdu.Router

using ...Bukdu.Routing
using ...Bukdu: ApplicationController

struct URLPathError <: Exception
    msg
end

function url_path(verb, C::Type{<:ApplicationController}, action)
    path = Routing.route_path(verb, C, action)
    if path isa Nothing
        URLPathError(string("failed to get the path: ", join([verb, C, action], ' ')))
    else
        path
    end
end

end # module Bukdu.Router.Helpers

end # module Bukdu.Router


text:
using JLD2, Random, LinearAlgebra, Statistics, CSV, DataFrames, FreqTables, Distributions,
FreqTables, Optim, HTTP, GLM, ForwardDiff

function PS4()
    # includes both the mlogit and lgwt functions
    include("functions.jl")

    #:::::::::::::::::::
    # question 1   :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    #:::::::::::::::::::

    url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2021/master/ProblemSets/PS4-mixture/nlsw88t.csv"
    df = CSV.read(HTTP.get(url).body, DataFrame)
    X = [df.age df.white df.collgrad]
    Z = hcat(df.elnwage1, df.elnwage2, df.elnwage3, df.elnwage4,
    df.elnwage5, df.elnwage6, df.elnwage7, df.elnwage8)
    y = df.occ_code

    startvals = [0.04037443688915296, 0.2439931891755024, -1.571319982008478,
    0.04332546800040571, 0.1468547553146026, -2.95910294680261, 0.1020574214746393,
     0.7473078162016722, -4.120049525011452, 0.037562865451966576, 0.6884890953253581,
     -3.6557702948516075, 0.0204543204452073, -0.3584022910222397, -4.376928748361732,
      0.1074636880063156, -0.5263752170984033, -6.199198099540245, 0.11688249286984193,
       -0.28705670558923196, -5.322248251070095, 1.30747948066456]

    # mlogit_with_Z included in the functions.jl script
    td = TwiceDifferentiable(theta -> mlogit_with_Z(theta, X, Z, y), startvals; autodiff = :forward)

           # run the optimizer
           theta_hat_optim_ad = optimize(td, startvals, LBFGS(), Optim.Options(g_tol = 1e-5, iterations=100_000, show_trace=true, show_every=1))

           theta_hat_mle_ad = theta_hat_optim_ad.minimizer
           # evaluate the Hessian at the estimates
           H  = Optim.hessian!(td, theta_hat_mle_ad)
           theta_hat_mle_ad_se = sqrt.(diag(inv(H)))
              println(
              DataFrame(Row =  1:length(theta_hat_mle_ad), Coef = theta_hat_mle_ad, se = theta_hat_mle_ad_se)
              )

    #:::::::::::::::::::
    # question 2   :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    #:::::::::::::::::::
# they make more sense here as they are positive

    #:::::::::::::::::::
    # question 3   :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    #:::::::::::::::::::
    # define distribution
    d = Normal(0,1) # mean=0, standard deviation=1

    # get quadrature nodes and weights for 7 grid points
    nodes, weights = lgwt(7,-4,4)

    # now compute the integral over the density and verify it's 1
    sum(weights.*pdf.(d,nodes))

    # now compute the expectation and verify it's 0
    sum(weights.*nodes.*pdf.(d,nodes))

    #:::::::::::::::::::
    # question 4   :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    #:::::::::::::::::::

    # b)
    #::::
    # define distribution
    d = Normal(0,2) # mean=0, standard deviation=1
    nstd = 5*sqrt(2)
    # get quadrature nodes and weights for 7 grid points
    nodes, weights = lgwt(7,-nstd,nstd)

    # now compute the expectation of x^2 where x~N(0,2)
    println(
    sum(weights.*nodes.*nodes.*pdf.(d,nodes))
    )
    # get quadrature nodes and weights for 10 grid points
    nodes, weights = lgwt(10,-nstd,5*nstd)
    println(
    sum(weights.*nodes.*nodes.*pdf.(d,nodes))
    )

    # The quadrature nodes and weights for 10 grid points yields better approximation as it is very close to 4, the variance

    # c)
    #::::
    NumObs = 10^6
    x = rand(Normal(0,2),NumObs)
    x_sq = x.*x

    # pdf of N(0,2)
          function f(y)
           1/(2*sqrt(2*pi)) * exp(-0.5 * (y/2)^2)
           end
    w = (2*nstd)/NumObs
    fuct1 = map(y -> f(y), x) .* x_sq
    MonteCarloEst1 = w * sum(fuct1) # result very close to 4
    println(MonteCarloEst1)

    fuct2 = map(y -> f(y), x) .* x
    MonteCarloEst2 = w * sum(fuct2) # result very close to 0
    println(MonteCarloEst2)

    fuct3 = map(y -> f(y), x)
    MonteCarloEst3 = w * sum(fuct3) #result Not too close to one;close if Numobs = 10^8
    println(MonteCarloEst3)
    # we can also just use pdf() to get the same results
    # pdf(d::MultivariateDistribution, x::AbstractArray)

    println(sum(
     w .*
    (x_sq .* pdf.(Normal(0,2),x))
    ))

    println(sum(
    w .*
    (x .* pdf.(Normal(0,2),x))
    ))

    println(sum(
    w .*
    pdf.(Normal(0,2),x)
    ))

    # D = 10^6 yields better approximation

    # D)
    #:::::
    # both methods are numerical methods for integral approximation; one can be mapped to the other

    #:::::::::::::::::::
    # question 5   :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    #:::::::::::::::::::

    # using Gauss-Legendre quadrature
    function mixedlogit(theta, gridPoints, X, Z, y)

            alpha = theta[1:end-1]
            gamma = theta[end]
            K = size(X,2)
            J = length(unique(y))
            N = length(y)
            bigY = zeros(N,J)
            for j=1:J
                bigY[:,j] = y.==j
            end
            bigAlpha = [reshape(alpha,K,J-1) zeros(K)]

            T = promote_type(eltype(X),eltype(theta))
            num   = zeros(T,N,J)
            dem   = zeros(T,N)

     # define distribution
     d = Normal(0,1) # mean=0, standard deviation=1
     # get quadrature nodes and weights for 7 grid points
     nodes, weights = lgwt(gridPoints,-4,4)

            for j=1:J
                num[:,j] = sum(
                (exp.(X*bigAlpha[:,j] .+ (Z[:,j] .- Z[:,J])*gamma)).*
                weights .* pdf.(d,nodes)
                )
                dem .+= num[:,j]
            end

            P = num./repeat(dem,1,J)

            loglike = -sum( bigY.*log.(P) )

            return loglike
        end

        # using Monte Carlo
        function mixedlogit(theta, NumObs, X, Z, y)

                alpha = theta[1:end-1]
                gamma = theta[end]
                K = size(X,2)
                J = length(unique(y))
                N = length(y)
                bigY = zeros(N,J)
                for j=1:J
                    bigY[:,j] = y.==j
                end
                bigAlpha = [reshape(alpha,K,J-1) zeros(K)]

                T = promote_type(eltype(X),eltype(theta))
                num   = zeros(T,N,J)
                dem   = zeros(T,N)

         # define distribution
         x = rand(Normal(0,1),NumObs) # mean=0, standard deviation=1
         #nstd = 4
         w = (8)/NumObs

                for j=1:J
                    num[:,j] = sum(
                    (exp.(X*bigAlpha[:,j] .+ (Z[:,j] .- Z[:,J])*gamma)).*
                    w .* pdf.(Normal(0,1),x)
                    )
                    dem .+= num[:,j]
                end

                P = num./repeat(dem,1,J)

                loglike = -sum( bigY.*log.(P) )

                return loglike
            end
end

#:::::::::::::::::::
# question 6   :::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
#:::::::::::::::::::
PS4()


text:
Dunit = 1.0
Eunit = 1.0
Tunit = 1.0
Funit = 1.0
Cunit = 1.0

kb = 8.6173303e-5*Eunit/Tunit
ϵ₀ = 8.85418782e-14*Funit/Dunit
q = 1.60217662e-19*Cunit

function defineConstants()
    global kb = 8.6173303e-5*Eunit/Tunit
    global ϵ₀ = 8.85418782e-14*Funit/Dunit
    global q = 1.60217662e-19*Cunit
end

function SemiUseUnits()
    global Dunit = u"cm"
    global Eunit = u"eV"
    global Tunit = u"K"
    global Funit = u"F"
    global Cunit = u"C"
    defineConstants()
end

function SemiUstrip()
    global Dunit = 1.0
    global Eunit = 1.0
    global Tunit = 1.0
    global Funit = 1.0
    global Cunit = 1.0
    defineConstants()
end


text:
using Documenter
using DataDeps

using Documenter

makedocs(
    format=Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    sitename = "DataDeps.jl",
    modules = [DataDeps],
    repo="https://github.com/oxinabox/DataDeps.jl/blob/{commit}{path}#L{line}",
    authors="Lyndon White",
)

deploydocs(
    repo = "github.com/oxinabox/DataDeps.jl.git",
)


text:
#=
Created on 13 August 2019
Last update: 14 July 2020

@author: Michiel Stock
michielfmstock@gmail.com

Functions to simulate species interaction networks with optimal transport.
=#

"""
    zerodiv(a, b)

Division which results 0 if both arguments are 0.
"""
zerodiv(a, b) = a == b == zero(a) ? zero(a) : a / b

function optimaltransport_(M::AbstractMatrix, a::AbstractVector; λ::Real=1.0)
    @assert length(a) == size(M, 1) throw(DimensionMismatch(
                "a ($(size(a))) does not match M ($(size(M)))"))
    K = exp.(λ * M)
    return Diagonal(a ./ sum(K, dims=2)[:]) * K
end

function optimaltransport_inplace(M::AbstractMatrix, a::AbstractVector, b::AbstractVector;
                        λ::Real=1.0, maxitter::Int=100, ϵ::Real=1e-5)
    Q = exp.(λ * M)
    iter = 0
    while iter < maxitter && maximum(abs.(sum(Q, dims=1)[:] - b)) > ϵ
        iter += 1
        lmul!(Diagonal(a ./ sum(Q, dims=2)[:]), Q)
        rmul!(Q, Diagonal(b ./ sum(Q, dims=1)[:]))
    end
    return Q
end

# version that is differentiable
function optimaltransport_(M::AbstractMatrix, a::AbstractVector, b::AbstractVector;
                        λ::Real=1.0, maxitter::Int=100, ϵ::Real=1e-5)
    Q = exp.(λ * M)
    iter = 0
    while iter < maxitter && maximum(abs.(sum(Q, dims=2)[:] - a)) > ϵ
        iter += 1
        Q = Diagonal(zerodiv.(a, sum(Q, dims=2)[:])) * Q
        Q = Q * Diagonal(zerodiv.(b, sum(Q, dims=1)))
    end
    return Q
end

function h(A, x, n=10)
    res = exp.(x)
    for i in 1:n
        res = A * x
    end
    return sum(res.^2)
end

Ps = [rand(10, 20) |> x -> x/sum(x) for i in 1:11]
M = randn(10, 20)
A = hcat([sum(P, dims=2) for P in Ps]...)
B = hcat([sum(P, dims=1)' for P in Ps]...)

"""Parallel version of optimal transport"""
function optimaltransport(M::AbstractMatrix, A::AbstractMatrix, B::AbstractMatrix;
                        λ::Real=1.0, maxitter::Int=10, ϵ::Real=1e-5)
    K = exp.(λ * M)
    U, V = A, B
    for itter in 1:maxitter
        U = A ./ (K * V)
        V = B ./ (K' * U)
    end
    return U, V
end

function par_kl(Ps, M, U, V)
    #Q = similar(first(Ps))
    K = exp.(M)
    l = 0.0
    for (i, P) in enumerate(Ps)
        Q = Diagonal(U[:,i]) * K * Diagonal(V[:,i])
        l += relative_entropy(P, Q)
    end
    return l
end

function optimaltransport(M::AbstractMatrix,
                a::Union{AbstractVector,Nothing}=nothing,
                b::Union{AbstractVector,Nothing}=nothing;
                λ::Real=1.0, maxitter::Int=1000, ϵ::Real=1e-5)
    if a == nothing && b == nothing
        return exp.(λ * M) |> x -> x / sum(x)
    elseif b == nothing
        return optimaltransport_(M, a, λ=λ)
    elseif a == nothing
        return optimaltransport_(M', b, λ=λ)'
    else
        return optimaltransport_(M, a, b, λ=λ, maxitter=maxitter, ϵ=ϵ)
    end
end

optimaltransport(a::AbstractVector, b::AbstractVector; λ::Real=0.0) = a * b'

max_stars_repo_name: PoisotLab/OTSIN

max_issues_repo_name: PoisotLab/OTSIN

max_forks_repo_name: PoisotLab/OTSIN


text:
# Convergence experiment (convergence_rates-pseudosamplers-gmm)
#
# Compares the Stein discrepancy convergence behavior of herding
# with that of independent random sampling. Make sure to run
# the matlab script at
#
# src/experiments/skh/MoG_experiments/script_test.m
#
# with all the desired parameters used here.

include("experiment_utils.jl")

# this controls the seed used for the approximating the wasserstein metric
@parseintcli seed "s" "seed" 7
# the dimension of the space
@parseintcli d "d" "dimension" 1
# the sampler ["FW", "FCFW", "IID", "QMC"]
@parsestringcli sampler "q" "sampler" "IID"
# the gap between the two modes
@parseintcli gap "g" "gap" 5
# approx wasserstein sample sizes
@parseintcli wasserstein_n "w" "wassersteinn" 5000

using MAT

using SteinDistributions: SteinGMM, SteinDiscrete, cdf, rand, gradlogdensity
using SteinDiscrepancy: gsd, wasserstein1d, approxwasserstein

# Select an optimization problem solver
#solver = "clp"
solver = "gurobi"
# the number of components
K = 2
# distribution sampled from
distname = "$K-gmm"
# get the results directory of the sample files
resultsdir = joinpath("src", "experiments", "skh", "MoG_experiments", "results")

# open dat matlab file
samplefile = @sprintf("matlab_skh_gmm_method=%s_d=%d_comps=%d_gap=%d.mat",
                       sampler, d, K, gap)
samplefilepath = joinpath(resultsdir, samplefile)
samples = matread(samplefilepath)
# get some params from the data
X = samples["X"]
n = size(X,1)
pointweights = samples["weights"]
gmmweights = vec(samples["pi_prob"])
gmmsigmas = samples["Sigma_mix"]
gmmmus = samples["mu_mix"]
target = SteinGMM(gmmmus, gmmsigmas, gmmweights)
# define gradlogp
function gradlogp(x::Array{Float64,1})
    gradlogdensity(target, x)
end
if sampler == "IID"
    # we use our own iid sampler so we can run multiple seeds
    X = rand(target, n)
end
# Sample sizes at which optimization problem will be solved
ns = vcat(10:10:min(1000,n))

@printf("Beginning optimization for %s, dimension=%d\n", sampler, d)
for i in ns
    sampleweights = vec(pointweights[i+1,1:i])
    # Compute Stein discrepancy for first i points in this trial
    res = gsd(points=X[1:i,:],
              weights=sampleweights,
              gradlogdensity=gradlogp,
              solver=solver)
    println("\tn = $(i), objective = $(res.objectivevalue)")
    wasserstein = Inf
    wasserstein_lb = Inf
    wasserstein_ub = Inf
    if d == 1
        (wasserstein, error) = wasserstein1d(X[1:i,:],
                                             weights=sampleweights,
                                             targetcdf=x::Float64 -> cdf(target, x))
        wasserstein_lb = wasserstein - error
        wasserstein_ub = wasserstein + error
    else
        (wasserstein_lb, wasserstein_ub) =
            @setseed seed approxwasserstein(points=X[1:i,:],
                                            weights=sampleweights,
                                            targetsamplegen=() -> rand(target, wasserstein_n),
                                            solver=solver)
        wasserstein = (wasserstein_lb + wasserstein_ub) / 2
    end
    # Package and save results
    instance_data = Dict{Any, Any}(
        "seed" => seed,
        "mus" => gmmmus,
        "sigmas" => gmmsigmas,
        "distweights" => gmmweights,
        "sampler" => sampler,
        "gap" => gap,
        "distname" => distname,
        "n" => i,
        "d" => d,
        "X" => res.points,
        "q" => res.weights,
        "g" => res.g,
        "gprime" => res.gradg,
        "objectivevalue" => res.objectivevalue,
        "operatorg" => res.operatorg,
        "wasserstein" => wasserstein,
        "wasserstein_lb" => wasserstein_lb,
        "wasserstein_ub" => wasserstein_ub,
        "edgetime" => res.edgetime,
        "solvetime" => res.solvetime,
        "wasserstein_n" => wasserstein_n,
    )

    save_json(
        instance_data;
        dir="convergence_rates-pseudosamplers-gmm",
        distname=distname,
        sampler=sampler,
        n=i,
        d=d,
        gap=gap,
        wassersteinn=wasserstein_n,
        seed=seed,
    )
end

@printf("COMPLETE!\n")


text:
function PID_factory(;name, kp=1.0, ki=0.0, kd=0.0)
    @parameters _kp _ki _kd t
    @variables e(t) Ie(t) de(t) u(t)
    D = Differential(t)

    eqs = [
        D(Ie) ~ e
        D(e) ~ de
        u ~ _ki*Ie+_kp*e+_kd*de
    ]
    ODESystem(eqs, t ;name, defaults=[u=>0.0,
                                      e=>0.0, de=>0.0, Ie=>0.0,
                                      _kp=>kp, _ki=>ki, _kd=>kd])
end

function lead_lag_factory(;name,k=1.0, zero=1.0, pole=0.0)
    @parameters t _zero _pole _k
    @variables input(t) dinput(t) output(t) doutput(t)
    D = Differential(t)
    eqs = [
        D(input) ~ dinput
        D(output) ~ doutput
        _pole*doutput + output ~ _k*(_zero*dinput + input)
    ]
    ODESystem(eqs, t ;name, defaults=[output=>0.0, doutput=>0.0,
                                      input=>0.0, dinput=>0.0,
                                      _zero=>zero, _pole=>pole, _k=>k])
end

function sensor(;name, bw=1)
    @parameters t a2 a1
    @variables input(t) output(t) doutput(t) ddoutput(t)
    D = Differential(t)
    eqs = [
        D(output) ~ doutput
        D(doutput) ~ ddoutput
        input ~ output+a1*doutput+a2*ddoutput
    ]
    ODESystem(eqs, t ;name, defaults=[output=>0.0, doutput=>0.0,
                                      ddoutput=>0.0, input=>0.0,
                                      a1=>2*0.7/(2π*bw), a2=>1/(2π*bw^2)])
end

function notch(;name,ω=1.0,gain=1.0,width=1.0)
    @parameters t i1 i2 o1 o2
    @variables input(t) dinput(t) ddinput(t) output(t) doutput(t) ddoutput(t)
    D = Differential(t)
    eqs = [
        D(output) ~ doutput
        D(doutput) ~ ddoutput
        D(input) ~ dinput
        D(dinput) ~ ddinput
        input+i1*dinput+i2*ddinput ~ output+o1*doutput+o2*ddoutput
    ]
    ODESystem(eqs, t ;name, defaults=[output=>0.0, doutput=>0.0, ddoutput=>0.0,
                                      input=>0.0, dinput=>0.0, ddinput=>0.0,
                                      o1=>width/(2π*ω), o2=>1/(2π*ω)^2,
                                      i1=>width*gain/(2π*ω), i2=>1/(2π*ω)^2])

end

function amp_factory(;name, a=1.0)
    @parameters t, _a
    @variables input(t) output(t)
    eqs = [
        output ~ _a*input
    ]
    ODESystem(eqs, t;name, defaults=[output=>0.0,
                                     input=>0.0,
                                     _a=>a])
end

export PID_factory, lead_lag_factory, amp_factory, sensor, notch


text:
using ElectronProgressBars
using Juno
using Logging

function demo(delay; title="delay=$delay", n = 3, m = 50)
    Juno.progress(name="$title (outer)") do id
        for i in 1:n
            Juno.progress(name="$title (inner)") do id
                for j in 1:m
                    @debug "$title (inner)" progress=j/m _id=id
                    sleep(delay)
                end
            end
            @debug "$title (outer)" progress=i/n _id=id
        end
    end
end

with_logger(ElectronProgressBars.get_logger()) do
    @sync begin
        @async demo(0.01)
        @async demo(0.02)
    end
end

with_logger(ElectronProgressBars.get_logger()) do
    @sync begin
        Threads.@spawn demo(0.01)
        Threads.@spawn demo(0.02)
    end
end


text:
mutable struct KalmanFiltered{T}
    filtered::Array{T}
    predicted::Array{T}
    error_cov::Array{T}
    pred_error_cov::Array{T}
    model::StateSpaceModel
    y::Array{T}
    u::Array{T}
    loglik::T
end

function Base.show(io::IO, filt::KalmanFiltered{T}) where T
    n = size(filt.y, 1)
    dx, dy = filt.model.nx, filt.model.ny
    println("KalmanFiltered{$T}")
    println("$n observations, $dx-D process x $dy-D observations")
    println("Negative log-likelihood: $(filt.loglik)")
end

function kalman_filter(y::Array{T}, model::StateSpaceModel{T}; u::Array{T}=zeros(size(y,1), model.nu)) where T

    @assert size(u,1) == size(y,1)
    @assert size(y,2) == model.ny
    @assert size(u,2) == model.nu

    function kalman_recursions(y_i::Vector{T}, u_i::Vector{T},
                                  C_i::Matrix{T}, D_i::Matrix{T}, W_i::Matrix{T},
                                  x_pred_i::Vector{T}, P_pred_i::Matrix{T})
        if !any(isnan, y_i)
            innov =  y_i - C_i * x_pred_i - D_i * u_i
            S = C_i * P_pred_i * C_i' + W_i  # Innovation covariance
            K = P_pred_i * C_i' / S # Kalman gain
            x_filt_i = x_pred_i + K * innov
            P_filt_i = (I - K * C_i) * P_pred_i
            dll = (dot(innov,S\innov) + logdet(S))/2
        else
            x_filt_i = x_pred_i
            P_filt_i = P_pred_i
            dll = 0
        end
        return x_filt_i, P_filt_i, dll
    end #kalman_recursions

    y = y'
    u = u'
    n = size(y, 2)
    x_pred = zeros(model.nx, n)
    x_filt = zeros(x_pred)
    P_pred = zeros(model.nx, model.nx, n)
    P_filt = zeros(P_pred)
    log_likelihood = n*model.ny*log(2pi)/2

    # first iteration
    A_1 = model.A(1)
    x_pred[:, 1] = model.x1
    P_pred[:, :, 1] = model.P1
    x_filt[:, 1], P_filt[:,:,1], dll = kalman_recursions(y[:, 1], u[:, 1],
                                            model.C(1), model.D(1), model.W(1),
                                            model.x1, model.P1)
    log_likelihood += dll

    for i=2:n
        A_i1 = model.A(i-1)
        x_pred[:, i] =  A_i1 * x_filt[:, i-1] + model.B(i-1) * u[:, i-1]
        P_pred[:, :, i] = A_i1 * P_filt[:, :, i-1] * A_i1' + model.V(i-1)
        x_filt[:, i], P_filt[:,:,i], dll = kalman_recursions(y[:, i], u[:, i],
                                                model.C(i), model.D(i), model.W(i),
                                                x_pred[:,i], P_pred[:,:,i])
        log_likelihood += dll
    end

    return KalmanFiltered(x_filt', x_pred', P_filt, P_pred, model, y', u', log_likelihood)
end

function loglikelihood(y::Array{T}, model::StateSpaceModel{T}; u::Array{T}=zeros(size(y,1), model.nu)) where T

    y, u = y', u'
    n = size(y, 2)

    @assert !any(isnan, u)
    @assert size(y, 1) == model.ny
    @assert size(u, 1) == model.nu
    @assert size(u, 2) == n

    y_notnan = .!isnan.(y)
    y = y .* y_notnan

    I0ny = zeros(model.ny, model.ny)
    Ksparse = issparse(model.A(1)) && issparse(model.V(1)) && issparse(model.C(1))

    x_pred_t, P_pred_t  = model.x1, model.P1
    ut, Ct, Wt          = zeros(model.nu), model.C(1), model.W(1)
    innov_t, innov_cov_inv_t  = zeros(model.ny), I0ny

    log_likelihood = n*model.ny*log(2pi)/2
    marginal_likelihood(innov::Vector, innov_cov::Matrix,
          innov_cov_inv::Matrix) =
              (dot(innov, innov_cov_inv * innov) + logdet(innov_cov))/2

    for t = 1:n

        # Predict using last iteration's values
        if t > 1

            At, But, Vt = model.A(t-1), model.B(t-1) * ut, model.V(t-1)
            AP_pred_t   = At * P_pred_t
            Kt          = AP_pred_t * Ct' * innov_cov_inv_t
            Ksparse && (Kt = sparse(Kt))

            x_pred_t = At * x_pred_t + But + Kt * innov_t
            P_pred_t = any(Wt .!= 0) ? AP_pred_t * (At - Kt * Ct)' + Vt : Vt
            P_pred_t = (P_pred_t + P_pred_t')/2

        end #if

        # Assign new values
        ut  = u[:, t]
        Ct  = model.C(t)
        Dut = model.D(t) * ut
        Wt  = model.W(t)

        # Check for and handle missing observation values
        missing_obs = !all(y_notnan[:, t])
        if missing_obs
            ynnt = y_notnan[:, t]
            I1, I2 = spdiagm(ynnt), spdiagm(!ynnt)
            Ct, Dut = I1 * Ct, I1 * Dut
            Wt = I1 * Wt * I1 + I2 * Wt * I2
        end #if

        # Compute nessecary filtering quantities
        innov_t           = y[:, t] - Ct * x_pred_t - Dut
        innov_cov_t       = Ct * P_pred_t * Ct' + Wt |> full
        nonzero_innov_cov = all(diag(innov_cov_t) .!= 0)
        innov_cov_inv_t   = nonzero_innov_cov ? inv(innov_cov_t) : I0ny

        nonzero_innov_cov && (log_likelihood += marginal_likelihood(innov_t, innov_cov_t, innov_cov_inv_t))

    end #for

    return log_likelihood

end #loglikelihood


text:
abstract type SpinHalf <: Model end

function getTotalSx(model::SpinHalf)::MPO
    ampo = OpSum()

    for j = 1:length(model.sites)
        ampo += 2.0, "Sx", j
    end

    return MPO(ampo, model.sites)
end

function getTotalSy(model::SpinHalf)::MPO
    ampo = OpSum()

    for j = 1:length(model.sites)
        ampo += 2.0, "Sy", j
    end

    return MPO(ampo, model.sites)
end

function getTotalSz(model::SpinHalf)::MPO
    ampo = OpSum()

    for j = 1:length(model.sites)
        ampo += 2.0, "Sz", j
    end

    return MPO(ampo, model.sites)
end


text:
struct SingleOrder{T <: AbstractFloat}
    hour::Int64                    # Hour for single order
    independent_volume::T          # Price independent order volume
    dependent_volumes::Vector{T}   # Price dependent order volume
    prices::Vector{T}              # Price dependent order prices

    function SingleOrder(h::Integer, x_i::AbstractFloat, x_d::AbstractVector, prices::AbstractVector)
        T = promote_type(typeof(x_i), eltype(x_d), eltype(prices), Float32)
        return new{T}(h, x_i, convert(Vector{T},x_d), convert(Vector{T}, prices))
    end
end
independent(order::SingleOrder) = order.independent_volume
dependent(order::SingleOrder) = order.dependent_volumes[2:end-1]
lowertechnical(order::SingleOrder) = order.dependent_volumes[1]
uppertechnical(order::SingleOrder) = order.dependent_volumes[end]

struct BlockOrder{T <: AbstractFloat}
    interval::Tuple{Int,Int}   # Time range for block order
    volume::T                  # Asked volume for block order
    price::T                   # Asked price for block order

    function BlockOrder(interval::Tuple{Int,Int}, volume::AbstractFloat, price::AbstractFloat)
        T = promote_type(typeof(volume), typeof(price), Float32)
        return new{T}(interval, volume, price)
    end
end

struct OrderStrategy{T <: AbstractFloat}
    horizon::Horizon                        # Planning horizon
    bidlevels::Vector{Vector{T}}            # Possible prices levels

    single_orders::Vector{SingleOrder{T}}   # Order orders each hour
    block_orders::Vector{BlockOrder{T}}     # All block orders

    function OrderStrategy(horizon::Horizon,
                           bidlevels::Vector{Vector{T}},
                           single_orders::Vector{SingleOrder{T}},
                           block_orders::Vector{BlockOrder{T}}) where T <: AbstractFloat
        return new{T}(horizon, bidlevels, single_orders, block_orders)
    end
end

function blockbidprice(bidlevels, hours_per_block, i)
    return mean([bidlevels[i][2:end-1] for i in hours_per_block])[i]
end

function OrderStrategy(model::AbstractHydroModel; variant = :rp)
    sp = model.internalmodel
    horizon = HydroModels.horizon(model)
    regulations = model.data.regulations
    bidlevels = model.data.bidlevels

    xᴵ = optimal_decision(sp, :xᴵ)
    xᴰ = optimal_decision(sp, :xᴰ)
    xᴮ = optimal_decision(sp, :xᴮ)

    # Check data length consistency
    hs = 1:nhours(horizon)
    hours_per_block = [collect(h:ending) for h in hs for ending in hs[h+regulations.blockminlength-1:end]]
    @assert nhours(horizon) == length(xᴵ) "Incorrect horizon of price independent single orders"
    @assert nhours(horizon) == JuMP.size(xᴰ)[2] "Incorrect horizon of price dependent single orders"
    # @assert length(hours_per_block) == JuMP.size(xᴮ)[2] "Incorrect horizon of block orders"
    # @assert all([length(bidlevels[t]) .== JuMP.size(xᴰ)[t] for t in hs]) "Incorrect number of possible orders in price dependent single orders"

    # Accumulate orders per hour
    single_orders = Vector{SingleOrder{eltype(bidlevels[1])}}(undef, nhours(horizon))
    for hour in 1:nhours(horizon)
        single_d = [xᴰ[order,hour] for order = 1:length(model.indices.bids)]
        single_orders[hour] = SingleOrder(hour, xᴵ[hour], single_d, bidlevels[hour])
    end

    block_orders = Vector{BlockOrder{eltype(bidlevels[1])}}()
    for (i,price) in enumerate(model.indices.blockbids)
        for (b,interval) in enumerate(model.indices.hours_per_block)
            ordervolume = xᴮ[i,b]
            if ordervolume > 1e-6
                price = blockbidprice(bidlevels, interval, i)
                push!(block_orders, BlockOrder(tuple(interval[1]-1,interval[end]), ordervolume, price))
            end
        end
    end

    return OrderStrategy(horizon,
                         bidlevels,
                         single_orders,
                         block_orders)
end

function singleorder(strategy::OrderStrategy, hour::Int64)
    if hour > 0 && hour <= nhours(strategy.horizon)
        return strategy.single_orders[hour]
    else
        throw(ArgumentError(string("Selected hour ", hour, " not within horizon 1 to ", nhours(strategy.horizon))))
    end
end
singleorders(strategy::OrderStrategy) = strategy.single_orders
blockorders(strategy::OrderStrategy) = strategy.block_orders
function volumes(strategy::OrderStrategy)
    fromBlock = (block,hour) -> (hour >= block.interval[1]+1 && hour <= block.interval[2]) ? block.volume : 0
    return [strategy.single_orders[hour].independent_volume + maximum(strategy.single_orders[hour].dependent_volumes) + sum([fromBlock(block,hour) for block in strategy.block_orders])
            for hour in 1:hours(strategy.horizon)]
end

function production(strategy::OrderStrategy{T},ρs::AbstractVector) where T <: AbstractFloat
    H = zeros(T,nhours(strategy.horizon))
    accepted_blocks = blockorders(strategy)[findall((o) -> o.price <= mean(ρs[o.interval[1]+1:o.interval[2]]),blockorders(strategy))]
    for h = 1:nhours(strategy.horizon)
        ρ = ρs[h]
        single = singleorder(strategy,h)
        H[h] = independent(single)
        for i = 2:length(single.prices)
            p1 = single.prices[i-1]
            p2 = single.prices[i]
            if p1 <= ρ < p2 || p1 < ρ <= p2
                H[h] += ((ρ - p1)/(p2 - p1))*single.dependent_volumes[i] + ((p2 - ρ)/(p2 - p1))*single.dependent_volumes[i-1]
                break
            end
        end
        blocks_in_hour = findall(o -> (o.interval[1]+1 <= h <= o.interval[2]),accepted_blocks)
        H[h] += sum([order.volume for order in accepted_blocks[blocks_in_hour]])
    end
    return H
end
totalproduction(strategy::OrderStrategy,ρ::AbstractVector) = sum(production(strategy,ρ))

revenue(strategy::OrderStrategy,ρs::AbstractVector) = production(strategy, ρs) .* ρs
totalrevenue(strategy::OrderStrategy,ρs) = production(strategy,ρs) ⋅ ρs

## Print / Plot routines ##
# ===================================================== #
function Base.show(io::IO, ::MIME"text/plain", order::SingleOrder)
    show(io,order)
end

function Base.show(io::IO, order::SingleOrder)
    formatter = (d) -> begin
        if abs(d) <= sqrt(eps())
            "0.0"
        elseif (log10(d) < -2.0 || log10(d) > 3.0)
            @sprintf("%.2e",d)
        else
            @sprintf("%.2f",d)
        end
    end
    if get(io, :multiline, false)
        print(io,@sprintf("Single Order"))
    else
        output = @sprintf("Single Order \nHour: %d \n",order.hour)
        output *= @sprintf("Price Independent Order Volume [MWh/h]: %s \n",formatter(order.independent_volume))
        output *= "Price Dependent Order Volumes [EUR/MWh -> MWh/h]: \n"
        for (price,volume) in zip(order.prices[2:end-1],order.dependent_volumes[2:end-1])
            output *= @sprintf(" %s -> %s \n",formatter(price),formatter(volume))
        end
    end
    print(io,chomp(output))
end

const KTH_colors = [RGB(25/255,84/255,166/255),
                    RGB(157/255,16/255,45/255),
                    RGB(98/255,146/255,46/255),
                    RGB(36/255,160/255,216/255),
                    RGB(228/255,54/255,62/255),
                    RGB(176/255,201/255,43/255),
                    RGB(216/255,84/255,151/255),
                    RGB(250/255,185/255,25/255),
                    RGB(101/255,101/255,108/255),
                    RGB(189/255,188/255,188/255),
                    RGB(227/255,229/255,227/255)]
const KTHBlue = KTH_colors[1]
const KTHRed = KTH_colors[2]
const KTHGreen = KTH_colors[3]
const KTHLBlue = KTH_colors[4]

@recipe f(order::SingleOrder{T}) where T <: AbstractFloat = (order,[])
@recipe f(order::SingleOrder{T},ρ::Real) where T <: AbstractFloat = (order,[ρ])
@recipe f(order::SingleOrder{T},ρs...) where T <: AbstractFloat = (order,[ρs...])
@recipe function f(order::SingleOrder{T}, ρs::AbstractVector) where T <: AbstractFloat
    hour = order.hour
    orderincrement = mean(abs.(diff(order.prices[2:end-1])))
    maxorder = max(order.independent_volume, maximum(order.dependent_volumes))

    line_v = []
    line_p = []
    interp_v = []
    interp_p = []
    for (v,volume) in enumerate(order.dependent_volumes)
        if v == 1 && volume == 0
            append!(line_v,fill(volume,2))
            append!(line_p,[0,order.prices[1]])
        end
        if v > 1
            previous = order.dependent_volumes[v-1]
            if !(abs(volume-previous) <= eps())
                # Horizontal strip
                append!(line_v,[previous,volume])
                append!(line_p,fill(order.prices[v],2))
                # Interpolation line
                if any(ρ -> (order.prices[v-1] <= ρ <= order.prices[v]),ρs)
                    push!(interp_v,[previous,volume])
                    push!(interp_p,[order.prices[v-1],order.prices[v]])
                end
            else
                append!(line_v,fill(volume,2))
                price = [order.prices[v-1],order.prices[v]]
                append!(line_p,price)
            end
        end
    end
    idx = findfirst(v -> v == -500, line_p)
    line_p[idx] = 0.
    idx = findfirst(v -> v == 3000, line_p)
    line_p[idx] = order.prices[end-1]+orderincrement

    v_outcomes = zeros(eltype(order.dependent_volumes),length(ρs))
    if !isempty(ρs)
        # Plot all resulting trading outcomes for the given prices
        for (i,ρ) in enumerate(ρs)
            if ρ < order.prices[1] || ρ > order.prices[end]
                v_outcomes[i] = 0
                continue
            end
            for p = 2:length(order.prices)
                p1 = order.prices[p-1]
                p2 = order.prices[p]
                if p1 <= ρ <= p2
                    v_outcomes[i] = ((ρ - p1)/(p2 - p1))*order.dependent_volumes[p] + ((p2 - ρ)/(p2 - p1))*order.dependent_volumes[p-1]
                end
            end
        end
    end

    # Plot attributes
    padding = maxorder/(2*length(order.dependent_volumes))
    xticks := LinRange(0, maxorder, length(order.prices))
    yticks := 0:orderincrement:order.prices[end-1]
    xlims := (-padding, maxorder+padding)
    ylims := (-orderincrement, order.prices[end-1]+orderincrement)
    formatter := (d) -> @sprintf("%.2f",d)

    title := "Order Curve"
    xlabel := "Order Volume [MWh/h]"
    ylabel := "Price [EUR/MWh]"
    legend := :topleft
    tickfontfamily := "sans-serif"
    guidefontfamily := "sans-serif"
    titlefontfamily := "sans-serif"
    legendfontfamily := "sans-serif"

    # Dashed line
    if !isempty(line_v)
        @series begin
            linestyle := :solid
            linewidth := 2
            linecolor := :black
            label := ""
            line_v, line_p
        end
    end

    # Price independent order
    @series begin
        markercolor --> KTHGreen
        seriestype := :scatter
        label := "Price Independent Order"

        [order.independent_volume],[0]
    end

    if any(v->v >= eps(),order.dependent_volumes)
        # Display the individual orders
        @series begin
            markercolor := KTHBlue
            seriestype := :scatter
            label := "Price Dependent Order"

            order.dependent_volumes[2:end-1], order.prices[2:end-1]
        end
    end
    # Interpolation line
    if !isempty(interp_v)
        @series begin
            linestyle := :dash
            linewidth := 1
            linecolor := :black
            label := ""
            interp_v,interp_p
        end
    end

    # Display the trading outcomes
    if !isempty(v_outcomes)
        @series begin
            markercolor := KTHGreen
            markershape := :diamond
            seriestype := :scatter
            if length(ρs) > 1
                label := "Trading Outcomes"
            else
                label := "Trading Outcome"
            end
            v_outcomes,ρs
        end
    end
end

@recipe f(orders::Vector{SingleOrder{T}}) where T <: AbstractFloat = (orders,[])
@recipe f(orders::Vector{SingleOrder{T}},ρ::Real) where T <: AbstractFloat = (orders,[ρ])
@recipe f(orders::Vector{SingleOrder{T}},ρs...) where T <: AbstractFloat = (orders,[ρs...])
@recipe function f(orders::Vector{SingleOrder{T}}, ρs::AbstractVector) where T <: AbstractFloat
    prices = orders[findmax([o.prices[end-1] for o in orders])[2]].prices
    orderincrement = mean(abs.(diff(prices[2:end-1])))
    ρ_min = !isempty(ρs) ? minimum(ρs) : []
    ρ_max = !isempty(ρs) ? maximum(ρs) : []
    δρ = !isempty(ρs) ? mean(abs.(diff(ρs))) : []
    δρ = !isempty(ρs) ? std(ρs) : []
    highest = !isempty(ρs) ? max(ρ_max, prices[end-1]) : prices[end-1]+2*orderincrement

    rect(x,y,w,h) = Shape(x .+ [0,w,w,0,0],y .+ [0,0,h,h,0])
    formatter = (d) -> begin
        if abs(d) <= sqrt(eps())
            text("0.0",font("sans-serif",:white))
        elseif (log10(d) < -2.0 || log10(d) > 3.0)
            text(@sprintf("%.1e",d),font("sans-serif",:white))
        elseif log10(d) > 2.0
            text(@sprintf("%.1f",d),font("sans-serif",:white))
        else
            text(@sprintf("%.2f",d),font("sans-serif",:white))
        end
    end

    independent_bars = Shape[]
    dependent_bars = Shape[]
    ordervolumes = []
    accepted = Int64[]
    rejected = Int64[]
    for order in orders
        if order.independent_volume >= sqrt(eps())
            push!(independent_bars,rect(order.hour-1,-orderincrement,1.0,orderincrement))
            push!(ordervolumes,(order.hour-0.5,-orderincrement/2,formatter(order.independent_volume)))
        end
        firstpos = findfirst(a -> a >= sqrt(eps()),order.dependent_volumes)
        if firstpos == 0
            continue
        elseif firstpos == 1
            firstpos += 1
        end
        if !any(a -> a >= sqrt(eps()),order.dependent_volumes)
            continue
        end
        start = !isempty(ρs) ? order.prices[1] : 0.0
        stop = 0.0
        splitfound = false
        if !isempty(ρs)
            ρ = ρs[order.hour]
            for i = 2:length(order.dependent_volumes)
                volume = order.dependent_volumes[i]
                previous = order.dependent_volumes[i-1]
                stop = order.prices[i]
                if !splitfound && start <= ρ
                    if start <= ρ <= stop
                        interp_volume = ((ρ - start)/(stop - start))*volume + ((stop - ρ)/(stop - start))*previous
                        stop = ρ
                        start = 0
                        if abs(interp_volume) <= sqrt(eps())
                            break
                        end
                        push!(accepted,length(dependent_bars)+1)
                        push!(dependent_bars,rect(order.hour-1,start,1.0,stop-start))
                        push!(ordervolumes,(order.hour-0.5,(start+stop)/2,formatter(interp_volume)))
                        start = stop
                        splitfound = true
                        break
                    else
                        if !(abs(volume-previous) <= sqrt(eps()))
                            start = stop
                        end
                    end
                end
            end
        else
            for i = 2:length(order.dependent_volumes)
                volume = order.dependent_volumes[i]
                stop,next = if i == length(order.dependent_volumes)
                    order.prices[end-1]+orderincrement,volume
                else
                    order.prices[i],order.dependent_volumes[i+1]
                end
                if !(abs(next-volume) <= sqrt(eps())) || i == length(order.dependent_volumes)
                    push!(dependent_bars,rect(order.hour-1,start,1.0,stop-start))
                    push!(ordervolumes,(order.hour-0.5,(start+stop)/2,formatter(volume)))
                    start = stop
                end
            end
        end

        if !isempty(ρs)
            stop = prices[end]
            push!(rejected,length(dependent_bars)+1)
            push!(dependent_bars,rect(order.hour-1,start,1.0,stop-start))
        end
    end

    height = highest+2*orderincrement
    # Independent Volume label
    push!(independent_bars,rect(25,-orderincrement,1.0,0.5*height))
    push!(ordervolumes,(25.5,0.25*height-orderincrement,text("Independent Volume [MWh]",font("sans-serif", :white, rotation = -90.))))
    # Dependent Volume label
    push!(accepted,length(dependent_bars)+1)
    push!(dependent_bars,rect(25,0.5*height-orderincrement,1.0,0.5*height))
    push!(ordervolumes,(25.5,0.75*height-orderincrement,text("Dependent Volumes [Mwh]",font("sans-serif", :white, rotation = -90.))))

    # Plot attributes
    xticks := collect(0:24)
    xlims := (-1,26)
    ylims --> (-orderincrement,highest+orderincrement)
    if !isempty(ρs)
        yticks := 0:mean([δρ,orderincrement]):highest+orderincrement
    else
        yticks := 0:orderincrement:highest
    end
    tickfontfamily := "sans-serif"
    guidefontfamily := "sans-serif"
    titlefontfamily := "sans-serif"
    legendfontfamily := "sans-serif"

    legend := :topleft
    annotations := ordervolumes
    yformatter := (d) -> @sprintf("%.2f",d)

    title := "Single Orders"
    xlabel := "Hour"
    ylabel := "Price [EUR/MWh]"

    # Always show the price independent orders
    @series begin
        seriestype := :shape
        if !isempty(ρs)
            seriescolor := KTHGreen
        else
            seriescolor := KTHLBlue
        end
        label := ""
        independent_bars
    end

    if !isempty(ρs)
        # Accepted orders
        @series begin
            seriestype := :shape
            seriescolor := KTHGreen
            label := "Accepted Orders"
            dependent_bars[accepted]
        end
        # Rejected orders
        @series begin
            seriestype := :shape
            seriescolor := KTHRed
            label := "Rejected Orders"
            dependent_bars[rejected]
        end
        # Show the price curve
        @series begin
            seriestype := :scatter
            seriescolor := KTHGreen
            label := "Market price"
            collect(0.5:1:23.5),ρs
        end
        @series begin
            seriestype := :path
            seriescolor := KTHGreen
            label := ""
            collect(0.5:1:23.5),ρs
        end
    else
        # Show all price dependent orders if no prices given
        @series begin
            seriestype := :shape
            seriescolor := KTHBlue
            label := ""
            dependent_bars
        end
    end
end

function Base.show(io::IO, ::MIME"text/plain", order::BlockOrder)
    show(io,order)
end

function Base.show(io::IO, order::BlockOrder)
    formatter = (d) -> begin
        if abs(d) <= sqrt(eps())
            "0.0"
        elseif (log10(d) < -2.0 || log10(d) > 3.0)
            @sprintf("%.2e",d)
        else
            @sprintf("%.2f",d)
        end
    end
    if get(io, :multiline, false)
        print(io,"Block Order")
    else
        output = "Block Order\n"
        output *= string("Interval: ",order.interval[1]," - ",order.interval[2],"\n")
        if order.volume >= eps()
                        output *= @sprintf("Order Volume [MWh/h]: %s \n",formatter(order.volume))
        else
            output *= "Order Volume [MWh/h]: 0\n"
        end
        output *= @sprintf("Price [EUR/MWh]: %s",formatter(order.price))
    end
    print(io,output)
end

@recipe f(order::BlockOrder{T}) where T <: AbstractFloat = order,[]
@recipe function f(order::BlockOrder{T},ρs::AbstractVector) where T <: AbstractFloat
    ρ_min = !isempty(ρs) ? minimum(ρs) : []
    ρ_max = !isempty(ρs) ? maximum(ρs) : []
    δρ = !isempty(ρs) ? mean(abs.(diff(ρs))) : []
    p_max = order.price .+ δρ
    highest = !isempty(ρs) ? max(p_max,ρ_max) : p_max
    annotation_font = font("sans-serif",:white)
    formatter = (d) -> begin
        if abs(d) <= sqrt(eps())
            text("0.0",font("sans-serif",:white))
        elseif (log10(d) < -2.0 || log10(d) > 3.0)
            text(@sprintf("%.2e",d),annotation_font)
        elseif log10(d) > 2.0
            text(@sprintf("%.1f",d),annotation_font)
        else
            text(@sprintf("%.2f",d),annotation_font)
        end
    end

    price = order.price
    width = order.interval[2]-order.interval[1]
    height = order.price/10
    color = KTHRed
    status = "Rejected Order"
    padding = !isempty(ρs) ? 1.5 : 2*(order.interval[2]-order.interval[1])/24

    if !isempty(ρs)
        ρ̅ = mean(ρs[order.interval[1]+1:order.interval[2]])
        if order.price <= ρ̅
            price = ρ̅
            color = KTHGreen
            status = "Accepted Order"
        else
            δρ *= 4
        end
    end

    # Plot attributes
    ylims := (0,2*price)
    xticks := [order.interval...]
    if !isempty(ρs)
        xlims := -1:1:24
        ylims := (ρ_min-δρ,highest+δρ)
        if δρ >= eps()
            yticks := ρ_min:δρ:highest+δρ
        else
            yticks := []
        end
    else
        xlims := (order.interval[1]-1, order.interval[2]+1)
        yticks := []
    end
    title := "Block Order"
    tickfontfamily := "sans-serif"
    guidefontfamily := "sans-serif"
    titlefontfamily := "sans-serif"
    legendfontfamily := "sans-serif"
    title := "Block Order"
    xlabel := "Hour"
    ylabel := !isempty(ρs) ? "Price [EUR/MWh]" : ""
    yformatter --> (d) -> @sprintf("%.2f",d)
    annotations --> [(order.interval[1] .+ padding, price,text(@sprintf("%.2f [EUR/MWh]",order.price), annotation_font)),
                     (order.interval[2] .- padding, price,text(@sprintf("%.2f [MWh/h]",order.volume), annotation_font))]

    if !isempty(ρs)
        # Accepted/Rejected orders
        @series begin
            seriestype := :shape
            seriescolor := color
            label := status
            Shape(order.interval[1] .+ [0.,width,width,0,0], price .- height/2 .+ [0.,0,height,height,0])
        end
        # Show the price curve
        @series begin
            seriestype := :scatter
            seriescolor := KTHGreen
            label := "Market price"
            collect(0.5:1:23.5), ρs
        end
        @series begin
            seriestype := :path
            seriescolor := KTHGreen
            label := ""
            collect(0.5:1:23.5), ρs
        end
    else
        @series begin
            seriestype := :shape
            seriescolor := KTHBlue
            label := ""
            Shape(order.interval[1] .+ [0,width,width,0,0], order.price .- height/2 .+ [0.,0,height,height,0])
        end
    end
end

@recipe f(orders::Vector{BlockOrder{T}}) where T <: AbstractFloat = orders,[]
@recipe function f(orders::Vector{BlockOrder{T}},ρs::AbstractVector) where T <: AbstractFloat
    if !isempty(ρs) && length(ρs) != 24
        throw(ArgumentError("Need to supply a price vector for the whole day when showing block orders"))
    end
    ρ_min = !isempty(ρs) ? minimum(ρs) : []
    ρ_max = !isempty(ρs) ? maximum(ρs) : []
    δρ = !isempty(ρs) ? std(ρs) : []
    p_max = maximum([order.price for order in orders]) .+ δρ
    orderincrement = mean(abs.(diff([order.price for order in orders])))
    increment = !isempty(ρs) ? mean([δρ,orderincrement]) : orderincrement
    highest = !isempty(ρs) ? max(p_max,ρ_max) : p_max
    levels = collect(1:1:length(orders))
    sorted = sort(orders,by=(order)->order.price)
    annotation_font = font("sans-serif",:white)
    formatter = (d) -> begin
        if abs(d) <= sqrt(eps())
            text("0.0",font("sans-serif",:white))
        elseif (log10(d) < -2.0 || log10(d) > 3.0)
            text(@sprintf("%.2e",d),annotation_font)
        elseif log10(d) > 2.0
            text(@sprintf("%.1f",d),annotation_font)
        else
            text(@sprintf("%.2f",d),annotation_font)
        end
    end
    plain_orders = Shape[]
    accepted_blocks = Shape[]
    rejected_blocks = []
    rejected_prices = []
    orderinfos = []

    rect(x,y,w,h) = Shape(x .+ [0,w,w,0,0],y .+ [0,0,h,h,0])
    if !isempty(ρs)
        accepted = findall(o -> o.price <= mean(ρs[o.interval[1]+1:o.interval[2]]),orders)
        if isempty(accepted)
            δρ *= 4
        end
        accepted_orders = [BlockOrder(order.interval,sum([same.volume for same in orders[accepted][findall(o -> o.interval == order.interval,orders[accepted])]]),order.price) for order in unique(o -> o.interval,orders[accepted])]
        rejected_orders = orders[setdiff(1:length(orders),accepted)]
        for order in accepted_orders
            width = order.interval[2]-order.interval[1]
            ρ̅ = mean(ρs[order.interval[1]+1:order.interval[2]])
            push!(accepted_blocks,rect(order.interval[1],ρ̅-increment/4,width, increment/2))
            push!(orderinfos,(order.interval[1]+0.5,ρ̅,formatter(ρ̅)))
            push!(orderinfos,(order.interval[2]-0.5,ρ̅,formatter(order.volume)))
        end
        for order in rejected_orders
            append!(rejected_blocks,[order.interval...])
            push!(rejected_blocks,NaN)
            append!(rejected_prices,fill(order.price,2))
            push!(rejected_prices,NaN)
        end
    else
        for (level,order) in zip(levels,sorted)
            width = order.interval[2]-order.interval[1]
            push!(plain_orders,rect(order.interval[1],level-0.5,width,1.0))
            push!(orderinfos,(order.interval[1]+0.5,level,formatter(order.price)))
            push!(orderinfos,(order.interval[2]-0.5,level,formatter(order.volume)))
        end
    end

    # Legend
    if !isempty(ρs)
        append!(orderinfos,[(20.5,highest+3*δρ/4,text("Price [EUR/MWh]",annotation_font)),(24.5,highest+3*δρ/4,text("Volume [MWh/h]",annotation_font))])
    else
        append!(orderinfos,[(20.5,length(orders)+2.5,text("Price [EUR/MWh]",annotation_font)),(24.5,length(orders)+2.5,text("Volume [MWh/h]",annotation_font))])
    end
    legendline = if !isempty(ρs)
        rect(19,highest+δρ/2,7,increment/2)
    else
        rect(19,length(orders)+2,7,1.0)
    end

    # Plot attributes
    xlims := (-1,26)
    xticks := 0:1:24
    if !isempty(ρs)
        ylims := (ρ_min-increment,highest+increment)
        if δρ >= eps()
            yticks := ρ_min:increment:highest+δρ
        else
            yticks := []
        end
    else
        ylims := (-1,length(orders)+3)
        yticks := []
    end
    yformatter := (d) -> @sprintf("%.2f",d)
    tickfontfamily := "sans-serif"
    guidefontfamily := "sans-serif"
    titlefontfamily := "sans-serif"
    legendfontfamily := "sans-serif"
    color_palette := KTH_colors
    title := "Block Orders"
    xlabel := "Hour"
    ylabel := "Price [EUR/MWh]"
    annotations --> orderinfos

    # Legend
    if !isempty(ρs)
        @series begin
            seriestype := :shape
            seriescolor := KTHGreen
            label := ""
            legendline
        end
    else
        @series begin
            seriestype := :shape
            seriescolor := KTHBlue
            label := ""
            legendline
        end
    end

    # Block orders
    if !isempty(ρs)
        # Display which orders are accepted/rejected
        # Accepted orders
        if !isempty(accepted_blocks)
            @series begin
                seriestype := :shape
                seriescolor := KTHGreen
                label := "Accepted Orders"
                accepted_blocks
            end
        end
        # Rejected orders
        if !isempty(rejected_blocks)
            @series begin
                linewidth := 4
                linecolor := KTHRed
                seriestype := :shape
                seriescolor := KTHRed
                label := "Rejected Orders"
                rejected_blocks,rejected_prices
            end
        end
        # Show the price curve
        @series begin
            seriestype := :scatter
            seriescolor := KTHGreen
            label := "Market price"
            collect(0.5:1:23.5),ρs
        end
        @series begin
            seriestype := :path
            seriescolor := KTHGreen
            label := ""
            collect(0.5:1:23.5),ρs
        end
    else
        # Show all orders
        @series begin
            seriestype := :shape
            seriescolor := KTHBlue
            label := ""
            plain_orders
        end
    end
end

function show(io::IO, ::MIME"text/plain", strategy::OrderStrategy)
    show(io,strategy)
end

function show(io::IO, strategy::OrderStrategy)
    if get(io, :multiline, false)
        print(io,"Order Strategy")
    else
        println(io,"Order Strategy")
        println(io,"Price levels:")
        Base.print_matrix(io,strategy.bidlevels)
    end
end

@recipe f(strategy::OrderStrategy{T}) where T <: AbstractFloat = strategy,[]
@recipe function f(strategy::OrderStrategy{T},ρs::AbstractVector) where T <: AbstractFloat
    if !isempty(ρs) && length(ρs) != 24
        throw(ArgumentError("Need to supply a price vector for the whole day when showing block orders"))
    end
    bidlevels = strategy.bidlevels
    orderincrement = mean(abs.(reduce(vcat, [diff(prices) for prices in bidlevels])))

    if !isempty(ρs)
        legend := :topleft
    end
    layout := (2,1)
    link := :x

    @series begin
        title := "Single Orders"
        subplot := 1
        strategy.single_orders,ρs
    end

    @series begin
        title := "Block Orders"
        subplot := 2
        strategy.block_orders,ρs
    end
end

function strategy(model::StochasticHydroModel; variant = :rp)
    status(model; variant = variant) == :Planned || error("Hydro model has not been planned yet")
    return OrderStrategy(model)
end


text:
#Ipopt.jl example - This is taken from the documentation available on:
# http://ipoptjl.readthedocs.org/en/latest/ipopt.html#

# Install Package if you don't have it
Pkg.add("Ipopt")

tic()
# Load the package
using Ipopt

# HS071
# min x1 * x4 * (x1 + x2 + x3) + x3
# st  x1 * x2 * x3 * x4 >= 25
#     x1^2 + x2^2 + x3^2 + x4^2 = 40
#     1 <= x1, x2, x3, x4 <= 5
# Start at (1,5,5,1)
# End at (1.000..., 4.743..., 3.821..., 1.379...)

## Set the number of arguments to solve for and bounds for these
## Aside for new Julia users: the '.0' after each number is important - it tells Julia
## that the vectors x_L, etc are of Float type.
## Note that Julia is a functional language, with methods (e.g. a function) working for the types
## for which it is defined (e.g. integer, Float). So, a function defined for an integer will not
## work if you enter a Float variable, and vice-versa.
n = 4
x_L = [1.0, 1.0, 1.0, 1.0]
x_U = [5.0, 5.0, 5.0, 5.0]

# Set the number of constraints and bounds on these - Note that equality
# constraints can be accommodated for by setting the upper bound = lower bound
m = 2
g_L = [25.0, 40.0]
g_U = [2.0e19, 40.0]

# Set the objective function - Returns the value of the objective function at the current solution x
# Note that the '::Vector{Float64}' below defines the type of x. This is a Julia-specific thing, and
# need not be defined. The function would be created even without including the type for x.
function eval_f(x::Vector{Float64})
  return x[1] * x[4] * (x[1] + x[2] + x[3]) + x[3]
end

# Set the constraint set - Sets the value of the constraint functions g at the current solution x
# Note that there is no need to create a separate vector g in which to save the constraints.
# The values of g are set 'in-place'
function eval_g(x::Vector{Float64}, g::Vector{Float64})
  g[1] = x[1]   * x[2]   * x[3]   * x[4]
  g[2] = x[1]^2 + x[2]^2 + x[3]^2 + x[4]^2
end

# Gradient of the objective function - Sets the value of the gradient of the objective function at the current solution x
# Note that there is no need to create a vector grad_f to fill in. The values of grad_f are set 'in-place'.
function eval_grad_f(x::Vector{Float64}, grad_f::Vector{Float64})
  grad_f[1] = x[1] * x[4] + x[4] * (x[1] + x[2] + x[3])
  grad_f[2] = x[1] * x[4]
  grad_f[3] = x[1] * x[4] + 1
  grad_f[4] = x[1] * (x[1] + x[2] + x[3])
end

# Gradient of the constraint set (Jacobian matrix). There are two ways of setting this:
# (1) Under the ':Structure' mode, indicate the rows and columns with non-zero values (i.e. tell Ipopt
# about the sparsity structure of the matrix)
# (2) Undet the ':Values' mode, tell it the actual values.
# Each non-zero value is indexed, starting from 1.
###########################################
# General specification is as follows:
###########################################
# function eval_jac_g(
#  x::Vector{Float64},         # Current solution
#  mode,                       # Either :Structure or :Values
#  rows::Vector{Int32},        # Sparsity structure - row indices
#  cols::Vector{Int32},        # Sparsity structure - column indices
#  values::Vector{Float64})    # The values of the Hessian

#  if mode == :Structure
    # rows[...] = ...
    # ...
    # cols[...] = ...
#  else
    # values[...] = ...
#  end
# end
#############################################
function eval_jac_g(x, mode, rows, cols, values)
  if mode == :Structure
    # Constraint (row) 1
    rows[1] = 1; cols[1] = 1
    rows[2] = 1; cols[2] = 2
    rows[3] = 1; cols[3] = 3
    rows[4] = 1; cols[4] = 4
    # Constraint (row) 2
    rows[5] = 2; cols[5] = 1
    rows[6] = 2; cols[6] = 2
    rows[7] = 2; cols[7] = 3
    rows[8] = 2; cols[8] = 4
  else
    # Constraint (row) 1
    values[1] = x[2]*x[3]*x[4]  # 1,1
    values[2] = x[1]*x[3]*x[4]  # 1,2
    values[3] = x[1]*x[2]*x[4]  # 1,3
    values[4] = x[1]*x[2]*x[3]  # 1,4
    # Constraint (row) 2
    values[5] = 2*x[1]  # 2,1
    values[6] = 2*x[2]  # 2,2
    values[7] = 2*x[3]  # 2,3
    values[8] = 2*x[4]  # 2,4
  end
end

# Hessian of the Lagrangian - this is an optional object. In terms of the Julia syntax, it is
# similar to the Jacobian.
# The general definition and types for the objects in the function below are:
# function eval_h(
# x::Vector{Float64},         # Current solution
# mode,                       # Either :Structure or :Values
# rows::Vector{Int32},        # Sparsity structure - row indices
# cols::Vector{Int32},        # Sparsity structure - column indices
# obj_factor::Float64,        # Lagrangian multiplier for objective
# lambda::Vector{Float64},    # Multipliers for each constraint
# values::Vector{Float64})    # The values of the Hessian
# if....else...end            # the function text goes here
# end                         # close the function

function eval_h(x, mode, rows, cols, obj_factor, lambda, values)
  if mode == :Structure
    # Symmetric matrix, fill the lower left triangle only
    idx = 1
    for row = 1:4
      for col = 1:row
        rows[idx] = row
        cols[idx] = col
        idx += 1
      end
    end
  else
    # Again, only lower left triangle
    # Objective
    values[1] = obj_factor * (2*x[4])  # 1,1
    values[2] = obj_factor * (  x[4])  # 2,1
    values[3] = 0                      # 2,2
    values[4] = obj_factor * (  x[4])  # 3,1
    values[5] = 0                      # 3,2
    values[6] = 0                      # 3,3
    values[7] = obj_factor * (2*x[1] + x[2] + x[3])  # 4,1
    values[8] = obj_factor * (  x[1])  # 4,2
    values[9] = obj_factor * (  x[1])  # 4,3
    values[10] = 0                     # 4,4

    # First constraint
    values[2] += lambda[1] * (x[3] * x[4])  # 2,1
    values[4] += lambda[1] * (x[2] * x[4])  # 3,1
    values[5] += lambda[1] * (x[1] * x[4])  # 3,2
    values[7] += lambda[1] * (x[2] * x[3])  # 4,1
    values[8] += lambda[1] * (x[1] * x[3])  # 4,2
    values[9] += lambda[1] * (x[1] * x[2])  # 4,3

    # Second constraint
    values[1]  += lambda[2] * 2  # 1,1
    values[3]  += lambda[2] * 2  # 2,2
    values[6]  += lambda[2] * 2  # 3,3
    values[10] += lambda[2] * 2  # 4,4
  end
end

# Creates and returns an IpoptProblem with the given options. Raises error if something goes wrong during construction.
# If you do not provide a callback for the Hessian, you must set the Hessian approximation option:
# addOption(prob, "hessian_approximation", "limited-memory")
prob = createProblem(n, x_L, x_U, m, g_L, g_U, 8, 10,
                     eval_f, eval_g, eval_grad_f, eval_jac_g, eval_h)

# Set starting solution
prob.x = [1.0, 5.0, 5.0, 1.0]

# Solve
status = solveProblem(prob)

# Display the solution
println(Ipopt.ApplicationReturnStatus[status])
println(prob.x)
println(prob.obj_val)

toc()

max_stars_repo_name: floswald/ucl-econ-julia

max_issues_repo_name: floswald/ucl-econ-julia

max_forks_repo_name: floswald/ucl-econ-julia


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -33.59531116254298 <= q <= 31.312008728451346)

                     add_NL_constraint(m, :(log(1 + exp(0.7464435003747418 + -0.7676428720221384*log(1 + exp(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4]))) + -0.5937285886485864*log(1 + exp(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4]))) + 0.70936536775501*log(1 + exp(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4]))) + 0.9647291374973435*log(1 + exp(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4]))) + 0.527810271934432*log(1 + exp(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))))) + log(1 + exp(-0.35636357579153266 + 0.8665784119884141*log(1 + exp(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4]))) + 0.03175667478129274*log(1 + exp(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4]))) + -0.7799384758999897*log(1 + exp(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4]))) + 0.4104070349373172*log(1 + exp(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4]))) + -0.7416668360909755*log(1 + exp(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))))) + log(1 + exp(0.5038932713674513 + 0.7223377417739489*log(1 + exp(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4]))) + -0.17110213406275054*log(1 + exp(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4]))) + 0.7726394930743403*log(1 + exp(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4]))) + 0.6118790711423934*log(1 + exp(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4]))) + 0.16250007841990932*log(1 + exp(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))))) + log(1 + exp(-0.7190092552996492 + 0.6952008118484776*log(1 + exp(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4]))) + 0.975536624348428*log(1 + exp(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4]))) + 0.03367235210941821*log(1 + exp(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4]))) + 0.8059867473566515*log(1 + exp(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4]))) + 0.5874603423029647*log(1 + exp(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))))) + log(1 + exp(-0.5715134646750415 + 0.6325466436773857*log(1 + exp(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4]))) + -0.2633742432221431*log(1 + exp(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4]))) + -0.3786600783719858*log(1 + exp(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4]))) + 0.5923927477083897*log(1 + exp(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4]))) + 0.9114166052727355*log(1 + exp(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
"Computes Gregory-type quadrature rules."
module Gregory

export integrate

"""
    B = calcBernoulli(n)

Return the `n`th Bernoulli number.  Follows the convention that `n=1`
Bernoulli number is equal to -1/2.
"""
function calcBernoulli(n)
    A = Vector{Rational{BigInt}}(undef, n + 1)
    for m = 0:n
        A[m + 1] = 1 // (m + 1)
        for j = m:-1:1
            A[j] = j * (A[j] - A[j + 1])
        end
    end
    if n == 1
        return -A[1]
    else
        return A[1]
    end
end

"""
    w = calcBoundaryWeights(numbnd, order)

Compute the boundary weights for a Gregory rule of order `order` that using
`numbnd` boundary points.
"""
function calcBoundaryWeights(numbnd::Int, order::Int)
    @assert( numbnd >= order-1,
             "number of boundary nodes must be larger than order-1")
    # Form the linear system to solve
    A = zeros(Rational{BigInt}, order-1, numbnd)
    b = zeros(Rational{BigInt}, order-1)
    for j = 1:order-1
        for i = 1:numbnd
            A[j,i] = j*(numbnd - (i-1))^(j-1)
        end
        b[j] = numbnd^j - ((-1)^j)*calcBernoulli(j)
    end
    w = A\b
    return w
end

"""
    val = integrate(data, numbnd, order[, h=1.0])

Integrate uniformly spaced `data` using a Gregory quadrature rule of order
`order` with `numbnd` boundary points.  The mesh spacing is defined by `h`.
"""
function integrate(data::Array{T,1}, numbnd::Int, order::Int;
                   h::T=one(T)) where {T<:Number}
    @assert( size(data,1) >= 2*numbnd,
             "number of data points must be larger than 2*numbnd")
    w = calcBoundaryWeights(numbnd, order)
    integral = zero(T)
    for i = 1:numbnd
        integral += w[i]*(data[i] + data[end-i+1])
    end
    for i = numbnd+1:size(data,1)-numbnd
        integral += data[i]
    end
    return integral*h
end

end # module Gregory

max_stars_repo_name: jehicken/high-order-and-chaos

max_issues_repo_name: jehicken/high-order-and-chaos

max_forks_repo_name: jehicken/high-order-and-chaos


text:
# Wrap everything with a module to allow redefition of type
module MC

"""
Composite type to represent a spin state
"""
mutable struct SpinState
    num_spins::Int
    s::Array{Int8,1}
    energy::Int
    tot_mag::Int
end

"""
Energy
"""
function energy(s)
    n = length(s)
    - sum((s[i] * s[ifelse(i == n, 1, i+1)] for i in 1:n))
end

"""
Total magnetization
"""
total_magnetization(s) = sum(s)

"""
Constructor
"""
function SpinState(s)
    ss = SpinState(length(s), copy(s), energy(s), total_magnetization(s))
    sanity_check(ss)
    ss
end

"""
Sanity check
"""
function sanity_check(ss)
    @assert energy(ss.s) == ss.energy
    @assert total_magnetization(ss.s) == ss.tot_mag
end

"""
Take an object of SpinState as an input and update it in place.
"""
function update!(ss, β, niters, rng)
    min_h = -2
    max_h = 2
    s = ss.s
    n = ss.num_spins
    prob = [1/(1+exp(-2*β*h)) for h in min_h:max_h]
    for iter in 1:niters, i in 1:n
        sl = s[ifelse(i == 1, n, i-1)]
        sr = s[ifelse(i == n, 1, i+1)]
        # h = -2, 0, 2
        h = sl + sr
        si_old = s[i]
        s[i] = ifelse(rand(rng) < prob[h-min_h+1], +1, -1)

        # Update observables with O(1) operations
        ss.energy += (si_old - s[i]) * h
        ss.tot_mag += (s[i] - si_old)
    end
end

end