text:
// Check if a portal can be spawned in any direction at the location of ACV_Main.
// All area effect clouds with tag ACV_Direction that would result in an invalid portal are killed.
// Of the remaining directions exactly one is retained with the following priority:
// * ACV_Up
// * ACV_Down
// * ACV_South
// * ACV_East
// * ACV_North
// * ACV_West
process ACV_validateDirections {
ACV_spaceValidation()
ACV_collisionValidation()
ACV_surfaceValidation()

/execute @e[type=area_effect_cloud,name=ACV_Up] ~ ~ ~ kill @e[type=area_effect_cloud,name=!ACV_Up,tag=ACV_Direction,dy=0]
/execute @e[type=area_effect_cloud,name=ACV_Down] ~ ~ ~ kill @e[type=area_effect_cloud,name=!ACV_Down,tag=ACV_Direction,dy=0]
/execute @e[type=area_effect_cloud,name=ACV_South] ~ ~ ~ kill @e[type=area_effect_cloud,name=!ACV_South,tag=ACV_Direction,dy=0]
/execute @e[type=area_effect_cloud,name=ACV_East] ~ ~ ~ kill @e[type=area_effect_cloud,name=!ACV_East,tag=ACV_Direction,dy=0]
/execute @e[type=area_effect_cloud,name=ACV_North] ~ ~ ~ kill @e[type=area_effect_cloud,name=!ACV_North,tag=ACV_Direction,dy=0]
/execute @e[type=area_effect_cloud,name=ACV_West] ~ ~ ~ kill @e[type=area_effect_cloud,name=!ACV_West,tag=ACV_Direction,dy=0]
}

install {
  /scoreboard objectives add ACV_HEIGHT_VALID dummy
}

uninstall {
  /scoreboard objectives remove ACV_HEIGHT_VALID
}

// Checks if there is enough space to create the portal (is there air or carpet below a non floor Portal),
// if there is not, all directions of the affected creation are killed except for ACV_Up.
process ACV_spaceValidation {
/scoreboard players set @e[type=area_effect_cloud,tag=ACV_Direction] ACV_HEIGHT_VALID 0
/scoreboard players set @e[type=area_effect_cloud,name=ACV_Up] ACV_HEIGHT_VALID 1
/execute @e[type=area_effect_cloud,name=!ACV_Up,tag=ACV_Direction] ~ ~ ~ detect ~ ~-1 ~ air * scoreboard players set @e[type=area_effect_cloud,tag=ACV_Direction,r=0,c=1] ACV_HEIGHT_VALID 1
/execute @e[type=area_effect_cloud,name=!ACV_Up,tag=ACV_Direction] ~ ~ ~ detect ~ ~-1 ~ carpet * scoreboard players set @e[type=area_effect_cloud,tag=ACV_Direction,r=0,c=1] ACV_HEIGHT_VALID 1
/kill @e[type=area_effect_cloud,tag=ACV_Direction,score_ACV_HEIGHT_VALID=0]
}

// Checks if there is already a portal of a different color or an item frame or painting that would collide,
// if there is, all directions of the affected creation are killed.
process ACV_collisionValidation {
/execute @e[type=area_effect_cloud,name=ACV_Main] ~ ~ ~ execute @e[type=painting,dy=-1,c=1] ~ ~ ~ kill @e[type=area_effect_cloud,tag=ACV_Direction,dy=1]
/execute @e[type=area_effect_cloud,name=ACV_Main] ~ ~ ~ execute @e[type=item_frame,name=!ACV_PortalDisplay,dy=-1,c=1] ~ ~ ~ kill @e[type=area_effect_cloud,tag=ACV_Direction,dy=1]
/execute @e[type=area_effect_cloud,name=ACV_Main,score_ACV_COLOR_min=0,score_ACV_COLOR=0] ~ ~ ~ execute @e[type=item_frame,name=ACV_PortalDisplay,score_ACV_COLOR=-1,dy=-1,c=1] ~ ~ ~ kill @e[type=area_effect_cloud,tag=ACV_Direction,dy=1]
/execute @e[type=area_effect_cloud,name=ACV_Main,score_ACV_COLOR_min=0,score_ACV_COLOR=0] ~ ~ ~ execute @e[type=item_frame,name=ACV_PortalDisplay,score_ACV_COLOR_min=1,dy=-1,c=1] ~ ~ ~ kill @e[type=area_effect_cloud,tag=ACV_Direction,dy=1]
/execute @e[type=area_effect_cloud,name=ACV_Main,score_ACV_COLOR_min=1,score_ACV_COLOR=1] ~ ~ ~ execute @e[type=item_frame,name=ACV_PortalDisplay,score_ACV_COLOR=0,dy=-1,c=1] ~ ~ ~ kill @e[type=area_effect_cloud,tag=ACV_Direction,dy=1]
/execute @e[type=area_effect_cloud,name=ACV_Main,score_ACV_COLOR_min=1,score_ACV_COLOR=1] ~ ~ ~ execute @e[type=item_frame,name=ACV_PortalDisplay,score_ACV_COLOR_min=2,dy=-1,c=1] ~ ~ ~ kill @e[type=area_effect_cloud,tag=ACV_Direction,dy=1]
/execute @e[type=area_effect_cloud,name=ACV_Main,score_ACV_COLOR_min=2,score_ACV_COLOR=2] ~ ~ ~ execute @e[type=item_frame,name=ACV_PortalDisplay,score_ACV_COLOR=1,dy=-1,c=1] ~ ~ ~ kill @e[type=area_effect_cloud,tag=ACV_Direction,dy=1]
/execute @e[type=area_effect_cloud,name=ACV_Main,score_ACV_COLOR_min=2,score_ACV_COLOR=2] ~ ~ ~ execute @e[type=item_frame,name=ACV_PortalDisplay,score_ACV_COLOR_min=3,dy=-1,c=1] ~ ~ ~ kill @e[type=area_effect_cloud,tag=ACV_Direction,dy=1]
/execute @e[type=area_effect_cloud,name=ACV_Main,score_ACV_COLOR_min=3,score_ACV_COLOR=3] ~ ~ ~ execute @e[type=item_frame,name=ACV_PortalDisplay,score_ACV_COLOR=2,dy=-1,c=1] ~ ~ ~ kill @e[type=area_effect_cloud,tag=ACV_Direction,dy=1]
/execute @e[type=area_effect_cloud,name=ACV_Main,score_ACV_COLOR_min=3,score_ACV_COLOR=3] ~ ~ ~ execute @e[type=item_frame,name=ACV_PortalDisplay,score_ACV_COLOR_min=4,dy=-1,c=1] ~ ~ ~ kill @e[type=area_effect_cloud,tag=ACV_Direction,dy=1]
}

install {
  // The scoreboard Surface Valid Stat is filled by a stat command
  /scoreboard objectives add ACV_SURF_VALID_S dummy
}

uninstall {
  /scoreboard objectives remove ACV_SURF_VALID_S
}

// Checks if both the upper and lower block is valid.
// They are if they are non air, and either listed by the Block_Option or the Block_Option is empty.
// Invalid directions are killed.
process ACV_surfaceValidation {
/execute @e[type=area_effect_cloud,name=ACV_Up] ~ ~ ~ summon area_effect_cloud ~ ~-1 ~ {CustomName:"ACV_LowerBlock",Tags:[ACV_Up],Duration:2147483647}
/execute @e[type=area_effect_cloud,name=ACV_Up] ~ ~ ~ summon area_effect_cloud ~ ~-1 ~ {CustomName:"ACV_UpperBlock",Tags:[ACV_Up],Duration:2147483647}

/execute @e[type=area_effect_cloud,name=ACV_Down] ~ ~ ~ summon area_effect_cloud ~ ~1 ~ {CustomName:"ACV_LowerBlock",Tags:[ACV_Down],Duration:2147483647}
/execute @e[type=area_effect_cloud,name=ACV_Down] ~ ~ ~ summon area_effect_cloud ~ ~1 ~ {CustomName:"ACV_UpperBlock",Tags:[ACV_Down],Duration:2147483647}

/execute @e[type=area_effect_cloud,name=ACV_South] ~ ~ ~ summon area_effect_cloud ~ ~-1 ~-1 {CustomName:"ACV_LowerBlock",Tags:[ACV_South],Duration:2147483647}
/execute @e[type=area_effect_cloud,name=ACV_South] ~ ~ ~ summon area_effect_cloud ~ ~ ~-1 {CustomName:"ACV_UpperBlock",Tags:[ACV_South],Duration:2147483647}

/execute @e[type=area_effect_cloud,name=ACV_East] ~ ~ ~ summon area_effect_cloud ~-1 ~-1 ~ {CustomName:"ACV_LowerBlock",Tags:[ACV_East],Duration:2147483647}
/execute @e[type=area_effect_cloud,name=ACV_East] ~ ~ ~ summon area_effect_cloud ~-1 ~ ~ {CustomName:"ACV_UpperBlock",Tags:[ACV_East],Duration:2147483647}

/execute @e[type=area_effect_cloud,name=ACV_North] ~ ~ ~ summon area_effect_cloud ~ ~-1 ~1 {CustomName:"ACV_LowerBlock",Tags:[ACV_North],Duration:2147483647}
/execute @e[type=area_effect_cloud,name=ACV_North] ~ ~ ~ summon area_effect_cloud ~ ~ ~1 {CustomName:"ACV_UpperBlock",Tags:[ACV_North],Duration:2147483647}

/execute @e[type=area_effect_cloud,name=ACV_West] ~ ~ ~ summon area_effect_cloud ~1 ~-1 ~ {CustomName:"ACV_LowerBlock",Tags:[ACV_West],Duration:2147483647}
/execute @e[type=area_effect_cloud,name=ACV_West] ~ ~ ~ summon area_effect_cloud ~1 ~ ~ {CustomName:"ACV_UpperBlock",Tags:[ACV_West],Duration:2147483647}

/scoreboard players add @e[type=area_effect_cloud,name=ACV_LowerBlock] ACV_SURF_VALID_S 0
/scoreboard players add @e[type=area_effect_cloud,name=ACV_UpperBlock] ACV_SURF_VALID_S 0
/stats entity @e[type=area_effect_cloud,name=ACV_LowerBlock] set SuccessCount @e[type=area_effect_cloud,name=ACV_LowerBlock,r=0,c=1] ACV_SURF_VALID_S
/stats entity @e[type=area_effect_cloud,name=ACV_UpperBlock] set SuccessCount @e[type=area_effect_cloud,name=ACV_UpperBlock,r=0,c=1] ACV_SURF_VALID_S

ACV_loadChunks()

/testforblocks 0 1 1 6 2 1 0 1 2
conditional: /scoreboard players tag @e[type=area_effect_cloud,name=ACV_LowerBlock] add ACV_SurfaceValid
conditional: /scoreboard players tag @e[type=area_effect_cloud,name=ACV_UpperBlock] add ACV_SurfaceValid

/execute @e[type=area_effect_cloud,name=ACV_LowerBlock] ~ ~ ~ testforblocks ~ ~ ~ ~ ~ ~ 0 1 1
/execute @e[type=area_effect_cloud,name=ACV_LowerBlock,score_ACV_SURF_VALID_S_min=1] ~ ~ ~ scoreboard players tag @e[type=area_effect_cloud,name=ACV_LowerBlock,r=0,c=1] add ACV_SurfaceValid
/execute @e[type=area_effect_cloud,name=ACV_LowerBlock] ~ ~ ~ testforblocks ~ ~ ~ ~ ~ ~ 1 1 1
/execute @e[type=area_effect_cloud,name=ACV_LowerBlock,score_ACV_SURF_VALID_S_min=1] ~ ~ ~ scoreboard players tag @e[type=area_effect_cloud,name=ACV_LowerBlock,r=0,c=1] add ACV_SurfaceValid
/execute @e[type=area_effect_cloud,name=ACV_LowerBlock] ~ ~ ~ testforblocks ~ ~ ~ ~ ~ ~ 2 1 1
/execute @e[type=area_effect_cloud,name=ACV_LowerBlock,score_ACV_SURF_VALID_S_min=1] ~ ~ ~ scoreboard players tag @e[type=area_effect_cloud,name=ACV_LowerBlock,r=0,c=1] add ACV_SurfaceValid
/execute @e[type=area_effect_cloud,name=ACV_LowerBlock] ~ ~ ~ testforblocks ~ ~ ~ ~ ~ ~ 3 1 1
/execute @e[type=area_effect_cloud,name=ACV_LowerBlock,score_ACV_SURF_VALID_S_min=1] ~ ~ ~ scoreboard players tag @e[type=area_effect_cloud,name=ACV_LowerBlock,r=0,c=1] add ACV_SurfaceValid
/execute @e[type=area_effect_cloud,name=ACV_LowerBlock] ~ ~ ~ testforblocks ~ ~ ~ ~ ~ ~ 4 1 1
/execute @e[type=area_effect_cloud,name=ACV_LowerBlock,score_ACV_SURF_VALID_S_min=1] ~ ~ ~ scoreboard players tag @e[type=area_effect_cloud,name=ACV_LowerBlock,r=0,c=1] add ACV_SurfaceValid
/execute @e[type=area_effect_cloud,name=ACV_LowerBlock] ~ ~ ~ testforblocks ~ ~ ~ ~ ~ ~ 5 1 1
/execute @e[type=area_effect_cloud,name=ACV_LowerBlock,score_ACV_SURF_VALID_S_min=1] ~ ~ ~ scoreboard players tag @e[type=area_effect_cloud,name=ACV_LowerBlock,r=0,c=1] add ACV_SurfaceValid
/execute @e[type=area_effect_cloud,name=ACV_LowerBlock] ~ ~ ~ testforblocks ~ ~ ~ ~ ~ ~ 6 1 1
/execute @e[type=area_effect_cloud,name=ACV_LowerBlock,score_ACV_SURF_VALID_S_min=1] ~ ~ ~ scoreboard players tag @e[type=area_effect_cloud,name=ACV_LowerBlock,r=0,c=1] add ACV_SurfaceValid
/execute @e[type=area_effect_cloud,name=ACV_LowerBlock] ~ ~ ~ testforblocks ~ ~ ~ ~ ~ ~ 0 2 1
/execute @e[type=area_effect_cloud,name=ACV_LowerBlock,score_ACV_SURF_VALID_S_min=1] ~ ~ ~ scoreboard players tag @e[type=area_effect_cloud,name=ACV_LowerBlock,r=0,c=1] add ACV_SurfaceValid
/execute @e[type=area_effect_cloud,name=ACV_LowerBlock] ~ ~ ~ testforblocks ~ ~ ~ ~ ~ ~ 1 2 1
/execute @e[type=area_effect_cloud,name=ACV_LowerBlock,score_ACV_SURF_VALID_S_min=1] ~ ~ ~ scoreboard players tag @e[type=area_effect_cloud,name=ACV_LowerBlock,r=0,c=1] add ACV_SurfaceValid
/execute @e[type=area_effect_cloud,name=ACV_LowerBlock] ~ ~ ~ testforblocks ~ ~ ~ ~ ~ ~ 2 2 1
/execute @e[type=area_effect_cloud,name=ACV_LowerBlock,score_ACV_SURF_VALID_S_min=1] ~ ~ ~ scoreboard players tag @e[type=area_effect_cloud,name=ACV_LowerBlock,r=0,c=1] add ACV_SurfaceValid
/execute @e[type=area_effect_cloud,name=ACV_LowerBlock] ~ ~ ~ testforblocks ~ ~ ~ ~ ~ ~ 3 2 1
/execute @e[type=area_effect_cloud,name=ACV_LowerBlock,score_ACV_SURF_VALID_S_min=1] ~ ~ ~ scoreboard players tag @e[type=area_effect_cloud,name=ACV_LowerBlock,r=0,c=1] add ACV_SurfaceValid
/execute @e[type=area_effect_cloud,name=ACV_LowerBlock] ~ ~ ~ testforblocks ~ ~ ~ ~ ~ ~ 4 2 1
/execute @e[type=area_effect_cloud,name=ACV_LowerBlock,score_ACV_SURF_VALID_S_min=1] ~ ~ ~ scoreboard players tag @e[type=area_effect_cloud,name=ACV_LowerBlock,r=0,c=1] add ACV_SurfaceValid
/execute @e[type=area_effect_cloud,name=ACV_LowerBlock] ~ ~ ~ testforblocks ~ ~ ~ ~ ~ ~ 5 2 1
/execute @e[type=area_effect_cloud,name=ACV_LowerBlock,score_ACV_SURF_VALID_S_min=1] ~ ~ ~ scoreboard players tag @e[type=area_effect_cloud,name=ACV_LowerBlock,r=0,c=1] add ACV_SurfaceValid
/execute @e[type=area_effect_cloud,name=ACV_LowerBlock] ~ ~ ~ testforblocks ~ ~ ~ ~ ~ ~ 6 2 1
/execute @e[type=area_effect_cloud,name=ACV_LowerBlock,score_ACV_SURF_VALID_S_min=1] ~ ~ ~ scoreboard players tag @e[type=area_effect_cloud,name=ACV_LowerBlock,r=0,c=1] add ACV_SurfaceValid

/execute @e[type=area_effect_cloud,name=ACV_UpperBlock] ~ ~ ~ testforblocks ~ ~ ~ ~ ~ ~ 0 1 1
/execute @e[type=area_effect_cloud,name=ACV_UpperBlock,score_ACV_SURF_VALID_S_min=1] ~ ~ ~ scoreboard players tag @e[type=area_effect_cloud,name=ACV_UpperBlock,r=0,c=1] add ACV_SurfaceValid
/execute @e[type=area_effect_cloud,name=ACV_UpperBlock] ~ ~ ~ testforblocks ~ ~ ~ ~ ~ ~ 1 1 1
/execute @e[type=area_effect_cloud,name=ACV_UpperBlock,score_ACV_SURF_VALID_S_min=1] ~ ~ ~ scoreboard players tag @e[type=area_effect_cloud,name=ACV_UpperBlock,r=0,c=1] add ACV_SurfaceValid
/execute @e[type=area_effect_cloud,name=ACV_UpperBlock] ~ ~ ~ testforblocks ~ ~ ~ ~ ~ ~ 2 1 1
/execute @e[type=area_effect_cloud,name=ACV_UpperBlock,score_ACV_SURF_VALID_S_min=1] ~ ~ ~ scoreboard players tag @e[type=area_effect_cloud,name=ACV_UpperBlock,r=0,c=1] add ACV_SurfaceValid
/execute @e[type=area_effect_cloud,name=ACV_UpperBlock] ~ ~ ~ testforblocks ~ ~ ~ ~ ~ ~ 3 1 1
/execute @e[type=area_effect_cloud,name=ACV_UpperBlock,score_ACV_SURF_VALID_S_min=1] ~ ~ ~ scoreboard players tag @e[type=area_effect_cloud,name=ACV_UpperBlock,r=0,c=1] add ACV_SurfaceValid
/execute @e[type=area_effect_cloud,name=ACV_UpperBlock] ~ ~ ~ testforblocks ~ ~ ~ ~ ~ ~ 4 1 1
/execute @e[type=area_effect_cloud,name=ACV_UpperBlock,score_ACV_SURF_VALID_S_min=1] ~ ~ ~ scoreboard players tag @e[type=area_effect_cloud,name=ACV_UpperBlock,r=0,c=1] add ACV_SurfaceValid
/execute @e[type=area_effect_cloud,name=ACV_UpperBlock] ~ ~ ~ testforblocks ~ ~ ~ ~ ~ ~ 5 1 1
/execute @e[type=area_effect_cloud,name=ACV_UpperBlock,score_ACV_SURF_VALID_S_min=1] ~ ~ ~ scoreboard players tag @e[type=area_effect_cloud,name=ACV_UpperBlock,r=0,c=1] add ACV_SurfaceValid
/execute @e[type=area_effect_cloud,name=ACV_UpperBlock] ~ ~ ~ testforblocks ~ ~ ~ ~ ~ ~ 6 1 1
/execute @e[type=area_effect_cloud,name=ACV_UpperBlock,score_ACV_SURF_VALID_S_min=1] ~ ~ ~ scoreboard players tag @e[type=area_effect_cloud,name=ACV_UpperBlock,r=0,c=1] add ACV_SurfaceValid
/execute @e[type=area_effect_cloud,name=ACV_UpperBlock] ~ ~ ~ testforblocks ~ ~ ~ ~ ~ ~ 0 2 1
/execute @e[type=area_effect_cloud,name=ACV_UpperBlock,score_ACV_SURF_VALID_S_min=1] ~ ~ ~ scoreboard players tag @e[type=area_effect_cloud,name=ACV_UpperBlock,r=0,c=1] add ACV_SurfaceValid
/execute @e[type=area_effect_cloud,name=ACV_UpperBlock] ~ ~ ~ testforblocks ~ ~ ~ ~ ~ ~ 1 2 1
/execute @e[type=area_effect_cloud,name=ACV_UpperBlock,score_ACV_SURF_VALID_S_min=1] ~ ~ ~ scoreboard players tag @e[type=area_effect_cloud,name=ACV_UpperBlock,r=0,c=1] add ACV_SurfaceValid
/execute @e[type=area_effect_cloud,name=ACV_UpperBlock] ~ ~ ~ testforblocks ~ ~ ~ ~ ~ ~ 2 2 1
/execute @e[type=area_effect_cloud,name=ACV_UpperBlock,score_ACV_SURF_VALID_S_min=1] ~ ~ ~ scoreboard players tag @e[type=area_effect_cloud,name=ACV_UpperBlock,r=0,c=1] add ACV_SurfaceValid
/execute @e[type=area_effect_cloud,name=ACV_UpperBlock] ~ ~ ~ testforblocks ~ ~ ~ ~ ~ ~ 3 2 1
/execute @e[type=area_effect_cloud,name=ACV_UpperBlock,score_ACV_SURF_VALID_S_min=1] ~ ~ ~ scoreboard players tag @e[type=area_effect_cloud,name=ACV_UpperBlock,r=0,c=1] add ACV_SurfaceValid
/execute @e[type=area_effect_cloud,name=ACV_UpperBlock] ~ ~ ~ testforblocks ~ ~ ~ ~ ~ ~ 4 2 1
/execute @e[type=area_effect_cloud,name=ACV_UpperBlock,score_ACV_SURF_VALID_S_min=1] ~ ~ ~ scoreboard players tag @e[type=area_effect_cloud,name=ACV_UpperBlock,r=0,c=1] add ACV_SurfaceValid
/execute @e[type=area_effect_cloud,name=ACV_UpperBlock] ~ ~ ~ testforblocks ~ ~ ~ ~ ~ ~ 5 2 1
/execute @e[type=area_effect_cloud,name=ACV_UpperBlock,score_ACV_SURF_VALID_S_min=1] ~ ~ ~ scoreboard players tag @e[type=area_effect_cloud,name=ACV_UpperBlock,r=0,c=1] add ACV_SurfaceValid
/execute @e[type=area_effect_cloud,name=ACV_UpperBlock] ~ ~ ~ testforblocks ~ ~ ~ ~ ~ ~ 6 2 1
/execute @e[type=area_effect_cloud,name=ACV_UpperBlock,score_ACV_SURF_VALID_S_min=1] ~ ~ ~ scoreboard players tag @e[type=area_effect_cloud,name=ACV_UpperBlock,r=0,c=1] add ACV_SurfaceValid

/execute @e[type=area_effect_cloud,name=ACV_LowerBlock] ~ ~ ~ detect ~ ~ ~ piston_extension * scoreboard players tag @e[type=area_effect_cloud,name=ACV_LowerBlock,r=0,c=1] add ACV_SurfaceValid
/execute @e[type=area_effect_cloud,name=ACV_UpperBlock] ~ ~ ~ detect ~ ~ ~ piston_extension * scoreboard players tag @e[type=area_effect_cloud,name=ACV_UpperBlock,r=0,c=1] add ACV_SurfaceValid

/execute @e[type=area_effect_cloud,name=ACV_LowerBlock] ~ ~ ~ detect ~ ~ ~ air * scoreboard players tag @e[type=area_effect_cloud,name=ACV_LowerBlock,r=0,c=1] remove ACV_SurfaceValid
/execute @e[type=area_effect_cloud,name=ACV_UpperBlock] ~ ~ ~ detect ~ ~ ~ air * scoreboard players tag @e[type=area_effect_cloud,name=ACV_UpperBlock,r=0,c=1] remove ACV_SurfaceValid

/execute @e[type=area_effect_cloud,name=ACV_LowerBlock,tag=ACV_Up,tag=!ACV_SurfaceValid] ~ ~1 ~ kill @e[type=area_effect_cloud,name=ACV_Up,dy=0]
/execute @e[type=area_effect_cloud,name=ACV_UpperBlock,tag=ACV_Up,tag=!ACV_SurfaceValid] ~ ~1 ~ kill @e[type=area_effect_cloud,name=ACV_Up,dy=0]

/execute @e[type=area_effect_cloud,name=ACV_LowerBlock,tag=ACV_Down,tag=!ACV_SurfaceValid] ~ ~-1 ~ kill @e[type=area_effect_cloud,name=ACV_Down,dy=0]
/execute @e[type=area_effect_cloud,name=ACV_UpperBlock,tag=ACV_Down,tag=!ACV_SurfaceValid] ~ ~-1 ~ kill @e[type=area_effect_cloud,name=ACV_Down,dy=0]

/execute @e[type=area_effect_cloud,name=ACV_LowerBlock,tag=ACV_South,tag=!ACV_SurfaceValid] ~ ~1 ~1 kill @e[type=area_effect_cloud,name=ACV_South,dy=0]
/execute @e[type=area_effect_cloud,name=ACV_UpperBlock,tag=ACV_South,tag=!ACV_SurfaceValid] ~ ~ ~1 kill @e[type=area_effect_cloud,name=ACV_South,dy=0]

/execute @e[type=area_effect_cloud,name=ACV_LowerBlock,tag=ACV_East,tag=!ACV_SurfaceValid] ~1 ~1 ~ kill @e[type=area_effect_cloud,name=ACV_East,dy=0]
/execute @e[type=area_effect_cloud,name=ACV_UpperBlock,tag=ACV_East,tag=!ACV_SurfaceValid] ~1 ~ ~ kill @e[type=area_effect_cloud,name=ACV_East,dy=0]

/execute @e[type=area_effect_cloud,name=ACV_LowerBlock,tag=ACV_North,tag=!ACV_SurfaceValid] ~ ~1 ~-1 kill @e[type=area_effect_cloud,name=ACV_North,dy=0]
/execute @e[type=area_effect_cloud,name=ACV_UpperBlock,tag=ACV_North,tag=!ACV_SurfaceValid] ~ ~ ~-1 kill @e[type=area_effect_cloud,name=ACV_North,dy=0]

/execute @e[type=area_effect_cloud,name=ACV_LowerBlock,tag=ACV_West,tag=!ACV_SurfaceValid] ~-1 ~1 ~ kill @e[type=area_effect_cloud,name=ACV_West,dy=0]
/execute @e[type=area_effect_cloud,name=ACV_UpperBlock,tag=ACV_West,tag=!ACV_SurfaceValid] ~-1 ~ ~ kill @e[type=area_effect_cloud,name=ACV_West,dy=0]

/kill @e[type=area_effect_cloud,name=ACV_LowerBlock]
/kill @e[type=area_effect_cloud,name=ACV_UpperBlock]
}


text:
######################################################################

`eta/nonempty_subsets` := (A::set) -> `if`(nops(A)=1,A,FAIL);

`gamma/nonempty_subsets` := (A::set,B::set) -> (p) -> proc(U,V)
 map(op,map(u -> V[u],U));
end;

ext: mpl


text:
`is_element/BW1` := eval(`is_element/prime_simplex`);
`is_equal/BW1` := eval(`is_equal/prime_simplex`);
`is_leq/BW1` := NULL;
`random_element/BW1` := eval(`random_element/prime_simplex`);
`list_elements/BW1` := NULL;
`count_elements/BW1` := NULL;
`phi/nonempty_subsets/BW1` := eval(`phi/nonempty_subsets/prime_simplex`);


text:
interface(screenwidth=9999):
if not (NewSLO :: `module`) then
  WARNING("loading NewSLO failed");
  `quit`(3);
end if;

with(Hakaru):
with(NewSLO):

#####################################################################
#
# Dirichlet conjugacy tests
#
#####################################################################

(dice_index, dice_index_t) := Concrete("examples/dice_index.hk"):
dice_index := value(eval(dice_index)):

# BUG: naming the var `lda' with a lowercase causes "Error: recursive
# assignment", but *not* inside Concrete (in the evaluation of the resulting
# expression!). `lda' doesn't seem to be a special name.
(LDA, LDA_t) := Concrete("examples/lda2.hk"):
LDA := value(eval(LDA)):

(gmm_gibbs, gmm_gibbs_t) := Concrete("examples/gmm_gibbs.hk"):
gmm_gibbs := value(eval(gmm_gibbs)):

(naive_bayes_gibbs, naive_bayes_gibbs_t) := Concrete("examples/naive_bayes_gibbs.hk"):
naive_bayes_gibbs := value(eval(naive_bayes_gibbs)):

#####################################################################
TestEfficient( dice_index, dice_index_t, KB:-empty, label="dice index" ):
TestEfficient( LDA, LDA_t, KB:-empty, label="LDA" ):
TestEfficient( gmm_gibbs, gmm_gibbs_t, KB:-empty, label="gmm gibbs" ):
TestEfficient( naive_bayes_gibbs, naive_bayes_gibbs_t, KB:-empty, label="naive bayes gibbs" ):

max_stars_repo_name: vmchale/hakaru

max_issues_repo_name: vmchale/hakaru

max_forks_repo_name: vmchale/hakaru


text:
# A maple code for the dual Leipniz-Hopf(LH) algebra
# Description: This code includes functions to compute
# the orverlapping product and the conjugation in the dual
# Leipniz-Hopf algebra
# Author: S. Kaji (Yamaguchi University, Japan)
# Version history: initial version -- 2 Feb. 2014.
#                  current version -- 2 Jan. 2015.
# Global variables:
# T[]: basis for LH-algebra
# S[]: basis for dual LH-algebra
# P[]: Cartan-Serre basis for Steenrod algebra
# Q[]: Cartan-Serre basis for dual Steenrod algebra
# x[]: Milnor basis for dual Steenrod algebra

with(LinearAlgebra[Modular]);

## the dual Leibniz-Hopf algebra
# overlapping shuffle product
os := proc (f::polynom) local ff, L, t; global p;
    ff := expand(f) mod p;
    if degree(ff) < 2 then
        return ff
    elif op(0, ff) = `+` then
        return expand(map(os, ff)) mod p
    elif op(0, ff) = `*` then
        if type(op(1, ff), 'numeric') then
            return op(1, ff)*os(ff/op(1, ff)) mod p;
        elif nops(ff) = 2 and type(op(1, ff), 'indexed') and type(op(2, ff), 'indexed') then
            L := osw(op(1, op(1, ff)), op(1, op(2, ff)));
            return `mod`(add(S[l], l in L), p)
        else
            return `mod`(expand(os(os(op(1, ff))*os(mul(op(i, ff), i = 2 .. nops(ff))))), p)
        end if;
    elif op(0, ff) = `^` and type(op(1, ff), 'indexed') then
        t := floor(log[p](op(2, ff)));
        if 0 < t then
            return `mod`(os(S[p^t*op(1, op(1, ff))]*op(1, ff)^(op(2, ff)-p^t)), p)
        else L := osw(op(1, op(1, ff)), op(1, op(1, ff)));
            return os(add(S[l], l in L)*op(1, ff)^(op(2, ff)-2))
        end if;
    else
        print(ff, op(0, ff), "Error: no rule to compute the product")
    end if;
end proc;

#overlapping shuffle of two words
osw := proc (v::list, w::list) local L, i; option remember;
    if nops(w) = 0 then
        return [v]
    end if;
    if nops(v) = 0 then
        return [w]
    end if;
    L := [seq([w[1], op(x)], x in osw(v, w[2 .. -1]))];
    for i to nops(v) do
        L := [op(L), seq([op(v[1 .. i]), w[1], op(x)], x in osw(v[i+1 .. -1], w[2 .. -1]))];
        L := [op(L), seq([op(v[1 .. i-1]), v[i]+w[1], op(x)], x in osw(v[i+1 .. -1], w[2 .. -1]))]
    end do;
    return L;
end proc;

# pi image of milnor's dual basis
xi := proc (n) local i, x; global p;
    if type(n, 'integer') then
        if n = 0 then
            return 1
        else
            return S[[seq(p^(n-i), i = 1 .. n)]]
        end if
    elif type(n, 'list') then
        x := 1;
        for i to nops(n) do
            x := x*xi(i)^n[i]
        end do;
        return x
    else print("Error: xi")
    end if
end proc;

# convert a polynomial in x to one in S
x2xi := proc(f::polynom) local s,ind;
    ind := [];
    for s in indets(f) do
        if op(0,s)=`x` and type(s,'indexed') then
            ind := [op(ind),op(1,s)]
        end if;
    end do;
    return eval(f,[seq(x[i]=xi(i),i=ind)]);
end proc:

# conjugation for dual LH-algebra
chi := proc (f::polynom) global p;
    if degree(f) < 1 then
        return f
    elif type(f, 'indexed') then
        return add((-1)^nops(op(1, f))*S[x], x in coarsening(ListTools[Reverse](op(1, f))))
    elif op(0, f) = `+` then
        return `mod`(map(chi, f), p)
    elif op(0, f) = `*` then
        if type(op(1, f), 'numeric') then
            return `mod`(op(1, f)*chi(f/op(1, f)), p)
        else
            return `mod`(chi(op(1, f))*chi(mul(op(i, f), i = 2 .. nops(f))), p)
        end if
    elif op(0, f) = `^` then
        return `mod`(chi(op(1, f))^op(2, f), p)
    else
        print("cannot compute chi:", f, op(0, f))
    end if
end proc;

# conjugation for LH-algebra
chi2 := proc (f::polynom) global p;
    if degree(f) < 1 then
        return f
    elif type(f, 'indexed') then
        return add((-1)^nops(op(1, f))*T[x], x in refinement(ListTools[Reverse](op(1, f))))
    elif op(0, f) = `+` then
        return `mod`(map(chi2, f), p)
    elif op(0, f) = `*` then
        if type(op(1, f), 'numeric') then
            return `mod`(op(1, f)*chi2(f/op(1, f)), p)
        else
            return `mod`(chi2(op(1, f))*chi2(mul(op(i, f), i = 2 .. nops(f))), p)
        end if
    elif op(0, f) = `^` then
        return `mod`(chi2(op(1, f))^op(2, f), p)
    else
        print("cannot compute chi2:", f, op(0, f))
    end if
end proc;

# retraction onto the dual steenrod algebra (taking admissble terms)
retract := proc (f::polynom) local ff, i, w; global p;
    ff := expand(f);
    if 1 < degree(ff) then
        print(ff, "input must be linear")
    elif op(0, ff) = `+` then
        return map(retract, ff)
    elif type(ff, 'indexed') then
        if op(0,ff) != `S` then return ff; end if;
        w := op(1, ff);
        for i to nops(w)-1 do
            if w[i] < p*w[i+1] then
                return 0
            end if
        end do;
        return Q[w]
    elif op(0, ff) = `*` and type(op(1, ff), 'numeric') then
        return `mod`(op(1, ff)*retract(ff/op(1, ff)), p)
    else
        print("cannot compute retract:",ff, op(0, ff));
        return 0
    end if
end proc;

# projection onto the steenrod algebra (applying Adem relations)
projection := proc (f::polynom) local ff, i, w, t, j; global p;
    ff := expand(f) mod p;
    if 1 < degree(ff) then
        print(ff, "input must be linear")
    elif type(ff,'numeric') then
        return ff;
    elif op(0, ff) = `+` then
        return map(procname, ff) mod p
    elif type(ff, 'indexed') then
        if op(0,ff) != `T` and op(0,ff) != `P` then return ff; end if;
        w := op(1, ff);
        for i to nops(w)-1 do
            if w[i] < p*w[i+1] then
                t := ((-1)^(w[i]))*binomial((p-1)*w[i+1]-1,w[i])
                      *P[[op(w[1..i-1]),w[i]+w[i+1],op(w[i+2..-1])]];
                for j to floor(w[i]/p) do
                    t:=t + ((-1)^(w[i]+j))*binomial((p-1)*(w[i+1]-j)-1,w[i]-p*j)
                      *P[[op(w[1..i-1]),w[i]+w[i+1]-j,j,op(w[i+2..-1])]];
                end do;
                return procname(t mod p);
            end if
        end do;
        return P[w]
    elif op(0, ff) = `*` and type(op(1, ff), 'numeric') then
        return `mod`(op(1, ff)*procname(ff/op(1, ff)), p)
    else
        print("cannot compute projection:",ff, op(0, ff));
        return 0
    end if
end proc;

## Combinatorics on words
# ordered partition of an integer
ordpart := proc (n::integer) local L, V; option remember;
    if n = 0 then
        return []
    elif n = 1 then
        return [[1]]
    end if;
    L := ordpart(n-1);
    V[1] := [seq([1, op(v)], v in L)];
    V[2] := [seq([v[1]+1, op(v[2 .. -1])], v in L)];
    return [op(V[1]), op(V[2])]
end proc;
# ordered partition of a list
ordpartlist := proc (L::list) local V, u, i, w, pos;
    V := [];
    for u in ordpart(nops(L)) do
        w := [];
        pos := 1;
        for i in u do
            w := [op(w), L[pos .. pos+i-1]];
            pos := pos+i
        end do;
        V := [op(V), w]
    end do;
    return V
end proc;

# coarsening of a word
coarsening := proc (w::list) local L, p, i, l, v; option remember;
    L := [];
    for p in ordpart(nops(w)) do
        v := [];
        l := 1;
        for i to nops(p) do
            v := [op(v), add(w[j], j = l .. l+p[i]-1)];
            l := l+p[i]
        end do;
        L := [op(L), v]
    end do;
    return L
end proc;

# refinement of a word
refinement := proc (w::list) local L, l, v, U; option remember;
    if nops(w)=1 then
        return ordpart(w[1]);
    else
        U := [];
        L := refinement(w[2..-1]);
        for v in ordpart(w[1]) do
            for l in L do
                U := [op(U), [op(v),op(l)]];
            end do;
        end do;
        return U;
    end if;
end proc:

# list monomials of degree "totaldeg" consisting of "gens",
# whose degrees are specified by "degs"
list_monom := proc (gens::list, degs::(list(integer)), totaldeg::integer) local v, L, V;
    L := wpartition(degs, totaldeg);
    V := [];
    for v in L do
        V := [op(V), mul(gens[i]^v[i], i = 1 .. nops(v))]
    end do;
    return V
end proc;

# weighted partition; express "totaldeg" as linear combinations of "degs"
wpartition := proc (degs::(list(integer)), totaldeg::integer) local i, L, V;
    L := [];
    if totaldeg = 0 then
        return [[seq(0, i = 1 .. nops(degs))]]
    elif nops(degs) = 0 then
        return []
    elif nops(degs) = 1 then
        if `mod`(totaldeg, degs[1]) = 0 then
            return [[totaldeg/degs[1]]]
        else
            return []
        end if
    end if;
    for i from floor(totaldeg/degs[-1]) by -1 to 0 do
        V := wpartition(degs[1 .. -2], totaldeg-i*degs[-1]);
        L := [op(L), seq([op(v), i], v in V)]
    end do;
    return L
end proc;

## Enumerating basis and the conjugation invariants
# enumerate admissible sequences with lex order
admseq := proc (totaldeg::integer) local n, L; global p;
    n := floor(log[p]((p-1)*totaldeg+1));
    L := wpartition([seq((p^i-1)/(p-1), i = 1 .. n)], totaldeg);
    return map(ex2adm, L)
end proc;

# convert excess representation to admissible sequence
ex2adm := proc (ex::list) local n; global p;
    n := nops(ex);
    while ex[n] = 0 do
        n := n-1
    end do;
    return [seq(add(p^(i-j)*ex[i], i = j .. nops(ex)), j = 1 .. n)]
end proc;

# enumerate pi image of milnor monomials of degree "deg"
milnor_images := proc (deg::integer) local L, V, i, s, n; global p;
    V := [];
    n := floor(log[p]((p-1)*deg+1));
    L := list_monom([seq(x[i], i = 1 .. n)], [seq((p^i-1)/(p-1), i = 1 .. n)], deg);
    for s in L do
        V := [op(V), s = os(eval(s, [seq(x[i] = xi(i), i = 1 .. n)]))]
    end do;
    return V
end proc;

# polynomial associated to a list
pp := proc (L::list) local x, w;
    x := 0;
    for w in L do
        x := x+mul(S[u], u in w)
    end do;
    return x
end proc;

# linear equation solver
# Solves equations L in terms of indeterminants A
sol := proc (L::list, A::list) local X, V, i, j, v, n, Y,Z, rp, d, R, U, r; global p;
    n := nops(L);
    X := Matrix(n, n+nops(A));
    for i to n do
        X[i, nops(A)+i] := 1;
        for j to nops(A) do
            X[i, j] := coeff(op(2, L[i]), A[j])
        end do
    end do;
    X := Mod(p, X, integer[]);
    U := Mod(p, X, integer[]);
    Y, rp, d := RowEchelonTransform(p, U, true, true, false, false);
    r := nops(rp); U := U[1 .. n, 1 .. n];
    Z := LinearAlgebra[IdentityMatrix](n, compact = false);
    for i to LinearAlgebra[Dimension](Y) do
        LinearAlgebra[RowOperation](Z, [i, Y[i]], inplace = true)
    end do;
    Z := Mod(p, Z, integer[]);
    X := Multiply(p, Z, X);
    X := Multiply(p, U, X);
    V := [];
    for v in convert(X, listlist) do
        V := [op(V), add(v[i]*A[i], i = 1 .. nops(A)) = add(v[nops(A)+i]*op(1, L[i]), i = 1 .. nops(L))]
    end do;
    return V
end proc;

# compute conjugation invariants in degree "deg" in the dual Steenrod algebra
conjinvQ := proc (deg::integer) local L, V, U, v, w; global p;
    L := milnor_images(deg);
    V := [];
    for v in L do
        V := [op(V), op(1, v) = `mod`(chi(op(2, v))-op(2, v), p)]
    end do;
    U := [];
    for v in sol(V, [seq(S[w], w=admseq(deg))]) do
        if op(1,v)=0 then
            U:=[op(U), [op(2, v),retract(os(x2xi(op(2,v))))]];
        end if;
    end do:
    return U;
end proc;
# compute conjugation invariants in degree "deg" in the dual LH algebra
conjinvS := proc (deg::integer) local L, V, U, v, w; global p;
    L := ordpart(deg);
    V := [];
    for v in L do
        V := [op(V), S[v] = `mod`(chi(S[v])-S[v], p)]
    end do;
    U:=[];
    for v in sol(V, [seq(S[w], w=L)]) do
        if op(1, v) = 0 then
            U:=[op(U),srt(op(2, v))]
        end if;
    end do;
    return U
end proc;
# compute conjugation invariants in degree "deg" in the LH algebra
conjinvT := proc (deg::integer) local L, V, U, v, w; global p;
    L := ordpart(deg);
    V := [];
    for v in L do
        V := [op(V), T[v] = `mod`(chi2(T[v])-T[v], p)]
    end do;
    U:=[];
    for v in sol(V, [seq(T[w], w=L)]) do
        if op(1, v) = 0 then
            U:=[op(U),srt(op(2, v))]
        end if;
    end do;
    return U
end proc;
# compute conjugation invariants in degree "deg" in the Steenrod algebra
conjinvP := proc (deg::integer) local L, V, U, v; global p;
    L := admseq(deg);
    V := [];
    for v in L do
        V := [op(V), P[v] = projection(chi2(T[v]) - T[v]) mod p]
    end do;
    U:=[];
    for v in sol(V, [seq(P[w],w=L)]) do
        if op(1, v) = 0 then
            U:=[op(U),srt(op(2, v))]
        end if;
    end do;
    return U
end proc;

# sort terms in reverse lex order
srt := proc (f::polynom) local L;
    if degree(f) <> 1 then
        print("degree must be one", f);
        return f
    end if;
    L := sort([seq(ListTools[Reverse](op(1, s)), s in indets(f))]);
    L := [seq(S[ListTools[Reverse](l)], l in L)];
    return
    sort(f, L)
end proc;

# duality between S_I and T^J
# Note that S_J is conjugation invariant <=> So is T^K in non-dual LH-algebra
dual := proc (f::{polynom,list}) local J,K,k,Sym;
    if op(0, f) = `+` or type(f,'list') then
        return map(dual, f);
    elif op(0, f) = `*` then
        if type(op(1, f), 'numeric') then
            return op(1, f)*dual(f/op(1, f));
        else
            print("cannot compute dual :",f);
        end if;
    elif type(f, 'indexed') then
        if op(0,f)=`S` then
            Sym := T;
        else
            Sym := S;
        end if;
        k := 0;
        J := op(1,f);
        K := [];
        while nops(J)>0 do
            k := k+1;
            if J[1]=1 then
                J := J[2..-1];
                if nops(J)=0 then
                    K := [op(K),k];
                end if;
            else
                J[1] := J[1]-1;
                K := [op(K),k];
                k := 0;
            end if;
        end do;
        return Sym[K];
    else
            print("cannot compute dual:",f);
    end if;
end proc:

max_stars_repo_name: shizuo-kaji/Leibniz-Hopf-algebra

max_issues_repo_name: shizuo-kaji/Leibniz-Hopf-algebra

max_forks_repo_name: shizuo-kaji/Leibniz-Hopf-algebra


text:
(*
 Copyright (C) 2017 M.A.L. Marques

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*)

(* type: gga_exc *)

cac  := 1.467:
mtau := 4.5:
bcgp_pt := t -> t*sqrt((mtau + t)/(mtau + cac*t)):
arg := 0.5:
crg := 0.16667:
drg := 0.29633:
fbeta_modRasoltGeldart := rs -> (1.0 + arg*rs*(1+crg*rs)) / (1.0 + arg*rs*(1+drg*rs)):

$define gga_c_pbe_params
$include "gga_c_pbe.mpl"

(* override definition of tp *)
tp := (rs, z, xt) -> bcgp_pt(tt(rs, z, xt)):

(* override definition of mbeta *)
mbeta := (rs, t) -> params_a_beta*fbeta_modRasoltGeldart(rs):


text:
MoreThuente := module()

  description "More-Thuente";

  # Module defined as a package (i.e.) collection of procedures
  option package, load = ModuleLoad, unload = ModuleUnLoad;

  export lines, linesearch, esatta, AB_list, doplot, doplot2;

  local ModuleLoad,
        ModuleUnLoad,
        _debug,
        min_quadratic,
        min_cubic,
        xtol, ftol, gtol, stpmin, stpmax, maxfev,
        stx, fx, dx, sty, fy, dy, stp, fp, dp, brackt,
        max_iter, A_list, B_list;

  uses LinearAlgebra;

  ModuleUnLoad := proc()
    printf("PotraShi unload\n");
    NULL;
  end proc:

  ModuleLoad := proc()
    printf("PotraShi load\n");
    _debug   := true;
    rho      := 0.25;
    sigma    := 0.75;
    tau1     := 0.1;
    tau2     := 0.49;
    tau3     := 2.5;
    Jbig     := 2;
    max_iter := 40;
    NULL;
  end proc;

  # Explicitly call ModuleLoad here so the type is registered when this
  # code is cut&pasted in.  ModuleLoad gets called when the module is
  # read from the repository, but the code used to define a module
  # (like the command below) is not called at library read time.
  ModuleLoad();

  AB_list := proc()
    A_list, B_list;
  end proc;

  min_quadratic := proc ( a, fa, Dfa, b, fb )
    local t, delta; # A+B*(x-a)+C*(x-a)^2
    # p  = fa + Dfa*(x-a)+t*(x-a)^2
    # p' = Dfa+2*t*(x-a)
    ##t     := ((fb-fa)/(b-a)-Dfa)/(b-a);
    ##delta := -Dfa/(2*t);
    t     := (fb-fa)/(b-a)-Dfa;
    delta := max(tau1,min(1-tau1,-Dfa/(2*t)));
    a+(b-a)*delta;
  end proc;

  lines := proc( f0, df0, L )
    [ [0,f0], [L,f0+L*rho*df0] ], [ [0,f0], [L,f0+L*sigma*df0] ];
  end proc;

  pqr := proc()
    theta := 3*(fx - fp)/(stp - stx) + dx + dp;
    s     := max([abs(theta),abs(dx),abs(dp)]);
    gamma := s*sqrt(max(0,(theta/s)^2 - (dx/s)*(dp/s)));
    if stp > stx then
      gamma = -gamma;
    end;
    p := (gamma - dp) + theta;
    q := (gamma + (dx - dp)) + gamma;
    r := p/q;
  end proc:
Skip to content
Search or jump to…
Pull requests
Issues
Marketplace
Explore

@ebertolazzi
JuliaNLSolvers
/
LineSearches.jl
Public
6
69
26
Code
Issues
17
Pull requests
1
Actions
Projects
Wiki
Security
Insights
LineSearches.jl/src/morethuente.jl
@pkofod
pkofod Pkm/fix lse (#146)
…
Latest commit 0f5571e on Aug 18, 2020
 History
 6 contributors
@anriseth@pkofod@antoine-levitt@femtocleaner@mohamed82008@cortner
658 lines (594 sloc)  20.1 KB

#   Translation of Matlab version by John Myles White
#   Translation of minpack subroutine cvsrch
#   Dianne O'Leary   July 1991
#
#     **********
#
#     Subroutine cvsrch
#
#     The purpose of cvsrch is to find a step which satisfies
#     a sufficient decrease condition and a curvature condition.
#     The user must provide a subroutine which calculates the
#     function and the gradient.
#
#     At each stage the subroutine updates an interval of
#     uncertainty with endpoints stx and sty. The interval of
#     uncertainty is initially chosen so that it contains a
#     minimizer of the modified function
#
#          f(x + stp * s) - f(x) - f_tol * stp * (gradf(x)' * s).
#
#     If a step is obtained for which the modified function
#     has a nonpositive function value and nonnegative derivative,
#     then the interval of uncertainty is chosen so that it
#     contains a minimizer of f(x + stp * s).
#
#     The algorithm is designed to find a step which satisfies
#     the sufficient decrease condition
#
#           f(x + stp * s) <= f(x) + f_tol * stp * (gradf(x)' * s),
#
#     and the curvature condition
#
#           abs(gradf(x + stp * s)' * s)) <= gtol * abs(gradf(x)' * s).
#
#     If f_tol is less than gtol and if, for example, the function
#     is bounded below, then there is always a step which satisfies
#     both conditions. If no step can be found which satisfies both
#     conditions, then the algorithm usually stops when rounding
#     errors prevent further progress. In this case stp only
#     satisfies the sufficient decrease condition.
#
#     The subroutine statement is
#
#        subroutine cvsrch(df,n,x,f,s,stp,f_tol,gtol,x_tol,
#                          alphamin,alphamax,maxfev,info,nfev,wa)
#
#     where
#
# df is the name of the user-supplied subroutine which
#    calculates the function and the gradient.  df must
#    be declared in an external statement in the user
#    calling program, and should be written as follows.
#
#    function [f,g] = df(n,x) (Matlab)
#                     (10/2010 change in documentation)
#                     (derived from Fortran subroutine df(n,x,f,g))
#    integer n
#    f
#    x(n),g(n)
#
#    Calculate the function at x and
#    return this value in the variable f.
#    Calculate the gradient at x and
#    return this vector in g.
#
#  n is a positive integer input variable set to the number
#   of variables.
#
# x is an Abstractarray of length n. On input it must contain the
#   base point for the line search. On output it contains
#    x + stp * s.
#
# f is a variable. On input it must contain the value of f
#    at x. On output it contains the value of f at x + stp * s.
#
# g is an Abstractarray of length n. On input it must contain the
#    gradient of f at x. On output it contains the gradient
#    of f at x + stp * s.
#
# s is an input Abstractarray of length n which specifies the
#    search direction.
#
# stp is a nonnegative variable. On input stp contains an
#    initial estimate of a satisfactory step. On output
#    stp contains the final estimate.
#
#  f_tol and gtol are nonnegative input variables. Termination
#    occurs when the sufficient decrease condition and the
#    directional derivative condition are satisfied.
#
# x_tol is a nonnegative input variable. Termination occurs
#    when the relative width of the interval of uncertainty
#   is at most x_tol.
#
# alphamin and alphamax are nonnegative input variables which
#   specify lower and upper bounds for the step.
#
# maxfev is a positive integer input variable. Termination
#    occurs when the number of calls to df is at least
#    maxfev by the end of an iteration.
#
# info is an integer output variable set as follows:
#
#   info = 0  Improper input parameters.
#
#   info = 1  The sufficient decrease condition and the
#              directional derivative condition hold.
#
#   info = 2  Relative width of the interval of uncertainty
#            is at most x_tol.
#
#   info = 3  Number of calls to df has reached maxfev.
#
#   info = 4  The step is at the lower bound alphamin.
#
#   info = 5  The step is at the upper bound alphamax.
#
#   info = 6  Rounding errors prevent further progress.
#              There may not be a step which satisfies the
#              sufficient decrease and curvature conditions.
#              Tolerances may be too small.
#
#    nfev is an integer output variable set to the number of
#         calls to df.
#
#     Argonne National Laboratory. MINPACK Project. June 1983
#     Jorge J. More', David J. Thuente
#
#     **********

# Returns x, f, stp, info, nfev
# TODO: Decide whether to update x, f, g and info
#       or just return step and nfev and let existing code do its job

"""
The line search implementation from:
  Moré, Jorge J., and David J. Thuente
    Line search algorithms with guaranteed sufficient decrease.
    ACM Transactions on Mathematical Software (TOMS) 20.3 (1994): 286-307.
"""
@with_kw struct MoreThuente{T}
    f_tol::T = 1e-4 # c_1 Wolfe sufficient decrease condition
    gtol::T = 0.9   # c_2 Wolfe curvature condition (Recommend 0.1 for GradientDescent)
    x_tol::T = 1e-8
    alphamin::T = 1e-16
    alphamax::T = 65536.0
    maxfev::Int = 100
end

function (ls::MoreThuente)(df::AbstractObjective, x::AbstractArray{T},
                           s::AbstractArray{T}, alpha::Real, x_new::AbstractArray{T},
                           ϕ_0, dϕ_0) where T
    ϕdϕ = make_ϕdϕ(df, x_new, x, s)
    ls(ϕdϕ, alpha, ϕ_0, dϕ_0)
end

(ls::MoreThuente)(ϕ, dϕ, ϕdϕ, alpha, ϕ_0, dϕ_0) = ls(ϕdϕ, alpha, ϕ_0, dϕ_0)

# TODO: Should we deprecate the interface that only uses the ϕdϕ argument?
function (ls::MoreThuente)(ϕdϕ,
                           alpha::T,
                           ϕ_0,
                           dϕ_0) where T
    @unpack f_tol, gtol, x_tol, alphamin, alphamax, maxfev = ls

    iterfinitemax = -log2(eps(T))
    info = 0
    info_cstep = 1 # Info from step

    zeroT = convert(T, 0)

    #
    # Check the input parameters for errors.
    #

    if  alpha <= zeroT || f_tol < zeroT || gtol < zeroT ||
        x_tol < zeroT || alphamin < zeroT || alphamax < alphamin || maxfev <= zeroT
        throw(LineSearchException("Invalid parameters to MoreThuente.", 0))
    end

    if dϕ_0 >= zeroT
        throw(LineSearchException("Search direction is not a direction of descent.", 0))
    end

    #
    # Initialize local variables.
    #

    bracketed = false
    stage1 = true
    nfev = 0
    finit = ϕ_0
    dgtest = f_tol * dϕ_0
    width = alphamax - alphamin
    width1 = 2 * width

    # Keep this across calls

    #
    # The variables stx, fx, dgx contain the values of the step,
    # function, and directional derivative at the best step.
    # The variables sty, fy, dgy contain the value of the step,
    # function, and derivative at the other endpoint of
    # the interval of uncertainty.
    # The variables alpha, f, dg contain the values of the step,
    # function, and derivative at the current step.
    #

    stx = zeroT
    fx = finit
    dgx = dϕ_0
    sty = zeroT
    fy = finit
    dgy = dϕ_0

    # START: Ensure that the initial step provides finite function values
    # This is not part of the original FORTRAN code
    if !isfinite(alpha)
        alpha = one(T)
    end
    stmin = stx
    stmax = alpha + 4 * (alpha - stx) # Why 4?
    alpha = max(alpha, alphamin)
    alpha = min(alpha, alphamax)

    f, dg = ϕdϕ(alpha)
    nfev += 1 # This includes calls to f() and g!()
    iterfinite = 0
    while (!isfinite(f) || !isfinite(dg)) && iterfinite < iterfinitemax
        iterfinite += 1
        alpha = alpha/2

        f, dg = ϕdϕ(alpha)
        nfev += 1 # This includes calls to f() and g!()

        # Make stmax = (3/2)*alpha < 2alpha in the first iteration below
        stx = (convert(T, 7)/8)*alpha
    end
    # END: Ensure that the initial step provides finite function values

    while true
        #
        # Set the minimum and maximum steps to correspond
        # to the present interval of uncertainty.
        #

        if bracketed
            stmin = min(stx, sty)
            stmax = max(stx, sty)
        else
            stmin = stx
            stmax = alpha + 4 * (alpha - stx) # Why 4?
        end

        #
        # Ensure stmin and stmax (used in cstep) don't violate alphamin and alphamax
        # Not part of original FORTRAN translation
        #
        stmin = max(alphamin,stmin)
        stmax = min(alphamax,stmax)

        #
        # Force the step to be within the bounds alphamax and alphamin
        #

        alpha = max(alpha, alphamin)
        alpha = min(alpha, alphamax)

        #
        # If an unusual termination is to occur then let
        # alpha be the lowest point obtained so far.
        #

        if (bracketed && (alpha <= stmin || alpha >= stmax)) ||
            nfev >= maxfev-1 || info_cstep == 0 ||
            (bracketed && stmax - stmin <= x_tol * stmax)
            alpha = stx
        end

        #
        # Evaluate the function and gradient at alpha
        # and compute the directional derivative.
        #
        f, dg = ϕdϕ(alpha)
        nfev += 1 # This includes calls to f() and g!()

        if isapprox(dg, 0, atol=eps(T)) # Should add atol value to MoreThuente
            return alpha, f
        end

        ftest1 = finit + alpha * dgtest

        #
        # Test for convergence.
        #

        # What does info_cstep stand for?

        if (bracketed && (alpha <= stmin || alpha >= stmax)) || info_cstep == 0
            info = 6
        end
        if alpha == alphamax && f <= ftest1 && dg <= dgtest
            info = 5
        end
        if alpha == alphamin && (f > ftest1 || dg >= dgtest)
            info = 4
        end
        if nfev >= maxfev
            info = 3
        end
        if bracketed && stmax - stmin <= x_tol * stmax
            info = 2
        end
        if f <= ftest1 && abs(dg) <= -gtol * dϕ_0
            info = 1
        end

        #
        # Check for termination.
        #

        if info != 0
            return alpha, f
        end

        #
        # In the first stage we seek a step for which the modified
        # function has a nonpositive value and nonnegative derivative.
        #

        if stage1 && f <= ftest1 && dg >= min(f_tol, gtol) * dϕ_0
            stage1 = false
        end

        #
        # A modified function is used to predict the step only if
        # we have not obtained a step for which the modified
        # function has a nonpositive function value and nonnegative
        # derivative, and if a lower function value has been
        # obtained but the decrease is not sufficient.
        #

        if stage1 && f <= fx && f > ftest1
            #
            # Define the modified function and derivative values.
            #
            fm = f - alpha * dgtest
            fxm = fx - stx * dgtest
            fym = fy - sty * dgtest
            dgm = dg - dgtest
            dgxm = dgx - dgtest
            dgym = dgy - dgtest
            #
            # Call cstep to update the interval of uncertainty
            # and to compute the new step.
            #
            stx, fxm, dgxm,
            sty, fym, dgym,
            alpha, fm, dgm,
            bracketed, info_cstep =
                cstep(stx, fxm, dgxm, sty, fym, dgym,
                      alpha, fm, dgm, bracketed, stmin, stmax)
            #
            # Reset the function and gradient values for f.
            #
            fx = fxm + stx * dgtest
            fy = fym + sty * dgtest
            dgx = dgxm + dgtest
            dgy = dgym + dgtest
        else
            #
            # Call cstep to update the interval of uncertainty
            # and to compute the new step.
            #
            stx, fx, dgx,
            sty, fy, dgy,
            alpha, f, dg,
            bracketed, info_cstep =
                cstep(stx, fx, dgx, sty, fy, dgy,
                      alpha, f, dg, bracketed, stmin, stmax)
        end

        #
        # Force a sufficient decrease in the size of the
        # interval of uncertainty.
        #

        if bracketed
            if abs(sty - stx) >= (convert(T, 2)/3) * width1
                alpha = stx + (sty - stx)/2
            end
            width1 = width
            width = abs(sty - stx)
        end
    end # while
end # function

# Translation of minpack subroutine cstep
# Dianne O'Leary   July 1991
#
# Subroutine cstep
#
# The purpose of cstep is to compute a safeguarded step for
# a linesearch and to update an interval of uncertainty for
# a minimizer of the function.
#
# The parameter stx contains the step with the least function
# value. The parameter stp contains the current step. It is
# assumed that the derivative at stx is negative in the
# direction of the step. If bracketed is set true then a
# minimizer has been bracketed in an interval of uncertainty
# with endpoints stx and sty.
#
# The subroutine statement is
#
# subroutine cstep(stx, fx, dgx,
#                  sty, fy, dgy,
#                  stp, f, dg,
#                  bracketed, alphamin, alphamax, info)
#
# where
#
# stx, fx, and dgx are variables which specify the step,
#   the function, and the derivative at the best step obtained
#   so far. The derivative must be negative in the direction
#   of the step, that is, dgx and stp-stx must have opposite
#   signs. On output these parameters are updated appropriately
#
# sty, fy, and dgy are variables which specify the step,
#   the function, and the derivative at the other endpoint of
#   the interval of uncertainty. On output these parameters are
#   updated appropriately
#
# stp, f, and dg are variables which specify the step,
#   the function, and the derivative at the current step.
#   If bracketed is set true then on input stp must be
#   between stx and sty. On output stp is set to the new step
#
# bracketed is a logical variable which specifies if a minimizer
#   has been bracketed. If the minimizer has not been bracketed
#   then on input bracketed must be set false. If the minimizer
#   is bracketed then on output bracketed is set true
#
# alphamin and alphamax are input variables which specify lower
#   and upper bounds for the step
#
# info is an integer output variable set as follows:
#   If info = 1,2,3,4,5, then the step has been computed
#   according to one of the five cases below. Otherwise
#   info = 0, and this indicates improper input parameters
#
# Argonne National Laboratory. MINPACK Project. June 1983
# Jorge J. More', David J. Thuente

function cstep(stx::Real, fx::Real, dgx::Real,
               sty::Real, fy::Real, dgy::Real,
               alpha::Real, f::Real, dg::Real,
               bracketed::Bool, alphamin::Real, alphamax::Real)

   T = promote_type(typeof(stx), typeof(fx), typeof(dgx), typeof(sty), typeof(fy), typeof(dgy), typeof(alpha), typeof(f), typeof(dg), typeof(alphamin), typeof(alphamax))
   zeroT = convert(T, 0)
   info = 0

   #
   # Check the input parameters for error
   #

   if (bracketed && (alpha <= min(stx, sty) || alpha >= max(stx, sty))) ||
     dgx * (alpha - stx) >= zeroT || alphamax < alphamin
       throw(ArgumentError("Minimizer not bracketed"))
   end

   #
   # Determine if the derivatives have opposite sign
   #

   sgnd = dg * (dgx / abs(dgx))

   #
   # First case. A higher function value.
   # The minimum is bracketed. If the cubic step is closer
   # to stx than the quadratic step, the cubic step is taken,
   # else the average of the cubic and quadratic steps is taken
   #

   if f > fx
      info = 1
      bound = true
      theta = 3 * (fx - f) / (alpha - stx) + dgx + dg
      # Use s to prevent overflow/underflow of theta^2 and dgx * dg
      s = max(abs(theta), abs(dgx), abs(dg))
      gamma = s * sqrt((theta / s)^2 - (dgx / s) * (dg / s))
      if alpha < stx
          gamma = -gamma
      end
      p = gamma - dgx + theta
      q = gamma - dgx + gamma + dg
      r = p / q
      alphac = stx + r * (alpha - stx)
      alphaq = stx + (dgx / ((fx - f) / (alpha - stx) + dgx)) / 2 * (alpha - stx)
      if abs(alphac - stx) < abs(alphaq - stx)
         alphaf = alphac
      else
         alphaf = (alphac + alphaq) / 2
      end
      bracketed = true

   #
   # Second case. A lower function value and derivatives of
   # opposite sign. The minimum is bracketed. If the cubic
   # step is closer to stx than the quadratic (secant) step,
   # the cubic step is taken, else the quadratic step is taken
   #

elseif sgnd < zeroT
      info = 2
      bound = false
      theta = 3 * (fx - f) / (alpha - stx) + dgx + dg
      # Use s to prevent overflow/underflow of theta^2 and dgx * dg
      s = max(abs(theta), abs(dgx), abs(dg))
      gamma = s * sqrt((theta / s)^2 - (dgx / s) * (dg / s))

      if alpha > stx
         gamma = -gamma
      end
      p = gamma - dg + theta
      q = gamma - dg + gamma + dgx
      r = p / q
      alphac = alpha + r * (stx - alpha)
      alphaq = alpha + (dg / (dg - dgx)) * (stx - alpha)
      if abs(alphac - alpha) > abs(alphaq - alpha)
         alphaf = alphac
      else
         alphaf = alphaq
      end
      bracketed = true

   #
   # Third case. A lower function value, derivatives of the
   # same sign, and the magnitude of the derivative decreases.
   # The cubic step is only used if the cubic tends to infinity
   # in the direction of the step or if the minimum of the cubic
   # is beyond alpha. Otherwise the cubic step is defined to be
   # either alphamin or alphamax. The quadratic (secant) step is also
   # computed and if the minimum is bracketed then the the step
   # closest to stx is taken, else the step farthest away is taken
   #

   elseif abs(dg) < abs(dgx)
      info = 3
      bound = true
      theta = 3 * (fx - f) / (alpha - stx) + dgx + dg
      # Use s to prevent overflow/underflow of theta^2 and dgx * dg
      s = max(abs(theta), abs(dgx), abs(dg))
      #
      # The case gamma = 0 only arises if the cubic does not tend
      # to infinity in the direction of the step
      #
      # # Use NaNMath in case s == zero(s)
      gamma = s * sqrt(NaNMath.max(zero(s), (theta / s)^2 - (dgx / s) * (dg / s)))

      if alpha > stx
          gamma = -gamma
      end
      p = gamma - dg + theta
      q = gamma + dgx - dg + gamma
      r = p / q
      if r < zeroT && gamma != zeroT
         alphac = alpha + r * (stx - alpha)
     elseif alpha > stx
         alphac = alphamax
      else
         alphac = alphamin
      end
      alphaq = alpha + (dg / (dg - dgx)) * (stx - alpha)
      if bracketed
         if abs(alpha - alphac) < abs(alpha - alphaq)
            alphaf = alphac
         else
            alphaf = alphaq
         end
      else
         if abs(alpha - alphac) > abs(alpha - alphaq)
            alphaf = alphac
         else
            alphaf = alphaq
         end
      end

   #
   # Fourth case. A lower function value, derivatives of the
   # same sign, and the magnitude of the derivative does
   # not decrease. If the minimum is not bracketed, the step
   # is either alphamin or alphamax, else the cubic step is taken
   #

   else
      info = 4
      bound = false
      if bracketed
         theta = 3 * (f - fy) / (sty - alpha) + dgy + dg
         # Use s to prevent overflow/underflow of theta^2 and dgy * dg
         s = max(abs(theta), abs(dgy), abs(dg))
         gamma = s * sqrt((theta / s)^2 - (dgy / s) * (dg / s))

         if alpha > sty
             gamma = -gamma
         end
         p = gamma - dg + theta
         q = gamma - dg + gamma + dgy
         r = p / q
         alphac = alpha + r * (sty - alpha)
         alphaf = alphac
     elseif alpha > stx
         alphaf = alphamax
      else
         alphaf = alphamin
      end
   end

   #
   # Update the interval of uncertainty. This update does not
   # depend on the new step or the case analysis above
   #

   if f > fx
      sty = alpha
      fy = f
      dgy = dg
   else
      if sgnd < zeroT
         sty = stx
         fy = fx
         dgy = dgx
      end
      stx = alpha
      fx = f
      dgx = dg
   end

   #
   # Compute the new step and safeguard it
   #

   alphaf = min(alphamax, alphaf)
   alphaf = max(alphamin, alphaf)
   alpha = alphaf
   if bracketed && bound
      if sty > stx
         alpha = min(stx + (convert(T, 2)/3) * (sty - stx), alpha)
      else
         alpha = max(stx + (convert(T, 2)/3) * (sty - stx), alpha)
      end
   end

   return stx, fx, dgx, sty, fy, dgy, alpha, f, dg, bracketed, info
end

  esatta := proc( fun )
    local res, alpha;
    #res := Optimization[Minimize](
    res := Optimization[NLPSolve](
      fun(alpha),
      {alpha >= 0, alpha <= 3},
      #initialpoint={alpha=1},
      iterationlimit=50
    ):
    subs(res[2],alpha);
  end proc:

  doplot := proc( f_in, df_in, alpha, amax)
    local f0, df0, x, AA, BB, CC, DD;
    f0  := evalf(f_in(0));
    df0 := evalf(df_in(0));

    AA := plot(f_in(x),x=0..amax);
    BB := plot([[0,f0],[amax,f0+amax*rho*df0]],color="LimeGreen");
    CC := plot([[0,f0],[amax,f0+amax*sigma*df0]],color="blue");
    DD := plot([[alpha,f_in(alpha)]],color="red",style = point);
    display(AA,BB,CC,DD);
  end;

  doplot2 := proc( f_in, df_in, alpha, amax, A_list, B_list )
    local i, f0, df0, x, AA, BB, CC, DD, EE;
    f0  := evalf(f_in(0));
    df0 := evalf(df_in(0));

    AA := plot(f_in(x),x=0..amax);
    BB := plot([[0,f0],[amax,f0+amax*rho*df0]],color="LimeGreen");
    CC := plot([[0,f0],[amax,f0+amax*sigma*df0]],color="blue");
    DD := plot([[alpha,f_in(alpha)]],color="red",style = point);
    EE := plot([seq([[A_list[i],-i/5],[B_list[i],-i/5]],i=1..nops(A_list))],color="black");
    display(AA,BB,CC,DD,EE);
  end;

end module:

max_stars_repo_name: ebertolazzi/NLtoolbox

max_issues_repo_name: ebertolazzi/NLtoolbox

max_forks_repo_name: ebertolazzi/NLtoolbox


text:
     {   ProdPlan.mpl   }

TITLE
    Production_Planning;

INDEX
    product  =  1..5
    month    =  (Jan,Feb,Mar,Apr,May,Jun)
    machine  =  (Grind,Vdrill,Hdrill,Boring,Planing) : 2

DATA
    price[product]          =        SPARSEFILE(price.dat);
    time[machine,product]   =  1/100 DATAFILE(time.dat);

    MaxSales[month,product] =        DATAFILE(sales.dat);

    DaysPerMonth   =  6*4;
    HoursPerDay    =  2*8;
    HoursPerMonth  =  HoursPerDay * DaysPerMonth;

    NrOfMachines[machine]   =  (4,2,3,3,3);

    Repairs[machine,month]  =  [Grind,   Jan:  1
                                Grind,   May:  1
                                Vdrill,  Apr:  1
                                Vdrill,  May:  1
                                Hdrill,  Feb:  2
                                Hdrill,  Jun:  1
                                Boring,  Mar:  1
                                Boring,  May:  1
                                Planing, Jun:  1]

    MachineHours[machine,month] :=

         HoursPerMonth * DENSE(NrOfMachines - Repairs) ;

    InventoryCost   =   0.50 ;

DECISION VARIABLES
    Sale[product,month]
    Prod[product,month]
    Invt[product,month] ;

MODEL

   MAX   Contribution   =

           SUM(product,month: price*Sale - InventoryCost*Invt) ;

SUBJECT TO

   Capacity[month,machine] -> Cap :  SUM(product: time*Prod)  <  MachineHours ;

   Balance[product,month] -> Bal  :  Invt  =  Invt[month-1] + Prod - Sale ;

BOUNDS
    MaximumSales  :  Sale  <  MaxSales;
    InvtCapacity  :  Invt[month<Jun]  <  100 ;
    InvtCapacity  :  Invt[month=Jun]  =   50 ;

END

max_stars_repo_name: BjarniMax/mpl-models

max_issues_repo_name: BjarniMax/mpl-models

max_forks_repo_name: BjarniMax/mpl-models


text:
(*
 Copyright (C) 2017 M.A.L. Marques

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*)

(* type: lda_exc *)

$define xc_dimensions_1d

eta0 := -log(sqrt(2*Pi)) - 3/4:
eta1 := 0.359933:
eps0 := -Pi^2/360:
eps1 := 0.00714:

kappa := 0.3083:

c0 := kappa*eta0:
c1 := 4*kappa*eta0 + kappa*sqrt(kappa)*eta1:
c2 := 5*eps0 + eps1/kappa:
c3 := eps1:

t := rs -> (sqrt(1 + 4*kappa*rs) - 1)/(2*kappa*rs):

f := (rs, z) ->
  t(rs)^2*(c0*(1 - t(rs))^3 + c1*t(rs)*(1 - t(rs))^2 + c2*t(rs)^2*(1 - t(rs)) + c3*t(rs)^3):


text:
# This code is for the cohomology of F(A,R^3), or F(A,R^{2n+1}) for
# any n > 0.  This is concentrated in even degrees so we do not
# need to worry about sign rules or gradings.

with(Groebner):

`vars/poisson` := proc(A::set)
 local a,b;
 [seq(seq(u[a,b],b in A),a in A)];
end:

`rels/poisson` := proc(A::set)
 local n,i,j,k,uf;
 n := nops(A);
 uf := (i,j) -> u[A[i],A[j]];

 return
  [seq(uf(i,i),i = 1..n),
   seq(seq(uf(i,j)+uf(j,i),j=i+1..n),i=1..n-1),
   seq(seq(uf(i,j)^2,j=1..n),i=1..n),
   seq(seq(seq(
    uf(i,j)*uf(j,k) + uj(j,k)*uf(k,i) + uf(k,i)*uf(i,j),
    k=j+1..n),j=i+1..n),i=1..n)
  ];
end:

`groebner_basis/poisson` := proc(A::set)
 Basis(`rels/poisson`(A),tdeg(`vars/poisson`(A)));
end:

`basis/poisson` := proc(A::set)
 local n,a,x,v,B,U,V;

 n := nops(A);
 if n <= 1 then
  return [1];
 fi;
 a := A[n];
 B := A minus {a};
 U := `basis/poisson`(A minus {a});
 V := [1,seq(u[x,a],x in B)];
 return map(op,[seq(U *~ v,v in V)]);
end:


text:
program IfstTC;	{sample25t}
var
 ch:char;int:integer;boolx,booly:boolean;
begin
boolx:=true;booly:=false;ch:='a';int:=66;
write(integer(ch));write(integer(int));write(integer(boolx),integer(booly));
writeln;
writeln(char(ch),char(int));writeln(boolean(ch),integer(int),integer(boolx),
integer(booly));end.

max_forks_repo_name: Hiroya-W/lang-processing


text:
fun:=proc()
    uses Grid;
    local node,i,N;
    node:=MyNode();
    N:=NumNodes();
    if node<>0 then
        Threads:-Sleep(N-node);
        print(node);
        Send(0,evalf(Pi-3+node,100));
    else
        for i from 1 to N-1 do
            print(Receive(i));
        end do;
    end if;
end proc:
Grid:-Launch(fun);


text:
input := FileTools:-Text:-ReadFile("AoC-2021-9-input.txt" ):

ogrid := ((parse~)@StringTools:-Explode)~(subs(""=NULL,StringTools:-Split(input))):
gridw := nops(ogrid[1]);
gridl := nops(ogrid);
grid := [ [ 9 $ (gridw+2) ], map(l->[9,l[],9], ogrid)[],  [ 9 $ (gridw+2) ] ]:

lowpoints := NULL:
risklevel := 0:
for i from 2 to gridl+1 do
    for j from 2 to gridw+1 do
        if      grid[i][j] < grid[i-1][j]
            and grid[i][j] < grid[i+1][j]
            and grid[i][j] < grid[i][j-1]
            and grid[i][j] < grid[i][j+1]
        then
            lowpoints := lowpoints, [i,j];
            risklevel := risklevel + grid[i][j]+1;
        end if;
    end do;
end do:

answer1 := risklevel;

counted := table(sparse=0):
basinsize := proc( input )
global grid, counted;
local nb, x, y, v;

    x, y, v := op(input);

    if v = 9 or counted[input] <> 0 then
        return 0;
    end if;

    nb := [[x, y-1], [x, y+1], [x-1, y], [x+1,y]];
    nb := remove( n->(min(n)<2 or n[1]>gridl+2 or n[2]>gridw+2), nb);
    nb := map(n->[n[1], n[2], grid[n[1]][n[2]]], nb);
    #nb := remove(n->n[3]=9 or n[3] < v, nb); # for non-9 bounded basins

    counted[ input ] := 1;

    if nb = [] then
        return 1;
    else
        return add(map(thisproc, nb))+1;
    end if;
end proc:

bs := map(l->basinsize([l[], grid[l[1]][l[2]]]), [lowpoints]):
answer2 := mul( sort(bs)[-3..-1] );


text:
use DocumentTools in
	if parse(GetProperty("Refine", value)) then
		SetProperty("Refine", value, false):
		SetProperty("UsingUpTo", enabled, false):
		SetProperty("MaxPermutations", enabled, false):
		SetProperty("Permutations", enabled, false):
	fi:
end use;

max_stars_repo_name: iliailmer/sian-web-app

max_issues_repo_name: iliailmer/sian-web-app

max_forks_repo_name: iliailmer/sian-web-app


text:
# The goal here is to set up a simplicial complex K and an embedding of
# |K| in R^3 such that the image is a version of the standard "pair of
# pants space", or in other words a sphere with three discs removed.
# More specifically, we start with the unit sphere in R^3, and remove
# three discs of radius 1/2 centred at points equally spaced around
# the plane perpendicular to [1,1,1].  This has an evident linear action
# of the group G = S_3 x C_2, and everything will be equivariant with
# respect to this.  There is a fundamental domain for the action on S^2
# that is a spherical triangle with one vertex at [1,1,1]/sqrt(3) and
# the other two vertices separated by Pi/3 on the equator perpendicular
# to [1,1,1].  One boundary circle meets this domain in a quarter
# circle, and the other boundary circles do not intersect it at all.
# We start by setting up a triangulation of the fundamental domain
# such that this quarter circle is essentially a subcomplex.  We
# can then propagate this equivariantly and remove the boundary
# discs to obtain a triangulated surface of the required homeomorphism
# type.  We can glue a standard horizontal cylinder to each boundary
# circle to connect to other pieces and so construct various surfaces.
# However, this is ugly because the tangent spaces do not line up
# along the glued circles.  To fix this, we construct a map
# f : S^2 -> R^3 with the property the f is the identity on our
# boundary circles but acts near those circles in such a way as to
# line up the tangent spaces with a horizontal cylinder.

lib_read("brent.mpl"):

make_pants_complex := proc(verbose := false)
 local plot_opts,dp,G,T,U,V,xx,i,j,k,lm,m,v,g,gi,gx,E,R,f2,sphere_retract,
  xa,xb,xc,ua,rb,rels,F,a,n,s,t,u;
 global pants_group,pants_sphere_sector_complex,pants_sphere_complex,
  pants_complex,condensed_pants_complex;

 plot_opts := scaling=constrained, axes=none:
 dp := `dot/R`(3);

 if verbose then printf("Making group\n"); fi;

 # This is the group S_3 x C_2, which will act isometrically on everything.
 G := table():
 pants_group := eval(G);
 G["elements"] := [seq(seq([op(u),s],u in combinat[permute](3)),s in [1,-1])]:
 G["id"] := [1,2,3,1]:
 G["o"] := (u,v) -> [u[v[1]],u[v[2]],u[v[3]],u[4] * v[4]]:
 G["inv"] := proc(u)
  local ui,i;
  ui := table():
  for i from 1 to 3 do ui[u[i]] := i; od:
  return [ui[1],ui[2],ui[3],u[4]];
 end:

 if verbose then printf("Making skeleton\n"); fi;
 T := table():
 pants_sphere_sector_complex := eval(T):

 T["u"] := table([
   1 = [ 2,-1,-1] /~ sqrt(6),
   2 = [ 1, 1,-2] /~ sqrt(6),
   3 = [-1, 2,-1] /~ sqrt(6),
   4 = [-2, 1, 1] /~ sqrt(6),
   5 = [-1,-1, 2] /~ sqrt(6),
   6 = [ 1,-2, 1] /~ sqrt(6),
 ##
   7 = [ 1, 1, 1] /~ sqrt(3),
   8 = [-1,-1,-1] /~ sqrt(3),
 ##
   9 = [ 0, 1,-1] /~ sqrt(2),
  10 = [-1, 1, 0] /~ sqrt(2),
  11 = [-1, 0, 1] /~ sqrt(2),
  12 = [ 0,-1, 1] /~ sqrt(2),
  13 = [ 1,-1, 0] /~ sqrt(2),
  14 = [ 1, 0,-1] /~ sqrt(2)
 ]):

 T["u_index"] := table():
 for i from 0 to 13 do
  T["u_index"][T["u"][i]] := i;
  T["u_index"][simplify(T["u"][i])] := i;
 od:

 xx := [x[1],x[2],x[3]];

 T["act_R3"] := table():
 for g in G["elements"] do
  gi := G["inv"](g);
  gx := [x[gi[1]],x[gi[2]],x[gi[3]]];
  if g[4] = -1 then
   gx := expand(gx -~ (2/3 * dp([1,1,1],gx)) *~ [1,1,1]);
  fi;
  T["act_R3"][g] := unapply(gx,x);
 od:

 T["act_u_table"] := table():
 T["act_u"] := table():
 for g in G["elements"] do
  T["act_u_table"][g] :=
   table([seq(i = T["u_index"][simplify(T["act_R3"][g](T["u"][i]))],i=0..13)]);
  T["act_u"][g] := (i) -> T["act_u_table"][g][i];
 od:

 T["equator"] := unapply(cos(t) *~ T["u"][1] +~ sin(t) *~ T["u"][9],t):
 T["equator_plot"] :=
  spacecurve(T["equator"](t),t=0..2*Pi,colour=green,plot_opts):
 T["equator_sector_plot"] :=
  spacecurve(T["equator"](t),t=0..Pi/3,colour=green,plot_opts):

 T["longitude"] := table():
 T["longitude_plot"] := table():
 for i from 1 to 6 do
  T["longitude"][i] :=
   unapply(expand(sin(t) *~ T["u"][7] +~ cos(t) *~ T["u"][i]),t);
  T["longitude_plot"][i] :=
   spacecurve(T["longitude"][i](t),t=-Pi/2..Pi/2,colour=blue,plot_opts):
 od:
 T["longitude_sector_plot"] :=
  display(
   spacecurve(T["longitude"][1](t),t=0..Pi/2,colour=blue,plot_opts),
   spacecurve(T["longitude"][2](t),t=0..Pi/2,colour=blue,plot_opts),
   plot_opts
  ):

 T["boundary_circle"] := table():
 T["boundary_circle_plot"] := table():
 for i from 1 to 3 do
  T["boundary_circle"][i] :=
   unapply(expand((sqrt(3)/2) *~ T["u"][2*i-1] +~ (cos(t)/2) *~ T["u"][2*i+7] +~ (sin(t)/2) *~ T["u"][7]),t);
  T["boundary_circle_plot"][i] :=
   spacecurve(T["boundary_circle"][i](t),t=0..2*Pi,colour=red,plot_opts)
 od:
 T["boundary_circle_sector_plot"] :=
   spacecurve(T["boundary_circle"][1](t),t=0..Pi/2,colour=red,plot_opts):

 T["wire_sphere_plot"]  := sphere(colour=grey,style=wireframe):
 T["solid_sphere_plot"] := sphere(colour=grey,style=patchnogrid):
 T["wire_sector_plot"]  :=
  plot3d(cos(t) *~ T["u"][7] +~ sin(t) *~ (cos(u) *~ T["u"][1] +~ sin(u) *~ T["u"][9]),
	 t = 0 .. Pi/2,u = 0 .. Pi/3,
	 colour=grey,style=wireframe,scaling=constrained, axes=none):
 T["solid_sector_plot"]  :=
  plot3d(cos(t) *~ T["u"][7] +~ sin(t) *~ (cos(u) *~ T["u"][1] +~ sin(u) *~ T["u"][9]),
	 t = 0 .. Pi/2,u = 0 .. Pi/3,
	 colour=grey,style=patchnogrid,scaling=constrained, axes=none):

 T["line_plot"] := display(
  T["equator_plot"],
  seq(T["longitude_plot"][i],i=1..6),
  seq(T["boundary_circle_plot"][i],i=1..3)
 ):

 T["sphere_plot"] := display(
  T["solid_sphere_plot"],T["line_plot"],
  plot_opts
 ):

 T["line_sector_plot"] := display(
  T["equator_sector_plot"],
  T["longitude_sector_plot"],
  T["boundary_circle_sector_plot"],
  plot_opts
 ):

 T["sector_plot"] := display(
  T["solid_sector_plot"],T["line_sector_plot"],
  plot_opts
 ):

 lm := 4:
 m := 2 ^ lm:
 T["m"] := m:

 if verbose then printf("Making simplices\n"); fi;

 v := (i,j) -> [i/m,j/m,(m-i-j)/m]:

 T["vertices"] := [seq(seq(v(i,j),j=0..m - i),i=0..m)]:

 T["edges"] := [
  seq(seq([v(i,j),v(i+1,j)],j=0..m-1-i),i=0..m-1),
  seq(seq([v(i,j),v(i,j+1)],j=0..m-1-i),i=0..m-1),
  seq(seq([v(i+1,j),v(i,j+1)],j=0..m-1-i),i=0..m-1),
 NULL
 ]:

 T["faces"] := [
  seq(seq([v(i,j),v(i+1,j),v(i,j+1)],j=0..m-1-i),i=0..m-1),
  seq(seq([v(i,j+1),v(i+1,j),v(i+1,j+1)],j=0..m-2-i),i=0..m-2)
 ]:

 T["max_simplices"] := T["faces"]:

 if verbose then printf("Making explicit embedding\n"); fi;

 unassign('i');
 T["f0"] := unapply(i[1] *~ T["u"][1] +~ i[2] *~ T["u"][2] +~ i[3] *~ T["u"][7],i):
 T["f1"] := proc(i) local x; x := T["f0"](i); return x /~ sqrt(add(x[j]^2,j=1..3)); end:
 T["embedding0"] := table():
 T["embedding1"] := table():
 for i in T["vertices"] do
  T["embedding0"][i] := T["f0"](i);
  T["embedding1"][i] := T["f1"](i);
 od:

 E := table():
 for i in T["vertices"] do E[i] := NULL; od:

 for i from 0 to m do
   E[[m-i,i,0] /~ m] := T["equator"](i/m * Pi/3);
   E[[0,m-i,i] /~ m] := T["longitude"][2](i/m * Pi/2);
 od:
 for i from 0 to m/2 do
   E[[m-i,0,i]       /~ m] := T["longitude"][1]((2*i/m) * Pi/6);
   E[[m/2-i,0,m/2+i] /~ m] := T["longitude"][1]((1+4*i/m) * Pi/6);
   E[[m/2,m/2-i,i]   /~ m] := T["boundary_circle"][1]((2*i/m) * Pi/2);
 od:

 sphere_retract := (u) -> evalf(u /~ sqrt(add(u[i]^2,i=1..3))):

 `midpoint_extend/simplicial_complex`(E,4,sphere_retract):
 T["embedding2"] := eval(E):
 T["embedding"] := eval(E):
 `plot/simplicial_complex`(T):
 `surface_plot/simplicial_complex`(T):

 T["proj"] := table():
 for k from 1 to 3 do
  T["proj"][k] := unapply(expand(xx - dp(x,T["u"][2*k-1]) *~ T["u"][2*k-1]),x);
 od:

 T["proj_a"] := unapply([dp(xx,T["u"][ 1]),dp(xx,T["u"][12])],x):
 T["proj_b"] := unapply([dp(xx,T["u"][ 5]),dp(xx,T["u"][ 7])],x):
 T["proj_c"] := unapply([dp(xx,T["u"][ 1]),dp(xx,T["u"][ 7])],x):
 T["proj_d"] := unapply([dp(xx,T["u"][12]),dp(xx,T["u"][ 7])],x):

 unassign('s');

 T["unproj_a"] := unapply(s[1] *~ T["u"][ 1] +~ s[2] *~ T["u"][12], s);
 T["unproj_b"] := unapply(s[1] *~ T["u"][ 5] +~ s[2] *~ T["u"][ 7], s);
 T["unproj_c"] := unapply(s[1] *~ T["u"][ 1] +~ s[2] *~ T["u"][ 7], s);
 T["unproj_d"] := unapply(s[1] *~ T["u"][12] +~ s[2] *~ T["u"][ 7] +~ (sqrt(3)/2) *~ T["u"][1], s);

 T["g"] := unapply(expand(mul(dp(x,T["u"][2*j-1]) - sqrt(3)/2,j=1..3)),x):
 T["h"] := unapply(expand((2 + dp(x,T["u"][7]) ^ 2)*sqrt(3)/4),x):

 # Here we define a particular polynomial function k(x), and then we
 # define f(x) in terms of k(x).  Most of the desired properties hold
 # independently of the choice of k(x).  Our particular choice is
 # designed to straighten out the tangent spaces on the boundary
 # circles and ensure that the width of he surface in the centre
 # is the same as the width of the circles.
 T["k"] := (x) -> (20 + 5 * sqrt(2) * x[1] - 15 * x[1]^2 + 5 * sqrt(2) * x[1]^3 +
		   (4 * x[1] - 3 * sqrt(2)) * (x[2] + x[3]) +
		   (8 * sqrt(2) * x[1] - 26) * x[2] * x[3] +
		   (-5) * sqrt(2) * (x[2] + x[3]) * x[2] * x[3])/27:
 T["f2"] := unapply(simplify(
   xx +~
    (T["g"](x)/T["h"](x)) *~ (
      T["k"]([x[1],x[2],x[3]]) *~ T["proj"][1](x) +~
      T["k"]([x[2],x[3],x[1]]) *~ T["proj"][2](x) +~
      T["k"]([x[3],x[1],x[2]]) *~ T["proj"][3](x))),
   x):

 for i in T["vertices"] do
  T["embedding3"][i] := evalf(T["f2"](T["embedding2"][i]));
 od:
 T["embedding"] := eval(T["embedding4"]);

 if verbose then printf("Making implicit embedding\n"); fi;

 # We now start to define an alternative embedding, using the surface
 # where the following function is zero.
 T["l"] := unapply(
  12*sqrt(3)*T["g"](xx) + (1-dp(xx,xx))*(27+6*(x[1]*x[2]+x[2]*x[3]+x[3]*x[1])+
  8*sqrt(3)*T["g"](xx)-6*dp(xx,xx)),x):
 T["l_min_root"] := fsolve(expand(T["l"]((-t) *~ T["u"][1])),t=0.68..0.70);
 T["l_max_root"] := fsolve(expand(T["l"](t *~ T["u"][1])),t=1.42..1.43);
 xa := (r * cos(t)) *~ T["u"][1] +~ (r * sin(t)) *~ T["u"][9] +~ u *~ T["u"][7]:
 ua := solve(simplify(expand(T["l"](xa))),u)[1]:
 ua := sqrt(factor(ua^2)):
 xa := subs(u = ua,xa):
 T["chart_a"] := unapply(xa,r,t):
 xb := (r * cos(t)) *~ T["u"][9] +~ (r * sin(t)) *~ T["u"][7] +~ u *~ T["u"][1]:
 rb := solve(simplify(expand(T["l"](xb))),r)[3]:
 rb := 1/sqrt(factor(expand(rationalize(1/rb^2)))):
 xb := subs(r = rb,xb):
 T["chart_b"] := unapply(xb,u,t):
 xc := (a[0] - a[1] * r^2 - a[2] * r^4) *~ T["u"][1] +~ r *~ (cos(t) *~ T["u"][9] +~ sin(t) *~ T["u"][7]):
 rels := expand([coeffs(collect(simplify(expand(convert(series(T["l"](xc),r=0,5),polynom,r))),r),r)]):
 xc := expand(evalf(subs(a[0] = T["l_max_root"],subs(solve([rels[2],rels[3]],{a[1],a[2]}),xc)))):
 T["chart_c"] := unapply(xc,r,t):

 T["dl"] := unapply([seq(diff(T["l"](xx),x[i]),i=1..3)],x);

 T["retract"] := proc(x)
  local g,y,sols;
  g := evalf(T["dl"](x));
  g := g /~ sqrt(add(g[i]^2,i=1..3));
  y := x +~ t *~ g;
  sols := [fsolve(evalf(expand(T["l"](y))),t)];
  sols := select(t -> abs(t) = min(map(abs,sols)),sols);
  if sols = [] then error("No solutions"); fi;
  return(evalf(subs(t = sols[1],y)));
 end:

 T["step"] := proc(q,d0,proj,unproj)
  local xx,st0,mx,dx,g,h,sol;
  xx  := eval(unproj([s,t]));
  st0 := eval(proj(q));
  mx := evalf(T["l"](xx));
  dx := (s - st0[1])^2 + (t - st0[2])^2 - d0^2;
  g := evalf(subs({s = st0[1],t = st0[2]},[-diff(mx,t),diff(mx,s)]));
  g := (d0 / sqrt(g[1]^2+g[2]^2)) *~ g;
  h := st0 +~ g;
  sol := fsolve({mx,dx},{s = h[1], t = h[2]});
  return evalf(subs(sol,xx));
 end:

 F[0] := (E) -> proc(d0)
  local i;
  E[[1,0,0]] := evalf(T["l_max_root"] *~ T["u"][1]);
  for i from 1 to 8 do
   E[[16-i,i,0] /~ 16] := T["step"](E[[17-i,i-1,0]/~16],d0,T["proj_a"],T["unproj_a"]):
  od:
  return table(["err" = evalf(sqrt(3)/2 - dp(T["u"][1],E[[1/2,1/2,0]]))]);
 end:

 F[1] := (E) -> proc(d0)
  global p; local i;
  E[[1/2,1/2,0]] := T["u"][14];
  for i from 1 to 8 do
   E[[8-i,8+i,0]/~16] := T["step"](E[[9-i,7+i,0]/~16],d0,T["proj_a"],T["unproj_a"]):
  od:
  return table(["err" = evalf(dp(T["u"][10],E[[0,1,0]]))]);
 end:

 F[2] := (E) -> proc(d0)
  global p; local i;
  E[[0,1,0]] := evalf(T["l_min_root"] *~ T["u"][2]):
  for i from 1 to 16 do
   E[[0,16-i,i]/~16] := T["step"](E[[0,17-i,i-1]/~16],d0,T["proj_b"],T["unproj_b"]):
  od:
  return table(["err" = evalf(-dp(T["u"][9],E[[0,0,1]]))]);
 end:

 F[3] := (E) -> proc(d0)
  global p; local i;
  E[[0,0,1]] := T["u"][7] /~ 2:
  for i from 1 to 8 do
   E[[i,0,16-i]/~16] := T["step"](E[[i-1,0,17-i]/~16],d0,T["proj_c"],T["unproj_c"]):
  od:
  return table(["err" = evalf(dp(E[[1/2,0,1/2]],T["u"][1]) - sqrt(3)/2)]);
 end:

 F[4] := (E) -> proc(d0)
  global p; local i,e;
  E[[1/2,0,1/2]] := (sqrt(3)/2) *~ T["u"][1] +~ (1/2) *~ T["u"][7];
  for i from 1 to 8 do
   E[[8+i,0,8-i]/~16] := T["step"](E[[7+i,0,9-i]/~16],d0,T["proj_c"],T["unproj_c"]):
  od:
  e := - evalf(dp(E[[1,0,0]],T["u"][7]));
  E[[1,0,0]] := evalf(T["l_max_root"] *~ T["u"][1]);
  return table(["err" = e]);
 end:

 F[5] := (E) -> proc(d0)
  global p; local i,e;
  for i from 1 to 8 do
   E[[8,8-i,i]/~16] := T["step"](E[[8,9-i,i-1]/~16],d0,T["proj_d"],T["unproj_d"]):
  od:
  e := - evalf(dp(E[[1/2,0,1/2]],T["u"][9]));
  E[[1/2,0,1/2]] := (sqrt(3)/2) *~ T["u"][1] +~ (1/2) *~ T["u"][7];
  return table(["err" = e]);
 end:

 E := table():
 for i in T["vertices"] do E[i] := NULL; od:
 T["F_step"] := table([0 = 0.10, 1 = 0.06, 2 = 0.05, 3 = 0.10, 4 = 0.10, 5 = 0.09]):
 T["F_error"] := table():
 for i from 0 to 5 do
  T["F_step"][i] :=
   brent_fsolve(F[i](E),T["F_step"][i],T["F_step"][i] + 0.01,false,false,10.^(-10))[1];
  T["F_error"][i] := F[i](E)(T["F_step"][i]);
 od:

 `midpoint_extend/simplicial_complex`(E,4,T["retract"]):

 T["embedding4"] := eval(E);

 if verbose then printf("Extending equivariantly\n"); fi;

 U := table():
 pants_sphere_complex := eval(U);
 U["act_i"] := table():
 U["act_i"][[1,2,3, 1]] := (j) -> [ j[1], j[2], j[3], j[4]]:
 U["act_i"][[2,3,1, 1]] := (j) -> [-j[2],-j[3], j[1], j[4]]:
 U["act_i"][[3,1,2, 1]] := (j) -> [ j[3],-j[1],-j[2], j[4]]:
 U["act_i"][[2,1,3, 1]] := (j) -> [ j[3], j[2], j[1], j[4]]:
 U["act_i"][[3,2,1, 1]] := (j) -> [-j[2],-j[1], j[3], j[4]]:
 U["act_i"][[1,3,2, 1]] := (j) -> [ j[1],-j[3],-j[2], j[4]]:
 U["act_i"][[1,2,3,-1]] := (j) -> [ j[1], j[2], j[3],-j[4]]:
 U["act_i"][[2,3,1,-1]] := (j) -> [-j[2],-j[3], j[1],-j[4]]:
 U["act_i"][[3,1,2,-1]] := (j) -> [ j[3],-j[1],-j[2],-j[4]]:
 U["act_i"][[2,1,3,-1]] := (j) -> [ j[3], j[2], j[1],-j[4]]:
 U["act_i"][[3,2,1,-1]] := (j) -> [-j[2],-j[1], j[3],-j[4]]:
 U["act_i"][[1,3,2,-1]] := (j) -> [ j[1],-j[3],-j[2],-j[4]]:

 U["vertex_convert"] := (i) -> [i[1],i[2],0,i[3]]:
 U["edge_convert"]   := (e) -> sort(map(U["vertex_convert"],e)):
 U["face_convert"]   := (f) -> sort(map(U["vertex_convert"],f)):

 U["embedding0"] := table():
 U["embedding1"] := table():
 U["embedding2"] := table():
 U["embedding3"] := table():
 U["embedding4"] := table():
 U["vertices"] := map(U["vertex_convert"],T["vertices"]):
 for i in T["vertices"] do
  U["embedding0"][U["vertex_convert"](i)] := evalf(T["embedding0"][i]);
  U["embedding1"][U["vertex_convert"](i)] := evalf(T["embedding1"][i]);
  U["embedding2"][U["vertex_convert"](i)] := evalf(T["embedding2"][i]);
  U["embedding3"][U["vertex_convert"](i)] := evalf(T["embedding3"][i]);
  U["embedding4"][U["vertex_convert"](i)] := evalf(T["embedding4"][i]);
 od:

 for i in U["vertices"] do
  for g in G["elements"] do
   U["embedding0"][U["act_i"][g](i)] := evalf(T["act_R3"][g](U["embedding0"][i]));
   U["embedding1"][U["act_i"][g](i)] := evalf(T["act_R3"][g](U["embedding1"][i]));
   U["embedding2"][U["act_i"][g](i)] := evalf(T["act_R3"][g](U["embedding2"][i]));
   U["embedding3"][U["act_i"][g](i)] := evalf(T["act_R3"][g](U["embedding3"][i]));
   U["embedding4"][U["act_i"][g](i)] := evalf(T["act_R3"][g](U["embedding4"][i]));
  od:
 od:
 U["embedding"] := eval(U["embedding4"]);

 U["vertices"] := sort(map(op,[indices(U["embedding0"])])):
 U["edges"] :=
  {seq(seq(sort(map(U["act_i"][g],U["edge_convert"](e))),
    e in T["edges"]),g in G["elements"])}:
 U["edges"] := sort([op(U["edges"])]):
 U["faces"] :=
  {seq(seq(sort(map(U["act_i"][g],U["face_convert"](f))),
    f in T["faces"]),g in G["elements"])}:
 U["faces"] := sort([op(U["faces"])]):

 R := T["act_R3"][[1,2,3,-1]];
 f2 := T["f2"];

 U["smooth_plot"] :=
  display(
   plot3d(  f2(sqrt(1-r^2) *~ T["u"][7] +~ (r * cos(t)) *~ T["u"][1] +~ (r * sin(t)) *~ T["u"][9]),
          r=0..0.6,t=0..2*Pi,style=patchnogrid),
   plot3d(R(f2(sqrt(1-r^2) *~ T["u"][7] +~ (r * cos(t)) *~ T["u"][1] +~ (r * sin(t)) *~ T["u"][9])),
          r=0..0.6,t=0..2*Pi,style=patchnogrid),
   plot3d(  f2(cos(s) *~ T["u"][7] +~ (sin(s) * cos(t)) *~ T["u"][1] +~ (sin(s) * sin(t)) *~ T["u"][9]),
          s=0.2*Pi..0.8*Pi,t=0..2*Pi,style=patchnogrid,numpoints=4000),
   plot_opts
  );

 U["chart_plot"] := display(
  plot3d(xa,r=0..0.5,t=0..2*Pi,style=patchnogrid),
  plot3d(T["act_R3"][[1,2,3,-1]](xa),r=0..0.5,t=0..2*Pi,style=patchnogrid),
  plot3d(xb,u=0..1.4,t=0..2*Pi,style=patchnogrid),
  plot3d(T["act_R3"][[2,3,1,1]](xb),u=0..1.4,t=0..2*Pi,style=patchnogrid),
  plot3d(T["act_R3"][[3,1,2,1]](xb),u=0..1.4,t=0..2*Pi,style=patchnogrid),
  plot3d(xc,r=0..0.3,t=0..2*Pi,style=patchnogrid),
  plot3d(T["act_R3"][[2,3,1,1]](xc),r=0..0.3,t=0..2*Pi,style=patchnogrid),
  plot3d(T["act_R3"][[3,1,2,1]](xc),r=0..0.3,t=0..2*Pi,style=patchnogrid),
  scaling=constrained,axes=none):

 if verbose then printf("Trimming\n"); fi;

 V := table():
 pants_complex := eval(V);
 V["vertices"] := select(i -> max(i[1],-i[2],i[3]) <= 1/2,U["vertices"]):
 V["edges"] := select(e -> {op(e)} minus {op(V["vertices"])} = {}, U["edges"]):
 V["faces"] := select(f -> {op(f)} minus {op(V["vertices"])} = {}, U["faces"]):
 V["max_simplices"] := V["faces"]:
 V["embedding0"] := table():
 V["embedding1"] := table():
 V["embedding2"] := table():
 V["embedding3"] := table():
 V["embedding4"] := table():
 for i in V["vertices"] do
  V["embedding0"][i] := U["embedding0"][i];
  V["embedding1"][i] := U["embedding1"][i];
  V["embedding2"][i] := U["embedding2"][i];
  V["embedding3"][i] := U["embedding3"][i];
  V["embedding4"][i] := U["embedding4"][i];
 od:
 V["embedding"] := V["embedding4"]:

 V["normal"] := table():
 for i in V["vertices"] do
  n := evalf(T["dl"](V["embedding"][i]));
  n := n /~ dp(n,n);
  V["normal"][i] := n;
 od:

 V["chart_plot"] :=
  display(
   plot3d(xa,r=0..0.5,t=0..2*Pi,style=patchnogrid),
   plot3d(T["act_R3"][[1,2,3,-1]](xa),r=0..0.5,t=0..2*Pi,style=patchnogrid),
   plot3d(xb,u=0..sqrt(3.)/2,t=0..2*Pi,style=patchnogrid),
   plot3d(T["act_R3"][[2,3,1,1]](xb),u=0..sqrt(3.)/2,t=0..2*Pi,style=patchnogrid),
   plot3d(T["act_R3"][[3,1,2,1]](xb),u=0..sqrt(3.)/2,t=0..2*Pi,style=patchnogrid),
  scaling=constrained,axes=none):

 condensed_pants_complex :=
  `condense/simplicial_complex`(pants_complex);

 if verbose then printf("Generating javascript\n"); fi;
 `set_javascript/simplicial_complex`(condensed_pants_complex);
end:

`base/pants` := table([
 "centre" = [0,0,0],
 "circle_centres" = [[sqrt(3)/2,0,0],[-sqrt(3)/4,3/4,0],[-sqrt(3)/4,-3/4,0]],
 "top" = [0,0,1/2]
]):

`outline/pants` := proc(P)
 local c,z,t;
 c := P["centre"];
 z := 2 *~ (P["top"] -~ P["centre"]);
 display(
  seq(line(P["centre"],u,colour=green),u in P["circle_centres"]),
  seq(spacecurve(u +~ (cos(t)/sqrt(3)) *~ cross_product(u -~ c,z) +~ (sin(t)/2) *~ z,t=0..2*Pi,colour=red),
       u in P["circle_centres"]),
  scaling=constrained,axes=none
 );
end:

`rough_plot/pants` := proc(P)
 local c,z,t,u,v,w,X;
 c := P["centre"];
 z := 2 *~ (P["top"] -~ P["centre"]);
 X := NULL;
 for u in P["circle_centres"] do
  v := evalf((u -~ c) *~ (2/sqrt(3)));
  w := cross_product(v,z);
  X := X,plot3d(c +~ (s * abs(cos(t))/4 + (1-s) * sqrt(3)/2) *~ v
                +~ (cos(t)/2) *~ w +~ (sin(t)/2) *~ z,
                s=0..1,t=0..2*Pi,style=patchnogrid,scaling=constrained,axes=none);
 od:
 return(display(X));
end:

`plot/pants` := proc(P)
 local dp,h,c,b,u,v,w,p,f,X;

 dp := (u,v) -> add(u[i] * v[i],i=1..3);
 h := (x) -> evalf(x /~ sqrt(dp(x,x)));
 c := evalf(P["centre"]);
 b := map(x -> x -~ c,P["circle_centres"]);
 u := h(b[1]);
 v := h(b[2] -~ b[3]);
 w := h(P["top"] -~ c);
 p := pants_sphere_sector_complex["u"];
 f := unapply(c +~ dp(x,p[1]) *~ u +~ dp(x,p[9]) *~ v +~ dp(x,p[7]) *~ w,x);
 X := `map/simplicial_complex`(f,condensed_pants_complex);
 return `surface_plot/simplicial_complex`(X);
end:

`map/pants` := proc(f,P)
 return table([
  "centre" = f(P["centre"]),
  "circle_centres" = map(f,P["circle_centres"]),
  "top" = f(P["top"])
 ]);
end:

`sprout/pants` := proc(P,i)
 local Q,u,c;
 c := P["centre"];
 u := P["circle_centres"][i] -~ c;
 Q := table():
 Q["centre"] := c +~ 2 *~ u;
 Q["circle_centres"] := map(v -> 2 *~ c + 2 *~ u -~ v,P["circle_centres"]);
 Q["top"] := P["top"] -~ P["centre"] +~ Q["centre"];
 return eval(Q);
end:

`pants_tube/base` := table([
  "centre"    = [0,0,0],
  "direction" = [1,0,0],  # unit vector
  "bend"      = [0,1,0],  # unit vector
  "length"    = 1,
  "curvature" = 0
]):

`outline/pants_tube` := proc(T)
 local c,d,b,l,k,w,x,y,z,s,t0,t;
 c := T["centre"];
 d := T["direction"];
 b := T["bend"];
 l := T["length"];
 k := T["curvature"];
 w := cross_product(d,b);
 x := c -~ (`C/pants_tube`(t*l*k/2)*t^2*l^2*k/8) *~ b +~ (`S/pants_tube`(t*l*k/2)*t*l/2) *~ d;
 y := cos(t*l*k/2) *~ b +~ sin(t*l*k/2) *~ d;
 z := x +~ (cos(s)/2) *~ y +~ (sin(s)/2) *~ w;
 return display(
  spacecurve(x,t = -1..1,colour=cyan),
  spacecurve(subs(t=-1,z),s=0..2*Pi,colour=orange),
  spacecurve(subs(t= 1,z),s=0..2*Pi,colour=orange),
  scaling=constrained,axes=none
 );
end:

`S/pants_tube` := unapply(convert(series(sin(x)/x,x=0,11),polynom,x),x);
`C/pants_tube` := unapply(convert(series(2*(1-cos(x))/x^2,x=0,11),polynom,x),x);
`A/pants_tube` := unapply(convert(series(arcsin(x)/x,x=0,11),polynom,x),x);

`plot/pants_tube` := proc(T)
 local c,d,b,l,k,w,x,y,z,s,t;
 c := T["centre"];
 d := T["direction"];
 b := T["bend"];
 l := T["length"];
 k := T["curvature"];
 w := cross_product(d,b);
 x := c -~ (`C/pants_tube`(t*l*k/2)*t^2*l^2*k/8) *~ b +~ (`S/pants_tube`(t*l*k/2)*t*l/2) *~ d;
 y := cos(t*l*k/2) *~ b +~ sin(t*l*k/2) *~ d;
 z := x +~ (cos(s)/2) *~ y +~ (sin(s)/2) *~ w;
 return display(
  plot3d(z,t=-1..1,s=0..2*Pi,style=patchnogrid),
  scaling=constrained,axes=none
 );
end:

`map/pants_tube` := proc(f,T)
 local U;
 U := table():
 U["centre"] := f(T["centre"]);
 U["direction"] := f(T["centre"] +~ T["direction"]) -~ U["centre"];
 U["bend"]      := f(T["centre"] +~ T["bend"]     ) -~ U["centre"];
 U["length"]    := T["length"];
 U["curvature"] := T["curvature"];
 return eval(U);
end:

`joiner/pants_tube` := proc(c1,n1,c2,n2_)
 local dp,d,s,c0,ca,sa,n2,err,p,q,k,l,b,c,T,alpha,m,pic;
 dp := `dot/R`(3);

 d := evalf(c2 -~ c1);
 s := sqrt(dp(d,d));           # straight line length
 d := d /~ s;                  # unit vector c1 -> c2
 c0 := evalf((c1 +~ c2) /~ 2); # straight line midpoint

 ca := -evalf(dp(n1,d));       # cos(alpha/2)
 sa := sqrt(1 - ca * ca);      # sin(alpha/2)

 if nargs >= 4 then
  n2 := n2_;
 else
  n2 := n1 +~ (2 * p) *~ d;
 fi;

 k := 2 * sa/s;       # curvature
 q := s/2 * sqrt((1-ca)/(1+ca));

 err := evalf(cross_product(d,n2 -~ n1));
 err := sqrt(dp(err,err));
 if err > 10.^(-6) then
  error("Non-matching circles");
 fi;

 b := -~ evalf(n1 +~ n2);

 if sqrt(dp(b,b)) < 10.^(-4) then
  # Tube will be straight.  Set b to be an arbitrary vector
  # orthogonal to d.

  b := [1,0,0]; p := abs(d[1]);
  if abs(d[2]) < p then b := [0,1,0]; p := abs(d[2]); fi;
  if abs(d[3]) < p then b := [0,0,1]; p := abs(d[3]); fi;
  b := b -~ dp(d,b) *~ d;
  b := b /~ sqrt(dp(b,b));
  k := 0;
  c := c0;
  l := s;
 else
  b := b /~ sqrt(dp(b,b));
  c := evalf(c0 +~ q *~ b);
#  l := s * `A/pants_tube`(s * k / 2);
  alpha := 2 * arccos(ca);
  l := alpha / k;
  m := c -~ b /~ k;
  pic := display(
   point(c1,colour=black),
   line(c1,c2,colour=black),
   arrow(c1,0.3*~n1,cross_product(b,d),0.05,0.1,0.025,colour=magenta),
   arrow(c2,0.3*~n2,cross_product(b,d),0.05,0.1,0.025,colour=magenta),
   line(c0,c,colour=orange),
   line(m,c1,colour=purple),
   line(m,c2,colour=purple),
   spacecurve(m +~ (cos(t)/k) *~ b +~ (sin(t)/k) *~ d,t=-alpha/2..alpha/2,colour=cyan),
   scaling=constrained,axes=none
  );
 fi;

 T := table([
  "c1" = evalf(c1), "c2" = evalf(c2), "n1" = evalf(n1), "n2" = evalf(n2), "d" = d, "s" = s, "c0" = c0, "ca" = ca, "sa" = sa, "q" = q, "alpha" = alpha,
  "centre" = c,
  "direction" = d,
  "bend" = b,
  "length" = l,
  "curvature" = k,
  "pic" = pic
 ]);

 return eval(T);
end:

`pants_joiner/pants_tube` := proc(P1,i1,P2,i2)
 local dp,a1,a2,c1,c2,n1,n2;
 dp := `dot/R`(3);
 a1 := P1["centre"];
 a2 := P2["centre"];
 c1 := P1["circle_centres"][modp(i1-1,3)+1];
 c2 := P2["circle_centres"][modp(i2-1,3)+1];
 n1 := a1 -~ c1;
 n2 := a2 -~ c2;
 n1 := n1 /~ sqrt(dp(n1,n1));
 n2 := n2 /~ sqrt(dp(n2,n2));
 return `joiner/pants_tube`(c1,n1,c2,n2);
end:


text:
input := FileTools:-Text:-ReadFile("AoC-2021-14-input.txt" ):

split := StringTools:-Split(input,"\n"):
polymer := split[1];
rules := table(map(`=`@op,
    map(s->StringTools:-Split(s)[[1,3]], split[3..-1]))):

# initial count
paircount := table(sparse=0):
for i to length(polymer)-1 do
   paircount[polymer[i..i+1]] += 1;
end do:

for i to 40 do
    update := indices(paircount, 'pairs');
    paircount := table(sparse=0):
    lettercount := table(sparse=0):
    for u in update do
        b := lhs(u); c := rules[lhs(u)];
        paircount[ cat(b[1],c) ] += rhs(u);
        paircount[ cat(c,b[2]) ] += rhs(u);
        lettercount[ b[1] ] += rhs(u);
        lettercount[ b[2] ] += rhs(u);
        lettercount[ c ] += 2*rhs(u);
    end do:
    lettercount[polymer[1]] += 1; lettercount[polymer[-1]] += 1;
    if i = 10 or i = 40 then
        t:=[entries(lettercount, 'nolist')];
        print( i=(t[max[index](t)] - t[min[index](t)]) / 2);
    end if;
end do:


text:
(*
 Copyright (C) 2017 M.A.L. Marques

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*)

(* nop scaling *)
s_scaling_0 := s -> s:

(* original scaling of Heyd *)
strans :=  8.3:
smax   :=  8.572844:
sconst := 18.79622316:
s_scaling_1 := s -> my_piecewise3(
  s < strans, s,
  smax - sconst/s^2
):

(* first version of the scaling by TM Henderson, apparently used by Gaussian *)
s_scaling_2 := s -> my_piecewise3(
  s < 1,  s,
  my_piecewise3(s > 15, smax, m_max(m_min(s, 15), 1) - log(1 + exp(m_max(m_min(s, 15), 1) - smax)))
):

(* second version of the scaling by TM Henderson *)
s_scaling_3 := s -> s - (1 - exp(-s))*log(1 + exp(s - smax)):

(* appendix of JCP 128, 194105 (2008) *)
s_p :=  [0.615482, 1.136921, -0.449154, 0.0175739*8.572844]:
s_q :=  [1.229195, -0.0269253, 0.313417, -0.0508314, 0.0175739]:

s_scaling_4 := s->
  s*(1 + s^3*add(s_p[i]*s^i, i=1..4))/(1 + s^3*add(s_q[i]*s^i, i=1..5)):


text:
# prepare lib
restart:
libname:=libname,"..":
with(MapleCodeReader):
# prepare data
ReadCode("../InvClassify/Basic.mpl");
setLogLevel(3);
vv:=[d(x), d(t), u*d(u), x*d(x)+2*t*d(t), 2*t*d(x)-x*u*d(u), 4*t*x*d(x)+4*t^2*d(t)-(x^2+2*t)*u*d(u)]:
As,A,eqs:=getTransMatAndPDE(vv):
# run test
ReadCode("../InvClassify/Closure.mpl");
map(findClosure,[seq(1..4)],[{2,3},{2},{3},{1,4}]);
getClosure(A);
getClosure(A,[a[2]=0,a[4]=0,a[6]=0]);
getClosure(A,[a[2]=0,a[4]=0,a[6]=0,a[1]=0,a[5]=0]);


text:
with(LinearAlgebra): with(ArrayTools): with(QuantumChemistry):
loaddata := readdata(cat(cdir,"_temp.rdm"), 8):
Flatten := proc(x) local n, a, i, j, k, l; `local`(a, n, i, j, k, l); `description`("convert chemists to numpy and flatten an array to form a matrix"); n := Size(x); a := Matrix(1 .. n[1]*n[2], 1 .. n[3]*n[4], datatype = float[8]); for i to round(n[1]) do for j to round(n[2]) do for k to round(n[3]) do for l to round(n[4]) do a[(i - 1)*n[1] + j, (k - 1)*n[3] + l] := x[i, j, k, l]; end do; end do; end do; end do; return a; end proc:
New := Array(1 .. 4, 1 .. 4, 1 .. 4, 1 .. 4, datatype = float[8]):
for i in loaddata[3 .. ()] do
    New[round(i[1]), round(i[2]), round(i[3]), round(i[4])] := i[5]:
end do:
pure := Purify2RDM(New, spin_free = false, electron_number = 4, conv_tol = 0.10000000):
ExportMatrix(cat(cdir, "_temp_purified.csv"), Flatten(pure[rdm2])):

max_stars_repo_name: damazz/HQCA

max_issues_repo_name: damazz/HQCA

max_forks_repo_name: damazz/HQCA


text:
(*
 Copyright (C) 2017 M.A.L. Marques

 This Source Code Form is subject to the terms of the Mozilla Public
 License, v. 2.0. If a copy of the MPL was not distributed with this
 file, You can obtain one at http://mozilla.org/MPL/2.0/.
*)

(* type: work_gga_x *)

$define gga_x_pbe_sol_params
$include "gga_x_pbe.mpl"

cc := 100:
c1 := 0.5217:

f1 := s -> f0_pbe(s)*(cc - s^4) + c1*s^3.5*(1 + s^2):
f2 := s -> cc + s^6:

f  := x -> f1(X2S*x)/f2(X2S*x):

max_stars_repo_name: rdietric/lsms

max_issues_repo_name: rdietric/lsms

max_forks_repo_name: rdietric/lsms


text:
`is_element/graphs` := (V::set) -> proc(E)
 if not(`is_element/autorel`(V)(E)) then
  return false;
 fi;

 if not(`is_symmetric/autorel`(V)(E)) then
  return false;
 fi;

 if not(`is_irreflexive/autorel`(V)(E)) then
  return false;
 fi;

 return true;
end:

`is_leq/graphs` := (V::set) -> (E0,E1) -> evalb(E0 minus E1 = {}):

`list_elements/graphs` := proc(V::set)
 local n,E0,EE,i,j;

 n := nops(V);
 E0 := {seq(seq([V[i],V[j]],j=i+1..n),i=1..n-1)};
 EE := `list_elements/subsets`(E0);
 EE := map(E -> E union `op/autorel`(V)(E),EE);
 return EE;
end:

`count_elements/graphs` := (V::set) -> 2^(nops(V)*(nops(V)-1)/2);

`random_element/graphs` := (V::set) -> proc()
 local n,E0,E,i,j;

 n := nops(V);
 E0 := {seq(seq([V[i],V[j]],j=i+1..n),i=1..n-1)};
 E := `random_element/subsets`(E0)();
 E := E union `op/autorel`(V)(E);
 return E;
end:

`neighbour_table/graphs` := (V) -> proc(E)
 local N,v,e;

 N := table();
 for v in V do N[v] := {}; od:
 for e in E do
  N[e[1]] := {op(N[e[1]]),e[2]};
 od:

 return eval(N);
end:

`component_relation/graphs` := (V) -> (E) ->
 `transitive_closure/autorel`(V)(
  E union `id/autorel`(V)
 );

`components/graphs` := (V) -> (E) ->
 `block_partition/equiv`(`component_relation/graphs`(V)(E));

`is_connected/graphs` := (V) -> (E) ->
 evalb(`component_relation/graphs`(V)(E) = `top/autorel`(V));

`is_forest/graphs` := (V) -> proc(E)
 local CC,V0,E0;

 CC := `components/graphs`(V)(E);

 for V0 in CC do
  E0 := select(e -> member(e[1],V0),V0);

  if nops(E0) <> nops(V0) - 1 then
   return false;
  fi;
 od:

 return true;
end:

`is_tree/graphs` := (V) -> proc(E)
 `is_connected/graphs`(V)(E) and
  nops(E) = nops(V) - 1;
end:

`skeleton/graphs` := (V::set) -> proc(E)
 local n,V0,E0,ix,i;

 n := nops(V);
 V0 := {seq(i,i=1..n)};
 ix := table():
 for i from 1 to n do ix[V[i]] := i; od;
 E0 := map(e -> [ix[e[1]],ix[e[2]]],E);
 return [V0,E0];
end:

`is_element/paths` := (V) -> (E) -> proc(p)
 local i,n;

 if not type(p,list) then
  return false;
 fi;

 n := nops(p) - 1;

 if n < 0 then
  return false;
 fi;

 if {op(p)} minus V <> {} then
  return false;
 fi;

 for i from 1 to n do
  if not(member([p[i],p[i+1]],E)) then
   return false;
  fi;
 od;

 return true;
end:

`is_small_element/paths` := (V) -> (E) -> (n::nonnegint) -> proc(p)
 return evalb(`is_element/paths`(V)(E)(p) and (nops(p) = n+1));
end:

`random_small_element/paths` := (V) -> (E) -> (n::nonnegint) -> proc()
 local v,p,e,i;

 if E = {} then
  if V = {} or n > 0 then
   return FAIL;
  else
   return [random_element_of(V)()];
  fi;
 fi;

 v := random_element_of(E)()[1];
 p := v;
 for i from 1 to n do
  e := random_element_of(select(e -> member(v,e),E))();
  v := op({op(e)} minus {v});
  p := p,v;
 od;
 return [p];
end:

`list_small_elements/paths` := (V) -> (E) -> proc(n::nonnegint)
 option remember;
 local P,Q,p,v,F,e;

 if n = 0 then
  return map(v -> [v],V);
 else
  P := `list_small_elements/paths`(V)(E)(n-1);
  Q := NULL;
  for p in P do
   v := p[1];
   F := select(e -> e[2] = v,E);
   Q := Q,seq([e[1],op(p)],e in F);
  od:
  return [Q];
 fi;
end:

`is_leq/paths` := NULL;
`count_elements/paths` := NULL:
`count_small_elements/paths` := NULL:

`length/paths` := (V) -> (E) -> (p) -> nops(p) - 1;

`is_trail/paths` := proc(p)
 local n,e,i;
 n := nops(p) - 1;
 e := {seq([p[i],p[i+1]],i=1..n)};
 return evalb(nops(e) = n);
end:

`is_cycle/paths` := proc(p)
 local n,i;
 n := nops(p) - 1;
 return evalb(p[1] = p[n+1] and nops({seq(p[i],i=1..n)}) = n);
end:

`is_small_element/trails` := (V) -> (E) -> (n) -> proc(p)
 if not `is_small_element/paths`(V)(E)(n)(p) then
  return false;
 fi;

 return `is_trail/paths`(p);
end:

`is_element/trails` := (V) -> (E) -> proc(p)
 if not `is_element/paths`(V)(E)(p) then
  return false;
 fi;

 return `is_trail/paths`(p);
end:

`list_small_elements/trails` := (V) -> (E) -> proc(n)
 option remember;
 local E0,E1,P,Q,p,v,F,i,e;

 if n = 0 then
  return map(v -> [v],V);
 else
  P := `list_small_elements/trails`(V)(E)(n-1);
  Q := NULL;
  for p in P do
   E0 := {seq([p[i],p[i+1]],i=1..n-1),seq([p[i+1],p[i]],i=1..n-1)};
   E1 := E minus E0;
   v := p[1];
   F := select(e -> e[2] = v,E1);
   Q := Q,seq([e[1],op(p)],e in F);
  od:
  return [Q];
 fi;
end:

`random_small_element/trails` := (V) -> (E) -> (n) -> proc(num_tries := 10)
 local k,p,E0,E1,v,F,e,i;

 k := num_tries;

 if n = 0 then
  if V = {} then
   return FAIL;
  else
   return [random_element_of(V)()];
  fi;
 elif n = 1 then
  if E = {} then
   return FAIL;
  else
   return random_element_of(E)();
  fi;
 fi;

 while k > 0 do
  p := `random_small_element/trails`(V)(E)(n-1)();
  k := k - 1;
  if p <> FAIL then
   E0 := {seq([p[i],p[i+1]],i=1..n-1),seq([p[i+1],p[i]],i=1..n-1)};
   E1 := E minus E0;
   v := p[1];
   F := select(e -> e[2] = v,E1);
   if F <> {} then
    e := random_element_of(F)();
    return [e[1],op(p)];
   fi;
  fi;
 od;

 return FAIL;
end:

`is_leq/trails` := NULL:
`count_elements/trails` := NULL:
`count_small_elements/trails` := NULL:

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