{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE StandaloneDeriving #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE DeriveTraversable #-}
{-# LANGUAGE DeriveGeneric #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE FunctionalDependencies #-}
{-# LANGUAGE TypeSynonymInstances #-}
{-# LANGUAGE LambdaCase #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE OverloadedStrings #-}

{-# OPTIONS_GHC -Wno-name-shadowing -Wno-unused-matches -Wno-unused-local-binds -Wno-unused-top-binds #-}

-- | Dot-Bid Auctions. A dot-bid auction is an auction where goods are sold at
-- integral quantities (often even unit quantities, although this is not a
-- requirement), and bidders are generally interested in a given quantity of
-- any of a set of alternative goods. That is, bidders want to express their
-- interest in buying N units of either good A at price P, or good B at price
-- Q, or good C at price R, etc., but only one of those options. There is a
-- fixed supply of each good available, and the algorithm will find the lowest
-- price at which this supply (or more) will be demanded completely.
--
-- Examples of goods for which this auction type is suitable include:
--
-- - FM radio frequencies: each bidder is typically only interested in one
--   frequency, but they will want to place bids on more than one of them,
--   because their preferred one might be so popular that they lose the
--   auction; but putting in two separate bids could have them end up with two
--   frequencies when they only need one.
-- - A government buying bad loans against a choice of collaterals. Here the
--   situation is reversed, as the auctioneer is *buying* loans; however, the
--   mechanism is exactly the same.
--
-- Throughout the code, we are referencing the paper
-- \"Finding a valid allocation\" (Baldwin, Goldberg, Klemperer 2017)
-- (sometimes referred to as \"validityAllocation\"); much
-- of the terminology used here is derived from that paper.
--
-- If you are not interested in the implementation, and just want to run an
-- auction, the following will be sufficient:
--
-- - The 'Auction' and 'AuctionResult' data structures
-- - The 'runAuction' / 'runAuctionM' functions
--
-- The REJECT good receives special treatment in the paper; it is a virtual
-- good that represents rejecting bids, that is, allocating on REJECT means
-- that we are rejecting bids, and finding a nonzero value in the REJECT
-- element of the resulting endowment means that not all the bidder's bids
-- have been honored.
--
-- Just like in the paper, we extend all inputs (bids, target bundle, reserve
-- price) with an additional REJECT coordinate, which we maintain throughout
-- the entire processing chain, to remove it again only at the very last step.
-- In other words, the problem is presented without the REJECT good on the
-- outside, but the REJECT coordinate is consistently present on the inside.
-- So, for an auction of 4 goods, we will see 4-dimensional vectors in the
-- input and output, but 5-dimensional vectors in all the internal functions.
--
-- The paper makes frequent use of sets; however, the operations performed on
-- them are often such that it is more efficient to represent them as vectors
-- where each component is either 0 or 1. Represented like this, set operations
-- become vector operations, e.g.:
--
-- @
-- set union -> component-wise maximum
-- set intersection -> component-wise minimum
-- set membership test -> component lookup + nonzero test
-- @
module ProductMixAuction.DotBids
(
-- * Basic Types
  Price (..)
, priceToFrac2
, unsafePriceFromFrac2
, priceFromFrac2
, Units (..)
, Good
, Bundle
, PriceVector
, emptyBundle
, isSubBundle
, priceOf
, addToPrice
, mulUnitPrice

-- * Helper functions for dealing with Vectors and Lists
, modifyAt
, generate
, extend
, reduce
, zero0
, dimension
, (!?)
, argmin, argmax
, opt
, exceedsUpperBound
, exceedsLowerBound
, unitWhereLower
, powersetL
, powerset
, eSet
, ePowerset
, eXs, eXs1, eX

-- * Bidders and Bid Agents
--
-- We are making a distinction between bidders and bid agents here.
--
-- A bidder is a real entity participating in the auction. For the inner
-- workings of the algorithms, we are however introducing an additional bid
-- agent representing the auctioneer; the intuition for this is that the
-- auctioneer will buy back any unsold goods up to the total supply (which we
-- express by adding suitable bids to the auction ourselves).
--
-- In order to model this, we allow the caller to pick a type for the bidder
-- (denoted as type variable @b@ throughout this module), which generally only
-- needs to have instances for 'Eq', 'Ord', and 'Show'; and we introduce the
-- 'BidAgentID' type, which can additionally represent the auctioneer.

, BidAgentID (..)

-- * Bids

, DotBidOf (..)
, dotBidQuantity
, dotBidPrices
, dotBidAgentID
, DotBid

-- * Calculating Demand

, maximumPrices
, bidsByAgents
, replaceAgentBids
, separateBidsAcrossGoodsMin
, separateBidsAcrossGoodsMax
, bundleDemandedMin
, bundleDemandedMax
, bundlesDemandedApprox
, setAuctioneerBids
, auctioneerBids
, auctioneerBidOn

-- * Calculating marginality
--
-- A bid is considered *marginal* between a set of goods G at a target price p
-- iff both the following are true:
--
-- - the difference between the bid price and the target price is maximal on
--   all goods in G.
-- - G has 2 or more elements.
--
-- And by extension, we say a bid is marginal on good g if there exists a set G
-- such that the bid is marginal on G, and g ∈ G.
--
-- Further yet, we say that a bid is marginal if it is marginal on any good g
-- (which also implies that it must be marginal on at least one other good).
--
-- The opposite, of course, is bids that are not marginal; we call these
-- non-marginal bids; a bid where the price difference between the bid price
-- and the target price is maximal on exactly one good is said to be
-- non-marginally accepted on that good. In theory, a bid can have all prices
-- below the target price, which would make it non-marginal but not accepted on
-- any goods; however, the REJECT good is always 0 for both the bids and the
-- target price, so if all other prices are below the target price, the bid
-- ends up being non-marginally accepted on the REJECT good. We call such bids
-- non-marginally rejected.

, hasNoMarginalBids
, marginalBidsOn
, marginalBids
, isMarginalOn
, nonMarginallyAcceptedBids
, nonMarginallyAcceptedBidsOn
, isNonMarginallyAcceptedOn
, anyMarginalOn
, biddersMarginalOn

-- * Finding an equilibrium price

-- ** Calculating g, the function to be minimised

, calcUG, calcU, calcG
, vpart

-- ** Minimizing g

, validTarget
, findBundleNaive
, findBundle
, findBundle_GreedyUpMinimal
, findBundle_GreedyUpMinimalM
, findBundle_TwoPhaseMinMin
, findBundle_TwoPhaseMinMinM
, findBundle_TwoPhaseMinMax
, findBundle_TwoPhaseMinMaxM
, findBundle_TwoPhaseM

, minimalMinimizer
, minimalMinimizerBy
, maximalMinimizer
, maximalMinimizerBy
, allPriceVectors
, upperBidBounds
, upperPriceBounds

, findIK
, minminCandidates'
, minminCandidates
, maxminCandidates'
, maxminCandidates

, jumpDistUp
, jumpDistDown

, minimalNeighbour
, minimalNeighbourBF
, minimalNeighbourFW
, minimalNeighbourBoth

-- * Helper functions for constructing bids, bundles, and price vectors.
--
-- Note that since price vector and bundle types have IsList instances, just
-- using 'fromList' or the @OverloadedLists@ extension makes these somewhat
-- superfluous.

, mkBid
, mkBid1
, mkPV1
, mkBundle1
, mkPosBid1
, mkNegBid1
, mkBid2
, mkPV2
, mkBundle2
, mkPosBid2
, mkNegBid2
, mkBid3
, mkPV3
, mkBundle3
, mkPosBid3
, mkNegBid3

-- * Generating sets of valid bids

, generateValidBids
, generateValidBidsW
, generateHorribleValidBidsW

-- * Allocating goods to bidders
-- ** Representing the Allocation Problem
, AllocationProblemOf (..)
, AllocationProblem

, apBids
, apResidualSupply
, apActiveGoods
, apPrice

-- ** Operations on Allocation Problems
, isVacuous
, totalGoods

-- ** Preparing Allocation Problems from bids and target bundles
, mkAllocationProblem
, mkAllocationProblemEx

-- ** Generating random Allocation Problems
, genAllocationProblem
, genHorribleAllocationProblem

-- ** Invariants
, allocationProblemTotalGoodsInvariant
, allocationProblemNoNegativeSupplyInvariant
, allocationProblemNoNegativeEndowmentsInvariant
, ppmMarginalityInvariant

, CheckInvariantException(..)
, check
, checkAllocationProblemInvariants

-- ** The Obvious Refinement
--
-- The Obvious Refinement takes care of all bids that can be trivially resolved
-- due to being non-marginally accepted or rejected. Any bid that is
-- non-marginally accepted can, by definition of the problem, be honored in
-- full, while any non-marginally rejected bid must be rejected in full. This
-- means that we can allocate the sum of quantities of all bids non-marginally
-- accepted on each good to their respective bidders on that good, and we can
-- remove all non-marginal bids from the allocation problem.
, allocateObvious
, allocateObviousL

-- ** Key Lists and Link Goods
--
-- A /Key List/ is the set of all goods that form a maximally connected
-- subgraph in the marginal bids graph from a given good and bidder.
--
-- A /Link Good/ is a good in the marginal bids graph that has edges for more
-- than one bidder; that is, there are bids marginal on this good from two or
-- more bidders.
--
-- Key lists and link goods are important, because we can unambiguously
-- allocate to all goods in key lists that contain zero or one link goods (see
-- below)
, findKeyList
, isLinkGood

-- ** Key List Allocation
, allocateByKeyList
, allocateByKeyListL

-- ** Jiggles
, jiggleBids
, jiggleBid
, jiggle
, jiggleL
, moveBid
, unjiggleBid

-- ** Projection Preserving Marginality
, ppm
, ppmBids
, ppmBid
, ppmPrice

-- ** Unraveling

, UnravelState (..)
, unravelUnambiguous
, unravelUnambiguousL

-- ** Priority-based allocation

, allocatePriorityBased
, allocatePriorityBasedL

-- * Running a full auction

-- ** Representing incoming auction data
, Auction (..)

, auctionReservePrice
, auctionSupply
, auctionBids

, IncomingBid (..)
, incomingBidQuantity
, incomingBidPrices

-- ** Representing outgoing auction result data

, AuctionResult (..)

, auctionResultPrices
, auctionResultEndowment
, auctionResultRemainingSupply

-- ** Running an auction
, runAuction
, runAuctionM

-- ** Converting to and from 'AllocationProblem's
, resultFromAP
, prepareAP
, prepareBids
, prepareAgentBids
, prepareAgentBid

-- ** Generating random 'Auction's
, genAuction
)
where

import Data.Default.Class
import Data.List
import Data.Ord
import qualified Data.Map.Strict as Map
import Data.Map.Strict (Map)
import qualified Data.Vector as V
import Data.Vector (Vector)
import Test.QuickCheck (Gen, shuffle, elements)
import Test.QuickCheck.Gen (choose)
import Control.Monad (forM_, forM, when, unless, foldM)
import Control.Monad.Identity (runIdentity)
import Data.Maybe (fromMaybe)
import Control.Monad.State (MonadState, StateT, execState, execStateT, put, get, runStateT)
import Control.Monad.Trans (lift)
-- import Debug.Trace as D
import GHC.Exts (IsList (..))
import qualified Data.Set as Set
import Data.Set (Set)
import Data.Monoid ( (<>) )
import Data.Bool (bool)
import Control.Lens ( use, over, zoom, (^.), (.=), (%=), (<>=), _1, _2, Lens )
import Control.Lens.TH (makeLenses)
import Control.Exception (throw, Exception (..))
import Data.Aeson (FromJSON, ToJSON, toJSON, parseJSON, FromJSONKey, ToJSONKey (..), (.:))
import Data.Aeson.TH (deriveJSON, fieldLabelModifier)
import qualified Data.Aeson as JSON
import qualified Data.Aeson.Encoding as JSON
import qualified Data.Csv as CSV
import GHC.Generics (Generic)

import ProductMixAuction.DotBids.Frac2
import ProductMixAuction.LinearAlgebra ( Dimension
                                       , (<+>)
                                       , (<->)
                                       , (<.>)
                                       , vmax
                                       , vmin
                                       , vscale
                                       , vlte
                                       , dotBy
                                       , zero
                                       , isZero
                                       , ComponentWise (..)
                                       )
import ProductMixAuction.FujishigeWolfe (findMin, fromBitMask, SubmodularFunction)
import System.Random


-- | Represents a price supplied by or presented to the user.
newtype Price = Price { _Price :: Integer }
  deriving (Enum, Eq, Integral, Ord, Random, Read, Real, Num, ToJSON, FromJSON, CSV.ToField, CSV.FromField, Generic)

instance Default Price

instance Show Price where
  show = show . _Price

priceToFrac2 :: Price -> Frac2
priceToFrac2 (Price x) = toFrac2 x

unsafePriceFromFrac2 :: Frac2 -> Price
unsafePriceFromFrac2 x = Price $ unsafeFromFrac2 x

priceFromFrac2 :: Frac2 -> Maybe Price
priceFromFrac2 x = Price <$> fromFrac2 x

-- | In dot-bids auctions, quantities are always integers (discrete
-- rather than continuous) but may be negative.
newtype Units = Units { _Units :: Integer }
  deriving (Eq, Num, Ord, Enum, Real, Integral, ToJSON, FromJSON, CSV.ToField, CSV.FromField, Generic)

instance Default Units where
  def = 1

instance Show Units where
  show (Units u) = show u

-- | Goods are numbered from 0 to N-1, where N is the dimension of the auction.
type Good = Int

-- | A bundle of goods (i.e. a quantity for each good).
type Bundle = Vector Units

-- | A vector in price space (i.e. a price for each good).
type PriceVector = Vector Price

-- | The bundle containing no goods.
emptyBundle :: Dimension -> Bundle
emptyBundle = zero

-- | Test whether the first bundle is a sub-bundle of the second
-- (i.e. it contains less of each good).
isSubBundle :: Bundle -> Bundle -> Bool
isSubBundle = vlte

-- | Find the price of a single good.
priceOf :: Good -> Vector a -> a
priceOf k p = p !? k

-- | Increment the price of a single good.
addToPrice :: (Num price) => Good -> price -> Vector price -> Vector price
addToPrice good p = modifyAt good (+ p)

-- | Modify the i-th component of a vector.
modifyAt :: Int -> (a -> a) -> V.Vector a -> V.Vector a
modifyAt i k v =
  -- TODO: this isn't very efficient with immutable vectors.
  v V.// [(i, k (v V.! i))]

mulUnitPrice :: (Num a) => Units -> a -> a
mulUnitPrice (Units u) p = fromIntegral u * p

-- | @generate n f@ generates a list-like data structure containing @n@ items,
-- such that item @n@ equals @f n@.
generate :: IsList l
         => Int
         -> (Int -> Item l)
         -> l
generate n f = fromList [ f n_ | n_ <- [0..n-1] ]

-- | Prepend a 0th element of value 0 to a vector.
-- This is effectively an implementation of the ι function described in
-- __Algorithm 1__, step 2.
extend :: Num a => Vector a -> Vector a
extend = V.cons 0

-- | Dual to 'extend': removes the first element from a vector.
reduce :: Vector a -> Vector a
reduce = V.drop 1

-- | Explicitly set the 0th element of a vector to 0. Several algorithms
-- require this step
zero0 :: Num a => Vector a -> Vector a
zero0 = (V.// [(0,0)])

-- | Get the dimension of a vector. Right now, it is just an alias for
-- 'length'.
dimension :: Vector a -> Int
dimension = length

-- | Get the n-th element from a vector, 0-based. This used to be important
-- when the codebase still had 1-based vectors with implicit 0th elements,
-- but now that this distinction has been removed, this is just an alias for
-- 'V.!'.
(!?) :: Vector a -> Int -> a
(!?) = (V.!)

argmin :: Ord b => (a -> b) -> [a] -> a
argmin p = minimumBy (comparing p)

argmax :: Ord b => (a -> b) -> [a] -> a
argmax p = maximumBy (comparing p)

-- | Identifies a bid agent. A bid agent can be either a real 'Bidder', or
-- the 'Auctioneer', who will buy all the unsold goods in the auction back at
-- any price.
data BidAgentID b
  = Auctioneer -- ^ The auctioneer, used to identify auctioneer bids.
  | Bidder b -- ^ A regular, "real" bidder.
  deriving (Eq, Ord, Show, Generic)

instance (ToJSON b, ToJSONKey b) => ToJSONKey (BidAgentID b) where
  toJSONKey =
    case JSON.toJSONKey of
      JSON.ToJSONKeyText f g ->
        JSON.ToJSONKeyText
          (\case
            Auctioneer -> ""
            Bidder b -> f b)
          (\case
            Auctioneer -> JSON.text ""
            Bidder b -> g b)
      JSON.ToJSONKeyValue f g ->
        JSON.ToJSONKeyValue
          (\case
            Auctioneer -> JSON.Null
            Bidder b -> f b)
          (\case
            Auctioneer -> JSON.null_
            Bidder b -> g b)

instance ToJSON b => ToJSON (BidAgentID b) where
  toJSON Auctioneer = JSON.Null
  toJSON (Bidder i) = toJSON i

instance (FromJSON b, FromJSONKey b) => FromJSONKey (BidAgentID b) where
  fromJSONKey = case JSON.fromJSONKey of
    JSON.FromJSONKeyText f ->
      JSON.FromJSONKeyText
        (\case
            "" -> Auctioneer
            s  -> Bidder (f s)
        )
    JSON.FromJSONKeyValue f ->
      JSON.FromJSONKeyValue
        (\case
            JSON.Null -> pure Auctioneer
            v         -> Bidder <$> f v
        )
    _ -> error "FromJSONKey BidAgentID: case not implemented"


instance FromJSON b => FromJSON (BidAgentID b) where
  parseJSON JSON.Null = pure Auctioneer
  parseJSON v = Bidder <$> parseJSON v

instance CSV.ToField b => CSV.ToField (BidAgentID b) where
  toField Auctioneer = ""
  toField (Bidder i) = CSV.toField i

instance CSV.FromField b => CSV.FromField (BidAgentID b) where
  parseField "" = pure Auctioneer
  parseField i = Bidder <$> CSV.parseField i

-- | A dot-bid is represented by a number of units (which may be
-- positive or negative) and a vector of prices for each good.
data DotBidOf b a =
  DotBid
    { _dotBidQuantity :: Units
    , _dotBidPrices   :: Vector a
    , _dotBidAgentID  :: BidAgentID b
    }
  deriving (Eq, Functor)

makeLenses ''DotBidOf
deriveJSON JSON.defaultOptions { fieldLabelModifier = drop 7 } ''DotBidOf

instance (CSV.ToField b, CSV.ToField a) => CSV.ToRecord (DotBidOf b a) where
  toRecord (DotBid q p a) =
    CSV.record $
      CSV.toField a :
      CSV.toField q :
      map CSV.toField (toList . reduce $ p)

type DotBid b = DotBidOf b Price

instance (Show b, Show a) => Show (DotBidOf b a) where
  show (DotBid q ps agent) =
    unlines
      [ "DotBid"
      , "{ _dotBidQuantity = " ++ show q
      , ", _dotBidPrices = " ++ show ps
      , ", _dotBidAgentID = " ++ show agent
      , "}"
      ]

-- | Calculate the maximum price bid on each good.
maximumPrices :: (Ord price, Num price)
              => Dimension -> [DotBidOf b price] -> Vector price
maximumPrices dim = foldl' (\ pv -> vmax pv . (^. dotBidPrices)) (zero dim)

-- | Sort a list of bids into a 'Map' indexed by bid agent. We will be using
-- this in algorithms that require us to only look at the set of bids from one
-- specific bid agent.
bidsByAgents :: forall a b
              . (Eq b, Ord b)
             => [DotBidOf b a] -> Map (BidAgentID b) [DotBidOf b a]
bidsByAgents bids =
  foldl insertBid Map.empty bids
  where
    insertBid :: Map (BidAgentID b) [DotBidOf b a] -> DotBidOf b a -> Map (BidAgentID b) [DotBidOf b a]
    insertBid m bid =
      Map.insertWith (++) (bid ^. dotBidAgentID) [bid] m

-- | @replaceAgentBids agentID bidsToKeep bids@ removes all bids for agent
-- @agentID@ from the list of @bids@, and then inserts all the bids in
-- @bidsToKeep@. The bids in @bidsToKeep@ should all be from the bidder
-- indicated by @agentID@, however this is not enforced.
replaceAgentBids :: Eq b
                 => BidAgentID b -> [DotBidOf b a] -> [DotBidOf b a] -> [DotBidOf b a]
replaceAgentBids agentID bidsToKeep =
  (++ bidsToKeep) . (filter ((/= agentID) . _dotBidAgentID))

-- | __Algorithm 2__ from validityAllocation: Separate bids across goods
-- at p, minimizing demand.
separateBidsAcrossGoodsMin :: [DotBid b] -> PriceVector -> [[DotBid b]]
separateBidsAcrossGoodsMin bids p =
  map
    (\i ->
      filter
        (\(DotBid { _dotBidPrices = v }) ->
          and
            [ (v !? i) - (p !? i) > (v !? k) - (p !? k)
            | k <- [0..pred i]
            ]
          &&
          and
            [ (v !? i) - (p !? i) >= (v !? k) - (p !? k)
            | k <- [0..dim-1]
            ]
        )
        bids
    )
    [0..dim-1]
  where
    dim = dimension p

-- | __Algorithm 3__ from validityAllocation: Separate bids across goods
-- at p, maximizing demand.
separateBidsAcrossGoodsMax :: [DotBid b] -> PriceVector -> [[DotBid b]]
separateBidsAcrossGoodsMax bids p =
  map
    (\i ->
      filter
        (\(DotBid { _dotBidPrices = v }) ->
          and
            [ (v !? i) - (p !? i) > (v !? k) - (p !? k)
            | k <- [succ i..dim-1]
            ]
          &&
          and
            [ (v !? i) - (p !? i) >= (v !? k) - (p !? k)
            | k <- [0..dim-1]
            ]
        )
        bids
    )
    [0..dim-1]
  where
    dim = dimension p

-- | __Algorithm 4__ from validityAllocation: Calculate a bundle demanded
-- at p, minimizing demand (via __Algorithm 2__ / 'separateBidsAcrossGoodsMin')
bundleDemandedMin :: [DotBid b] -> PriceVector -> Bundle
bundleDemandedMin bids p =
  fromList $
    map
      (sum . map _dotBidQuantity)
      bidsPerGood
  where
    bidsPerGood = separateBidsAcrossGoodsMin bids p

-- | __Algorithm 4__ from validityAllocation: Calculate a bundle demanded
-- at p, maximizing demand (via __Algorithm 2__ / 'separateBidsAcrossGoodsMax')
bundleDemandedMax :: [DotBid b] -> PriceVector -> Bundle
bundleDemandedMax bids p =
  fromList $
    map
      (sum . map _dotBidQuantity)
      bidsPerGood
  where
    bidsPerGood = separateBidsAcrossGoodsMax bids p

-- | Calculates an overestimation of bundles demanded at a given price point.
--
-- Not all bundles returned are guaranteed to be actually demanded at the price
-- point, but any bundle that is actually demanded is guaranteed to be in the
-- returned set.
--
-- The approximation is calculated based on marginality:
-- - bids that are non-marginally accepted at the price point are always
--   selected on the good that they are accepted on
-- - for each bid that is marginal on any goods (including REJECT), bundles are
--   explored that select the bid on each of the goods that it is marginal on
-- - any other bid (neither non-marginally accepted nor marginal) is ignored
bundlesDemandedApprox :: Eq b
                      => [DotBid b] -> PriceVector -> [Bundle]
bundlesDemandedApprox bids0 p = do
  go baseBundle0 mBids
  where
    dim = dimension p
    nmaBids = nonMarginallyAcceptedBids bids0 p
    mBids = nub . concat $ marginalBids p bids0
    baseBundle0 = fromList (map (sum . map _dotBidQuantity) nmaBids)

    go :: Bundle -> [DotBid b] -> [Bundle]
    go baseBundle [] = [baseBundle]
    go baseBundle (bid:bids) = do
      let diffBundles = (zero dim:) . concat $
            map
              (\i -> if isMarginalOn i p bid then [vscale (_dotBidQuantity bid) (eX dim i)] else [])
              [0..dim-1]
      diffBundle <- diffBundles
      go (baseBundle <+> diffBundle) bids

-- | Checks if the price vector has no marginal bids.
hasNoMarginalBids :: [DotBid b] -> PriceVector -> Bool
hasNoMarginalBids bids p = all null (marginalBids p bids)

-- | Get all bids marginal on a given good at a given price vector.
-- See __Algorithm 5__ in validityAllocation.pdf.
marginalBidsOn :: (Num a, Ord a) => Good -> Vector a -> [DotBidOf b a] -> [DotBidOf b a]
marginalBidsOn i p bids =
  filter (isMarginalOn i p) bids

-- | Get all bids marginal on any good at the given price vector.
marginalBids :: (Num a, Ord a) => Vector a -> [DotBidOf b a] -> [[DotBidOf b a]]
marginalBids p bids =
  map (\good -> marginalBidsOn good p bids) [0..dim-1]
  where
    dim = dimension p

-- | Check if a bid is marginal on a given good at a given price vector.
-- See __Algorithm 5__ in validityAllocation.pdf.
isMarginalOn :: (Num a, Ord a) => Good -> Vector a -> DotBidOf b a -> Bool
isMarginalOn i p bid =
  and
    [ (v !? i) - (p !? i) >= (v !? k) - (p !? k)
    | k <- [0..dim-1]
    ]
  &&
  or
    [ (v !? i) - (p !? i) == (v !? k) - (p !? k)
    | k <- [0..dim-1], k /= i
    ]
  where
    dim = dimension p
    v = bid ^. dotBidPrices
    -- w = bid ^. dotBidQuantity

-- | Calculate u and g, as shown in Lemma 3 in section 3.3 of "LP for dot-bits,
-- and usage of the Ellipsoid Method".
--
-- For a given set of bids and a target bundle, minimising g gives the price
-- vector at which the target bundle is demanded.
--
-- Returns the value of u in the first component, and the value of g in the
-- second component.
calcUG :: forall a b. (Num a, Ord a) => [DotBidOf b a] -> Vector Units -> Vector a -> (a, a)
calcUG bids t p =
  let
    n = length p
    vjs :: Vector (Int, [DotBidOf b a])
    vjs = fromList
            [ (j, vpart j p bids)
            | j <- [0 .. n - 1]
            ]
    x :: Vector Units
    x = V.map (\ (_, vj) -> sum $ map _dotBidQuantity vj) vjs
    u = sum
          [ mulUnitPrice (_dotBidQuantity b) (priceOf j (_dotBidPrices b))
          | (j, vj) <- V.toList vjs
          , b <- vj
          ]
    g = u - (dotBy mulUnitPrice (x <-> t) p)
  in
    (u, g)

-- | Calculate u only.
calcU :: (Num a, Ord a) => [DotBidOf b a] -> Vector Units -> Vector a -> a
calcU bids t p = fst (calcUG bids t p)

-- | Calculate g only.
calcG :: (Num a, Ord a) => [DotBidOf b a] -> Vector Units -> Vector a -> a
calcG bids t p = snd (calcUG bids t p)

-- | Given a good j and a price vector prices, this computes
-- the set @V_j(prices)@ from Definition 5 in Section 3.3 from
-- a list of dot bids.
--
-- A dot bid will end up in (at most) one of the V_j sets,
-- where j indicates the good which is demanded by this bid.
-- If multiple goods are equally close to the bid, then we
-- put the bid in the set for the good with the lowest index.

-- TODO: It seems computationally wasteful to compute all
-- the sets individually if we're actually talking about a
-- partitioning. A better algorithm would perhaps just go
-- through the list of dot bids once and put every one in
-- the correct set.
--
-- TODO: The document additionally defines
-- a rejection set for situations where the differences are
-- all 0 and lower. We don't, but allow 0-difference
-- coordinates.
--
-- Discussed with Elizabeth a different with an older
-- version, which had @vjpj >= 0@. The new version is
-- important if we use vpart in Lemma 2 for the computation
-- of the bundle, because it will give us the "better" bundle.
--
-- Problematic points are the ones where diagonal lines leave
-- to the upper right, because the "lower left" corner may
-- belong to more than one bundle.
--
vpart :: forall a b. (Num a, Ord a) => Good -> Vector a -> [DotBidOf b a] -> [DotBidOf b a]
vpart j prices = filter $ \ bid ->
  let
    n =  dimension prices
    v :: Good -> a
    v i = (bid ^. dotBidPrices) !? i
    p :: Good -> a
    p i = prices !? i
    vjpj = v j - p j
    vkpks =
      all
        (\ k -> let vkpk = v k - p k in vjpj >= vkpk && (vjpj /= vkpk || k >= j))
        [1 .. n - 1]
  in
    vjpj > 0 && vkpks

-- | Check that the given target bundle is a sub-bundle of the bundle
-- demanded at the origin.  If this is not the case, the search
-- algorithms may return an invalid result or fail to terminate.
validTarget :: [DotBid b] -> Bundle -> Bool
validTarget bids t =
    t `isSubBundle` bundleDemandedMin bids (zero $ dimension t)

-- | Find prices that cause the given bundle to be demanded.
--
-- This uses a brute-force search for the minimum of g across all
-- prices.
findBundleNaive :: [DotBid b] -> Bundle -> PriceVector
findBundleNaive bids t = argmin (calcG bids t) all_ps
  where
    all_ps = allPriceVectors (zero dim) (maximumPrices dim bids)
    dim = dimension t

-- | @minimalMinimizer p xs@ finds the smallest @x@ in @xs@ for which @p x@ is
-- minimal.
minimalMinimizer :: (Ord a, Ord b) => (a -> b) -> [a] -> a
minimalMinimizer = minimalMinimizerBy id

-- | @minimalMinimizerBy f p xs@ finds @x'@ in @xs'@ such that @xs'@ is the set of
-- @x@ for which @p x@ is minimal, and @f x'@ is minimal in @xs'@.
minimalMinimizerBy :: (Ord a, Ord b) => (c -> a) -> (c -> b) -> [c] -> c
minimalMinimizerBy f p =
    foldl1' (\a b -> if (p a, f a) < (p b, f b) then a else b)

-- | @maximalMinimizer p xs@ finds the smallest @x@ in @xs@ for which @p x@ is
-- maximal.
maximalMinimizer :: (Ord a, Ord b) => (a -> b) -> [a] -> a
maximalMinimizer = maximalMinimizerBy id

-- | @maximalMinimizerBy f p xs@ finds @x'@ in @xs'@ such that @xs'@ is the set of
-- @x@ for which @p x@ is minimal, and @f x'@ is maximal in @xs'@.
maximalMinimizerBy :: (Ord a, Ord b) => (c -> a) -> (c -> b) -> [c] -> c
maximalMinimizerBy f p =
    foldl1' (\a b -> if (p a, f b) < (p b, f a) then a else b)

-- | List all the prices contained in the hyperrectangle with given
-- bottom-left and top-right corners.
allPriceVectors :: PriceVector -> PriceVector -> [PriceVector]
allPriceVectors p1 p2 =
  -- TODO: surely we can do this more efficiently, even for the stupid
  -- search algorithm?
  go 1 (all_k 0 p1)
  where
    dim = dimension p1

    go k xs
      | k < dim   = go (k+1) (concatMap (all_k k) xs)
      | otherwise = xs

    all_k k p
      | p !? k <= p2 !? k =
          p : all_k k (addToPrice k 1 p)
      | otherwise =
          []


-- | Find prices that cause the given bundle to be demanded.
-- The 'validTarget' predicate should be satisfied to make sure this doesn't
-- loop.
findBundle :: (Eq b, Show b)
           => [DotBid b] -> Bundle -> PriceVector
findBundle = findBundle_GreedyUpMinimal

-- | Get the upper limit of bid prices per dimension, like 'upperPriceBounds',
-- but for 'DotBid's, not 'PriceVector's.
upperBidBounds :: Ord a => [DotBidOf b a] -> Maybe (Vector a)
upperBidBounds = upperPriceBounds . map _dotBidPrices

-- | Get the upper limit of prices per dimension. E.g.,
-- for bids @[ (1, 10), (5, 5), (9, 2) ]@, the upper limit is @(9, 10)@.
--
-- Empty input returns 'Nothing'.
upperPriceBounds :: Ord a => [Vector a] -> Maybe (Vector a)
upperPriceBounds [] =
    Nothing
upperPriceBounds xs =
    Just $ foldl1' vmax xs

exceedsUpperBound :: (Foldable f, ComponentWise f, Ord a) => f a -> f a -> Bool
x `exceedsUpperBound` b =
  or $ componentWise (>) x b

exceedsLowerBound :: (Foldable f, ComponentWise f, Ord a) => f a -> f a -> Bool
x `exceedsLowerBound` b =
  or $ componentWise (<) x b

-- | @unitWhereLower a b@ gets a vector v such that
-- @v[n] = 1@ if @a[n] < b[n]@, else @0@.
unitWhereLower :: (Ord a, ComponentWise f, Functor f, Num b) => f a -> f a -> f b
unitWhereLower a b =
  fmap (bool 1 0) $ componentWise (<) a b

-- | __Algorithm 10__ from validityAllocation.pdf: Identify bids non-marginally
-- accepted on each good (including REJECT)
nonMarginallyAcceptedBids :: [DotBid b] -> PriceVector -> [[DotBid b]]
nonMarginallyAcceptedBids bids p =
  map (\i -> nonMarginallyAcceptedBidsOn i bids p) [0..dim-1]
  where
    dim = dimension p

-- | Identify bids non-marginally accepted on a given good.
nonMarginallyAcceptedBidsOn :: (Num a, Ord a) => Good -> [DotBidOf b a] -> Vector a -> [DotBidOf b a]
nonMarginallyAcceptedBidsOn i bids p =
  filter (isNonMarginallyAcceptedOn i p) bids

-- | Check if a bid is non-marginally accepted on a given good.
isNonMarginallyAcceptedOn :: (Num a, Ord a) => Good -> Vector a -> DotBidOf b a -> Bool
isNonMarginallyAcceptedOn i p bid =
  and
    [ (v !? i) - (p !? i) > (v !? k) - (p !? k)
    | k <- [0..dim-1], k /= i
    ]
  where
    dim :: Int
    dim = dimension p
    v = bid ^. dotBidPrices

-- | Add auctioneer bids to an auction, replacing any existing auctioneer bids.
-- Needed for __Algorithm 1__.
setAuctioneerBids :: (Num a, Eq b) => Bundle -> [DotBidOf b a] -> [DotBidOf b a]
setAuctioneerBids t =
  replaceAgentBids Auctioneer (auctioneerBids t)

-- | Calculate auctioneer bids for a bundle.
-- Needed for __Algorithm 1__.
auctioneerBids :: (Num a) => Bundle -> [DotBidOf b a]
auctioneerBids t =
  filter ((/= 0) . (^. dotBidQuantity)) .
  map (\good -> auctioneerBidOn good t) $
    [1..dimension t-1]

-- | Calculate the auctioneer bid for a given good and bundle.
-- Needed for __Algorithm 1__.
auctioneerBidOn :: (Num a) => Good -> Bundle -> DotBidOf b a
auctioneerBidOn i t =
  mkBid w v Auctioneer
  where
    v = fromList [ if k == i || k == 0 then 0 else -1 | k <- [0..dim-1]]
    w = (t !? i) + 1
    dim = dimension t

-- | __Algorithm 11__ from validityAllocation.pdf: Identify I and K as in
-- Proposition 4.4.
-- Needed for Algorithms 12 and 13 (optimized minimal/maximal minimizers)
findIK :: (Num a, Ord a) => [DotBidOf b a] -> Vector a -> Bundle -> (Vector a, Vector a)
findIK bids p t =
  let (i_,k_) = unzip $ do
        i <- [0..dim-1]
        let nmaBids = nonMarginallyAcceptedBidsOn i bids p
            mBids = marginalBidsOn i p bids
            nmaSum = sum (map _dotBidQuantity nmaBids)
            mSum = sum (map _dotBidQuantity mBids)
        if nmaSum > t !? i
          then return (1, 0) -- add to I
          else if nmaSum + mSum <= t !? i
            then return (0, 1) -- add to K
            else return (0, 0) -- add to neither
  in (fromList i_, fromList k_)
  where
    dim = dimension p

-- | __Algorithm 12__ from the validityAllocation.pdf document.
-- This doesn't implement all of __Algorithm 12__, rather it just filters a given
-- list of candidates according to the conditions outlined in the
-- last section (such that...). Finding the actual minimal minimizer
-- is already handled in the crawler algorithm itself, and other
-- filters are applied there already, so we do not duplicate this
-- here.
minminCandidates' :: (Num a, Ord a) => [Vector a] -> [DotBidOf b a] -> Vector a -> Bundle -> [Vector a]
minminCandidates' candidates bids p t =
  filter useful candidates
  where
    (i, k) = findIK bids p t
    useful x =
      (reduce i <.> reduce x == reduce i) &&
      (isZero $ reduce k <.> reduce x)

minminCandidates :: (Num a, Ord a) => [DotBidOf b a] -> Vector a -> Bundle -> [Vector a]
minminCandidates bids p t =
  minminCandidates' candidates bids p t
  where
    -- We only need candidates with the 0th coordinate (REJECT) set to 0,
    -- because the inclusion condition never holds for x[0] /= 0.
    candidates = map extend $ eXs (dim - 1)
    dim = dimension p

-- | __Algorithm 13__ from the validityAllocation.pdf document.
-- This doesn't implement all of __Algorithm 13__, rather it just filters a given
-- list of candidates according to the conditions outlined in the
-- last section (such that...). Finding the actual minimal minimizer
-- is already handled in the crawler algorithm itself, and other
-- filters are applied there already, so we do not duplicate this
-- here.
maxminCandidates' :: (Num a, Ord a) => [Vector a] -> [DotBidOf b a] -> Vector a -> Bundle -> [Vector a]
maxminCandidates' candidates bids p t =
  filter useful candidates
  where
    (i, k) = findIK bids p t
    l = fmap (\case { 0 -> 0; _ -> 1 }) p
    kl = k <.> l
    useful x =
      (kl <.> x == kl) &&
      (isZero $ i <.> x)


maxminCandidates :: (Num a, Ord a) => [DotBidOf b a] -> Vector a -> Bundle -> [Vector a]
maxminCandidates bids p t =
  maxminCandidates' candidates bids p t
  where
    dim = dimension p

    -- Generated candidates need to include an explicit 0th coordinate
    -- Unlike minminCandidates, where we have to consider both x[0] = 1 and
    -- x[0] = 0, we skip the x[0] = 1 case here, because those would make
    -- p[0] - x[0] < 0, thus violating the p - x >= 0 requirement.
    candidates = eXs dim

-- | Find the powerset of a set, implemented using lists.
powersetL :: [a] -> [[a]]
powersetL [] = [[]]
powersetL (x:xs) = do
  ys <- powersetL xs
  [x:ys,ys]

-- | Find the powerset of a set.
powerset :: Ord a => Set a -> [Set a]
powerset base = map Set.fromList $ powersetL (Set.toList base)

-- | Create a unit vector from a set of coordinates.
eSet :: Num a => Good -> Set Good -> Vector a
eSet dim goods =
  fromList
    [ if good `Set.member` goods then 1 else 0
    | good <- [0..dim-1]
    ]

-- | Find the powerset of the set of N goods, represented as a list of vectors
-- such that for each good in the set, the corresponding vector's coordinate
-- is set to 1, while for each good not in the set, the vector's coordinate is
-- set to 0.
ePowerset :: Num a => Good -> Set Good -> [Vector a]
ePowerset dim goods = map (eSet dim) (powerset goods)

-- | Find prices that cause the given bundle to be demanded.
--
-- This implementation uses the GreedyUpMinimal algorithm described in
-- section 3 of "LP for dot-bids, and usage of the Ellipsoid Method", also
-- described as __Algorithm 8__ in validityAllocation.
findBundle_GreedyUpMinimal :: (Num a, Ord a, Show a, Integral a, Eq b, Show b)
                           => [DotBidOf b a]
                           -> Bundle
                           -> Vector a
findBundle_GreedyUpMinimal bids t =
  runIdentity $ findBundle_GreedyUpMinimalM
    (const $ return ())
    minimalNeighbour
    bids t

-- | The opt function frequently used in the validityAllocation paper: given
-- vectors p and v, find the goods on which p - v is maximal, and return
-- the maximal value and a vector representing the set of goods on which the
-- difference is maximal.
--
-- Informally, this gets us the best price relative to the target price.
opt :: (Num a, Ord a) => Vector a -> Vector a -> (a, Vector a)
opt p v =
  ( optVal
  , fmap (\x -> if x == optVal then 1 else 0) d
  )
  where
    optVal = V.maximum $ d
    d = (v <-> p)

-- | λ p x from the validityAllocation.pdf paper, used in __Algorithm 14__.
jumpDistUp :: forall a b. (Num a, Ord a) => [DotBidOf b a] -> Vector a -> Vector a -> a
jumpDistUp bids p x =
  execState go maxPrice
  where
    maxPrice :: a
    maxPrice = maximum (map (V.maximum . _dotBidPrices) bids)

    dim = dimension p
    go = forM_ bids $ \bid -> do
      let v = _dotBidPrices bid
          (a, optPV) = opt p v
      when (optPV <.> x == optPV) $ do
        forM_ [0..dim-1] $ \k -> do
          lambda <- get
          let lambda' = a - (v !? k - p !? k)
          when ((x !? k == 0) && (lambda' < lambda)) $
            put lambda'

-- | ν p x from the validityAllocation.pdf paper, used in __Algorithm 15__.
jumpDistDown :: (Num a, Ord a) => [DotBidOf b a] -> Vector a -> Vector a -> a
jumpDistDown bids p x =
  let minCandidates = fmap snd . filter ((/= 0) . fst) . toList $ componentWise (,) x p
  in if null minCandidates
    then 0
    else execState go (minimum minCandidates)
  where
    dim = dimension p
    go = do
      forM_ bids $ \bid -> do
        let v = _dotBidPrices bid
            (a, optPV) = opt p v
        when (isZero $ optPV <.> x) $ do
          forM_ [0..dim-1] $ \k -> do
            nu <- get
            let nu' = a - ((v !? k) - (p !? k))
            when ((x !? k /= 0) && (nu' < nu)) $
              put nu'

-- | Flavor of 'findBundle_GreedyUpMinimal' with tracing, in an arbitrary
-- 'Monad'.
findBundle_GreedyUpMinimalM ::
      forall a b m
    . (Monad m, Ord a, Num a, Show a, Eq b, Show b)
   => (String -> m ()) -- ^ trace function
   -> ([DotBidOf b a] -> Vector a -> Bundle -> Vector a) -- ^ find local minimum
   -> [DotBidOf b a] -- ^ bids
   -> Bundle -- ^ target bundle
   -> m (Vector a)
findBundle_GreedyUpMinimalM trace localMin bids' t = do
  -- Step 0: Find a vector p◦ ∈ dom g such that {p : p ≥ p◦ } contains some
  --         minimiser of g. In practice, this would be some price below all the
  --         bids. Set p := p◦.
  -- starting from Zero is good enough; alternatively, one could find the
  -- minimum of all bids along all dimensions, and start from there instead,
  -- however the jump optimization performed later on will amount to the same
  -- thing anyway.
  loop $ zero dim

  where
    bids :: [DotBidOf b a]
    bids = setAuctioneerBids t bids'

    dim :: Dimension
    dim = dimension t

    upperBound :: Vector a
    upperBound = fromMaybe (zero dim) (upperBidBounds bids)

    loop :: Vector a -> m (Vector a)
    loop p = do
      trace $ "Loop: " ++ show p
      trace $ "g = " ++ show (calcG bids t p)
      -- Step 1: Find the minimal minimizer X ⊆ [N] of g(p + e X ).
      let x1 = localMin bids p t
      trace $ "Found minimal minimizer eX: " ++ show x1
      let dist = jumpDistUp bids p x1
      trace $ "Jump by " ++ show dist
      let xscaled = vscale dist x1
      let x = if xscaled `exceedsUpperBound` (upperBound <-> p)
                then x1
                else xscaled
      let p' = p <+> x
      -- Step 2: If X = ∅, then output p and stop.
      if isZero x || p' `exceedsUpperBound` upperBound
          then do
            trace $ "Found solution: " ++ show p
            return p
          -- Step 3: Set p := p + eX and go to Step 1.
          else do
            loop p'

minimalNeighbour :: (Num a, Ord a, Integral a, Show a, Show b)
                 => [DotBidOf b a]
                 -> (Vector a)
                 -> Bundle
                 -> Vector a
minimalNeighbour = minimalNeighbourFW

minimalNeighbourBoth :: forall a b
                      . (Num a, Ord a, Integral a, Show a, Show b)
                     => [DotBidOf b a]
                     -> (Vector a)
                     -> Bundle
                     -> Vector a
minimalNeighbourBoth bids p t
  | mnFW == mnBF = mnFW
  | otherwise = error $ unlines
      [ "FW and brute force do not agree."
      , "Bids:"
      , unwords . lines $ show bids
      , "p: " ++ show p
      , "t: " ++ show t
      , "BF says: " ++ show mnBF ++ "(g = " ++ show (calcG bids t (p <+> mnBF)) ++ ")"
      , "FW says: " ++ show mnFW ++ "(g = " ++ show (calcG bids t (p <+> mnFW)) ++ ")"
      ]
  where
    mnFW = minimalNeighbourFW bids p t
    mnBF = minimalNeighbourBF bids p t

-- | Brute-force way of finding a minimal neighbour.
minimalNeighbourBF :: (Num a, Ord a)
                   => [DotBidOf b a]
                   -> Vector a
                   -> Bundle
                   -> Vector a
minimalNeighbourBF bids p t =
  let
    dim = dimension t
    candidates = minminCandidates bids p t
  in
    if null candidates then
      zero dim
    else
      minimalMinimizerBy sum (calcG bids t . (<+> p)) $ candidates

-- | Find minimal neighbour via Fujishige-Wolfe algorithm.
minimalNeighbourFW :: forall a b
                    . (Num a, Ord a, Integral a)
                   => [DotBidOf b a]
                   -> (Vector a)
                   -> Bundle
                   -> Vector a
minimalNeighbourFW bids p t =
  let
    dim = dimension t
    offset :: a
    offset = calcG bids t p
    sf :: SubmodularFunction
    sf v = fromIntegral (calcG bids t (p <+> fromBitMask v) - offset)
    p' = findMin dim sf
  in
    fromBitMask p'

-- | Find prices that cause the given bundle to be demanded.
--
-- This implementation uses the TwoPhaseMinMin algorithm described in
-- section 3 of "LP for dot-bids, and usage of the Ellipsoid Method". The
-- same algorithm is also described in validityAllocation as __Algorithm 9__.
findBundle_TwoPhaseMinMin :: (Eq b)
                          => PriceVector -> [DotBid b] -> Bundle -> PriceVector
findBundle_TwoPhaseMinMin p0 bids t =
  runIdentity $ findBundle_TwoPhaseMinMinM (const $ return ()) p0 bids t

-- This implementation uses the TwoPhaseMinMax algorithm described in
-- section 3 of "LP for dot-bids, and usage of the Ellipsoid Method".
findBundle_TwoPhaseMinMax :: (Eq b)
                          => PriceVector -> [DotBid b] -> Bundle -> PriceVector
findBundle_TwoPhaseMinMax p0 bids t =
  runIdentity $ findBundle_TwoPhaseMinMaxM (const $ return ()) p0 bids t

-- | Flavor of findBundle_TwoPhaseMinMin with tracing, in an arbitrary 'Monad'.
findBundle_TwoPhaseMinMinM :: forall m b
                            . (Eq b, Monad m)
                           => (String -> m ()) -> PriceVector -> [DotBid b] -> Bundle -> m PriceVector
findBundle_TwoPhaseMinMinM = findBundle_TwoPhaseM maximalMinimizerBy

-- | Flavor of findBundle_TwoPhaseMinMin with tracing, in an arbitrary 'Monad'.
findBundle_TwoPhaseMinMaxM :: forall m b
                            . (Eq b, Monad m)
                           => (String -> m ()) -> PriceVector -> [DotBid b] -> Bundle -> m PriceVector
findBundle_TwoPhaseMinMaxM = findBundle_TwoPhaseM minimalMinimizerBy

-- | Base implementation of the TwoPhase algorithms, with tracing.
findBundle_TwoPhaseM :: forall m b
                      . (Eq b, Monad m)
                     => ( (Vector Price -> Price)
                          -> (Vector Price -> Price)
                          -> [PriceVector]
                          -> Vector Price
                        )
                     -> (String -> m ())
                     -> PriceVector -> [DotBid b] -> Bundle -> m PriceVector
findBundle_TwoPhaseM minimizerDown trace p0 bids' t =
  -- Step 0: Find a vector p◦ ∈ dom g. Set p := p◦.
  -- (Note that in our implementation, we delegate the choice of p0 to the
  -- caller)
  loopUp $ vmin p0 upperBound

  where
    bids :: [DotBid b]
    bids = setAuctioneerBids t bids'

    dim :: Dimension
    dim = dimension t

    upperBound :: PriceVector
    upperBound = fromMaybe (zero dim) (upperBidBounds bids)

    lowerBound :: PriceVector
    lowerBound = zero dim

    -- Up Phase
    loopUp :: PriceVector -> m PriceVector
    loopUp p = do
      trace $ "-- Loop up: " ++ show p
      trace $ "g(p): " ++ show (calcG bids t p)

      -- Step 1: Find the minimal minimizer X ⊆ [N] of g(p + eX).
      let candidates =
            filter
              (\x -> not (x `exceedsUpperBound` (upperBound <-> p)))
              (minminCandidates bids p t)

      let x1 = minimalMinimizerBy sum (calcG bids t . (<+> p)) $ candidates
      trace $ "Found minimal minimizer eX: " ++ show x1
      let dist = jumpDistUp bids p x1
      trace $ "Jump by " ++ show dist
      let xscaled = vscale dist x1
      let x = if xscaled `exceedsUpperBound` (upperBound <-> p)
                then x1
                else xscaled
      let p' = p <+> x

      -- Step 2: If X = ∅, then go to the Down Phase
      if isZero x
          then do
            trace "Entering Down Phase"
            loopDown p
          -- Step 3: Set p := p + eX and go to Step 1.
          else do
            loopUp p'

    -- Down Phase
    loopDown :: PriceVector -> m PriceVector
    loopDown p = do
      trace $ "-- Loop down: " ++ show p
      trace $ "g(p): " ++ show (calcG bids t p)

      -- Step 1: Find the maximal minimizer X ⊆ [N] of g(p - eX).
      let candidates :: [PriceVector]
          candidates =
            map (zero0) $
            filter
              (\x -> not ((reduce $ p <-> x) `exceedsLowerBound` reduce lowerBound))
              (maxminCandidates bids p t)
      trace $ "maxminCandidates: " ++ show (maxminCandidates bids p t)
      trace $ "candidates: " ++ show candidates
      let x1 = if null candidates then
                 zero dim
               else
                 minimizerDown sum (calcG bids t . (p <->)) $ candidates
      trace $ "Found minimal minimizer eX: " ++ show x1
      let dist = jumpDistDown bids p x1
      trace $ "Jump by " ++ show dist
      let xscaled = vscale dist x1
      let x = if (p <-> xscaled) `exceedsLowerBound` lowerBound
                then x1
                else xscaled
      let p' = p <-> x

      -- Step 2: If X = ∅, then output p and stop.
      if isZero x
          then do
            trace $ "Found solution: " ++ show p
            return p
          -- Step 3: Set p := p - eX and go to Step 1.
          else do
            loopDown p'


-- | eX, from section 3 of "LP for dot-bids, and usage of the Ellipsoid
-- Method", used in several of the algorithms described there.
--
-- For X ⊂ [N] write eX for the vector defined by (eX)i = 1 iff i ∈ X, (eX)i =
-- 0 otherwise.
--
-- For efficiency reasons, instead of calculating individual eX for a given
-- subset X of [N], and considering the use cases, we implement instead a
-- function eXs, that enumerates *all* the possible eX for a given N.
eXs :: forall f a
     . ( IsList (f a)
       , Item (f a) ~ a
       , Num a
       )
    => Dimension -> [f a]
eXs dim = map fromList $ go dim
  where
    go :: Dimension -> [[a]]
    go 0 =
      [[]]
    go n = do
      h <- [0, 1]
      t <- go (pred n)
      return (h:t)

-- | A variation of eXs that does not include the zero vector.
eXs1 :: Dimension -> [PriceVector]
eXs1 = tail . eXs

eX :: ( IsList (f price)
      , Item (f price) ~ price
      , Num price
      )
   => Dimension -> Dimension -> f price
eX dim i = fromList [ if k == i then 1 else 0 | k <- [0..dim-1] ]

mkBid :: Units -> Vector a -> BidAgentID b -> DotBidOf b a
mkBid = DotBid

mkBid1 :: Units -> Price -> BidAgentID b -> DotBid b
mkBid1 q p = DotBid q (mkPV1 p)

mkPV1 :: Price -> PriceVector
mkPV1 p = extend $ V.singleton p

mkBundle1 :: Units -> Bundle
mkBundle1 u = extend $ V.singleton u

mkPosBid1 :: Price -> BidAgentID b -> DotBid b
mkPosBid1 = mkBid1 1

mkNegBid1 :: Price -> BidAgentID b -> DotBid b
mkNegBid1 = mkBid1 (-1)

mkBid2 :: Units -> (Price, Price) -> BidAgentID b -> DotBid b
mkBid2 q ps = DotBid q (mkPV2 ps)

mkPV2 :: (Price, Price) -> PriceVector
mkPV2 (p1, p2) = extend $ V.fromList [p1, p2]

mkBundle2 :: (Units, Units) -> Bundle
mkBundle2 (u1, u2) = extend $ V.fromList [u1, u2]

mkPosBid2 :: (Price, Price) -> BidAgentID b -> DotBid b
mkPosBid2 = mkBid2 1

mkNegBid2 :: (Price, Price) -> BidAgentID b -> DotBid b
mkNegBid2 = mkBid2 (-1)

mkBid3 :: Units -> (Price, Price, Price) -> BidAgentID b -> DotBid b
mkBid3 q ps = DotBid q (mkPV3 ps)

mkPV3 :: (Price, Price, Price) -> PriceVector
mkPV3 (p1, p2, p3) = extend $ V.fromList [p1, p2, p3]

mkBundle3 :: (Units, Units, Units) -> Bundle
mkBundle3 (u1, u2, u3) = extend $ V.fromList [u1, u2, u3]

mkPosBid3 :: (Price, Price, Price) -> BidAgentID b -> DotBid b
mkPosBid3 = mkBid3 1

mkNegBid3 :: (Price, Price, Price) -> BidAgentID b -> DotBid b
mkNegBid3 = mkBid3 (-1)

-- | Monadically generate a vector with 1-based indexing.
genM :: Monad m => Int -> (Int -> m a) -> m (V.Vector a)
genM size fun = V.generateM size (\ i -> fun (i + 1))

-- | Generate a vector with 1-based indexing.
gen :: Int -> (Int -> a) -> V.Vector a
gen size fun = V.generate size (\ i -> fun (i + 1))

-- | Lookup a 1-based index in a vector.
(!) :: V.Vector a -> Int -> a
v ! i = v V.! (i - 1)

-- | The call @sel b k@ yields all possibilities to select
-- @k@ elements out of @[1 .. b]@.
--
sel :: Int -> Int -> [[Int]]
sel b k = selhelper k [1 .. b]
  where
    selhelper :: Int -> [Int] -> [[Int]]
    selhelper 0 _  = [[]]
    selhelper n xs = [ (i : is) | (i, ys) <- neTails xs, is <- selhelper (n - 1) ys ]

    neTails :: [a] -> [(a, [a])]
    neTails []       = []
    neTails (x : xs) = (x, xs) : neTails xs

-- | This implements the previous version of the valid bids generation
-- procedure as shown in section 3.1 in "Valid Combinations of Bids", where
-- all bid weights are 1. For a generalized version as described in a newer
-- version of the same paper, which produces randomly weighted bids, see
-- 'generateValidBidsW'.
generateValidBids :: Integer -> Int -> Int -> BidAgentID b -> Gen ([PriceVector], [DotBid b])
generateValidBids = generateValidBidsW 1

-- | This implements the valid bids generation procedure as shown in section
-- 3.1 in "Valid Combinations of Bids". This flavor produces randomly weighted
-- bids; set @maxWeight@ to 1 to disable this behavior, or use
-- 'generateValidBids' instead.
generateValidBidsW :: Integer -> Integer -> Int -> Int -> BidAgentID b -> Gen ([PriceVector], [DotBid b])
generateValidBidsW maxWeight maxPrice n b agent = do -- b is chosen in Step 1, but we take it as an input
  w <- choose (1, maxWeight)
  a <- genM n (const $ Price <$> choose (1, maxPrice))
  a0 <- Price <$> choose (1, maxPrice)
  c <- genM n (const $ Price <$> choose (1, maxPrice)) -- Step 2
  generateValidBidsW' w a0 a c n b agent

-- | A variation on the 'generateValidBidsW' function that intentionally
-- generates bids sets that \"line up with each other in a horrible way\".
generateHorribleValidBidsW :: Integer -> Integer -> Int -> BidAgentID b -> Gen ([PriceVector], [DotBid b])
generateHorribleValidBidsW maxWeight _maxPrice n agent = do -- b is chosen in Step 1, but we take it as an input
  w <- choose (1, maxWeight)
  let a = gen n (const 1)
      a0 = 1
      c = gen n (const 1)
  b <- choose (n-2, n)
  generateValidBidsW' w a0 a c n b agent

generateValidBidsW' ::
  Integer -> Price -> PriceVector -> PriceVector -> Int -> Int
  -> BidAgentID b
  -> Gen ([PriceVector], [DotBid b])
generateValidBidsW' w a0 a c n b agent = do
  sigmaMap <- V.fromList <$> shuffle [1 .. n]
  let sigma j = sigmaMap ! j -- Step 3
  -- Steps 4 and 5, vs corresponds to V
  vs <-
    genM b $ \ i ->
    genM n $ \ j ->
    if      sigma j == i then pure (a ! j) -- Step 4 b
    else if sigma j <= b then pure 0       -- Step 4 a
    else    elements [0, a ! j]            -- Step 4 c
  -- Step 6 a
  let js i1 i2 = [ j | j <- [1 .. n], vs ! i1 ! j /= vs ! i2 ! j ] -- contains the indices where v_i1 and v_i2 differ
  let e is = gen n (\ i -> if i `elem` is then 1 else 0)
  -- Step 6 b
  let x i1 i2 = vmax (vs ! i1) (vs ! i2) <+> vscale a0 (e (js i1 i2))
  -- Step 8
  let mixed =
        map
          (\ (k, is) ->
            DotBid
              (Units ((-1) ^ (k - 1)))
              (extend $ foldl' vmax (gen n (const 0)) (map (vs !) is))
              agent
          )
          (concatMap (\ k -> map ((,) k) (sel b k)) [1 .. b] :: [(Int, [Int])])
  -- Step 9
  let xs = map (\ [i1, i2] -> DotBid (Units 1) (extend $ x i1 i2) agent) (sel b 2)
  -- Step 10
  let bids =
        map
          (\ (DotBid (Units q) v a_) -> DotBid (Units $ q * w) (v <+> extend c) a_)
          (mixed ++ xs)
  pure (map extend $ toList vs, bids)

-- | An allocation problem is defined as:
--
-- 1. An equilibrium price vector, found using one of the 'findBundle'
--    functions (named @p@ in the paper).
-- 2. Bids, including auctioneer bids (named @V@ in the paper).
-- 3. Bid weights (implemented as the '_dotBidQuantity' field of the 'DotBidOf'
--    data structure; named @w@ in the paper).
-- 4. Such that for each bidder (@j ∈ J@ in the paper), the set of bids for
--    that bidder is P-valid.
-- 5. An endowment map, associating each bidder with a bundle allocated to that
--    bidder (@m@ in the paper).
-- 6. A residual supply (@r@ in the paper), initialized to the target bundle
--    (@t@ in the paper) but, over the course of resolving the problem, to be
--    reduced to all-zero if possible.
-- 7. Additionally, we track a set of active goods, that is, all the goods that
--    haven't been resolved fully.
--
-- Allocation problems are polymorphic on a bidder ID type, @b@, and a price
-- type, @a@. The latter is necessary because some parts of the algorithm
-- require fractional prices (specifically, half-unit prices).
data AllocationProblemOf b a =
  AllocationProblem
    { _apPrice :: Vector a
    , _apBids :: [DotBidOf b a]
    , _apEndowment :: Map (BidAgentID b) Bundle
    , _apResidualSupply :: Bundle
    , _apActiveGoods :: Set Good
    }
    deriving (Show, Eq, Functor)

makeLenses ''AllocationProblemOf

type AllocationProblem b = AllocationProblemOf b Price

-- | An allocation problem is considered vacuous if any of the following hold
-- true:
--
-- 1. there are no more bids to consider
-- 2. there is no residual supply left to allocate
-- 3. there are no more active goods in the active goods set
--
-- When an allocation problem is vacuous, we can assume that no further
-- processing is desired.
isVacuous :: AllocationProblemOf b a -> Bool
isVacuous ap =
  null (ap ^. apBids) || isZero (ap ^. apResidualSupply) || null (ap ^. apActiveGoods)

-- | Create an initial allocation problem from a set of bids and a target
-- bundle (initial supply).
mkAllocationProblem :: (Eq b, Show b)
                    => [DotBid b] -> Bundle -> AllocationProblem b
mkAllocationProblem bids t =
  runIdentity $
    mkAllocationProblemEx
      (const $ return ())
      (zero dim)
      (\trace -> findBundle_GreedyUpMinimalM trace minimalNeighbour)
      bids
      t
  where
    dim = dimension t

-- | Fully configurable preparation of an allocation problem.
mkAllocationProblemEx :: (Eq b, Monad m)
                      => (String -> m ()) -- ^ Trace logger
                      -> PriceVector -- ^ Auctioneer reserve price
                      -> ((String -> m ()) -> [DotBid b] -> Bundle -> m PriceVector) -- ^ Equilibrium finder
                      -> [DotBid b] -- ^ Bids
                      -> Bundle -- ^ Available supply (target bundle)
                      -> m (AllocationProblem b)
mkAllocationProblemEx trace pMin findEquilibrium bids t = do
  p <- findEquilibrium trace bids' t
  return $ AllocationProblem
    p
    bids'
    Map.empty
    r
    (Set.fromList [0..dim-1])
  where
    -- Move all bids such that the auctioneer's minimum price point is mapped
    -- to the origin
    bids' = setAuctioneerBids t $ map (over dotBidPrices (<-> pMin)) bids
    r = t V.// [(0, r0)]
    r0 = (sum $ map _dotBidQuantity bids') - (V.sum . reduce $ t)
    dim = dimension t

-- | Generate a valid allocation problem with multiple bidders.
genAllocationProblem :: Integer -> Integer -> Int -> Int -> Int -> Gen (AllocationProblem Int)
genAllocationProblem maxWeight maxPrice n b maxNumAgents = do
  bids <- mconcat
          . fmap snd
          <$> forM
                [1..maxNumAgents]
                (generateValidBidsW maxWeight maxPrice n b . Bidder)
  bundle <- extend . fromList <$> forM [1..n] (const $ Units . fromIntegral <$> choose (0, maxWeight))
  return $ mkAllocationProblem bids bundle

-- | Generate a \"horrible\" (but valid) allocation problem with multiple
-- bidders.
genHorribleAllocationProblem :: Integer -> Integer -> Int -> Int -> Gen (AllocationProblem Int)
genHorribleAllocationProblem maxWeight maxPrice n maxNumAgents = do
  bids <- mconcat
          . fmap snd
          <$> forM
                [1..maxNumAgents]
                (generateHorribleValidBidsW maxWeight maxPrice n . Bidder)
  bundle <- extend . fromList <$> forM [1..n] (const $ Units . fromIntegral <$> choose (0, maxWeight))
  return $ mkAllocationProblem bids bundle

-- | Apply the Obvious Refinement defined in section 6.2 of the
-- validityAllocation paper.
--
-- Loosely following __Algorithm 16__.
allocateObvious :: forall a b m. (Num a, Ord a, Monad m, Show a, Ord b, Show b)
                => StateT (AllocationProblemOf b a) m ()
allocateObvious = allocateObviousL (const $ return ())

-- | Flavor of 'allocateObvious' with monadic trace logging.
allocateObviousL :: forall a b m. (Num a, Ord a, Show a, Ord b, Show b, Monad m)
                 => (String -> m ()) -> StateT (AllocationProblemOf b a) m ()
allocateObviousL trace' = do
  get >>= trace . show
  agentBids <- bidsByAgents <$> use apBids
  activeGoods <- use apActiveGoods
  -- dim <- dimension <$> use apPrice

  -- Clear the ActiveGoods list before going into the per-agent processing.
  -- This is required because in the final result, the list of active goods
  -- must contain exactly those goods for which there are any marginal goods,
  -- however, we don't want to do the marginality calculation twice, so we
  -- need to do it per agent, and then accumulate the results across
  -- per-agent iterations. In order to do that, we buffer the original set of
  -- active goods, clear the active goods set in the allocation problem we
  -- pass around, and mappend any goods we find to still be active as we
  -- encounter them.
  apActiveGoods .= Set.empty
  mapM_ (processAgent activeGoods) $ Map.toList agentBids

  where
    trace = lift . trace'

    -- | Apply the Obvious Refinement for one agent (j ∈ J)
    -- processAgent :: Set Good -> (BidAgentID, [DotBid b]) -> State (AllocationProblemOf b a) ()
    processAgent :: Set Good -> (BidAgentID b, [DotBidOf b a]) -> StateT (AllocationProblemOf b a) m ()
    -- processAgent _ (Auctioneer, _) = return ()
    processAgent activeGoods (agentID, bids) = do
      ap0 <- get
      trace $ show agentID
      p <- use apPrice
      let dim = dimension p
      let nmaBids = [ if good `Set.member` activeGoods then
                        nonMarginallyAcceptedBidsOn good bids p
                      else
                        []
                    | good <- [0..dim-1]
                    ]
          -- The amount we want to allocate for each good is the sum of the
          -- weights of all the bids non-marginally accepted on that good for
          -- the current bidder.
          bundleToAllocate = fromList $
            map (sum . map _dotBidQuantity) nmaBids
          -- We only need to keep the bids for this bidder that are marginal on
          -- any good: the other bids are either non-marginally accepted (which
          -- means we are allocating them in this step), or non-marginally
          -- rejected (which means we can just throw them away).
          (activeGoods', mBids) = flip execState (Set.empty, []) $ do
            forM_ bids $ \bid -> do
              marginal <- foldM
                (\marginal good -> do
                  if isMarginalOn good p bid then do
                    _1 %= Set.insert good
                    return True
                  else
                    return marginal
                )
                False
                (Set.toList activeGoods)
              when marginal $ _2 %= (bid:)

      trace $ "p: " ++ show p
      trace $ "NMA: " ++ show nmaBids
      trace $ "M: " ++ show mBids
      trace $ "To Allocate: " ++ show bundleToAllocate
      use apResidualSupply >>= trace . ("Residual supply before: " ++) . show
      apBids %= replaceAgentBids agentID mBids
      apEndowment %= Map.insertWith (<+>) agentID bundleToAllocate
      apResidualSupply %= (<-> bundleToAllocate)
      apActiveGoods <>= activeGoods'
      use apResidualSupply >>= trace . ("Residual supply after: " ++) . show
      get >>= checkAllocationProblemInvariants ap0

allocationProblemTotalGoodsInvariant :: AllocationProblemOf b a -> AllocationProblemOf b a -> Bool
allocationProblemTotalGoodsInvariant ap1 ap2 =
  totalGoods ap1 == totalGoods ap2

allocationProblemNoNegativeSupplyInvariant :: (Num a, Ord a) => AllocationProblemOf b a -> Bool
allocationProblemNoNegativeSupplyInvariant ap =
  all (>= 0) (ap ^. apResidualSupply)

allocationProblemNoNegativeEndowmentsInvariant :: (Num a, Ord a) => AllocationProblemOf b a -> Bool
allocationProblemNoNegativeEndowmentsInvariant ap =
  all (all (>= 0)) (ap ^. apEndowment)

ppmMarginalityInvariant :: AllocationProblem b -> AllocationProblem b -> Bool
ppmMarginalityInvariant ap0 ap1 =
  counts ap0 == counts ap1
  where
    counts :: AllocationProblem b -> (Int, Int)
    counts ap =
      let p = ap ^. apPrice
          bids = ap ^. apBids
          numMarginal = length . concat $ marginalBids p bids
          numNonMarginal = length . concat $ nonMarginallyAcceptedBids bids p
      in
        (numMarginal, numNonMarginal)

data CheckInvariantException = CheckInvariantException String
  deriving (Eq, Show)

instance Exception CheckInvariantException

check :: Monad m => String -> Bool -> m ()
check msg False = throw $ CheckInvariantException msg
check _ True = return ()

checkAllocationProblemInvariants :: (Num a, Ord a, Monad m)
                                  => AllocationProblemOf b a
                                  -> AllocationProblemOf b a
                                  -> m ()
checkAllocationProblemInvariants pre post = do
  check "total goods constant" $
    allocationProblemTotalGoodsInvariant pre post
  check "no negative supply" $
    allocationProblemNoNegativeSupplyInvariant post
  check "no negative endowments" $
    allocationProblemNoNegativeEndowmentsInvariant post

-- | The total number of units per good in a given allocation problem, i.e.,
-- the sum of all allocated goods plus the residual supply.
totalGoods :: AllocationProblemOf b a -> Bundle
totalGoods ap =
  foldl' (<+>) (ap ^. apResidualSupply) (Map.elems (ap ^. apEndowment))

-- | Determine whether a good is a link good at the specified price point, as
-- per Definition 6.6 and __Algorithm 17__ in the validityAllocation paper.
--
-- A good is considered a link good iff the number of distinct bidders marginal
-- on this good at the given price point exceeds 1, or, in terms of the
-- Marginal Goods Graph, iff it has edges labelled by more than one bidder.
isLinkGood :: (Num a, Ord a, Ord b) => [DotBidOf b a] -> Vector a -> Good -> Bool
isLinkGood bids p good =
  Set.size bidders >= 2
  where
    bidders = biddersMarginalOn bids p good

-- | Check if any bids are marginal on the given good.
anyMarginalOn :: (Num a, Ord a) => Good -> Vector a -> [DotBidOf b a] -> Bool
anyMarginalOn good p bids =
  any (isMarginalOn good p) bids

-- | List all unique bidders that are marginal on a given good at a given
-- price point.
biddersMarginalOn :: (Num a, Ord a, Ord b) => [DotBidOf b a] -> Vector a -> Good -> Set (BidAgentID b)
biddersMarginalOn bids p good =
  Set.fromList [ agent | (agent, abids) <- agentBids, anyMarginalOn good p abids ]
  where
    agentBids = toList $ bidsByAgents bids

-- | Find a key list within a set of bids, starting at a given good that is
-- assumed to be marginal for at least one bid at the given price vector p.
--
-- The given bids are assumed to all belong to the same bidder.
--
-- A key list is a set of at least two goods and exactly one bidder, such that
-- the goods form the nodes of a maximally connected subgraph of the marginal
-- bids graph for the given bids and price vector.
--
-- Unlike __Algorithm 18__ described in validityAllocation, we return only the
-- goods that make up a key list, not the bidder, which is implicitly assumed
-- to be the only bidder in the list of given bids.
--
-- Also unlike __Algorithm 18__, we return a singleton set if the given good is not
-- part of a key list; consumers who wish to treat this case specially will
-- have to perform the relevant checks themselves.
findKeyList :: (Num a, Ord a) => [DotBidOf b a] -> Vector a -> Good -> Set Good
findKeyList bids p good =
  -- We are not following __Algorithm 18__ to the letter; particularly, rather than
  -- keeping track of the marginal bids to consider, we track two sets of goods.
  --
  -- The first set of goods is the "to do list": these are the goods we have
  -- found in previous steps, so we know they are part of the key list, but we
  -- haven't expanded further from them, so we do not know yet whether they have
  -- any neighbors that are also part of the key list.
  --
  -- The second set of goods is the tentative output: this set contains all the
  -- goods that we have already visited, so we know that 1) they are in the key
  -- list, 2) all their neighbors that are also in the key list are either in the
  -- "to do" set or in the "output" set themselves.
  --
  -- With these two sets in hand, our algorithm's iteration becomes simple:
  -- - If the todo list is empty, terminate and return the output set.
  --   Otherwise:
  -- - For each pair of (i: good in the todo list, j: good in neither the todo
  --   list nor the output list), check if there are any bids that are marginal
  --   on both i and j.
  -- - If there are, then add j to the new todo list.
  -- - Move i from the todo list to the output list.
  -- - Repeat.
  go (Set.singleton good) Set.empty
  where
    dim = dimension p
    go :: Set Good -> Set Good -> Set Good
    go goods output = case Set.toList goods of
        [] ->
          output
        (x:_xs) ->
          let
            -- Find all the goods y for which...
            found = Set.fromList
              [ y
              | y <- [0..dim-1]
              -- ...y is not the good we're currently looking at,
              , y /= x
              -- ...y is not in the todo set,
              , not (y `Set.member` goods)
              -- ...y is not in the output yet,
              , not (y `Set.member` output)
              -- ...and a bid exists such that...
              , bid <- bids
              -- ...that bid is marginal on the good we're currently looking at...
              , isMarginalOn x p bid
              -- ...and also marginal on the good y we picked
              , isMarginalOn y p bid
              ]
          in
            go (Set.delete x goods <> found) (Set.insert x output)

-- | Allocate to key list goods, if possible, as described in __Algorithm 19__ in
-- validityAllocation.
allocateByKeyList :: (Num a, Ord a, Show a, Ord b, Show b, Monad m)
                  => BidAgentID b -> Set Good -> StateT (AllocationProblemOf b a) m ()
allocateByKeyList =
  allocateByKeyListL (const $ return ())

-- | 'allocateByKeyList' with tracing in an arbitrary 'Monad'.
allocateByKeyListL :: forall a b m
                    . (Num a, Ord a, Show a, Ord b, Show b, Monad m)
                   => (String -> m ())
                   -> BidAgentID b -> Set Good -> StateT (AllocationProblemOf b a) m ()
allocateByKeyListL trace' agent keyGoods = do
  bids <- use apBids
  p <- use apPrice
  let dim = dimension p
  let linkGoods = Set.filter (isLinkGood bids p) keyGoods
  trace $ "Link goods: " ++ show linkGoods
  -- More than 1 link good means we cannot easily decide this, jiggling is
  -- required.
  unless (Set.size linkGoods >= 2) $ do
    -- Proceed.
    rr <- use apResidualSupply
    let r i = rr !? i
        r' = rr <-> deltaM

        agentBids = filter belongsToAgent bids

        -- All the bids from this bidder that are marginal on at least one good
        -- from the key list, as well as those that are non-marginally accepted
        -- on a good from the key list.
        acceptableBids =
          filter
            isAcceptable
            agentBids

        unacceptableBids =
          filter
            (not . isAcceptable)
            agentBids

        remainingBids = replaceAgentBids agent unacceptableBids bids

        belongsToAgent bid = bid ^. dotBidAgentID == agent

        isAcceptable bid =
          or [ isMarginalOn good p bid || isNonMarginallyAcceptedOn good p bid
             | good <- toList keyGoods
             ]

        deltaM = fromList $ map (\i ->
            if not (i `Set.member` keyGoods) then
              0
            else if i `Set.member` linkGoods then
              sum (map _dotBidQuantity acceptableBids) -
              sum [ if ii == i then 0 else r ii | ii <- toList keyGoods ]
            else
              r i
          )
          [0..dim-1]

        nonLinkKeyGoods = keyGoods `Set.difference` linkGoods

    trace $ "Residual supply: " ++ show rr
    trace $ "Key goods: " ++ show keyGoods
    trace $ "Link goods: " ++ show linkGoods
    trace $ "Acceptable agent bids: " ++ show acceptableBids
    trace $ "Sum of acceptable bid weights: " ++ show (sum (map _dotBidQuantity acceptableBids))
    trace $ "Total residual supply on non-link key goods: " ++
      (show . sum . map r . toList $ nonLinkKeyGoods)
    trace $ "Total residual supply on key goods: " ++
      (show . sum . map r . toList $ keyGoods)
    trace $ "Total supply allocated: " ++
      (show . V.sum $ deltaM)
    trace $ "Allocation to " ++ show agent ++ ": " ++ show deltaM
    apResidualSupply .= r'
    apEndowment %=
      Map.insertWith
        (<+>)
        agent
        deltaM
    apBids .= remainingBids
    apActiveGoods %= (`Set.difference` nonLinkKeyGoods)
  where
    trace :: forall s. String -> StateT s m ()
    trace = lift . trace'


-- | A jiggle, as per definition 7.2 in the validityAllocation paper.
jiggleBids :: (Eq b)
           => BidAgentID b -> Good -> [DotBidOf b Frac2] -> [DotBidOf b Frac2]
jiggleBids agentID good bids =
  [ if bid ^. dotBidAgentID == agentID then
      jiggleBid good bid
    else
      bid
  | bid <- bids
  ]

-- | Jiggle an individual bid.
jiggleBid :: Good -> DotBidOf b Frac2 -> DotBidOf b Frac2
jiggleBid good bid =
  over dotBidPrices (<+> vscale 0.5 (eX dim good)) bid
  where
    dim = dimension (bid ^. dotBidPrices)

-- | A single (i*, j*) jiggle, as described in __Algorithm 20__ in
-- validityAllocation.
jiggle :: (Monad m, Show b, Eq b)
       => Good -> BidAgentID b -> StateT (AllocationProblemOf b Frac2) m ()
jiggle = jiggleL (const $ return ())

-- | 'jiggle' with tracing in an arbitrary 'Monad'.
jiggleL :: (Monad m, Show b, Eq b)
        => (String -> m ()) -> Good -> BidAgentID b -> StateT (AllocationProblemOf b Frac2) m ()
jiggleL trace' good agent = do
  bids <- use apBids
  p <- use apPrice
  r <- use apResidualSupply
  activeGoods <- use apActiveGoods
  let dim = dimension p
      -- 1.-3. Set Vj*' and w'
      bids' = jiggleBids agent good bids

      -- 4. Use __Algorithm 11__ to find I and K
      (i,k) = findIK bids' p r
      -- 5. Be ready to calculate g given Vj',w' and r
      g' = calcG bids' r
      g0 = g' p
      exsMin = map (vscale 0.5)
        . filter (\e -> (e !? good) /= 0)
        . ePowerset dim
        . Set.insert good
        $ activeGoods
      -- 6. Find the minimal minimizer of g(p + ½eX)
      candidatesMin = exsMin
      -- candidatesMin =
      --   minminCandidates'
      --     exsMin
      --     bids'
      --     p
      --     r
      xMin =
        if null candidatesMin then
          zero dim
        else
          minimalMinimizerBy sum (g' . (p <+>)) candidatesMin
      pMin = p <+> xMin
      gMin = g' pMin
      -- 7. Find the maximal minimizer of g(p - ½eX)
      exsMax = map (vscale 0.5)
        . ePowerset dim
        . Set.delete good
        $ activeGoods
      candidatesMax = exsMax
      -- candidatesMax =
      --   maxminCandidates'
      --     exsMax
      --     bids'
      --     p
      --     r
      xMax =
        if null candidatesMax then
          zero dim
        else
          maximalMinimizerBy sum (g' . (p <->)) candidatesMax
      pMax = p <-> xMax
      gMax = g' pMax
      -- 6.-8.
      -- If g(p+½eX) < g(p) then p' = p+½eX and terminate;
      -- If g(p-½eX) < g(p) then p' = p-½eX and terminate;
      -- Else keep p.
      p' = if gMin < g0 then
             pMin
           else if gMax < g0 then
             pMax
           else
             p

  trace $ "jiggle candidates at " ++ show p ++ " for good " ++ show good
  trace $ "p = " ++ show p ++ ": g(p) = " ++ show g0
  trace $ "pMin = " ++ show pMin ++ ": g(pMin) = " ++ show gMin
  trace $ "pMax = " ++ show pMax ++ ": g(pMax) = " ++ show gMax
  trace $ "bids:"
  mapM_ (trace . show) bids'
  -- error $ show p ++ ", " ++ show pMin ++ ", " ++ show pMax
  --       ++ "\n"
  --       ++ show g0 ++ ", " ++ show gMin ++ ", " ++ show gMax
  --       ++ "\n"
  --       ++ show activeGoods
  --       ++ "\n"
  --       ++ show good

  apPrice .= p'
  apBids .= bids'
  where
    trace = lift . trace'


-- | The Projection Preserving Marginality as described in __Algorithm 21__ of
-- the validityAllocation paper, applied to all bids in an allocation problem.
ppm :: (Num a, Ord a, Monad m) => StateT (AllocationProblemOf b a) m ()
ppm = do
  p <- use apPrice
  apBids %= ppmBids p

-- | The Projection Preserving Marginality as described in __Algorithm 21__ of
-- the validityAllocation paper, applied to a list of bids.
ppmBids :: (Num a, Ord a) => Vector a -> [DotBidOf b a] -> [DotBidOf b a]
ppmBids p = map (ppmBid p)

-- | The Projection Preserving Marginality as described in __Algorithm 21__ of
-- the validityAllocation paper, applied to a single bid.
ppmBid :: (Num a, Ord a) => Vector a -> DotBidOf b a -> DotBidOf b a
ppmBid p =
  over dotBidPrices (ppmPrice p)

-- | The Projection Preserving Marginality as described in __Algorithm 21__ of
-- the validityAllocation paper, applied to a price vector.
ppmPrice :: (Num a, Ord a) => Vector a -> Vector a -> Vector a
ppmPrice p v =
  componentWise step optPV v
  where
    (_optval, optPV) = opt p v

    -- R can be picked from all { R ∈ ℝ | R > 1 }; we pick the smallest
    -- possible integer, which is 2.
    r = 2

    step3 opti vi =
      if opti /= 0 then vi else (-r)
    step4 opti vi =
      if opti /= 0 then vi + r else (-r)

    step =
      if (optPV !? 0) == 0 then step4 else step3

allocateJiggleL :: (Monad m, Show b, Eq b, Ord b, Show b)
                => (String -> m ()) -> Good -> BidAgentID b -> StateT (AllocationProblemOf b Frac2) m ()
allocateJiggleL trace' good agent = do
  -- Remember the original price, because we need it to un-jiggle in the last
  -- step.
  ap0 <- get
  trace "allocateJiggle"
  traceState "@1: before"
  get >>= checkAllocationProblemInvariants ap0
  p <- use apPrice
  -- let dim = dimension p
  -- 1. Apply __Algorithm 20__ (jiggle)
  jiggleL trace' good agent
  traceState "@2: after jiggle"
  get >>= checkAllocationProblemInvariants ap0
  -- 2. Apply __Algorithm 16__ (obvious refinement)
  allocateObviousL trace'
  traceState "@3: after obvious allocation"
  get >>= checkAllocationProblemInvariants ap0
  -- 3. Apply __Algorithm 21__ (projection preserving marginality)
  ppm
  traceState "@4: after ppm"
  get >>= checkAllocationProblemInvariants ap0
  -- 4. Un-jiggle.
  apPrice .= p
  apBids %= map (unjiggleBid good agent)
  traceState "@5: final"
  get >>= checkAllocationProblemInvariants ap0
  where
    trace = lift . trace'
    traceState label = do
      trace label
      use apBids >>= mapM_ (trace . show)
      use apBids >>= trace . ("#bids: " ++) . show . length
      use apResidualSupply >>= trace . show

moveBid :: (Num a) => Vector a -> DotBidOf b a -> DotBidOf b a
moveBid p =
  over dotBidPrices (<+> p)

unjiggleBid :: (Num a, Fractional a, Show b, Eq b) => Good -> BidAgentID b -> DotBidOf b a -> DotBidOf b a
unjiggleBid good agent bid =
  if bid ^. dotBidAgentID == agent then
    moveBid p bid
  else
    bid
  where
    p = vscale (-0.5) $ eX dim good
    dim = dimension (bid ^. dotBidPrices)

transformState :: Monad m => (s -> t) -> (t -> s) -> StateT s m a -> StateT t m a
transformState fwd rev action = do
  s <- rev <$> get
  (a, s') <- lift (runStateT action s)
  put $ fwd s'
  return a

data UnravelState b a =
  UnravelState
    { _usAP :: AllocationProblemOf b a
    , _usGoodsQueue :: [Good]
    }
    deriving (Show)

makeLenses ''UnravelState

pop :: MonadState s m => Lens s s [a] [a] -> m (Maybe a)
pop lens = do
  items <- use lens
  case items of
    x:xs -> do
      lens .= xs
      return (Just x)
    _ -> return Nothing

-- | __Algorithm 23__ from the validityAllocation paper. Unravel all the
-- unambiguous allocations in an allocation problem.
unravelUnambiguous :: forall a b m. (Num a, Ord a, Show a, Eq b, Ord b, Show b, Monad m)
                   => StateT (AllocationProblemOf b a) m ()
unravelUnambiguous = unravelUnambiguousL (const $ return ())

-- | Flavor of 'unravelUnambiguous' with monadic trace logging.
unravelUnambiguousL :: forall a b m. (Num a, Ord a, Show a, Eq b, Ord b, Show b, Monad m)
                    => (String -> m ())
                    -> StateT (AllocationProblemOf b a) m ()
unravelUnambiguousL trace' = do
  goodsQueue <- Set.toList . Set.delete 0 <$> use apActiveGoods
  trace $ "Goods to process: " ++ show goodsQueue
  transformState _usAP (\ap -> UnravelState ap goodsQueue) go
  where
    trace :: forall s. String -> StateT s m ()
    trace = lift . trace'

    go :: StateT (UnravelState b a) m ()
    go = do
      -- Unlike the algorithm described in the paper, we remove i from L here
      -- unconditionally, because all the paths in the original algorithm end
      -- up removing it eventually anyway.
      trace "*** Next unraveNext unravel step"
      get >>= trace . show
      pop usGoodsQueue >>= \case
        Nothing -> do
          trace "Done unraveling"
          return ()
        Just good -> do
          trace $ "Unraveling good " ++ show good
          findMarginalAgents good >>= \case
            -- (c) if there exists unique j such that V_ j i p /= ∅ then proceed
            -- 2. Identify some v ∈ V_ j i p /= ∅
            [(agent, (v:_))] -> do
              trace $ "Found 1 marginal agent on good " ++ show good ++ ": " ++ show agent
              -- 3. Apply __Algorithm 18__
              p <- use (usAP . apPrice)
              bids <- use (usAP . apBids)
              let agentBids = filter ((== agent) . (^. dotBidAgentID)) bids
              trace $ "Agent bids: " ++ show agentBids
              let keyGoods = findKeyList agentBids p good
              -- 4. Use __Algorithm 17__ to search for link goods
              let linkGoods = Set.filter (isLinkGood bids p) keyGoods
              trace $ "Key list: " ++ show keyGoods
              trace $ "Link goods: " ++ show linkGoods
              case Set.toList linkGoods of
                -- 5. If 0 link goods are found, apply __Algorithm 19__
                [] -> do
                  trace $ "No link goods found, allocating by key list"
                  zoom usAP $ allocateByKeyListL trace' agent keyGoods
                  go
                -- 7. If exactly one link good is found, let this be i'.
                [linkGood] -> do
                  trace $ "Exactly one link good found, allocating by key list"
                  -- (a) Apply __Algorithm 19__.
                  -- (b) Reset A' to be the result of this.
                  zoom usAP $ allocateByKeyListL trace' agent keyGoods
                  -- (c) Remove i from the head of L and prepend i' to the
                  -- head of L, removing it from its position later in L if it
                  -- occurs later in L.
                  usGoodsQueue %= filter (/= linkGood)
                  usGoodsQueue %= (linkGood:)
                  go

                -- 6. As soon as a second link good has been found, remove i
                -- from the head of L and go to the head of this loop.
                _ -> do
                  trace $ "Found multiple link goods, skip"
                  go

            -- (a) / (b)
            _ -> do
              trace "Cannot unravel on this good (#agents /= 1)"
              go

    findMarginalAgents :: (Num a, Ord a, Monad m)
                       => Good
                       -> StateT (UnravelState b a) m [(BidAgentID b, [DotBidOf b a])]
    findMarginalAgents good = do
      bids <- use (usAP . apBids)
      p <- use (usAP . apPrice)
      let marginalAgentBids =
            bidsByAgents .
            marginalBidsOn good p $
            bids
      let marginalAgents =
            filter (not . null . snd) .
            Map.toList $
            marginalAgentBids
      return marginalAgents

data PrioState b a =
  PrioState
    { _prioAP :: AllocationProblemOf b a
    , _prioList :: [(Good, BidAgentID b)]
    }
    deriving (Show)

makeLenses ''PrioState

-- | Run a 'MonadState' action, but track the complete state before and after
-- the action, and return a boolean indicating whether they are different
-- (as per 'Eq').
detectStateChanges :: (Eq s, MonadState s m)
                  => m a -> m Bool
detectStateChanges action = do
  s0 <- get
  _  <- action
  s1 <- get
  return $ s0 /= s1

-- | __Algorithm 24__ from the validityAllocation paper: priority-based
-- procedure to solve an allocation problem.
allocatePriorityBased :: (Ord prio, Eq b, Ord b, Show b, Monad m)
                      => (Good -> BidAgentID b -> prio) -- ^ Priority function
                      -> StateT (AllocationProblem b) m ()
allocatePriorityBased =
  allocatePriorityBasedL (const $ return ())

-- | Flavor of 'allocatePriorityBased' with monadic trace logging.
allocatePriorityBasedL :: forall prio b m
                        . (Ord prio, Eq b, Ord b, Show b, Monad m)
                       => (String -> m ()) -- ^ Trace logger
                       -> (Good -> BidAgentID b -> prio) -- ^ Priority function
                       -> StateT (AllocationProblem b) m ()
allocatePriorityBasedL trace' getPrio = do
  -- Initialise Priority to be a list ordering all elements of Q' × J in order
  -- of priority
  trace $ "** Original situation"
  ap0 <- get
  trace . show $ ap0
  activeGoods <- use apActiveGoods
  agents <- nub . map (^. dotBidAgentID) <$> use apBids
  let prios =
        sortBy
          (\a b -> comparing (uncurry getPrio) b a)
          [ (i,j)
          | i <- toList activeGoods
          , j <- agents
          ]
  transformState _prioAP (\ap -> PrioState ap prios) go
  get >>= checkAllocationProblemInvariants ap0
  where
    trace :: forall s. String -> StateT s m ()
    trace = lift . trace'

    go = do
      -- 1. Apply first __Algorithm 16__ and then __Algorithm 21__
      trace $ "** First obvious allocation"
      zoom prioAP $ do
        allocateObviousL trace'
        ppm
      trace . show =<< get
      r <- use $ prioAP . apResidualSupply
      unless (isZero r) step2

    filterPrioList = do
      activeGoods <- use (prioAP . apActiveGoods)
      prioList %= filter ((`Set.member` activeGoods) . fst)

    step2 = do
      trace $ "** Unravel"
      zoom prioAP (unravelUnambiguousL trace')
      trace . show =<< get
      activeGoods <- use (prioAP . apActiveGoods)
      r <- use $ prioAP . apResidualSupply
      unless (Set.null activeGoods || isZero r) step4

    step4 = do
      trace "** Jiggle"
      filterPrioList
      nextPairMay <- pop prioList
      case nextPairMay of
        Nothing -> do
          trace "*** Priority list empty, exiting."
          return ()
        Just (good, agent) -> do
          trace $ "*** Jiggle on " ++ show (good, agent)
          anythingChanged <-
            zoom prioAP .
            detectStateChanges .
            transformState (fmap unsafePriceFromFrac2) (fmap priceToFrac2) $
              allocateJiggleL trace' good agent
          if anythingChanged then do
            trace . show =<< get
            step2
          else do
            trace "*** Nothing changed."
            step4

-- | A raw bid, as provided by the user.
data IncomingBid =
  IncomingBid
    { _incomingBidQuantity :: Units
    , _incomingBidPrices :: PriceVector
    }
    deriving (Show, Eq, Generic)

makeLenses ''IncomingBid

instance ToJSON IncomingBid where
  toJSON (IncomingBid w v) = toJSON (w, v)

instance FromJSON IncomingBid where
  parseJSON x = do
    (w, v) <- parseJSON x
    return $ IncomingBid w v

instance CSV.ToRecord IncomingBid where
  toRecord (IncomingBid q p) =
    CSV.record $
      CSV.toField q :
      map CSV.toField (toList p)

instance CSV.FromRecord IncomingBid where
  parseRecord xs =
    IncomingBid
      <$> CSV.index xs 0
      <*> CSV.parseRecord (V.drop 1 xs)

-- | An auction, as provided by the user.
-- Before the algorithms can process this, preprocessing as defined in
-- __Algorithm 1__ is required, which is provided by 'prepareAP'
data Auction b =
  Auction
    { _auctionReservePrice :: PriceVector
    , _auctionSupply :: Bundle
    , _auctionBids :: Map b [IncomingBid]
    }
    deriving (Show, Eq)

makeLenses ''Auction

instance (Ord b, FromJSONKey b) => FromJSON (Auction b) where
  parseJSON (JSON.Object o) =
    Auction
      <$> o .: "ReservePrice"
      <*> o .: "Supply"
      <*> o .: "Bids"
  parseJSON _ = fail "Invalid JSON for Auction"

instance ToJSONKey b => ToJSON (Auction b) where
  toJSON (Auction price supply bids) =
    JSON.object
      [ ("ReservePrice", toJSON price)
      , ("Supply", toJSON supply)
      , ("Bids", toJSON bids)
      , ("Goods", toJSON [1..dimension supply])
      ]

data AuctionResult b =
  AuctionResult
    { _auctionResultPrices :: PriceVector
    , _auctionResultEndowment :: Map b Bundle
    , _auctionResultRemainingSupply :: Bundle
    }
    deriving (Show, Eq)

makeLenses ''AuctionResult

instance (ToJSON b, ToJSONKey b, Eq b, Show b) => ToJSON (AuctionResult b) where
  toJSON (AuctionResult prices endowment supply) =
    JSON.object
      [ ("Prices", toJSON prices)
      , ("Endowment", toJSON endowment)
      , ("RemainingSupply", toJSON supply)
      ]

-- | Run a full auction, using the priority-based algorithm (__Algorithm 24__
-- in validityAllocation).
runAuction :: forall b
            . (Show b, Ord b)
           => Auction b
           -> AuctionResult b
runAuction auction =
  runIdentity $ runAuctionM (const $ return ()) auction

-- | Flavor of 'runAuction' with monadic trace logging.
runAuctionM :: forall b m
             . (Show b, Ord b, Monad m)
            => (String -> m ())
            -> Auction b
            -> m (AuctionResult b)
runAuctionM trace auction = do
  (reservePrice, ap) <- prepareAP trace auction
  ap' <- execStateT (allocatePriorityBasedL trace (,)) ap
  return $ resultFromAP reservePrice ap'

-- | Convert a solved allocation problem to an auction result.
-- This covers the following:
--
-- - Reverting reserve price normalization
-- - Removing auctioneer endowments
-- - Removing REJECT coordinates from all endowments, prices, and bundles
-- - Collecting all remaining supply as well as supply allocated to the
--   auctioneer into a single rejected bundle
-- - Carrying over good labels from the original auction
resultFromAP :: Ord b
             => PriceVector -> AllocationProblem b -> AuctionResult b
resultFromAP reservePrice ap =
  AuctionResult
    (reduce $ (ap ^. apPrice) <+> reservePrice)
    (fmap reduce . removeAuctioneer $ ap ^. apEndowment)
    rejectedBundle

  where
    rejectedBundle = reduce $
      (ap ^. apResidualSupply) <+>
      (fromMaybe (zero dim) . Map.lookup Auctioneer $ ap ^. apEndowment)

    dim = dimension (ap ^. apResidualSupply)

removeAuctioneer :: (Ord b, Eq b) => Map (BidAgentID b) a -> Map b a
removeAuctioneer bids =
  Map.fromList [ (b, v) | (Bidder b, v) <- Map.toList bids ]

-- | Prepare an auction for processing, creating an allocation problem.
-- This includes:
--
-- - Adding a 0th element to every price and bundle (the REJECT good)
-- - Applying reserve price normalization
-- - Adding an auctioneer bidder and auctioneer bids
-- - Finding an equilibrium price
prepareAP :: (Ord b, Show b, Monad m)
          => (String -> m ()) -> Auction b -> m (PriceVector, AllocationProblem b)
prepareAP trace_ auction = do
  let reservePrice = (extend $ auction ^. auctionReservePrice)
  ap <- mkAllocationProblemEx
          trace_
          reservePrice
          (\trace -> findBundle_GreedyUpMinimalM trace minimalNeighbour)
          (prepareBids $ auction ^. auctionBids)
          (extend $ auction ^. auctionSupply)
  return (reservePrice, ap)

-- | Prepare a list of bids for processing.
prepareBids :: Ord b => Map b [IncomingBid] -> [DotBid b]
prepareBids =
  concat . map (uncurry prepareAgentBids) . Map.toList

-- | Prepare a list of bids from one agent for processing.
prepareAgentBids :: b -> [IncomingBid] -> [DotBid b]
prepareAgentBids agent = map (prepareAgentBid agent)

-- | Prepare one agent + bid pair for processing.
prepareAgentBid :: b -> IncomingBid -> DotBid b
prepareAgentBid agent incoming =
  DotBid
    { _dotBidQuantity = incoming ^. incomingBidQuantity
    , _dotBidPrices = extend $ incoming ^. incomingBidPrices
    , _dotBidAgentID = Bidder agent
    }

unprepareBid :: DotBid b -> IncomingBid
unprepareBid (DotBid quantity prices _) =
  IncomingBid quantity (reduce prices)

-- | Generate a valid auction with multiple bidders.
genAuction :: (Ord b, Eq b)
           => Integer -- ^ Maximum weight per bid (1 for unit-bids only)
           -> Integer -- ^ Maximum price for bids
           -> Int -- ^ Number of goods in the auction / dimensionality
           -> Int -- ^ The @b@ parameter from the paper
           -> [b] -- ^ Bidder names. List length implies number of bidders.
           -> Gen (Auction b)
genAuction maxWeight maxPrice n b bidderNames = do
  bidsList <- forM
    bidderNames $ \agent -> do
      bbids <- snd <$> generateValidBidsW maxWeight maxPrice n b (Bidder agent)
      return (agent, map unprepareBid bbids)
  let bids = Map.fromList bidsList

  supply <- fromList <$> forM [1..n] (const $ Units . fromIntegral <$> choose (0, maxWeight))
  reservePrice <- fromList <$> forM [1..n] (const $ Price <$> choose (0, maxPrice))
  return $ Auction reservePrice supply bids


-- A list of concrete tasks:
--
-- - try to reimplement bundleDemanded based on Lemma 2,
--   in such a way that it additionally indicates whether the bundle returned
--   is unique
-- - then hasNoMarginalBids becomes trivial in terms of bundleDemanded
--   ideally implement TwoPhaseMinMin
--
-- Done:
-- - reimplement findBundle crawling upwards from 0;
-- - create a traced version of findBundle
--

-- Notes from the call with Elizabeth 2017-07-11:
--
-- The function g is constant within a region (including its borders)
-- as determined by the bids when called with the bundle demanded in that
-- region. It also is minimal.
--
-- Outside of the region, g is still piecewise linear. The rate at which it
-- is increasing in a particular direction depends on the difference between
-- the target bundle and the bundle corresponding to the region we're calling
-- it in.
--
-- gTest implements a small part of this, but more is needed.
--
-- Furthermore, bundleDemanded can be replaced by Lemma 2. The calculation
-- in Lemma 2 gives us a bundle that is demanded at the current price.
--
-- The calculation of V_j says something about the borders. If any of the
-- inequalities specifying the region is weak, then we're in fact in a border.
-- The implication requiring j to be sufficiently large for an equality is
-- actually an arbitrary way to count the borders as parts of some regions,
-- breaking ties according to the order of goods.
--
--
-- Notes about ellipsoid.pdf (version of July 17)
--
-- Intro to Section 1 ok
-- Observation 1 -> End of Section 1.0 cut
-- Bottom of page 3 (Section 1.1 is fine) up to end of Section 2
--
-- Section 2 is ok, but irrelevant
--
-- Section 3, in particular the three algos, are all relevant.
-- We want to start at 0 and crawl up.
--
-- TwoPhaseMinMin seems the superior algorithm, because it always
-- finds the lower left corner.
--
-- The other relevant thing is Corollary 2 on page 14.
--
-- In particular Corollary 2.2 gives the optimised case when
-- in a UDR.
--
-- Corollary 2.1 and 2.3 are further optimisations, and require
-- knowing all demanded bundles.
--
-- 2.2 can be used if all inequalities are strict in the formula
-- to compute a demanded Bundle (Lemma 2 at the end, I think; see
-- the other notes above).

-- Notes about allocation phase (call with Elizabeth 2017-08-02):
--
-- [Clarified Lemma 2, see comments above.]
--
-- We know:
-- g is minimized at the point,
-- know the bundle,
-- don't know how to allocate the bundle to the bidders.
--
-- In principle, this works by tweaking the prices the bidders
-- bid for slightly and then re-running the algorithm.
--
-- For an individual bid that is non-marginal, there is only one
-- thing we can do (i.e., if that bid is accepted, we have to
-- allocate one unit of one particular good to that bidder).
--
-- The critical bids are the ones that are marginal, because for
-- them we have to decide which of the goods they're getting, or
-- if they're getting a good at all.
--
-- In principle, there are two approaches:
--
-- 1. This is an iterative approach. We first allocate all the
-- non-marginal bids. Then we're left with a residual problem
-- and only marginal bids. We then select some of these bids and
-- tweak their prices a bit (effectively by 0.5 in one direction
-- or the other). Then we re-run the algorithm. Elizabeth will
-- have proved that at least one bid that was marginal becomes
-- now non-marginal. So we can allocate all non-marginal goods
-- again. Then return to the original prices, and select other
-- bids to tweak. Repeat until everything is allocated.
--
-- 2. If we want to avoid returning to the original prices, we
-- could tweak by increasingly small amounts (0.5, 0.25, 0.125, ...).
-- We could then, in theory, perform multiple such tweaks at once
-- and would have to run the algorithm less often.
--
-- However, the fractions we work with could get extremely small,
-- and the algorithms we're working with, in particular submodular
-- function minimisation, would have to be able to properly deal
-- with all these small fractions.
--
-- It is currently unclear whether this would work.
--
-- So Andres should perhaps check the feasibility of approach 2.
-- On the other hand, approach 1 seems uncritical and not very
-- difficult, so we could also just go for that one.

{-
vim: sw=2 ts=2
-}