text:
module senstools
export OSM!, SCM!
#######################################################################
#   Function for performing different types of sensitivity analysis   #
#######################################################################

type SCMdata
  nopars            :: Int
  parameters        :: Array{Float64, 1}
  names             :: Array{String, 1}
  time              :: Array{Float64, 1}
  sensitivities     :: Array{Float64, 2}
  twonorms          :: Array{Float64, 1}
  ordering          :: Array{Int, 1}
  subset            :: Array{Int, 1}
  nochosen          :: Array{Int, 1}
  twonormfactor     :: Array{Float64, 1}
  normval           :: Array{Float64, 1}
  precision         :: Array{Float64, 1}
  correlation_limit :: Array{Float64, 1}
end

type OSMdata
  parameters    :: Array{Float64, 1}
  names         :: Array{String, 1}
  time          :: Array{Float64, 1}
  sensitivities :: Array{Float64, 2}
  E             :: Array{Float64,1}
  subset        :: Array{Int, 1}
  available     :: Array{Int, 1}
  angles        :: Array{Float64, 1}
  I             :: Array{Float64, 1}
  tolerance     :: Float64
end

type Morrisdata
  level         :: Int
  points        :: Int
  no_parameters :: Int
  f             :: Function
  print         :: Bool
  names         :: Array{String, 1}
  mu
  s2
end

type Soboldata
  points        :: Int              # Number of points to use (M)
  no_parameters :: Int              # Number of parameters
  f             :: Function         # Function to be investigated. Must take array of parameters as only input.
  print         :: Bool             # Print results?
  names         :: Array{String, 1} # Names of parameters
  Si                                # Resulting Si values
  STi                               # Resulting STi values
end

function SCM!(
  Data :: SCMdata)

  """ Strucural Correlation Matrix Method from Olufsen et. al REF"""

  run(`figlet -w 800 Structural Correlations`)
  @printf("\n\nSTATUS: Started Structural Correlation investigation.\n")

  SCM_twonorm_order!(
    Data :: SCMdata);

  SCM_select_from_twonorm!(
    Data :: SCMdata);

  flag = 1
  while flag > 0
    flag = SCM_max_correlation_reduction!(
      Data :: SCMdata)
  end

  SCM_print_results(
    Data :: SCMdata)

  return int32(0);
end

function SCM_twonorm_order!(
  Data :: SCMdata)

  @printf("STATUS: Calculating norms for sensitivities...\n")
  Data.twonorms[1:Data.nopars]   = float64([norm(Data.sensitivities[:,j],Data.normval[1]) for j in 1:Data.nopars]) :: Array{Float64,1};
  Data.twonorms[1:Data.nopars]  /= maximum(Data.twonorms);
  Data.ordering[1:Data.nopars]   = sortperm(Data.twonorms)[end:-1:1]                                               :: Array{Int,1};

  return int32(0);
end

function SCM_select_from_twonorm!(
  Data :: SCMdata)

  #Data.nochosen[1] = length(Data.subset[Data.twonorms[Data.subset].>(1.0/Data.twonormfactor[1])])
  Data.subset = Data.subset[Data.twonorms[Data.subset].>(1.0/Data.twonormfactor[1])]
  Data.nochosen[1] = length(Data.subset)

  return int32(0)
end

function SCM_max_correlation_reduction!(
  Data :: SCMdata;
  condlimit = 1e10)

  @printf("STATUS: Calculating correlations...\n")
  @printf("----------------------------------------------\n")

  subset = Data.subset[1:Data.nochosen[1]]
  STS = Data.sensitivities[:,subset]'*Data.sensitivities[:,subset]

  ###########################################
  #  Check if we can calculate the inverse  #
  ###########################################

  c = try inv(STS)
    catch error("STS can not be inverted. Perhaps too many parameters have been chosen,
    or the ones chosen have linearly dependent sensitivities.")
    end

  # The (STS)^-1 * STS shoulde be I. Consdire the difference from the identity matrix.
  I = c*STS-eye(c)

  # If the largest error is larger than the precision we want, return an error.
  if maximum(abs(I)) > Data.precision[1]
    error("Too large imprecission in covariance matrix. Perhaps too many parameters have been chosen,
    or the ones chosen have linearly dependent sensitivities.\n")
  end

  C = zeros(c);

  #C = [c[i,j]/sqrt(c[i,i]*c[j,j]) for i in 1:Data.nochosen[1], j in 1:Data.nochosen[1]]-eye(c)
  for i in 1:Data.nochosen[1]
    for j in 1:Data.nochosen[1]
      C[i,j] = c[i,j]/sqrt(c[i,i]*c[j,j]);
    end
  end
  C -= eye(C);

  if maximum(C) > Data.correlation_limit[1]
    (i,j) = max_index(abs(C))

    @printf("CORRELATION FOUND:\t %2.2f\n", C[i,j])
    @printf("%s (%2.2f) and %s (%2.2f)\n",
      Data.names[Data.subset[i]],
      Data.twonorms[Data.subset[i]],
      Data.names[Data.subset[j]],
      Data.twonorms[Data.subset[j]])

    if Data.twonorms[Data.subset[i]] > Data.twonorms[Data.subset[j]]
      removed = j
    else
      removed = i
    end
    @printf("Removed parameter:\t %s\n", Data.names[Data.subset[removed]])
    @printf("----------------------------------------------\n")
    splice!(Data.subset,removed)
    Data.nochosen[1] -= 1

    return int32(1)

  else
    @printf("----------------------------------------------\n")
    @printf("No correlation larger than %2.2f have been identified.\n",Data.correlation_limit[1])
    @printf("----------------------------------------------\n")
    (i,j) = max_index(abs(C))
    return int32(0)
  end
end

function SCM_print_results(
  Data :: SCMdata)

  @printf("\n\nSTATUS: Structural Correlations analysis finished.\n")
  @printf("Found parameters are:\n")
  @printf("----------------------------------------------\n")
  @printf("2-norm\t Name\n")
  @printf("----------------------------------------------\n")
  #for i in 1:Data.nochosen[1]
  for i in sortperm(Data.twonorms[Data.subset])[end:-1:1]
    @printf("%2.4f\t %s\n", Data.twonorms[Data.subset[i]], Data.names[Data.subset[i]])
  end
  @printf("----------------------------------------------\n")
end

function max_index(
  Ain :: Array{Float64,2})

  A = copy(Ain)
  B = zeros(Int,size(A))

  # Order each column
  for i in 1:size(A,1)
    B[:,i] = sortperm(A[:,i])
    A[:,i] = A[B[:,i],i]
  end
  # Save index for largest row in each column
  maxrows = copy(B[end,:])

  # Transpose, and repeat
  A = A'
  B = zeros(Int,size(A))
  for i in 1:size(A,1)
    B[:,i] = sortperm(A[:,i])
    A[:,i] = A[B[:,i],i]
  end

  # Last component of the last column of B is the column number we are looking for
  j = copy(B[end,end])
  i = maxrows[j]

  return (i,j) :: (Int, Int);
end

function OSM_calcE!(
  Data :: OSMdata)

  """ Function to be called by OSM! to calculate the E-rank
  of the parameters. """

    #p = numpy.size(basis,1)
    #STS =  numpy.dot(numpy.transpose(basis),basis)
    #[l, Q] = numpy.linalg.eig(STS);

  STS = Data.sensitivities'*Data.sensitivities
  (l, Q) = eig(STS)

    #E = abs(Q)
    #E = numpy.array(numpy.dot(numpy.abs(Q),numpy.abs(l)))/numpy.sum(numpy.abs(l))
    #E = E[0]

  Data.E = abs(Q)*abs(l)/sum(abs(l))

  return int(0)
end

function OSM_select_par!(
  index :: Int,
  Data :: OSMdata;
  I = NaN :: Float64,
  angle = NaN :: Float64)
  """ Function to be called by OSM! to select a parameter for the subset. """

  push!(Data.subset, Data.available[index])
  push!(Data.I,I)
  push!(Data.angles,angle)
  splice!(Data.available,index)

  return int(0)
end

function OSM_calculate_I_and_select!(
  Data :: OSMdata)

  d = zeros(Float64,size(Data.available))
  for (i,j) in enumerate(Data.available)
    a = Data.sensitivities[:,Data.subset]\Data.sensitivities[:,j]
    s = Data.sensitivities[:,Data.subset]*a

    ########################################################
    #  Calculate the I-index for the remaining parameters  #
    ########################################################

    # Consinder the easy options first

    if norm(Data.sensitivities[:,j]) < Data.tolerance
      # If the column has near 0 2-norm, set d[j] = 0
      d[i] = 0.0

    elseif norm(a) < Data.tolerance
      # If the closests fit has all zero components, set d[j] = 1
      d[i] = 1.0

    else
      # Calculate the projection length and denominator
      proj_length = dot(s,Data.sensitivities[:,j])
      denom = norm(s)*norm(Data.sensitivities[:,j])

      if abs(proj_length/denom) > 1.0
        # This could happen with numerical error, but would only be slightly larger
        d[i] = 0
      else
        d[i] = sin(acos(proj_length/denom))
      end

    end
  end

  I = Data.E[Data.available].*d

  #####################################
  #  Now for selecting the parameter  #
  #####################################

  if norm(I) == 0.0
    # If all I-values are 0
    println("All I-values are 0")

    if norm(Data.E[Data.available]) == 0
      # If E-value of remaining parameters are 0, select first from the list
      selected_element = sortperm(Data.available)[1]

    else
      # Select available parameter with largest E-value
      selected_element = sortperm(E[Data.available])[end]

    end

  else
    # Otherwise, select parameter with largest I-value
    selected_element = sortperm(I)[end]

  end

  OSM_select_par!(
    selected_element :: Int,
    Data :: OSMdata;
    I=I[selected_element],
    angle=d[selected_element])

  return int(0)
end

function OSM_print_results(
  Data :: OSMdata)

  @printf("Results for OSM analysis:\n")
  @printf("---------------------------------------------------------------------------------------------------\n")
  @printf("Rank\t\t I\t\t E\t\t d\t\t Name\n")
  @printf("---------------------------------------------------------------------------------------------------\n")
  for (i,j) in enumerate(Data.subset)
    @printf("%i\t\t %1.0e \t\t %1.0e \t\t %1.0e \t\t %s\n", i, Data.I[i], Data.E[j], Data.angles[i], Data.names[j])
  end
  @printf("---------------------------------------------------------------------------------------------------\n")
end

function OSM!(
  Data :: OSMdata)

  ###############
  #  Load data  #
  ###############

  run(`figlet -w 800 Orthogonal Sensitivities`)

  #@printf("    ____       _   _                                   _    _____                _ _   _       _ _   _
   #/ __ \\     | | | |                                 | |  / ____|              (_) | (_)     (_) | (_)
  #| |  | |_ __| |_| |__   ___   __ _  ___  _ __   __ _| | | (___   ___ _ __  ___ _| |_ ___   ___| |_ _  ___  ___
  #| |  | | '__| __| '_ \\ / _ \\ / _` |/ _ \\| '_ \\ / _` | |  \\___ \\ / _ \\ '_ \\/ __| | __| \\ \\ / / | __| |/ _ \\/ __|
  #| |__| | |  | |_| | | | (_) | (_| | (_) | | | | (_| | |  ____) |  __/ | | \\__ \\ | |_| |\\ V /| | |_| |  __/\\__ \\
   #\\____/|_|   \\__|_| |_|\\___/ \\__, |\\___/|_| |_|\\__,_|_| |_____/ \\___|_| |_|___/_|\\__|_| \\_/ |_|\\__|_|\\___||___/
                               #__/ |
                              #|___/                                                                             \n")
  # Calculate E values
  OSM_calcE!(
    Data :: OSMdata)

  # If subset is empty, select the parameter with largest E-value
  if length(Data.subset) == 0
    first = sortperm(Data.E)[end]
    OSM_select_par!(
      first :: Int,
      Data :: OSMdata)
    end

  while length(Data.available)>0
    # As long as there are still parameters to add to the subset
    OSM_calculate_I_and_select!(
      Data :: OSMdata)
  end

  # Print results
  OSM_print_results(
    Data :: OSMdata)

  return int(0)
end

function Morris!(
  Data :: Morrisdata)

  """
  [Mus, S2] = morris_series_parallel(l,r,p,response)

  Calculates elementary effects using Morris method.

  Parameters
  ----------
  l : The number of points per dimension on the parameter grid.
  r : The number of evaluations of elementary effects for each parameter.
  p : The number of parameters.
  response : a function that takes a numpy.array of length p with parameter-values uniformly distributed between 0 and 1.

  Returns
  -------
  Mus : Sampling mean (the mean of the ABSOLUTE value of the elementary effects).
  S2 : Sampling variance.

  References
  ----------
  The code is presented on the algorithm presented by Smith in the book "Uncertainty Quantification: Theory, Implementation, and Applications" page 331-334.

  Comments
  --------
  The code could probably be improved by implementing methods for selecting good initial points on the parameter grid.

  Author
  ------
  Christian Haargaard Olsen <chaarga@ncsu.edu>

  """

  l = Data.level
  r = Data.points
  p = Data.no_parameters
  response = Data.f

  run(`figlet -w 800 Morris`)

  @printf("\n\nCalculating Morris measures mean and variance for l=%i, r=%i.\n",l,r)

  # For l even, this delta prevent Type2 errors
  delta = l/(2.0*(l-1))

  #Sample base parameter values
  q = Morris_qvalue(l,p)

  ###############################################################################
  #                       Do steps and calculate indices                        #
  ###############################################################################

  #Determine size of output
  resp = response(vec(q))
  n_t = length(resp)

  #Containers for responses and d values
  d   = zeros(r,p,n_t)
  uis = zeros(p,n_t)
  ui  = zeros(p,n_t)
  s2i = zeros(p,n_t)

  for i in 1:r
    # generate initial parameter configuration - the point to start in the grid
    q = Morris_qvalue(l,p)

    results = zeros(p+1,n_t)

    # Generate the permutation matrix - See Smith. 2013
    A = zeros(p+1,p)
    J = ones(p+1,p)
    for j in 1:size(A,1)
      for k in 1:size(A,2)
        if j>k
          A[j,k] = 1
        end
      end
    end

    D = diagm(floor(2*rand(p)-1)*2+1)
    P = eye(p)
    P = P[shuffle([1:p]),:]

    B = (J.*q' + delta/2*((2*A-J)*D + J))*P

    # In case the permutation takes any values to be larger than 1, find the indexes that overshoots, and subtract 2*delta instead
    over = -floor(-mod(B,1).*floor(B/1))
    B = B-2*delta*over

    # Find the order of steps, so the elementary effects can be assigned to the correct parameter
    shift = diff(B)
    shiftindex = [indmax(abs(shift[i,:])) for i in 1:p]
    deltas = sum(shift,1)

    #parallel code
    y0 = zeros(p+1,n_t)
    for j in 1:p+1
      y0[j,:] = response(vec(B[j,:]))
    end

    # Calculate the elementary effects, and assign them to the correct parameter
    for j in 1:p
      d[i,shiftindex[j],:] = (y0[j+1,:]-y0[j,:])/deltas[j]
    end
  end

  for i in 1:p
    for j in 1:n_t
      uis[i,j] = sum(abs(d[:,i,j]))/r
      ui[i,j] = sum(d[:,i,j])/r
      s2i[i,j] = sum((d[:,i,j]-ui[i,j]).^2)/(r-1)
    end
  end

  Data.mu = uis
  Data.s2 = sqrt(s2i)

  if Data.print
    Morris_print(Data)
  end

  return int32(0)
end

function Morris_print(
  Data :: Morrisdata)

  if Data.print
    munorm = [norm(Data.mu[i,:]) for i in 1:Data.no_parameters]
    s2norm = [norm(Data.s2[i,:]) for i in 1:Data.no_parameters]
    @printf("\n\nResults for Morris elementary effects analysis\n")
    @printf("---------------------------------------------------------------------------------------------------\n")
    @printf("Name \t \t μ \t \t σ2\n")
    @printf("---------------------------------------------------------------------------------------------------\n")
    for i in 1:Data.no_parameters
      @printf("%s\t\t %1.0e\t\t %1.0e.\n", Data.names[i], munorm[i], s2norm[i])
    end
    @printf("---------------------------------------------------------------------------------------------------\n")
  end
  @printf("See data structure fields .mu and .s2 for results.\n\n\n")
end

function Morris_qvalue(
  l :: Int,
  n=1 :: Int)
  # Draw parameter values from the proper grid points
  return floor(rand(n)*(l))/(l-1)
end

function Sobol!(
  Data :: Soboldata)

  """
  [Si, STi] = sobol_weirs(M,no_p,response)

  Calculates first order and total Sobol Indices of a function returning an array, using the method produced by Weirs et al., 2010.

  Parameters
  ----------
  M : The number of points in parameter space used.
  no_p: Number of parameters.
  response : A function accepting as input a numpy.array of parameter values.

  Returns
  -------
  Si : First order indices
  STi : Total indices

  References
  ----------
  The code is based on the method presented in the paper "Sensitivity analysis techniques applied to a system of hyperbolic conservation laws" by V. Gregory Weirs, James R. Kamm, Laura P. Swiler, Stefano Tarantola, Marco Ratto, Brian M. Adams, William J. Rider, Michael S. Eldred in Reliability Engineering and System Safety, vol. 107, 2012, pages 157-170.
  """

  M = Data.points
  run(`figlet -w 800 Sobol`)

  @printf("\n\nCalculating first order Sobol Indices Si and STi using M=%i points.\n", M)

  no_p     = Data.no_parameters
  response = Data.f

  # Draw random parameter values ~ U(0,1)
  A = rand(M,no_p)
  B = rand(M,no_p)

  # Asses the size of the output.
  A0   = copy(A)
  resp = response(vec(A0[1,:]))
  n_t  = length(resp)

  # Container for model output
  yA  = zeros(M,n_t)
  yB  = zeros(M,n_t)
  yAB = zeros(M,n_t)
  yBA = zeros(M,n_t)

  # Some parallel magic
  # from http://stackoverflow.com/questions/9742739/how-do-i-make-processes-able-to-write-in-an-array-of-the-main-program/9849971
  # pool = Pool(processes=cpu_count())

  # Populate rows of yA and yB with responses corresponding to parameters from rows of
  # A and B
  for i in 1:M
    yA[i,:] = response(vec(A[i,:]))
    yB[i,:] = response(vec(B[i,:]))
  end

  yC          = [yA,yB]
  yC2         = mean(yC,1).^2
  denominator = sum(yC.*yC,1)-2*M*yC2

  Data.Si     = zeros(no_p,n_t)
  Data.STi    = zeros(no_p,n_t)

  # For each parameter
  for i in 1:no_p
    Bi = copy(B)
    Bi[:,i] = copy(A[:,i])
    Ai = copy(A)
    Ai[:,i] = copy(B[:,i])

    # For parameter set
    for j=1:M
      yAB[j,:] = response(vec(Ai[j,:]))
      yBA[j,:] = response(vec(Bi[j,:]))
    end

    Data.Si[i,:]  = 2*sum(yA.*(yBA -yB),1)./denominator
    Data.STi[i,:] = sum((yA-yAB).*(yA-yAB),1)./denominator

    @printf("Done for parameter: %i/%i.\n", i, no_p)
  end

  Sobol_print(Data)
  return 0
end

function Sobol_print(
  Data :: Soboldata)

  if Data.print
    Sinorm = [norm(Data.Si[i,:]) for i in 1:Data.no_parameters]
    STinorm = [norm(Data.STi[i,:]) for i in 1:Data.no_parameters]
    @printf("\n\nResults for Sobol indices analysis\n")

    @printf("---------------------------------------------------------------------------------------------------\n")
    @printf("Name \t \t Si \t \t STi\n")
    @printf("---------------------------------------------------------------------------------------------------\n")
    for i in 1:Data.no_parameters
      @printf("%s\t\t %1.0e\t\t %1.0e.\n", Data.names[i], Sinorm[i], STinorm[i])
    end
    @printf("---------------------------------------------------------------------------------------------------\n")
  end
  @printf("See data structure fields .Si and .STi for results.\n\n\n")
end

function ishigami_function(
  x :: Array{Float64,1};
  a = 7.0 :: Float64,
  b = 0.05 :: Float64)
  """
  y = ishigami_function(x; a=7.0, b=0.05)

  The Ishigami function of Ishigami & Homma (1990) is used as an example for uncertainty
  and sensitivity analysis methods, because it exhibits strong nonlinearity and nonmonotonicity.
  It also has a peculiar dependence on x3, as described by Sobol' & Levitan (1999).
  (Description from http://www.sfu.ca/~ssurjano/ishigami.html)

  Parameters
  ----------
  x is a vector with three components with numbers between 0 and 1. They are scaled to -pi:pi before
  they are used in the function expression.
  a = 7.0 as suggested by Sobol & Levitan (1999)
  b = 0.05 as suggested by Sobol & Levitan (1999)

  Returns
  -------
  Returns function output for the function
      f(x) = sin(x[1]) + a*sin(x[2])^2 + b*x[3]^4*sin(x[1]).

  Comments
  --------
  Analytical values for a=7, b=0.05 for Si should be [0.222, 0.697, 0.000] (Smith, 2014.)

  References
  ----------
  The code is presented on the algorithm presented by Smith in the book "Uncertainty Quantification:
  Theory, Implementation, and Applications" page 331-334.

  Author
  ------
  Christian Haargaard Olsen <christianhaargaard@gmail.com>
  """

  # scale parameters to interval -pi:pi
  x = x*2.0*pi-pi
  y = (
    + sin(x[1])
    + a*(sin(x[2]))^2
    + b*x[3]^4*sin(x[1])
    )
  #return [y,y] # To test the tools for function returning more than one output eg. time-dependent
  return y
end

function sobols_function(
  x :: Array{Float64,1};
  a1 = 78.0:: Float64,
  a2 = 12.0 :: Float64,
  a3 = 0.50 :: Float64,
  a4 = 2.00 :: Float64,
  a5 = 97.0 :: Float64,
  a6 = 33.0 :: Float64,
  p = 6 :: Int,
  )
  """
  y = sobols_function(x; a=7.0, b=0.05)

  Parameters
  ----------
  x is a vector with three components with numbers between 0 and 1. They are scaled to -pi:pi before
  they are used in the function expression.
  a = 7.0 as suggested by Sobol & Levitan (1999)
  b = 0.05 as suggested by Sobol & Levitan (1999)

  Returns
  -------
  Returns function output for the function
      f(x) = sin(x[1]) + a*sin(x[2])^2 + b*x[3]^4*sin(x[1]).

  Comments
  --------
  Numerical values from Smith (2014) are
      |   x1    |   x2    |    x3   |    x4   |   x5    |   x6
    -------------------------------------------------------------
    μ | 5.60e-2 | 2.77e-1 |  1.76e0 |  1.19e0 | 3.50e-2 | 9.90e-2
    σ | 3.84e-2 | 1.90e-1 |  1.68e0 |  1.13e0 | 3.74e-2 | 2.38e-2
    -------------------------------------------------------------

  References
  ----------
  The code is presented on the algorithm presented by Smith in the book "Uncertainty Quantification:
  Theory, Implementation, and Applications" page 331-334.

  Author
  ------
  Christian Haargaard Olsen <christianhaargaard@gmail.com>
  """

  # scale parameters to interval -pi:pi
  y = (
    (abs(4*x[1]-2.0)+a1)/(1.0+a1)*
    (abs(4*x[2]-2.0)+a2)/(1.0+a2)*
    (abs(4*x[3]-2.0)+a3)/(1.0+a3)*
    (abs(4*x[4]-2.0)+a4)/(1.0+a4)*
    (abs(4*x[5]-2.0)+a5)/(1.0+a5)*
    (abs(4*x[6]-2.0)+a6)/(1.0+a6)
    )
  #return [y,y] # To test the tools for function returning more than one output eg. time-dependent
  return y
end
end


text:
## Model definition

""" `Model` objects assume that there is a pure-kwarg constructors (such as provided by
Parameters.jl) """
abstract type Model end

# Models should implement these functions (or rely on the defaults)
transition_mat(m::Model, inputs, i) = Identity()
transition_noise(m::Model, inputs, i) = no_noise()
observation_mat(m::Model, inputs, i) = Identity()
function observation_noise end  # necessary for all models
function observation end        # optional; you can pass `observations=...` instead
function initial_state end      # optional; you can pass `initial_state=...` instead

# More exotic model functions
skip_observation(m::Model, inputs, i) = false
nobservations(m::Model, inputs) = size(observation_mat(m, inputs, 1), 1) #fails if skipped
function labels end             # for pretty-printing
parameters(m::M) where M<:Model = fieldnames(M)

################################################################################

predicted_state(state_prior::Gaussian, transition_mat, transition_noise::Gaussian) =
    Gaussian(transition_mat * mean(state_prior) + mean(transition_noise),
             transition_mat * cov(state_prior) * transition_mat' + cov(transition_noise))

predicted_state(prev_state::Gaussian, m::MiniKalman.Model, inputs, t::Int) =
    predicted_state(prev_state, transition_mat(m, inputs, t),
                    transition_noise(m, inputs, t))

function predicted_observation(state::Gaussian, observation_mat, observation_noise)
    # Prediction step (μ_(t|t-1), Σ_(t|t-1))
    C = observation_mat
    Du, R = observation_noise     # Du := Dₜuₜ
    μ, Σ = state
    ŷ = C * μ + Du               # Murphy forgot the Du term in his book
    S = C * Σ * C' + R
    return Gaussian(ŷ, S)
end

predicted_observation(state::Gaussian, m::MiniKalman.Model, inputs, t::Int) =
    predicted_observation(state, observation_mat(m, inputs, t),
                          observation_noise(m, inputs, t))

function kfilter_obs(state::Gaussian,
                     observation, observation_mat, observation_noise::Gaussian)
    # Following Machine Learning, A probabilistic perspective, by Kevin Murphy, 18.3.1.2,
    # with the exception that I'm setting the forcing term to 0, but allowing noise terms
    # with means (without loss of generality)

    C = observation_mat
    y = observation
    μ, Σ = state

    ŷ, S = predicted_obs =
        predicted_observation(state, observation_mat, observation_noise)

    # Filter
    K = Σ * C' / S         # Kalman gain matrix
    r = y - ŷ
    filtered_state = Gaussian(μ + K*r,
                              (I - K*C) * Σ)
    ll = logpdf(predicted_obs, y)  # log-likelihood
    return (filtered_state, ll, predicted_obs)
end

""" Perform one step of Kalman filtering, for online use. We assume equations:

```julia
    current_state = transition_mat * previous_state + transition_noise
    observation = observation_mat * current_state + observation_noise
```

and return this tuple:

```
   (current_state::Gaussian,   # the posterior `P(state|observation)`
    ll,                        # the log-likelihood for that step
    predicted_obs::Gaussian)   # `P(observation|state)`
```

To add an "input" term, pass it as `transition_noise = Gaussian(input_term, noise_cov)`.
"""
function kfilter(prev_state::Gaussian, m::MiniKalman.Model, inputs, t::Int,
                 observations=nothing)
    state = predicted_state(prev_state, m, inputs, t)
    if skip_observation(m, inputs, t)
        # Optimization note: returning ll=0.0 might make this function type-unstable when
        # computing gradients, for gradient descent.
        # We can't have `kfilter_obs` run when `skip_observation` is true, because
        # bad matrices (singular / containing missing) can trigger exceptions.
        (state, 0.0, predicted_observation(state, m, inputs, t))
    else
        kfilter_obs(state,
                    observations===nothing ? observation(inputs, t) : observations[t],
                    observation_mat(m, inputs, t), observation_noise(m, inputs, t))
    end
end

function prediction_type(m::Model, inputs, any_state::Gaussian)
    # This juggling is necessary to get the proper `SVector/Vector` output
    # type. It's lame, but I don't see how else to do it (except maybe using `accumulate`,
    # but IIRC there were some performance concerns)
    fake_obs_mat = hcat((mean(any_state) for _ in 1:nobservations(m, inputs))...)'
    fake_obs = fake_obs_mat * mean(any_state)
    return Gaussian{typeof(fake_obs), typeof(fake_obs*fake_obs')}
end

""" (Internal function) return three vectors, appropriate for storing
states, likelihoods, and P(observation). """
function output_vectors(m::Model, inputs, observations=nothing; length=length(inputs))
    # Note: this function was split off for our own internal purposes (switching states)
    # It could be merged back into `kalman_filter!`, I suppose.
    state = make_full(initial_state(m))
    filtered_states = Vector{typeof(state)}(undef, length)
    predicted_obs = Vector{prediction_type(m, inputs, state)}(undef, length)
    lls = Vector{Float64}(undef, length)
    return (filtered_states, lls, predicted_obs)
end

""" Perform Kalman filtering and store the results in the `filtered_states, lls,
predicted_obs` vectors, for the given `steps`. """
function kalman_filter!(filtered_states::AbstractVector, lls::AbstractVector,
                        predicted_obs::AbstractVector,
                        m::Model, inputs, observations=nothing;
                        steps::AbstractRange=1:length(filtered_states),
                        initial_state=(steps[1]==1 ? initial_state(m) :
                                       filtered_states[steps[1]-1]))
    state = make_full(initial_state)  # we need make_full so that the state does
                                      # not change type during iteration
    for t in steps
        state, lls[t], predicted_obs[t] = kfilter(state, m, inputs, t, observations)
        filtered_states[t] = state
    end
end

function kalman_filter(m::Model, inputs, observations=nothing;
                       initial_state=initial_state(m),
                       steps=1:length(observations===nothing ? inputs : observations))
    out_vecs = output_vectors(m, inputs, observations, length=length(steps))
    kalman_filter!(out_vecs..., m, inputs, observations; steps=steps,
                   initial_state=initial_state)
    return out_vecs
end

function log_likelihood(m::Model, inputs, observations=nothing;
                        initial_state=MiniKalman.initial_state(m),
                        steps=1:length(observations===nothing ? inputs : observations))
    # Since this is in the inner loop of `optimize`, we make sure it's non-allocating
    ll_sum = 0.0
    state = make_full(initial_state)
    for t in steps
        state, ll, _ = kfilter(state, m, inputs, t, observations)
        ll_sum += ll
    end
    return ll_sum
end

kalman_filtered(args...; kwargs...) = kalman_filter(args...; kwargs...)[1]  # convenience
kalman_predicted(args...; kwargs...) = kalman_filter(args...; kwargs...)[3]  # convenience

""" Convenience; returns a vector of the total likelihood up to each step. """
cumulative_log_likelihood(args...; kwargs...) =
    cumsum(kalman_filter(args...; kwargs...)[2])

ext: jl


text:
"""Calculate the number of grains on square `square`."""
function on_square(square)
    square < 1 && throw(DomainError("Square cannot be negative"))
    square > 64 && throw(DomainError("There are only 64 squares on a chessboard"))
    _on_square(square)
end

# "Follow the question" implementation
function _on_square(square)
    square == 1 && return Int128(1)
    2*on_square(square-1)
end

"""Calculate the total number of grains after square `square`."""
function total_after(square)
    square < 1 && throw(DomainError("Square cannot be negative"))
    sum(on_square.(1:square))
end

# "Let's do some maths" implementation
function _on_square_quick(square)
    Int128(2)^(square-1)
end

# I had two options:
# 1. Try to remember how sums of geometric series work
# 2. use SymPy
# ```
# using SymPy
# @vars(k, n)
# summation(2^(k-1),(k,(1,n)))
# ```
function total_after_quick(square)
    Int128(2)^(square+1)÷2-1
end


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

using Test

Test.@testset "Package BAT" begin
    include("utils/test_utils.jl")
    include("rngs/test_rngs.jl")
    include("distributions/test_distributions.jl")
    include("parameters/test_parameters.jl")
    include("statistics/test_statistics.jl")
    include("densities/test_densities.jl")
    include("samplers/test_samplers.jl")
    include("io/test_io.jl")
    include("plotting/test_plotting.jl")
end


text:
module ESValues

using GLM
using StatsBase

# package code goes here
include("iterators.jl")

export esvalues, ESValuesEstimator

typealias MaskType Float64

type ESValuesEstimator{T}
    x
    f::Function
    X
    link
    featureGroups::Vector{Vector{Int64}}
    weights
    nsamples::Int64
    varyingInds::Vector{Int64}
    varyingFeatureGroups::Vector{Vector{Int64}}
    data::Matrix{T}
    maskMatrix::Matrix{MaskType}
    kernelWeights::Vector{MaskType}
    y::Vector{Float64}
    ey::Vector{Float64}
    lastMask::Vector{Float64}
    P::Int64
    N::Int64
    M::Int64
    nsamplesAdded::Int64
    nsamplesRun::Int64
    fx::Float64
    fnull::Float64
end

"Designed to determine the ES values (importance) of each feature for f(x)."
function esvalues(x, f::Function, X, link=identity; featureGroups=nothing, weights=nothing, nsamples=0)
    esvalues(ESValuesEstimator(f, X, link; featureGroups=featureGroups, weights=weights, nsamples=nsamples), x)
end

function esvalues(e::ESValuesEstimator, x)

    @assert length(x) == e.P "Provided 'x' length must match the data matrix features count ($(length(x)) != $(e.P))!"
    e.x = x

    # find the feature groups we will test. If a feature does not change from its
    # current value then we know it doesn't impact the model
    e.varyingInds = varying_groups(e.x, e.X, e.featureGroups)
    e.varyingFeatureGroups = e.featureGroups[e.varyingInds]
    e.M = length(e.varyingFeatureGroups)

    # find f(x) and E_x[f(x)]
    e.fx = e.f(x)[1]
    e.fnull = sum(vec(e.f(e.X)) .* e.weights)

    # if no features vary then there no feature has an effect
    if e.M == 0
        return e.fx,zeros(length(e.featureGroups)),zeros(length(e.featureGroups))

    # if only one feature varies then it has all the effect
    elseif e.M == 1
        fx = mean(e.f(x))
        fnull = sum(vec(e.f(e.X)) .* e.weights)
        φ = zeros(length(e.featureGroups))
        φ[e.varyingInds[1]] = e.link(e.fx) - e.link(e.fnull)
        return e.fnull,φ,zeros(length(e.featureGroups))
    end

    # pick a reasonable number of samples if the user didn't specify how many they wanted
    if e.nsamples == 0
        e.nsamples = 2*e.M+1000
    end
    if e.M <= 30 && e.nsamples > 2^e.M-2
        e.nsamples = 2^e.M-2
    end
    @assert e.nsamples >= min(2*e.M, 2^e.M-2) "'nsamples' must be at least 2 times the number of varying feature groups!"

    # add the singleton samples
    allocate!(e)
    for (m,w) in take(drop(eskernelsubsets(collect(1:e.M), ones(e.M)), 2), 2*e.M)
        addsample!(e, x, m, w)
    end
    run!(e)

    # if there might be more samples then enumarate them
    if length(e.y) >= 2*e.M

        # estimate the variance of each ES value estimate
        variances = zeros(e.M)
        for i in 1:2:2*e.M
            variances[div(i+1,2)] = var([e.y[i] - e.fnull, e.fx - e.y[i+1]])
        end

        # now add the rest of the samples giving priority to ES values with high estimated variance
        for (m,w) in take(drop(eskernelsubsets(collect(1:e.M), variances), 2*e.M+2), e.nsamples-(2*e.M))
            addsample!(e, x, m, w)
        end
        run!(e)
    end

    # solve then expand the ES values vector to contain the non-varying features as well
    vφ,vφVar = solve!(e)
    φ = zeros(length(e.featureGroups))
    φ[e.varyingInds] = vφ
    φVar = zeros(length(e.featureGroups))
    φVar[e.varyingInds] = vφVar

    # return the Shapley values along with variances of the estimates
    e.fnull,φ,φVar
end

function ESValuesEstimator{T}(f::Function, X::Matrix{T}, link=identity; featureGroups=nothing, weights=nothing, nsamples=0)
    P,N = size(X)

    # give default values to omitted arguments
    weights != nothing || (weights = ones(N))
    weights ./= sum(weights)
    featureGroups != nothing || (featureGroups = Array{Int64,1}[Int64[i] for i in 1:size(X)[1]])
    featureGroups = convert(Array{Array{Int64,1},1}, featureGroups)
    @assert length(weights) == N "Provided 'weights' must match the number of representative data points (size(X)[2])!"

    ESValuesEstimator(
        zeros(1),
        f,
        X,
        link,
        featureGroups,
        weights,
        nsamples,
        Int64[],
        Vector{Int64}[],
        zeros(T, 1, 1),
        zeros(MaskType, 1, 1),
        zeros(MaskType, 1),
        zeros(Float64, 1),
        zeros(Float64, 1),
        zeros(Float64, 1),
        P,
        N,
        0,
        0,
        0,
        0.0,
        0.0
    )
end

function allocate!{T}(e::ESValuesEstimator{T})
    e.data = zeros(T, e.P, e.nsamples * e.N)
    e.maskMatrix = zeros(MaskType, e.M-1, e.nsamples)
    e.kernelWeights = zeros(MaskType, e.nsamples)
    e.y = zeros(Float64, e.nsamples * e.N)
    e.ey = zeros(Float64, e.nsamples)
    e.lastMask = zeros(Float64, e.nsamples)
end

function addsample!(e::ESValuesEstimator, x, m, w)
    offset = e.nsamplesAdded * e.N
    e.nsamplesAdded += 1
    for i in 1:e.N
        for j in 1:e.M
            for k in e.varyingFeatureGroups[j]
                if m[j] == 1.0
                    e.data[k,offset+i] = x[k]
                else
                    e.data[k,offset+i] = e.X[k,i]
                end
            end
        end
    end
    e.maskMatrix[:,e.nsamplesAdded] = m[1:end-1] - m[end]
    e.lastMask[e.nsamplesAdded] = m[end]
    e.kernelWeights[e.nsamplesAdded] = w
end

function run!(e::ESValuesEstimator)
    e.y[e.nsamplesRun*e.N+1:e.nsamplesAdded*e.N] = e.f(e.data[:,e.nsamplesRun*e.N+1:e.nsamplesAdded*e.N])

    # find the expected value of each output
    for i in e.nsamplesRun+1:e.nsamplesAdded
        eyVal = 0.0
        for j in 1:e.N
            eyVal += e.y[(i-1)*e.N + j] * e.weights[j]
        end
        e.ey[i] = eyVal
        e.nsamplesRun += 1
    end
end

function solve!(e::ESValuesEstimator)

    # adjust the y value according to the constraints for the offset and sum
    eyAdj = e.link.(e.ey) .- e.lastMask*(e.link(e.fx) - e.link(e.fnull)) - e.link(e.fnull)

    # solve a weighted least squares equation to estimate φ
    tmp = e.maskMatrix .* e.kernelWeights'
    tmp2 = inv(tmp*e.maskMatrix')
    w = tmp2*(tmp*eyAdj)
    wlast = (e.link(e.fx) - e.link(e.fnull)) - sum(w)
    φ = [w; wlast]

    yHat = e.maskMatrix'w
    φVar = var(yHat .- eyAdj) * diag(tmp2)
    φVar = [φVar; maximum(φVar)] # since the last weight is inferred we use a pessimistic guess of its variance

    # a finite sample adjustment based on how much of the weight is left in the sample space
    fractionWeightLeft = 1 - sum(e.kernelWeights)/sum([(e.M-1)/(s*(e.M-s)) for s in 1:e.M-1])

    φ,φVar*fractionWeightLeft
end

"Identify which feature groups vary."
function varying_groups(x, X, featureGroups)
    varying = zeros(length(featureGroups))
    for (i,inds) in enumerate(featureGroups)
        varying[i] = sum(vec(sum(x[inds] .== X[inds,:],1) .!= length(inds)))
    end
    find(varying)
end

end # module


text:
struct PoincareInvariant1st{ET,DT,TT,ΘT}
    equ::ET
    Θ::ΘT
    Δt::TT
    nloop::Int
    ntime::Int
    nsave::Int
    nt::Int
    I::Vector{DT}
    J::Vector{DT}
    L::Vector{DT}
end

function PoincareInvariant1st(f_equ::Function, f_loop::Function, Θ::ΘT, Δt::TT, d::Int, nloop::Int, ntime::Int, nsave::Int=1, DT=Float64) where {TT,ΘT}

    if get_config(:verbosity) > 1
        println()
        println("First Euler-Poincaré Integral Invariant")
        println("=======================================")
        println()
        println(" nloop = ", nloop)
        println(" ntime = ", ntime)
        println(" nsave = ", nsave)
        println(" Δt    = ", Δt)
        println()
    end

    # compute initial conditions
    q₀ = zeros(DT, (d, nloop))

    for i in 1:nloop
        q₀[:,i] .= f_loop(i/nloop)
    end

    # initialise euation
    equ = f_equ(q₀)

    # create arrays for results
    nt = div(ntime, nsave)

    I = zeros(DT, nt+1)
    J = zeros(DT, nt+1)
    L = zeros(DT, nt+1)

    PoincareInvariant1st{typeof(equ),DT,TT,ΘT}(equ, Θ, Δt, nloop, ntime, nsave, nt, I, J, L)
end

function evaluate_poincare_invariant(pinv::PoincareInvariant1st, sol::Solution)
    p = zeros(sol.q.d)
    g = zeros(sol.q.d)
    v = zeros(size(sol.q.d,1), size(sol.q.d,3))
    γ = zeros(size(sol.q.d,1), size(sol.q.d,3))

    compute_one_form(sol.t.t, sol.q.d, p, pinv.Θ)

    if isdefined(sol, :λ)
        compute_correction(sol.t.t, sol.q.d, sol.λ.d, g, pinv.equ.g)
    end

    for i in 1:size(sol.q.d,2)
        compute_velocity(sol.q.d[:,i,:], v)
        pinv.I[i] = compute_loop_integral(p[:,i,:], v)

        if isdefined(sol, :p)
            pinv.J[i] = compute_loop_integral(sol.p.d[:,i,:], v)
        end

        if isdefined(sol, :λ)
            compute_velocity(sol.λ.d[:,i,:], γ)
            pinv.L[i] = compute_loop_integral(p[:,i,:] .- pinv.Δt .* g[:,i,:], v .- pinv.Δt .* γ)
        end
    end

    (pinv.I, pinv.J, pinv.L)
end

function CommonFunctions.write_to_hdf5(pinv::PoincareInvariant1st, sol::Solution, output_file::String)
    # h5open(output_file, isfile(output_file) ? "r+" : "w") do h5
    h5open(output_file, "w") do h5

        write(h5, "t", sol.t.t)
        write(h5, "I", pinv.I)

        isdefined(sol, :p) ? write(h5, "J", pinv.J) : nothing
        isdefined(sol, :λ) ? write(h5, "L", pinv.L) : nothing

    end
end


text:
"""
    SLab.Machine

Machine definition.
"""
mutable struct Machine

    """machine id"""
    id::Int64

    """machine type"""
    T::Int64

    """machine name"""
    name::String

end

"""
    SLab.Operation

Operation definition.
"""
mutable struct Operation

    """parent job id"""
    job_id::Int64

    """operation id"""
    id::Int64

    """compatible machine type"""
    C::Int64

    """process time in sec."""
    tau::Int64

    """start time"""
    S::Int64

    """machine id to process this operation"""
    E::Int64

    """note"""
    note::String

    function Operation(job_id::Int64, id::Int64, C::Int64, tau::Int64, note::String)
        new(job_id, id, C, tau, 0, 0, note)
    end

end

"""get end time of an operation"""
function end_time(op::Operation)::Int64
    if op.S == 0
        return 0
    else
        return op.S + op.tau
    end
end

"""
    SLab.TCMB

Time constraint by mutual boundaries
"""
mutable struct TCMB

    """operation id 1"""
    operation_id_1::Int64

    """boundary of operation 1 (begin or end)"""
    boundary_1::Symbol

    """operation id 2"""
    operation_id_2::Int64

    """boundary of operation 2 (begin or end)"""
    boundary_2::Symbol

    """maximum difference in time (sec.) between two boundaries"""
    alpha::Int64

end

"""
    SLab.Job

Job definition, containing a series of operations, dependency graph among operations, and TCMB.
"""
mutable struct Job

    """job id"""
    id::Int64

    """list of operations"""
    operations::Vector{Operation}

    """operation dependency graph"""
    dag::Vector{Tuple{Int64,Int64}}

    """time constraint by mutual boundaries"""
    constraints::Vector{TCMB}

end

"""
    SLab.Batch

To schedule multiple jobs simultaneously, combine those jobs into a single Batch.
"""
mutable struct Batch

    """list of contained jobs"""
    jobs::Vector{Job}

    """all operations in the jobs in a single vector"""
    operations::Vector{Operation}

    """conversion table from (job_id, op_id_in_job) to op_id_in_batch"""
    operation_indices::Dict{Tuple{Int64,Int64},Int64}

    function Batch(jobs::Vector{Job})
        operations = Vector{Operation}()
        operation_indices = Dict{Tuple{Int64,Int64},Int64}()
        for job in jobs
            for op in job.operations
                push!(operations, op)
                operation_indices[job.id, op.id] = length(operations)
            end
        end
        new(jobs, operations, operation_indices)
    end

end

"""get number of operations in a batch"""
function get_N(batch::Batch)::Int64
    return length(batch.operations)
end

"""get operation dependency graph in a batch"""
function get_P(batch::Batch)::Matrix{Bool}
    P = zeros(Bool, get_N(batch), get_N(batch))
    for job in batch.jobs
        for edge in job.dag
            P[batch.operation_indices[job.id, edge[1]], batch.operation_indices[job.id, edge[2]]] = 1
        end
    end
    return P
end

"""get compatible machine types of operations in a batch"""
function get_C(batch::Batch)::Vector{Int64}
    return [op.C for op in batch.operations]
end

"""get processing times for operations in a batch"""
function get_tau(batch::Batch)::Vector{Int64}
    return [op.tau for op in batch.operations]
end

"""get TCMB in a batch"""
function get_constraints(batch::Batch)::Vector{TCMB}
    constraints = Vector{TCMB}()
    for job in batch.jobs
        for con in job.constraints
            push!(
                constraints,
                TCMB(
                    batch.operation_indices[job.id, con.operation_id_1],
                    con.boundary_1,
                    batch.operation_indices[job.id, con.operation_id_2],
                    con.boundary_2,
                    con.alpha,
                ),
            )
        end
    end
    return constraints
end

"""set start time and machine id for each operation in a batch"""
function set_schedule(batch::Batch, S::Vector{Int64}, E::Vector{Int64})
    for a = 1:get_N(batch)
        batch.operations[a].S = S[a]
        batch.operations[a].E = E[a]
    end
end


text:
using ITensors
using ITensors.NDTensors
using Test

@testset "ITensor (Empty)" begin
  @testset "ITensor set elements" begin
    i = Index(2; tags="i")

    E = ITensor(i', dag(i))

    @test conj(E) == E
    @test 1.2 * E == E

    @test hassameinds(E, (i', i))
    @test order(E) == 2
    @test E[i' => 1, i => 1] == 0

    E[i' => 1, i => 2] = 2.3

    @test E[i' => 1, i => 1] == 0
    @test E[i' => 2, i => 1] == 0
    @test E[i' => 1, i => 2] == 2.3
    @test E[i' => 2, i => 2] == 0
  end

  @testset "ITensor (Empty) convert to complex" begin
    i = Index(2; tags="i")
    E = ITensor(i', dag(i))
    @test eltype(E) == NDTensors.EmptyNumber

    Ec = complex(E)
    @test eltype(Ec) == Complex{NDTensors.EmptyNumber}
    Ec[1, 1] = 2.3
    @test eltype(Ec) == ComplexF64

    Ec = complex(E)
    @test eltype(Ec) == Complex{NDTensors.EmptyNumber}
    Ec[1, 1] = 2.3f0
    @test eltype(Ec) == ComplexF32

    E2 = copy(E)
    E2c = complex!(E2)
    @test eltype(E2c) == Complex{NDTensors.EmptyNumber}
  end

  @testset "ITensor set elements (QN)" begin
    i = Index(QN(0) => 2, QN(1) => 2; tags="i")

    E = ITensor(i', dag(i))

    @test hassameinds(E, (i', i))
    @test order(E) == 2
    @test isnothing(flux(E))
    @test E[i' => 1, i => 3] == 0

    E[i' => 3, i => 2] = 2.3

    @test flux(E) == QN(1)

    @test E[i' => 1, i => 1] == 0
    @test E[i' => 2, i => 1] == 0
    @test E[i' => 3, i => 2] == 2.3
    @test E[i' => 2, i => 3] == 0
    @test_throws ErrorException E[i' => 2, i => 3] = 3.2
  end

  @testset "ITensor()" begin
    i = Index(QN(0) => 2, QN(1) => 2; tags="i")

    E = ITensor()

    @test isnothing(flux(E))
    @test order(E) == 0
    @test_throws MethodError E[i' => 1, i => 3] = 0

    A = randomITensor(i', dag(i))
    E += A

    @test norm(E - A) < 1E-8
  end
end

nothing


text:
# This file was generated, do not modify it. # hide
knnc.K = 3
fit!(classif, rows=train)
ŷ = predict_mode(classif, rows=test)
accuracy(ŷ, y[test]) |> r3

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
function SAG(X, y, x0, f, g, eps, λ, maxiter, stepsize, num_data)

    x_k = x0
    x_old = x_k
    oval = Inf
    fval = Inf
    # Arrays of grad_inf and function_values
    val_array = zeros(0)

    # M is my matrix of derivatives (column i = f_i'(x))
    M = zeros(size(x0,1),num_data)
    sum = zeros(size(x0,1))
    for i in 1:num_data
        # compute f' for each data_point at x0
        x_i = X[i:i,:]'
        gval_idx = g( x_i ,y[i] ,x0 ,λ/num_data )
        # fill i-th column of M
        M[:,i:i] = gval_idx
        sum += gval_idx
    end
    sum = sum / num_data

    @printf("Initial derivatives are initialized.\n")
    @printf("Stepsize was choosen to α = %f.\n", stepsize )
    @printf("iter \t\t j-pick \t\t fval \t\t\t diff \t\t opt\n")

    @time for i in 1:maxiter

        # Step 1:
        # pick a j uniformly at random. j ∈ [1,n]
        j = rand(1:num_data)
        x_j = X[j:j,:]'
        old_gval_j = M[:,j:j]

        # Step 2:
        # Take phi_j_k_plus = xk and store f'_j(phi_j_k_plus) in the table
        phi_j_k_plus = x_k
        gval_j = g(x_j ,y[j] , phi_j_k_plus, λ / num_data )
        stoppingCriteria = norm(gval_j,Inf)
        # update gval_j in table
        M[:,j:j] = gval_j

        # Step 3:
        # Update x using f'_j(phi_j_k_plus),
        x_old = x_k
        # update the sum
        sum = sum - old_gval_j / num_data
        sum = sum + gval_j / num_data

        # two different updates
        #x_k = x_k - stepsize * ( (gval_j-old_gval_j)/num_data +  sum )
        x_k = x_k - stepsize * sum

        oval = fval
        fval = f(x_k)
        append!( val_array, fval )
        @printf("%i \t\t %i \t\t %f \t\t %f \t\t %e \n",i,j,fval,(oval-fval),stoppingCriteria)
        end
    return (x_k, "iter", maxiter, val_array)
end

max_stars_repo_name: zagazao/num-opt

max_issues_repo_name: zagazao/num-opt

max_forks_repo_name: zagazao/num-opt


text:
abstract type AbstractDemand end
abstract type AbstractGeneration end
abstract type AbstractStorage end
abstract type AbstractConverter end
abstract type AbstractGrid end
abstract type AbstractDesigner end
abstract type AbstractController end

mutable struct GlobalParameters
    ns::Int64 # number of [nh, ny] scenarios
    Δh::Int64 # operations time step in hours
    nh::Int64 # number of operation stages
    Δy::Int64 # investment time step in years
    ny::Int64 # number of investment stages
    τ::Float64 # discount rate
    renewable_share::Float64 # share of renewables [0,1]

    GlobalParameters(nh, ny, ns;
                Δh = 1,
                Δy = 1,
                τ = 0.045,
                renewable_share = 0.) = renewable_share == 1. ? new(ns, Δh, nh, Δy, ny, τ, 0.9999) : new(ns, Δh, nh, Δy, ny, τ, renewable_share)
end

mutable struct Microgrid
    parameters::GlobalParameters
    demands::Vector{AbstractDemand}
    generations::Vector{AbstractGeneration}
    storages::Vector{AbstractStorage}
    converters::Vector{AbstractConverter}
    grids::Vector{AbstractGrid}

    Microgrid(; parameters = GlobalParameters(8760, 20, 1)) = new(parameters)
end

# Add node to microgrid
# TODO : add asset by asset without initializing everything...
function add!(mg::Microgrid, assets...)
    # Build and preallocate
    mg.demands = [preallocate!(a, mg.parameters.nh, mg.parameters.ny, mg.parameters.ns) for a in assets if a isa AbstractDemand]
    mg.generations = [preallocate!(a, mg.parameters.nh, mg.parameters.ny, mg.parameters.ns) for a in assets if a isa AbstractGeneration]
    mg.storages = [preallocate!(a, mg.parameters.nh, mg.parameters.ny, mg.parameters.ns) for a in assets if a isa AbstractStorage]
    mg.converters = [preallocate!(a, mg.parameters.nh, mg.parameters.ny, mg.parameters.ns) for a in assets if a isa AbstractConverter]
    mg.grids = [preallocate!(a, mg.parameters.nh, mg.parameters.ny, mg.parameters.ns) for a in assets if a isa AbstractGrid]
end

# Preallocate abstract designer
function preallocate!(mg::Microgrid, designer::AbstractDesigner)
    designer.decisions = (generations = [zeros(mg.parameters.ny, mg.parameters.ns) for a in mg.generations],
                          storages = [zeros(mg.parameters.ny, mg.parameters.ns) for a in mg.storages],
                          converters = [zeros(mg.parameters.ny, mg.parameters.ns) for a in mg.converters])
end
# Preallocate abstract controller
function preallocate!(mg::Microgrid, controller::AbstractController)
    controller.decisions = (converters = [zeros(mg.parameters.nh, mg.parameters.ny, mg.parameters.ns) for a in mg.converters],
                            storages = [zeros(mg.parameters.nh, mg.parameters.ny, mg.parameters.ns) for a in mg.storages])
end

# Find if the datatype is in a mg field
function isin(field::Vector, type::DataType)
    # Return true if the datatype is in the field and its index
    bool, idx = false, NaN
    for (k, a) in enumerate(field)
        if a isa type
            bool, idx = true, k
        end
    end
    return bool, idx
end


text:
include("nanop.jl")

using NCDatasets
using Formatting
using Statistics

Re  = 6371000.0 #  m
c_p  = 1004.0   #  J / kg K
g   = 9.81      #  m / s^2
Lw  = 2.26e6    #  J / kg

linestyle = ["r-", "g-", "b-"]

Dataset("domain_nc_files/domain.lnd.fv4x5_gx3v7.091218.nc", "r") do ds
    global lat, lon, lev, ilev, rec

    lat = replace(ds["yc"][1,:], missing=>NaN)
    lon = replace(ds["xc"][:,1], missing=>NaN)

end

using PyPlot
plt[:rcParams]["font.size"] = 40
plt[:rcParams]["axes.label.size"] = 40

fig, ax = plt[:subplots](1, 1, figsize=(12,8))
#ax[:plot]([-90, 90], [0, 0], linewidth=2, color="#cccccc")

casename = "lowres_SSM_NK"

Dataset("extract_nc_files/$casename.ocn.h.ma.fv4x5.0001-0020.nc", "r") do ds

    t_beg = 240 - 10*12 + 1
    t_end = 240

    cnt_JJA = 0
    cnt_DJF = 0

    mld_JJA = zeros(length(lon), length(lat))
    mld_DJF = zeros(length(lon), length(lat))

    for t = t_beg:t_end
        m = mod(t - 1, 12) + 1

        if m in [11, 12, 1]
            mld_DJF += replace(ds["mld"][:, :, t], missing=>NaN)
            cnt_DJF += 1
        elseif m in [6, 7, 8]
            mld_JJA += replace(ds["mld"][:, :, t], missing=>NaN)
            cnt_JJA += 1
        end

    end

    mld_JJA ./= cnt_JJA
    mld_DJF ./= cnt_DJF

    mld_JJA = nanmean(mld_JJA, dims=(1,))[1, :]
    mld_DJF = nanmean(mld_DJF, dims=(1,))[1, :]

    ax[:plot](lat, mld_DJF, "k-", linewidth=2,  label="DJF")
    ax[:plot](lat, mld_JJA, "k--", linewidth=2,  label="JJA")

end

ax[:legend](loc="lower center", fontsize=15)
ax[:grid](linestyle="--")
ax[:set_xlabel]("Latitude [deg]", size=15)
ax[:set_ylabel]("Mixed-layer Depth [m]", size=15)
ax[:set_xlim](-90, 90)
ax[:set_ylim](-5, 800)

ax[:invert_yaxis]()

ax[:set_title]("lowres_SSM_NK 11-20 year MLD average in DJF(solid) and JJA(dashed).")

fig[:savefig]("$(casename)_MLD.png", dpi=200)

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
mutable struct AlexyaApp
    title::String # The window title
    width::Int # Canvas width
    height::Int # Canvas height
    framerate::Real # FPS
    layout::Function # Layout manager
    window::GtkWindow # The GtkWindow widget
    canvas::GtkCanvas # The GtkCanvas widget
    events::Dict{Symbol, Function}
    loop::Bool
    framecount::Int
    fps::Real # Current FPS
    widgets::Vector{<:GtkWidget}
    drawing::Drawing # Luxor Drawing object

    AlexyaApp() = new()
end

function destroy(app::AlexyaApp)
    Gtk.destroy.([app.widgets..., app.canvas, app.window])
end

const CURRENT_ALEXYA_APP = Array{AlexyaApp, 1}()

function current_app()
    try
        return CURRENT_ALEXYA_APP[1]
    catch
        throw(error("There is no current Alexya Application."))
    end
end

current_app(field::Symbol) = getfield(current_app(), field)

"""
        addwidget(widget::GtkWidget)

Adds the widget `widget` in the current window.

# Examples

```julia
julia> slider = Slider(1:100)
GtkSlider...

julia> addwidget(slider)
GtkSlider...
```
"""
function addwidget(widget::GtkWidget)
    widgets = current_app(:widgets)
    push!(widgets, widget)
    return widget
end

"""
        init(title [, width = 400, height = 400])

Initialize `Alexya` by creating a new `AlexyaApp` global object.

The `title` parameter is the title of the `window`, the `width`
and `height` parameters are not necessarily the same as the
window, but they determine the size of the canvas. The size
of the window and canvas can be controlled by the current
`layout`, and the behavior of the `addwidget` function or
`@create` macro are also controlled by the `layout`.

The default `layout` is `canvasonly`, wich means that
the canvas will be the only widget on the window and
uses of the `addwidget` function or `@create` macro
will not work.

You can set the current `layout` by calling the `uselayout` function
or the `@layout` macro. The `layout` is a just a function that
determines how the widgets (including the canvas) are added to the
window. You can also create your own layouts.

# Examples

```julia
julia> init("My beautiful animation", 800, 600)
AlexyaApp...

julia> @layout aside(:h, 200)

julia> @layout aside(:v, 300)

julia> @layout MyLayout

julia> uselayout() do canvas, widgets
            Box(:v) do
                canvas,
                widgets...
            end # Box
        end

```
"""
function init(title::String, width::Int = 400, height::Int = 400)
    @debug "Application initialized."
    app = AlexyaApp()
    app.title = title
    app.width = width
    app.height = height
    app.framerate = 65.0
    app.layout = canvasonly
    app.widgets = GtkWidget[]
    app.events = Dict{Symbol, Function}(:appclosed => appclosed)
    app.loop = true

    app.canvas = GtkCanvas()
    app.window = Window((), title, width, height)
    app.drawing = Drawing(width, height, :image)

    empty!(CURRENT_ALEXYA_APP)
    push!(CURRENT_ALEXYA_APP, app)

    return app
end

"""
        start(; [ async = false])

Starts the current `AlexyaApp`.

This function will setup all the events, apply the current `layout` and start the canvas drawing loop.
The `async` parameter determines if the loop should be asynchronous.
"""
function start(; async::Bool = false)
    app = current_app()
    layout = app.layout
    widgets = app.widgets
    window = app.window
    canvas = app.canvas

    add!(window, layout(canvas, widgets...))
    showall(window)
    setevent.([Event(T, callback) for (T, callback) in app.events])

    if async
        @async while app.loop == true
            Gtk.draw(canvas)
            sleep(1/app.framerate)
        end
    else
        while app.loop == true
            Gtk.draw(canvas)
            sleep(1/app.framerate)
        end
    end
end
"""
        framerate(fps::Real)

Sets the current `framerate`.

The default `framerate` is approximaly `60` frames/second.

!!! note

    The framerate is dynamic, wich means that you can
    set the framerate inside the `update` function or anywhere else.

"""
function framerate(fps::Real)
    app = current_app()
    app.framerate = fps
end

"""
        noLoop()

This function will prevent the drawing loop from start, wich will call the `update`
function just one time (or zero if you call on the `setup` function).
"""
function noLoop()
    app = current_app()
    app.loop = false
end

"""
        use(type::Symbol, func::Function)

`Alexya` works with events, and event listeners can be attached to the current
`application` with the `use` function or the `@use` macro.

Each event has a unique type (a `symbol`) that determines how it is fired.

!!! note

    The `@use` macro automatically determines the type of event
    by the name of the function that is passed as argument.

| Event        | Is fired when        | Usage                                                   |
|:-------------|:---------------------|:--------------------------------------------------------|
|`setup`       | The `app` starts     | Do something when the app is starting                   |
|`update`      | Every frame          | Draw every frame on the canvas using Luxor              |
|`appclosed`   | Window closed        | Do something when the window is closed                  |
|`keypress`    | Keyboard key pressed | Do something when a key is pressed                      |
|`mousepress`  | Canvas is clicked    | Do something when the canvas is clicked                 |
|`mouserelease`| Canvas is clicked    | Do something after the canvas is clicked                |
|`mousemotion` | Mouse is moving      | Do something when the mouse is moving inside the canvas |

# Examples

```julia
@use function setup()
    println("My application has started!")
end

@use function update()
    background("black")
    # Luxor drawing here
end

@use function keypress(event)
    @info "Key press event" key=event.keyval
end

@use function mousepress(event)
    @info "Mouse clicked" pos=event.pos
end

@use function mouserelease(event)
    @info "Mouse button release" pos=event.pos
end

@use function mousemotion(event)
    @info "Cursor moving" pos=event.pos
end

@use function appclosed()
    @info "Window closed"
end
```
"""
function use(type::Symbol, func::Function)
    app = current_app()
    app.events[type] = func
end

"""
        @use event

Adds a `event` listener on the current `application`.

The `@use` macro should be use in a function declaration expression.

The event `type` is automatically determined by the name of the function.
"""
macro use(ex)
    if ex isa Symbol
        name = QuoteNode(ex)
        :( use($name, $(esc(ex))) )
    elseif Meta.isexpr(ex, :function)
        name = ex.args[1].args[1]
        quote
            $(esc(ex))
            use($(QuoteNode(name)), $(esc(name)))
        end
    end
end

"""
        @window

Returns the current `window`.
"""
macro window()
    :( current_app(:window) )
end

"""
        @canvas

Returns the current canvas.
"""
macro canvas()
    :( current_app(:canvas) )
end

"""
        @width

Returns the current canvas `width`.
"""
macro width()
    :( current_app(:width) )
end

"""
        @height

Returns the current canvas `height`.
"""
macro height()
    :( current_app(:height) )
end

"""
        @framecount

Returns the current `frame count`.
"""
macro framecount()
    :( current_app(:framecount) )
end

"""
        @framerate

Returns the current `framerate`.
"""
macro framerate()
    :( current_app(:fps) )
end

"""
        @create widget

Adds a `widget` in the current `application`.
"""
macro create(ex)
    :( addwidget($(esc(ex))) )
end

"""
        @init title [width = 400, height = 400]

Shorthand macro for the `init` function.
"""
macro init(title)
    :( init($(esc(title)), 400, 400) )
end
"""
        @init title width height

Shorthand macro for the `init` function.
"""
macro init(title, width, height)
    :( init($(esc(title)), $(esc(width)), $(esc(height))) )
end


text:
export @interface

macro interface(ex)
    return interfacem(__module__, __source__, ex)
end

function interfacem(__module__::Module, __source__::LineNumberNode, ex::Symbol)
    esc(quote
        export $ex
        Core.@__doc__ $(Expr(:function, ex))
    end)
end


text:
using CSV
include("oeid.jl");
include("idanalysis.jl");

prefix = "stergiopulos1999_data/"
AoP = CSV.read(prefix*"apressure.csv", DataFrame, header=false, types=[Float64, Float64])
AoQ = CSV.read(prefix*"aflow.csv", DataFrame, header=false, types=[Float64, Float64])
rename!(AoP, ["timestamp", "vb_va"])
rename!(AoQ, ["timestamp", "vb_va"])
size(AoP) == size(AoQ) || @warn "Pressure and flow waveforms not equal in size."
AoQ[!,"vb_va"] = AoQ[!,"vb_va"].*(60/1000) # convert from mL/s to L/min
pltPQ = plotPQ([AoP], [AoQ], "stergiopulos1999", size=(800,400))

h = 0.005 # sampling rate

# input beats
# FOH fill to match initial and final value of beat sequence
smooth!(x,n=5) = append!(x,(collect(1:n-1).*(x[1]-x[end])/n).+x[end])
y = repeat(smooth!(AoP[:,:vb_va]),20)
u = repeat(smooth!(AoQ[:,:vb_va]),20)

# with parameters from paper
# θ = Rp,C,Rc,L
θ1999mLs, _ = params_stergiopulos1999total()
mLs2Lmin(θ) = Diagonal([1000/60, 60/1000, 1000/60, 1000/60])*θ
θ1999 = mLs2Lmin(θ1999mLs) # convert impedance units to match mmHg pressures and L/min flows
x01999 = solvex0(G(θ1999,h),u)
y1999, x1999, t, plt = plot_result(θ1999,x01999,y,u,h)
plt

# first ID attempt
Rp = mean(y)/mean(u);
θ0 = [Rp, 1, 1, 1]
x00 = solvex0(G(θ0,h), u)
θid = [1,1,1,1]
x0id = false
x0solve = true
id1999 = IdConditions(θ0,x00,θid,x0id,x0solve)
@time phat, res = identify(θ0,x00,u,y,h,θid=θid,x0id=x0id,x0solve=x0solve);
θhat, x0hat = deal(phat, id1999)
θ0, x00
x0hatsolve = solvex0(G(θhat,h),u)
hess = Optim.trace(res)[end].metadata["h(x)"]
svd(hess).U
_, _, _, plt = plot_result(θhat,x0hatsolve,y,u,h)
plot(plt, xlims=[0,1])

# check hessian and MSE with original parameters
@time phat1999,res1999 = identify(θ1999,solvex0(G(θ1999,h),u),u,y,h,θid=θid,x0id=x0id,x0solve=x0solve);
hess1999 = res1999.trace[1].metadata["h(x)"]
cond(hess1999)
svd(hess1999)

# multiple initialization, identify 4 element windkessel model
θ0s, θhats, x0hats, res = id_multiple_init(u, y, x00, θ0, 100, θid=θid,x0id=x0id,x0solve=x0solve, paramdist=[false,true,true,true]);

results4, nanratio = multiinitresults(θhats, x0hats, res);
results4[1,:hess]
results4[1,:hessSVD].U
best1999_4 = results4[1,[1:7...,11,12]]

θbest4 = convert(Array, results4[1, 1:length(θ0)]);
[θ0, round.(θbest4, sigdigits=3)]
x0best4 = Array(results4[1,["x01", "x02"]]);
[x00, round.(x0best4, sigdigits=3)]
_, _, t, plt_best = plot_result(θbest4, x0best4, y, u, h)
plot(plt_best, size=(1200,400))
# savefig("results/best_fit_1999_4_$(Dates.format(now(), "yyyymmdd")).png")
# CSV.write("results/results_1999_4_$(Dates.format(now(), "yyyymmdd")).csv",results4)

# plot parameter histograms before and after fit
θlabels = ["Rp" "C" "Rc" "L"]
fithistograms(θ0s,θhats,θlabels)
# savefig("./results/param_conv_" * Dates.format(now(), "yyyymmdd-HHMMSS") * ".png")

yhat4, = lsim(G(θbest4,h),u,h.*(1:length(u)),x0best4)
yhat4 = yhat4[:]
w = exp10.(range(-3, stop=3, length=500))
mag4, phase4, w = bode(G(θbest4,h,discretize=false),w)
mag4plt = plot(w,mag4[:],xaxis=:log, yaxis=:log, label="WK4",title="1999");
phase4plt = plot(w,phase4[:], xaxis=:log, label="WK4");
plot(mag4plt, phase4plt, layout=(2,1))

# multiple initialization, identify 3 element windkessel model
θ0 = copy(θbest4[1:3])
x00 = solvex0(G3(θ0,h), u)
θid = [1,1,1]
x0id = false
x0solve = true
θ0s, θhats, x0hats, res = id_multiple_init(u, y, x00, θ0, 100, θid=θid,x0id=x0id,x0solve=x0solve, G_=G3, paramdist=[false,true,true]);

results3, nanratio = multiinitresults(θhats, x0hats, res);
θbest3 = Array(results3[1,1:3])
x0best3 = [results3[1,"x01"]]
results3[1,:hess]
results3[1,:hessSVD].U
best1999_3 = results3[1,["Rp", "C", "Rc", "x01", "MSE","hesscond"]]
# CSV.write("results/results_1999_3_$(Dates.format(now(), "yyyymmdd")).csv",results3)
# results3 = CSV.read("results/results_1999_3_20210309.csv", DataFrame)

yhat3, = lsim(G3(θbest3,h),u,h.*(1:length(u)),x0best3)
yhat3 = yhat3[:]
mag3, phase3, w = bode(G3(results3[1,1:3],h,discretize=false),w)
gr()
mag3plt = plot(w,mag3[:],xaxis=:log, yaxis=:log, label="WK3",title="1999");
phase3plt = plot(w,phase3[:], xaxis=:log, label="WK3");
plot(mag3plt, phase3plt, layout=(2,1))

# multiple initialization, identify 2 element windkessel model
θ0 = copy(θbest4[1:2])
x00 = solvex0(G2(θ0,h), u)
θid = [1,1]
θ0s, θhats, x0hats, res = id_multiple_init(u, y, x00, θ0, 100, θid=θid,x0id=x0id,x0solve=x0solve, G_=G2, paramdist=[false,true]);

results2, nanratio = multiinitresults(θhats, x0hats, res);
θbest2 = Array(results2[1,1:2])
x0best2 = [results2[1,"x01"]]
results2[1,:hess]
results2[1,:hessSVD].U
best1999_2 = results2[1,["Rp", "C", "x01", "MSE","hesscond"]]
# CSV.write("results/results_1999_2_$(Dates.format(now(), "yyyymmdd")).csv",results2)
# results2 = CSV.read("results/results_1999_2_20210309", DataFrame)

yhat2, = lsim(G2(θbest2,h),u,h.*(1:length(u)),x0best2)
yhat2 = yhat2[:]
mag2, phase2, w = bode(G2(θbest2,h,discretize=false),w)
gr()
mag2plt = plot(w,mag2[:],xaxis=:log, yaxis=:log, label="WK2",title="1999");
phase2plt = plot(w,phase2[:], xaxis=:log, label="WK2");
plot(mag2plt, phase2plt, layout=(2,1))

y = AoP.vb_va
u = AoQ.vb_va
v, = pvloop(y,u,h)
y1999 = y1999[1:length(y),1]
# load results
# results4 = CSV.read("results/results_1999_4_20210309.csv", DataFrame)
# svd1999 = eval(Meta.parse(results4[1,:hessSVD]))
# x0best4 = solvex0(G(θbest4,h),u)
# θbest4 = Array(results4[1,1:4])
# yhat4 = simulate(θbest4,x0best4,u,h)[:,1]
svd1999 = results4[1,:hessSVD]
plot([y,y1999,yhat4],xlims=(0,300))
plot([v,v,v],[y,y1999,yhat4])

# export bode data for pgfplots
comment= *(
    "# fit to data: published data stergiopulos1999total doi = {10.1152/ajpheart.1999.276.1.H81}, Figure 4A type A beat.\n",
    "# mag4,phase4: 4-element parallel Windkessel model\n",
    "# θhat = [Rc,L,C,Rp]\n",
    "# θhat = $(θbest4)\n",
    "# mag3,phase3: 3-element parallel Windkessel model\n",
    "# θhat = [Rc,C,Rp]\n",
    "# θhat = $(θbest3)\n",
    "# mag2,phase2: 2-element parallel Windkessel model\n",
    "# θhat = [C,Rp]\n",
    "# θhat = $(θbest2)\n")
filename = "bodesterg1999.csv"
open(filename, "w") do f
    write(f, comment)
end
CSV.write(filename, DataFrame(w = w,
mag4 = mag4[:], phase4 = phase4[:],
mag3 = mag3[:], phase3 = phase3[:],
mag2 = mag2[:], phase2 = phase2[:]),header=true,append=true)

# export data for pgfplots
beatstart = 1
beatend = length(AoP.vb_va)
comment = *("# data: from published data stergiopulos1999total doi = {10.1152/ajpheart.1999.276.1.H81}, Figure 4A type A beat.\n",
"# t: time in seconds\n",
"# q: published measured Aortic flow in L/min, parsed using https://automeris.io/WebPlotDigitizer\n",
"# p: published measured Aortic pressure in mmHg, parsed using https://automeris.io/WebPlotDigitizer\n",
"# phat1999: modelled Aortic pressure in mmHg, using parameters published in paper θhat1999\n",
"# phat4: modelled Aortic pressure in mmHg\n",
"# v: measured volume through the aorta in ml\n",
"# epsilon4: error signal of the model p-phat4\n",
"# 4-element parallel Windkessel model\n",
"# θhat = [Rp,C,Rc,L]\n",
"# θhat1999 = $(θ1999); solved for x0hat1 = $(x01999)\n",
"# θhat = $(θbest4); x0hat2 = $(x0best4)\n",
"# as well as the corresonding phat and epsilon values for\n",
"# 3-element parallel Windkessel model\n",
"# θhat = [Rc,C,Rp]\n",
"# θhat = $(θbest3); x0hat = $(x0best3)\n",
"# 2-element parallel Windkessel model\n",
"# θhat = [C,Rp]\n",
"# θhat = $(θbest2); x0hat = $(x0best2)\n")
filename = "pqvsterg1999.csv"
open(filename, "w") do f
    write(f, comment)
end
CSV.write(filename, DataFrame(t = round.(collect(0:beatend-beatstart)*h, digits=5), q = u[beatstart:beatend],
p = y[beatstart:beatend],
phat1999 = y1999[beatstart:beatend],
phat4 = yhat4[beatstart:beatend],
phat3 = yhat3[beatstart:beatend],
phat2 = yhat2[beatstart:beatend],
v = v[beatstart:beatend].-v[beatstart],
epsilon1999 = y[beatstart:beatend].-y1999[beatstart:beatend],
epsilon4 = y[beatstart:beatend].-yhat4[beatstart:beatend],
epsilon3 = y[beatstart:beatend].-yhat3[beatstart:beatend],
epsilon2 = y[beatstart:beatend].-yhat2[beatstart:beatend]),
header=true,append=true)

max_stars_repo_name: hpigot/windkessel-id

max_issues_repo_name: hpigot/windkessel-id

max_forks_repo_name: hpigot/windkessel-id


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

mutable struct IoK8sApiSettingsV1alpha1PodPresetSpec <: SwaggerModel
    env::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1EnvVar} } # spec name: env
    envFrom::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1EnvFromSource} } # spec name: envFrom
    selector::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1LabelSelector } # spec name: selector
    volumeMounts::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1VolumeMount} } # spec name: volumeMounts
    volumes::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1Volume} } # spec name: volumes

    function IoK8sApiSettingsV1alpha1PodPresetSpec(;env=nothing, envFrom=nothing, selector=nothing, volumeMounts=nothing, volumes=nothing)
        o = new()
        validate_property(IoK8sApiSettingsV1alpha1PodPresetSpec, Symbol("env"), env)
        setfield!(o, Symbol("env"), env)
        validate_property(IoK8sApiSettingsV1alpha1PodPresetSpec, Symbol("envFrom"), envFrom)
        setfield!(o, Symbol("envFrom"), envFrom)
        validate_property(IoK8sApiSettingsV1alpha1PodPresetSpec, Symbol("selector"), selector)
        setfield!(o, Symbol("selector"), selector)
        validate_property(IoK8sApiSettingsV1alpha1PodPresetSpec, Symbol("volumeMounts"), volumeMounts)
        setfield!(o, Symbol("volumeMounts"), volumeMounts)
        validate_property(IoK8sApiSettingsV1alpha1PodPresetSpec, Symbol("volumes"), volumes)
        setfield!(o, Symbol("volumes"), volumes)
        o
    end
end # type IoK8sApiSettingsV1alpha1PodPresetSpec

const _property_map_IoK8sApiSettingsV1alpha1PodPresetSpec = Dict{Symbol,Symbol}(Symbol("env")=>Symbol("env"), Symbol("envFrom")=>Symbol("envFrom"), Symbol("selector")=>Symbol("selector"), Symbol("volumeMounts")=>Symbol("volumeMounts"), Symbol("volumes")=>Symbol("volumes"))
const _property_types_IoK8sApiSettingsV1alpha1PodPresetSpec = Dict{Symbol,String}(Symbol("env")=>"Vector{IoK8sApiCoreV1EnvVar}", Symbol("envFrom")=>"Vector{IoK8sApiCoreV1EnvFromSource}", Symbol("selector")=>"IoK8sApimachineryPkgApisMetaV1LabelSelector", Symbol("volumeMounts")=>"Vector{IoK8sApiCoreV1VolumeMount}", Symbol("volumes")=>"Vector{IoK8sApiCoreV1Volume}")
Base.propertynames(::Type{ IoK8sApiSettingsV1alpha1PodPresetSpec }) = collect(keys(_property_map_IoK8sApiSettingsV1alpha1PodPresetSpec))
Swagger.property_type(::Type{ IoK8sApiSettingsV1alpha1PodPresetSpec }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiSettingsV1alpha1PodPresetSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiSettingsV1alpha1PodPresetSpec }, property_name::Symbol) =  _property_map_IoK8sApiSettingsV1alpha1PodPresetSpec[property_name]

function check_required(o::IoK8sApiSettingsV1alpha1PodPresetSpec)
    true
end

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


text:
struct JspritJavaTypes
    JCoord::DataType
    JLocation::DataType
    JLocationBuilder::DataType
    JTW::DataType
    JAbstractJob::DataType
    JJob::DataType
    JPickup::DataType
    JPickupBuilder::DataType
    JDelivery::DataType
    JDeliveryBuilder::DataType
    JShipment::DataType
    JShipmentBuilder::DataType
    JService::DataType
    JServiceBuilder::DataType
    JVTypeI::DataType
    JVType::DataType
    JVTypeBuilder::DataType
    JAbstractVehicle::DataType
    JVehicleI::DataType
    JVehicle::DataType
    JVehicleBuilder::DataType
    JCostsI::DataType
    JCosts::DataType
    JCostsBuilder::DataType
    FleetSize::DataType
    JVrp::DataType
    JVrpBuilder::DataType
    JAlgorithm::DataType
    JSP::DataType
    JSPBuilder::DataType
    JSolution::DataType
    JSolutions::DataType
    JVehicleRoute::DataType
    JTourActivity::DataType
    JJobActivity::DataType
    JStartActivity::DataType
    JEndActivity::DataType
    JCollection::DataType
    JList::DataType
    JMap::DataType
    function JspritJavaTypes()
        jcoord = JavaCall.@jimport com.graphhopper.jsprit.core.util.Coordinate
        jlocation = JavaCall.@jimport com.graphhopper.jsprit.core.problem.Location
        jlocationbuilder = JavaCall.@jimport com.graphhopper.jsprit.core.problem.Location$Builder
        jtw = JavaCall.@jimport com.graphhopper.jsprit.core.problem.solution.route.activity.TimeWindow
        jabsjob = JavaCall.@jimport com.graphhopper.jsprit.core.problem.AbstractJob
        jjob = JavaCall.@jimport com.graphhopper.jsprit.core.problem.job.Job
        jpickup = JavaCall.@jimport com.graphhopper.jsprit.core.problem.job.Pickup
        jpickupbuilder = JavaCall.@jimport com.graphhopper.jsprit.core.problem.job.Pickup$Builder
        jdelivery = JavaCall.@jimport com.graphhopper.jsprit.core.problem.job.Delivery
        jdeliverybuilder = JavaCall.@jimport com.graphhopper.jsprit.core.problem.job.Delivery$Builder
        jshipment = JavaCall.@jimport com.graphhopper.jsprit.core.problem.job.Shipment
        jshipmentbuilder = JavaCall.@jimport com.graphhopper.jsprit.core.problem.job.Shipment$Builder
        jservice = JavaCall.@jimport com.graphhopper.jsprit.core.problem.job.Service
        jservicebuilder = JavaCall.@jimport com.graphhopper.jsprit.core.problem.job.Service$Builder
        jvtypei = JavaCall.@jimport com.graphhopper.jsprit.core.problem.vehicle.VehicleType
        jvtypeimpl = JavaCall.@jimport com.graphhopper.jsprit.core.problem.vehicle.VehicleTypeImpl
        jvtypebuilder = JavaCall.@jimport com.graphhopper.jsprit.core.problem.vehicle.VehicleTypeImpl$Builder
        jabsvehicle = JavaCall.@jimport com.graphhopper.jsprit.core.problem.AbstractVehicle
        jvehiclei = JavaCall.@jimport com.graphhopper.jsprit.core.problem.vehicle.Vehicle
        jvehicleimpl = JavaCall.@jimport com.graphhopper.jsprit.core.problem.vehicle.VehicleImpl
        jvehiclebuilder = JavaCall.@jimport com.graphhopper.jsprit.core.problem.vehicle.VehicleImpl$Builder
        jcostsi = JavaCall.@jimport com.graphhopper.jsprit.core.problem.cost.VehicleRoutingTransportCosts
        jcosts = JavaCall.@jimport com.graphhopper.jsprit.core.util.TimeDependentVehicleRoutingTransportCosts
        jcostsbuilder = JavaCall.@jimport com.graphhopper.jsprit.core.util.TimeDependentVehicleRoutingTransportCosts$Builder
        fleetsize = JavaCall.@jimport com.graphhopper.jsprit.core.problem.VehicleRoutingProblem$FleetSize
        jvrp = JavaCall.@jimport com.graphhopper.jsprit.core.problem.VehicleRoutingProblem
        jvrpbuilder = JavaCall.@jimport com.graphhopper.jsprit.core.problem.VehicleRoutingProblem$Builder
        jalgorithm = JavaCall.@jimport com.graphhopper.jsprit.core.algorithm.VehicleRoutingAlgorithm
        jsp = JavaCall.@jimport com.graphhopper.jsprit.core.algorithm.box.Jsprit
        jspbuilder = JavaCall.@jimport com.graphhopper.jsprit.core.algorithm.box.Jsprit$Builder
        jsolution = JavaCall.@jimport com.graphhopper.jsprit.core.problem.solution.VehicleRoutingProblemSolution
        jsolutions = JavaCall.@jimport com.graphhopper.jsprit.core.util.Solutions
        jvehicleroute = JavaCall.@jimport com.graphhopper.jsprit.core.problem.solution.route.VehicleRoute
        jtouractivity = JavaCall.@jimport com.graphhopper.jsprit.core.problem.solution.route.activity.TourActivity
        jjobactivity = JavaCall.@jimport com.graphhopper.jsprit.core.problem.solution.route.activity.TourActivity$JobActivity
        jstartactivity = JavaCall.@jimport com.graphhopper.jsprit.core.problem.solution.route.activity.Start
        jendactivity = JavaCall.@jimport com.graphhopper.jsprit.core.problem.solution.route.activity.End
        jcollection = JavaCall.@jimport java.util.Collection
        jlist = JavaCall.@jimport java.util.List
        jmap = JavaCall.@jimport java.util.Map
        new(
            jcoord,
            jlocation,
            jlocationbuilder,
            jtw,
            jabsjob,
            jjob,
            jpickup,
            jpickupbuilder,
            jdelivery,
            jdeliverybuilder,
            jshipment,
            jshipmentbuilder,
            jservice,
            jservicebuilder,
            jvtypei,
            jvtypeimpl,
            jvtypebuilder,
            jabsvehicle,
            jvehiclei,
            jvehicleimpl,
            jvehiclebuilder,
            jcostsi,
            jcosts,
            jcostsbuilder,
            fleetsize,
            jvrp,
            jvrpbuilder,
            jalgorithm,
            jsp,
            jspbuilder,
            jsolution,
            jsolutions,
            jvehicleroute,
            jtouractivity,
            jjobactivity,
            jstartactivity,
            jendactivity,
            jcollection,
            jlist,
            jmap
        )
    end
end


text:
function residuals_SS_bounded(x; par::Pars_v3)
    Y_guess = x[1]
    w_guess = x[2]
    w_P_guess = x[3]
    w_G_guess = x[4]
    #p_D_guess = x[5]^2

    @set! par.Y = Y_guess
    @set! par.w = w_guess
    @set! par.w_P = w_P_guess
    @set! par.w_G = w_G_guess
    @set! par.p_D = find_p_D(par)

    sol_BS = solve_firm_problem(par; type = :BS)
    sol_LS = solve_firm_problem(par; type = :LS)
    sol_LU = solve_firm_problem(par; type = :LU)
    sol_BU = solve_firm_problem(par; type = :BU)

    l_P_BS, l_G_BS, l_BS, D_E_BS, D_I_BS, D_BS = get_solution(sol_BS; par, type = :BS)
    l_P_LS, l_G_LS, l_LS, D_E_LS, D_I_LS, D_LS = get_solution(sol_LS; par, type = :LS)
    l_P_LU, l_G_LU, l_LU, D_E_LU, D_I_LU, D_LU = get_solution(sol_LU; par, type = :LU)
    l_P_BU, l_G_BU, l_BU, D_E_BU, D_I_BU, D_BU = get_solution(sol_BU; par, type = :BU)

    #D_BS = bundle(D_I_BS, D_E_BS; par)
    #D_LS = bundle(D_I_LS, D_E_LS; par)
    #D_LU = bundle(D_I_LU, D_E_LU; par)
    #D_BU = bundle(D_I_BU, D_E_BU; par)

    Y_BS = firm_output(l_P_BS, D_BS, l_BS; par, type = :BS)
    Y_LS = firm_output(l_P_LS, D_LS, l_LS; par, type = :LS)
    Y_LU = firm_output(l_P_LU, D_LU, l_LU; par, type = :LU)
    Y_BU = firm_output(l_P_BU, D_BU, l_BU; par, type = :BU)

    Y = aggregate_output(Y_BS, Y_LS, Y_LU, Y_BU; par)

    res = zeros(4)

    res[1] = (Y - Y_guess)
    res[2] = labor_market_clearing(l_BS, l_LS, l_LU, l_BU; par)
    res[3] = labor_market_clearing(l_P_BS, l_P_LS, l_P_LU, l_P_BU; par)
    res[4] = labor_market_clearing(l_G_BS, l_G_LS, l_G_LU, l_G_BU; par)

    # l_G = [l_G_BS, l_G_LS, l_G_LU, l_G_BU]
    # D_E = [D_E_BS, D_E_LS, D_E_LU, D_E_BU]
    # D_I = [D_I_BS, D_I_LS, D_I_LU, D_I_BU]

    # res[5] = data_market_clearing(l_G, D_E, D_I; par)
    @show diff = sum(res .^ 2)
    return sum(res .^ 2)
end

# function residuals_SS_single(x; par::Pars_v3)
#     Y_guess = x[1]^2
#     w_guess = x[2]^2
#     w_P_guess = x[3]^2
#     w_G_guess = x[4]^2
#     p_D_guess = x[5]^2

#     @set! par.Y = Y_guess
#     @set! par.w = w_guess
#     @set! par.w_P = w_P_guess
#     @set! par.w_G = w_G_guess
#     @set! par.p_D = p_D_guess

#     sol_BS = solve_firm_problem(par; type = :BS)
#     sol_LS = solve_firm_problem(par; type = :LS)
#     sol_LU = solve_firm_problem(par; type = :LU)
#     sol_BU = solve_firm_problem(par; type = :BU)

#     l_P_BS, l_G_BS, l_BS, D_E_BS, D_I_BS, D_BS = get_solution(sol_BS; par, type =:BS)
#     l_P_LS, l_G_LS, l_LS, D_E_LS, D_I_LS, D_LS = get_solution(sol_LS; par, type =:LS)
#     l_P_LU, l_G_LU, l_LU, D_E_LU, D_I_LU, D_LU = get_solution(sol_LU; par, type =:LU)
#     l_P_BU, l_G_BU, l_BU, D_E_BU, D_I_BU, D_BU = get_solution(sol_BU; par, type =:BU)

#     #D_BS = bundle(D_I_BS, D_E_BS; par)
#     #D_LS = bundle(D_I_LS, D_E_LS; par)
#     #D_LU = bundle(D_I_LU, D_E_LU; par)
#     #D_BU = bundle(D_I_BU, D_E_BU; par)

#     Y_BS = firm_output(l_P_BS, D_BS, l_BS; par, type = :BS)
#     Y_LS = firm_output(l_P_LS, D_LS, l_LS; par, type = :LS)
#     Y_LU = firm_output(l_P_LU, D_LU, l_LU; par, type = :LU)
#     Y_BU = firm_output(l_P_BU, D_BU, l_BU; par, type = :BU)

#     Y = aggregate_output(Y_BS, Y_LS, Y_LU, Y_BU; par)

#     res = zeros(5)

#     res[1] = (Y - Y_guess)
#     res[2] = labor_market_clearing(l_BS, l_LS, l_LU, l_BU; par)
#     res[3] = labor_market_clearing(l_P_BS, l_P_LS, l_P_LU, l_P_BU; par)
#     res[4] = labor_market_clearing(l_G_BS, l_G_LS, l_G_LU, l_G_BU; par)

#     l_G = [l_G_BS, l_G_LS, l_G_LU, l_G_BU]
#     D_E = [D_E_BS, D_E_LS, D_E_LU, D_E_BU]
#     D_I = [D_I_BS, D_I_LS, D_I_LU, D_I_BU]

#     res[5] = data_market_clearing(l_G, D_E, D_I; par)

#     return res
# end

"""
    Solving the model with a general solution algorithm does not work performantly.
"""
function find_steady_state(par::Pars_v3)
    f(x) = residuals_SS(x; par)
    initial_x = [par.Y, par.w, par.w_P, par.w_G]
    optimize(f, initial_x)
end

function constraint_data_selling(D_S)
    if D_S > 0
        return 0
    elseif D_S <= 0
        100000 * D_S^2
    else
        error("I should never be here")
    end
end

function wage_int_goods2(K; par)
    K_BS = K[1]
    K_LS = K[2]
    K_LU = K[3]
    K_BU = K[4]

    w_BS = weight_of_type(:BS; par)
    w_LS = weight_of_type(:LS; par)
    w_LU = weight_of_type(:LU; par)
    w_BU = weight_of_type(:BU; par)

    return par.Y^par.α_Y * par.α_L_hat * (w_BS * K_BS + w_LS * K_LS + w_LU * K_LU + w_BU * K_BU)
end

"""
    Wage for data generation
"""
function wage_gen_data2(Π, D; par)
    Π_BS = Π[1] # firm revenue
    Π_LS = Π[2] # firm revenue
    Π_LU = Π[3] # firm revenue
    Π_BU = Π[4] # firm revenue

    D_BS = D[1] # data bundle
    D_LS = D[2] # data bundle
    D_LU = D[3] # data bundle
    D_BU = D[4] # data bundle

    w_BS = weight_of_type(:BS; par) # weights
    w_LS = weight_of_type(:LS; par) # weights
    w_LU = weight_of_type(:LU; par) # weights
    w_BU = weight_of_type(:BU; par) # weights

    A_G_BS = A_G_of_type(:BS; par) # prod of data generation
    A_G_LS = A_G_of_type(:LS; par) # prod of data generation
    A_G_LU = A_G_of_type(:LU; par) # prod of data generation
    A_G_BU = A_G_of_type(:BU; par) # prod of data generation

    c_BS = Π_BS / D_BS * A_G_BS # Put some things together
    c_LS = Π_LS / D_LS * A_G_LS # Put some things together
    c_LU = Π_LU / D_LU * A_G_LU # Put some things together
    c_BU = Π_BU / D_BU * A_G_BU # Put some things together

    return par.α_K_hat * ((1 - par.γ_S) * w_BS * c_BS + (1 - par.γ_S) * w_LS * c_LS + (1 - par.γ_U) * w_LU * c_LU + (1 - par.γ_U) * w_BU * c_BU)
end

"""
    Return a vector containing all types
"""
function weights_of_types(par::Pars_v3)
    w_BS = weight_of_type(:BS; par)
    w_LS = weight_of_type(:LS; par)
    w_LU = weight_of_type(:LU; par)
    w_BU = weight_of_type(:BU; par)

    return [w_BS, w_LS, w_LU, w_BU]
end


text:
"""
`AllenCahnProblem1D` - initialize the problem data structure
"""
function AllenCahnProblem1D(a, b, ϵ, nₓ, nₜ, t_max, left_bc::TLBC, right_bc::TRBC, u₀::TU) where {TLBC <: NeumannBC, TRBC <: NeumannBC, TU <: Function}
    @assert a < b "Cannot construct: Invalid interval!"
    @assert ϵ > 0 "Cannot construct: ϵ must be positive!"
    @assert (nₓ > 0 && nₜ > 0) "Cannot construct: nₓ and nₜ must both be positive!"

    Δx = (b - a)/(nₓ + 1)
    x = [a + i*Δx for i ∈ 0:nₓ + 1]
    Δt = t_max/nₜ
    f(x) = x - x^3
    u = u₀.(x)
    return AllenCahnProblem1D(a, b, ϵ, nₓ, Δx, x, nₜ, t_max, Δt, left_bc, right_bc, spzeros(nₓ + 2, nₓ + 2), zeros(nₓ + 2), f, u)
end

"""
`AllenCahnProblem1D` - initialize the problem data structure
"""
function AllenCahnProblem1D(a, b, ϵ, nₓ, nₜ, t_max, left_bc::TLBC, right_bc::TRBC, u₀::TU) where {TLBC <: PeriodicBC, TRBC <: PeriodicBC, TU <: Function}
    @assert a < b "Cannot construct: Invalid interval!"
    @assert ϵ > 0 "Cannot construct: ϵ must be positive!"
    @assert (nₓ > 0 && nₜ > 0) "Cannot construct: nₓ and nₜ must both be positive!"

    Δx = (b - a)/(nₓ + 1)
    x = [a + i*Δx for i ∈ 0:nₓ]
    Δt = t_max/nₜ
    f(x) = x - x^3
    u = u₀.(x)
    return AllenCahnProblem1D(a, b, ϵ, nₓ, Δx, x, nₜ, t_max, Δt, left_bc, right_bc, spzeros(nₓ + 1, nₓ + 1), zeros(nₓ + 1), f, u)
end


text:
"""
    AbstractScenario

Abstract supertype for scenario objects.
"""
abstract type AbstractScenario end
"""
    Probability

A type-safe wrapper for `Float64` used to represent probability of a scenario occuring.
"""
mutable struct Probability
    π::Float64
end
"""
    probability(scenario::AbstractScenario)

Return the probability of `scenario` occuring.

Is always defined for scenarios created through @scenario. Other user defined scenario types must implement this method to generate a proper probability. The default behaviour is to assume that `scenario` has a `probability` field of type [`Probability`](@ref)

See also: [`Probability`](@ref)
"""
probability(scenario::AbstractScenario)::Float64 = scenario.probability.π
"""
    probability(scenarios::Vector{<:AbstractScenario})

Return the probability of that any scenario in the collection `scenarios` occurs.
"""
probability(scenarios::Vector{<:AbstractScenario}) = sum(probability.(scenarios))
"""
    set_probability!(scenario::AbstractScenario, probability::AbstractFloat)

Set the probability of `scenario` occuring.

Is always defined for scenarios created through @scenario. Other user defined scenario types must implement this method.
"""
function set_probability!(scenario::AbstractScenario, π::AbstractFloat)
    scenario.probability.π = π
end
function Base.zero(::Type{S}) where S <: AbstractScenario
    error("zero not implemented for scenario type: ", S)
end
function Base.show(io::IO, scenario::S) where S <: AbstractScenario
    print(io, "$(S.name.name) with probability $(probability(scenario))")
    scenariotext(io, scenario)
    return io
end
"""
    scenariotext(io::IO, scenario::AbstractScenario)

Custom printout called when printing `scenario`.
"""
function scenariotext(io::IO, scenario::AbstractScenario)
    return io
end
function Base.getindex(ξ::AbstractScenario, field::Symbol)
    return getfield(ξ, field)
end
"""
    ExpectedScenario{S <: AbstractScenario}

Wrapper type around an `AbstractScenario`. Should for convenience be used as the result of a call to `expected`.

See also [`expected`](@ref)
"""
struct ExpectedScenario{S <: AbstractScenario} <: AbstractScenario
    scenario::S

    (::Type{ExpectedScenario})(scenario::AbstractScenario) = new{typeof(scenario)}(scenario)
end
function Base.show(io::IO, scenario::ExpectedScenario{S}) where S <: AbstractScenario
    print(io, "Expected scenario of type $(S.name.name)")
    scenariotext(io, scenario.scenario)
    return io
end
"""
    expected(scenarios::Vector{<:AbstractScenario})

Return the expected scenario out of the collection `scenarios` in an [`ExpectedScenario`](@ref) wrapper.

This is defined through classical expectation: sum([probability(s)*s for s in scenarios]), and is always defined for scenarios created through @scenario, if the requested fields support it.

Otherwise, user-defined scenario types must implement this method for full functionality.

See also [`ExpectedScenario`](@ref)
"""
function expected(::Vector{S}) where S <: AbstractScenario
    error("Expectation not implemented for scenario type: ", S)
end
"""
    Scenario

Conveniece type that adheres to the [`AbstractScenario`](@ref) abstraction. Useful when uncertain parameters are a finite set of scalar values.
"""
struct Scenario{T} <: AbstractScenario
    probability::Probability
    data::T

    function Scenario(data::T; probability::AbstractFloat = 1.0) where T
        return new{T}(Probability(probability), data)
    end

    function Scenario(; probability::AbstractFloat = 1.0, kw...)
        data = values(kw)
        NT = typeof(data)
        return new{NT}(Probability(probability), data)
    end
end
function Base.zero(::Type{Scenario{NT}}) where NT <: NamedTuple
    return Scenario(NamedTuple{Tuple(NT.names)}(zero.(NT.types)); probability = 1.0)
end
function scenariotext(io::IO, scenario::Scenario{NT}) where NT <: NamedTuple
    for (k,v) in pairs(scenario.data)
        print(io, "\n  $k: $v")
    end
    return io
end
function scenariotext(io::IO, scenario::Scenario)
    print(io, " and underlying data:\n\n")
    print(io, scenario.data)
    return io
end

Scenarios{S <: AbstractScenario} = Vector{S}
ScenarioTypes{N} = NTuple{N, Union{DataType, UnionAll}}

function expected(scenarios::Vector{<:Scenario{NT}}) where NT <: NamedTuple
    isempty(scenarios) && return StochasticPrograms.ExpectedScenario(zero(Scenario{NT}))
    expected = reduce(scenarios) do s₁, s₂
        keys(s₁.data) == keys(s₂.data) || error("Iconsistent scenarios. $(keys(s₁)) and $(keys(s₂)) do not match.")
        Scenario(NamedTuple{Tuple(keys(s₁.data))}([probability(s₁)*x + probability(s₂)*y for (x,y) in zip(values(s₁.data), values(s₂.data))]); probability = 1.0)
    end
    return StochasticPrograms.ExpectedScenario(expected)
end


text:
import Voronoi

using Random
include("Draw.jl")

function randf(start, finish, n)
   v = rand(n)
   return map(x -> start + x*(finish-start), v)
end

function demoFortune()
   WIDTH = 300.0
   HEIGHT = 100.0
   n = 2
   for i in 1:1000
      #println("Random seed: ", i)
      #Random.seed!(i)
      #points = convert(Array{Tuple{Real, Real}}, collect(zip(randf(1, WIDTH-1, n), randf(1, HEIGHT-1, n))))
      #points = convert(Array{Tuple{Real, Real}}, [(90.0, 90.0), (90.0, 50), (90, 30)])
#      points = convert(Array{Tuple{Real, Real}}, [(90.0, 90.0), (80.0, 50), (70.0, 55.0), (30.0, 85.0), (30, 45), (35.0, 44.0),
#                                                  (1, 90),
#                                                  (5, 90), (10, 90), (15, 90),
#                                                  (10, 55), (10,60), (10, 20), (80, 20), (30, 20),
#                                                  (10, 1), (20, 1), (30, 1),
#                                                 ])
#      points = convert(Array{Tuple{Real, Real}}, [
#                                                  (100*0.28, 100*0.8), (100*0.28, 100*0.6), (100*0.28, 100*0.4), (100*0.28, 100*0.2),
#                                                  (100*0.8,  100*0.8), (100*0.8,  100*0.6), (100*0.8,  100*0.4), (100*0.8,  100*0.2),
#                                                  (100*1.32, 100*0.8), (100*1.32, 100*0.6), (100*1.32, 100*0.4), (100*1.32, 100*0.2),
#                                                  (100*1.84, 100*0.8), (100*1.84, 100*0.6), (100*1.84, 100*0.4), (100*1.84, 100*0.2),
#                                                  (100*2.36, 100*0.8), (100*2.36, 100*0.6), (100*2.36, 100*0.4), (100*2.36, 100*0.2),
#                                                  (100*2.84, 100*0.8), (100*2.84, 100*0.6), (100*2.84, 100*0.4), (100*2.84, 100*0.2)
#                                                 ])
      points = convert(Array{Tuple{Real, Real}}, [ (100*0.2219056800797507, 100*0.7000658172436546), (100*0.49169960385155936, 100*0.6636482987025588), (100*0.7618642131758435, 100*0.7432736698473899), (100*0.8846773855975993, 100*0.6317667493965284) ])

      # test for repeated points
      repeat = false
      for i in 1:length(points)
         for j in 1:length(points)
            if j == i
               continue
            end

            if points[i][1] == points[j][1] && points[i][2] == points[j][2]
               repeat = true
            end
         end
      end
      if repeat
         println("REPEATED POINTS!")
         continue
      end

      V, T, Q = Voronoi.Fortune.init(points, WIDTH, HEIGHT)
      command = nothing
      ly = HEIGHT
      Draw.init(WIDTH, HEIGHT)

      while (event = Voronoi.EventQueue.pop(Q)) != nothing
         Voronoi.Fortune.handleEvent(V, T, Q, event) # multiple dispatch decides if it's a site event or circle event

         ly = Voronoi.EventQueue.coordinates(event)[2] # sweep line height
         Draw.fortuneIteration(V, T, Q, points, ly)
         if command != "a"
            println("Press Return for a step or enter \"a\" to animate until the end.")
            command = readline(stdin)
         else
            sleep(0.001)
         end
      end

      Voronoi.Intersect.intersect(V, Voronoi.Intersect.Rectangle(WIDTH, HEIGHT))
      Draw.init(WIDTH, HEIGHT)
      Draw.voronoiDiagram(V)

      println("Press Return to finish.")
      readline(stdin)
      n += 1
   end
end

demoFortune()

max_stars_repo_name: arthursantana/coveringWithCircles

max_issues_repo_name: arthursantana/coveringWithCircles

max_forks_repo_name: arthursantana/coveringWithCircles


text:
using Diderot: Arc, Node

@testset "struct equality" begin
    @test Arc(1, true, 1.0) == Arc(1, true, 1.0)
    @test Arc(1, true, 1.0) != Arc(0, true, 1.0)
    @test Arc(1, true, 1.0) != Arc(1, false, 1.0)
    @test Arc(1, true, 1.0) != Arc(1, true, 1.1)

    @test Arc(BitSet(1:3), 1, 0.0) == Arc(BitSet(1:3), 1, 0.0)

    N = Node{Int,Int,Int}
    @test N(1, Arc(1, 1, 1), true) == N(1, Arc(1, 1, 1), true)
    @test N(1, Arc(1, 1, 1), true) != N(1, nothing, true)
    @test N(1, Arc(1, 1, 1), true) != N(1, Arc(1, 1, 1), false)
    @test N(1, nothing, true) == N(1, nothing, true)
end