text:
# Advanced APIs

export
   askSetMode, askHungryStatus, askSleepyStatus,
   askPawPadsStatus,
   askBodyTouchedStatus,
   askPostureStatus, askChangePosture, askMoveHead, askStay,
   askBitingStatus, askPlayBone, askPlayDice, askPlayMotion,
   askNameCalledStatus, askVoiceCommandStatus,
   askApproachPerson, askChasePerson, askFindPerson,
   askExplore, askMoveToPosition,
   askTurnAround, askMoveAlongCircle, askMoveDirection,
   askMoveForward, askMoveSideways,
   askFoundObjectsStatus, askApproachObject, askChaseObject,
   askFindObject, askGetCloseToObject,
   askKickObject, askReleaseObject

## Status related APIs

"""
    askSetMode(
      modeName="NORMAL";
      enqueue=false,
      target_deviceID=nothing,
      target_nickname=nothing,
      timeoutLimit=10)

asks `set_mode`.

- `askSetMode` should be one of the following strings:
  - `NORMAL`
  - `DEVELOPMENT`

This method is equivalent to `askAction("set_mode", Dict(ModeName=>modeName, Enqueue=>enqueue))`

"""
function askSetMode(
   modeName="NORMAL";
   enqueue=false,
   target_deviceID=nothing,
   target_nickname=nothing,
   timeoutLimit=10)
   arguments = Dict("ModeName"=> modeName)
   if enqueue
      arguments["Enqueue"]=enqueue
   end
   askAction("set_mode", arguments,
      target_deviceID=target_deviceID, target_nickname=target_deviceID, timeoutLimit=timeoutLimit)
end

"""
    askHungryStatus(;
      target_deviceID=nothing,
      target_nickname=nothing,
      timeoutLimit=10)

asks `hungry_status`.

This method is equivalent to `askAction("hungry_status")`
"""
function askHungryStatus(;
   target_deviceID=nothing,
   target_nickname=nothing,
   timeoutLimit=10)
   askAction("hungry_status",
      target_deviceID=target_deviceID, target_nickname=target_deviceID, timeoutLimit=timeoutLimit)
end

"""
    askSleepyStatus(;
      target_deviceID=nothing,
      target_nickname=nothing,
      timeoutLimit=10)

asks `sleepy_status`.

This method is equivalent to `askAction("sleepy_status")`

"""
function askSleepyStatus(;
   target_deviceID=nothing,
   target_nickname=nothing,
   timeoutLimit=10)
   askAction("sleepy_status",
      target_deviceID=target_deviceID, target_nickname=target_deviceID, timeoutLimit=timeoutLimit)
end

"""
    askBodyTouchedStatus(;
      target_deviceID=nothing,
      target_nickname=nothing,
      timeoutLimit=10)

asks `body_touched_status`.

This method is equivalent to `askAction("body_touched_status")`

"""
function askBodyTouchedStatus(;
   target_deviceID=nothing,
   target_nickname=nothing,
   timeoutLimit=10)
   askAction("body_touched_status",
      target_deviceID=target_deviceID, target_nickname=target_deviceID, timeoutLimit=timeoutLimit)
end

"""
    askPawPadsStatus(;
      target_deviceID=nothing,
      target_nickname=nothing,
      timeoutLimit=10)

asks `paw_pads_status`.

This method is equivalent to `askAction("paw_pads_status")`

"""
function askPawPadsStatus(;
   target_deviceID=nothing,
   target_nickname=nothing,
   timeoutLimit=10)
   askAction("paw_pads_status",
      target_deviceID=target_deviceID, target_nickname=target_deviceID, timeoutLimit=timeoutLimit)
end

## Posture related APIs

"""
    askPostureStatus(;
      target_deviceID=nothing,
      target_nickname=nothing,
      timeoutLimit=10)

asks `posture_status`.

This method is equivalent to `askAction("posture_status")`

"""
function askPostureStatus(;
   target_deviceID=nothing,
   target_nickname=nothing,
   timeoutLimit=10)
   askAction("posture_status",
      target_deviceID=target_deviceID, target_nickname=target_deviceID, timeoutLimit=timeoutLimit)
end

"""
    askChangePosture(
      finalPosture="sit";
      enquee=false,
      target_deviceID=nothing,
      target_nickname=nothing,
      timeoutLimit=10)

asks `change_posture`.

- `finalPosture` should be one of the following constants:
  - `back`
  - `crouch`
  - `down`
  - `down_and_lengthen_behind`
  - `down_and_shorten_behind`
  - `sit_and_raise__both_hands`
  - `sit`
  - `sleep`
  - `stand`
  - `stand_straight`

This method is equivalent to `askAction("change_posture", Dict(FinalPosture=>finalPosture, Enqueue=>enqueue))`

"""
function askChangePosture(
   finalPosture="sit";
   enqueue=false,
   target_deviceID=nothing,
   target_nickname=nothing,
	timeoutLimit=10)
   arguments = Dict("FinalPosture"=> finalPosture)
   if enqueue
      arguments["Enqueue"]=enqueue
   end
   askAction("change_posture", arguments,
      target_deviceID=target_deviceID, target_nickname=target_deviceID, timeoutLimit=timeoutLimit)
end

"""
    askMoveHead(;
      azimuth=0,
      elevation=0,
      velocity=0,
      enqueue=false,
      target_deviceID=nothing,
      target_nickname=nothing,
      timeoutLimit=10)

asks `move_head`.

- `azimuth` : -80 -- 80 (degree)
- `elevation` : -40 -- 40 (degree)
- `velocity` : 10 -- 80 (deg/sec)

This method is equivalent to `askAction("move_head", Dict(Duration=>duration, Enqueue=>enqueue))`

"""
function askMoveHead(;
   azimuth=0,
   elevation=0,
   velocity=0,
   enqueue=false,
   target_deviceID=nothing,
   target_nickname=nothing,
	timeoutLimit=10)
   arguments = Dict(
      "Azimuth" => azimuth,
      "elevation" => elevation,
      "Velocity" => velocity)
   if enqueue
      arguments["Enqueue"]=enqueue
   end
   askAction("move_head", arguments,
      target_deviceID=target_deviceID, target_nickname=target_deviceID, timeoutLimit=timeoutLimit)
end

"""
    askStay(;
      duration=60,
      enqueue=false,
      target_deviceID=nothing,
      target_nickname=nothing,
      timeoutLimit=10)

asks `stay` for `duration` seconds.

- `duration` : 0--360 (seconds)

This method is equivalent to `askAction("stay", Dict(Duration=>duration, Enqueue=>enqueue))`

"""
function askStay(;
   duration=60,
   enqueue=false,
   target_deviceID=nothing,
   target_nickname=nothing,
	timeoutLimit=10)
   arguments = Dict("Duration"=> duration)
   if enqueue
      arguments["Enqueue"]=enqueue
   end
   askAction("stay", arguments,
      target_deviceID=target_deviceID, target_nickname=target_deviceID, timeoutLimit=timeoutLimit)
end

## Play related APIs

"""
    askBitingStatus(;
      target_deviceID=nothing,
      target_nickname=nothing,
      timeoutLimit=10)

asks `biting_status`.

This method is equivalent to `askAction("biting_status")`

"""
function askBitingStatus(;
   target_deviceID=nothing,
   target_nickname=nothing,
   timeoutLimit=10)
   askAction("biting_status",
      target_deviceID=target_deviceID, target_nickname=target_deviceID, timeoutLimit=timeoutLimit)
end

"""
    askPlayBone(
      category="holdMouth";
      enqueue=false,
      target_deviceID=nothing,
      target_nickname=nothing,
      timeoutLimit=10)

asks `play_bone`.

This method is equivalent to `askAction("play_bone", Dict(Category=>holdMouth, Enqueue=>enqueue))`

"""
function askPlayBone(
   category="holdMouth";
   enqueue=false,
   target_deviceID=nothing,
   target_nickname=nothing,
   timeoutLimit=10)
   arguments = Dict("Category"=> category)
   if enqueue
      arguments["Enqueue"]=enqueue
   end
   askAction("play_bone", arguments,
      target_deviceID=target_deviceID, target_nickname=target_deviceID, timeoutLimit=timeoutLimit)
end

"""
    askPlayDice(
      category="holdMouthDice";
      enqueue=false,
      target_deviceID=nothing,
      target_nickname=nothing,
      timeoutLimit=10)

asks `change_posture`.

- `category` should be one of the following constants:
  - `holdMouthDice`
  - `rollDiceLeft`
  - `rollDiceRight`
  - `rollDicePush`
  - `rollDicePull`
  - `stackDice`

This method is equivalent to `askAction("play_dice", Dict(Category=>category, Enqueue=>enqueue))`

"""
function askPlayDice(
   category="holdMouthDice";
   enqueue=false,
   target_deviceID=nothing,
   target_nickname=nothing,
   timeoutLimit=10)
   arguments = Dict("Category"=> category)
   if enqueue
      arguments["Enqueue"]=enqueue
   end
   askAction("play_dice", arguments,
      target_deviceID=target_deviceID, target_nickname=target_deviceID, timeoutLimit=timeoutLimit)
end

"""
    askPlayMotion(
      category="agree";
      enqueue=false,
      mode="NONE",
      target_deviceID=nothing,
      target_nickname=nothing,
      timeoutLimit=10)

asks `play_motion`.

- `category` should be one of the following constants. Optional `Mode` parameters other than "NONE" are also listed.
  - `agree`
  - `bad`
  - `bark`
  - `beckon`
    - `Mode`: `BODY_LEFT`
    - `Mode`: `BODY_RIGHT`
  - `belch`
  - `bentBack`
  - `breath`
  - `curious`
  - `dance`
  - `drawInOnesChin`
  - `dreaming`
  - `friendly`
  - `friendlyPolite`
  - `gasp`
  - `halfAsleep`
  - `handsUp`
  - `happyOrHot`
  - `heading`
    - `Mode`: `NONE`
    - `Mode`: `SPACE_LEFT`
    - `Mode`: `SPACE_RIGHT`
  - `helloIloveYou`
  - `highFive`
    - `Mode`: `NONE`
    - `Mode`: `BODY_LEFT`
    - `Mode`: `BODY_RIGHT`
  - `imFriendly`
  - `jiggle`
  - `kiss`
  - `lickBody`
    - `Mode`: `BODY_LEFT`
    - `Mode`: `BODY_RIGHT`

  - `lookAroundHead`
  - `lookOver`
  - `marking`
    - `Mode`: `BOY`
    - `Mode`: `GIRL`
  - `nodHead`
  - `openMouth`
  - `overJoyed`
  - `paw`
    - `Mode`: `BODY_LEFT`
    - `Mode`: `BODY_RIGHT`
  - `peace`
    - `Mode`: `SPACE_LEFT`
    - `Mode`: `SPACE_RIGHT`
  - `perceive`
  - `playBiting`
  - `prettyPlease`
  - `ready`
  - `relax`
  - `restless`
  - `rubBack`
  - `scratchFloor`
  - `scratchHead`
    - `Mode`: `HIND_LEFT`
    - `Mode`: `HIND_RIGHT`
  - `shakeHead`
  - `shakeHipsBehind`
  - `shudder`
  - `sideKick`
    - `Mode`: `FRONT_LEFT`
    - `Mode`: `FRONT_RIGHT`
  - `sideUp`
    - `Mode`: `BODY_LEFT`
    - `Mode`: `BODY_RIGHT`
  - `sing`
  - `sleepTalking`
  - `sneeze`
  - `sniff`
  - `sniffDown`
  - `sniffUp`
  - `startled`
  - `startledLittle`
  - `stretch`
  - `swing`
  - `touched`
    - `Mode`: `SPACE_CENTER`
  - `waiting`
  - `washFace`
  - `whine`
  - `wiggleEar`
    - `Mode` : `BODY_BOTH`
  - `woof`
  - `yap`
  - `yawn`

This method is equivalent to `askAction("play_motion", Dict(Category=>category, Mode=>mode, Enqueue=>enqueue))`

"""
function askPlayMotion(
   category="agree";
   enqueue=false,
   mode="NONE",
   target_deviceID=nothing,
   target_nickname=nothing,
   timeoutLimit=10)
   arguments = Dict(
      "Category"=> category,
      "Mode"=>mode)
   if enqueue
      arguments["Enqueue"]=enqueue
   end
   askAction("play_motion", arguments,
      target_deviceID=target_deviceID, target_nickname=target_deviceID, timeoutLimit=timeoutLimit)
end

## Voice related APIs

"""
    askNameCalledStatus(;
      target_deviceID=nothing,
      target_nickname=nothing,
      timeoutLimit=10)

asks `name_called_status`.

This method is equivalent to `askAction("name_called_status")`
"""
function askNameCalledStatus(;
   target_deviceID=nothing,
   target_nickname=nothing,
   timeoutLimit=10)
   askAction("name_called_status",
      target_deviceID=target_deviceID, target_nickname=target_deviceID, timeoutLimit=timeoutLimit)
end

"""
    askVoiceCommandStatus(;
      target_deviceID=nothing,
      target_nickname=nothing,
      timeoutLimit=10)

asks `voice_command_status`.

This method is equivalent to `askAction("voice_command_status")`
"""
function askVoiceCommandStatus(;
   target_deviceID=nothing,
   target_nickname=nothing,
   timeoutLimit=10)
   askAction("voice_command_status",
      target_deviceID=target_deviceID, target_nickname=target_deviceID, timeoutLimit=timeoutLimit)
end

## Person related APIs

"""
    askApproachPerson(
      distanceFromTarget=1;
      enqueue=false,
      target_deviceID=nothing,
      target_nickname=nothing,
      timeoutLimit=10)

asks `approach_person`.

- `distanceFromTarget` : 0.4 -- 2 (meter)

This method is equivalent to `askAction("approach_person", Dict(DistanceFromTarget=>distanceFromTarget, Enqueue=>enqueue))`

"""
function askApproachPerson(
   distanceFromTarget=1;
   enqueue=false,
   target_deviceID=nothing,
   target_nickname=nothing,
   timeoutLimit=10)
   arguments = Dict("DistanceFromTarget"=> distanceFromTarget)
   if enqueue
      arguments["Enqueue"]=enqueue
   end
   askAction("approach_person", arguments,
      target_deviceID=target_deviceID, target_nickname=target_deviceID, timeoutLimit=timeoutLimit)
end

"""
    askChasePerson(
      chasingDurationMsec=15*1000;
      enqueue=false,
      target_deviceID=nothing,
      target_nickname=nothing,
      timeoutLimit=10)

asks `chase_person`.

- `chasingDurationMsec` : 0 -- 360000 (milliseconds)

This method is equivalent to `askAction("chase_person", Dict(ChasingDurationMsec=>chasingDurationMsec, Enqueue=>enqueue))`

"""
function askChasePerson(
   chasingDurationMsec=15*1000;
   enqueue=false,
   target_deviceID=nothing,
   target_nickname=nothing,
   timeoutLimit=10)
   arguments = Dict("ChasingDurationMsec"=> chasingDurationMsec)
   if enqueue
      arguments["Enqueue"]=enqueue
   end
   askAction("chase_person", arguments,
      target_deviceID=target_deviceID, target_nickname=target_deviceID, timeoutLimit=timeoutLimit)
end

"""
    askFindPerson(;
      target_deviceID=nothing,
      target_nickname=nothing,
      timeoutLimit=10)

asks `find_person`.

This method is equivalent to `askAction("find_person")`
"""
function askFindPerson(;
   target_deviceID=nothing,
   target_nickname=nothing,
   timeoutLimit=10)
   askAction("find_person",
      target_deviceID=target_deviceID, target_nickname=target_deviceID, timeoutLimit=timeoutLimit)
end

### Position related APIs

"""
    askExplore(
      duration=60;
      enqueue=false,
      target_deviceID=nothing,
      target_nickname=nothing,
      timeoutLimit=10)

asks `explore`.

- `duration` : 0 -- 360 (seconds)

This method is equivalent to `askAction("explore", Dict(Duration=>duration, Enqueue=>enqueue))`

"""
function askExplore(
   duration=60;
   enqueue=false,
   target_deviceID=nothing,
   target_nickname=nothing,
	timeoutLimit=10)
   arguments = Dict("Duration"=> duration)
   if enqueue
      arguments["Enqueue"]=enqueue
   end
   askAction("explore", arguments,
      target_deviceID=target_deviceID, target_nickname=target_deviceID, timeoutLimit=timeoutLimit)
end

"""
    askMoveToPosition(
      targetType="charging_station";
      enqueue=false,
      target_deviceID=nothing,
      target_nickname=nothing,
      timeoutLimit=10)

asks `move_to_position`.

- `targetType` : one of the following constants:
  - `charging_station`
  - `greeting_spot`
  - `toilet`

This method is equivalent to `askAction("move_to_position", Dict(TargetType=>targetType, Enqueue=>enqueue))`
"""
function askMoveToPosition(
   targetType="charging_station";
   enqueue=false,
   target_deviceID=nothing,
   target_nickname=nothing,
	timeoutLimit=10)
   arguments = Dict("TargetType"=> targetType)
   if enqueue
      arguments["Enqueue"]=enqueue
   end
   askAction("move_to_position", arguments,
      target_deviceID=target_deviceID, target_nickname=target_deviceID, timeoutLimit=timeoutLimit)
end

### Moving related APIs

"""
    askTurnAround(;
      turnAngle=0,
      turnSpeed=1,
      enqueue=false,
      target_deviceID=nothing,
      target_nickname=nothing,
      timeoutLimit=10)

asks `turn_around`.

- `turnSpeed` : one of 0 (slow), 1, 2 (fast)
- `turnAngle` : -180 -- 180, counterclockwise with respect to the current aibo angle.

This method is equivalent to `askAction("turn_around", Dict(TurnSpeed=> turnSpeed, TurnAngle=> turnAngle, Enqueue=>enqueue))`
"""
function askTurnAround(;
   turnAngle=0,
   turnSpeed=1,
   enqueue=false,
   target_deviceID=nothing,
   target_nickname=nothing,
	timeoutLimit=10)
   arguments = Dict(
      "TurnSpeed"=> turnSpeed,
      "TurnAngle"=> turnAngle)
   if enqueue
      arguments["Enqueue"]=enqueue
   end
   askAction("turn_around", arguments,
      target_deviceID=target_deviceID, target_nickname=target_deviceID, timeoutLimit=timeoutLimit)
end

"""
    askMoveAlongCircle(;
      movingAngle=360,
      radius=1,
      walkSpeed=1,
      enqueue=false,
      target_deviceID=nothing,
      target_nickname=nothing,
      timeoutLimit=10)

asks `move_along_circle`.

- `walkSpeed`   : one of 0 (slow), 1, 2 (fast)
- `radius`      : 0.5 -- 3 (meter)
- `movingAngle` : 0 -- 1080 (degree)

This method is equivalent to `askAction("move_along_circle", Dict(WalkSpeed=> walkSpeed, Radius=> radius, MovingAngle=>movingAngle, Enqueue=>enqueue))`

"""
function askMoveAlongCircle(;
   movingAngle=360,
   radius=1,
   walkSpeed=1,
   enqueue=false,
   target_deviceID=nothing,
   target_nickname=nothing,
	timeoutLimit=10)
   arguments = Dict(
      "WalkSpeed"=> walkSpeed,
      "Radius"=> radius,
      "MovingAngle"=> movingAngle)
   if enqueue
      arguments["Enqueue"]=enqueue
   end
   askAction("move_along_circle", arguments,
      target_deviceID=target_deviceID, target_nickname=target_deviceID, timeoutLimit=timeoutLimit)
end

"""
    askMoveDirection(;
      targetAngle=0,
      targetDistance=1,
      walkSpeed=1,
      enqueue=false,
      target_deviceID=nothing,
      target_nickname=nothing,
      timeoutLimit=10)

asks `move_direction`.

- `walkSpeed`      : one of 0 (slow), 1, 2 (fast)
- `targetDistance` : 0 -- 6 (meter)
- `targetAngle`    : -180 -- 180, counterclockwise with respect to the current aibo angle.

This method is equivalent to `askAction("move_direction", Dict(WalkSpeed=> walkSpeed, TargetDistance=> targetDistance, TargetAngle=> targetAngle, Enqueue=>enqueue))`
"""
function askMoveDirection(;
   targetAngle=0,
   targetDistance=1,
   walkSpeed=1,
   enqueue=false,
   target_deviceID=nothing,
   target_nickname=nothing,
	timeoutLimit=10)
   arguments = Dict(
      "WalkSpeed"=> walkSpeed,
      "TargetDistance"=> targetDistance,
      "TargetAngle"=> targetAngle)
   if enqueue
      arguments["Enqueue"]=enqueue
   end
   askAction("move_direction", arguments,
      target_deviceID=target_deviceID, target_nickname=target_deviceID, timeoutLimit=timeoutLimit)
end

"""
    askMoveForward(;
      walkDistance=1,
      walkSpeed=1,
      enqueue=false,
      target_deviceID=nothing,
      target_nickname=nothing,
      timeoutLimit=10)

asks `move_forward`.

- `walkSpeed`    : one of 0 (slow), 1, 2 (fast)
- `walkDistance` : -6 -- 6 (meter)

This method is equivalent to `askAction("move_forward", Dict(WalkSpeed=> walkSpeed, WalkDistance=> walkDistance, Enqueue=>enqueue))`
"""
function askMoveForward(;
   walkDistance=1,
   walkSpeed=1,
   enqueue=false,
   target_deviceID=nothing,
   target_nickname=nothing,
	timeoutLimit=10)
   arguments = Dict(
      "WalkSpeed"=> walkSpeed,
      "WalkDistance"=> walkDistance)
   if enqueue
      arguments["Enqueue"]=enqueue
   end
   askAction("move_forward", arguments,
      target_deviceID=target_deviceID, target_nickname=target_deviceID, timeoutLimit=timeoutLimit)
end

"""
    askMoveSideways(;
      walkDistance=0.5,
      walkSpeed=1,
      enqueue=false,
      target_deviceID=nothing,
      target_nickname=nothing,
      timeoutLimit=10)

asks `move_sideways`.

- `walkSpeed`    : one of 0 (slow), 1, 2 (fast)
- `walkDistance` : -6 -- 6 (meter)

This method is equivalent to `askAction("move_sideways", Dict(WalkSpeed=> walkSpeed, WalkDistance=> walkDistance, Enqueue=>enqueue))`
"""
function askMoveSideways(;
   walkDistance=0.5,
   walkSpeed=1,
   enqueue=false,
   target_deviceID=nothing,
   target_nickname=nothing,
	timeoutLimit=10)
   arguments = Dict(
      "WalkSpeed"=> walkSpeed,
      "WalkDistance"=> walkDistance)
   if enqueue
      arguments["Enqueue"]=enqueue
   end
   askAction("move_sideways", arguments,
      target_deviceID=target_deviceID, target_nickname=target_deviceID, timeoutLimit=timeoutLimit)
end

### Object related APIs

"""
    askFoundObjectsStatus(;
      target_deviceID=nothing,
      target_nickname=nothing,
      timeoutLimit=10)

asks `found_objects_status`.

This method is equivalent to `askAction("found_objects_status")`
"""
function askFoundObjectsStatus(;
   target_deviceID=nothing,
   target_nickname=nothing,
   timeoutLimit=10)
   askAction("found_objects_status",
      target_deviceID=target_deviceID, target_nickname=target_deviceID, timeoutLimit=timeoutLimit)
end

"""
    askApproachObject(
      targetType="pinkball";
      enqueue=false,
      target_deviceID=nothing,
      target_nickname=nothing,
      timeoutLimit=10)

asks `approach_object`.

- `targetType` : one of the following constants:
  - `aibo`
  - `aibone`
  - `dice`
  - `pinkball`

This method is equivalent to `askAction("approach_object", Dict(TargetType=>targetType, Enqueue=>enqueue))`
"""
function askApproachObject(
   targetType="pinkball";
   enqueue=false,
   target_deviceID=nothing,
   target_nickname=nothing,
	timeoutLimit=10)
   arguments = Dict("TargetType"=> targetType)
   if enqueue
      arguments["Enqueue"]=enqueue
   end
   askAction("approach_object", arguments,
      target_deviceID=target_deviceID, target_nickname=target_deviceID, timeoutLimit=timeoutLimit)
end

"""
    askChaseObject(
      targetType="pinkball";
      chasingDurationMsec=30*1000,
      enqueue=false,
      target_deviceID=nothing,
      target_nickname=nothing,
      timeoutLimit=10)

asks `chase_object`.

- `targetType` : one of the following constants:
  - `aibo`
  - `aibone`
  - `dice`
  - `pinkball`
- `chasingDurationMsec` : 0 -- 360000 (milliseconds)

This method is equivalent to `askAction("chase_object", Dict(TargetType=>targetType, ChasingDurationMsec=> chasingDurationMsec, Enqueue=>enqueue))`
"""
function askChaseObject(
   targetType="pinkball";
   chasingDurationMsec=30*1000,
   enqueue=false,
   target_deviceID=nothing,
   target_nickname=nothing,
	timeoutLimit=10)
   arguments = Dict(
      "TargetType"=> targetType,
      "ChasingDurationMsec"=> chasingDurationMsec)
   if enqueue
      arguments["Enqueue"]=enqueue
   end
   askAction("chase_object", arguments,
      target_deviceID=target_deviceID, target_nickname=target_deviceID, timeoutLimit=timeoutLimit)
end

"""
    askFindObject(
      targetType="pinkball";
      enqueue=false,
      target_deviceID=nothing,
      target_nickname=nothing,
      timeoutLimit=10)

asks `find_object`.

- `targetType` : one of the following constants:
  - `aibo`
  - `aibone`
  - `dice`
  - `pinkball`

This method is equivalent to `askAction("find_object", Dict(TargetType=>targetType, Enqueue=>enqueue))`
"""
function askFindObject(
   targetType="pinkball";
   enqueue=false,
   target_deviceID=nothing,
   target_nickname=nothing,
	timeoutLimit=10)
   arguments = Dict("TargetType"=> targetType)
   if enqueue
      arguments["Enqueue"]=enqueue
   end
   askAction("find_object", arguments,
      target_deviceID=target_deviceID, target_nickname=target_deviceID, timeoutLimit=timeoutLimit)
end

"""
    askGetCloseToObject(
      targetType="pinkball";
      distance=0.2,
      enqueue=false,
      target_deviceID=nothing,
      target_nickname=nothing,
      timeoutLimit=10)

asks `get_close_to_object`.
Invoke this method after `askApproachObject`.

- `targetType` : one of the following constants:
  - `aibo`
  - `aibone`
  - `dice`
  - `pinkball`

- `distance` : 0.1 -- 0.3 (meter)

This method is equivalent to `askAction("get_close_to_object", Dict(TargetType=>targetType, Distance=>distance, Enqueue=>enqueue))`
"""
function askGetCloseToObject(
   targetType="pinkball";
   distance=0.2,
   enqueue=false,
   target_deviceID=nothing,
   target_nickname=nothing,
	timeoutLimit=10)
   arguments = Dict(
      "TargetType"=> targetType,
      "Distance"=> distance)
   if enqueue
      arguments["Enqueue"]=enqueue
   end
   askAction("get_close_to_object", arguments,
      target_deviceID=target_deviceID, target_nickname=target_deviceID, timeoutLimit=timeoutLimit)
end

"""
    askKickObject(
      targetType="pinkball";
      kickMotion="kick",
      enqueue=false,
      target_deviceID=nothing,
      target_nickname=nothing,
      timeoutLimit=10)

asks `kick_object`.

- `targetType` : one of the following constants:
  - `pinkball`

- `kickMotion` : one of the following constants:
  - `kick`
  - `heading`

This method is equivalent to `askAction("kick_object", Dict(TargetType=>targetType, KickMotion=>kickMotion, Enqueue=>enqueue))`
"""
function askKickObject(
   targetType="pinkball";
   kickMotion="kick",
   enqueue=false,
   target_deviceID=nothing,
   target_nickname=nothing,
	timeoutLimit=10)
   arguments = Dict(
      "TargetType"=> targetType,
      "KickMotion"=> kickMotion)
   if enqueue
      arguments["Enqueue"]=enqueue
   end
   askAction("kick_object", arguments,
      target_deviceID=target_deviceID, target_nickname=target_deviceID, timeoutLimit=timeoutLimit)
end

"""
    askReleaseObject(
      targetType="aibone";
      enqueue=false,
      target_deviceID=nothing,
      target_nickname=nothing,
      timeoutLimit=10)

asks `release_object`.

- `targetType` : one of the following constants:
  - `aibone`
  - `dice`

This method is equivalent to `askAction("release_object", Dict(TargetType=>targetType, Enqueue=>enqueue))`
"""
function askReleaseObject(
   targetType="aibone";
   enqueue=false,
   target_deviceID=nothing,
   target_nickname=nothing,
	timeoutLimit=10)
   arguments = Dict("TargetType"=> targetType)
   if enqueue
      arguments["Enqueue"]=enqueue
   end
   askAction("release_object", arguments,
      target_deviceID=target_deviceID, target_nickname=target_deviceID, timeoutLimit=timeoutLimit)
end


text:
module Theta

using LinearAlgebra
using NLopt
using GSL: sf_gamma_inc

include("riemann_matrix.jl")
include("eval_theta.jl")

include("lattice.jl")
include("ellipsoid.jl")
include("radius.jl")

include("characteristics.jl")
include("siegel_transform.jl")

include("schottky4.jl")
include("accola.jl")
include("fgsm.jl")
include("schottky_null.jl")

export theta
export RiemannMatrix, random_siegel
export siegel_transform, symplectic_transform
export theta_char, even_theta_char, odd_theta_char, check_azygetic
export schottky_genus_4, random_nonschottky_genus_4
export accola_chars, accola, random_nonaccola
export fgsm_chars, fgsm, random_nonfgsm
export hessian, schottky_null, min_even_theta_char

end

ext: jl


text:
export tabs,
       pages,
       menu,
       submenu,
       dropdown,
       dropdownmenu,
       menuitem,
       icon,
       iconbutton,
       toolbar
# Higher-order layouts: e.g. tabs, pages

# Icons and icon button
@api icon => (Icon <: Tile) begin
    doc("An Icon.")
    typedarg(
        icon::AbstractString="",
        doc=md"""The name of the icon. Valid icons can be found in the
        [Polymer iron-icon documentation](https://www.polymer-project.org/0.5/components/iron-icons/demo.html)"""
    )
    kwarg(
        url::AbstractString="",
        doc="Optionally you can specify a url to the png/svg of the icon. icon field will be ignored if url is non-empty."
    )
end
render(i::Icon, state) = begin
    if i.url == ""
        Elem("iron-icon", attributes=@d(:icon => i.icon))
    else
        Elem("iron-icon", attributes=@d(:src => i.url))
    end
end

@api iconbutton => (IconButton <: Widget) begin
    doc("A button with an inset icon.")
    typedarg(
        icon::AbstractString="",
        doc=md"""The name of the icon. Valid icons can be found in the
        [Polymer iron-icon documentation](https://www.polymer-project.org/0.5/components/iron-icons/demo.html)"""
    )
    kwarg(
        url::AbstractString="",
        doc="Optionally you can specify a url to the png/svg of the icon. icon field will be ignored if url is non-empty."
    )
end
render(i::IconButton, state) = begin
    if i.url == ""
        Elem("paper-icon-button", attributes=@d(:icon => i.icon))
    else
        Elem("paper-icon-button", attributes=@d(:src => i.url))
    end
end

wrapbehavior(w::IconButton) =
    clickable(w)

abstract Selection <: Widget

@api pages => (Pages <: Selection) begin #FIXME: Why is this a widget?
    doc("A set of pages. Only one selected page will be visible at any given time.")
    arg(pages::TileList, doc="Pages.")
    kwarg(selected::Integer=1, doc="Index of the currently visible page.")
end

render(ps::Pages, state) =
    Elem("iron-pages",
        map(t -> Elem("div", render(t, state)), ps.pages.tiles),
        attributes = @d(:selected=>ps.selected-1))

@api tabs => (Tabs <: Selection) begin
    doc("A horizontal tab bar.")
    arg(tabs::TileList, doc="The tabs.")
    kwarg(selected::Integer=1, doc="Index of the currently selected tab.")
end

render(tabs::Tabs, state) =
    Elem("paper-tabs",
        map(t -> Elem("paper-tab", render(t, state)), tabs.tabs.tiles),
        attributes = @d(:selected=>tabs.selected-1))

wrapbehavior(t::Tabs) = selectable(t)

# Menus

@api menu => (Menu <: Selection) begin
    doc("A menu.")
    arg(
        items::TileList,
        doc="Menu items. Some of the items can also be sub-menus"
    )
    kwarg(
        multi::Bool=false,
        doc="Can multiple items be selected? If set, output signal contains a vector of indices"
    )
    kwarg(selected::Integer=1, doc="Index of the currently selected menu item.")
end

render(m::Menu, state) =
    Elem("paper-menu",
        map(x -> render(wrapitem(x), state), m.items.tiles),
        attributes = @d(:selected=>m.selected-1, :multi=>boolattr(m.multi)))

wrapbehavior(m::Menu) = selectable(m, multi=m.multi)

@api submenu => (SubMenu <: Selection) begin
    doc("A submenu. Must be put inside a menu.")
    arg(label::Tile, doc="The title of the sub menu.")
    curry(items::TileList, doc="Sub menu items.")
    kwarg(
        multi::Bool=false,
        doc="Can multiple items be selected? If set, output signal contains a vector of indices"
    )
    kwarg(
        selected::Integer=1,
        doc="Index of the currently selected sub menu item."
    )
end

render(m2::SubMenu, state) =
    Elem("paper-submenu",
        [addclasses(render(wrapitem(m2.label), state), "menu-trigger"),
         addclasses(render(menu(m2.items, multi=m2.multi,
                                selected=m2.selected), state), "menu-content")])

wrapbehavior(m::SubMenu) = selectable(m, multi=m.multi)

@api menuitem => (MenuItem <: Tile) begin
    arg(tile::Tile)
    kwarg(disabled::Bool=false)
end

render(i::MenuItem, state) =
    Elem("paper-item", render(i.tile, state),
        attributes=@d(:disabled=>boolattr(i.disabled)))

wrapitem(x::MenuItem) = x
wrapitem(x::SubMenu) = x
wrapitem(x) = menuitem(x)

# Toolbar

@api toolbar => (Toolbar <: Tile) begin
    doc("A toolbar.")
    arg(
        tiles::TileList,
        doc="""Contents of the toolbar. Will be laid out in an hbox.
               Use of icon button is recommended inside a toolbar."""
    )
end

render(t::Toolbar, state) =
    render(t.tiles, "paper-toolbar", state)

# @api dropdown => (Dropdown <: Widget) begin
#     doc(md"A dropdown. For a dropdown menu use `dropdownmenu`")
#     arg(tile::Tile, doc="Contents of the dropdown.")
#     kwarg(
#         halign::Side{Horizontal}=right,
#         doc=md"""Horizontal alingment with respect to container.
#              Valid values are `left` and `right`."""
#     )
#     kwarg(
#         valign::Side{Vertical}=top,
#         doc=md"""Vertical alingment with respect to container.
#              Valid values are `top` and `bottom`."""
#     )
# end
# render(d::Dropdown, state) =
#     Elem("iron-dropdown",
#         render(d.tile, state),
#         attributes=@d(
#             "horizontal-align"=>lowercase(name(d.halign)),
#             "vertical-align"=>lowercase(name(d.valign))))

@api dropdownmenu => (DropdownMenu <: Selection) begin
    doc("A dropdown menu.")
    arg(label::AbstractString="", doc="The label.")
    arg(menu::Tile, doc="The menu (any Selection widget).")
    kwarg(
        halign::Side{Horizontal}=left,
        doc=md"""Horizontal alingment with respect to container.
             Valid values are `left` and `right`."""
    )
    kwarg(
        valign::Side{Vertical}=top,
        doc=md"""Vertical alingment with respect to container.
             Valid values are `top` and `bottom`."""
    )
    kwarg(
        floatinglabel::Bool=true,
        doc="Whether to show a floating label at the top."
    )
    kwarg(
        opened::Bool=false,
        doc="If set to true dropdown menu will be opened."
    )
    kwarg(
        disabled::Bool=false,
        doc="If set to true, the dropdown menu is disabled."
    )
end

render(dm::DropdownMenu, state) = begin
    # paper-dropdown-menu spec requires these classes

    Elem("paper-dropdown-menu",
        addclasses(render(dm.menu, state), "dropdown-content"),
        attributes = @d(
            "label"=>dm.label,
            "opened"=>boolattr(dm.opened),
            "no-label-float"=>boolattr(!dm.floatinglabel),
        )
    )
end
wrapbehavior(d::DropdownMenu) = selectable(d, selector=".dropdown-content")


text:
include("basic_dqn.jl")
include("dqn.jl")
include("prioritized_dqn.jl")
include("rainbow.jl")


text:
module MaxHelpingHandReskinnableStarTrackSpinner

using ..Ahorn, Maple

@pardef ReskinnableStarTrackSpinner(x1::Integer, y1::Integer, x2::Integer=x1 + 16, y2::Integer=y1, speed::String="Normal", startCenter::Bool=false, spriteFolder::String="danger/MaxHelpingHand/starSpinner", particleColors::String="EA64B7|3EE852,67DFEA|E85351,EA582C|33BDE8") =
    Entity("MaxHelpingHand/ReskinnableStarTrackSpinner", x=x1, y=y1, nodes=Tuple{Int, Int}[(x2, y2)], speed=speed, startCenter=startCenter, spriteFolder=spriteFolder, particleColors=particleColors)

const placements = Ahorn.PlacementDict(
    "Star (Track, $(uppercasefirst(speed)), Reskinnable) (max480's Helping Hand)" => Ahorn.EntityPlacement(
        ReskinnableStarTrackSpinner,
        "line",
        Dict{String, Any}(
            "speed" => speed
        )
    ) for speed in Maple.track_spinner_speeds
)

Ahorn.editingOptions(entity::ReskinnableStarTrackSpinner) = Dict{String, Any}(
    "speed" => Maple.track_spinner_speeds
)

Ahorn.nodeLimits(entity::ReskinnableStarTrackSpinner) = 1, 1

function Ahorn.selection(entity::ReskinnableStarTrackSpinner)
    startX, startY = Ahorn.position(entity)
    stopX, stopY = Int.(entity.data["nodes"][1])

    return [Ahorn.Rectangle(startX - 8, startY - 8, 16, 16), Ahorn.Rectangle(stopX - 8, stopY - 8, 16, 16)]
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::ReskinnableStarTrackSpinner, room::Maple.Room)
    startX, startY = Ahorn.position(entity)
    stopX, stopY = entity.data["nodes"][1]
    spriteFolder = get(entity.data, "spriteFolder", "danger/MaxHelpingHand/starSpinner")

    Ahorn.drawSprite(ctx, "$(spriteFolder)/idle0_00", stopX, stopY)
    Ahorn.drawArrow(ctx, startX, startY, stopX, stopY, Ahorn.colors.selection_selected_fc, headLength=10)
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::ReskinnableStarTrackSpinner, room::Maple.Room)
    startX, startY = Ahorn.position(entity)
    spriteFolder = get(entity.data, "spriteFolder", "danger/MaxHelpingHand/starSpinner")

    Ahorn.drawSprite(ctx, "$(spriteFolder)/idle0_00", startX, startY)
end

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
### Py
#
# We use pickle to serialise Python objects to bytes.

function serialize_py(s, x::Py)
    if pyisnull(x)
        serialize(s, nothing)
    else
        b = pyimport("pickle").dumps(x)
        serialize(s, pybytes_asvector(b))
    end
end

function deserialize_py(s)
    v = deserialize(s)
    if v === nothing
        pynew()
    else
        pyimport("pickle").loads(pybytes(v))
    end
end

function Serialization.serialize(s::AbstractSerializer, x::Py)
    Serialization.serialize_type(s, Py, false)
    serialize_py(s, x)
end

Serialization.deserialize(s::AbstractSerializer, ::Type{Py}) = deserialize_py(s)

### PyException
#
# Traceback objects are not serialisable, but Exception objects are (because pickle just
# ignores the __traceback__ attribute), so we serialise a PyException by just serialising
# its (normalised) `v` field, from which we can recover the type and traceback.
#
# This means the user can install something like "tblib" to enable pickling of tracebacks
# and for free this enables serializing PyException including the traceback.

function Serialization.serialize(s::AbstractSerializer, x::PyException)
    Serialization.serialize_type(s, PyException, false)
    serialize_py(s, x.v)
end

Serialization.deserialize(s::AbstractSerializer, ::Type{PyException}) = PyException(deserialize_py(s))


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

struct ResponseCollector
    collectors::Dict{Bytes,Channel{Transaction_Res_All}}
    transact_result_channel::Channel{Proto.Transaction_Server}
    access_lock::ReentrantLock
end

function ResponseCollector(transact_result_channel::Channel{Proto.Transaction_Server})
    collectors = Dict{Bytes,Channel{Transaction_Res_All}}()
    access_lock = ReentrantLock()
    resp_col = ResponseCollector(collectors, transact_result_channel, access_lock)
    @async response_worker(resp_col)
    return resp_col
end

"""
    push!(resp_collector::ResponsCollector, request::T) where {T<:Proto.ProtoType}
Function is meant to give back the result_channel in which the results for one request
will be collected.
Attention! Don't put a new Id manually on the ResponsCollector. It wouldn't be thread safe
"""
function Base.push!(resp_collector::ResponseCollector, req_id::Bytes)
    res_channel = Channel{Transaction_Res_All}(10)
    try
        lock(resp_collector.access_lock)
        resp_collector.collectors[req_id] = res_channel
    finally
        unlock(resp_collector.access_lock)
    end
    return res_channel
end

"""
    delete!(resp_collector::ResponsCollector, id::Bytes)
The function will close the collecting result channel and remove
this from the the resonse collector.
Attention! Don't remove a result_channel manually from the Dictionary.
This will not be thread safe.
"""
function Base.delete!(resp_collector::ResponseCollector, id::Bytes)
    try
        lock(resp_collector.access_lock)
        close(resp_collector.collectors[id])
        delete!(resp_collector.collectors, id)
    finally
        unlock(resp_collector.access_lock)
    end
end

function response_worker(response_collector::ResponseCollector)
    resp_chan = response_collector.transact_result_channel
    collectors = response_collector.collectors
    while isopen(resp_chan)
        yield()
        try
            if isready(resp_chan)
                req_result = take!(resp_chan)
                tmp_result = _process_Transaction_Server(req_result)
                put!(collectors[tmp_result.req_id], tmp_result)
            end
        catch ex
            @info "response_worker shows an error \n
            $ex"
        end
    end
    @debug "response_collector is Done"
end

function _process_Transaction_Server(input::Proto.Transaction_Server)
    kind_of_respond = Proto.which_oneof(input,:server)
    result = getproperty(input, kind_of_respond)
    return result
end

function close(res_collector::ResponseCollector)
    #close the resul channels which are open
    safe_close.(values(res_collector.collectors))
end


text:
function euclidian(x::Vector{T},y::Vector{T}) where T<:Real
  nsz = length(x)
  dist = zeros(T,nsz,nsz)

  for i in 1:nsz, j in 1:nsz
    if i<=j
      xd = (x[i]-x[j])^2
      yd = (y[i]-y[j])^2
      dist[i,j] = dist[j,i] = round(sqrt(xd+yd),RoundNearestTiesUp)
    end
  end
  return dist
end

function manhattan(x::Vector{T},y::Vector{T}) where T<:Real
  nsz = length(x)
  dist = zeros(T, nsz, nsz)

  for i in 1:nsz, j in 1:nsz
    if i<=j
      xd = abs(x[i] - x[j])
      yd = abs(y[i] - y[j])
      dist[i, j] = dist[j, i] = round(xd + yd, RoundNearestTiesUp)
    end
  end
  return dist
end

function max_norm(x::Vector{T},y::Vector{T}) where T<:Real
  nsz = length(x)
  dist = zeros(T, nsz, nsz)

  for i in 1:nsz, j in 1:nsz
    if i<=j
      xd = abs(x[i] - x[j])
      yd = abs(y[i] - y[j])
      dist[i, j] = dist[j, i] = max(round(xd, RoundNearestTiesUp), round(yd, RoundNearestTiesUp))
    end
  end
  return dist
end

function att_euclidian(x::Vector{T},y::Vector{T}) where T<:Real
  nsz = length(x)
  dist = zeros(T,nsz,nsz)

  for i in 1:nsz, j in 1:nsz
    if i<=j
      xd = (x[i]-x[j])^2
      yd = (y[i]-y[j])^2
      dist[i,j] = dist[j,i] = ceil(sqrt((xd+yd)/10.0))
    end
  end
  return dist
end

function ceil_euclidian(x::Vector{T},y::Vector{T}) where T<:Real
  nsz = length(x)
  dist = zeros(T,nsz,nsz)

  for i in 1:nsz, j in 1:nsz
    if i<=j
      xd = (x[i]-x[j])^2
      yd = (y[i]-y[j])^2
      dist[i,j] = dist[j,i] = ceil(sqrt(xd+yd))
    end
  end
  return dist
end

function geo(x::Vector{T},y::Vector{T}) where T<:Real
  PI = 3.141592
  RRR = 6378.388
  nsz = length(x)
  dist = zeros(T,nsz,nsz)
  degs = trunc.(hcat(x,y))
  mins = hcat(x,y).-degs
  coords = PI.*(degs.+(5.0.*(mins./3.0)))./180.0
  lat = coords[:,1]
  lon = coords[:,2]

  for i in 1:nsz, j in 1:nsz
    if i<=j
      q1 = cos(lon[i]-lon[j])
      q2 = cos(lat[i]-lat[j])
      q3 = cos(lat[i]+lat[j])
      dij = RRR.*acos(0.5.*((1.0.+q1).*q2.-(1.0.-q1).*q3)).+1.0
      dist[i,j] = dist[j,i] = floor(dij)
    end
  end
  return dist
end


text:
function rollout(mdp::Union{POMDP,MDP}, s, d::Int64,
                 cost::Float64, weight::Float64, c_reduction::Symbol,
                 action_distrib_fn::Function)
    if d == 0 || isterminal(mdp, s)
        return cost, weight
    else
        p_action = action_distrib_fn(mdp, s)
        action = rand(p_action)

        (sp, r) = @gen(:sp, :r)(mdp, s, action, Random.GLOBAL_RNG)
        new_cost = update_cost(cost, r, c_reduction)
        return rollout(mdp, sp, d - 1, cost, weight, c_reduction,
                       action_distrib_fn)
    end
end

function estimate_value(mdp::Union{POMDP,MDP}, state, depth::Int,
                        cost::Float64, weight::Float64, c_reduction::Symbol,
                        action_distrib_fn::Function)
    return rollout(mdp, state, depth, cost, weight, c_reduction,
                   action_distrib_fn)
end

function next_action(gen::UniformActionGenerator, mdp::Union{POMDP,MDP}, s, snode::PISStateNode)
    rand(gen.rng, actions(mdp, s))
end

init_Q(f::Function, mdp::Union{MDP,POMDP}, s, a) = f(mdp, s, a)
init_Q(n::Number, mdp::Union{MDP,POMDP}, s, a) = convert(Float64, n)

init_N(f::Function, mdp::Union{MDP,POMDP}, s, a) = f(mdp, s, a)
init_N(n::Number, mdp::Union{MDP,POMDP}, s, a) = convert(Int, n)


text:
#############################################################################
# quantum_entropy returns -tr(X*log(X)) where X is a positive semidefinite.
# Note this function uses logarithm base e, not base 2, so return value is in
# units of nats, not bits.
#
# Quantum entropy is concave. This function implements the semidefinite
# programming approximation given in the reference below.  Parameters m and k
# control the accuracy of this approximation: m is the number of quadrature
# nodes to use and k the number of square-roots to take. See reference for
# more details.
#
# Implementation uses the expression
#   H(X) = -trace( D_{op}(X||I) )
# where D_{op} is the operator relative entropy:
#   D_{op}(X||Y) = X^{1/2}*logm(X^{1/2} Y^{-1} X^{1/2})*X^{1/2}
#
# All expressions and atoms are subtypes of AbstractExpr.
# Please read expressions.jl first.
#
#REFERENCE
#   Ported from CVXQUAD which is based on the paper: "Lieb's concavity
#   theorem, matrix geometric means and semidefinite optimization" by Hamza
#   Fawzi and James Saunderson (arXiv:1512.03401)
#############################################################################

struct QuantumEntropy <: AbstractExpr
    head::Symbol
    id_hash::UInt64
    children::Tuple{AbstractExpr,}
    size::Tuple{Int, Int}
    m::Integer
    k::Integer

    function QuantumEntropy(X::AbstractExpr, m::Integer, k::Integer)
        children = (X,)
        n = size(X)[1]
        if size(X) != (n, n)
            throw(DimensionMismatch("X must be square"))
        end
        return new(:quantum_entropy, hash(children), children, (1, 1), m, k)
    end
end

function sign(atom::QuantumEntropy)
    return Positive()
end

function monotonicity(atom::QuantumEntropy)
    return (NoMonotonicity(),)
end

function curvature(atom::QuantumEntropy)
    return ConcaveVexity()
end

function evaluate(atom::QuantumEntropy)
    X = evaluate(atom.children[1])
    return quantum_entropy(X)
end

const MatrixOrConstant = Union{AbstractMatrix, Constant}

function quantum_entropy(X::AbstractExpr, m::Integer=3, k::Integer=3)
    #println("quantum_entropy general case")
    return QuantumEntropy(X, m, k)
end

function quantum_entropy(X::MatrixOrConstant, m::Integer=0, k::Integer=0)
    #println("quantum_entropy constant X")
    return -quantum_relative_entropy(X, Matrix(1.0*I, size(X)))
end

function conic_form!(atom::QuantumEntropy, unique_conic_forms)
    if !has_conic_form(unique_conic_forms, atom)
        X = atom.children[1]
        m = atom.m
        k = atom.k
        n = size(X)[1]
        eye = Matrix(1.0*I, n, n)

        conic_form!(X ⪰ 0, unique_conic_forms)

        is_complex = sign(X) == ComplexSign()
        if is_complex
            τ = ComplexVariable(n, n)
        else
            τ = Variable(n, n)
        end
        conic_form!(τ in RelativeEntropyEpiCone(X, eye, m, k), unique_conic_forms)

        # It's already a real mathematically, but need to make it a real type.
        τ = real(-tr(τ))
        cache_conic_form!(unique_conic_forms, atom, minimize(τ))
    end
    return get_conic_form(unique_conic_forms, atom)
end


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

mutable struct IoK8sApiAppsV1beta2DaemonSetUpdateStrategy <: SwaggerModel
    rollingUpdate::Any # spec type: Union{ Nothing, IoK8sApiAppsV1beta2RollingUpdateDaemonSet } # spec name: rollingUpdate
    type::Any # spec type: Union{ Nothing, String } # spec name: type

    function IoK8sApiAppsV1beta2DaemonSetUpdateStrategy(;rollingUpdate=nothing, type=nothing)
        o = new()
        validate_property(IoK8sApiAppsV1beta2DaemonSetUpdateStrategy, Symbol("rollingUpdate"), rollingUpdate)
        setfield!(o, Symbol("rollingUpdate"), rollingUpdate)
        validate_property(IoK8sApiAppsV1beta2DaemonSetUpdateStrategy, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        o
    end
end # type IoK8sApiAppsV1beta2DaemonSetUpdateStrategy

const _property_map_IoK8sApiAppsV1beta2DaemonSetUpdateStrategy = Dict{Symbol,Symbol}(Symbol("rollingUpdate")=>Symbol("rollingUpdate"), Symbol("type")=>Symbol("type"))
const _property_types_IoK8sApiAppsV1beta2DaemonSetUpdateStrategy = Dict{Symbol,String}(Symbol("rollingUpdate")=>"IoK8sApiAppsV1beta2RollingUpdateDaemonSet", Symbol("type")=>"String")
Base.propertynames(::Type{ IoK8sApiAppsV1beta2DaemonSetUpdateStrategy }) = collect(keys(_property_map_IoK8sApiAppsV1beta2DaemonSetUpdateStrategy))
Swagger.property_type(::Type{ IoK8sApiAppsV1beta2DaemonSetUpdateStrategy }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAppsV1beta2DaemonSetUpdateStrategy[name]))}
Swagger.field_name(::Type{ IoK8sApiAppsV1beta2DaemonSetUpdateStrategy }, property_name::Symbol) =  _property_map_IoK8sApiAppsV1beta2DaemonSetUpdateStrategy[property_name]

function check_required(o::IoK8sApiAppsV1beta2DaemonSetUpdateStrategy)
    true
end

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


text:
################################
###    Summary statistics    ###
################################
"""
	poisson_fixation(observed_values,λds, λdn)

Divergence sampling from Poisson distribution. The expected neutral and selected fixations are subset through their relative expected rates ([`fix_neut`](@ref), [`fix_neg_b`](@ref), [`fix_pos_sim`](@ref)). Empirical values are used are used to simulate the locus *L* along a branch of time *T* from which the expected *Ds* and *Dn* raw count estimated given the mutation rate (``\\mu``). Random number generation is used to subset samples arbitrarily given the success rate ``\\lambda`` in the distribution.

```math
\\mathbb{E}[D_N] = X \\in Poisson\\left(\\lambda = D \\times \\left[\\frac{\\mathbb{E}[D_+] + \\mathbb{E}[D_-]}{\\mathbb{E}[D_+] + \\mathbb{E}[D_-] + \\mathbb{E}[D_0]}\\right]\\right)
```
```math
\\mathbb{E}[D_S] = X \\in Poisson\\left(\\lambda = D \\times \\left[\\frac{\\mathbb{E}[D_0]}{\\mathbb{E}[D_+] + \\mathbb{E}[D_-] + \\mathbb{E}[D_0]}\\right]\\right)
```
# Arguments
 - `observed_values::Array`: Array containing the total observed divergence.
 - ` λds::Float64`: expected neutral fixations rate.
 - ` λdn::Float64`: expected selected fixations rate.
# Returns
 - `Array{Int64,1}` containing the expected count of neutral and selected fixations.

"""
function poisson_fixation(;observed_values::Vector{Int64}, λds::Vector{Float64}, λdn::Vector{Float64},λweak::Vector{Float64},λstrong::Vector{Float64})

	ds = @. λds / (λds + λdn) * observed_values
	dn = @. λdn / (λds + λdn) * observed_values
	dweak = @. λweak / (λds + λdn) * observed_values
	dstrong = @. λstrong / (λds + λdn) * observed_values

	sampled_ds     = pois_rand.(ds)
	sampled_dn     = pois_rand.(dn)
	sampled_weak   = pois_rand.(dweak)
	sampled_strong = pois_rand.(dstrong)

	alphas = @. [sampled_weak/sampled_dn sampled_strong/sampled_dn (sampled_weak+sampled_strong)/sampled_dn]

	out = alphas,sampled_dn, sampled_ds
	return out
end

"""
	poisson_polymorphism(observed_values,λps,λpn)

Polymorphism sampling from Poisson distributions. The total expected neutral and selected polimorphism are subset through the relative expected rates at the frequency spectrum ([`fix_neut`](@ref), [`sfs_neut`](@ref),). Empirical sfs are used to simulate the locus *L* along a branch of time *T* from which the expected *Ps* and *Pn* raw count are estimated given the mutation rate (``\\mu``). Random number generation is used to subset samples arbitrarily from the whole sfs given each frequency success rate ``\\lambda`` in the distribution.

The success rate managing the Poisson distribution by the observed count each frequency.  We considered both sampling variance and process variance is affecting the number of variable alleles we sample from SFS. This variance arises from the random mutation-fixation process along the branch. To incorporate this variance we do one sample per frequency-bin and use the total sampled variation and the SFS for the summary statistics.

```math
\\mathbb{E}[P_N] = \\sum_{x=0}^{x=1} X \\in Poisson\\left(\\lambda = SFS_{(x)} \\times \\left[\\frac{\\mathbb{E}[P_{+(x)}] + \\mathbb{E}[P_{-(x)}]}{\\mathbb{E}[P_{+(x)}] + \\mathbb{E}[P_{-(x)}] + \\mathbb{E}[P_{0(x)}]}\\right]\\right)
```

```math
\\mathbb{E}[P_S] = \\sum_{x=0}^{x=1} X \\in Poisson\\left(\\lambda = SFS_{(x)} \\times \\left[\\frac{\\mathbb{E}[P_{0(x)}]}{\\mathbb{E}[P_{+(x)}] + \\mathbb{E}[P_{-(x)}] + \\mathbb{E}[P_{0(x)}]}\\right]\\right)
```

# Arguments
 - `observed_values::Array{Int64,1}`: Array containing the total observed divergence.
 - ` λps::Array{Float64,1} `: expected neutral site frequency spectrum rate.
 - ` λpn::Array{Float64,1} `: expected selected site frequency spectrum rate.
# Returns
 - `Array{Int64,2}` containing the expected total count of neutral and selected polymorphism.

"""
function poisson_polymorphism(;observed_values::Vector{Float64}, λps::Matrix{Float64}, λpn::Matrix{Float64})

	# Neutral λ;
	λ1 = @. λps / (λps + λpn) * observed_values
	# Selected λ;
	λ2 = @. λpn / (λps + λpn) * observed_values

	# Replace negative and 0 rates values. Strange behaviour depending on θ and Γ values
	# Relative rates output NaN values due to 0 divisons.

	replace!(λ1,NaN=>1)
	replace!(λ2,NaN=>1)

	sampled_ps = pois_rand.(λ1)
	sampled_pn = pois_rand.(λ2)

	return (sampled_pn, sampled_ps)
end

"""
	sampled_alpha(observed_values,λds, λdn)

Ouput the expected values from the Poisson sampling process. Please check [`poisson_fixation`](@ref) and [`poisson_polymorphism`](@ref) to understand the samplingn process. α(x) is estimated through the expected values of Dn, Ds, Pn and Ps.

# Arguments
 - `param::parameters`: Array containing the total observed divergence.
 - `d::Array`: observed divergence.
 - `afs::Array`: observed polymorphism.
 - ` λdiv::Array{Float64,2}`: expected fixations rate.
 - ` λdiv::Array{Float64,2}`: expected site frequency spectrum rates.
# Returns
α_summaries,exp_dn,exp_ds,exp_pn,exp_ps,ssAlpha
 - `Array{Int64,2}` containing α(x) values.
 - `Array{Int64,1}` expected non-synonymous divergence.
 - `Array{Int64,1}` expected synonymous divergence.
 - `Array{Int64,1}` expected non-synonymous polymorphism.
 - `Array{Int64,1}` expected synonymous polymorphism.
 - `Array{Int64,1}` expected synonymous polymorphism.
 - `Array{Int64,1}` expected synonymous polymorphism.
 - `Array{Int64,2}` containing α(x) binned values.

	sampled_from_rates(gammaL,gammaH,ppos_l,ppos_h,observedData,nopos)

"""
function sampled_from_rates(models::SubArray,fs::Vector{Float64},d::Vector{Int64},neut::SubArray,sel::SubArray,dsdn::SubArray,output::String)

	ds             = dsdn[:,1]
	dn             = dsdn[:,2]
	dweak          = dsdn[:,3]
	dstrong        = dsdn[:,4]
	gn             = abs.(models[:,4])
	sh             = round.(models[:,end-1],digits=5)

	## Outputs
	alphas, exp_dn, exp_ds = poisson_fixation(observed_values=d,λds=ds,λdn=dn,λweak=dweak,λstrong=dstrong);

	exp_pn, exp_ps         = poisson_polymorphism(observed_values=fs,λps=permutedims(neut),λpn=permutedims(sel));

	## Alpha from expected values. Used as summary statistics
	α_summaries = @. round(1 - ((exp_ds/exp_dn) * (exp_pn/exp_ps)'),digits=5);

	expected_values = hcat(round.(alphas,digits=5),gn,sh,α_summaries);

	expected_values = filter_expected(expected_values);

	w(x,name) = CSV.write(name,DataFrame(x,:auto),delim='\t',header=false);
	w(expected_values,output)
end

"""
	summary_statistics(param::parameters,rates::JLD2.JLDFile,analysis_folder::String,summstat_size::Int64,replicas::Int64,bootstrap::Bool)

Estimate summary statistics using observed data and analytical rates. *analysis_folder* will check for the SFS and divergence file and will be used to output summary statistics

# Arguments
 - `param::parameters` : Mutable structure containing the models
 - `h5_file::String` : HDF5 containing solved models, fixation and polymorphic rates
 - `sfs::Vector{Matrix{Float64}}`: SFS matrix.
 - `divergence::Vector{Matrix{Float64}}` : divergence matrix.
 - `summstat_size::Int64` : Number of summary statistics
 - `replicas::Int64` : Number of bootstrap replicas
 - `bootstrap::Bool` : Boolean to perform or not bootstrapping
# Output
 - Summary statistics to ABC inference
"""
function summary_statistics(;param::parameters,h5_file::String,sfs::Vector,divergence::Vector,analysis_folder::String,summstat_size::Int64,bootstrap::Union{Bool,Int64}=false,threads::Bool=false)

	## Opening files
	if length(sfs) > 1 & bootstrap != false
		throw(ArgumentError("You have more than one SFS and divergence file. Please be sure you have on set of files to bootstrap manually your data."))
	end

	assertion_params(param);

	α,sfs,divergence = data_to_poisson(sfs,divergence,param.dac,bootstrap);

	if any(0 .∈ sfs) | any(0 .∈ divergence)
		throw(ArgumentError("Your SFS contains 0 values at the selected DACs or the divergence is 0. Please consider to bin the SFS and re-estimate the rates using the selected bin as sample the new sample size."))
	end

	# Open rates
	string_cutoff = "cutoff=[" * string(param.cutoff[1]) * "," * string(param.cutoff[end]) * "]";
	h         = jldopen(h5_file);
	tmp       = h[string(param.N)* "/" * string(param.n) * "/" * string_cutoff]

	# Subset index
	idx    = sample(1:size(tmp["models"],1),summstat_size,replace=false);
	models = view(Array(tmp["models"]),idx,:);
	dsdn   = view(Array(tmp["dsdn"]),idx,:);

	# Filtering polymorphic rate by dac
	n    = hcat(map(x -> view(tmp["neut"][x],:),param.dac)...);
	s    = hcat(map(x -> view(tmp["sel"][x],:),param.dac)...);

	neut = view(n,idx,:);
	sel  = view(s,idx,:);

	# Making summaries
	summ_output = analysis_folder .* "/summstat_" .* string.(1:size(sfs,1)) .* ".txt"

	if threads
		ThreadsX.map((x,y,z) -> sampled_from_rates(models,x,y,neut,sel,dsdn,z),sfs,divergence,summ_output);
	else
		progress_pmap((x,y,z) -> sampled_from_rates(models,x,y,neut,sel,dsdn,z),sfs,divergence,summ_output);
	end

	w(x,name) = CSV.write(name,DataFrame(x,:auto),delim='\t',header=false);
	α_output = analysis_folder * "/alphas_" .* string.(1:size(sfs,1)) .* ".txt";

	w.(α, α_output);

	return(summ_output);
end

function filter_expected(x::Matrix{Float64})

	replace!(x, -Inf=>NaN)
	x = x[vec(.!any(isnan.(x),dims=2)),:]
	x = x[(x[:,3] .< 1),:]
	x = x[(x[:,1] .> 0),:]

	return(x)
end

function pol_correction!(sfs_all::Vector{Matrix{Float64}},sfs_in::Vector{Matrix{Float64}};column::Vector{Int}=[2])

	pn_all = map(x->sum(x[:,column]),sfs_all)

	sfs_pn = map(x -> sum(x[:,column]),sfs_in)

	ratio_all_sub = map( (x,y) -> x/y,pn_all,sfs_pn)

	map( (x,y) -> y[:,column] .= x * y[:,column],ratio_all_sub,sfs_in)
end


text:
using BenchmarkTools

function foo1()
    s = 0
    for i in 1:1e5
        s += i
    end
end

@benchmark foo1()


text:
using DelimitedFiles

filename = "./log-conv-CNS-1D-nosmooth-CFL.5"
file = readdlm(filename)

num_scheme = 2
Narr = [2;3]
Karr = [50;100;200;400;800]
block_size = 4
num_lines = size(file,1)
num_block = div(num_lines,block_size)
L1_err_dict = Dict()
for b in 1:num_block
    l = 1+block_size*(b-1)
    N = parse(Int,file[l,3][1])
    K = file[l,6]
    L1_err = file[l+1,4]
    @show N,K
    @show L1_err
    if !haskey(L1_err_dict,(N,K))
        L1_err_dict[(N,K)] = [L1_err]
    else
        push!(L1_err_dict[(N,K)],L1_err)
    end
end

num_row = length(Karr)-1
num_col = 2*length(Narr)
conv_table = zeros(num_row,num_col)
for i = 1:num_row
    for j = 1:num_col
        K = Karr[i]
        N = Narr[div(j-1,num_scheme)+1]
        s = mod1(j,num_scheme)
        conv_table[i,j] = log2(L1_err_dict[(N,K)][s]/L1_err_dict[(N,2*K)][s])
    end
end

err_table = zeros(length(Karr),2*length(Narr))
for i = 1:length(Karr)
    for j = 1:length(Narr)
        K = Karr[i]
        N = Narr[j]
        err_table[i,2*(j-1)+1] = L1_err_dict[(N,K)][1]
        err_table[i,2*(j-1)+2] = L1_err_dict[(N,K)][2]
    end
end

max_stars_repo_name: yiminllin/ESDG-PosLimit

max_issues_repo_name: yiminllin/ESDG-PosLimit

max_forks_repo_name: yiminllin/ESDG-PosLimit


text:
using InfiniteLinearAlgebra, BlockBandedMatrices, BlockArrays, BandedMatrices, InfiniteArrays, FillArrays, LazyArrays, Test, MatrixFactorizations, LinearAlgebra, Random
import InfiniteLinearAlgebra: qltail, toeptail, tailiterate , tailiterate!, tail_de, ql_X!,
                    InfToeplitz, PertToeplitz, TriToeplitz, InfBandedMatrix,
                    rightasymptotics, QLHessenberg, ConstRows, PertConstRows, BandedToeplitzLayout, PertToeplitzLayout
import Base: BroadcastStyle
import BlockArrays: _BlockArray
import BlockBandedMatrices: isblockbanded, _BlockBandedMatrix
import MatrixFactorizations: QLPackedQ
import BandedMatrices: bandeddata, _BandedMatrix, BandedStyle
import LazyArrays: colsupport, ApplyStyle, MemoryLayout, ApplyLayout, LazyArrayStyle

@testset "∞-Toeplitz and Pert-Toeplitz" begin
    A = BandedMatrix(1 => Fill(2im,∞), 2 => Fill(-1,∞), 3 => Fill(2,∞), -2 => Fill(-4,∞), -3 => Fill(-2im,∞))
    @test A isa InfToeplitz
    @test MemoryLayout(typeof(A.data)) == ConstRows()
    @test MemoryLayout(typeof(A)) == BandedToeplitzLayout()
    V = view(A,:,3:∞)
    @test MemoryLayout(typeof(bandeddata(V))) == ConstRows()
    @test MemoryLayout(typeof(V)) == BandedToeplitzLayout()

    @test BandedMatrix(V) isa InfToeplitz
    @test A[:,3:end] isa InfToeplitz

    A = BandedMatrix(-2 => Vcat(Float64[], Fill(1/4,∞)), 0 => Vcat([1.0+im,2,3],Fill(0,∞)), 1 => Vcat(Float64[], Fill(1,∞)))
    @test A isa PertToeplitz
    @test MemoryLayout(typeof(A)) == PertToeplitzLayout()
    V = view(A,2:∞,2:∞)
    @test MemoryLayout(typeof(V)) == PertToeplitzLayout()
    @test BandedMatrix(V) isa PertToeplitz
    @test A[2:∞,2:∞] isa PertToeplitz
end

@testset "Algebra" begin
    @testset "BandedMatrix" begin
        A = BandedMatrix(-3 => Fill(7/10,∞), -2 => 1:∞, 1 => Fill(2im,∞))
        @test A isa BandedMatrix{ComplexF64}
        @test A[1:10,1:10] == diagm(-3 => Fill(7/10,7), -2 => 1:8, 1 => Fill(2im,9))

        A = BandedMatrix(0 => Vcat([1,2,3],Zeros(∞)), 1 => Vcat(1, Zeros(∞)))
        @test A[1,2] == 1

        A = BandedMatrix(-3 => Fill(7/10,∞), -2 => Fill(1,∞), 1 => Fill(2im,∞))
        Ac = BandedMatrix(A')
        At = BandedMatrix(transpose(A))
        @test Ac[1:10,1:10] ≈ (A')[1:10,1:10] ≈ A[1:10,1:10]'
        @test At[1:10,1:10] ≈ transpose(A)[1:10,1:10] ≈ transpose(A[1:10,1:10])

        A = BandedMatrix(-1 => Vcat(Float64[], Fill(1/4,∞)), 0 => Vcat([1.0+im],Fill(0,∞)), 1 => Vcat(Float64[], Fill(1,∞)))
        @test MemoryLayout(typeof(view(A.data,:,1:10))) == ApplyLayout{typeof(hcat)}()
        Ac = BandedMatrix(A')
        At = BandedMatrix(transpose(A))
        @test Ac[1:10,1:10] ≈ (A')[1:10,1:10] ≈ A[1:10,1:10]'
        @test At[1:10,1:10] ≈ transpose(A)[1:10,1:10] ≈ transpose(A[1:10,1:10])

        A = BandedMatrix(-2 => Vcat(Float64[], Fill(1/4,∞)), 0 => Vcat([1.0+im,2,3],Fill(0,∞)), 1 => Vcat(Float64[], Fill(1,∞)))
        Ac = BandedMatrix(A')
        At = BandedMatrix(transpose(A))
        @test Ac[1:10,1:10] ≈ (A')[1:10,1:10] ≈ A[1:10,1:10]'
        @test At[1:10,1:10] ≈ transpose(A)[1:10,1:10] ≈ transpose(A[1:10,1:10])

        A = _BandedMatrix(Fill(1,4,∞),∞,1,2)
        @test A*A isa ApplyArray
        @test (A^2)[1:10,1:10] == (A*A)[1:10,1:10] == (A[1:100,1:100]^2)[1:10,1:10]
        @test (A^3)[1:10,1:10] == (A*A*A)[1:10,1:10] == (A[1:100,1:100]^3)[1:10,1:10]
    end

    @testset "BlockTridiagonal" begin
        A = BlockTridiagonal(Vcat([fill(1.0,2,1),Matrix(1.0I,2,2),Matrix(1.0I,2,2),Matrix(1.0I,2,2)],Fill(Matrix(1.0I,2,2), ∞)),
                            Vcat([zeros(1,1)], Fill(zeros(2,2), ∞)),
                            Vcat([fill(1.0,1,2),Matrix(1.0I,2,2)], Fill(Matrix(1.0I,2,2), ∞)))

        @test A isa InfiniteLinearAlgebra.BlockTriPertToeplitz
        @test isblockbanded(A)

        @test A[Block.(1:2),Block(1)] == A[1:3,1:1] == reshape([0.,1.,1.],3,1)

        @test BlockBandedMatrix(A)[1:100,1:100] == BlockBandedMatrix(A,(2,1))[1:100,1:100] == BlockBandedMatrix(A,(1,1))[1:100,1:100] == A[1:100,1:100]

        @test (A - I)[1:100,1:100] == A[1:100,1:100]-I
        @test (A + I)[1:100,1:100] == A[1:100,1:100]+I
        @test (I + A)[1:100,1:100] == I+A[1:100,1:100]
        @test (I - A)[1:100,1:100] == I-A[1:100,1:100]
    end

    @testset "Fill" begin
        A = _BandedMatrix(Ones(1,∞),∞,-1,1)
        @test 1.0 .* A isa BandedMatrix{Float64,<:Fill}
        @test_skip Ones(∞) .* A
        @test 2.0 .* A isa BandedMatrix{Float64,<:Fill}
        @test A .* 2.0 isa BandedMatrix{Float64,<:Fill}
        @test Eye(∞)*A isa BandedMatrix{Float64,<:Fill}
        @test A*Eye(∞) isa BandedMatrix{Float64,<:Fill}
    end

    @testset "Banded Broadast" begin
        A = _BandedMatrix((1:∞)',∞,-1,1)
        @test 2.0 .* A isa BandedMatrix{Float64,<:Adjoint}
        @test A .* 2.0 isa BandedMatrix{Float64,<:Adjoint}
        @test Eye(∞)*A isa BandedMatrix{Float64,<:Adjoint}
        @test A*Eye(∞) isa BandedMatrix{Float64,<:Adjoint}
        A = _BandedMatrix(Vcat((1:∞)',Ones(1,∞)),∞,0,1)
        @test 2.0 .* A isa BandedMatrix
        @test A .* 2.0 isa BandedMatrix
        @test Eye(∞) * A isa BandedMatrix
        @test A * Eye(∞) isa BandedMatrix
        b = 1:∞
        @test BroadcastStyle(typeof(b)) isa LazyArrayStyle{1}
        @test BroadcastStyle(typeof(A)) isa BandedStyle
        @test BroadcastStyle(LazyArrayStyle{1}(), BandedStyle()) isa LazyArrayStyle{2}
        @test BroadcastStyle(LazyArrayStyle{2}(), BandedStyle()) isa LazyArrayStyle{2}
        @test bandwidths(b .* A) == (0,1)

        @test colsupport(b.*A, 1) == 1:1
        @test Base.replace_in_print_matrix(b.*A, 2,1,"0.0") == " ⋅ "
        @test bandwidths(A .* b) == (0,1)
        @test A .* b' isa BroadcastArray
        @test bandwidths(A .* b') == bandwidths(A .* b')
        @test colsupport(A .* b', 3) == 2:3

        A = _BandedMatrix(Ones{Int}(1,∞),∞,0,0)'
        B = _BandedMatrix((-2:-2:-∞)', ∞,-1,1)
        C = Diagonal( 2 ./ (1:2:∞))
        @test A*(B*C) isa MulMatrix
        @test bandwidths(A*(B*C)) == (-1,1)
    end

    @testset "Triangle OP recurrences" begin
        # mortar((n -> 1:n).(1:∞))
    end
    # Multivariate OPs Corollary (3)
    # n = 5
    # BlockTridiagonal(Zeros.(1:∞,2:∞),
    #         (n -> Diagonal(((n+2).+(0:n)))/ (2n + 2)).(0:∞),
    #         Zeros.(2:∞,1:∞))
end

include("test_hessenbergq.jl")
include("test_infql.jl")
include("test_infqr.jl")


text:
"""
    CatmullRom

Interpolate using the centripetal parameterization of Catmull-Rom splines.
This is the form that does not crimp, provides taughtness, and looks good.
"""
module CatmullRom

export catmullrom,               # populates with points placed between those given
       catmullrom_by_arclength,  # populates with points placed between those given relative to arclength
       close_seq!,               # ensures a point sequence is closed, precisely
       floatvecs                 # obtain a vector of vectors of floats

using LinearAlgebra: dot, norm, normalize
using Polynomials: Poly, polyval

include("utilities.jl")

# The sorts of sequences understood to hold point coordinates
#     defines `Points`m `npoints(Points)`, `ncoords(Points)`
include("presence/pointsequences.jl")

# suggest outermost two points for Catmull-Rom spline sequence
const ReflectionScale = 1.0
const ArcpointsMin = 8
const ArcpointsMax = 64

include("fewpoints/outside.jl")

# centripetal Catmull-Rom compution
include("centripetal/catmullrom.jl")
include("centripetal/arcbased.jl")

# extrapolation using low-order interpolatory fits
include("fewpoints/twopoints.jl")    # linear fit, reflection
include("fewpoints/threepoints.jl")  # thiele3 rational fit, general quadratic fit
include("fewpoints/fourpoints.jl")   # thiele4 rational fit (preferred)

# myunzip, anglesep
include("presence/support.jl")

end # CatmullRom


text:
# test the MDP module

include("./script.jl")
import .MDPSimulator

function test1()
    # graph:
    # a <-> b <-> c <-> d <-> e
    # terminals: a (10) and e (1)
    # actions: left or right
    states = ["a", "b", "c", "d", "e"]
    actions = ["L", "R"]
    terminals = ["a", "e"]
    terminal_vals = [10.0, 1.0]
    gamma = 0.2
    graph = MDPSimulator.createMDPGraph(states, actions, terminals, terminal_vals, gamma)
    MDPSimulator.addMDPGraphTransition!(graph, "b", "L", "a", prob=1.0)
    MDPSimulator.addMDPGraphTransition!(graph, "b", "R", "c", prob=1.0)
    MDPSimulator.addMDPGraphTransition!(graph, "c", "L", "b", prob=1.0)
    MDPSimulator.addMDPGraphTransition!(graph, "c", "R", "d", prob=1.0)
    MDPSimulator.addMDPGraphTransition!(graph, "d", "L", "c", prob=1.0)
    MDPSimulator.addMDPGraphTransition!(graph, "d", "R", "e", prob=1.0)
    MDPSimulator.runMDP!(graph, 100)
    MDPSimulator.displayMDP(graph)
end

function test2()
    # graph:
    # A <-> C
    # A <-> B
    # C <-> B
    # terminals: None
    # actions: C or CC
    states = ["A", "B", "C"]
    actions = ["C", "CC"]
    gamma = 0.5
    terminals::Array{String} = []
    terminal_vals::Array{AbstractFloat} = []
    graph = MDPSimulator.createMDPGraph(states, actions, terminals, terminal_vals, gamma)
    MDPSimulator.addMDPGraphTransition!(graph, "A", "C", "B", prob=1.0, reward=0.0)
    MDPSimulator.addMDPGraphTransition!(graph, "A", "CC", "C", prob=1.0, reward=-2.0)
    MDPSimulator.addMDPGraphTransition!(graph, "B", "C", "A", prob=0.4, reward=-1.0)
    MDPSimulator.addMDPGraphTransition!(graph, "B", "C", "C", prob=0.6, reward=2.0)
    MDPSimulator.addMDPGraphTransition!(graph, "B", "CC", "A", prob=0.6, reward=2.0)
    MDPSimulator.addMDPGraphTransition!(graph, "B", "CC", "C", prob=0.4, reward=-1.0)
    MDPSimulator.addMDPGraphTransition!(graph, "C", "C", "A", prob=0.6, reward=2.0)
    MDPSimulator.addMDPGraphTransition!(graph, "C", "C", "B", prob=0.4, reward=2.0)
    MDPSimulator.addMDPGraphTransition!(graph, "C", "CC", "A", prob=0.4, reward=2.0)
    MDPSimulator.addMDPGraphTransition!(graph, "C", "CC", "B", prob=0.6, reward=0.0)
    MDPSimulator.runMDP!(graph, 100)
    MDPSimulator.displayMDP(graph)
end

function test3()
    # graph:
    # A -> B,C
    # B -> D
    # C -> B,E
    # D -> F
    # E -> D,F
    # terminal: None
    # action: g
    states = ["A", "B", "C", "D", "E", "F"]
    actions = ["g"]
    gamma = 0.5
    terminals::Array{String} = []
    terminal_vals::Array{AbstractFloat} = []
    graph = MDPSimulator.createMDPGraph(states, actions, terminals, terminal_vals, gamma)
    MDPSimulator.addMDPGraphTransition!(graph, "A", "g", "B", prob=0.5, reward=1.0)
    MDPSimulator.addMDPGraphTransition!(graph, "A", "g", "C", prob=0.5, reward=1.0)
    MDPSimulator.addMDPGraphTransition!(graph, "B", "g", "D", prob=1.0, reward=1.0)
    MDPSimulator.addMDPGraphTransition!(graph, "C", "g", "B", prob=0.5, reward=1.0)
    MDPSimulator.addMDPGraphTransition!(graph, "C", "g", "E", prob=0.5, reward=1.0)
    MDPSimulator.addMDPGraphTransition!(graph, "D", "g", "F", prob=1.0, reward=1.0)
    MDPSimulator.addMDPGraphTransition!(graph, "E", "g", "D", prob=0.5, reward=1.0)
    MDPSimulator.addMDPGraphTransition!(graph, "E", "g", "F", prob=0.5, reward=1.0)
    MDPSimulator.addMDPGraphTransition!(graph, "F", "g", "F", prob=1.0, reward=0.0)
    MDPSimulator.runMDP!(graph, 100)
    MDPSimulator.displayMDP(graph)
end

println("Test 1")
test1()
println()

println("Test 2")
test2()
println()

println("Test 3")
test3()
println()

max_stars_repo_name: teamclouday/MDPSimulator

max_issues_repo_name: teamclouday/MDPSimulator

max_forks_repo_name: teamclouday/MDPSimulator


text:
module KmerAnalysisMakie

export
    spectracn

using Observables
using AbstractPlotting

import KmerAnalysis: KmerFrequencySpectra

"""
    _find_plottable_subset(spec::KmerFrequencySpectra{1})

Utility method that automatically finds the highest k-mer frequency in a spectra
with a non-zero count, and returns a view of the spectra data that only goes up
to that detected frequeny.

Currently this is useful to have Makie primitives like `scatter` and `barplot`
not plot a large portion of the graph where there's basically nothing interesting.
"""
function _find_plottable_subset(spec::KmerFrequencySpectra{1})
    return view(spec.data, firstindex(spec.data):findprev(x -> x > 0, spec.data, lastindex(spec.data)))
end

"""
    _find_plottable_subset(spec::KmerFrequencySpectra{2})

Utility method that automatically finds the highest k-mer frequency in both of
a 2D spectras dimensions with a non-zero count, and returns a sub-view of the
spectra data that only goes up to those detected frequencies.

Currently this is useful to have Makie primitives like `heatmap` and `contour`
not plot a large portion of the graph where there's basically nothing interesting.
"""
function _find_plottable_subset(spec::KmerFrequencySpectra{2})
    mat = spec.data
    furthest_row = 0
    furthest_col = 0
    for col in 1:257
        for row in 1:257
            v = mat[row, col]
            furthest_row = ifelse((row > furthest_row) & !iszero(v), row, furthest_row)
            furthest_col = ifelse((col > furthest_col) & !iszero(v), col, furthest_col)
        end
    end
    return view(mat, 1:furthest_row, 1:furthest_col)
end

include(joinpath("type_recipes", "kmer_spectra.jl"))
include(joinpath("full_recipes", "spectracn.jl"))

end # module


text:
#
# structure factor test
# @author Tobias Weber <tobias.weber@tum.de>
# @date nov-2017
# @license GPLv2
#

include("sgs.jl")
include("slens.jl")

function structfact(atoms, G)
	F = 0 + 0im

	for atom in atoms
		F += exp(2.*pi * 1im * dot(G, atom))
	end

	return F
end

function has_pos(allpos, pos)
	for curpos in allpos
		if(curpos == pos)
			return true
		end
	end

	return false
end

function generate_pos(pos, trafos)
	allpos = []

	for trafo in trafos
		rot = trafo[1:3, 1:3]
		trans = trafo[1:3, 4:4]

		newpos = rot * pos + trans

#		for idx in range(1, length(newpos))
#			while(newpos[idx] < 0.)
#				newpos[idx] += 1.
#			end
#			while(newpos[idx] > 1.)
#				newpos[idx] -= 1.
#			end
#		end

		if !has_pos(allpos, newpos)
			push!(allpos, newpos)
		end
	end

	#println(allpos)
	return allpos
end

atom_1 = "Mn"
atom_2 = "Si"
atompos_1 = [ 0.14; 0.14; 0.14 ]
atompos_2 = [ 0.85; 0.85; 0.85 ]
sgroup = "P2_13"

slen_1 = slens[atom_1]["coh"]
slen_2 = slens[atom_2]["coh"]
trafos = sgs[sgroup]["trafos"]

allatompos_1 = generate_pos(atompos_1, trafos)
allatompos_2 = generate_pos(atompos_2, trafos)

G = [ 1.; 1.; 0. ]

F_1 = structfact(allatompos_1, G) * slen_1
F_2 = structfact(allatompos_2, G) * slen_2
F = F_1 + F_2

println(atom_1, " atoms: ", allatompos_1)
println(atom_2, " atoms: ", allatompos_2)

println("G = ", G, ", F = ", F, ", |F| = ", sqrt(F*conj(F)))

max_stars_repo_name: t-weber/libcrystal

max_issues_repo_name: t-weber/libcrystal

max_forks_repo_name: t-weber/libcrystal


text:
# Julia wrapper for header: /usr/include/curl/curl.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

@c Cint curl_strequal (Ptr{UInt8}, Ptr{UInt8}) libcurl
@c Cint curl_strnequal (Ptr{UInt8}, Ptr{UInt8}, size_t) libcurl
@c CURLFORMcode curl_formadd (Ptr{Ptr{Void}}, Ptr{Ptr{Void}}) libcurl
@c Cint curl_formget (Ptr{Void}, Ptr{Void}, curl_formget_callback) libcurl
@c Void curl_formfree (Ptr{Void},) libcurl
@c Ptr{UInt8} curl_getenv (Ptr{UInt8},) libcurl
@c Ptr{UInt8} curl_version () libcurl
@c Ptr{UInt8} curl_easy_escape (Ptr{CURL}, Ptr{UInt8}, Cint) libcurl
@c Ptr{UInt8} curl_escape (Ptr{UInt8}, Cint) libcurl
@c Ptr{UInt8} curl_easy_unescape (Ptr{CURL}, Ptr{UInt8}, Cint, Ptr{Cint}) libcurl
@c Ptr{UInt8} curl_unescape (Ptr{UInt8}, Cint) libcurl
@c Void curl_free (Ptr{Void},) libcurl
@c CURLcode curl_global_init (Cint,) libcurl
@c CURLcode curl_global_init_mem (Cint, curl_malloc_callback, curl_free_callback, curl_realloc_callback, curl_strdup_callback, curl_calloc_callback) libcurl
@c Void curl_global_cleanup () libcurl
@c Ptr{Void} curl_slist_append (Ptr{Void}, Ptr{UInt8}) libcurl
@c Void curl_slist_free_all (Ptr{Void},) libcurl
@c time_t curl_getdate (Ptr{UInt8}, Ptr{time_t}) libcurl
@c Ptr{CURLSH} curl_share_init () libcurl
@c CURLSHcode curl_share_setopt (Ptr{CURLSH}, CURLSHoption) libcurl
@c CURLSHcode curl_share_cleanup (Ptr{CURLSH},) libcurl
@c Ptr{curl_version_info_data} curl_version_info (CURLversion,) libcurl
@c Ptr{UInt8} curl_easy_strerror (CURLcode,) libcurl
@c Ptr{UInt8} curl_share_strerror (CURLSHcode,) libcurl
@c CURLcode curl_easy_pause (Ptr{CURL}, Cint) libcurl
@c Ptr{CURL} curl_easy_init () libcurl
@c CURLcode curl_easy_setopt (Ptr{CURL}, CURLoption) libcurl
@c CURLcode curl_easy_perform (Ptr{CURL},) libcurl
@c Void curl_easy_cleanup (Ptr{CURL},) libcurl
@c CURLcode curl_easy_getinfo (Ptr{CURL}, CURLINFO) libcurl
@c Ptr{CURL} curl_easy_duphandle (Ptr{CURL},) libcurl
@c Void curl_easy_reset (Ptr{CURL},) libcurl
@c CURLcode curl_easy_recv (Ptr{CURL}, Ptr{Void}, size_t, Ptr{size_t}) libcurl
@c CURLcode curl_easy_send (Ptr{CURL}, Ptr{Void}, size_t, Ptr{size_t}) libcurl
@c Ptr{CURLM} curl_multi_init () libcurl
@c CURLMcode curl_multi_add_handle (Ptr{CURLM}, Ptr{CURL}) libcurl
@c CURLMcode curl_multi_remove_handle (Ptr{CURLM}, Ptr{CURL}) libcurl
@c CURLMcode curl_multi_fdset (Ptr{CURLM}, Ptr{fd_set}, Ptr{fd_set}, Ptr{fd_set}, Ptr{Cint}) libcurl
@c CURLMcode curl_multi_perform (Ptr{CURLM}, Ptr{Cint}) libcurl
@c CURLMcode curl_multi_cleanup (Ptr{CURLM},) libcurl
@c Ptr{CURLMsg} curl_multi_info_read (Ptr{CURLM}, Ptr{Cint}) libcurl
@c Ptr{UInt8} curl_multi_strerror (CURLMcode,) libcurl
@c CURLMcode curl_multi_socket (Ptr{CURLM}, curl_socket_t, Ptr{Cint}) libcurl
@c CURLMcode curl_multi_socket_action (Ptr{CURLM}, curl_socket_t, Cint, Ptr{Cint}) libcurl
@c CURLMcode curl_multi_socket_all (Ptr{CURLM}, Ptr{Cint}) libcurl
@c CURLMcode curl_multi_timeout (Ptr{CURLM}, Ptr{Cint}) libcurl
@c CURLMcode curl_multi_setopt (Ptr{CURLM}, CURLMoption) libcurl
@c CURLMcode curl_multi_assign (Ptr{CURLM}, curl_socket_t, Ptr{Void}) libcurl


text:
using Test
using Convex: MOI
import SCS

# https://optimization.mccormick.northwestern.edu/index.php/Quadratic_programming#Numerical_example
Q = Float64[[3 1]; [1 1]]
p = Float64[1, 6]
G = Float64[[-2 -3]; [-1 0]; [0 -1]]
h = Float64[-4, 0, 0]
A = Float64[1 1]
b = Float64[1.5]

@testset "qp.jl: Solve QP Test 1" begin
    x, λ, ν = solve_qp(Q, p, G, h, A, b)
    @test size(x) == (2,)
    @test size(λ) == (3,)
    @test size(ν) == (1,)
    @test x ≈ [0.5, 1.0] atol = 1e-5
    @test λ ≈ [4.0, 0.0, 0.0] atol = 2e-5
    @test ν ≈ [4.5] atol = 3e-5
end

@testset "qp.jl: Solve QP Test 2" begin
    optimizer = SCS.Optimizer()
    MOI.set(optimizer, MOI.RawParameter("rho_x"), 1e-5)
    x, λ, ν = solve_qp(Q, p, G, h, A, b, optimizer, verbose = false)
    @test size(x) == (2,)
    @test size(λ) == (3,)
    @test size(ν) == (1,)
    @test x ≈ [0.5, 1.0] atol = 1e-10
    @test λ ≈ [4.0, 0.0, 0.0] atol = 5e-10
    @test ν ≈ [4.5] atol = 5e-10
end

@testset "qp.jl: Solve QP Test 3" begin
    _, λ, ν = solve_qp(Q, p, G, h, Float64[], Float64[])
    @test isempty(ν)
    @test !isempty(λ)
end

@testset "qp.jl: Solve QP Test 4" begin
    _, λ, ν = solve_qp(Q, p, Float64[], Float64[], A, b)
    @test isempty(λ)
    @test !isempty(ν)
end

@testset "qp.jl: Solve Batch QP" begin
    Qs = reshape(repeat(Q, 1, 5), 2, 2, 5)
    ps = repeat(p, 1, 5)
    Gs = reshape(repeat(G, 1, 5), 3, 2, 5)
    hs = repeat(h, 1, 5)
    As = reshape(repeat(A, 1, 5), 1, 2, 5)
    bs = repeat(b, 1, 5)
    X, λs, νs = solve_qp_batch(Qs, ps, Gs, hs, As, bs)
    @test size(X) == (2, 5)
    @test size(λs) == (3, 5)
    @test size(νs) == (1, 5)
    @test X ≈ repeat([0.5, 1.0], 1, 5) atol = 5e-6
    @test λs ≈ repeat([4.0, 0.0, 0.0], 1, 5) atol = 5e-5
    @test νs ≈ repeat([4.5], 1, 5) atol = 1e-4
end


text:
# example of neural net functions in Econometrics.jl
using Econometrics, Glob

# generate draws from linear regression model, and
# fitted coefficients from correct model, plus
# quadratic and cubic models (irrelevant regressors)
# and 5 pure noise statistics
function make_simdata(reps=100000)
    n = 30
    simdata = zeros(reps, 35+6)
    for rep = 1:reps
        # draw the regressors
        x = randn(n,4)
        z = [ones(n,1) x]
        # draw the parameters from prior
        b = randn(5,1)
        sig = exp.(randn(1))
        # generate dependent variable
        y = z*b + sig.*randn(n,1)
        # linear model
        bhat1 = z\y
        uhat = y-z*bhat1
        sighat1 = sqrt.(uhat'*uhat/(n-size(z,2)))
        # quadratic model
        z = [ones(n,1) x 0.1*x.^2.0]
        bhat2 = z\y
        uhat = y-z*bhat2
        sighat2 = sqrt.(uhat'*uhat/(n-size(z,2)))
        # cubic model
        z = [ones(n,1) x 0.1*x.^2.0 0.01*x.^3.0]
        bhat3 = z\y
        uhat = y-z*bhat3
        sighat3 = sqrt.(uhat'*uhat/(n-size(z,2)))
        # pure noise
        z = randn(1,5)
        # assemble:
        simdata[rep,:] = [b' log.(sig) bhat1' log.(sighat1) bhat2' log.(sighat2) bhat3' log.(sighat3) z]
    end
    return simdata
end

# fit a neural net to the linear model data, and show influence of statistics
function main()
    data = make_simdata()
    noutputs = 6
    trainsize = 80000
    savefile = "olsnet"
    layerconfig = [20, 10, 0, 0]
    epochs = 30
    TrainNet(data, trainsize, noutputs, layerconfig, 512, epochs, savefile)
    Y = data[trainsize+1:end,1:noutputs]
    olsfit = data[trainsize+1:end,(noutputs+1):(2*noutputs)]
    title = "linear regression example"
    params = ["constant", "x1","x2","x3","x4","x5"]
    fit = AnalyzeNet(savefile, epochs, data, trainsize, noutputs, title=title, params=params, doplot=true);
    rm.(glob("olsnet-*"))
end
main();


text:
struct EnumerativeBatchQs <: MultiObjectiveBatchQs
    model::SVDD.OCClassifier
    inf_measure::SequentialPQs
    rep_measure::Function
    div_measure::Function
    normalization::Function
    k::Int
    λ_inf::Float64
    λ_rep::Float64
    λ_div::Float64

    function EnumerativeBatchQs(model::SVDD.OCClassifier, informativeness::SequentialPQs; representativeness::Symbol, diversity::Symbol,
        k::Int, λ_inf::T1 where T1<:Real=1, λ_rep::T2 where T2 <: Real=1, λ_div::T3 where T3 <: Real=1)::EnumerativeBatchQs
        (k < 1) && throw(ArgumentError("Invalid batch size k=$(k)."))

        representativeness_measure = get_rep_measure(representativeness)
        diversity_measure = get_div_measure(diversity)
        return new(model, informativeness, representativeness_measure, diversity_measure, min_max_normalization, k, λ_inf, λ_rep, λ_div)
    end
end

"""
Select best batch with weighted sum of requirements.

Enumerative selection: Compute weighted sum for all possible candidates.
Warning: This may take a while.
"""
function select_batch(qs::EnumerativeBatchQs, x::Array{T, 2}, labels::Dict{Symbol, Vector{Int}}, candidate_indices::Vector{Int})::Vector{Int} where T <: Real
    num_observations = length(candidate_indices)
    if num_observations <= qs.k
        return candidate_indices
    end
    #informativeness
    inf_scores_normalized = qs.normalization(qs_score(qs.inf_measure, x, labels)[candidate_indices])

    # representativeness
    rep_scores_normalized = qs.normalization(qs.rep_measure(qs.model, x, labels, candidate_indices))

    best_batch = Vector{Int}()
    best_score = -Inf

    for batch in subsets(1:num_observations, qs.k)
        inf_score = sum(inf_scores_normalized[batch])/qs.k
        rep_score = sum(rep_scores_normalized[batch])/qs.k
        div_score = qs.div_measure(qs.model, x, candidate_indices[batch])
        combined_score = qs.λ_inf * inf_score + qs.λ_rep * rep_score + qs.λ_div * div_score

        if combined_score > best_score
            best_score = combined_score
            best_batch = candidate_indices[batch]
        end
    end
    return best_batch
end


text:
"""
    perm,invp = find_ped_order(pedlist [,start])

Traverses the pedigree and finds a permutation so that the ancestors precede their progeny.
The resulting vector `perm` associates the current ID with the new ID i.e., `perm[newid]=oldid`.
The vector `invp` has the inverse permutation i.e. `invp[oldid]=newid`.
By default, the function takes each animal as the starting animal.
When supplying an optional vector `start` with starting IDs, the tracing starts only with these progeny.
In this case, `perm` and `invp` will include `0` for non-tracked animals.
"""
function find_ped_order(pedlist::Matrix{Ti},start::Vector{Ti}=zeros(Ti,0)) where Ti<:Integer
   n::Ti = size(pedlist,2)
   perm = zeros(Ti,n)
   invp = zeros(Ti,n)
   newid = zero(Ti)
   if length(start)>0
      for k=1:length(start)
         i = start[k]
         if i<1 || i>n
            error("starting ID out of range")
         end
         newid = find_ped_order_within_family!(pedlist,i,newid,0,invp,perm)
      end
   else
      for i=1:n
         newid = find_ped_order_within_family!(pedlist,i,newid,0,invp,perm)
      end
   end
   return perm,invp
end

# assign new (permutaed) id using recursion within family
# perm[newid]  = origid
# invp[origid] = newid
function find_ped_order_within_family!(pedlist::Matrix{Ti},aid::Ti,lastaid::Ti,gen::Ti,invp::Vector{Ti},perm::Vector{Ti}) where Ti<:Integer
   sid = pedlist[1,aid]
   did = pedlist[2,aid]
   n::Ti = size(pedlist,2)
   newaid::Ti = lastaid

   # already registered in invp
   if invp[aid] > 0
      return newaid
   end
   # too many generations
   if gen>64
      error("too many generations back")
   end
   # not registered in invp
   if (sid<1||sid>n) && (did<1||did>n)
      # both unknown: increment new id

   elseif sid<1||sid>n
      # sire unknown: track the dam first; then register this animal
      newaid = find_ped_order_within_family!(pedlist,did,newaid,gen+1,invp,perm)
   elseif did<1||did>n
      # dam unknown: track the sire first; then register this animal
      newaid = find_ped_order_within_family!(pedlist,sid,newaid,gen+1,invp,perm)
   else
      # both known: just track both parents
      newaid = find_ped_order_within_family!(pedlist,sid,newaid,gen+1,invp,perm)
      newaid = find_ped_order_within_family!(pedlist,did,newaid,gen+1,invp,perm)
   end
   # not registered yet
   if invp[aid] < 1
      newaid = newaid + 1
      invp[aid] = newaid
      perm[newaid] = aid
   end
   return newaid
end

# a function return the renumbered ID for a parent
function get_new_parent_code(oldpid::Ti,n::Ti,pedlist::Matrix{Ti},invp::Vector{Ti}) where Ti<:Integer
   if oldpid<1
      # missing parent
      return 0
   elseif oldpid>n
      # unknown parent group
      return oldpid
   else
      # renumbered parent ID
      return invp[oldpid]
   end
end

# apply the permutation to pedigree
# perm[newid]  = origid
# invp[origid] = newid
"""
    permute_ped!(invp, pedlist [,idtable])

This function "renumbers" the animals by a permutation vector.
It updates the code in `pedlist` according to the permutation list in `invp` as `invp[oldid]=newid`, which is generated with some functions like `find_ped_order`.
The function also replaces IDs in `idtable` (optional).
A code for unknown-parent group, which is usually greater than the maximum ID of real animals, remains in the new pedigree.
"""
function permute_ped!(invp::Vector{Ti},pedlist::Matrix{Ti},idtable::Dict{String,Ti}=Dict{String,Ti}()) where Ti<:Integer
   # incomplete permutation vector
   n::Ti = length(invp)
   m::Ti = count((invp.>=1) .& (invp.<=n))
   if length(invp)!=size(pedlist,2)
      error("dimension mismatch between pedlist and invp")
   end
   if m!=n
      error("incomplete permutation vector(s); use ExtractPedigree() to extract the subset pedigree")
   end

   # swap the entry
   tmplist = zeros(Ti,2,n)
   for oldaid=1:n
      # animal
      newaid = invp[oldaid]
      # sire
      oldsid = pedlist[1,oldaid]
      newsid = get_new_parent_code(oldsid,n,pedlist,invp)
      # dam
      olddid = pedlist[2,oldaid]
      newdid = get_new_parent_code(olddid,n,pedlist,invp)
      # update
      tmplist[1:2,newaid] = [newsid; newdid]
   end
   pedlist[1:2,1:n] = tmplist[1:2,1:n]

   # change the ID table (optional)
   for id in keys(idtable)
      oldaid = idtable[id]
      if oldaid < 1
         newaid = 0
      elseif oldaid > n
         newaid = oldaid
      else
         newaid = invp[oldaid]
      end
      idtable[id] = newaid
   end
end

# extract the pedigree based on the permutation vector
# perm[newid]  = origid
# invp[origid] = newid
"""
    newlist,newtable = extract_ped(invp, pedlist ,idtable)
    newlist          = extract_ped(invp, pedlist)

Extract the code in `pedlist` according to the permutation list in `invp` as `invp[oldid]=newid`, which is generated with some functions like `find_ped_order`.
With optional `idtable`, it will also return with the new table.
The function re-assigns the code to animals in `newlist`.
The behavior is similar to `permute_ped`, but this function supports reducing the amount of pedigree.
"""
function extract_ped(invp::Vector{Ti},pedlist::Matrix{Ti},idtable::Dict{String,Ti}=Dict{String,Ti}()) where Ti<:Integer
   # animals in the full pedigree
   n::Ti = length(invp)
   # animals in the subset
   m::Ti = count((invp.>=1) .& (invp.<=n))
   if length(invp)!=size(pedlist,2)
      error("dimension mismatch between pedlist and invp")
   end

   # create a new subset of pedigree
   newlist = zeros(Ti,2,m)
   for oldaid=1:n
      # animal
      newaid = invp[oldaid]
      if newaid>=1 && newaid<=m
         # sire
         oldsid = pedlist[1,oldaid]
         newsid = get_new_parent_code(oldsid,n,pedlist,invp)
         # dam
         olddid = pedlist[2,oldaid]
         newdid = get_new_parent_code(olddid,n,pedlist,invp)
         # adjustment needed for unknown parent groups
         # convert oldupg>n to newupg>m by newid-n+m
         newsid = newsid>n ? newsid-n+m : newsid
         newdid = newdid>n ? newdid-n+m : newdid
         # update
         newlist[1:2,newaid] = [newsid, newdid]
      end
   end

   # change the ID table
   if length(idtable)>0
      newtable::Dict{String,Ti} = Dict{String,Ti}()
      for id in keys(idtable)
         # conversion with unknown parent groups
         oldaid = idtable[id]
         if oldaid < 1
            newaid = 0
         elseif oldaid > n
            newaid = oldaid-n+m
         else
            newaid = invp[oldaid]
         end
         newtable[id] = newaid
      end
      return newlist,newtable
   else
      return newlist
   end
end

"""
    invs = invsubidx(subidx [,n])

Returns the "inverse permutation" of sun-indices (`subidx[newid]=oldid`) resulting in `invs[oldid]=newid`.
The `subid` is similar to the permutation vector `perm`, but `subid` doesn't have to include all the index; `length(subid)` can be smaller than the maximum ID.
Additional argument `n` specifies the maximum ID in the pedigree.
"""
function invsubidx(subidx::Vector{T},n::Integer) where{T}
   invs = zeros(T,n)
   for i=1:length(subidx)
      if subidx[i]>0
         invs[subidx[i]] = i
      end
   end
   return invs
end

"""
    r = get_upg_range(pedlist)

    Finds the range of IDs for unknown parent groups.
"""
function get_upg_range(pedlist)
   return (size(pedlist,2)+1):maximum(pedlist)
end

"""
    perm,invp,subpedlist = subset_ped(pedlist,idlist)

Creates a subset of pedigree list `pedlist` with two permutation vectors: `perm[newid]=oldid` and invp[oldid]=newid`.
The subset includes the animals listed in `idlist` and their ancestors.
The subset pedigree is chronologically sorted.

This function is a wrapper for `find_ped_order` and `extract_ped`.
See the functions fro details.
"""
function subset_ped(pedlist::Matrix{Ti},idlist::Vector{Ti}=zeros(Ti,0)) where Ti<:Integer
   # extract pedigree
   # perm[newid]  = origid
   # invp[origid] = newid
   perm,invp = find_ped_order(pedlist,idlist)
   subpedlist = extract_ped(invp,pedlist)
   # new pedigree
   #n = Int(size(subpedlist,2))
   #nupg = Int(maximum(subpedlist)) - n
   return perm,invp,subpedlist
end

"""
    mgslist = mgs_ped(pedlist)
    mgslist = mgs_ped(pedlist, males)

Create a pedigree list with sire and maternal-grandsire (MGS) only.
The returned list `mgslist` has the same dimentions as `pedlist`, but the 2nd row has the MGS code.
The UPG code in `pedlist` will be removed.
The Boolean vector `males` indicates male (`true`) or female (`false`); by default, all animals will be treated as male.
"""
function mgs_ped(pedlist, males=fill(true,size(pedlist,2)))
   mgslist = similar(pedlist)
   mgslist .= 0
   n = size(pedlist,2)
   for i=1:n
      if males[i]
         s = pedlist[1,i]
         d = pedlist[2,i]
         mgslist[1,i] = s
         if d>=1 && d<=n
            mgslist[2,i] = pedlist[1,d]
         end
      end
   end
   return mgslist
end


text:
const SpinHalfSite = Union{TagType"S=1/2", TagType"SpinHalf"}

function siteinds(::SpinHalfSite,
                  N::Int; kwargs...)
  conserve_qns = get(kwargs,:conserve_qns,false)
  conserve_sz = get(kwargs,:conserve_sz,conserve_qns)
  if conserve_sz
    return [Index(QN("Sz",+1)=>1,QN("Sz",-1)=>1;tags="Site,S=1/2,n=$n") for n=1:N]
  end
  return [Index(2,"Site,S=1/2,n=$n") for n=1:N]
end

function state(::SpinHalfSite,
               st::AbstractString)
  if st == "Up" || st == "↑"
    return 1
  elseif st == "Dn" || st == "↓"
    return 2
  end
  throw(ArgumentError("State string \"$st\" not recognized for SpinHalf site"))
  return 0
end

function op(::SpinHalfSite,
            s::Index,
            opname::AbstractString;kwargs...)::ITensor
  Up = s(1)
  UpP = s'(1)
  Dn = s(2)
  DnP = s'(2)

  Op = emptyITensor(s',dag(s))

  if opname == "S⁺" || opname == "Splus" || opname == "S+"
    Op[UpP, Dn] = 1.
  elseif opname == "S⁻" || opname == "Sminus" || opname == "S-"
    Op[DnP, Up] = 1.
  elseif opname == "Sˣ" || opname == "Sx"
    Op[UpP, Dn] = 0.5
    Op[DnP, Up] = 0.5
  elseif opname == "iSʸ" || opname == "iSy"
     Op[UpP, Dn] = 0.5
     Op[DnP, Up] = -0.5
  elseif opname == "Sʸ" || opname == "Sy"
     Op = complex(Op)
     Op[UpP, Dn] = -0.5*im
     Op[DnP, Up] = 0.5*im
  elseif opname == "Sᶻ" || opname == "Sz"
     Op[UpP, Up] = 0.5
     Op[DnP, Dn] = -0.5
  elseif opname == "projUp"
     Op[UpP, Up] = 1.
  elseif opname == "projDn"
    Op[DnP, Dn] = 1.
  elseif opname == "Up" || opname == "↑"
    pU = emptyITensor(s)
    pU[Up] = 1.
    return pU
  elseif opname == "Dn" || opname == "↓"
    pD = emptyITensor(s)
    pD[Dn] = 1.
    return pD
  else
    throw(ArgumentError("Operator name '$opname' not recognized for SpinHalfSite"))
  end
  return Op
end


text:
# rawview
AAFixed{T<:FixedPoint,N} = AbstractArray{T,N}
function Base.showarg(io::IO, A::MappedArray{T,N,AA,typeof(reinterpret)}, toplevel=false) where {T<:Integer,N,AA<:AAFixed}
    print(io, "rawview(")
    Base.showarg(io, parent(A), false)
    print(io, ')')
    toplevel && print(io, " with eltype ", T)
end

# normedview
AAInteger{T<:Integer,N} = AbstractArray{T,N}
function Base.showarg(io::IO, A::MappedArray{T,N,AA,F,typeof(reinterpret)}, toplevel=false) where {T<:FixedPoint,N,AA<:AAInteger,F}
    print(io, "normedview(")
    ColorTypes.showcoloranttype(io, T)
    print(io, ", ")
    Base.showarg(io, parent(A), false)
    print(io, ')')
    toplevel && print(io, " with eltype ", T)
end

function Base.showarg(io::IO, r::Base.ReinterpretArray{T}, toplevel) where {T<:Colorant}
    print(io, "reinterpret(")
    ColorTypes.colorant_string_with_eltype(io, T)
    print(io, ", ")
    Base.showarg(io, parent(r), false)
    print(io, ')')
end

function Base.showarg(io::IO, r::Base.ReinterpretArray{T}, toplevel) where {T<:FixedPoint}
    print(io, "reinterpret(")
    ColorTypes.showcoloranttype(io, T)
    print(io, ", ")
    Base.showarg(io, parent(r), false)
    print(io, ')')
end